On the average Gamma-Ray Burst X-ray flaring activity

Gamma-ray burst X-ray flares are believed to mark the late time activity of the central engine. We compute the temporal evolution of the average flare luminosity $$ in the common rest frame energy band of 44 GRBs taken from the large \emph{Swift} 5-years data base. Our work highlights the importance of a proper consideration of the threshold of detection of flares against the contemporaneous continuous X-ray emission. In the time interval $30 \rm{s}\propto t^{-2.7\pm 0.1}$; this implies that the flare isotropic energy scaling is $E_{\rm{iso,flare}}\propto t^{-1.7}$. The decay of the continuum underlying the flare emission closely tracks the average flare luminosity evolution, with a typical flare to steep-decay luminosity ratio which is $L_{\rm{flare}}/L_{\rm{steep}}=4.7$: this suggests that flares and continuum emission are deeply related to one another. We infer on the progenitor properties considering different models. According to the hyper-accreting black hole scenario, the average flare luminosity scaling can be obtained in the case of rapid accretion ($t_{\rm{acc}}\ll t$) or when the last $\sim 0.5 M_{\sun}$ of the original $14 M_{\sun}$ progenitor star are accreted. Alternatively, the steep $\propto t^{-2.7}$ behaviour could be triggered by a rapid outward expansion of an accretion shock in the material feeding a convective disk. If instead we assume the engine to be a rapidly spinning magnetar, then its rotational energy can be extracted to power a jet whose luminosity is likely to be between the monopole ($L\propto e^{-2t}$) and dipole ($L\propto t^{-2}$) cases. In both scenarios we suggest the variability, which is the main signature of the flaring activity, to be established as a consequence of different kinds of instabilities.

(ii) Lifetime. CEs are required to be long-lived and in particular to be active up to t ∼ 10 4 − 10 5 s after the end of the prompt emission. The probability of a CE to re-start decays with time.
(iii) Variability. CEs are required to store and release energy in the form of erratic, short-lived episodes of emission giving rise to flares whose duration linearly grows with time.
Ideas on how to revive the CE have been explored by a number of authors: flares could be the result of a two-stage stellar collapse leading to core fragmentation and subsequent accretion (King et al. 2005). Alternatively, the fragmentation of the outer part of an hyper-accreting disk into rings of material as a result of gravitational instabilities could potentially lead to large-amplitude variations of the CE output (Perna, Armitage & Zhang 2006). However, Piro & Pfhal (2007) noticed that tidal disruption of a fragment and accretion by the central black hole might be too rapid to account for the durations of observed flares. Proga & Zhang (2006) proposed that magnetic flux accumulated around the accretor could stop and then restart the CE release of energy possibly giving rise to flares. Dai et al. (2006) invoked the presence of a differentially rotating millisecond pulsar where magnetic reconnection events are responsible for the observed flare emission. Magnetic reconnection is also suggested by Giannios (2006) in the framework of the internal plus external shock scenario with strongly magnetised ejecta: in this case, no late-time CE activity is needed and flares come from the revival of MHD instabilities triggered by the deceleration of the magnetised ejecta in the external medium.
The entire list of flare models basically refers to two different CEs: an accretion disk around a newly formed black hole (e.g. MacFadyen & Woosley 1999) or a rapidly rotating magnetar (e.g. Usov 1992;Thompson 1994;Wheeler et al. 2000; Thompson et al. 2004). In magnetar models the ultimate source of energy is represented by the magnetar spin and flares are likely to be connected to episodes of magnetic energy dissipation. In the case of accretion models, the source of energy powering the flare emission is possibly of gravitational origin. However, magnetic effects could still play the key role, deeply affecting the accretion dynamics and release of energy (e.g. Fan et al. 2005;Proga & Zhang 2006;Lei et al. 2009;Zhang & Dai 2009). X-ray flares mark the late time activity of the CE after the main power emission and provide information about the way the CE power is progressively declining. In the following we show that this information turns out to be of fundamental importance in constraining existing CE models (accretion or magnetar).
The main goal of this work is to compute the evolution of the average luminosity of the flaring component with time, taking advantage from the large Swift 5-years database. A variety of physical mechanisms which might properly account for the new episodes of energy release are discussed in the context of both accretion and magnetar CE models. A previous attempt in this direction was made by Lazzati et al. (2008) (L08 hereafter): from a sample of 9 long GRBs and 1 short GRB (which was the largest sample with known redshift available at that time) these authors showed that the average flare luminosity follows a power-law decay L L08 ∝ t −1.5±0.2 and concluded that accretion onto a compact object could reasonably account for the L decay. In the present work we aim at presenting a substantial update of the L08 sample, including 44 GRBs with redshift; the flare luminosity of each GRB is computed in the rest frame energy band which is common to the entire sample: this assures L to be obtained from strictly comparable properties of each GRB. Moreover, no assumption is made on the functional form of the flare temporal profile. Finally, since flares are detected as sharp features superimposed to the X-ray afterglow, it is of fundamental importance to understand the associated observational biasses: we compute the temporal evolution of the flare detection threshold for each GRB, and quantitatively discuss the influence of the flare detection threshold on the obtained L temporal evolution. This work is organised as follows: the sample selection and data reduction is detailed in Sect. 2 while the computation of the flare and underlying continuum components of each event is presented in Sect. 3. Section 4 describes the method which leads to the definition of a flare detection threshold of each GRB. A luminosity-time plane avoidance region is presented in Sect. 5. Results are discussed in Sect. 6: possible physical interpretations in the context of both accretion and magnetar models of the central GRB engine are critically analysed in Sect. 6.1.1 and 6.1.2. Constraints to physical scenarios able to account for the erratic behaviour of the flare emission are derived in Sect. 6.2. Conclusions are drawn in Sec. 7.
The phenomenology of the bursts is presented in the rest frame unless otherwise stated. Uncertainties are quoted at the 68% confidence level (c.l.): a warning is added if it is not the case. Standard cosmological quantities have been adopted: H0 = 70 Km s −1 Mpc −1 , ΩΛ = 0.7, ΩM = 0.3.

