A common stochastic process rules gamma–ray burst prompt emission and X–ray flares
Abstract
Prompt –ray and early X–ray afterglow emission in gamma–ray bursts (GRBs) are characterized by a bursty behavior and are often interspersed with long quiescent times. There is compelling evidence that X–ray flares are linked to prompt –rays. However, the physical mechanism that leads to the complex temporal distribution of –ray pulses and X–ray flares is not understood. Here we show that the waiting time distribution (WTD) of pulses and flares exhibits a power–law tail extending over 4 decades with index and can be the manifestation of a common time–dependent Poisson process. This result is robust and is obtained on different catalogs. Surprisingly, GRBs with many () –ray pulses are very unlikely to be accompanied by X–ray flares after the end of the prompt emission ( Gaussian confidence). These results are consistent with a simple interpretation: an hyperaccreting disk breaks up into one or a few groups of fragments, each of which is independently accreted with the same probability per unit time. Prompt –rays and late X–ray flares are nothing but different fragments being accreted at the beginning and at the end, respectively, following the very same stochastic process and likely the same mechanism.
Subject headings:
gammaray: bursts, waiting time distribution1. Introduction
The first electromagnetic messenger of a gamma–ray burst (GRB) is the so–called –ray prompt emission, followed by the early X–ray afterglow on a timescale from minutes to hours. Long duration (– s) GRBs are nowadays known to be associated with the core collapse of some kind of massive stars rid of hydrogen envelopes (see Woosley & Bloom 2006; Hjorth & Bloom 2012 for reviews). Prompt –rays (with energies in the keV–MeV range) are observed within a given GRB as a sequence of pulses (typically a few up to several dozens). In addition, for a sizable fraction of GRBs the following decaying X–ray emission, which marks the end of the –rays, is characterized by the presence of X–ray flares which are sometimes observed as late as s (Burrows et al., 2005; Chincarini et al., 2007; Falcone et al., 2007; Curran et al., 2008; Bernardini et al., 2011). Although mounting evidence exists that X–ray flares, like –ray pulses, result from the GRB inner engine activity rather than from external shocks (Lazzati & Perna, 2007; Chincarini et al., 2010; Margutti et al., 2010), key questions remain unanswered: what radiation process(es)? What information on the inner engine can we extract? Is there a common process ruling inner engine activity across several decades in time?
As a matter of fact, both emissions represent a temporal point process, i.e. a time series characterized by the discrete occurrence of impulsive events superposed on a continuum. Intense bursting periods are often interspersed with relatively long (several up to tens of seconds) intervals with very low activity, compatible with the detector background, which are often referred to as quiescent times (QTs; RamirezRuiz & Merloni 2001; Nakar & Piran 2002; Quilligan et al. 2002; Drago & Pagliara 2007). The study of the waiting time distribution (WTD), i.e. of how time intervals between adjacent peaks distribute, provides clues on the nature of the stochastic process.
In particular, it reveals the degree of memory and correlation and constrains the physical process responsible for the discontinuous and bursty release of energy.
Processes showing similar on–off intermittency, or, equivalently, bursty behavior or clusterization, can be found in many fields (Platt et al., 1993). The corresponding WTDs often show power–law tails at long waiting times (WTs), whose index depends on the degree of clusterization of the time series. Examples encompass the aftershock sequence observed in earthquakes, described by Omori’s law (Utsu, 1961), neuronal firing activity, as well as a wide range of dynamical systems of human activity, such as mail and email exchanges (Oliveira & Barabási, 2005; Eckmann et al., 2004), phone calls (Karsai et al. 2012 and references therein) all the way to violent conflicts (Picoli et al., 2014). These processes are often modeled and interpreted in the context of self–organized criticality (SOC), where a nonlinear dynamical system reaches a stable critical point in which continuous energy input is released intermittently through avalanches and in a scale–free way. SOC naturally predicts power–laws in energy and WT distributions. See Aschwanden et al. (2014) for a recent review on the many areas displaying SOC behavior.
In astrophysics WTDs are studied for many different kinds of sources, such as outbursting magnetars (Göǧüş et al., 1999; Göǧüş et al., 2000; Gavriil et al., 2004), flare stars (Arzner & Güdel, 2004), and particularly the activity of the Sun throughout its cycle. WTDs of solar X–ray flares exhibit power–law tails with indices in the range – across several decades (Boffetta et al., 1999; Wheatland, 2000), depending on the class of flares and flux thresholds. Related bursty emission from the Sun such as coronal mass ejections (CMEs) are found to show very similar WTDs, whose index ranges from to in low to high–activity periods of the solar cycle (Wheatland, 2003). Likewise, WTDs of solar radio storms (Eastwood et al., 2010), of solar energetic particle and of solar electron events show very similar power–law indices (Li et al., 2014). Such power–law tailed WTDs are usually interpreted as the consequence of a time–varying Poisson process produced by SOC systems, in which the energy input rate is intermittent and directly affects the degree of clusterization of flares (Aschwanden & McTiernan, 2010; Li et al., 2014). In this model, the bursty energy release is the result of avalanches produced in active regions where magnetic flux is twisted by the moving footpoints, leading to a series of independent magnetic reconnection events and consequent plasma acceleration. Alternatively to SOC, interpretations in the context of fully developed MHD turbulence have also been proposed to explain the bursty dynamics and the power–law WTD: the intermittent character is the result of a nonlinear dynamics which makes the convective motion of the fluid and magnetic field swing between laminar and turbulent regimes repeatedly and chaotically (Boffetta et al., 1999; Lepreti et al., 2004).
The WTD between adjacent peaks in GRB –ray prompt emission profiles was found to be described by a lognormal –which implies some degree of memory– (Li & Fenimore, 1996) with an excess at long values due to QTs (Nakar & Piran, 2002; Quilligan et al., 2002; Drago & Pagliara, 2007). However, when the peak detection efficiency is carefully taken into account, it is found that the intrinsic WTD at short values is also compatible with an exponential, that is what is expected for a constant Poisson (thus memoryless) process (Baldeschi & Guidorzi, 2015). On the other side of the distribution long QTs could either mark the inner engine temporarily switching off, or result from modulation of the relativistic wind of shells (RamirezRuiz et al., 2001), or they could be due to a different physical mechanism from that of short WTs (e.g., Drago et al. 2008).
In spite of the impressive data that are routinely being acquired in the Swift era, little progress has been reported on WTDs in GRBs. Recently, energy and WT distributions for GRB X–ray flares have been shown to have power–law tails very similar to what is observed for solar X–ray flares. In particular, the WTD has a power–law index of (Wang & Dai, 2013). These results were interpreted as evidence for SOC possibly driven by magnetic reconnection episodes triggered in magnetized shells emitted by differentially rotating millisecond pulsars or, alternatively, by an hyperaccreting disk around a black hole (Popham et al., 1999).
Yet, there are several crucial issues which can be tackled with WTDs: is there additional evidence for a link between prompt –rays and late X–ray flares? To what extent do QTs differ from short WTs? Is it possible to provide a common description of short WTs, QTs, X–ray flares? What about rest–frame properties? Is there evidence for memory in GRB engines? What can be inferred on GRB engines through the WTD study?
In this paper we address these issues through the analysis of the WTD of GRB prompt peaks for three independent data sets: Swift/BAT, CGRO/BATSE, and Fermi/GBM. For the Swift GRBs which have also been promptly observed with XRT, we present a joint analysis of –ray peaks and X–ray flares merged together. Section 2 describes the data sample selection and how we modeled the WTDs. Here we deliberately did not consider the energy distribution of peaks and flares, because even though our peak search algorithm identified moderately overlapped pulses, estimating their energy would require specific assumptions on their temporal structure. We therefore decided to postpone it for future investigation. The results, their implications and interpretation are reported in Sections 3 and 4, respectively. Hereafter, uncertainties on best fit parameters are given at 90% confidence, unless stated otherwise.
2. Data analysis
We searched all long–duration –ray light curves with mepsa^{1}^{1}1http://www.fe.infn.it/u/guidorzi/new_guidorzi_files/code.html (Guidorzi, 2015, 2014), a peak search algorithm designed and calibrated to this goal. The advantage of mepsa compared with analogous algorithms such as the one by Li & Fenimore (1996) (LF) is twofold:

