Nonstationary Fast-driven, Self-organized Criticality in Solar Flares

Waiting times Δt, the inter-event time intervals between two subsequent events of a Poissonian point process, are expected to exhibit an exponential function in the case of a stationary random process. The time series sample may consist of time intervals observed in statistically independent events and sampled at different locations and times. Thus, the probability distribution function p(Δt) is defined by

$$\begin{eqnarray}&&p({\rm{\Delta }}t)={\lambda }_{0}\ {\exp }^{-{\lambda }_{0}{\rm{\Delta }}t},\end{eqnarray} \tag{ 1 }$$

where λ₀ represents the mean event occurrence rate and the distribution is normalized to unity, i.e., ${\int }_{0}^{\infty }p({\rm{\Delta }}t)d{\rm{\Delta }}t=1$. A random process can be called a stationary Poisson process when the average flaring rate λ₀ is time independent and stays constant as a function of time.

A more general approach of waiting time distributions is the concept of inhomogeneous or nonstationary Poisson processes, where the mean flaring rate λ(t) becomes a function of time itself (e.g., Scargle 1998; Wheatland et al. 1998, 2000; Litvinenko & Wheatland 2001; Wheatland & Litvinenko 2002; Jaynes 2003; Sivia & Skilling 2006). Applying Bayesian statistics, a time series can be subdivided into Bayesian blocks, during which the occurrence rate λ_i is assumed to be piecewise stationary during a time interval [t_i, t_{i + 1}],

$$\begin{eqnarray}p(t,{\rm{\Delta }}t)=\left\{\begin{array}{cc}{\lambda }_{1}\ {\exp }^{-{\lambda }_{1}{\rm{\Delta }}t} & \mathrm{for}\ {t}_{1}\lt t\lt {t}_{2}\\ {\lambda }_{2}\ {\exp }^{-{\lambda }_{2}{\rm{\Delta }}t} & \mathrm{for}\ {t}_{2}\lt t\lt {t}_{3}\\ \ \ ....... & \\ {\lambda }_{n}{\exp }^{-{\lambda }_{n}{\rm{\Delta }}t} & \mathrm{for}\ {t}_{n}\lt t\lt {t}_{n+1}\end{array}\right.\end{eqnarray} \tag{ 2 }$$

The summation of the piecewise Bayesian blocks over discrete time intervals can be converted into a continuous integral function,

$$\begin{eqnarray}&&p({\rm{\Delta }}t)={\int }_{0}^{T}\lambda (t)\ p(t,{\rm{\Delta }}t)\ {dt},\end{eqnarray} \tag{ 3 }$$

where the probability p(t, Δt) in each Bayesian block is weighted by the number of events λ(t). The total duration of the time series is T, and the normalization is given by the total number of events, i.e., $N={\int }_{0}^{T}\lambda (t)\ {dt}$. Inserting the time-dependent probability $p(t,{\rm{\Delta }}t)=\lambda (t)\ {\exp }^{-\lambda (t){\rm{\Delta }}t)}$ into Equation (3) yields

$$\begin{eqnarray}&&p({\rm{\Delta }}t)=\displaystyle \frac{{\int }_{0}^{T}\lambda {\left(t\right)}^{2}\ {\exp }^{-\lambda (t){\rm{\Delta }}t}{dt}}{{\int }_{0}^{T}\lambda (t){dt}}.\end{eqnarray} \tag{ 4 }$$

Following Wheatland et al. (1998, 2000), we substitute the time variable t with the event occurrence rate λ, by defining the function f(λ) = (1/T) dt(λ)/dλ, which is equivalent to f(λ) dλ = dt/T,

$$\begin{eqnarray}&&p({\rm{\Delta }}t)=\displaystyle \frac{{\int }_{0}^{\infty }f(\lambda )\ {\lambda }^{2}\ {\exp }^{-\lambda {\rm{\Delta }}t}d\lambda }{{\int }_{0}^{\infty }\lambda f(\lambda )d\lambda }.\end{eqnarray} \tag{ 5 }$$