SAMPLE SELECTION AND DATA REDUCTION
We select the GRBs which visually show flaring activity superimposed to the smooth X-ray afterglow, good data quality and redshift, detected by the Swift X-ray telescope (XRT, Burrows et al. 2005) from the beginning of the mission up to the end of February, 2010. The sample comprises 44 GRBs (Table 1), with redshift ranging from 0.3 to 6.3 1 . XRT data have been processed with the heasoft package v. 6.6.1 and corresponding calibration files: standard filtering and screening criteria have been applied. Piled-up Window Timing (WT) data have been corrected following the prescriptions by Romano et al. (2006), while piled-up Photon Counting (PC) data have been extracted from an annular region whose inner radius has been derived comparing the observed to the nominal point spread function (PSF, Moretti et al. 2005;Vaughan et al. 2006). The background is estimated from a source-free portion of the sky and then subtracted. The 0.3-10 keV background subtracted, PSF and vignetting corrected light-curve of each GRB has been rebinned so as to assure a minimum signal-to-noise (SN) equal to 4. The count-rate light-curves are calibrated into luminosity light-curves using a time dependent count-to-flux conversion factor as described in Margutti et al. (2010a). This procedure produces luminosity curves where the possible spectral evolution of the source is properly taken into account. Each 0.3-10 keV XRT light-curve is calibrated in the common rest frame energy band defined by the redshift distribution of the sample which turns out to be 2.2-14.4 keV. This allows us to make a direct comparison between the X-ray afterglows of different bursts while avoiding extrapolation of the signal to an unobserved energy band 2 .

CONTINUUM AND EXCESS ESTIMATION
We produced a software to automatically identify the smooth continuum underlying the X-ray afterglow of GRBs with superimposed flaring activity. The procedure is based on the χ 2 statistics and can be described as a two-step process: (i) Identify the continuum afterglow component underlying the flaring activity; (ii) Calculate the flaring contribution as excess with respect to the continuum.
A first blind fit of the entire X-ray afterglow (continuum plus flares) is done in log-log units using a power-law, smoothly joint broken power-law or double broken power-law models. If the P-value (probability of obtaining a result at least as extreme as the one that is actually observed) associated to this fit is lower than 5%, the data point with the largest positive residual is removed from the light-curve and a new fit is performed. This process is repeated until a P-value > 5% is obtained. The F-test is used to choose between the different nested models when necessary. The best fitting model satisfying the P-value condition is identified with the underlying continuum associated to a particular GRB X-ray afterglow (red dot-dashed line in Fig. 1 for the case of GRB 060607A) Figure 1. X-ray afterglow of GRB 060607A in the rest frame energy band 2.2-14.4 keV. Red dot-dashed line: best estimate of the continuum underlying the flaring emission obtained as detailed in Sec. 3. Blue solid line: 2σ threshold of flare detection (i.e. the minimum peak luminosity that a flare must have to be detected as statistically significant fluctuation superimposed to the X-ray continuum), calculated as described in Sec. 4. In this case, for t 1000 s, the threshold is found to be ∼ 50% the continuum. Inset: residuals with respect to the continuum. and is subtracted from the original light-curve. The resulting residuals and respective errors (both the statistical uncertainty associated to the original light-curve bins and the one coming from the continuum estimation 3 are properly taken into account and propagated) constitute the candidate X-ray flaring component associated to a GRB: this is shown in the inset of Fig. 1 for GRB 060607A. Note that no particular flare functional shape is assumed in this analysis. This method allows us to account for small variations superimposed to the continuum: in C10 flares were instead visually identified and then fitted with a specific profile.
The average GRB 2.2-14.4 keV rest frame flaring component for any time interval ti − t f is computed as: where N is the number of GRB displaying a positive flaring component at a minimum 2σ significance during ti − t f ; Li is the 2.2-14.4 keV luminosity of the flaring component of the i th burst evaluated at at t = (ti + t f )/2. Linear interpolation is used when necessary. The uncertainty affecting L is found through standard propagation of the uncertainties affecting each flaring component, as determined in the previous paragraph. The typical uncertainty on each GRB flaring component data point is of the order of 30% the value of the GRB afterglow light-curve before the flare subtraction. An example is shown in Fig. 1, inset. The use of the logarithmic mean in Eq. 1 prevents L from being dominated by a single bright event (the use of the median has been proven to lead to very similar results); instead, the linear mean of the excesses leads to L values biassed towards the bright end of the luminosity distribution of the excesses at any time t, and is therefore discarded. GRB 050904 shows an extended In the time interval t < 1000 s (t > 1000 s) the best fit reads L/erg s −1 = 10 54.5±0.1 (t/s) −2.7±0.1 ( L/erg s −1 = 10 51.4±0.1 (t/s) −1.2±0.1 ). Panel (b): average re-normalised luminosity curve: the luminosity has been re-scaled for the average prompt 15-350 keV luminosity observed during the T 90 . The xaxis is expressed in T 90 units. In the time interval 1 t * 100 the re-normalised flare luminosity follows a power-law decay with best fitting model: flaring activity between 500 s and 3.1 ks: given the high luminosity of the flaring component and the peculiarity of a flaring emission covering the entire plateau phase duration, GRB 050904 is not included in the calculation above. This guarantees that L is not biassed by the contribution of a single peculiar burst: however, the inclusion of this event would not change the main conclusions of this work. The flare luminosity curve of GRB 050904 is investigated and discussed in Sec. 6. The result of the application of Eq. 1 is shown in Fig. 2, panel (a): in the time interval 15 − 30 s the average light-curve of the flaring component is likely to track the prompt gamma-ray emission, while for 30 < t < 1000 s the average light-curve is well represented by a powerlaw with best-fitting index α1 = −2.7 ± 0.1. For t > 1000 s the best-fitting power-law index reads α2 = −1.2 ± 0.1. We note that a simple power-law modelling of the entire 15 s < t < 10 6 s temporal window yields a statistically unacceptable fit (χ 2 /dof = 1271.4/50, P-value≪ 10 −10 ) that systematically overestimates the computed average flare luminosity for t > 1 ks. The two data points at t ∼ 10 3 s largely above the detection threshold (green thick area, see Sec. 4 for details) are mainly due to the flaring contribution of GRB 081028, GRB 090809 and GRB 090417B: these 3 events show a large flare at this rest frame epoch. The exclusion of these 3 bursts from Eq. 1 would bring the two data points slightly above the threshold. The shown error bars account for the 1σ uncertainty affecting the average luminosity value: however, the dispersion of the sample luminosity values at any given time t could be remarkably different. For this reason a yellow line-fillled area marking the 90% sample dispersion has been added to Fig. 2, both panels.
The GRB prompt γ-ray emission is known to span different orders of magnitude in luminosity. We correct for the different prompt luminosity characterising different events re-normalising the flaring component by the average 15-350 keV luminosity observed during the BAT T90 (interval of time of emission of 90% of the BAT fluence). The re-normalised flaring luminosity L * ≡ L/Lprompt in the time interval t * i − t * f is computed as: where L * i ≡ Li/ Li,γ , being Li,γ the average 15-350 keV prompt luminosity during the T90 of the i th burst. For each burst, the time t * is expressed in T90 units 4 . Both the BAT T90 and the fluence information are retrieved from the burst BAT-refined circulars 5 . Figure 2, panel (b), shows the resulting re-normalised average flaring luminosity curve: for t * 1 the XRT caught the X-ray counterpart of the prompt γ-ray emission; in the time interval 1 t * 100 the renormalised flaring activity follows a power-law decay with best-fitting index α = −1.2 ± 0.1.