it has a lower false positive rate. This is particularly true for the time intervals in which the signal drops to background between two adjacent activity periods: in the best cases the LF false positive rate is 3–5 bin, while the mepsa one is 1–2 bin (Guidorzi, 2015);

it has a higher true positive rate, especially at low signal to noise ratios (4–5).
2.1. Sample selection
2.1.1 Swift/BAT data
We started with the GRBs detected by Swift/BAT in burst mode from January 2005 to September 2014, collecting 825 GRBs. We extracted the mask–weighted light curves in the 15–150 keV energy band with a uniform bin time of 64 ms following the standard procedure recommended by the BAT team^{2}^{2}2http://swift.gsfc.nasa.gov/analysis/threads/bat_threads.html and applied mepsa. We then imposed a minimum threshold of significance, which ensures a very low false positive rate ( bin; Guidorzi 2015) and selected the GRBs with at least two peaks. We then removed from our sample the short duration GRBs (both with and without extended emission) by crosschecking with the classification provided in the BAT catalog (Sakamoto et al., 2011), as they will the subject of future investigation. Since this catalog does not include GRBs from 2010, for these GRBs we used the values as published in the BAT refined GCN circulars regularly published by the BAT team and set a conservative minimum threshold of s. A couple of GRBs detected by BAT exhibited a very long duration which could not be covered entirely in burst mode: in one case we used the WIND/Konus light curve for GRB 091024 (Virgili et al., 2013), while in the case of GRB 130925A we used the peak times as they have been obtained by Evans et al. (2014) from the corresponding Konus light curve. Finally, we ended up with a sample of 418 long GRBs with at least two significant () peaks each, totaling 2000 peaks and 1582 WTs. Hereafter, we refer to this sample as the BAT set.
2.1.2 Cgro/BATSE data
We took the concatenated 64–ms burst data distributed by the BATSE team ^{3}^{3}3ftp://cossc.gsfc.nasa.gov/compton/data/batse/ascii_data/64ms/. For each curve we interpolated the background by fitting with polynomials of up to fourth degree as suggested by the BATSE team (e.g., Guidorzi 2005). Like in the BAT case, we applied mepsa to an initial sample of 2024 light curves in the full passband. We applied the same selection on the peak significance and selected the long GRBs by requiring s, where was taken from the GRB catalog^{4}^{4}4http://gammaray.msfc.nasa.gov/batse/grb/catalog/current/tables/duration_table.txt (Paciesas et al., 1999). We ended up with a sample of 1089 long GRBs with at least two – significant peaks. Overall we collected 7649 peaks and 6560 WTs. Hereafter, this will be referred to as the BATSE sample.
We also applied the same selection procedure to the light curves corresponding to the sum of the two softest energy channels (1 and 2) and to the sum of the two hardest channels (3 and 4), respectively within the 25–110 keV and 110 keV bands. We collected 1065 and 922 GRBs, with 5156 and 4912 WTs, respectively. These two groups will be hereafter referred to as BATSE12 and BATSE34 sets, respectively.
2.1.3 Fermi/GBM data
We selected 586 long GRBs detected with Fermi/GBM (Meegan et al., 2009) from July 2008 to December 2013. We extracted the light curves of the two brightest GBM units in the energy band 8–1000 keV with 64 ms resolution and subtracted the background through interpolation with a polynomial of up to third degree. We selected the long GRBs by imposing s, where was taken from the official catalog.^{5}^{5}5http://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermigbrst.html We restricted to time intervals whose median values range from to s with reference to the trigger time. This corresponds to the time interval covered by the time tagged event (TTE) data type in trigger mode (Paciesas et al., 2012; Gruber et al., 2014). Before s and after 600 s the time resolution is that of CTIME data, s. In most cases we did not consider intervals s, because interpolationestimated background often becomes critical and the required effort for a proper estimate is beyond our scope (Gruber et al., 2011; Fitzpatrick et al., 2012). We did not consider GRBs showing prolonged activity beyond this time interval. We then applied the same selection criteria as for the previous sets. Finally, by visual inspection we removed phosphorescence spikes due to high–energy particles (Meegan et al., 2009), by comparing the same profiles in different GBM units. We ended up with a final sample of 2383 peaks out of 544 GRBs with at least two significant peaks. The total number of WTs is 1839.
2.1.4 Swift/XRT X–ray flares
We considered the catalog of 498 X–ray flare candidates detected with Swift/XRT obtained by Swenson & Roming (2014). This was extracted from 680 XRT light curves from January 2005 to December 2012 with a method based on the identification of breakpoints in the residuals of the fitted piecewise power–law light curves: these points mark sudden changes in the mean value due to unfitted features. The optimal set of breakpoints was then found by minimizing the residual sum of squares against piecewise constant functions. To counter the effect of overfitting with unnecessary breakpoints, they made use of the Bayesian Information Criterion (see Swenson et al. 2013; Swenson & Roming 2014 for further details). In this catalog each candidate is assigned a confidence value. We conservatively selected the subsample with a minimum confidence of 90%, ending up with 205 X–ray flare candidates.
We separately merged each X–ray flare catalog with the BAT one by joining, for each GRB, the sequence of –ray peak times and flare peak times into a unique sequence of temporal peaks. In doing this, every peak which was seen in both instruments was tagged as a –ray peak and not considered any more as an X–ray flare. Analogously to the requirements for the previous sets, we selected those GRBs with at least two (either or X–ray) peaks, so as to have at least one WT. We ended up with a sample of 1098 (954 – and 144 X–ray) peaks in 244 GRBs (01/2005 – 12/2012). We hereafter refer to this joint set as BATX sample.
Finally, we selected the subsample with known redshift, so as to derive the WTD in the source rest–frame. This was done by simply correcting for cosmic dilation and thus dividing the observed WTs by the corresponding . Unlike the width of a given pulse, which is affected not only by cosmic dilation but also by the energy passband shift, for their nature WTs are affected by the latter to a much lesser extent. We found 359 WTs in 94 GRBs with known redshift. The subset with known redshift will be referred to as BATXz.
As an independent check, we in parallel considered the X–ray flare catalogs of Chincarini et al. (2010) and Bernardini et al. (2011), which respectively include 113 early–time ( s) X–ray flares from April 2005 to March 2008 and 36 late–time ( s) flares from April 2005 to December 2009. However, due to lower statistics, we hereafter focus on the BATX sample.
2.2. Waiting time distribution modeling
In physics a Poisson process is usually assumed to be characterized by a constant expected rate. The WTD of this process is exponential with e–folding ,
(1) 
where is the constant mean rate and is the mean WT. A time–varying Poisson process is characterized by a variable mean rate : the process is locally Poisson, but the expected rate changes with time as described by . According to this definition, the resulting process is the combination of two different processes at play and is often referred to as a “Cox process” (e.g., Cox & Isham 1980):
 (a)

