A MACROSCOPIC DESCRIPTION OF A GENERALIZED SELF-ORGANIZED CRITICALITY SYSTEM: ASTROPHYSICAL APPLICATIONS

We start with geometric parameters, such as a length scale L, a Euclidean area A, and a Euclidean volume V, embedded in a Euclidean space with a dimension of d =1, 2, or 3. Euclidean means space-filling here, while inhomogeneous structures are described by a fractal dimension D_d, which also depends on the Euclidean dimension d.

SOC phenomena (like avalanches on a sandpile) can be triggered by the infall of a single sand grain, and thus the causal consequence of a tiny input or disturbance can have an unpredictable magnitude of the outcome or nonlinear response of an SOC system. Henceforth, the geometric size of an SOC avalanche can cover a considerable range L₁ ⩽ L ⩽ L₂ from the size L₁ of a single sand grain to the finite size L₂ of the SOC system. If only next-neighbor interactions are allowed in an SOC system, such as in the Bak–Tang–Wiesenfeld (BTW) model (Bak et al. 1987), a continuous distribution of length scales L of avalanches is expected when sampled over a long time. Naturally, small avalanches have a higher probability to occur than large ones, because they can happen simultaneously at different places of a sandpile, while a large system-wide avalanche can occur only once at a time. So, we can ask the question about the probability distribution function (PDF), N(L)dL, of avalanches with size L to occur in an SOC system. In order to solve this problem, we proceed in the same way as the PDF of random processes is derived.

The simplest statistical distribution is obtained from rolling dice, by enumerating all possible outcomes. The PDF of outcomes of rolling one dice, two dice, and three dice is shown in Figure 1, the classical binomial distribution that approaches a Gaussian normal distribution (Figure 1) for a large number of dice, with possible outcomes of n ⩽ x ⩽ 6n for six-sided dice, while the PDF N(x)dx is a Gaussian function centered at x = n(6 + 1)/2.

Zoom In Zoom Out Reset image size

Going to the statistical probability distributions of avalanches with size L, we use the same method by enumerating all possible states with size L that can occur in an SOC system with finite size L₂. The case with a Euclidean space dimension of d = 2 is illustrated in Figure 2, where we use logarithmic bins with size x = 1, 2, 4, 8, 16. In a system with finite size L₂ = 16, one avalanche of this maximum size L₂ is possible in a given time interval, and thus N(x = 16) = 1. For a bin with half the size, L = L₂/2 = 8 we have four possible areas with a length scale of x = L₂/2, and thus N(x = 8) = 4. Proceeding to quarter bins, L = L₂/4, we have 16 possible areas with size L = L₂/4, and thus N(x = 4) = 16 = 2⁴, and so forth. Obviously, the probability distribution scales as N(x)dx = (L₂/L)² for Euclidean dimension d = 2. We can easily imagine the probabilities for the other Euclidean dimensions d = 1, which is N(x)dx = (L₂/L), and for d = 3, which is N(x)dx = (L₂/L)³. Therefore, we obtain a generalized probability distribution of length scales L according to

$$\begin{equation} N(L) dL \propto L^{-d} \ {\it dL}, \end{equation} \tag{ 1 }$$

which we call the scale-free probability conjecture (Aschwanden 2012a), being related to packing rules (e.g., sphere packing, or dense packing) in geometric aggregation problems. A similar approach of using geometric scaling laws was also pioneered for earthquakes (Main & Burton 1984). The term scale-free is generally used to express that no special scale is present in a statistical distribution, unlike the first moment or center value of a Gaussian (normal) distribution, or the e-folding value in an exponential distribution. Our scale-free probability does not require that all possible avalanches in an SOC system have to occur simultaneously, or in any particular sequential order. They just represent the expected distribution of a statistically representative sample, similar to the rolling of dice. For instance, using n = 10²⁶ dice to mimic the number of atoms per cm³, there is no way to execute all possible rolls, but we expect for any statistically representative subset of possible outcomes a Gaussian distribution. Similarly, we expect a length distribution N(L) according to Equation (1) for any statistically representative subset of avalanches occurring in an SOC system. We expect that Equation (1) has universal validity in SOC systems, because it is only based on a statistical argument of random processes on all scales, without any other constraints given by specific physical parameters or the dynamic behavior of an SOC system. This scale-free probability conjecture (Equation (1)) may also occur in other nonlinear systems, such as in turbulence. We may be able to discriminate between the two systems by the sparseness of avalanches (in slowly driven SOC systems) and the space-filling of structures (in turbulent media).

Zoom In Zoom Out Reset image size

In the following, we are going to derive size distributions of SOC avalanches by using geometric scaling laws, which is a standard approach that has been applied in a number of previous works (e.g., Bak et al. 1988; Robinson 1994; Munoz et al. 1999; Biham et al. 2001).

Besides the length scale L, other geometric parameters are the Euclidean area A or the Euclidean volume V. The simplest definition of an area A as a function of a length scale L is the square-dependence,

$$\begin{equation} A \propto L^2, \end{equation} \tag{ 2 }$$

which applies also to circular areas, A∝πr², or more complicated solid areas, differing only by a constant factor for self-similar geometric shapes. A direct consequence of this simple geometric scaling law is that the statistical probability distribution of avalanche areas is directly coupled to the scale-free probability distribution of length scales (Equation (1)), and can be computed by substitution of L(A)∝A^1/2 (Equation (2)), into the distribution of Equation (1), N(L) = N(L[A]) = L[A]^−d = (A^1/2)^−d = A^−d/2, and with the derivative dL/dA∝A^−1/2,

$$\begin{equation} N(A) {\it dA} \propto N(L[A]) \left| {dL \over dA} \right| dA \propto A^{-(1+d)/2}\ {\it dA}. \end{equation} \tag{ 3 }$$

Thus, we expect an area distribution N(A) depending on the dimensionality d = 2, 3 of the SOC system,

$$\begin{equation} N(A) {\it dA} \propto A^{-\alpha _A} {\it dA}, \qquad {\rm where}\ \left\lbrace \begin{array}{@{}ll}\alpha _A=1.5 & {\rm for}\ d=2 \\ \alpha _A=2.0 & {\rm for}\ d=3 , \\ \end{array} \right. \, \end{equation} \tag{ 4 }$$

which should also have universal validity for SOC systems. In spatially resolved astrophysical observations, such as of the Sun or magnetosphere, a length scale L or area A are the only directly measurable geometric parameters, while a volume V is generally derived from the observed area of an SOC event.

Similar to the area, we can derive the geometric scaling for volumes V, which simply scales with the cubic power in three-dimensional space,

$$\begin{equation} V \propto L^3, \end{equation} \tag{ 5 }$$

which represents a cube but differs only by a constant factor for a sphere, i.e., V = (4π/3)r³. Consequently, we can also derive the probability distribution N(V)dV of volumes V directly from the scale-free probability conjecture (Equation (1)), where the definition of Equation (5) demands d = 3. Substituting L∝V^1/3 into N(L[V])∝L[V]^d∝V^−d/3 and the derivative dL/dV = V^−2/3, we obtain for d = 3,

$$\begin{equation} N(V) dV \propto N(L[V]) \left| {dL \over dV} \right| dV \propto V^{-\alpha _V}\ dV \propto V^{-5/3}\ dV. \end{equation} \tag{ 6 }$$

Thus, a power law slope of α_V = 5/3 is predicted in three-dimensional Euclidean space, which applies also to the Euclidean volume of a time-integrated SOC avalanche in lattice simulations. However, since avalanches have a fractal geometry, it is the time-integrated fractal volume that is equivalent to the number of active pixels in a lattice simulation, rather than the Euclidean volume.