FLARE DETECTION THRESHOLD
Flares are always detected as emission components superimposed on top of the overall afterglow decay: this means that whatever the origin of the continuum emission is, the sensitivity to temporal fluctuations is degraded as the flux level -and consequently the statistics-of the underlying continuum decays with time (see e.g. Morris 2008). This naturally requires flares at late times to be longer in duration and/or have a higher flux contrast with respect to the continuum in order to be detected as statistically significant fluctuations. It is therefore of primary importance to understand at which level the evolution with time of the average flaring component portrayed in Fig. 2 is a by-product of the decay of the underlying X-ray afterglow which acts as a time-variable detection threshold.
To this end a set of simulations is run for each GRB of our sample: the aim is to determine the flux and luminosity detection threshold of the flaring component at any time t given the observed X-ray afterglow statistics. The continuum contribution to each afterglow and associated uncertainty have been calculated in Sect. 3: we evaluate the continuum at the original light-curve central time bins and combine the systematic uncertainty coming from the continuum evaluation procedure with a statistical uncertainty equal to the relative error on the original luminosity bins of that particular burst. A grid of 545 flare peak time values t pk equally spaced in logarithmic units between 1 s and 10 7 s is generated together with a grid of 171 flare peak luminosity values in the interval 10 40 − 10 57 erg s −1 . For each t pk , we generate 171 flares with growing peak luminosity values: a Norris et al. (2005) profile specified in terms of amplitude A, width w and peak time t pk is assumed. The Norris 2005 profile is defined as f ≡ f (A, τ1, τ2, ts) where τ1 and τ2 are two flare shape parameters related to the rise and decay phases, respectively, while ts defines the pulse starting time. We take advantage of previous studies on the flare phenomenology to reduce the flare parameter space and write τ1 ∼ 4w/3; τ2 ∼ w/3; ts ∼ t pk − 2w/3; w ∼ 10 −0.94 t (1.12±0.18) pk (see C10 and M10 for details). The signal of the fake flare is then integrated in the original light-curve bins and a statistical error is assigned equal to the relative uncertainty affecting the original luminosity curve at the same time. Finally, the flare contribution is added to the continuum and the uncertainties propagated into a final fake GRB afterglow. To simulate the entire procedure, the fake afterglow is then processed by the automatic excess estimation routine described in the previous section: a new continuum is determined together with the associated flaring component, as before. If the newly determined flaring component contains a positive fluctuation with a minimum 2σ significance, the peak luminosity of the fake flare is recorded as 2σ flare detection threshold at t pk and the following peak time value is considered. The entire procedure is repeated for each t pk , every time starting from the minimum peak luminosity value of the grid. The result is shown in Fig. 1 for GRB 060607A taken as an example (blue solid line): due to the degrading quality of the signal, a growing flare-to-continuum flux contrast is required for the flare to be detected at later times.
The 2σ detection luminosity curve of the different GRBs are then combined to produce the average detection threshold L th (t) in the rest frame band 2.2-14.4 keV. In strict analogy with Eq. 1, L th (t) is defined as: being N the number of GRBs displaying a significant positive flaring component in the time period ti − t f ; L i,th is the luminosity of the 2σ detection threshold of the i th GRB at time t = (ti + t f )/2. The re-normalised average 2σ threshold L * th (t * ) is computed as: where L * i,th ≡ L i,th / Li,γ . As before, N collects the GRBs with significant excesses in t * i − t * f . The result is portrayed in Fig. 2: in both panels a thick green area marks the 2σ detection threshold region. From panel (a) it is clear that the −1.2 slope is mainly due to the selection of bright fluctuations able to overshine the threshold. The fact that L(t) tracks the L th (t) temporal behaviour suggests that the sample of significant excesses detected in the time period 2-3000 ks are not representative of the real population of flares at those times. The bright end of the flare luminosity distribution has been sampled and a biassed average flare luminosity has been re-constructed. Notably, the decay of the detection threshold at t ∼ 10 4 s allows us to appreciate a steeper decay (α ∼ −2.5) of 3σ significant positive fluctuations extending to t ∼ 6 × 10 4 s. No more than 2σ significant flaring activity can be quoted for (1 < t < 6) × 10 5 s due to the flatter decay (and actually enhancement) of L th (t) . The unbiassed average flare luminosity curve is therefore likely to be steeper than −1.2 at late times. At early times the situation is different: for t 300 s the L(t) ∝ t −2.7 decay is not affected by threshold effects; the extension of the same temporal law up to t ∼ 1000 s suggests that L th (t) is likely to play a minor role in the time interval (0.3 < t < 1) ks, as well. Finally, Fig. 2, panel (b) shows that the re-normalised average flare luminosity tracks the threshold starting from t * ∼ 100. L * (t * ) is consistent with the −1.2 power-law decay up to t * ∼ 800. For t * 1000 the re-normalised activity significance is 2σ.

LUMINOSITY-TIME PLANE: THE FLARE AVOIDANCE REGION
The 2σ detection threshold luminosity curves calculated in Sec. 4 define for each GRB, for each time t, the minimum peak luminosity which a flare must have to be detected as statistically significant positive fluctuation above the X-ray afterglow (see Fig. 1, blue solid line). We combine this information with the residual luminosity curves derived in Sec. 3 to draw an upper limit to the X-ray flaring luminosity associated to each afterglow at any time t, as follows: for each event we substitute the luminosity value of non significant fluctuations with the respective 2σ upper limits from the corresponding threshold luminosity curve, keeping the significant fluctuations unchanged. In this way we account for the presence of undetectable flaring activity. The resulting curve divides the luminosity-time plane into two halves: the real flaring luminosity curve is expected to lie below this line at a minimum 2σ c.l.: we refer to this region as the X-ray flare-area associated to a particular GRB. This procedure is repeated for each burst of the sample for both luminosity and re-normalised luminosity units. The different flare-areas are then superimposed and a luminosity-time plane contour plot is created (Fig. 3): the different colours refer to the different number of GRB flare-areas superimposed on a particular (ti, Li) pixel. In both panels the flaring luminosity of the 90% of the GRBs of our sample lies below the red solid line 6 . We refer to the area above the red line as the luminosity-time plane flare avoidance region. Note that a direct implication of Fig. 3, panel (a) is that in the time period 60 s t 400 s the flare luminosity function decay is steeper than ∝ t ∼−1.8 at 90% confidence.