We now make a special choice for the flaring rate distribution f(λ) that contains (i) a reciprocal relationship f(λ) ∝ λ⁻¹ for small flaring rates λ ≲ λ₀ and (ii) contains an exponential drop-off at large flaring rates λ ≳ λ₀ (see also Equation 5.2.16 in Aschwanden 2011a),

$$\begin{eqnarray}&&f(\lambda )={\lambda }^{-1}\ \exp \left(-\displaystyle \frac{\lambda }{{\lambda }_{0}}\right).\end{eqnarray} \tag{ 6 }$$

The scale-free range of λ < λ₀ is visualized in Figure 2 (left panel, solid line), together with the exponential component (Figure 2, left panel, dashed line). The scale-free property with the scaling f(λ) ∝ λ⁻¹ is easy to understand, because the number of events f(λ) is proportional to the mean waiting time $\langle {\rm{\Delta }}t\rangle$, which in turn is reciprocal to the mean flaring rate $\langle \lambda \rangle$, e.g., $f(\lambda )\propto \langle {\lambda }^{-1}\rangle \propto \langle {\rm{\Delta }}t\rangle$, and thus is universally valid for every waiting time distribution. In addition, the exponential term in Equation (6) essentially produces an upper boundary of the reciprocal function at λ ≳ λ₀. The expression given in Equation (6) also fulfills the normalization ${\int }_{0}^{\infty }\lambda \ f(\lambda )\ d\lambda ={\lambda }_{0}$. The waiting time distribution (Equation (5)) can then be written as

$$\begin{eqnarray}&&p({\rm{\Delta }}t)={\int }_{0}^{\infty }\left(\displaystyle \frac{\lambda }{{\lambda }_{0}}\right)\exp \left(-\displaystyle \frac{\lambda }{{\lambda }_{0}}[1+{\lambda }_{0}{\rm{\Delta }}t]\right)\ d\lambda ,\end{eqnarray} \tag{ 7 }$$

which, by defining a = −(1 + λ₀Δt)/λ₀, corresponds to the integral $\int {{xe}}^{{ax}}{dx}=({e}^{{ax}}/{a}^{2})({ax}-1)$ and becomes ${\int }_{0}^{\infty }{{xe}}^{{ax}}{dx}=1/{a}^{2}$ when integrated over $[0\lt x\lt \infty ]$, yielding the solution $p({\rm{\Delta }}t)=1/({a}^{2}{\lambda }_{0}^{2}$), and we obtain for the waiting time distribution,

$$\begin{eqnarray}&&p({\rm{\Delta }}t)=\displaystyle \frac{{\lambda }_{0}}{{\left(1+{\lambda }_{0}{\rm{\Delta }}t\right)}^{2}}.\end{eqnarray} \tag{ 8 }$$

Note that this waiting time distribution contains no free variables, except for the normalization constant λ₀. Thus, this model predicts universally a power-law slope of α_Δt = 2 for any waiting time distribution. The only underlying assumption is the reciprocality of flaring rates and waiting times, which naturally emerges from the property of scale-freeness in SOC models (Aschwanden & McTiernan 2010).

Zoom In Zoom Out Reset image size

A comparison of stationary and nonstationary waiting time distributions is shown in Figure 3, as well as for a slow-driven and fast-driven SOC model (Figure 3). Note the reciprocal relationship between the flare occurrence rate (y-axis in Figure 3) and the waiting time (x-axis in Figure 3), differing by a factor of 10². A parametric set of theoretically predicted waiting times with various values of ${\lambda }_{0}=0.02,...,0.12$ is shown in Figure 4(c).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is customary to perform Monte Carlo simulations of waiting time distributions N(Δt)dΔt or occurrence frequency distributions N(x)dx by random generator values x = x₁, x₂, ..., x_n that have a prescribed function of their frequency distribution. Examples for exponential and power-law distributions are given in Section 7.1.4 of Aschwanden (2011a). The normalization is given by the integral of the probability function p(x),

$$\begin{eqnarray}&&{\int }_{0}^{\infty }p(x)\ {dx}=1.\end{eqnarray} \tag{ 9 }$$