Since all the assumptions made so far are universal, such as the scale-free probability conjecture (Equation (1)) and the geometric scaling laws A∝L² and V∝L³, the resulting predicted occurrence frequency distributions of N(A) (Equation (3)) and N(V) (Equation (6)) are universal too, and power law functions are predicted from this derivation from first principles, which is consistent with the property of universality in theoretical SOC definitions.

“Fractals in nature originate from self-organized critical dynamical processes” (Bak & Chen 1989). Fractal geometries were pioneered in the context of self-similar structures before the advent of SOC models (Mandelbrot 1977, 1983, 1985) and have been applied to spatio-temporal SOC structures extensively (e.g., Bak et al. 1987, 1988; Bak & Chen 1989; Ito & Matsuzaki 1990; Feder & Feder 1991; Rinaldo et al. 1993; Erzan et al. 1995; Barabasi & Stanley 1995). Since the fractal geometry is a postulate of SOC processes invoked by the first pioneers of SOC, it is appropriate to approximate spatial structures of SOC avalanches by a fractal dimension. The simplest fractal is the Hausdorff dimension D_d, which is a monofractal and depends on the Euclidean space dimension d = 1, 2, 3. The Hausdorff dimension D₃ for the three-dimensional Euclidean space (d = 3) is

$$\begin{equation} D_3 = {\log {V_f(t)} \over \log {(L)} }, \end{equation} \tag{ 7 }$$

and analogously for the two-dimensional Euclidean space (d = 2),

$$\begin{equation} D_2 = {\log {A_f(t)} \over \log {(L)} }, \end{equation} \tag{ 8 }$$

with A_f(t) and V_f(t) being the fractal area and volume of an SOC avalanche during an instant of time t. These fractal dimensions can be determined by a box-counting method, where the area fractal D₂ can readily be obtained from images from the real world, while the volume fractal D₃ is generally not available unless one obtains three-dimensional data (or by numerical simulations).

A good approximation for the expected fractal dimension D_d is the mean value of the smallest possible fractal dimension D_{d, min} ≈ 1 and the largest possible fractal dimension D_{d, max} = d. The minimum possible fractal dimension is near the value of 1 because the next-neighbor interactions in SOC avalanches require some continuity between active nodes in a lattice simulation of a cellular automaton, while smaller fractal dimensions D_d < 1 are too sparse to allow an avalanche to propagate via next-neighbor interactions. Thus, the mean value of a fractal dimension is expected to be (Aschwanden 2012a)

$$\begin{equation} D_d \approx {D_{d,{\rm min}} + D_{d,{\rm max}} \over 2} = {(1+d) \over 2}. \end{equation} \tag{ 9 }$$

Thus, we expect fractal dimensions of D₃ ≈ (1 + 3)/2 = 2.0 for the three-dimensional space, and D₂ ≈ (1 + 2)/2 = 1.5 for the two-dimensional space. This conjecture of the mean value of the fractal dimension D_d has been numerically tested with cellular automaton simulations for Eucledian dimensions d = 1, 2, 3 and the following mean values were found: D₁ = 1.00  ±  0.00 (Aschwanden 2012a); then, D₂ = 1.58  ±  0.02 (Charbonneau et al. 2001), D₂ = 1.58 ± 0.03 (McIntosh et al. 2002), D₂ = 1.60 ± 0.17, 1.62 ± 0.18 (Aschwanden 2012a) for the two-dimensional case, for which D₂ = 1.5 is predicted, and D₃ = 1.78 ± 0.01 (Charbonneau et al. 2001; McIntosh et al. 2002), D₃ = 1.94 ± 0.27, 1.97 ± 0.29 (Aschwanden 2012a) for the three-dimensional case, for which D₃ = 2.0 is predicted. Thus, the mean value defined in Equation (9) is a reasonably accurate prediction based on the standard (BTW) cellular automaton model.

This relationship (Equation (9)) allows also a scaling between the fractal dimensions of the two-dimensional and three-dimensional Euclidean space,

$$\begin{equation} {D_3 \over D_2} \approx {(1+3) \over (1+2)} = {4 \over 3}. \end{equation} \tag{ 10 }$$

An extensive discussion of measuring the fractal geometry in SOC systems is given in Aschwanden (2011a, Chapter 8) and McAteer (2013). Fractals are measurable from the spatial structure of an avalanche at a given instant of time. Therefore, they enter the statistics of time-evolving SOC parameters, such as the observed flux per time unit, which is proportional to the number of instantaneously active nodes in a lattice-based SOC avalanche simulation.

The next important step is to include time scales, which together with the geometric scaling laws define the spatio-temporal evolution of SOC events. We model an SOC event simply as an instability that is triggered when a local threshold is exceeded. The universal behavior of any instability is an initial nonlinear growth phase and a subsequent saturation phase. We model the saturation phase with a diffusive function, as shown in Figure 3 (upper panel),

$$\begin{equation} r(t) = \kappa (t - t_0)^{\beta /2}, \end{equation} \tag{ 11 }$$

where t₀ is the onset time of the instability, κ is the diffusion coefficient, and β is the spreading exponent. A value of β ≳ 0 corresponds to logistic growth with an upper limit of the spatial volume (Aschwanden 2011a, 2012b), β ≈ 0.5 corresponds to subdiffusion, β = 1 to classical diffusion, β ≈ 1.5 to hyper diffusion or Lévy flight, and β = 2 to linear expansion.

Zoom In Zoom Out Reset image size