DISCUSSION
In the previous sections we showed that: • For 15 s < t 30 s the average flare luminosity function is likely to track the prompt γ-ray emission; • In the time interval 30 s t 1000 s we find L ∝ t −2.7±0.1 . No threshold-related argument can be invoked to explain this scaling.
• For t 1000 s L closely follows the temporal behaviour of the threshold of detection: L ∝ t −1.2±0.1 . For this reason we believe the un-biassed average flare luminosity function at this epoch to be steeper.
Notably, a similar L pk − t pk scaling was obtained by Chincarini et al. (2010) for t < 1000 s once the intrinsic scatter of the relation was properly accounted for: L pk ∝ t −2.7±0.5 pk , where L pk and t pk are the flare peak luminosity 6 In the time interval 60 s < t < 2 × 10 5 s the 90% contour of Fig. 3 and flare peak time as obtained from a Norris et al. (2005) profile fit. In the following we therefore focus our attention on the 30 s < t < 1000 s interval of time, where the flare luminosity function L ∝ t −2.7 . In Sect. 6.1 we first investigate the source of the discrepancy between our results and those reported in L08. Possible physical interpretations in the context of both accretion and magnetar models of the central GRB engine are discussed in Sect. 6.1 and 6.2.
The GRB average X-ray flare luminosity function was first analysed by L08 starting from a sample of 24 flares coming from 9 long GRBs and 1 short GRB. The mean luminosity was found to decline with a much shallower power-law in time: L L08 ∝ t −1.5±0.2 . This scaling was obtained averaging the contribution of single flares of assumed square shape over time-scales longer than the flare duration. We investigate the source of the discrepancy ( L L08 ∝ t −1.5 vs. L ∝ t −2.7 of our work) below. We follow the prescriptions of L08 and apply their method to the sample of early time flares of C10 supplemented by late time flares with known redshift: the enlarged sample includes 55 flares detected in 29 GRBs. We do not assume a flare square shape and used instead the best fitting parameters of a Norris et al. (2005) profile: however, the choice of a particular flare profile is unlikely to bias the final result, as already noted by L08. The interval of times of integration of the mean luminosity function have been estimated from their Fig. 1. Fig-ure 4 shows the final result: the average flare luminosity function is best modelled by a power-law with best fitting index α = −2.6 ± 0.1. The shallower t −1.5 scaling is also represented for comparison. Errors are computed following the standard theory of error propagation. This result confirms the findings of Sect. 3 with a completely independent method. The reason for the discrepancy lies elsewhere.
We further check the consistency between the L08 method and ours, measuring the flaring activity evolution of GRB 050904. In this case we find L ∝ t −0.95±0.05 which is consistent with L ∝ t −1.0±0.3 quoted by L08. The flaring evolution of GRB 050904 is much flatter than the average behaviour ∝ t −2.7 . Prompted by this result we investigate the temporal evolution of the flaring component in the subsample of GRBs of Table 1 showing multiple episodes of flaring emission. The subsample contains 10 elements (black dots in Fig. 5): the mean flare function decay index is −1.8. This suggests that multiple flare GRBs have a flatter than average flare luminosity function. Interestingly, this subsample also shows a flatter than average steep decay: the median temporal decay index is found to be αsteep = 1.7 to be compared to αsteep = 2.5 of the remaining 34 GRBs of the original sample. This topic is further explored in Sec. 6.2. The average number of flares per GRB in the L08 sample is n f = 2.4 while our enlarged sample contains a wider population of single-flare GRBs: n f = 1.9. In particular, 60% of the GRBs of L08 contains more than one flare, while only 25% of the enlarged sample shows more than one episode of activity. The discussion above leads to the conclusion that the choice of a sample biassed towards multiple-flare GRBs (which was the wider sample of flares with known redshift available at the time of writing), led L08 to determine a flatter mean luminosity temporal scaling. In agreement with this picture, in the two GRBs with multiple flares for which L08 determined the flare function evolution they obtained α = −1.0 ± 0.4 (GRB 050803) and α = −1.0 ± 0.3 (GRB 050904): both values are flatter -even if consistentthan the average −1.5 ± 0.2 decay index derived from their entire sample.
While we benefit from a higher statistics which makes our sample more representative of the entire bursts population, the same choice of GRBs with easily recognisable flares in their afterglows is possibly introducing a bias which is difficult to quantify. The analysis of the flare avoidance region of Sect. 5, which properly accounts for the fraction of undetectable flares is the safest approach in this respect: this analysis indicates the average flare luminosity decay to be steeper than ∝ t −1.8 at the 90% c.l. for t 400 s. We finally investigate if the different prompt luminosity (and consequently afterglow luminosity) of the 44 GRBs of Table 1, could bias the flare L temporal scaling. Renormalising the flaring component of each GRB by the average prompt luminosity during the burst T90 duration 7 the L/Lprompt is best fitted by a power-law decay with index α ∼ −2.5 for t < 1000 s.
The ∝ t ∼−2.7 behaviour of the average flare luminosity has been obtained in 4 completely independent ways: (i) Equation 1 applied to 44 GRBs; 7 Note that this is different from what displayed in Fig. 2, where the x-axis is in T 90 units.
(ii) Equation 1 with re-normalised flaring contribution; (iii) L pk vs. t pk relation from C10; (iv) Application of the L08 method to the enlarged sample of 29 long GRBs with known redshift fitted with a Norris et al. (2005) profile.
It is important to underline that the ∝ t −2.7 behaviour represents the average scaling of the GRB flaring activity: the decay of the flaring component in a particular GRB can considerably differ from the average relation, as shown by Fig. 5 (see also Fig. 2, upper panel, where the yellow hatched area marks the 90% sample dispersion). Any model aiming at explaining the flare phenomenology is required to provide an explanation for the observed scatter, as well. Different physical interpretations are discussed below.