at a given time events are generated according to a Poisson process with rate and, as such, are statistically independent;
 (b)

the expected rate is itself a function of time, which can vary either randomly or deterministically as time passes.
To derive the corresponding WTD, one may approximate as a piecewise constant function in a number of adjacent time intervals and treat it as a sequence of several Poisson processes with rate . Following Aschwanden & McTiernan (2010) and references therein, the resulting WTD is
(2) 
where
(3) 
is proportional to the expected number of WTs in interval where . In the continuous limit, Eq. (2) becomes
(4) 
where is the total duration. When is either unknown or hard to treat, it is possible to define such that , that is the fraction of time during which the expected rate lies within the range . Equation (4) becomes,
(5) 
We adopted the model for provided by Li et al. (2014) in their eq. (5), which has been proposed to fit the WTD obtained for solar X–ray flares and solar energetic particle events,
(6) 
with and free parameters and is a normalization term (). This model generalizes several other models which had been put forward in the same context (Wheatland, 2000; Aschwanden & McTiernan, 2010). The mean rate is
(7) 
where is the gamma function. From Equations (56) the corresponding WTD is
(8) 
and it is normalized like a probability density function (pdf), i.e. . There are only two free parameters, , which determines the level of clusterization, and the characteristic WT at which the WTD breaks: at , Eq. (8) becomes a power–law with an index of .
Equation (6) naturally gives rise to clusterization, i.e. time intervals characterized by an intense activity with a high rate of peaks (high ), interspersed with quiescent periods, during which the rate drops significantly (low ). The larger , the shallower the power–law regime at long WTs, and the more clustered the time profile (Aschwanden & McTiernan, 2010; Aschwanden et al., 2014). The details of how clustered the time profile looks like, in particular how the shot rate varies with time, are directly described by Eq. (6). At a given average rate , by increasing the variance of increases correspondingly, i.e. the shot rate varies more wildly. This means deviating more and more from the constant–rate case, thus enhancing the clustering character. Figure 1 illustrates the difference between a time–varying process like (, ) in Eq. (6) and a constant one sharing the same mean rate over a 100–s time window. The temporal sequence of events for the former is evidently more clustered than that of the latter and, in spite of the typical fluctuations of a Poisson point process, tracks the behavior of .
It is worth nothing that, in general, in a Poisson process individual events are independent of each other and, as such, have no memory of the events that occurred earlier, regardless of whether the expected rate is constant or time–dependent. The difference instead lies in the observed average rates as a function of time, so not on (a) but on (b): depending on whether varies either in a deterministic way, or randomly with/without memory, the average rate inherits the corresponding degree of correlation.
The WTDs we wanted to model are characterized by rare long WTs, where the count statistics is so low that one cannot use a simple minimization to fit the expected distribution of Eq. (8) to the counts collected in each bin. On the other side, merging the bins so as to have enough counts loses information and resolution. We therefore devised a log–likelihood based on Poisson statistics, which is essentially the C statistic (Cash, 1979) and holds exactly even in the low count regime. We used it in the context of a Bayesian Markov Chain Monte Carlo approach. The details are reported in Appendix A.
3. Results
Table 1 reports the best fit parameters for all of the WTDs we considered. In all cases the model of Eq. (8) provides an acceptable description. The lowest confidence level is that of BATSE (%), still comparable with nominal 5% usually adopted as a threshold. The BATSE sample is the largest (6560 WTs), so the high statistical sensitivity is likely to enhance small deviations from the model.
Figure 2 displays the WTDs for the –ray peak samples only. Apart from the GBM, whose power–law index is significantly steeper, the BAT and BATSE samples are fit with comparable indices, and , respectively. This is remarkable, given the different kind of detectors, energy passbands, different GRB populations each instrument is mostly sensitive to (Band, 2006). The soft and hard BATSE samples have the same index, showing no significant dependence on the energy channel.
Sample  Size  PL index  CL  