The corresponding velocity v(t) of an expanding SOC avalanche is shown in Figure 3 (second panel), which monotonically decreases with time and is obtained from the time derivative of r(t) (Equation (11),

$$\begin{equation} v(t) = {d r(t) \over dt} = {\kappa \beta \over 2}\ (t - t_0)^{\beta /2-1}. \end{equation} \tag{ 12 }$$

What spatio-temporal scaling law do we expect from this macroscopic description of an SOC avalanche? A spatial scale L could be defined from the maximum size of the avalanche at the end time T = (t − t₀), and thus we expect from Equation (11) the statistical spatio-temporal scaling law

$$\begin{equation} L \propto \kappa \ T^{\beta /2}. \end{equation} \tag{ 13 }$$

Substituting this scaling law L(T) into the PFD of length scales (Equation (1)), we expect a power law distribution of time scales,

$$\begin{equation} N(T) dT = N(L[T]) {dL\over dT} dT = T^{-[1+(d-1)\beta /2]} = T^{-\alpha _T}, \end{equation} \tag{ 14 }$$

with the power law slope of α_T = 1 + (d − 1)β/2, which has a value of α_T = 1 + β = 2.0 for three-dimensional Euclidean space (d = 3) and classical diffusion (β = 1). This power law slope for avalanche time scales is a prediction of universal validity, since it is only based on the scale-free probability conjecture (Equation (1)), N(L)∝L^−d, and the diffusive nature (or random-walk statistics) of the saturation phase.

The spatio-temporal scaling law (Equation (13)), based on random-walk or a diffusion process, is used here as a simple approximation in an empirical way. Diffusive transport has been applied to SOC theory and SOC phenomena in a number of previous studies, e.g., by using the spreading exponents to determine the critical points of systems with multiple absorbing states (Grassberger & de la Torre 1979), as a discretized diffusion process using the Langevin equation (Wiesenfeld et al. 1989; Zhang 1989; Foster et al. 1977; Medina et al. 1989), in terms of classical (Lawrence 1991) and anomalous diffusion of magnetic flux events (Lawrence & Schrijver 1993), in deriving spatio-temporal scaling laws with mean-field theory and branching theory (Vespignani & Zapperi 1998), as a continuum limit of a fourth-order hyper-diffusive system (Liu et al. 2001; Charbonneau et al. 2001), or in terms of a diffusion entropy description (Grigolini et al. 2002).

In numerical SOC simulations, such as in lattice-based cellular automaton models of the BTW type (Bak et al. 1987), energy is dissipated in every node that exceeds a threshold temporarily, and thus the energy that is dissipated during an SOC avalanche is proportional to the total number of all active nodes, summed over space at each instant of time. If we count these active nodes at a given time interval, we have a quantity that is proportional to the instantaneous energy dissipation rate, which has the unit of energy per time. In the real world, we observe a signal from an SOC avalanche in the form of an intensity flux f(t) (e.g., seismic waves from earthquakes, hard X-ray (HXR) flux from solar flares, or the amount of lost dollars per day in the stockmarket). Let us assume that this intensity flux is proportional to the volume of active nodes, which corresponds to the instantaneous fractal volume V_f(t) of an SOC avalanche in our spatio-temporal SOC model (Figure 4), also called fractal-diffusive (FD-SOC) model (Aschwanden 2012a),

$$\begin{equation} f(t) \propto V_f(t) \propto r(t)^{D_d}, \end{equation} \tag{ 15 }$$

which is shown in Figure 3 (third panel) for β = 0.1, 0.5, and 1. The flux time profile f(t) is expected to fluctuate substantially in real data or in lattice simulations, because the fractal dimension can vary in the range of D_{d, min} ≈ 1 and D_max = d, while we use only the mean value D_d = (D_min + D_max)/2 (Equation (9)) in our macroscopic model. Occasionally, the instantaneous fractal dimension may reach its maximum value, i.e., D_d(t) ≲ d, which defines an expected upper limit f_max(t) of

$$\begin{equation} f_{{\rm max}}(t) \propto V(t) \propto r(t)^d. \end{equation} \tag{ 16 }$$

Zoom In Zoom Out Reset image size

Integrating the time-dependent flux f(t) over the time interval [0, t] yields the total dissipated energy e(t) up to time t (using Equation (11)),

$$\begin{eqnarray} e(t) \propto \int _{0}^t V_f(t) dt &=& \int _{0}^t r^{D_d}(t) dt = \int _{0}^t \kappa ^{D_d}_d (t-t_0)^{D_d \beta /2} dt\nonumber\\ &=& {\kappa ^{D_d} \over D_d \beta /2+1} (t-t_0)^{D_d \beta /2 + 1}, \end{eqnarray} \tag{ 17 }$$

which is a monotonically increasing quantity with time (Figure 3, bottom panel).

From this time-dependent evolution of an SOC avalanche, we can characterize at the end time t a time duration T = (t − t₀), a spatial scale L = r(t = t₀ + T), an expected flux or energy dissipation rate F = f(t = t₀ + T), an expected peak flux or peak energy dissipation rate P = f_max(t = t₀ + T), and a dissipated energy E = e(t = t₀ + T), which is proportional to the avalanche size S in BTW models. Thus, we have the following scaling relations between the different SOC parameters and the length scale L (using Equations (15)–(17)),

$$\begin{equation} F \propto L^{D_d} \propto T^{D_d \beta /2}, \end{equation} \tag{ 18 }$$

$$\begin{equation} P \propto L^d \propto T^{d \beta /2}, \end{equation} \tag{ 19 }$$

$$\begin{equation} E \propto S \propto L^{D_d+2/\beta } \propto T^{D_d \beta /2+1}. \end{equation} \tag{ 20 }$$

An alternative notation for the diffusive spreading exponent β used in literature is D_T = 2/β, so that the spatio-temporal scaling law (Equation (13)) reads as $T \propto L^{D_T}$ and the energy scaling law (Equation (20)) as $E = S \propto V_f T \propto L^{D_d+D_T}$, which can be expressed as $S \propto L^{D_S}$ with the exponent D_S = D₃ + D_T. Slight variations of this scaling law have been inferred from observations in different wavelengths, such as D_S = D_A/2 + D_T for magnetic events (Equation (18) in Uritsky et al. 2013), which seems to be equally consistent with observations as our generalized (wavelength-independent) FD-SOC model (see EIT and MDI events from Uritsky et al. 2013 in Table 1).

Finally, we want to quantify the occurrence frequency distributions of the (smoothed) energy dissipation rate N(F), the peak flux N(P), and the dissipated energy N(E), which all can readily be obtained by substituting the scaling laws (Equations (18)–(20)) into the fundamental length scale distribution (Equation (1)), yielding

$$\begin{equation} N(F) dF = N(L[F]) \left| {dL \over dF} \right| dF \propto F^{-[1+(d-1)/D_d]} \ dF, \end{equation} \tag{ 21 }$$

$$\begin{equation} N(P) dP = N(/[P]) \left| {dL \over dP} \right| dP \propto P^{-[2-1/d]} \ dP, \end{equation} \tag{ 22 }$$

$$\begin{equation} N(E) dE = N(L[E]) \left| {dL \over dE} \right| dE \propto E^{-[1+(d-1)/(D_d+2/\beta)]} \ dE.\quad \end{equation} \tag{ 23 }$$

Thus, this derivation from first principles predicts power law functions for all parameters L, T, F, P, E, and S, which are the hallmarks of SOC systems. In summary, if we denote the occurrence frequency distributions N(x) of a parameter x with a power law distribution with power law index α_x,

$$\begin{equation} N(x) dx \propto x^{-\alpha _x} \ dx, \end{equation} \tag{ 24 }$$

we have the following power law coefficients α_x for the parameters x = L, T, F, P, E, and S,

$$\begin{equation} \begin{array}{r@{\ }l}\alpha _L &= d \\ \alpha _T &= 1+(d-1)\beta /2 \\ \alpha _F &= 1+(d-1)/D_d \\ \alpha _P &= 1+(d-1)/d \\ \alpha _E = \alpha _S &= 1+(d-1)/(D_d+2/\beta) .\\ \end{array} \end{equation} \tag{ 25 }$$

If we restrict the case to three-dimensional Euclidean space (d = 3), as is almost always the case for real-world data, the predicted power law indexes are

$$\begin{equation} \begin{array}{r@{\ }l}\alpha _L &= 3 \\ \alpha _T &= 1+\beta \\ \alpha _F &= 1+2/D_3 \\ \alpha _P &= 1+2/3 \\ \alpha _E = \alpha _S &= 1+1/(D_3/2+1/\beta).\\ \end{array} \end{equation} \tag{ 26 }$$

Restricting to classical diffusion (β = 1) and an estimated mean fractal dimension of D₃ ≈ (1 + 3)/2 = 2, we have the following absolute predictions

$$\begin{equation} \begin{array}{r@{\ }l}\alpha _L &= 3 \\ \alpha _T &= 2 \\ \alpha _F &= 2 \\ \alpha _P &= 5/3 \\ \alpha _E = \alpha _S &= 3/2 .\\ \end{array} \end{equation} \tag{ 27 }$$

The FD-SOC model predicts a power law distribution $N(T) \propto T^{-\alpha _T}$ of event durations T with a slope of α_T = [1 + (d − 1)β/2] (Equation (25)) that derives directly from the scale-free probability conjecture N(L)∝L^−d (Equation (1)) and the random walk (diffusive) transport (L∝T^β/2; Equation (13)). For classical diffusion (β = 1) and space dimension d = 3 the predicted power law is α_T = 2. From this time scale distribution, we can also predict the waiting time distribution with a simple probability argument. If we define a waiting time Δt as the time interval between the start time of two subsequent events so that no two events overlap with each other temporally, the waiting time cannot be shorter than the time duration of the intervening event, i.e., Δt_i ⩾ (t_{i + 1} − t_i). Let us consider the case of non-intermittent, contiguous flaring, but with no time overlap between subsequent events. In this case, the waiting times are identical with the event durations, and therefore their waiting time distributions are equal also, reflecting the same statistical probabilities,

$$\begin{equation} N(\Delta t) d\Delta t \propto N(T) dT \propto \Delta t^{-\alpha _{\Delta t}} d\Delta t, \end{equation} \tag{ 28 }$$

with the power law slope,

$$\begin{equation} \alpha _{\Delta t} = \alpha _T = 1 + (d-1)\beta /2. \end{equation} \tag{ 29 }$$

This statistical argument is true regardless what the order of subsequent event durations is, so it fulfills the Abelian property. Now we relax the contiguity condition and subdivide the time series into blocks with contiguous flaring (with intervals Δt ≈ T), interrupted by arbitrarily long quiet periods Δt = Δt^q when no events occur (Figure 5). The distribution of quiet periods Δt^q may be drawn from a random process, which has an exponential distribution

$$\begin{equation} N(\Delta t^q) d \Delta t^q \propto \exp {(-\Delta t^q / \Delta t^q_0)} \ d \Delta t^q. \end{equation} \tag{ 30 }$$

If we define a maximum event duration T₂ and assume that this is also approximately a lower limit for the quiet time intervals, i.e., $\Delta t^q_{{\rm min}} \approx T_2$, then we expect a power law distribution with a slope of α_Δt = α_T for the range of waiting times that are shorter than the maximum flare duration Δt ≲ T₂, with an exponential cutoff at Δt ≳ T₂. The contributions of waiting times from the subset of contiguous time blocks will still be identical, while those time intervals from the intervening quiet periods add some longer random waiting times. The predicted power law slope of short waiting times (T₁ ≲ Δt ≲ T₂) is then α_Δt = 2.0 for classical diffusion β = 1 and space dimension d = 3. Interestingly, this predicted slope is identical to that of nonstationary Poisson processes in the limit of intermittency (Aschwanden & McTiernan 2010). At the same time, this waiting time model predicts also clustering of events during active periods, and thus event statistics with memory and persistence, as it was demonstrated recently for CME events using Weibull distributions (Telloni et al. 2014).

Zoom In Zoom Out Reset image size

We can define a mean waiting time 〈Δt〉 from the total duration of the observing period T_obs and the number of observed events n_obs,

$$\begin{equation} \langle \Delta t \rangle = {T_{{\rm obs}} \over n_{{\rm obs}}}. \end{equation} \tag{ 31 }$$

From the distribution of event durations T, we have an inertial range of time scales [T₁, T₂], over which we observe a power law distribution, $N(T) \propto T^{-\alpha _T}$, with the corresponding number of events [N₁, N₂], so that we can define a nominal power law slope of α_T = log (N₂/N₁)/log (T₂/T₁). If the mean waiting time of an observed time series becomes shorter than the upper limit of time scales T₂ during very busy periods, we start to see time-overlapping events, a situation we call “event pile-up” or “pulse pile-up”. In such a case, we expect that the waiting time distribution starts to be modified, because the time durations of the long events are underestimated (by some automated detection algorithm), so that the nominal power law slope that is expected with no pulse pile-up, α_Δt = log (N₂/N₁)/log (T₂/T₁), has to be modified by replacing the lower time scale T₁ with the mean waiting time 〈Δt〉,

$$\begin{equation} \alpha _{\Delta t}^{{\rm pileup}} \ = \ \alpha _{\Delta t} \times \left\lbrace \begin{array}{@{}ll}1 & {\rm for}\ \langle \Delta t \rangle \ > \ T_2 \\ \log (T_2)/\log {\langle \Delta t \rangle } & {\rm for}\ \langle \Delta t \rangle \ \le \ T_2. \\ \end{array} \right. \end{equation} \tag{ 32 }$$

As a consequence, the measurements of event durations must suffer from the same pile-up effect, and a similar correction is expected for the time scale distribution N(T),

$$\begin{equation} \alpha _T^{{\rm pileup}} \ = \ \alpha _T \times \left\lbrace \begin{array}{@{}ll}1 & {\rm for}\ \langle \Delta t \rangle \ > \ T_2 \\ \log (T_2)/\log {\langle \Delta t \rangle } & {\rm for}\ \langle \Delta t \rangle \ \le \ T_2. \\ \end{array} \right. \end{equation} \tag{ 33 }$$

Thus, the predicted waiting time distribution has a slope of α_T = 2 in the slowly driven limit but can be steeper in the strongly driven limit. For instance, the waiting time distributions of solar flares correspond to the slowly driven limit during the minima of the solar 11 yr cycle, while the power law slopes indeed steepen during the maxima of the solar cycle (Aschwanden & Freeland 2012), when the flare density becomes so high that the slowly driven limit, and thus the separation of time scales, is violated.

Our fractal-diffusive SOC model developed so far has universal validity because it is entirely derived from statistical probabilities and fractal-diffusive transport. The predicted scaling laws and occurrence frequency distributions derived above do not depend on any specific physical parameter of an SOC phenomenon. Using real-world observations, however, some physical scaling laws are involved between the observables and the spatio-temporal parameters used so far. For instance, the strength of an earthquake is measured in magnitudes of the Gutenberg–Richter scale (Richter 1958), which may be related to the observed earthquake rupture area by some mechanical scaling law that determines the statistics (Main & Burton 1984). For solar flares, the observed fluxes in soft X-rays (SXRs) or HXRs are related to the physical parameters of electron temperatures, densities, and pressures of heated plasma, as it can be derived for the equilibrium point between heating and cooling (e.g., Rosner et al. 1978). Other scaling laws used in solar physics include, for instance, relationships between magnetic energies and the reduced MHD equations (Longcope & Sudan 1992), or the magnetic reconnection geometry (Craig 2001), or between the heating rate and the magnetic field strength (Schrijver et al. 2004). Such physical scaling laws allow us to derive the power law slope of the frequency distribution of both the observables and the physical parameters, which is examined elsewhere (e.g., Aschwanden et al. 2013).

The predicted frequency distributions for energies and fluxes derived in Section 2.5, are strictly only valid for systems where the assumption of proportionality between the flux and the instantaneous fractal volume is fulfilled, i.e., $F(t) \propto V_f(t) \propto r(t)^{D_d}$ (Equation (15)), because the scaling of observables depends then on geometric parameters only, which can be derived entirely from statistical probabilities, in terms of the scale-free probability conjecture (Equation (1)).

Without specializing on a particular physical mechanism of a given SOC system, we can give some general rules on how to derive the power law function of physical parameters. The simplest situation is a two-parameter correlation or scaling law, where a physical parameter x is related to the geometric length scale by a power law function with index γ,

$$\begin{equation} L \propto x^\gamma. \end{equation} \tag{ 34 }$$

Inserting this scaling law into the fundamental length scale distribution N(L)∝L^−d (Equation (1)) and using the derivative dL/dx = x^{γ − 1} then directly yields the occurrence frequency distribution N(x),

$$\begin{equation} N(x) dx = N(L[x]) {dL \over dx} dx = x^{-[1+(d-1)\gamma ]} dx. \end{equation} \tag{ 35 }$$

Also common is a three-parameter correlation or scaling law, such as in terms of two physical parameters x and y and the length scale L, i.e.,

$$\begin{equation} L \propto x^\gamma y^\delta, \end{equation} \tag{ 36 }$$

in terms of power law functions with exponents γ and δ. The probability distribution for one of the physical parameters, say y, can then be written as

$$\begin{eqnarray} N(y) dy &=& \int N(L, x[L,y]) dL \ dy\nonumber\\ &=& \int N(L) N(x[L,y]) dL \ dy \propto y^{a_x \delta / \gamma } \ dy,\qquad \end{eqnarray} \tag{ 37 }$$

after the integration over the variable L is carried out. Thus, the resulting distribution $N(y) dy \propto y^{-\alpha _y} dy$ has a power law slope of α_y = α_xδ/γ. The power law solution is strictly valid only for complete sampling of the parameters, which in reality is often not possible due to limited statistics, instrumental sensitivity limits, and data noise. This leads to truncation effects and finite-size effects, which can be simulated with Monte Carlo simulations or analytically calculated (see the Appendix in Aschwanden et al. 2013, for examples).

Besides physical scaling laws that are specific to a particular physical mechanism of an SOC system, there are also instrument-dependent scaling laws that are not universal and depend on the specific instrument used in an observation of SOC phenomena. If there is a nonlinear scaling between the observable and the geometric volume of an SOC avalanche, we cannot expect to measure the same power law slope of an observable with different instruments. In order to make observed frequency distributions obtained with different instruments compatible, it is often advisable to reduce the observable parameters to physical parameters using a well-established instrument calibration. For astrophysical observations in SXRs and extreme ultraviolet (EUV), for instance, an instrument-independent physical quantity is the differential emission measure distribution, which can be inverted from observed fluxes in different wavelengths (e.g., Aschwanden et al. 2013).

The original BTW model revealed the generic scale invariance of simulated or observed SOC parameters, which ideally exhibits power law functions for the occurrence frequency distributions, possibly related to the 1/f noise of power spectra (Bak et al. 1987). The property of a power law shape became the hallmark of SOC phenomena, but it was recognized that this is a necessary but not a satisfactory condition, since other phenomena (such as turbulence or percolation) produce power laws also.

Our fractal-diffusive SOC model (FD-SOC) derives the PFD based on a statistical probability argument, which leads to a power law function of spatial and geometric scales. The additional assumption of fractal-diffusive transport leads to a power law function of temporal scales. Further, we define the size of an avalanche from the time-integrated fractal volume that participates in an avalanche, and consequently, we obtain also power law distributions for the size or total dissipated energy of avalanches. Since all these assumptions are of a statistical nature and do not depend on any physical parameters of an SOC system, the predictions of the PDFs of spatial, temporal, and energy SOC parameters have universal applicability, regardless of the physical process that is involved in the nonlinear energy dissipation process. The prediction of a pure power law function for the size distributions at all scales is also called universality in theoretical SOC models (e.g., Sethna et al. 2001) and is fulfilled in the macroscopic description of our FD-SOC model by design (as a consequence of the scale-free probability conjecture; Equation (1)). However, we should be aware that this simple FD-SOC model provides only a first-order prediction, while additional effects (such as truncation, incomplete sampling, or finite-size effects) may modify the observed size distributions into broken power laws, double power laws, or other power law–like distribution functions. However, similar effects occur also in cellular automaton simulations.

SOC systems are expected to exhibit the spatio-temporal correlations (Jensen 1998) of an SOC state variable ${\bf B({\bf r},t)}$,

$$\begin{equation} C({\bf r}, t) = \langle B({\bf r}_0, t_0) B({\bf r}_0+{\bf r}, t_0+t) \rangle - \langle {\bf B}({\bf r}_0, t_0) \rangle ^2. \end{equation} \tag{ 38 }$$

Such correlations are absent in systems with random noise. In our FD-SOC model, however, the random structure of the background in a state near criticality is episodically disturbed by an avalanche event, which carves out a “hole” with a size L during a time scale T, which represents a major disturbance in the form of a spatially and temporally coherent structure, which can be restored to the critical state only gradually, for slowly driven SOC systems. Naturally, large avalanches leave their footprints behind and produce spatio-temporal correlations during the local restoration time. The correlation is best for large avalanches with similar shapes. The time profiles of avalanches in our FD-SOC system are self-similar to some extent, since they are characterized by a common fractal dimension D_d (Equation (15)), diffusion constant κ, and diffusive spreading exponent β (Equation (11)). We visualize the spatial correlations with a cartoon in Figure 6, which shows coherent disturbances as deviations from the critical state in large avalanches occurring in sandpiles and in solar flares.

Zoom In Zoom Out Reset image size

In our two-component model of waiting times (Section 2.6), an observed time series consist of quiet intervals Δt^q > T₂ with no avalanching (which have a random distribution), and active intervals with contiguous flaring (which have a power law distribution like the event durations N(Δt) = N(T)). This dual behavior is also called intermittency and has the consequence that the combined waiting time distribution has both a power law range (T₁ ⩽ Δt ⩽ T₂) and an exponential cutoff (Δt ⩾ T₂). Consequently, we expect spatio-temporal correlations (Equation (38)) during the intermittently active periods only, while they are expected to be absent during the quiet time intervals. Avalanching during active periods is also expected to exhibit persistence and memory, while no memory is expected during quiet time intervals. This property seems to be more consistent with observations (e.g., Telloni et al. 2014), but is different from the pure random (Poisson) statistics of the original BTW model. However, it reconciles related debates about the functional shape of the waiting time distributions (e.g., Boffetta et al. 1999; Lepreti et al. 2000).

Classical SOC systems operate in the limit of slow driving, which implies a separation between the duration of an avalanche and the waiting time interval between two subsequent avalanches. Numerically, the separation of time scales is simply realized by allowing only one single disturbance of an SOC system at a time, which triggers an avalanche (with duration T) or not, while the next disturbance is not initiated after a waiting time Δt > T, in the case of an avalanche.

In our FD-SOC model, the energy dissipation rate (de/dt) grows monotonically after a triggering disturbance, which exceeds the system-wide threshold value (de/dt)_crit until the spatial diffusion stops after time T, due to a lack of unstable nodes among the next-neighbors of an instantaneous avalanche area or volume. Therefore, the energy dissipation rate during an avalanche exceeds the threshold value during the entire duration of an avalanche. Energy conservation between the slowly driven energy input rate and the intermittent avalanching output rate can therefore only be obtained with sufficiently long waiting times Δt during which the energy loss of an avalanche is restored. This requires a balance of the long-term averages of the energy input and output rates, i.e.,

$$\begin{equation} \langle (de/dt)_{{\rm in}} \rangle \langle \Delta t \rangle \approx \langle (de/dt)_{{\rm out}} \rangle \langle T \rangle. \end{equation} \tag{ 39 }$$

Since 〈(de/dt)_in〉 ⩽ (de/dt)_crit ≪ 〈(de/dt)_out〉, it follows that 〈Δt〉 ≫ 〈T〉, which warrants a separation of the time scales, i.e., the waiting time Δt and the avalanche duration T.

The resulting time profile of the energy dissipation rate (de/dt) of an SOC system is then necessarily highly intermittent due to the long waiting times in between subsequent avalanches. In addition, the time profile is strongly fluctuating during an avalanche, according to $f(t) \propto r(t)^{D_d}$ (Equation (15)), since the fractal dimension D_d(t) can fluctuate in the entire range between the minimum and maximum value as a function of time, i.e., 1 ≲ D_d(t) ⩽ d. However, for the scaling laws in the FD-SOC model (Equations (18)–(20)), we can replace the fluctuating value of D_d(t) with a constant mean value 〈D_d(t)〉 = (1 + d)/2 and obtain the same size distributions.

How does our fractal-diffusive SOC model reinforce self-organized criticality? In classical SOC models, criticality is obtained by a slowly driven input of energy which restores the energy losses of avalanches until the system-wide critical threshold is reached (more or less) and new avalanches can be triggered by a local excess of the critical threshold. In our FD-SOC model, the time evolution of an avalanche has a generic shape that is given by fractal-diffusive transport, while the energy balance between energy input (disturbances) and output (avalanches) is not explicitly reinforced, unlike cellular automaton models which iterate a mathematical redistribution rule to drive the dynamics of an SOC system and are designed to conserve energy. Instead, self-organization of the FD-SOC model is constrained by statistical probability only, which does not need to be self-tuning to produce a particular functional form of a size distribution, because there is only one statistical distribution with maximum likelihood, which is a power law distribution function of spatial scales (according to our scale-free probability conjecture). So, we can say that the FD-SOC model gravitates around the statistically most likely state, like entropy in self-contained statistical systems without external influence. This may be a more general definition of SOC then originally proposed by Per Bak and coworkers but explains the concept of self-organization by the most general principle of maximum statistical likelihood. This should not surprise us, since the entire evolution of our universe followed maximum statistical likelihood, from the initial big bang expansion all the way to the bio-chemical evolution of life, forming complexity out of simple structures based on processes that are driven by statistical likelihood (e.g., Mendel's law in genetics).

If we mount a large container with a gel-like surface below a circular plate that holds Per Bak's sandpile, we would record impact craters from each sandpile avalanche in the viscous gel and could infer the avalanche sizes from the diameters of the impact craters (see experimental setup of sandpile experiment conducted by Held et al. 1990). Similarly, the Moon was targeted by many impacting meteors and meteorites, especially during an intense bombardment in the final sweep-up of debris at the end of the formation of the solar system between 4.6 and 4.0 billion years ago (e.g., Neukum et al. 2001). The sizes of lunar craters were measured with the first lunar spacecraft (Ranger 7, 8, 9) in the early 1960's, and a cumulative power law distribution with sizes in the range of L ≈ 10⁰–10^4.5 cm was found, with a power law slope of $\alpha ^{{\rm cum}}_L \approx 2.0$ for the cumulative distribution (Cross 1966), which corresponds to a value of $\alpha _L=\alpha ^{{\rm cum}}_L+1=3.0$ for the differential size distribution. This quite accurate result (for a size distribution covering a range of over four orders of magnitude) corresponds exactly to our prediction of the scale-free probability conjecture, N(L)∝L⁻³ (Equation (1)). A similar value of α_L = 2.75 was found for the size distribution of meteorites and space debris from man-made rockets and satellites (Figure 3.11 in Sornette 2004). The formation of the sizes of meteors and meteorites may have been controlled by a nonlinear process that includes a combination of self-gravity, gravitational disturbances, collisions, depletions, fragmentation, and captures of incoming new bodies in the solar system (e.g., Ivanov 2001). The Moon acts as a target that records the sizes of impacting meteorites that were produced by an SOC process, similar to the gel-filled plate under Bak's sandpile.

The origin of the asteroid main belt is believed to be associated with a time period of intense collisional evolution shortly after the formation of the planets (e.g., Botke et al. 2005). The asteroids are a leftover of the planetesimals that were either too small to form a planet by self-gravitation, or they orbited in an unstable region of the solar system that constantly got disturbed by the largest planets Jupiter and Saturn.

In Table 1, we compile some values of measured size distributions of asteroids, given as power law slopes α_L of the differential size distributions (related to the slope $\alpha ^{{\rm cum}}_L$ of the cumulative size distribution by $\alpha _L=\alpha ^{{\rm cum}}_L+1$, which includes values in the range of α_L ≈ 2.3–4.0, obtained from the Palomar Leiden Survey (Van Houten et al. 1970), the Spacewatch Surveys (Jedicke & Metcalfe 1998), the Sloan Digital Sky Survey (Ivezić et al. 2001), and the Subaru Main-Belt Asteroid Survey (Yoshida et al. 2003; Yoshida & Nakamura 2007). These values of the power law slopes agree within ≈25% with our theoretical prediction of α_L = 3.0, but the statistical range of sizes covers less than two decades, and thus incomplete sampling of small sizes is likely to limit the accuracy.

The Saturn ring extends over a range of 7000–80,000 km above Saturn's equator and has a mass of 3 × 10¹⁹ kg, consisting of a myriad of small particles with sizes in the range from 1 mm to 20 m (Zebker et al. 1985; French & Nicholson 2000). The particle size distribution was measured in eight different ring regions with Voyager I radio occultation measurements (Zebker et al. 1985). These size distributions were found to have slightly different power law slopes in each ring zone, with values of α_L = 2.74–3.03 for ring A, α_L = 2.79 for the Cassini division, and α_L = 3.05–3.22 for ring C (Zebker et al. 1985). Averaging the values from all eight zones, we find α_L = 2.89 ± 0.16, which is remarkably close to the prediction α_L = 3.0 of the scale-free probability conjecture (Equation (1)). Thus, the fragmentation of Saturn ring particles is consistent with the statistics of SOC avalanches, and the process of collisional fragmentation driven by celestial mechanics can be considered as a self-organizing system that is constantly driven toward the collisional instability threshold. An instability occurs by a collision of particles. If the system has a too low density, no collisions occur and the system is subcritical, while a too high density of particles would result into an excessive collision rate that would destroy the structure of the Saturn ring. Hence, the long-lived Saturn ring can be considered as an SOC system that self-tunes to a critical collisional limit that maintains its shape and conserves its (kinetic) energy, similar to Bak's SOC sandpile that maintains its slope and conserves the potential energy.

The Earth's magnetosphere displays a number of phenomena that have been associated with SOC models (Table 1), such as active and quiet substorms and auroral events (Lui et al. 2000; Uritsky et al. 2001, 2002, 2006; Kozelov et al. 2004; Klimas et al. 2010), substorm flow bursts (Angelopoulos et al. 1999), auroral electron (AE-index) bursts (Takalo 1993; Takalo et al. 1999), upper auroral (AU-index) bursts (Freeman et al. 2000; Chapman & Watkins 2001), or outer radiation belt electron events (Crosby et al. 2005). The power law indexes of observed size distributions of these phenomena are listed in Table 1.

Accurate measurements, using the same definition of time-integrated avalanche sizes as in the BTW model (Bak et al. 1987; Charbonneau et al. 2001) and in this paper, were carried out for auroral events in UV by Uritsky et al. (2002), and in visible light by Kozelov et al. (2004), yielding size distribution with power law slopes of α_A ≈ 1.7–2.0, α_T ≈ 2.0–2.5, α_P ≈ 1.66–2.0, and α_E = 1.4–1.7, which agree well with the predictions of the FD-SOC model (α_A = 2.0, α_T = 2.0, α_P ≈ 1.67, α_E = 1.5) (Table 1). The earlier reported lower values for the power law slopes of auroral fluences (Lui et al. 2000) are incompatible with recent observational results as well as with the FD-SOC model, because the auroral sizes were measured from snapshots taken in regular time intervals, rather than measured individually for each avalanche event (Uritsky et al. 2002). This case with contradicting statistical results measured from the same data is an example of a validation test using the FD-SOC model.

The number of electrons in the outer radiation belt (at 4–8 L-shell distances) is modulated by the solar wind, exhibiting size distributions of electron peak fluxes with power law slopes of α_P ≈ 1.5–2.1 (Crosby et al. 2005). The variation of the power law slope is mostly attributed to variations of the orbits of the microsatellites (STRV-1a and 1b) that record the electron bursts at different intersections of the radiation belt with the orbits. Nevertheless, the mean value averaged over different years and L-shell distances, α_P = 1.7 ± 0.2, is quite consistent with the theoretical prediction α_P = 1.67 of the FD-SOC model. The radiation belt can be considered as an SOC system, where the input is driven by solar wind electrons, which become trapped in the outer radiation belt, while magnetic variations modulate the untrapping of electrons by a self-organizing loss-cone angle, producing avalanches of electrons bursts.

Solar flares have been interpreted as an SOC phenomenon since 1991 (Lu & Hamilton 1991) and numerous studies have been performed to establish the size distributions of various solar flare parameters measured in HXRs, SXRs, EUV, and radio wavelengths. A representative selection of power law slopes from size distributions of solar flare length scales (α_L), flare areas (α_A), time durations (α_T), peak fluxes (α_P), and fluences or energies (α_E) is given in Table 1 (see references in footnote of Table 1). We note that most of the power law slopes measured in HXR, SXR, and EUV agree well with the theoretical predictions of our FD-SOC model, i.e., α_L = 3.0, α_A = 2.0, α_T = 2.0, α_P = 1.67, and α_E = 1.50, say typically within 5% to 10%. The remaining differences can be attributed to the different instrumental bias and the different analysis methods (threshold definition, preflare background subtraction, temperature bias) of the observations. Also the peak fluxes observed in radio wavelengths are commensurable with the predictions for incoherent emission mechanisms, such as gyrosynchrotron emission in microwave bursts. Only the solar energetic particles (SEPs) appear to have a flatter distribution than predicted, which has been interpreted in terms of a selection bias for large events (Cliver et al. 2012), or alternatively in terms of the geometric dimensionality of the SOC system (Kahler 2013). In summary, except for the SEP events, solar flares observed in almost all wavelengths are in agreement with the FD-SOC model and provide the strongest support for SOC models among all astrophysical phenomena.

What are the physical mechanisms in an SOC system that produce solar flares. The solar corona is considered to be a multi-component SOC system, where each active region or quiet Sun region represents a different SOC sandpile, with its own spatial (finite-size) boundary, lifetime, and flaring rate. Interestingly, the statistics of a single SOC system (one active region) seems not to be significantly different from the statistics of an ensemble of SOC systems (in the entire corona), except for a different largest-event cutoff (Kucera et al. 1997). The energy input comes ultimately from build-up of non-potential magnetic fields (with electric currents) that is driven by subphotospheric magneto-convection and magnetic flux emergence. The coronal SOC system is slowly driven by continuous emergence of magnetic flux, braiding, and stressing of the magnetic field. Once a local threshold for instability is exceeded (kink instability, torus instability, tearing mode instability, etc.), an avalanche of magnetic energy dissipation is triggered that ends in a fractal-diffusive phase. The dissipated energy can be converted into thermal energy of heated plasma (visible in SXRs and EUV), and into kinetic energy of accelerated particles (detectable in HXRs and in gyrosynchrotron emission in radio wavelengths). The fact that we measure similar power law slopes in all wavelengths (HXR, SXR, EUV, radio) implies that all converted energies are approximately proportional to the emitting volume, i.e., E∝V_fT. Only SEP events and coherent emission in radio wavelengths show a much flatter power law slope, which indicates a nonlinear scaling law E∝(V_fT)^γ or a selection bias for large events.

Stellar flares have been observed in small numbers during a few hours with the Hubble Space Telescope (Robinson et al. 1999), the Extreme Ultraviolet Explorer (EUVE; Audard et al. 2000; Kashyap et al. 2002; Güdel et al. 2003; Arzner & Güdel 2004; Arzner et al. 2007), and the XMM-Newton (Stelzer et al. 2007), which produced size distributions of flare energies (time-integrated EUV fluxes) with power law slopes in a range of α_E = 2.17 ± 0.25 (Table 1). These values are significantly steeper than derived for solar flare energies (α_E ≈ 1.5–1.6) but are expected for small samples near the exponential fall-off at the upper end of the size distribution (Aschwanden 2011a). In addition, since plasma cooling extends the SXR and EUV flux beyond the time interval of energy release, the fluence of the largest flares may be over estimated for the largest solar and most stellar flare events.

Much larger statistics of stellar flares became available recently from the Kepler mission: 373 flaring stars were identified in a search for white-light flares on ≈23,000 cool dwarfs in the Kepler Quarter 1 long cadence data (Walkowicz et al. 2011; Maehara et al. 2012); a total of 1547 superflares (several orders of magnitude larger than solar flares) were detected on 279 G-type (solar-like) stars (Notsu et al. 2013; Shibayama et al. 2013). The flare energies were estimated from the time-integrated bolometric luminosity in visible light. Similar energy size distributions were found as in earlier smaller samples (with EUVE), with power law slopes of α_E = 2.0 ± 0.2 for flares on all G-type stars, and α_E = 2.3 ± 0.3 for flares on slowly rotating G-type stars (Maehara et al. 2012; Shibayama et al. 2013). We show the size distribution for the total sample of 1538 stellar flares in Figure 7 (middle panel), which has a power law slope of α_E = 2.04 ± 0.13. From Kretzschmar (2011, Table 1 therein), we derive a scaling law between the bolometric fluence (total solar irradiance; which is equivalent to the bolometric energy E_b) and the GOES 1–8 Å peak flux P_x (Figure 7, top panel),

$$\begin{equation} E_b \propto P_x^{(0.78\pm 0.13)}. \end{equation} \tag{ 40 }$$

Using this scaling law we can derive the distribution of GOES peak fluxes of the stellar flares,

$$\begin{equation} N(P_x) {\it d P}_x \propto N(E_b[P_x]) {{\it dE}_b \over {\it d P}_x}\ {\it d P}_x \propto P_x^{-1.81\pm 0.12}\ {\it d P}_x, \end{equation} \tag{ 41 }$$

which is consistent with the size distribution of GOES fluxes directly obtained by applying the scaling law of Kretzschmar (2011) given in Equation (40), with a power law slope of α_P = 1.88 ± 0.09 (Figure 7, bottom). Interestingly, the so obtained peak flux α_P ≈ 1.88 agrees better with the theoretical prediction α_P = 1.67 of our FD-SOC model, than the bolometric fluence. This may indicate that the bolometric fluence is not an accurate proxy of the flare energy or flare volume, possibly due to a nonlinear scaling of the bolometric fluence with flare energies.

Zoom In Zoom Out Reset image size

Pulsars exhibit intermittent irregular radio pulses, besides the regular periodic pulses that are synchronized with their rotation period. The irregular pulses indicate some glitches in the positive spin-ups of the neutron star, possibly caused by sporadic unpinning of vortices that transfer momentum to the crust (Warzawski & Melatos 2008), which was interpreted as an SOC system (Young & Kenny 1996). A size distribution of the radio fluxes from the Crab pulsar exhibited a power law distribution with slopes in the range of α_P = 3.06–3.50 (Argyle & Gower 1972; Lundgren et al. 1995). A similar value of α_P = 2.8 ± 0.1 was found for PSR B1937+21 (Cognard et al. 1996). Statistical measurements of the size distribution of pulsar glitches obtained from about a dozen other pulsars yielded a large scatter of values in the range of α_P = −0.13, …, 2.4 (Melatos et al. 2008). The reasons for these inconsistent values may be rooted in the small-number statistics (N = 6, …, 30) and methodology (rank-order plots). If the more reliable values from the crab pulsar hold up (α_P ≈ 3.0), which are typical for size distribution of length scales (α_L = 3.0), physical models that predict a proportionality between peak fluxes P and length scales L should be considered.

Soft gamma-ray repeaters, detected at energies of >25 keV with the Compton Gamma Ray Observatory (CGRO) and the Rossi X-ray Timing Explorer, exhibited size distributions with fluences in the range of α_E = 1.43–1.76 (Gogus et al. 1999, 2000), which is quite consistent with the values measured from solar flares at the same energies and predicted by our FD-SOC model (α_E = 1.5). However, the physical mechanisms of soft gamma-ray repeaters are entirely different from solar flares, believed to originate from slowly rotating, extremely magnetized neutron stars that are located in supernova remnants (Kouveliotou et al. 1998, 1999), where neutron star crust fractures occur, driven by the stress of an evolving, ultrastrong magnetic field in the order of B ≳ 10¹⁴ G (Thompson & Duncan 1996). The fact that solar flares and soft gamma-ray repeaters exhibit the same energy size distribution, although the underlying physical processes are entirely different, supports the universal applicability of our FD-SOC model.

Cygnus X-1, the first galactic X-ray source that has been identified as a black-hole candidate, emits HXRs pulses with a time variability down to 1 ms, which is attributed to bremsstrahlung X-ray pulses from mass infalling toward the black hole and the resulting turbulence in the accretion disk. Observations with Ginga and Chandra exhibit complex 1/f noise spectra and size distributions of peak fluxes with very steep power law slopes of α_P ≈ 7.1 (Negoro et al. 1995; Mineshige & Negoro 1999), which have been interpreted in terms of SOC models applied to accretion disks (Takeuchi et al. 1995; Mineshige & Negoro 1999). Such steep values of the power law slope of peak fluxes are difficult to understand in terms of our standard FD-SOC model, which predicts α_P = 1.67. They exclude a linear scaling between the peak flux P and the emitting volume V covered by an X-ray pulse. Such a steep slope can only be produced by an extremely weak dependence of the X-ray peak flux P on the avalanche volume V, requiring a quenching mechanism that limits every fluctuation to almost the same level. The cellular automaton model of Mineshige & Negoro (1999), which can produce power law size distributions with such steep slopes of α_P ≈ 7, indeed prescribes a non-random distribution of time scales for large pulses (shots), where the occurrence of large pulses is suppressed for a certain period after each large pulse.

Blazars (BL Lacertae objects) are high-polarization quasars and optically violent variable stars, which exhibit a high degree of fluctuation in radio and X-ray emission due to their particular orientation with the jet axis almost coaligned with our line of sight. Light curves from GX 0109+224 were analyzed and found to exhibit a 1/f noise spectrum, i.e., P(ν)∝ν^−p with p = 1.57–2.05, and a size distribution of peak fluxes with a power law slope of α_P = 1.55, and have been interpreted in terms of an SOC model (Ciprini et al. 2003). This value is quite consistent with the prediction of our standard FD-SOC model (α_P = 1.67), which suggests that the peak flux P emitted (in optical and radio wavelengths) is proportional to the emitting volume V. The agreement between observations and the theoretical prediction supports the universal applicability of the FD-SOC model.

Cosmic rays are high-energetic particles that propagate through a large part of our universe and are accelerated by galactic and extragalactic magnetic fields. Cosmic-ray energy spectra span over a huge range of E = 10⁹–10²¹ eV, where the lower limit of ≈1 GeV corresponds to the largest energies that can be accelerated in solar flares and coronal mass ejections. This cosmic-ray energy spectrum exhibits an approximate power law function with a mean slope of α_E ≈ 3.0 (Figure 8, bottom right). A more detailed inspection reveals actually a broken power law with a slope of α_E1 ≈ 2.7 below the knee at E_knee ≈ 10¹⁶ eV, and a slope of α_E2 ≈ 3.3 above the knee. The two energy regimes are associated with the particle origin in galactic space (E ≲ E_knee) and extragalactic space (E ≳ E_knee).

Zoom In Zoom Out Reset image size

If we interpret a cosmic-ray energy spectrum as a size distribution of particle energies, we can apply our universal SOC model. The driver of the SOC system is a generation mechanism of seed populations of charged particles, which are mostly bound to astrophysical objects in a collisional plasma. A critical threshold is given by transitions of the particles from collisional to collisionless plasma (such as in the “run-away regime”), where a particle can freely be accelerated, either by Fermi first-order or diffusive shock acceleration. The subsequent particle transport combined with numerous acceleration steps during every passage of suitable electric fields or shock fronts represents the build-up of an avalanche, until the particle hits Earth's upper atmosphere where it is detected by a shower of secondary particles. If we could observe all end products of an avalanche, we would expect an energy spectrum of α_E = 1.5. In reality, the energy spectrum of cosmic rays is α_E ≈ 3.0, assuming that the detected energies are proportional to the avalanche volume. How can we explain this discrepancy? A power law index of this value is expected for the size distribution of length scales, N(L)∝L⁻³ (Equation (1)). Therefore, a similar energy spectrum of N(E)∝E⁻³ can only be produced if the energy E is proportional to the length scale L, requiring that the fractal volume $V_f(t) \propto r(t)^{D_3}$ has a fractal dimension of D₃ = 1.0. Such a scaling can be arranged if only a linear subvolume of an entire three-dimensional avalanche is observed, which is indeed the case for in situ detection at Earth, since the origin of the cosmic-ray avalanche is located far away. The situation is visualized in Figure 8. In solar flares, on the other hand, almost all energetic particles accelerated during a flare lose their energy in the chromosphere, and thus, we can detect the entire energy content of an avalanche event by remote-sensing. This is not possible in cosmic rays, because we cannot detect in situ all energy losses of cosmic-ray particles that originated isotropically from the same avalanche, in a remote place such as in a supernova or black hole. Another aspect that our FD-SOC model predicts is the random walk diffusion during its propagation, which is consistent with the current thinking of cosmic-ray particle transport.

Adopting the resulting scaling law between the energy E of a cosmic-ray particle and the Euclidian length scale L that a cosmic-ray particle has traveled at the time of detection, E∝L, which is expected for direct electric field acceleration in a voltage drop, as well as for any other particle acceleration mechanism with a fixed amount of energy extraction per distance increase, dE/dL ≈ const, we can even determine the distance to the origination site of the cosmic-ray particle. The cosmic-ray spectrum shown in Figure 8 has a knee at E_free ≈ 0.5 × 10¹⁶ eV, which marks the distance of the Earth to the center of our galaxy (L_gal ≈ 50 lt-yr or L_gal ≈ 5 × 10²² cm). The Euclidean distance where a cosmic-ray particle with a maximum energy of E_max ≈ 5 × 10²⁰ erg originated is then,

$$\begin{equation} L_{{\rm max}} \approx L_{{\rm gal}} \left({E_{{\rm max}} \over E_{{\rm knee}}} \right), \end{equation} \tag{ 42 }$$

which yields L_max ≈ 5 × 10²⁷ cm, which corresponds to about 10% of the size of our universe (R_uni ≈ 4 × 10²⁸ cm). Since the intergalactic and extragalactic magnetic fields have different field strengths, the diffusion coefficient of cosmic-ray particles is also expected to be different in these two regimes, which may explain the slightly different power law slopes below and above the galactic boundary L_gal and the related energy E_knee.

The lowest energies of the cosmic-ray spectrum are at E_min ≈ 10⁹ eV. Using the same linear scaling of energy with length scale,

$$\begin{equation} L_{{\rm min}} \approx L_{{\rm gal}} \left({E_{{\rm min}} \over E_{{\rm knee}}} \right), \end{equation} \tag{ 43 }$$

we estimate a distance of L_min ≈ 10¹⁶ cm or 200 AU, which is located somewhat outside of the termination shock of our heliosphere.