Accretion-powered flares scenario
In the case of a thin disk 8 with constant viscosity for t ≫ tα ν (where tα ν is the viscous time-scale) the accretion rate approaches the asymptotic regimeṁ ∝ t −1.25 (Franck, King & Raine 1985). A similar slope is obtained by Cannizzo et al. (1990) performing numerical simulations using the α-viscosity prescription of Shakura & Syunyaev 1973:ṁ ∝ t −1.2 . Recently Metzger et al. (2008) computed the viscous evolution of an efficiently ν−cooled, isolated ring of material under similar assumptions as the study of Cannizzo et al. (1990): they report a self-similar behaviouṙ m ∝ t −4/3 . The t −2.7±0.1 average flare luminosity power-law scaling obtained in this work is much steeper than the theoretical expectations listed above: an ad hoc non-linear relation between the luminosity and the outflow accretion rate, with the efficiency of conversion of mass accretion into luminosity decreasing as L/ṁ ∝ t −1.5 , is required to preserve the thin disk draining accretion scenario. Lee et al. (2009) demonstrated that the presence of powerful winds driven by the recombination of α−particles into nucleons and launched from the disk surface could lead to a substantial deviation from the t −4/3 regime, withṁ entering an exponential-like decay. This regime is unlikely to extend up to t ∼ 100−1000 s (see Lee et al. 2009, their Fig. 1). Alternatively, a temporal evolution of the flare opening angle causing later flares to be less beamed than earlier ones would restore the ∼ −1.5 slope 9 . However, most of the beaming evolution has been proven to take place at early times (see e.g. Morsony et al. 2007 and references therein) and is consequently unlikely to cause an overestimation of the flare luminosity steepness. We therefore consider alternative scenarios.
Fallback of the stellar envelope that did not reach the escape velocity during the explosion can provide the source of accreting material at late times: Chevalier (1989) showed that in this case the accretion rate is expected to decline asṁ ∝ t −5/3 which is too shallow to explain the observed t −2.7±0.1 flare luminosity scaling. The process was later investigated by MacFadyen et al. (2001): these authors simulated the fallback of the stellar envelope following the failure of the shockwave to unbind the star and recovered theṁRmin ∝ t −5/3 dependency (whereṁRmin is the radial fallback rate through their inner numerical boundary Rmin = 10 9 cm). However, departures from the t −5/3 scaling are expected if the infalling material passes through an accretion shock developing at R < 10 9 cm or/and if the mass fall-back rate decreases suddenly at some time t. These possibilities are discussed below. Kumar et al., (2008a,b) suggested that the prompt and X-ray afterglow reflect the modulation in the rate of central accretion of a rotating progenitor with a core-envelope structure onto a black hole. A thick disk configuration is assumed. According to their analytical treatment the conditionṁ ∝ t −3 can be basically achieved in two ways: first, in the case of rapid accretion, (tacc ≪ t),ṁ is expected to track the temporal evolution of the mass fall-back rateṁ fb . Assuming a 14 M ⊙ GRB progenitor star (model 16T1 of Woosley & Heger 2006), Kumar et al., (2008a,b) expectṁ fb ∝ t −3 (which directly translate intoṁ ∝ t −3 ) when the last ∼ 0.5 M ⊙ of the star near the surface is accreted. A second possibility arises when the stellar collapse leaves behind an accreting disk and no further mass is being added:ṁ fb = 0. The black hole accretion rate is found to scale as:ṁ ≈ṁ0t −4(s+1)/3 (Kumar et al. 2008b, their Eq. 43) 10 . The steepest decline allowed by these conditions iṡ m ∝ t −8/3 , suggestively close to the average flare luminosity t −2.7 temporal decay. This regime establishes for s = 1, when a convection-dominated accretion flow (CDAF, see e.g. Narayan et al. 2001 and references therein) sets in and the accretion rate decreases linearly with disk radius. In reality both the steeply falling density profile of the GRB progenitor and the loss of accreting matter in a wind (which is the signature of a CDAF or ADAF scenarios), are likely to play a major role in decreasing the fraction of mass which effectively makes it to the black hole. These two effect coupled together are able to account for Ljet ∝ t −3 or steeper, for s 0.5 (see Kumar et al. 2008b, their Fig. 4).
This scenario was further explored and confirmed by Lindner et al. (2010). Using axisymmetric hydrodynamical 2D simulations, Lindner et al. (2010) showed that a steady accretion rateṁ ∼ 0.1 − 0.2 M ⊙ s −1 giving rise to the prompt γ-ray phase is followed at t 40 s by a sudden and rapid power-law decline of the central accretion rate which is responsible for the prompt-to-steep decay X-ray afterglow transition. The steep decline is triggered by a rapid outward expansion of an accretion shock through the material feeding a convective thick accretion disk and lasts ∼ 500 s. During this phase the central engine is supposed to be active and the accretion rate is found to decay asṁ ∝ t −2.7 (Lindner et al. 2010). Again, we find the average flare luminosity function to show a very similar scaling (∝ t −2.7±0.1 for t < 1000 s), implying a linear relation between the luminosity and accretion rate: L ∼ ηaccṁ. The efficiency ηacc of mass rate to jet luminosity conversion depends on many parameters; among these, the black hole spin a is likely to play a major role. In Fig. 4 we use ηacc = 10 −3 which corresponds to a = 0.75  Table 1 with extended flaring emission (black dots). GRBs which show multiple X-ray flares but otherwise lack of the redshift measurement are indicated with red stars. These events are not part of the main sample of Table 1.
following the prescriptions of McKinney (2005). The main limitation of this study lies in that Lindner et al. (2010) do not simulate the fluid magnetic field which could alter thė m(t) dependence.

Magnetar models
Pulsars arise from core-collapse. A subset of them are inferred to host magnetic fields B ∼ 10 14 − 10 15 G. The rotational energy of a rapidly rotating magnetar can be extracted to power a strong relativistic jet (e.g. Bucciantini et al. 2009) that would give rise to the prompt emission. The energy release in radiation of this system at late times (100 s t 1000 s) is subject to great uncertainties (T. Thompson, private communication). However, the output of radiation is likely to be between two limiting cases: dipole radiation and monopole radiation. In the former case L ∝ (1+t/T ) −2 which implies L ∝ t −2 for t ≫ T , where T is the characteristic spin-down time-scale (e.g. Wheeler et al. 2000). In the latter, L ∝ e −2t (Thompson et al. 2004). Based on the braking indices of observed pulsars in the Galaxy, L ∝ t −2.3 is expected (Bucciantini et al. 2006). Additionally, the interaction with the stellar material can be responsible for shallower slopes detected in multiple-flare GRBs (Lazzati et al. 2008).