The total probability ρ(x) to have a value in the range of [0, x] is then the integral

$$\begin{eqnarray}&&\rho (x)={\int }_{0}^{x}p(x^{\prime} )\ {dx}^{\prime} .\end{eqnarray} \tag{ 10 }$$

Then, we invert the integral function ρ(x) and denote it by the analytical inverse function ρ⁻¹, so that

$$\begin{eqnarray}&&x={\rho }^{-1}(\rho )={\rho }^{-1}(\rho [x]),\end{eqnarray} \tag{ 11 }$$

yielding a transform that allows us to generate values x_i from a distribution of probability values ρ_i. There are many numerical random generator algorithms available that produce a random number ρ_i in a homogeneous range of [0, 1], which can then be used to generate values x_i with the mapping transform ${x}_{i}={\rho }^{-1}({\rho }_{i})$. The frequency distribution of these values x_i will then fulfill the prescribed function p(x).

In our case, we want to simulate the waiting time distribution function that is given by the probability function p(Δt) (Equation (8)),

$$\begin{eqnarray}&&p({\rm{\Delta }}t)=\displaystyle \frac{{\lambda }_{0}}{{\left(1+{\lambda }_{0}{\rm{\Delta }}t\right)}^{2}},\end{eqnarray} \tag{ 12 }$$

which fulfills the normalization

$$\begin{eqnarray}&&{\int }_{0}^{1}p({\rm{\Delta }}t)\ d{\rm{\Delta }}t=1.\end{eqnarray} \tag{ 13 }$$

The integral function ρ(Δt) of the probability function p(Δt) is then

$$\begin{eqnarray}&&\rho ({\rm{\Delta }}t)={\int }_{0}^{{\rm{\Delta }}t}\displaystyle \frac{{\lambda }_{0}}{{\left(1+{\lambda }_{0}{\rm{\Delta }}t^{\prime} \right)}^{2}}d{\rm{\Delta }}t^{\prime} =\displaystyle \frac{{\lambda }_{0}{\rm{\Delta }}t}{1+{\lambda }_{0}{\rm{\Delta }}t}.\end{eqnarray} \tag{ 14 }$$

The inversion of the probability function ρ(Δt) is then simply

$$\begin{eqnarray}&&{\rm{\Delta }}t=\displaystyle \frac{\rho }{{\lambda }_{0}(1-\rho )},\end{eqnarray} \tag{ 15 }$$

which can be used to simulate a set of waiting times Δt_i using random numbers ρ_i in the homogeneous range [0, 1],

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{i}=\displaystyle \frac{{\rho }_{i}}{{\lambda }_{0}(1-{\rho }_{i})},\quad \mathrm{for}\ [0\lt {\rho }_{i}\lt 1].\end{eqnarray} \tag{ 16 }$$

Such a simulation for N = 575 events and λ₀ = 1.7 is shown in Figure 4 (middle panel), along with the theoretical distribution function p(Δt) (Equation (8)). An example of an observed waiting distribution function is shown in Figure 4(a), sampled from Geostationary Operational Environmental Satellite (GOES) M- and X-class flares (histogram in Figure 4(a)), which matches the predicted distribution (solid curve in Figure 4(a)) according to Equation (8) within the statistical uncertainties. A parametric set of the same type of waiting time distributions (Equation (8)) is shown in Figure 4(c), where the normalization constant λ₀ is varied.

We analyzed the same data set of 170 solar flares presented in Aschwanden et al. (2014a), which includes all M- and X-class flares observed with the SDO (Pesnell et al. 2011) during the first 3.5 yr of the mission. This selection of events has a heliographic longitude range of [−45°, +45°], for which magnetic field modeling can be facilitated without too severe foreshortening effects near the solar limb. We use the 45 s line-of-sight magnetograms from the HMI/SDO and make use of all coronal extreme-ultraviolet (EUV) channels of the Atmospheric Imaging Assembly (AIA)/SDO (in the six wavelengths 94, 131, 171, 193, 211, 335 Å), which are sensitive to strong iron lines (Fe viii, ix, xii, xiv, xvi, xviii, xxi, xxiv) in the temperature range of T ≈ 0.6–16 MK. For most of the analysis, we analyzed images with a cadence of 6 minutes, but present one event with the full AIA time cadence of 12 s.