(s)  (%)  
BAT  1582  
BATSE  6560  
BATSE12  5156  
BATSE34  4912  
GBM  1839  
BATtrunc  1445  
BATX  854  
BATXz  359 
We investigated the reasons for the steeper WTD of the GBM set as follows: the dearth of long WTs is likely due to the shorter scanned time intervals, mostly from to s (Sect. 2.1.3). We therefore truncated the light curves of the Swift/BAT set and revised the WT selection accordingly. The results are reported in Table 1 as BATtrunc set, which includes 1445 WTs. Compared with the original BAT set, the WTD of the truncated data becomes steeper, from to , i.e. compatible with the GBM value within uncertainties. Hence we interpret the slightly steeper value of the GBM set as the result of shorter time profiles which disfavor long WTs.
Figure 3 displays the BATX set (squares) together with the corresponding best fit model. The power–law index is , i.e. compatible with BATSE sets within uncertainties. What is more, merging X–ray flares did not change the stochastic nature exhibited by the WTD, but extended its dynamical range by at least one order of magnitude with WTs as long as s. A common stochastic model is found to well describe the WTD observed across more than five orders of magnitude.
A similar result is obtained when one restricts to the known–redshift sample BATXz in the GRB rest–frame (circles in Fig. 3): here the power–law index is , i.e. somewhat shallower. The rest–frame characteristic time is significantly shorter because of cosmic dilation, s instead of the observer–frame values of 6–7 s.
We also searched for possible correlations between WTs and peak intensities and between WTs and peak fluences of adjacent pulses, but found none. Finally, we repeated the analysis for various subsets of GRBs, by requiring a minimum number of of peaks per GRB and found no significant difference.
3.1. –ray peaks vs. X–ray flares
Figure 4 shows the distribution of the number of –ray peaks per GRB for two different classes of GRBs, depending on whether their subsequent X–ray emission contains X–ray flares.
Surprisingly, it is found that almost all GRBs (23/25) with –ray peaks have no X–ray flares, although the two groups have comparable size, 131 and 163 GRBs with and without flares, respectively. The two distributions are unlikely to share a common population of events: a Kolmogorov–Smirnov (KS) test yields a mere % probability, i.e. they are different with (Gaussian) confidence. We visually inspected each of these – and X–ray light curves and found only one case of a flareless GRB, whose X–ray light curve exhibited some low–level flaring activity which did not pass the 90%–confidence threshold in the flare sample selection (Sect. 2.1.4). Therefore, GRBs with many pulses are unlikely to exhibit flares in the following declining X–ray emission.
4. Discussion
Our results may be summarized in four fundamental aspects:

–ray peaks and X–ray flares are compatible with being different aspects of the same stochastic process, which goes on after the end of the GRB itself and spans more than five orders of magnitude in time;

short (interpulse) and long (quiescent) WTs between –ray peaks are different realizations of the same stochastic process, the latters being only less frequent than the formers; hence a GRB with QTs and another without are by no means more different from each other than any other kind of GRBs are;

GRBs with several () –ray pulses are unlikely to exhibit X–ray flares after the end of the prompt emission;

–ray peaks and X–ray flares tend to cluster in much the same way that solar flares, energetic particle events, and CMEs do, even though the processes may be different.
The lognormal nature of the WTD originally claimed (Li & Fenimore, 1996) has recently been shown to be possibly an artifact of the peak detection algorithms in the short WT end ( s), where peaks significantly overlap and can hardly be separated (Baldeschi & Guidorzi, 2015). We found that the long value ( few s) tail needs no more to be described as the sum of a lognormal tail and a power–law excess due to the presence of QTs, that were interpreted as a different component (Nakar & Piran, 2002; Quilligan et al., 2002; Drago & Pagliara, 2007). This apparent diversity also concerns the so–called precursors (Lazzati, 2005; Burlon et al., 2008, 2009; Charisi et al., 2014), which are nothing but emission periods that are less intense than the following activity from which they are separated by a quiescent time. Our results (1. and 2.) show that all kinds of WTs, including precursors, can be described within a common stochastic process, and this holds all the way up to late X–ray flares, thus pointing towards a common mechanism, which keeps on working during and after the end of the prompt –ray emission, before the afterglow emission due to the interaction with the external medium takes over.
Another question concerns the break at low values in the WTD modeled in terms of the characteristic WT : is it an intrinsic property or is it entirely due to the low efficiency at short values of the peak detection algorithm? The capability of separating overlapping structures depends on a number of variables, such as the ratio between WT and peak widths, on intensities, and on temporal structures. While the drop at s is certainly due to the algorithm efficiency, the break itself modeled with is more complex: is shorter at harder energies (Table 1). A given pulse has a narrower temporal structure at harder energies (Fenimore et al., 1995), whereas in the softest energy channels there is a slow–varying component (Vetere et al., 2006). The presence of such soft component might hinder the peak identification in some case, so we examined the light curves in the harder channels. Visual inspection suggests that the paucity of subsecond WTs with respect to the power–law extrapolation is real and is unlikely to be a mere artifact of the peak identification process. In addition, minimum pulse widths observed in GRB profiles typically are in the range – s (Fenimore et al., 1995; Norris et al., 1996; Margutti et al., 2011). The MEPSA efficiency is above 10% for such pulse widths, for WTs s, and for measured signal to noise ratios (Guidorzi, 2015). It is therefore unlikely that the algorithm efficiency is entirely responsible for the observed exponential cutoff observed in the low end of the WTDs.
4.1. A simple toy model
We devised a very simple toy model to explore more in detail how a time–dependent Poisson process like the one of Eq. (6) could be obtained in a GRB engine. For the sake of clarity, suppose each pulse marks the accretion of a single fragment of an hyperaccreting disk. Actually, the idea behind this model is more general and only deals with the sequence of bursty emission episodes and their probability of occurring within a given time; however, hereafter we refer to the model of an hyperaccreting disk being fragmented as the source of the stochastic process which is responsible for the prompt –ray emission (Woosley, 1993; MacFadyen & Woosley, 1999) as well as for the subsequent X–ray flare activity (King et al., 2005; Perna et al., 2006; Kumar et al., 2008; Cannizzo & Gehrels, 2009; Geng et al., 2013). We briefly summarize the basic ingredients of the model, which are then thoroughly described in the remaining part of this Section:

a number of fragments are independently accreted with the same, constant, probability per unit time;

the number of available fragments is obviously decreasing with time; this naturally leads to a time–dependent Poisson process whose mean accretion rate decreases with time;

at the beginning, if the mean rate is too high (), accretion becomes inefficient and is suppressed by a factor of ;