The variability issue
Flares are observed as episodic, large-scale amplitude variations in the light-curves with a typical ∆t/t ≈ 0.1 (Chincarini et al. 2010): this automatically raises the question of the physical source of the flare variability. The presence of a correlation between the properties of the prompt and flare emission within the same GRB could in principle shed light on the mechanism powering the flaring activity. However we found no correlation between the prompt luminosity and the total energy and luminosity of the flaring component; the same is true if one were to consider the 15-150 keV prompt energy; no correlation has been found be- tween the prompt luminosity (or energy) and the peak time of the last flare; the peak time of the last flare is also not correlated with the prompt T90; the number of prompt pulses is not the key factor determining the number of flares to appear at t 50 s (Chincarini et al. 2007); finally, a hint for a prompt-fluence vs. flare-fluence correlation was recently reported by C10: however, the weak correlation is mostly due to the presence of two short GRBs. The conclusion is that, while flares do share with prompt pulses several key observational properties (a notable example is the lag-luminosity relation of M10), at the moment it is not possible to infer the flaring activity of a burst starting from its prompt emission.
Flares are under-represented in simple power law X-ray afterglows (M10). The X-ray afterglow morphology -flaring activity connection is here further explored starting from the findings of Sec. 6.1: multiple flare GRBs have on average flatter flare luminosity functions (i.e. ∝ t −α with α < 2.7). This, together with the observation that the flare detection threshold of Fig. 2, upper panel, closely tracks the temporal decay of the average flaring activity for t 300 s, brought us to consider the possibility of a correlation between the X-ray flare luminosity decay and the underlying continuum emission decay. The average continuum is shown in Fig. 6: the best fitting power-law index αsteep = 2.8 ± 0.1 is consistent with the flare function decay ∝ t −2.7±0.1 . Suggestively, this happens for t 1000 s: at later times the continuum is likely to be dominated by the shallow decay component instead of the steep decay emission (Figure 6, lower panel, shows a decreasing flare-to-continuum ratio around t ∼ 400 s due to the progressively increasing contribution of the shallow decay component to the continuum). During the first ∼ 1000 s the typical flare-to-continuum ratio is L flare /Lsteep = 4.7 (median value) albeit with a large scatter. Using the data from C10 we obtain a similar median value F flare /Fsteep = 4.3 (Fig. 6, upper panel, inset). The physical mechanism powering each flare emission is therefore required to release an average amount of energy E flare : where ∆t is the flare duration. The median value ∆t/t = 0.23 from C10 has been used. The average continuum vs. flare function temporal behaviour of Fig. 6 results from the presence of a correlation linking the steep decay flux evolution to the flaring activity within individual GRBs. Modelling the steep decay and the flaring activity of each burst with decaying power-laws of index αsteep and α flare , respectively, we obtain the result drawn in Fig. 5 11 : flares seems to be linked to the contemporaneous steep decay flux evolution in a way that causes flatter flare luminosity functions to be associated to more gradual steep decays.
Generically speaking, these observations are consistent with a model where a first physical mechanism is responsible for the steep decay continuum (without flares), while mechanism 2 powers X-ray flares. The presence of the αsteep vs. α flare relation suggests that the two mechanisms are in some way related and that the sporadic appearance of mechanism 2 is triggered by some properties of mechanism 1. Instabilities affecting mechanism 1 can in principle provide the source of episodic releases of energy manifesting as flares. If this is the case, instabilities are likely to be triggered by some physical quantity related to the decay of continuum flux at time t: this would explain why late-time flares are so rare (only ∼ 5% GRBs show clear flaring activity for t > 1000 s), but also the paucity of flares in GRBs with simple power-law decay X-ray afterglows. In those GRBs, the steep decay is likely to be hidden by a contemporaneous but physically different emission component (see Margutti et al. 2010a for a detailed analysis of the two emission components in GRB 081028). We speculate that, given the remarkable similarity between X-ray flares and prompt pulses, it is possible that both mechanisms also operate during the γ-ray prompt emission producing slowly varying (mechanism 1) and fast varying (mechanism 2) emission components. The presence of both long (several seconds) and short (≈ ms) variability time-scales in the same burst is a known feature of the GRB prompt emission (e.g. Norris et al. 1996, but also Vetere et al. 2006, Borgonovo et al. 2007. L08 concluded that the outflow interaction with the stellar envelope cannot produce flares out of a continuous flow. In the following we analyse various mechanisms able to produce variability in the central engine release of en-ergy in the context of accretion (Sect. 6.2.1) and magnetar models (Sect. 6.2.2).

Variability in accretion models
In the context of GRB accretion models, a flare corresponds to a sudden increase in the jet luminosity as a consequence of an abrupt change of the mass accretion rate 12 . Generally speaking, largeṁ variations arise if the disk develops a ringlike structure: flares would be associated to the accretion of blobs of material initially located at various radii which subsequently evolve on the viscous time scale. This class of models would naturally account for the observed duration-time scale correlation and duration-peak luminosity anticorrelation (C10) as shown by Perna, Armitage & Zhang (2006). Different classes of disk instabilities could lead to this scenario 13 .
Viscous instability arises for dṁ/dΣ < 0 where Σ is the disk surface density (see e.g. Franck, King & Raine 1985). When this condition is satisfied more material will be fed in those regions of the disk that are more dense while material will be removed from less dense areas, so that the disc will likely break up into rings. Viscous and thermal instabilities (the latter taking place for dQ − /dT < dQ + /dT being Q + and Q − the heating and cooling rates respectively) have been invoked to explain Dwarf Novae outbursts: in this case the instabilities give rise to a limit-cycle behaviour (see e.g. Cannizzo et al. 1998) and do not lead to a total disk break down. While theṁ required to explain the Dwarf Novae outbursts luminosity is completely different from the conditions expected to hold in GRB hyper-accreting disks, we note that outside-in instability bursts (see Franck, King & Raine 1985 and references therein) qualitatively share some observational properties with GRB X-ray flares/prompt pulses: they have rapid rise -slower decay profiles and rise first at longer wavelength (see M10). The stability conditions for GRB hyper-accreting disks have been thoroughly studied: the steady state disk model of Di Matteo et al. (2002) led the authors to conclude that the accretion flow is both thermally and viscously stable for a variety ofṁ values supposed to give rise to the GRB phenomenology. However, X-ray flares suggest that the GRB engine is long -lived: a time-dependent computation is therefore required. The timedependent studies of Janiuk et al. (2004) and Janiuk et al. (2007) revealed the thermal-viscous instability to be an intrinsic property of the innermost disk radii for torus densities ∼ 10 12 g cm −3 ; the GRB disk is not stabilised but rather breaks down into rings leading to several episodes of dramatic accretion on the viscous time-scale of each ring. However, large hyper-accretion rates are required for the instability to set in:ṁ 10 M ⊙ s −1 . This value is ∼ 2 order of magnitude higher than the typical accretion rate invoked to explain the prompt GRB phase which iṡ mprompt ∼ 0.1 − 0.2 M ⊙ s −1 (see Lindner et al. 2010 and references therein). 12 The efficiency ηacc ≡ ηacc(a) of the process is unlikely to undergo abrupt temporal changes (see e.g. Kumar et al. 2008a,b). 13 The importance of hydrodynamical instabilities in collapsar disks has been recently demonstrated by the studies of Taylor et al. (2010).
Both the black hole spin and the presence of largestructure magnetic fields can deeply modify the accreting torus structure: for a rapidly spinning black hole, the diskblack hole magnetic coupling transfer rotational energy from the black-hole to the inner disk where viscous-thermal instability is known to arise for a Schwarzschild black hole. The magnetic coupling results in reducing theṁ value for which thermal and viscous instabilities establish: Lei et al. (2009) showed that the disk becomes unstable in its inner region for accretion rates as low asṁ > 0.1 M ⊙ s −1 ; Janiuk & Yuan (2010) found that for a Kerr black hole instability in the inner edge may be effective forṁ ∼ 0.5 M ⊙ s −1 , the actual value depending on the viscosity parameter α. While a time-dependent analysis is required also in this case, these works prove that viscous and thermal instabilities could indeed arise forṁ ≈ṁprompt.
Alternatively, the ring-like structure of the accretion disk could be the outcome of a gravitational instability occurring in the disk outskirts. Gravitational instability arises when Q = Ωcs/πGΣ < 1 (Toomre 1964) and is known to lead to fragmentation in the outer regions of active galactic nuclei (see e.g. Shlosman et al. 1990) and of young stellar objects (see e.g. Adams & Lin 1993). Hyper-accreting GRB disks in the steady state regime have been proven to undergo gravitational instabilities whenṁ 2 M ⊙ s −1 and radii r 30 rg (rg ≡ 2GM/c 2 ) under the assumption of a viscous parameter α = 0.1. The threshold to gravitational instability lowers toṁ 0.2 M ⊙ s −1 , r 15 rg for α = 0.01 (Chen & Beloborodov 2007, but see also Di Matteo et al. 2002). This mechanism has been suggested as possible source of the GRB X-ray flares by Perna, Armitage & Zhang (2006) (see however Proga & Zhang 2006;Piro & Pfhal 2007 for a critical view).
Whatever the origin of the disk fragmentation 14 and/or ring-like structure is (i.e. viscous, thermal or gravitational), the total energy released by each flare is likely to be proportional to the total mass of the fragment m f which gives rise to the burst of accretion. Writing E flare ∼ L flare ∆t, where ∆t is a measure of the flare duration, since L flare ∝ t −2.7 and ∆t ∝ t (C10), the average energy of a flare is found to scale as E flare ∝ t −1.7 for 30 s < t < 1000 s. In particular, using the best-fitting average flare luminosity calculated in Sec. 4 and the flare width vs. peak time relation (w ∼ 0.2 × t pk ) from C10 we find: where E iso,flare refers to the isotropic energy emitted in the 1−10000 keV energy band. Equation 6 is valid in the time interval 40 s < t < 1000 s. A multiplicative factor of 2 has been applied to account for the limited 2.2-14.4 keV X-ray energy band in which L flare has been computed 15 . Our observations would therefore imply the fragments mass m f (t) ∝ E flare (t) to scale as: Since the arrival time t of a blob of material initially at a distance R is likely to scale as a positive power of the radius t ∝ R β with β > 0 (β = 3/2 in the case of advection dominated flows, Perna, Armitage & Zhang 2006, their Eq. 2), Eq. 7 suggests that fragmentation taking place at larger radii and/or later times gives rise to lower mass objects. This is qualitatively consistent with Perna, Armitage & Zhang (2006), their Eq. 5: however, the actual mass distribution of the fragments strongly depends on the local disk properties and can be addressed only via numerical simulations (R. Perna, private communication). Multiple flares GRBs would correspond to a flatter m f vs. t (and m f vs. R) dependence as implied by the results reported at the beginning of Sect. 6.2. Yet another possibility giving rise to variability in accretion models is a modulation of the mass fall back rateṁ fb as proposed by Kumar et al. 2008a,b: while the underlying physical mechanism has still to be understood, this mechanism is potentially able to give rise to sharp variations oḟ m provided that tacc ≪ t (as shown by Kumar et al. 2008b, their Fig. 6 and 7) and naturally accounts for the ∝ t −2.7 behaviour as discussed in Sect. 6.1.1. A possible realisation of this process is given by gravitational fragmentation of the fall back material into several bound objects (Rosswog 2007).
Magnetic fields can suppress disk fragmentation (Banerjee & Pudritz 2006), give rise to magneto-rotational instabilities (MRI) and strongly modify the dynamics of accretion both in magnetars and black-hole systems (see e.g. Zhang & Dai 2009;Proga & Begelman 2003 and references therein). In particular, a variable output may result from a modulation of the accretion rate by the magnetic-barrier and gravity (Narayan et al. 2003). As mass is being accreted on to the black-hole, the magnetic flux is accumulated in the inner region, causing the accretion through the torus to be repeatedly stopped and then restarted. When accretion resumes, an X-ray flare is produced (Proga & Zhang 2006). According to Spruit & Uzdensky (2005) the accumulated magnetic field B can support the gas against gravity until the radial magnetic force Fm is of the order of a few percent the gravitational force Fg. The force balance yields a minimum B ∝ṁ 1/2 . Since the magnetic flux captured by accretion also depends onṁ, then it is reasonable to expect the systems characterised by a flatterṁ(t) decay (manifested as a flatter afterglow steep decay) to be able to meet the Fg ∼ Fm requirement more easily and repeatedly (i.e. a number of flares are expected). This scenario would therefore naturally account for the α flare vs. αsteep relation of Fig.  5.

Variability in magnetic models
This section concentrate on scenarios where the variability completely depends on magnetic effects. Giannios (2006) suggested that flares can be powered by magnetic reconnection triggered by the deceleration with the external medium. However, the α flare vs. αsteep relation (Fig. 5), would require the steep decay to be intimately connected to the deceleration phase of the fireball: Swift observations instead favour an "internal" origin of the steep decay, being the steep de- . Flare peak flux to continuum ratio for the sample of early-time flares of C10. The distribution can be fitted by the superposition of two gaussian profiles with the following best fitting parameters: x 1 = 0.50 ± 0.03, σ 1 = 0.52 ± 0.03; x 2 = 2.10 ± 0.08, σ 2 = 0.36 ± 0.09. cay either interpreted as the tail of the prompt emission (e.g. Willingale et al. 2010 for a recent study) or the result of the prolonged engine activity (e.g. Kumar 2008a,b).
Alternatively, in the context of magnetar models differentially rotating millisecond pulsars can provide an engine able to repeatedly store and release energy. Differential rotation causes toroidal magnetic fields to be repeatedly wound up to ∼ 10 17 G and then pushed to and through the pulsar surface by buoyant forces: this allows the neutron star spin energy to be emitted in powerful bursts of pulsar wind (Kluzniak & Ruderman 1998;Ruderman et al. 2000), an extreme but transient realisation of the Usov (1992) pulsar. This mechanism has been invoked to explain the flare phenomenology by Dai et al. (2006). The re-windup time, i.e., the time between subbursts is anticorrelated with the pulsar differential rotation ∆Ω; however, the luminosity of the subbursts (manifesting as flares) is expected to roughly scale as ≈ Ω(t) 4 , being Ω the neutron star rotation rate (Ruderman et al. 2000, their Eq. 19 and 20). If the pulsar rotation Ω is also the ultimate source of energy which powers the steep decay, then it is not unreasonable to expect the correlation displayed in Fig. 6 and Fig. 5.

SUMMARY AND CONCLUSIONS
We analysed the X-ray afterglows of 44 long GRBs observed by Swift with the goal to determine the average flaring component in the common rest frame energy band which is 2.2 − 14.4 keV. No assumption on the flare functional shape has been made. Our work highlights the importance of the proper consideration of the threshold of detection of flares against the contemporaneous continuous X-ray emission. In particular we showed that: • The best fit 2.2-14.4 keV average flare luminosity curve is: L = 10 (54.5±0.1) (t/s) −2.7±0.1 (erg s −1 ) for 30 s < t < 1000 s (Fig. 2). A similar scaling was obtained by C10 for the subsample of 43 early-time flares with redshift: L pk ∝ t −2.7±0.5 pk . The application of the L08 method to this subsample of burst leads to L ∝ t −2.6±0.1 (Fig. 4).
• For t > 1000 s threshold effects related to the presence of a continuum X-ray emission underlying the erratic appearance of flares start to play an important role. The resulting L ∝ t −1.2±0.1 is biassed towards the bright end of the flare luminosity distribution at these times. The unbiassed L is likely to be steeper (Fig. 2).
• According to the flare avoidance region analysis which properly accounts for the fraction of undetectable flares at any time t, the power-law decay index of L is steeper than −1.8 at the 90% c.l. for t < 1000 s (Fig. 3).
• GRBs with multiple-flare have a flatter than average flare luminosity function: L ∝ t −α with 0.6 α 3. Parenthetically, this is probably the reason why L08 determined a flatter average flare luminosity function L L08 ∝ t −1.5 .
• The decay of the continuum closely tracks the decay of the average flare luminosity function (Fig. 6) for t < 1000 s, suggesting that the two components are deeply related to one another. In particular, within individual GRBs, the power-law decay index of the steep decay is positively correlated to the power-law decay index of the flaring component (Fig. 5). GRB 100212A is a show case in this respect (Grupe et al. 2010). The typical flare to steep-decay luminosity ratio is L flare /Lsteep = 4.7.
• As a result, the typical flare energy at time t < 1000 s obeys the following relation: E flare ∼ L flare (t)∆t ≈ tLsteep(t) These findings suggest a model where the steep decay is produced by some form of activity of the internal engine which would be required to be still alive at those times (see however Genet & Granot 2009 for a complementary view), while flares could be powered by instabilities affecting the physical source of energy which gives origin to the steep decay. This would explain the X-ray flares erratic behaviour. In this picture, the shallow decay phase would be due to a completely distinct component of emission which progressively hides both the steep decay and the X-ray flares as time proceeds.
The L ∝ t −2.7±0.1 has been analysed in the context of accretion and magnetar models of GRBs. In particular: • According to the hyper-accreting black hole scenario, the L ∝ t −2.7 scaling can be obtained in the case of rapid accretion (tacc ≪ t) or when the last ∼ 0.5M ⊙ of the original 14M ⊙ progenitor star are accreted (Kumar et al., 2008a,b). Alternatively, the steepṁ ∝ t −2.7 behaviour could be triggered by a rapid outward expansion of an accretion shock in the material feeding a convective disk (Lindner et al. 2010).
• The rotational energy of a rapidly spinning magnetar can be extracted to power a jet whose luminosity is likely to be between the monopole (L ∝ e −2t ) and dipole (L ∝ t −2 ) cases. L ∝ t ∼−2.3 is expected, based on the braking indices of observed pulsars in the Galaxy.
In both scenarios the variability, which is the main signature of the flaring activity, establishes as a consequence of different kinds of instabilities. In the case of accretion models, thermal, viscous or gravitational instabilities could either lead to disk breakdown or fragmentation. Our analysis constrains the mass of the accreting material to scale as m f (t) ∝ t −1.7 (Eq. 7). However, the presence of magnetic fields gives rise to MRI instabilities and strongly modifies the dynamics of accretion: the accumulation of magnetic flux during the accretion can repeatedly stop and restart the accretion process (Proga & Zhang 2006). This would account for the erratic flare emission while explaining the α flare vs. αsteep relation. Alternatively, differentially rotating millisecond pulsars provide a viable mechanism where the existence of the α flare vs. αsteep relation can be reasonably explained. We note that if the flare origin is linked to the magnetic energy dissipation, the flare emission is likely to be polarised (Fan et al. 2005), while a disk fragmentation origin is likely to be accompanied by detectable gravitational wave signal (Piro & Pfhal 2007). Both signals will be detectable in the near future.
Whatever the mechanism powering the X-ray flare emission is, it is extremely difficult to account for the late-time (t > 1000 s) flare activity displayed by some bursts using the L ∝ t −2.7 component: exceptional circumstances leading to the revival of the instabilities would be needed. An interesting possibility is offered by Fig. 7: while the lower edge of the distribution of the flux ratio is probably incomplete, this figure suggests the existence of two populations of Xray flares (see also C10, their figure 13). A first population with flux contrast ∼ 5 (which is the one responsible for the average flare and continuum behaviour of Fig. 6) and a second population of bright flares with a typical ∆F/F ∼ 100 but extending up to ∆F/F ∼ 1000. Fig. 6 directly links the X-ray flares to the underlying steep decay flux. It is therefore possible that at late times only the small fraction of flares belonging to the second population are able to overshine the contemporaneous shallow decay component. This would explain why late time flares are so rare. The detailed characterisation of the two flare populations is beyond the scope of the present work and is left for a future study.