The coronal magnetic field is modeled by using the line-of-sight magnetogram B_z(x, y) from HMI/SDO and (automatically detected) projected loop coordinates [x(s), y(s)] in each EUV wavelength of AIA. A full 3D magnetic field model B(x, y, z) is computed for each time interval and flare with a cadence of 6 minutes (0.1 hr). The total duration of a flare is defined by the GOES flare start and end times, including a margin of 0.5 hr before and after each flare. Such a margin (commensurable with the duration of the shortest flares; see Figure 8(d)) allows us to study the time evolution of a flare more comprehensively, because flare-related emission in soft X-rays often precedes the NOAA flare start time and extends past the NOAA flare end time.

The magnetic field is computed with the vertical-current approximation nonlinear force-free field (VCA-NLFFF) code, which is described for the original first version (Aschwanden 2013) and has been improved in accuracy in the second (Aschwanden et al. 2016b) and third (VCA3-NLFFF) versions (Aschwanden 2019b). Traditional nonlinear force-free field (NLFFF) codes have been found to produce large uncertainties in the horizontal (transverse to the line-of-sight) magnetic field components (e.g., Wheatland et al. 2000; Wiegelmann 2004; Wiegelmann et al. 2006, 2012), mostly due to the fact that the force-free magnetic field is extrapolated from photospheric magnetograms, although the photosphere is not force free (Metcalf et al. 1995). Improvements have been attempted by preprocessing of the line-of-sight magnetograms by additional constraints that minimize the force-freeness and net torque balance (Wiegelmann et al. 2012), by applying a magnetohydrostatic model (Zhu et al. 2013; Wiegelmann et al. 2017; Zhu & Wiegelmann 2018) or different magnetic helicity computation methods (Thalmann et al. 2019). In contrast, our method of the Vertical-Current Approximation (VCA) NLFFF code circumvents this problem, because the magnetic modeling is constrained by coronal loop structures, which are believed to be force free in the low plasma-beta corona (although not during flares, but before and after flares).

The main physical parameter that we are interested in here is the time evolution of the free energy, which is defined as the difference between the potential and nonpotential magnetic field, i.e., ${E}_{\mathrm{free}}(t)={E}_{{np}}(t)-{E}_{p}(t)$.

We show the time evolution of the free energy E_free(t) for 20 flare events (out of the 170 analyzed events) in Figures 5–7. We decompose the time profiles into pulses that consist of a rise time phase τ_rise = t_p − t_s and a decay time phase τ_decay = t_e − t_p. The peak times t_p are measured at the local maxima of the time evolution function E_free(t), and the starting times t_s and end times are derived from the local minima preceding and following each peak time. For clarity, we represent the decay phases of the pulses with gray areas in Figures 5 and 6. In Figure 5, we show relatively simple flare events with one single peak (n_p = 1) or two peaks (n_p = 2), while the 10 cases shown in Figure 6 were selected from the flare events with the longest duration, which exhibit from n_p = 5 to n_p = 13 peaks.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 7, we show the time profiles of the free energy E_free(t) with higher time resolution: the nominal resolution is 6 minutes (Figure 7(b)), an intermediate resolution is 1 minute (Figure 7(c)), and the full time resolution of AIA is 12 s (Figure 7(d)). The fluctuations visible at the highest cadence (Figure 7(d)) show a mean and standard deviation of E_free = (42 ± 7) × 10³⁰ erg, which indicates uncertainties of σ_E ≈ 7/42 ≈ 0.17. This uncertainty in the free energy includes numerical noise, mostly caused by the decomposition of unipolar magnetic charges from the HMI magnetograms and from the automated detection of coronal loops in the AIA images. Nevertheless, the time profiles shown in Figure 7 reveal about one to three significant pulses for this event, while Figure 6 shows 5–13 significant energy dissipation pulses per flare.

Zoom In Zoom Out Reset image size