for some (%) GRBs the reservoir of fragments is split into two separate groups sharing the same individual accretion probability per unit time, but with the second group becoming available only at later times (e.g., the late group could be identified with the outer part of the accretion disk).
Let us assume that the disk or the inner part of it has been split into a number of fragments, each of which has the same given probability of accreting per unit time, independently of the others. The probability for a given fragment to survive up to a given time is proportional to , where is the mean accretion time for each fragment. The total expected rate scales as the number of fragments that are still available, . The analogous to Eq. (6) is found as
(9) 
which corresponds to the case in Eq. (6) at .
Rather than a continuously varying, of this model is described more exactly by a piecewise Poisson process like the one of Eqs. (23), where is the expected rate when fragments are left over. All terms have equal weights ’s, since each piecewise constant process contributes one WT. The resulting WTD is thus given by Eq. (2),
(10) 
which can also be expressed as,
(11) 
where . In Figure 5 an example of such WTD is shown, with initial fragments, s. As time goes by, decreases and the e–folding of the individual exponential WTDs (thin solid) increase correspondingly. As a result, the total WTD (thick solid) show a power–law regime with index 2 at intermediate values of . At the total WTD is dominated by the initial exponential with e–folding , when fragments are all available. This agrees with the result that the intrinsic (i.e., corrected for the algorithm efficiency) WTD at short values is likely exponential, that is, compatible with a constant Poisson process (Baldeschi & Guidorzi, 2015).
Thus far, with reference to Eq. (6) our model implicitly assumed (see Eq. 9). However, in our attempt to reproduce the observed WTD with the piecewise constant process of Eq. (10) failed to model the observed break at . So we required that, when the expected rate becomes comparable or higher than , the number of observed WTs is suppressed by a factor of with respect to our model. This can be interpreted as if, when the number of fragments that can be readily accreted is such that the expected rate is , the overall process becomes inefficient and the rate is suppressed by . In other words, the number of WTs shorter than is smaller than what is expected from Eq. (9). This introduces some degree of memory in the initial stages of the accretion process: as long as the number of fragments ready to be accreted is too high () some of them are temporarily halted from accreting by some mechanism connected with the accretion rate itself. For instance, this self–regulating mechanism could be driven by the magnetic field (e.g., Proga & Zhang 2006; Uzdensky & MacFadyen 2006; Bernardini et al. 2013), which is known to have a complex role in accretion processes of utterly different objects such as T Tauri stars (Stephens et al., 2014). However, we cannot provide a more specific and physical justification for the exponential character of this self–quenching mechanism, which is therefore ad–hoc in its present formulation.
We assumed the logarithmic average and 1 scatter of the BATX WTD, s and a multiplicative scatter of , to generate the values for for each simulated burst. The number of generated peaks in each simulated curve was taken from the observed distribution and was augmented by 20% to ensure that the detected peaks were enough (since some are missed by the algorithm). The peak times for each simulated curve were randomly generated from an exponential distribution with e–folding , i.e. independently from each other. To mimic the drop in the peak detection efficiency at short WTs as well as the mechanism mentioned above about the suppression at high rates, we overlooked each peak occurring within s of the previous one, while through a binomial we assigned each a probability of being observed, where was set to the fitted value of the real WTD (Table 1).
To obtain a good match with the observed WTD over the same range we had to make a further assumption: we assumed that for a fraction of GRBs (%) randomly selected through a binomial, the disk is fragmented equally in two groups, the first of which is available for being accreted from the beginning (), while the second one becomes available from on, where is the common mean accretion time for each individual fragment from . The reason behind this is the observation of two similar bunches of –ray peaks interspersed with a long quiescent time (up to several tens of seconds) for a small fraction of GRBs. Physically, this could be the result of an outer part of the disk being accreted at later times with respect to the inner one or, more in general, delayed additional energy reservoir becoming available for late internal dissipation, with minimum variability timescale comparable with that of the early prompt emission (Fan & Wei, 2005; Lazzati & Perna, 2007; Troja et al., 2014). While the choice of the fraction of such GRBs and the duration of the quiescence period are somewhat arbitrary, the good match between simulated and real WTD does not depend crucially on them. Overall, the goal here is just to show the plausibility of the essential properties of this model, which can reproduce the observed properties in spite of the simple assumptions. We ended up with a set of 903 simulated WTs, whose distribution is compared with the real one in Fig. 6.
We further tested the consequences of this toy model by studying the distribution of the ratio between adjacent WTs. While WTDs describe how WTs distribute as a whole, losing information on their temporal sequence, the ratio distribution focuses on that. We therefore selected from the BATX as well as from the simulated sample the GRBs with peaks, so as to have at least two WTs, and derived the two distributions shown in Fig. 7. A KS test between the two sets yields 43% probability that they were drawn from a common population. Logarithmic mean and dispersion for the real (simulated) data are and ( and ). Similar results are obtained adopting other non–parametric tests, such as the Wilcoxon–Mann–Whitney or the more sensitive Epps–Singleton one (Epps & Singleton, 1986), respectively yielding 45% and 9% probability. Interestingly, simply replacing Eq. (9) with a constant Poisson process and applying the very same following steps, one ends up with a narrower and more zerocentered logarithmic ratio distribution, and , which is rejected with high confidence (p–value ) from a KS test. This means that for a constant process the ratio is, on average, closer to 1, and is less scattered around it than real data. The compatibility of the ratio distribution predicted by the toy model with the real one proves that the temporal sequence of WTs is compatible with an evolving Poisson process and is incompatible with a constant one on long timescales. In particular, X–ray flares are nothing but some of the last fragments that are left over and which are accreted on long timescales, when the rate decreases in a granular way, following the very same stochastic process ruling the accretion of the earliest ones. Hence, no correlation is to be expected between –ray prompt emission duration () and X–ray flares times, in agreement with observations (Liang et al., 2006). Finally, the result of Sect. 3.1 can be easily explained: GRBs with many –ray peaks accrete fragments rapidly with relatively short , so that at late time very few or none at all are left over for X–ray flares. The same result could be explained differently though: multi–peaked GRBs could have on average a brighter early X–ray afterglow continuum which overshines possible X–ray flares, which would then go unseen.
Overall, we assumed a direct connection between emission and observed times of the peaks. Within the context of internal shocks the observed time profile of both prompt –rays and late X–ray flares is strictly connected to the emission history (Kobayashi et al., 1997; Maxham & Zhang, 2009). Should this not be the case, little could be inferred about on emission times and potentially the times of individual accretion episodes from the study of the observed WTD. However, this connection becomes more complex due to the variety of Lorentz factors associated with the wind of shells colliding with each other. Even though the intrinsic duration of the GRB engine activity may differ by a factor of a few from the observed one (Gao & Mészáros, 2014), on average the temporal sequence of mutual collisions between randomly–assigned Lorentz–factor shells should track the emission time history. The nature of a given WTD is not altered as a whole when one passes from the emission to the observed times: in fact, each shell has a Lorentz factor which is in principle what can make the observed WTs very different from the emitted ones that is independent of the WTs preceding and following that shell. This statistical independence ensures that the observed WTD keeps memory on the emission time distribution. Only at late times, when the average Lorentz factor is expected to systematically decrease and the statistical independent character likely begins to fail, long WTs are likely to be affected as a consequence.
4.2. Solar activity: analogies and differences
It is remarkable and intriguing that WTDs of both solar eruptive events (X–ray flares, radio storms, high–energy particle events, CMEs) and of GRBs can be modeled with the same kind of timedependent Poisson process. The power–law characterization of the WTD heavy tail must not be overinterpreted from a mathematical viewpoint, since power–laws are, in general, what one ends up with when dealing with the sum of independent heavy–tailed variables. It works much in the same way that a normal distribution is the final outcome of the sum of independent finite–variance variables. In addition, claiming that data are power–law distributed is contrived whenever the explored range does not cover at least two decades (Stumpf & Porter, 2012). In this sense, invoking a SOC–driven mechanism for GRBs purely based on the power–law character of the WTD, and possibly of the energy distribution too, as suggested for X–ray flares from GRBs (Wang & Dai, 2013) or from other black hole systems (Wang et al., 2014), is a stretched interpretation of the data, as we explain below.
The same or very similar power–law indices imply that both processes have very similar degrees of clusterization, with analogous swings between intense and low–activity periods, apart from temporal rescaling (seconds for GRBs, hours for the Sun). However, one has to be careful extending this analogy to a common physical mechanism. Overall, there is a fundamental difference in terms of dynamical systems between GRB inner engines during core collapse and the Sun: for the latter, the regions where eruptive phenomena take place continuously receive energy, which is then released through avalanches, thus making the SOC interpretation plausible (although alternatives based on MHD turbulence seem equally compatible with observations). Instead, GRB engines are systems which start with a veryfarfromequilibrium configuration, evolve very fast using up all of the available energy, which no matter how much is limited. A GRB inner engine cannot return to its original configuration, it goes through an obviously irreversible evolution, whereas this is not the case for the solar active regions over sufficiently long timescales. For this reason, one needs not invoke SOC mechanisms related to accretion disks, in particular there is no need for a mechanism like the one proposed to explain fluctuations in black hole power spectra (Mineshige et al., 1994).
Therefore, a simple time–varying Poisson process explains the secular evolution of the mean rate of bursts/flares as well as the stochastic independent character of individual energy release episodes. This model disregards the physical origin of fragmentation and how energy is distributed among different fragments: thus, in principle it is compatible with various physical drivers, such as gravitational (Perna et al., 2006), or magnetorotational instability (Proga & Zhang, 2006).
5. Conclusions
In this paper we studied the waiting time distributions of GRB –ray pulses in three catalogs, CGRO/BATSE, Fermi/GBM, and Swift/BAT. For the latter, for the first time we merged –ray pulses and X–ray flares detected with Swift/XRT belonging to the same GRBs, and for a subsample the same analysis was carried out in the source rest–frame. We found that all WTDs can be described in terms of a common time–varying Poisson process that rules different waiting time intervals, which thus far in the literature have been treated differently: specifically, we showed that short WTs ( s), long quiescent times ( s) all the way up to late time X–ray flares are the manifestation of a common stochastic process. GRB WTDs exhibit heavy tails which are modeled with power–laws over 4–5 decades in time with indices in the range –, depending on the relative weight of late time events, such as X–ray flares, in each GRB sample. Because of the ubiquitous nature of power–laws (central limit theorem for heavy–tailed distributions), the character of WTDs must not be imbued with a mystical sense or overinterpreted as evidence for a universal process. In this sense, the similarity of the WTD power–law index with that of solar eruptive phenomena, such as flares and coronal mass ejections, proves nothing but a similar degree of clusterization in time. Nonetheless, it is remarkable that the WTD of –ray pulses and that of X–ray flares not only have compatible power–law indices but they join and extend the dynamical range for a common sample of GRBs. All this points to a common stochastic process ruling both phenomena. The unification under a common process of all different kinds of waiting times in GRBs (short interpulse times, long quiescent times, time intervals following precursors, time intervals between the end of the –ray prompt emission and subsequent X–ray flares) is a new result.
Another noteworthy result is that GRBs with many () –ray pulses are unlikely ( confidence in Gaussian units) to exhibit X–ray flares in their subsequent early X–ray emission. This result is naturally explained in the context of the time–varying Poisson process: many pulses observed in the prompt of a given GRB are indicative of a relatively short mean accretion time for a single disk fragment. Consequently, most of the available fragments are consumed during the prompt phase, with very few or none at all left over for the subsequent phase.
In the light of the irreversible evolution of GRB inner engines, the interpretation of a timevarying Poisson process appears to be simple and reasonable: the secular evolution of the expected rate of events is naturally linked to the energy reservoir being gradually used up, whereas the stochastically independent accretion of individual fragments is explained by the Poissonian character of the process.
Although self–organized criticality models naturally predict power–law tailed distributions of waiting times and energy, drawing upon this kind of dynamics for GRBs might be premature. Other equally plausible alternatives, such as fully developed MHD turbulence, can explain the same properties, as it was suggested for the solar case. Possible evidence for turbulence in GRBs has also been suggested from the analysis of power density spectra (Beloborodov et al., 1998, 2000; Guidorzi et al., 2012; Dichiara et al., 2013). Yet, we find that a simple time–varying Poisson process such as that of a system gradually using up all the available pieces already provides a remarkably accurate description.
The energy distribution, which was beyond the scope of this paper, will help to further constrain the stochastic process and possibly clarify whether more complex dynamical models, such as SOC or MHD turbulence, are to be considered.
Appendix A Log–likelihood to fit the distributions
Let the WTD consist of logarithmically spaced bins, each collecting counts. Let and be the lower and upper bounds of the –th bin (). Integrating Eq. (8) within this time interval yields the corresponding expected counts, , where we explicitly meant that it depends on the model parameters:
(A1) 
where , and is a function of both model parameters (Sect. 2.2). The probability of counts is ruled by the Poisson distribution where is the expected value,
(A2) 
where the dependence on model parameters is implicit through . The total probability is thus
(A3) 
where is the set of observed counts per bin (). The corresponding negative log–likelihood is therefore
(A4) 
We determine the best fit model parameters and their uncertainties in the Bayesian context. The posterior probability density function of the parameters for a given observed distribution , is (Bayes theorem)
(A5) 
where the first term in the numerator of the righthand side of eq. (A5) is the likelihood function of Eq. (A3), is the prior distribution of the model parameters, and the denominator is the normalization term. We assumed a uniform prior distribution, since no a priori information is available on the model parameters. The mode of the posterior probability of Eq. (A5) is therefore found by minimizing Eq. (A4).
We estimate the posterior density of the model parameters through a Markov Chain Monte Carlo (MCMC) algorithm such as the random–walk Metropolis–Hastings in the implementation of the ^{6}^{6}6http://cran.rproject.org/ package MHadaptive^{7}^{7}7 http://cran.rproject.org/web/packages/MHadaptive/index.html. (v.1.18). We initially approximate the posterior using a bivariate normal distribution centred on the mode and with covariance matrix obtained by minimizing Eq. (A4). For each WTD we generate sets of simulated model parameters and retain one every 5 MCMC iterations after excluding the first 1000. The remaining sets of parameters are therefore used to approximate the posterior density. Finally, once the best fit model parameters are determined, the bivariate posterior distribution of is sampled via MCMC simulations, which yield expected value and 90% confidence intervals for each of them.
As a matter of fact, since the bins in the low end of distribution are strongly affected by the poor efficiency of mepsa, these are to be ignored. In practice, one has to replace in Eq. (A1) with , where and are the first and last bins to be considered. In addition, the same in Eq. (A1) has to be further divided by a renormalizing factor so that it becomes,
(A6) 
For the WTDs discussed in the present paper we considered s ( s) in the observer (source) rest frame.
To assess the goodness of the fit for a given WTD, we use each set of simulated values for to generate as many synthetic WTDs from the the posterior predictive distribution. Hence for a given observed WTD, we directly calculate synthetic realizations of the same WTD. For each of these WTDs we calculate the negative log–likelihood with Eq. (A4) and derive a corresponding distribution of values, against which the value obtained from the real WTD is checked. This comparison directly provides a confidence level of modelling the observed WTD in terms of the best fit model of Eq. (8).
References
 Arzner & Güdel (2004) Arzner, K., & Güdel, M. 2004, ApJ, 602, 363
 Aschwanden et al. (2014) Aschwanden, M. J., et al. 2014, SSRv
 Aschwanden & McTiernan (2010) Aschwanden, M. J., & McTiernan, J. M. 2010, ApJ, 717, 683
 Baldeschi & Guidorzi (2015) Baldeschi, A., & Guidorzi, C. 2015, A&A, 573, L7
 Band (2006) Band, D. L. 2006, ApJ, 644, 378
 Beloborodov et al. (1998) Beloborodov, A. M., Stern, B. E., & Svensson, R. 1998, ApJ, 508, L25
 Beloborodov et al. (2000) Beloborodov, A. M., Stern, B. E., & Svensson, R. 2000, ApJ, 535, 158
 Bernardini et al. (2013) Bernardini, M. G., et al. 2013, ApJ, 775, 67
 Bernardini et al. (2011) Bernardini, M. G., Margutti, R., Chincarini, G., Guidorzi, C., & Mao, J. 2011, A&A, 526, A27
 Boffetta et al. (1999) Boffetta, G., Carbone, V., Giuliani, P., Veltri, P., & Vulpiani, A. 1999, Physical Review Letters, 83, 4662
 Burlon et al. (2009) Burlon, D., Ghirlanda, G., Ghisellini, G., Greiner, J., & Celotti, A. 2009, A&A, 505, 569
 Burlon et al. (2008) Burlon, D., Ghirlanda, G., Ghisellini, G., Lazzati, D., Nava, L., Nardini, M., & Celotti, A. 2008, ApJ, 685, L19
 Burrows et al. (2005) Burrows, D. N., et al. 2005, Science, 309, 1833
 Cannizzo & Gehrels (2009) Cannizzo, J. K., & Gehrels, N. 2009, ApJ, 700, 1047
 Cash (1979) Cash, W. 1979, ApJ, 228, 939
 Charisi et al. (2014) Charisi, M., Márka, S., & Bartos, I. 2014, in press, arXiv 1409.2491
 Chincarini et al. (2010) Chincarini, G., et al. 2010, MNRAS, 406, 2113
 Chincarini et al. (2007) Chincarini, G., et al. 2007, ApJ, 671, 1903
 Cox & Isham (1980) Cox, D. R., & Isham, V. 1980, Point Processes (London: Chapman & Hall)
 Curran et al. (2008) Curran, P. A., Starling, R. L. C., O’Brien, P. T., Godet, O., van der Horst, A. J., & Wijers, R. A. M. J. 2008, A&A, 487, 533
 Dichiara et al. (2013) Dichiara, S., Guidorzi, C., Amati, L., & Frontera, F. 2013, MNRAS, 431, 3608
 Drago & Pagliara (2007) Drago, A., & Pagliara, G. 2007, ApJ, 665, 1227
 Drago et al. (2008) Drago, A., Pagliara, G., & SchaffnerBielich, J. 2008, Journal of Physics G Nuclear Physics, 35, 014052
 Eastwood et al. (2010) Eastwood, J. P., Wheatland, M. S., Hudson, H. S., Krucker, S., Bale, S. D., Maksimovic, M., Goetz, K., & Bougeret, J.L. 2010, ApJ, 708, L95
 Eckmann et al. (2004) Eckmann, J.P., Moses, E., & Sergi, D. 2004, Proceedings of the National Academy of Science, 101, 14333
 Epps & Singleton (1986) Epps, T. W., & Singleton, K. J. 1986, J. Statist. Comput. Sim., 26, 177
 Evans et al. (2014) Evans, P. A., et al. 2014, MNRAS, 444, 250
 Falcone et al. (2007) Falcone, A. D., et al. 2007, ApJ, 671, 1921
 Fan & Wei (2005) Fan, Y. Z., & Wei, D. M. 2005, MNRAS, 364, L42
 Fenimore et al. (1995) Fenimore, E. E., in ’t Zand, J. J. M., Norris, J. P., Bonnell, J. T., & Nemiroff, R. J. 1995, ApJ, 448, L101
 Fitzpatrick et al. (2012) Fitzpatrick, G., McBreen, S., Connaughton, V., & Briggs, M. 2012, in Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, Vol. 8443, Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, 3
 Gao & Mészáros (2014) Gao, H., & Mészáros, P. 2014, arXiv:1411.2650, submitted
 Gavriil et al. (2004) Gavriil, F. P., Kaspi, V. M., & Woods, P. M. 2004, ApJ, 607, 959
 Geng et al. (2013) Geng, J. J., Wu, X. F., Huang, Y. F., & Yu, Y. B. 2013, ApJ, 779, 28
 Göǧüş et al. (1999) Göǧüş , E., Woods, P. M., Kouveliotou, C., van Paradijs, J., Briggs, M. S., Duncan, R. C., & Thompson, C. 1999, ApJ, 526, L93
 Göǧüş et al. (2000) Göǧüş, E., Woods, P. M., Kouveliotou, C., van Paradijs, J., Briggs, M. S., Duncan, R. C., & Thompson, C. 2000, ApJ, 532, L121
 Gruber et al. (2014) Gruber, D., et al. 2014, ApJS, 211, 12
 Gruber et al. (2011) Gruber, D., et al. 2011, A&A, 528, A15
 Guidorzi (2005) Guidorzi, C. 2005, MNRAS, 364, 163
 Guidorzi (2014) Guidorzi, C. 2014, MEPSA: Multiple Excess Peak Search Algorithm, ascl:1410.002
 Guidorzi (2015) Guidorzi, C. 2015, A&C, in press, arXiv:1501.01117
 Guidorzi et al. (2012) Guidorzi, C., Margutti, R., Amati, L., Campana, S., Orlandini, M., Romano, P., Stamatikos, M., & Tagliaferri, G. 2012, MNRAS, 422, 1785
 Hjorth & Bloom (2012) Hjorth, J., & Bloom, J. S. 2012, The GammaRay Burst  Supernova Connection 169
 Karsai et al. (2012) Karsai, M., Kaski, K., Barabási, A.L., & Kertész, J. 2012, Scientific Reports, 2, 397
 King et al. (2005) King, A., O’Brien, P. T., Goad, M. R., Osborne, J., Olsson, E., & Page, K. 2005, ApJ, 630, L113
 Kobayashi et al. (1997) Kobayashi, S., Piran, T., & Sari, R. 1997, ApJ, 490, 92
 Kumar et al. (2008) Kumar, P., Narayan, R., & Johnson, J. L. 2008, MNRAS, 388, 1729
 Lazzati (2005) Lazzati, D. 2005, MNRAS, 357, 722
 Lazzati & Perna (2007) Lazzati, D., & Perna, R. 2007, MNRAS, 375, L46
 Lepreti et al. (2004) Lepreti, F., Carbone, V., Giuliani, P., SorrisoValvo, L., & Veltri, P. 2004, P&SS, 52, 957
 Li et al. (2014) Li, C., Zhong, S. J., Wang, L., Su, W., & Fang, C. 2014, ApJ, 792, L26
 Li & Fenimore (1996) Li, H., & Fenimore, E. E. 1996, ApJ, 469, L115
 Liang et al. (2006) Liang, E. W., et al. 2006, ApJ, 646, 351
 MacFadyen & Woosley (1999) MacFadyen, A. I., & Woosley, S. E. 1999, ApJ, 524, 262
 Margutti et al. (2011) Margutti, R., Guidorzi, C., & Chincarini, G. 2011, International Journal of Modern Physics D, 20, 1969
 Margutti et al. (2010) Margutti, R., Guidorzi, C., Chincarini, G., Bernardini, M. G., Genet, F., Mao, J., & Pasotti, F. 2010, MNRAS, 406, 2149
 Maxham & Zhang (2009) Maxham, A., & Zhang, B. 2009, ApJ, 707, 1623
 Meegan et al. (2009) Meegan, C., et al. 2009, ApJ, 702, 791
 Mineshige et al. (1994) Mineshige, S., Ouchi, N. B., & Nishimori, H. 1994, PASJ, 46, 97
 Nakar & Piran (2002) Nakar, E., & Piran, T. 2002, MNRAS, 331, 40
 Norris et al. (1996) Norris, J. P., Nemiroff, R. J., Bonnell, J. T., Scargle, J. D., Kouveliotou, C., Paciesas, W. S., Meegan, C. A., & Fishman, G. J. 1996, ApJ, 459, 393
 Oliveira & Barabási (2005) Oliveira, J. G., & Barabási, A.L. 2005, Nature, 437, 1251
 Paciesas et al. (1999) Paciesas, W. S., et al. 1999, ApJS, 122, 465
 Paciesas et al. (2012) Paciesas, W. S., et al. 2012, ApJS, 199, 18
 Perna et al. (2006) Perna, R., Armitage, P. J., & Zhang, B. 2006, ApJ, 636, L29
 Picoli et al. (2014) Picoli, S., CastilloMussot, M. D., Ribeiro, H. V., Lenzi, E. K., & Mendes, R. S. 2014, Scientific Reports, 4, 4773
 Platt et al. (1993) Platt, N., Spiegel, E. A., & Tresser, C. 1993, Physical Review Letters, 70, 279
 Popham et al. (1999) Popham, R., Woosley, S. E., & Fryer, C. 1999, ApJ, 518, 356
 Proga & Zhang (2006) Proga, D., & Zhang, B. 2006, MNRAS, 370, L61
 Quilligan et al. (2002) Quilligan, F., McBreen, B., Hanlon, L., McBreen, S., Hurley, K. J., & Watson, D. 2002, A&A, 385, 377
 RamirezRuiz & Merloni (2001) RamirezRuiz, E., & Merloni, A. 2001, MNRAS, 320, L25
 RamirezRuiz et al. (2001) RamirezRuiz, E., Merloni, A., & Rees, M. J. 2001, MNRAS, 324, 1147
 Sakamoto et al. (2011) Sakamoto, T., et al. 2011, ApJS, 195, 2
 Stephens et al. (2014) Stephens, I. W., et al. 2014, Nature, 514, 597
 Stumpf & Porter (2012) Stumpf, M. P. H., & Porter, M. A. 2012, Science, 335, 665
 Swenson & Roming (2014) Swenson, C. A., & Roming, P. W. A. 2014, ApJ, 788, 30
 Swenson et al. (2013) Swenson, C. A., Roming, P. W. A., De Pasquale, M., & Oates, S. R. 2013, ApJ, 774, 2
 Troja et al. (2014) Troja, E., Piro, L., Vasileiou, V., Omodei, N., Burgess, J. M., Cutini, S., Connaughton, V., & McEnery, J. E. 2014, in press, arXiv 1411.1415
 Utsu (1961) Utsu, T. 1961, Geophysical Magazine, 30, 521
 Uzdensky & MacFadyen (2006) Uzdensky, D. A., & MacFadyen, A. I. 2006, ApJ, 647, 1192
 Vetere et al. (2006) Vetere, L., Massaro, E., Costa, E., Soffitta, P., & Ventura, G. 2006, A&A, 447, 499
 Virgili et al. (2013) Virgili, F. J., et al. 2013, ApJ, 778, 54
 Wang & Dai (2013) Wang, F. Y., & Dai, Z. G. 2013, Nature Physics, 9, 465
 Wang et al. (2014) Wang, F. Y., Dai, Z. G., Yi, S. X., & Xi, S. Q. 2014, ApJS, in press, arXiv 1411.4209
 Wheatland (2000) Wheatland, M. S. 2000, ApJ, 536, L109
 Wheatland (2003) Wheatland, M. S. 2003, Sol. Phys., 214, 361
 Woosley (1993) Woosley, S. E. 1993, ApJ, 405, 273
 Woosley & Bloom (2006) Woosley, S. E., & Bloom, J. S. 2006, ARAA, 44, 507