Statistics of timescales is given in Figure 8. The number of energy dissipation pulses per flare ranges from n_p = 1 to n_p = 13 (see Figures 5 and 6), as derived from the (slightly smoothed) time profiles of the free energy, E_free(t). Each of the pulses is characterized by the rise time (which can be interpreted as magnetic energy loading time by new flux emergence), τ_rise = 0.1–1.2 hr = 6–72 minutes (Figure 8(a)); the pulse decay time (which can be interpreted as magnetic energy dissipation time), τ_decay = 0.1–0.7 hr = 6–42 minutes (Figure 8(b)); and the total pulse duration τ_pulse = 0.2– 1.5 hr = 12–90 minutes (Figure 8(c)). The lower limit of τ_rise,min = τ_decay,min = 0.1 hr = 6 minutes is caused by the chosen cadence in the calculation of magnetic energies.

Zoom In Zoom Out Reset image size

The flare duration times have a range of τ_flare = 1.1–5.2 hr (Figure 8(d)), which is about an order of magnitude longer than the pulse rise or decay times. This difference can be explained by the fact that the timescale of magnetic energy dissipation, which is similar to the duration of hard X-ray emission, is generally shorter than the timescale of soft X-ray emission, which was used by NOAA to define the flare duration.

Finally, we also measure the waiting times of flare events, using all GOES M- and X-class flare events during the first 3.5 yr of the SDO mission (from 2010 June 12 to 2014 November 16), including those events near the limb for which magnetic modeling was not feasible. The range of waiting times derived from the starting time difference of these 575 flares covers Δt = 0.2–2000 hr (Figure 8(e)). Note that truncation effects due to solar rotation and the selected longitudinal range (±45°) are ignored in the waiting time statistics here, although it could affect the correct waiting time measurement for events near the east or west limb. The waiting time distribution forms a power-law distribution with a slope of α_Δt = 2.0 for timescales of Δt ≳ 1 hr (Figure 8(e)) and closely follows the predicted function derived theoretically (Equation (8)) for a normalization constant of λ₀ = 0.07 hr⁻¹. According to the definitions of waiting times Δt, energy storage times τ_storage ≈ τ_rise, and energy dissipation times τ_diss ≈ τ_decay given in Figure 1, most of the storage times (Figure 8(a)) are much shorter than the waiting times (Figure 8(e)) and thus are consistent with the fast-driven SOC model (Figure 1(b)), rather than with the slow-driven SOC model (Figure 1(a)).

If the magnetic free energy is the main energy input in solar flares, and the energy converted into the acceleration of (nonthermal) particles E_nth conveys the major energy output, we would expect some correlation between the free energy time profile E_free(t) and the hard X-ray flux time profile F_HXR(t), which most easily can be inferred from the time derivative of the GOES soft X-ray time profile, i.e., F_HXR = ∂F_HXR(t)/∂t, according to the Neupert effect (Neupert 1968; Dennis & Zarro 1993).

We juxtapose these two time profiles E_diss(t) and F_HXR(t) for 20 flare events in Figures 5 and 6, where the time profiles of the energy dissipation (inferred from the pulse decay time intervals marked with gray areas) and the GOES 1–8 Å flux (marked with hatched areas) are shown. While there are obvious correlations between the two time intervals in a number of single-pulse flares (e.g., event #53 in Figure 5(a), #187 in Figure 5(c)), in double-pulse flares (e.g., event #367 in Figure 5(i)), or in multipulse flares (e.g., event #54 in Figure 6(d), #150 in Figure 6(e)), we see also surprising cases where hard X-ray emission is detected for a single pulse only when a sequence of five magnetic energy pulses is present (e.g., event #171 in Figure 6(j), #219 in Figure 6(i)). Thus, we find both, well-correlated flare events as well as mismatching time profiles. This outcome of our study indicates that the simple-minded notion of magnetic energy dissipation with subsequent particle acceleration does not always fit the data.

We consider two different scenarios of the time evolution of energy dissipation in solar flares: the slow-driven self-organized criticality (SOC) model (Figure 1(a)) and the fast-driven SOC model (Figure 1(b)). The slow-driven SOC model corresponds to the model of cosmic transients proposed by Rosner & Vaiana (1978), while their time evolution can also be characterized by an exponential-growth model, a power-law growth model, or a logistic-growth model (Section 3 of Aschwanden 2011a). Besides the application to solar flare observations, slow external driving of photospheric motion is expected to lead to occasional relaxation events also, at random times, with random amplitudes (Longcope & Sudan 1992). The essential property of the slow SOC model is the exponential growth of energy build-up during the time interval between two subsequent flare events, which eventually creates a flare at a random time interval, and then relaxes into a more stable state than before. The exponential-growth function, together with Poissonian random statistics, leads to the prediction of a power-law function of the flare size distribution (Rosner & Vaiana 1978). Moreover, the monotonic growth of the free energy predicts a correlation between the flare size and the interflare (waiting) time interval. However, observational searches for such a correlation between the flare sizes and flare waiting times turned out to be negative (Lu 1995; Crosby et al. 1998; Wheatland 2000a; Moon et al. 2001; Lippiello et al. 2010). The only correlation found was that smaller active regions produce smaller flare sizes (Wheatland 2000b) and that small active regions produce deviations from power laws (Wheatland 2010). There are also the problems that large flares sometimes occur within shorter waiting times than the required energy build-up times of the Rosner–Vaiana model: sometimes a larger flare volume than available is required, or too many e-folding growth times are necessary (Lu 1995). Nevertheless, a correlation of the flare size with the time interval after a flare (rather than before) was claimed for a small sample of flare events in the same active region (Hudson 2019). In summary, none of the predictions of the slow-driven SOC model of Rosner & Vaiana (1978) could be confirmed by solar flare observations.

Most of the simulations of the (frequency occurrence) size distributions of SOC avalanches assume a separation of timescales, which means that the avalanche duration τ_flare or energy dissipation timescale τ_diss is much shorter than the waiting time between two subsequent avalanches, i.e., τ_flare ≪ Δt_wait. If the input rate (e.g., of sand grains dripped on a sandpile) is sufficiently slow, the statistical properties of avalanche sizes and durations are expected not to change (Pruessner 2012). However, the observed statistics of solar flares was found to violate the timescale separation during the solar cycle maximum era (Aschwanden 2011a, 2011b, 2011c; Aschwanden & Freeland 2012), when the flare duration exceeded the waiting times, i.e., τ_flare ≳ Δt_wait, which we call a fast-driven SOC system. Because the mean waiting time $\langle {\rm{\Delta }}{t}_{\mathrm{wait}}\rangle$ is defined by the total duration T of the observations, divided by the total number N_ev of events (or intervals),

$$\begin{eqnarray}&&\langle {\rm{\Delta }}{t}_{\mathrm{wait}}\rangle =\displaystyle \frac{T}{{N}_{\mathrm{ev}}},\end{eqnarray} \tag{ 17 }$$

the mean waiting time decreases reciprocally with the number of events, and thus becomes shorter for a faster input rate, as shown in Figure 3 for a fast driver that has a factor of 10² higher event number, but also a factor of 10² shorter mean waiting time. As the 10 examples in Figure 6 demonstrate, a number of N_peak = 5–13 flare peaks occur in large flares, which represent elementary flare substructures (Aschwanden et al. 1998) that we interpret as individual energy dissipation events in a fast-driven SOC system. Thus, we detect rapid fluctuations of the free energy E_free(t) before, during, and after large flares in a fast-driven SOC system (Figures 5 and 6), but the free energy does not monotonically increase between two subsequent flares (Figure 1(a)). Hence, the fast-driven SOC model (Figure 1(b)) is more consistent with the observations than the slow-driven SOC model (Figure 1(a)).

If the free (magnetic) energy that is dissipated during a solar flare would all be stored in the corona, we should see a negatively dropping step function of the free energy E_free(t) during the flare duration (Figure 1(a)). One of the most detailed studies on the time evolution of the free energy shows a gradual build-up of free energy over two days, culminating with an X2.2 GOES-class flare and a simultaneous downward step in the free energy (Sun et al. 2012; Aschwanden et al. 2014b). However, discrepancies up to a factor of ≲10 have been noticed in the decrease of free energy during flares, when the standard Wiegelmann-NLFFF code (with preprocessing) was employed in addition to our VCA-NLFFF code (Aschwanden et al. 2014b), which was reduced down to a factor of ≲3 in recent refined magnetic modeling (Aschwanden 2019b). Besides the expected step functions, we also observe in the present study a number of pulses in the free energy that have a short rise time and decay time, on the order of ${\tau }_{\mathrm{rise}}\approx {\tau }_{\mathrm{decay}}\,\approx {\tau }_{\mathrm{pulse}}/2\approx 0.2\pm 0.1\,\mathrm{hr}$ = 12 ± 6 minutes (Figures 8(a), (b), (c)).

A puzzling question is what mechanism causes the relatively short rise time of the free energy? One mechanism that we know produces an increase of the free energy is the helical twisting by vertical currents (as it is incorporated in the VCA-NLFFF code used here), but then the twisting with subsequent untwisting produces a time-symmetric pulse in the free energy without a net energy transfer. Another possible mechanism is the coronal illumination effect, where the twisted loops are not visible in the initial flare phase, but become detectable when chromospheric evaporation starts to fill up the flare loops (Aschwanden et al. 2014a). A third possibility is chromospheric energy injection into the corona produced by energy transferred from the turbulent convection zone and photosphere into the corona, e.g., via anomalous current dissipation (Rosner et al. 1978). Such a scenario, with the ultimate energy source in the convection zone rather than in the corona, can draw large amounts of free energy to generate a flare without requiring coronal storage. Magneto-convection as seen in photospheric granulation cells has typical spatial scales of ≈1000 km and turnover times of ≈7 minutes, which produces new emerging flux on timescales close to the observed pulse rise times of τ_pulse ≈ 12 ± 6 minutes (Figure 8(a)). In conclusion, the time evolution of the free energy E_free(t) provides crucial constraints on how and where the flare energy is stored.

From the waiting time distribution, we can learn whether an SOC system is stationary or nonstationary, which means whether or not the mean flaring rate is constant, as a function of time. In the original SOC concepts of Bak et al. (1987), it was assumed that individual avalanches are statistically independent events, and thus the waiting time distribution should form a Poissonian (or exponential) distribution function. If there is a deviation from a Poissonian distribution apparent, individual avalanche events could not be independent events, such as in sympathetic flares (Moon et al. 2002, 2003; Wheatland 2002, 2006; Wheatland & Craig 2006). However, when the flaring rate is not constant, the resulting waiting time distribution can be calculated by summing the partial waiting time distributions for each flaring rate (Wheatland et al. 1998; Wheatland & Glukhov 1998; Wheatland 2000c) as we summarize in Section 2 of this paper (and in Section 5 of Aschwanden 2011a). Waiting time distributions of solar flare data generally show a power-law distribution with a slope of α_Δw ≈ 2–3, (Wheatland et al. 1998; Moon et al. 2001; Wheatland 2003; Aschwanden & McTiernan 2010; Kanazir & Wheatland 2010), which is explained here with a model that is based on on the universal reciprocal relationship between the (time-varying) mean flaring rate and the (time-varying) waiting time, and predicts a slope of α_Δt = 2. In summary, the nonstationary Poissonian model provides the most natural explanation for the observed power-law-like waiting time distributions.

Besides the nonstationary Poissonian model of a fast-driven SOC model, some alternative interpretations have been explored, too, by other studies in the literature. An energy balance model in terms of a master equation between energy build-up and energy loss by dissipation of free energy has been proposed (Wheatland & Glukhov 1998; Wheatland & Litvinenko 2001, 2002; Wheatland 2008, 2009, 2009). Other approaches use scaling laws from magnetic reconnection processes (Litvinenko 1996; Wheatland & Craig 2003, 2006). Alternative functions for waiting time distributions were also tested, finding that lognormal and inverse Gaussian distribution functions are more likely to fit the observations than the exponential function (Kubo 2008).