Self-organized Criticality in Solar and Stellar Flares: Are Extreme Events Scale-free?

An ideal power-law distribution that represents a differential occurrence frequency distribution, also simply called a size distribution N(x) of events, as a function of some size parameter x, can be described by four parameters (${x}_{1},{x}_{2},{n}_{0},{\alpha }_{x}$),

$$\begin{eqnarray}&&N(x){dx}={n}_{0}\ {x}^{-{\alpha }_{x}}{dx},\qquad {x}_{1}\leqslant x\leqslant {x}_{2},\end{eqnarray} \tag{ 1 }$$

where x₁ and x₂ are the lower and upper bounds of the power-law inertial range, α_x is the power-law slope, and n₀ is the normalization constant,

$$\begin{eqnarray}&&{n}_{0}={n}_{\mathrm{ev}}(1-{\alpha }_{x}){\left[{\left({x}_{2}\right)}^{1-{\alpha }_{x}}-{\left({x}_{1}\right)}^{1-{\alpha }_{x}}\right]}^{-1},\end{eqnarray} \tag{ 2 }$$

with n_ev being the total number of events contained in the inertial range [x₁ ≤ x ≤ x₂].

The ideal power-law function (Equation (1)) can be generalized with an additional “shift” parameter x₀ (with respect to x), also called the Lomax distribution (Lomax 1954), generalized Pareto distribution (Hosking & Wallis 1987), or thresholded power-law size distribution (Aschwanden 2015),

$$\begin{eqnarray}&&N(x){dx}={n}_{0}{\left({x}_{0}+x\right)}^{-{\alpha }_{x}}{dx},\end{eqnarray} \tag{ 3 }$$

with the normalization constant n₀ for the range ${x}_{1}\leqslant x\leqslant {x}_{2}$,

$$\begin{eqnarray}&&{n}_{0}={n}_{\mathrm{ev}}(1-{\alpha }_{x}){\left[{\left({x}_{2}+{x}_{0}\right)}^{1-{\alpha }_{x}}-{\left({x}_{1}+{x}_{0}\right)}^{1-{\alpha }_{x}}\right]}^{-1}.\end{eqnarray} \tag{ 4 }$$

This additional parameter x₀ accommodates three different features: truncation effects due to incomplete sampling of events below some threshold (if x₀ > 0), incomplete sampling due to instrumental sensitivity limits (if x₀ > 0), or subtraction of event-unrelated background (if ${x}_{0}\lt 0$), as it is common in astrophysical data sets.

While Equations (1) and (3) represent differential size distributions N_diff(x), it is statistically more advantageous to employ cumulative size distributions N_cum(x), especially for small data sets and near the upper cutoff, where we deal with a small number of events per bin. We will use cumulative size distributions that include all events accumulated above some size x, such as for the ideal power-law function,

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\mathrm{cum}}(\gt x){dx} & = & {\int }_{x}^{{x}_{2}}{n}_{0}\ {x}^{-{\alpha }_{x}}{dx}\\ & = & 1+({n}_{\mathrm{ev}}-1)\left(\displaystyle \frac{{x}_{2}^{1-\alpha }-{x}^{1-\alpha }}{{x}_{2}^{1-\alpha }-{{x}_{1}}^{1-\alpha }}\right).\end{array}\end{eqnarray} \tag{ 5 }$$

or for the thresholded power-law size distribution,

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\mathrm{cum}}(\gt x) & = & {\int }_{x}^{{x}_{2}}{n}_{0}{\left(x+{x}_{0}\right)}^{-{\alpha }_{x}}{dx}\\ & = & 1+\,({n}_{\mathrm{ev}}-1)\left(\displaystyle \frac{{\left({x}_{2}+{x}_{0}\right)}^{1-{\alpha }_{x}}-{\left(x+{x}_{0}\right)}^{1-{\alpha }_{x}}}{{\left({x}_{2}+{x}_{0}\right)}^{1-{\alpha }_{x}}-{\left({x}_{1}+{x}_{0}\right)}^{1-{\alpha }_{x}}}\right).\end{array}\end{eqnarray} \tag{ 6 }$$

We show an example of a size distribution that contains the counts of events per bin n_cts(x) in Figure 1(a), the corresponding differential size distribution ${N}_{\mathrm{diff}}(x)={n}_{\mathrm{cts}}(x)/{\rm{\Delta }}x$ in Figure 1(b), and the corresponding cumulative size distribution N_cum(>x) in Figure 1(c). The event count histogram (Figure 1(a)) displays the inertial range (x₀, x₂), the minimum (x₁), and maximum value (x₂) of the size parameter x, with a peak in the event count histogram at x₀ (Figure 1(a)), which provides a suitable threshold definition, because incomplete sampling of small values $x\lt {x}_{0}$ is manifested by the drop of detected events on the left side of the peak x₀. This definition of a threshold x₀ has been proven to be very useful for characterizing data sets with incomplete or sensitivity-limited sampling (Aschwanden 2015). This provides also a definition for the inertial range (x₀, x₂), which we will quantify by the number of decades as a logarithmic ratio, i.e., $\mathrm{log}({x}_{2}/{x}_{0})$.

Zoom In Zoom Out Reset image size

An equivalent method to the cumulative size distribution is the rank-order plot. If the statistical sample is rather small, in the sense that no reasonable binning of a histogram can be done, either because we do not have multiple events per bin or because the number of bins is too small to represent a distribution function, we can create a rank-order plot. A rank-order plot is essentially an optimum adjustment to small statistics, by associating a single bin to every event. From an event list of a parameter ${x}_{i},i=1,\,\ldots ,\,{n}_{x}$, which is generally not sorted, we first have to generate a rank-ordered list by ordering the events according to increasing size,

$$\begin{eqnarray}&&{x}_{1}\leqslant {x}_{2}\,\leqslant \ldots \leqslant \,{x}_{j}\,\leqslant ...\leqslant {x}_{n},\quad j=1,\,\ldots ,\,n.\end{eqnarray} \tag{ 7 }$$

The bins are generally not equidistant, neither on a linear nor logarithmic scale, defined by the difference between subsequent values of the ordered x_j,

$$\begin{eqnarray}&&{\rm{\Delta }}{x}_{j}^{\mathrm{bin}}={x}_{j+1}^{\mathrm{bin}}-{x}_{j}^{\mathrm{bin}}.\end{eqnarray} \tag{ 8 }$$

In a rank-ordered sequence of n_x events, the probability for the largest value is $1/{n}_{x}$, for events that are larger than the second-largest event it is 2/n_x, and so forth, while events larger than the smallest event occur in this event list with a probability of unity. Thus, the cumulative frequency distribution is simply the reversed rank order,

$$\begin{eqnarray}&&{N}_{\mathrm{cum}}(\gt {x}_{j})=({n}_{x}+1-j),\qquad j=1,\,\ldots ,\,{n}_{x},\end{eqnarray} \tag{ 9 }$$

and the distribution varies from ${N}_{\mathrm{cum}}(\gt {x}_{1})={n}_{x}$ for j = 1 to ${N}_{\mathrm{cum}}(\gt {x}_{n})=1$ for $j={n}_{x}$. We can plot a cumulative frequency distribution with ${N}_{\mathrm{cum}}(\gt {x}_{j})$ on the y-axis versus the size x_j on the x-axis. The distribution is normalized to the number of events n_x,

$$\begin{eqnarray}&&{\int }_{{x}_{1}}^{{x}_{n}}N(x){dx}={N}_{\mathrm{cum}}(\gt {x}_{1})={n}_{x}.\end{eqnarray} \tag{ 10 }$$

We overlay rank-order plots in each of the cumulative size distributions in Figures 1–6 (shown with black diamond symbols), for the n_ee = 10 most extreme events in every size distribution. If there is a “Dragon-King” event, it should manifest itself by a significant offset from the best-fit cumulative size distribution at the upper end in the zone of most extreme events near x ≲ x₂, say by a factor of two at least. The example shown in Figure 1 exhibits a ratio of q_ee = 1.03 for the most extreme event, which is a very small offset from the best-fit cumulative size distribution that does not manifest a Dragon-King event.

In order to test the consistency of the size distribution of extreme events with the size distribution of all events, we perform least-square fits of the cumulative size distributions, such as shown in Figure 1(c). The cumulative distribution of an observed parameter x is shown in form of a black histogram with error bars N_cum(>x), while the best-fit cumulative size distribution is shown with a red curve (fitted over the range (x₀, x₂) above the threshold), and the extrapolated range below the threshold x₀ is shown with a dashed red curve. We see in Figure 1(c) that the reduced chi-square value is χ_all = 1.46 for all events, while the chi-square value is χ_ee = 0.56 in the range that contains the n_ee = 10 most extreme events, and thus is self-consistent between the two ranges, so that we can conclude that the extreme events do not show any significant deviation (χ ≳ 2) from the overall (best-fit) thresholded power-law distribution in this particular data set.

We estimate the goodness-of-fit from the uncertainties σ_cum(x_i) of the values ${N}_{i},i=1,\,\ldots ,\,{n}_{x}$ in the fitted size distribution with the standard reduced χ²-criterion,

$$\begin{eqnarray}&&{\chi }_{\mathrm{cum}}=\sqrt{\displaystyle \frac{1}{({n}_{i}-{n}_{\mathrm{par}})}\displaystyle \sum _{i=1}^{{n}_{x}}\displaystyle \frac{{\left[{N}_{\mathrm{cum},\mathrm{fit}}\left({x}_{i}\right)-{N}_{\mathrm{cum},\mathrm{obs}}\left({x}_{i}\right)\right]}^{2}}{{\sigma }_{\mathrm{cum},i}^{2}}},\end{eqnarray} \tag{ 11 }$$

where n_x is the number of bins, ${n}_{\mathrm{par}}=2$ is the number of free parameters (${n}_{0},a$) of the fitted size distribution, ${N}_{\mathrm{cum},\mathrm{obs}}(\gt {x}_{i})$ is the observed cumulative number of events in each bin, ${N}_{\mathrm{cum},\mathrm{fit}}(\gt {x}_{i})$ is the fitted cumulative number of events in each bin, and σ_cum(x_i) is the estimated uncertainty of the cumulative number of events,

$$\begin{eqnarray}&&{\sigma }_{\mathrm{cum},i}=\sqrt{{N}_{\mathrm{cum}}({x}_{i})-{N}_{\mathrm{cum}}({x}_{i+1})}.\end{eqnarray} \tag{ 12 }$$

This definition is slightly different from the estimate in the previous study (Aschwanden 2015), i.e., ${\sigma }_{\mathrm{cum},i}=\sqrt{({N}_{\mathrm{cum},i})}$, which represents a lower limit on the uncertainty only, since the counts of events in each bin are not independent of each other in a cumulative distribution function. Our new approach takes only independent events into account for the estimation of the uncertainty, which in a cumulative distribution is the increment $[{N}_{\mathrm{cum}}(\gt {x}_{i})-{N}_{\mathrm{cum}}(\gt {x}_{i+1})]$, while the remainder of the counts in the bins ${N}_{\mathrm{cum}}(\gt {x}_{i+1}),\,\ldots ,\,{N}_{\mathrm{cum}}(\gt {x}_{{n}_{x}})$ do not vary in the cumulative fit of bin ${N}_{\mathrm{cum}}(\gt {x}_{i})$, and therefore do not contribute to the uncertainty $\sigma ({x}_{i})$.

Besides the fits of the cumulative size distributions (Figure 1(c)), we also perform fits to the differential size distributions for consistency tests (Figure 1(b)). Since the differential size distribution is quantified in terms of counts per bin width, ${N}_{\mathrm{diff}}(x)={n}_{\mathrm{cts}}(x)/{\rm{\Delta }}x$, the corresponding χ²-criterion is,

$$\begin{eqnarray}&&{\chi }_{\mathrm{diff}}=\sqrt{\displaystyle \frac{1}{({n}_{x}-{n}_{\mathrm{par}})}\displaystyle \sum _{i=1}^{{n}_{x}}\displaystyle \frac{{\left[{N}_{\mathrm{fit},\mathrm{diff}}\left({x}_{i}\right)-{N}_{\mathrm{sim},\mathrm{diff}}\left({x}_{i}\right)\right]}^{2}}{{\sigma }_{\mathrm{diff},i}^{2}}},\end{eqnarray} \tag{ 13 }$$

and the uncertainty ${\sigma }_{\mathrm{diff},i}$ in terms of Poisson statistics is,

$$\begin{eqnarray}&&{\sigma }_{\mathrm{diff},i}=\sqrt{{n}_{\mathrm{cts}}({x}_{i})}/{\rm{\Delta }}{x}_{i}=\sqrt{{N}_{\mathrm{diff}}({x}_{i}){\rm{\Delta }}{x}_{i})}/{\rm{\Delta }}{x}_{i}.\end{eqnarray} \tag{ 14 }$$

An example of such a fit is shown in Figure 1(b). Note that the size distribution of the differential size distribution shows a straight power-law function even in the range of extreme events, while the cumulative size distribution shows an exponential-like drop-off towards the zero value at x₂, as a consequence of the integral function defined in Equation (6). We can compare now the results of the two power-law fits and find a power-law slope of α_diff,all = 1.75 ± 0.02 for the differential size distribution (Figure 1(b)), while the the power-law slope is ${\alpha }_{\mathrm{cum},\mathrm{all}}=1.73\pm 0.02$ for the cumulative size distribution (Figure 1(c)). Now, if we inspect the χ²-values in the range of the ten most extreme events, we find χ_cum,all = 1.46 for all events (above the threshold x₀), and a value of χ_cum,ee = 0.56 for the 10 most extreme events (Figure 1(c)). Therefore, not only the two methods of differential and cumulative distributions are self-consistent, but the two different ranges (of small-to-large and extreme events) lead to a self-consistent size distribution function also, so that we can conclude that extreme events belong to the same population as all other events, and that no significant deviation from the thresholded power-law distribution function is found for this data set, which was obtained from hard X-ray peak photon count rates in solar flares, using the Hard X-ray Burst Spectrometer (HXRBS) onboard the Solar Maximum Mission (SMM). The first size distributions of these solar flare data were published in Dennis (1985) and Crosby et al. (1993).

Alternatively to the chosen method of least-square fitting of differential or cumulative size distributions, other methods have been used to test the consistency of an ideal power-law function with observed size distributions, such as the maximum likelihood estimator, the Bayesian information criterion, the Anderson–Darling test, or the Kolmogorov–Smirnov statistic (Clauset et al. 2009). However, since no accommodation for (i) undersampling of small events, (ii) subtraction of event-unrelated background, or (iii) automated determination of a threshold is provided in those alternative methods, they are not suitable for modeling of astrophysical data (Aschwanden 2015).

The results of the data analysis are presented in Figures 2–6 and in Table 1. While we used the generic term “size” in the definition of occurrence frequency distributions, we should be aware that the size can be represented by any measure, magnitude, or physical quantity that can be measured in a set of events, but we should not expect to find the same power-law slope for different quantities, because there is often some nonlinear scaling involved. For the solar observations analyzed here we sample the following quantities: peak count rates of the hard X-ray photon fluxes of flares (P) (observed with the instruments HXRBS and the Burst and Transient Source Experiment (BATSE) on the Compton Gamma Ray Observatory (CGRO), and the Ramaty High Energy Solar Spectroscopic Imager (RHESSI) spacecraft), total counts (or fluences) of hard X-ray fluxes (C), flare durations (T), emission measure (EM) in soft X-rays and extreme ultra-violet (EUV), thermal flare energies (E_th), nonthermal flare energies (E_nth), dissipated magnetic energies (E_mag), flare volumes (V), flare areas (A), and flare length scales (L). The energetic parameters have been calculated in recent statistical studies on the global energetics of solar flares (Aschwanden et al. 2014, 2015, 2016).

Zoom In Zoom Out Reset image size

Table 1.  Results of Forward-fitting of Differential and Cumulative Size Distributions for Different Statistical

	BG	N	NT	IR	x0	αdiff	αcum	χdiff	χcum	χee	qee
Solar flares:
Peak HXRBS)	58	10856	7579	4.0	19	1.75 ± 0.02	1.73 ± 0.02	0.96	1.46	0.56	1.03
Peak HXRBS	41	6223	3979	3.8	31	1.69 ± 0.03	1.70 ± 0.02	1.18	1.86	1.12	1.09
Peak BATSE	0	7234	4996	3.3	501	1.91 ± 0.03	1.81 ± 0.02	3.18	4.39	0.78	1.11
Peak RHESSI	0	7996	4399	3.7	20	1.87 ± 0.03	1.86 ± 0.02	1.28	2.00	0.77	1.88
Peak RHESSI (A16)	0	290	149	1.7	822	2.44 ± 0.19	2.28 ± 0.13	0.70	1.07	0.81	1.45
Counts HXRBS	2	8129	3314	4.6	3100	1.72 ± 0.03	1.66 ± 0.02	1.54	2.68	0.65	1.40
Counts BATSE	0	2748	2748	4.8	2000	1.48 ± 0.03	1.50 ± 0.03	1.67	5.20	2.94	1.07
Counts RHESSI	0	11548	8081	4.9	4783	1.60 ± 0.02	1.64 ± 0.02	3.34	10.98	2.77	1.54
Counts RHESSI (A16)	0	289	147	1.8	$1.42\times {10}^{6}$	2.40 ± 0.19	2.00 ± 0.12	1.17	1.95	1.95	1.63
Duration HXRBS	2	11539	4341	2.2	126	2.61 ± 0.04	2.66 ± 0.02	4.43	5.45	1.58	0.90
Duration BATSE	0	7241	4562	2.1	37	2.32 ± 0.03	2.19 ± 0.03	5.65	5.81	0.00	1.77
Duration RHESSI	0	11445	8296	2.2	25	2.05 ± 0.02	1.90 ± 0.02	7.59	8.28	10.19	0.93
Duration (A14)	0	171	123	1.3	0	2.76 ± 0.24	2.73 ± 0.21	0.58	0.83	0.56	1.15
Duration (A16)	0	289	189	1.0	0	3.22 ± 0.22	3.33 ± 0.20	1.53	1.49	0.64	1.83
Length (A14)	0	171	112	0.9	26	5.13 ± 0.44	5.03 ± 0.38	1.10	0.52	0.55	1.25
Length (A15)	0	389	210	0.8	8	4.23 ± 0.27	4.45 ± 0.23	1.18	1.00	1.04	0.95
Area (A14)	0	171	100	1.5	899	2.84 ± 0.27	2.82 ± 0.22	0.71	1.05	0.73	1.75
Volume (A14)	0	171	96	2.2	30300	2.32 ± 0.22	2.09 ± 0.16	0.94	1.34	0.68	3.08
Volume (A15)	0	389	239	2.1	640	1.83 ± 0.12	1.86 ± 0.09	0.56	0.72	0.55	1.07
Magnetic energy (A14)	0	171	68	1.1	119	3.04 ± 0.34	2.83 ± 0.22	1.11	0.92	0.65	0.92
Nonthermal energy (A16)	0	192	89	2.6	0.32	2.50 ± 0.25	1.59 ± 0.11	1.28	2.34	2.08	1.87
Thermal energy (A15)	0	390	279	1.8	3	2.15 ± 0.13	2.16 ± 0.11	0.86	0.73	0.73	1.02
Emission measure (A15)	0	391	251	1.4	$6.5\times {10}^{8}$	2.97 ± 0.19	2.82 ± 0.14	0.86	0.58	0.49	1.00
Stellar flares:
Peak KEPLER	0	208	124	2.8	60255	1.66 ± 0.14	1.68 ± 0.12	0.81	1.26	1.45	3.55
Bolom. KEPLER	0	1537	727	1.5	3.8 × 1034	2.56 ± 0.09	2.55 ± 0.07	1.31	2.18	0.67	0.97
Terrestrial data:
Earthquakes	1	17425	9423	4.6	100	1.80 ± 0.02	1.85 ± 0.01	2.22	3.11	1.45	2.82
Fires	1	55853	38190	5.3	2	1.56 ± 0.01	1.53 ± 0.01	13.90	30.65	30.93	0.92
Cities	0	19445	12784	4.0	501	1.79 ± 0.02	1.79 ± 0.01	5.69	18.37	14.18	1.79
Blackouts	0	209	134	1.7	50000	1.98 ± 0.16	1.90 ± 0.13	0.80	0.83	0.63	3.35
Terrorism	1	2699	1949	2.4	2	2.48 ± 0.05	2.28 ± 0.04	4.86	2.25	0.43	6.15
Words	1	6609	4980	3.6	1	2.15 ± 0.03	1.99 ± 0.02	6.31	3.33	1.09	1.99
Surnames	2	2280	1663	2.2	15700	2.55 ± 0.06	2.50 ± 0.05	1.10	1.89	1.34	1.13
Weblinks	1	$9.4\times {10}^{7}$	$5.9\times {10}^{7}$	2.4	4	2.47 ± 0.03	2.39 ± 0.01	4.39	223	0.67	1.07

Download table as:  ASCIITypeset image

The results of forward-fitting differential and cumulative size distributions of different events are summarized in Table 1, which includes: background counts (BG), the number of all events (N), the number of events above threshold (NT), the logarithm of the inertial range $\mathrm{log}({x}_{2}/{x}_{1})$ (IR), the threshold value (x₀), the power-law slope from fitting differential (α_diff) and cumulative occurrence frequency distributions (α_cum), the χ²-values of differential (χ_diff) and cumulative size distributions (χ_cum), and the mean χ²-values in the range of the 10 most extreme events (χ_ee).

We determined the background level by minimizing the goodness-of-fit χ²-criterion as a function of the background level (see Figure 7 in Aschwanden et al. 2015). It turned out that only the HXRBS data needed background subtraction (by 41 and 58 cts s⁻¹), while the BATSE and RHESSI flare catalogs provide background-subtracted flux values.

The total number of events (N in Table 1) encompasses all sampled events, but only those (NT) above the threshold x₀ could be used in the fitting of size distributions. The inertial range is characterized by the logarithmic ratio $\mathrm{log}({x}_{2}/{x}_{0})$, which ranges from 2 to 5 decades. Those with the largest inertial ranges and largest number of events provide the most accurate fits of size distributions.

A first important test is the self-consistency between the two different methods, i.e., the fitting of the differential and cumulative size distributions. We can simply compare the resulting power-law slopes, which are labeled as α_diff and α_cum in Table 1. We find a satisfactory agreement between the two methods within the derived uncertainties of a few percent, as it can be seen for all data sets with large statistics (Table 1: $N\approx {10}^{3}\mbox{--}{10}^{4}$ events). This corroborates the self-consistency between the two methods. The high accuracy is mostly achieved by modeling the size distribution with a thresholded power-law distribution function (Equations (3) and (6)), which is more suitable to represent the data than an ideal power-law function.

A second important test is the comparison of the goodness-of-fit between the differential and cumulative size distributions, labeled as χ_diff and χ_cum in Table 1. Among the nine fitted size distributions of three solar observables (P, C, T) with three different instruments each (HXRBS, BATSE, RHESSI; shown in Figure 2), we find that only one data set is consistent with the fitted (differential and cumulative) size distributions, namely P of HXRBS, with χ_diff ≈ 2 and χ_cum ≲ 2, while the other eight data sets exhibit significant mismatches χ_diff ≈ χ_cum ≈ 3–13. This indicates that the large-number statistics of these data sets is sensitive to apparently little but significant deviations from power laws. These deviations can clearly be recognized in the poor values χ_diff,all and χ_cum,all in the plots of Figures 2(b)–(i). There is a tendency that peak count (P) size distributions match a power-law function closest, while the total counts (C) and the flare durations (T) deviate significantly. The reason for the latter property could be due to a violation of the principle of timescale separation that is required in self-organized criticality (SOC) models (Aschwanden et al. 2016), caused by confusion between overlapping short and long-duration flares.

The flare observables (P, C, T) are easiest to measure and have often been used in modeling of SOC models. However, since these observables may not be a good proxy for the representation of the energy contained in SOC avalanches, we also inspect the size distribution of spatial and flare parameters, such as ($L,A,V,{E}_{\mathrm{mag}},{E}_{\mathrm{th}},{E}_{\mathrm{nth}},\mathrm{EM}$), which have been determined under the scope of global energetics in solar flares in a series of ongoing studies (Aschwanden et al. 2014, 2015, 2016). The cumulative size distributions are presented in Figures 3 and 4. Since the number of events for which energetic parameters have been determined includes M and X Geostationary Operational Environmental Satellites (GOES) class flares only, we deal with small-number statistics in the order of $N\approx 100\mbox{--}300$ events per data set. Consequently, the resulting fits yield acceptable values for the goodness-of-fit, i.e., χ² ≲ 2, but these small data sets have insufficient sensitivity to measure deviations from power-law functions. Nevertheless, since they contain large flares only, they are suitable to detect the outliers of extreme events.

Zoom In Zoom Out Reset image size

The only obvious deviation from a power law is detected for the size distribution of nonthermal energies (Figure 4(b)), which is suspected to arise due to an instrumental irregularity, such as pulse pile-up in RHESSI data.

Zoom In Zoom Out Reset image size

We now investigate the so-called Dragon-King hypothesis (Sornette 2009; Sornette & Ouillon 2012). This hypothesis suggests that the most extreme events may belong to a different population of event size distributions than smaller or even large events. In this study we define three different parameters to characterize Dragon-King events: (i) a range of the largest extreme events in a data set, which we set arbitrarily to n_ee = 10 here; (ii) the ratio of the observed size ${x}_{\max ,\mathrm{obs}}={x}_{2}$ to the expected size ${x}_{\max ,\mathrm{fit}}$ of the most extreme event, based on the cumulative size distribution fit, i.e., ${q}_{\mathrm{ee}}={x}_{\max ,\mathrm{obs}}/{x}_{\max ,\mathrm{fit}}$, (which we may call the extreme event excess factor); and (iii) the goodness-of-fit χ_ee of the cumulative size distribution fit in the range of the largest n_ee = 10 events. These three values are given in the lower right corner of each panel in Figures 1–6 and Table 1.

As an example we consider Figure 1, where the goodness-of-fit in the extreme event part of the size distribution is i.e., χ_ee = 0.56, which compares well with the value obtained from fitting all events above the threshold, i.e., χ_all = 1.46, which does not yield any evidence for a different population. A second test is the extreme event excess factor, which for this event is q_ee = 1.03 and yields no evidence for a Dragon-King event either (Figure 1(c)).

We list all evaluated extreme event excess factors q_ee in Table 1 (also given in the panels of Figures 1–6). Among the solar flare parameters, all those factors vary in the range of q_ee ≈ 0.9–1.9, with a single outlier of q_ee = 3.08 (for one flare volume data set, see Figure 3(c)), out of 23 solar data sets. Therefore, we see no indication that solar extreme events belong to different populations of event distributions.

The estimation of an unknown tail distribution function can be derived with tools from extreme value theory, such as with the Pickands–Balkema–de Haan theorem (Balkema & de Haan 1974; Pickands 1975), the Fisher–Tippett–Gnedenko theorem (Fisher 1930), or as described in Gumbel (1958), which converges to three possible distributions (the Gumbel distribution, the Fréchet distribution, or the Weibull distribution). In our case, however, the cumulative form of the thresholded power-law size distribution (generalized Pareto) function is our preferred model of choice, due to its ability of taking incomplete sampling, background subtraction, and finite system size into account, and thus is known a priori, so that it does not need to be estimated with extreme value theory, as it would be needed in the case of an unknown tail distribution function. The theoretically defined distribution is then fitted to the observed distributions directly, and the extreme event excess factor q_ee is conveniently quantified in terms of standard deviations in common χ²-statistics.

Impulsive flaring with rapid increases in the brightness in UV or EUV has been observed for a number of so-called flare stars, such as AD Leo, AB Dor, YZ Cmi, EK Dra, or Eri. These types of stars include cool M dwarfs, brown dwarfs, A-type stars, giants, and binaries in the Hertzsprung–Russell diagram. Most of these stars are believed to have hot, soft X-ray emitting coronae, similar to our Sun (a G5 star), and thus magnetic reconnection processes are believed to operate in a similar way as on our Sun (see, e.g., review by Güdel 2004).

We complement the solar flare data sets with three data sets of stellar flares, two recorded with KEPLER and one observed with the EUVE telescope, with the cumulative size distributions shown in Figures 4(g), (h) and 5. The fit of a cumulative power-law distribution shown in Figure 4(g) produces a size distribution with a power-law slope of α = 1.68 ± 0.12 for the peak counts, which is similar to the slope observed in solar flares, i.e., α = 1.70 ± 0.02 with HXRBS (Figure 2(a)), α = 1.81 ± 0.02 with BATSE (Figure 2(b)), and α = 1.86 ± 0.02 with RHESSI (Figure 2(c)). The goodness-of-fit for extreme events is χ_ee = 1.45, which is similar to χ_all = 1.26 for all events (Figure 4(g)).

A second data set with stellar flare data from KEPLER, derived from the bolometric energy, is shown in Figure 4(h), yielding a power-law slope of ${\alpha }_{\mathrm{bol}}=2.55\pm 0.06$, which is similar to the value expected for flare areas ${\alpha }_{A}=2.33$ (Table 2). Since the largest contributions to the bolometric intensity come from the solar or stellar photospheric surface area, it makes sense that its size distribution scales with the flare area (rather than with a 3D volume).

Table 2.  Theoretically Predicted Power-law Slopes of Solar Flare Size Distributions

Parameter	Power-law Exponent	Power-law Exponent	Power-law Exponent	Data Set
	(General Expression)	$S=3,\beta =1,\gamma =1$,	of Observed Solar Flares
		${D}_{2}=3/2,{D}_{3}=2$	αx
Length L	${\alpha }_{L}=S$	${\alpha }_{L}=3.00$	${\alpha }_{L}=(4.45\pm 0.23)$	A15
Length L	${\alpha }_{L}=S$	${\alpha }_{L}=3.00$	${\alpha }_{L}=(5.03\pm 0.39)$	A14
Area A	${\alpha }_{A}=1+(S-1)/{D}_{2}$	${\alpha }_{A}=2.33$	${\alpha }_{A}=(2.82\pm 0.22)$	A14
Magnetic energy Emag	${\alpha }_{\mathrm{mag}}\approx {\alpha }_{A}$	${\alpha }_{A}=2.33$	${\alpha }_{\mathrm{mag}}=(2.83\pm 0.22$	A14, A16
Emission measure EM	${\alpha }_{\mathrm{EM}}\approx {\alpha }_{A}$	${\alpha }_{A}=2.33$	${\alpha }_{\mathrm{EM}}=(2.82\pm 0.14)$	A15
Bolom. energy Ebol	${\alpha }_{\mathrm{bol}}\approx {\alpha }_{A}$	${\alpha }_{A}=2.33$	${\alpha }_{\mathrm{bol}}=2.55\pm 0.07$	KEPLER
Volume V	${\alpha }_{V}=1+(S-1)/{D}_{3}$	${\alpha }_{V}=2.00$	${\alpha }_{V}=(2.09\pm 0.16)$	A14
Volume V	${\alpha }_{V}=1+(S-1)/{D}_{3}$	${\alpha }_{V}=2.00$	${\alpha }_{V}=(1.86\pm 0.09)$	A15
Thermal energy Eth	${\alpha }_{{E}_{\mathrm{th}}}\approx {\alpha }_{V}$	${\alpha }_{V}=2.00$	${\alpha }_{\mathrm{th}}=(2.16\pm 0.11)$	A15, A16
Time duration T	αT = 1 + (S − 1)β/2	αT = 2.00	αT = 2.66 ± 0.03	HXRBS
Time duration T	αT = 1 + (S − 1)β/2	αT = 2.00	αT = 2.19 ± 0.03	BATSE
Time duration T	αT = 1 + (S − 1)β/2	αT = 2.00	αT = 1.90 ± 0.02	RHESSI
Time duration T	αT = 1 + (S − 1)β/2	αT = 2.00	αT = (2.73 ± 0.21)	A14, A16
Peak flux P	αP = 1 + (S − 1)/(γS)	${\alpha }_{P}=1.67$	${\alpha }_{P}=1.70\pm 0.02$	HXRBS
Peak flux P	αP = 1 + (S − 1)/(γS)	${\alpha }_{P}=1.67$	${\alpha }_{P}=1.81\pm 0.02$	BATSE
Peak flux P	αP = 1 + (S − 1)/(γS)	${\alpha }_{P}=1.67$	${\alpha }_{P}=1.86\pm 0.02$	RHESSI
Peak flux P	αP = 1 + (S − 1)/(γS)	${\alpha }_{P}=1.67$	${\alpha }_{P}=(2.28\pm 0.13)$	A16
Peak flux P	αP = 1 + (S − 1)/(γS)	${\alpha }_{P}=1.67$	${\alpha }_{P}=(1.68\pm 0.12)$	KEPLER
Fluence C	${\alpha }_{C}=1+(S-1)/(\gamma {D}_{S}+2/\beta )$	${\alpha }_{C}=1.50$	${\alpha }_{C}=1.67\pm 0.02$	HXRBS
Fluence C	${\alpha }_{C}=1+(S-1)/(\gamma {D}_{S}+2/\beta )$	${\alpha }_{C}=1.50$	${\alpha }_{C}=1.50\pm 0.03$	BATSE
Fluence C	${\alpha }_{C}=1+(S-1)/(\gamma {D}_{S}+2/\beta )$	${\alpha }_{C}=1.50$	${\alpha }_{C}=1.64\pm 0.02$	RHESSI
Fluence C	${\alpha }_{C}=1+(S-1)/(\gamma {D}_{S}+2/\beta )$	${\alpha }_{C}=1.50$	${\alpha }_{C}=(2.00\pm 0.12)$	A16
Nonthermal energy Enth	${\alpha }_{{E}_{\mathrm{nth}}}\approx {\alpha }_{C}$	${\alpha }_{C}=1.50$	${\alpha }_{{E}_{\mathrm{nth}}}=(1.59\pm 0.12)$	A16

Note. The power-law slopes fitted to observed data sets (HXRBS, BATSE, RHESSI) are listed in the fourth column. Measurements with poor statistics and less reliable values are marked with parenthesis (...), obtained from the data sets A14, A15, and A16 (Aschwanden et al. 2014, 2015, 2016).

Download table as:  ASCIITypeset image

In Figure 5 we present cumulative size distribution fits of 12 cool (type F to M) stars, observed by Audard et al. (2000), with typically 5–15 flare events per star. For most of these stars, the flare intensities closely follow the fitted cumulative distribution all the way to the most extreme event, except for two cases (out of 12): HD 2726 (F2 V) in the top left panel and CN Leo, 1994 (M6 V) in the bottom middle panel, which could qualify as outliers (or Dragon-King events). The power-law slopes inferred in this sample vary strongly within a range of $\alpha \approx 2.0\pm 0.4$, which tend to be somewhat higher than the mean values of solar flare size distributions, but are not reliable in such small-number statistics.

Zoom In Zoom Out Reset image size

In summary, stellar flares fit cumulative size distributions at the upper cutoff (within an energy range of ≈0.8–1.4 decades), but the power-law slopes appear to be slightly steeper than observed in solar flares, and outliers of extreme events occur only rarely (in 2 out of 12 cases).

While we are interested in the investigation of deviations from ideal power laws, which is a second-order effect, we should consider first the physical justification of the existence of power-law size distributions in the first place. The original discoverers of SOC models noted power-law–like size distributions with power-law slopes of ${\alpha }_{x}\approx 2$ of the sizes x of SOC avalanches, as well as power-law slopes of α_T ≈ 2 for the durations T of SOC avalanches, obtained in cellular automaton simulations (e.g., Pruessner 2012). A more physical approach of modeling and predicting power-law slopes of size distributions observed from astrophysical SOC systems has been put forward in terms of the so-called fractal-diffusive (FD-SOC) model (Aschwanden 2013, 2015; Aschwanden et al. 2016). We juxtapose theoretical predictions of this model and mean observational values in Table 2.

The most fundamental assumption in the FD-SOC model is the scale-free probability conjecture, which directly predicts a power-law size distribution function of geometric length scales L from first principles, i.e.,

$$\begin{eqnarray}&&N(L){dL}={L}^{-S}\ {dL},\end{eqnarray} \tag{ 15 }$$

where S is defined as the Euclidean dimension of the SOC system (with possible values of $S=1,2,3$). From this fundamental power-law distribution, other power-law size distributions of any physical size parameter x that has a linear (x ∝ L), or a nonlinear relationship $x(L)\propto {L}^{p}$ to the length scale L (with power index p), can be derived. The following scaling laws have been used in the calculation of the power-law indices listed in Table 2:

$$\begin{eqnarray}&&A\propto {L}^{{D}_{2}}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&V\propto {L}^{{D}_{3}}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&L\propto {T}^{\beta /2}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&F\propto {V}^{\gamma }\propto {L}^{{D}_{S}\gamma }\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&P={F}_{\max }\propto {V}_{\max }^{\gamma }\propto {L}^{S\gamma }\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&E=F\ T\propto {L}^{(\gamma {D}_{S}+2/\beta )}\end{eqnarray} \tag{ 21 }$$

where D₂ is the area (A) fractal dimension, D₃ is the volume (V) fractal dimension, β is the spatiotemporal diffusion coefficient for the time duration (T) of a SOC avalanche with classical diffusion (β = 1), γ is the flux-volume scaling of observed fluxes F of astrophysical sources, V_max is the maximum fractal volume, which obeys the Euclidean dimension S of the SOC system, $P={F}_{\max }$ is the peak or maximum flux, and E = F T is the fluence, defined as the product of the flux F and time duration T.

The mean fractal dimension is estimated from the average of the minimum (${D}_{S,\min }\approx 1$) and maximum ${D}_{S,\max }=S$ values,

$$\begin{eqnarray}&&{D}_{S}\approx \displaystyle \frac{({D}_{S,\min }+{D}_{S,\max })}{2}=\displaystyle \frac{(1+S)}{2},\end{eqnarray} \tag{ 22 }$$

which yields ${D}_{2}=(1+2)/2=3/2$ and ${D}_{3}=(1+3)/2=\,2$ for SOC models that occupy either an Euclidean 2D or 3D dimensional space.

Since all relationships given in Equations (15)–(21) are expressed in terms of the variable L, the power-law indices α_x for the size distribution of the parameters $x\,=A,V,T,F,P,E$ can be derived by substitution of variables, i.e., $N(x){dx}=\ N[x(L)]\ | {dx}/{dL}]\ {dL}\propto {x}^{-{\alpha }_{x}}$, yielding the power-law indices α_x,

$$\begin{eqnarray}&&{\alpha }_{L}=S\approx 3\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{\alpha }_{A}=1+(S-1)/{D}_{2}\approx 7/3\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{\alpha }_{V}=1+(S-1)/{D}_{3}\approx 2\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{\alpha }_{T}=1+(S-1)\beta /2\approx 2\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\alpha }_{F}=1+(S-1)/(\gamma \ {D}_{S})\approx 2\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\alpha }_{P}=1+(S-1)/(\gamma \ S)\approx 5/3\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\alpha }_{E}=1+(S-1)/(\gamma \ {D}_{S}+2/\beta )\approx 3/2,\end{eqnarray} \tag{ 29 }$$

where the approximative numerical values are obtained by inserting $S=3,\beta =1,\gamma =1,{D}_{2}=1.5,{D}_{3}=2$.

In Table 2 we juxtapose these theoretical power-law slope values with the observed values. We see that there is mostly a good agreement (within a few percent) for solar flare size distributions with large statistics ($N\approx {10}^{3}\mbox{--}{10}^{4}$), such as for peak fluxes (α_P), total counts (α_C), and flare durations (α_T), while others with small-number statistics show larger deviations from the predictions, such as for length scales (α_L), areas (α_A), and volumes (α_V).

New results are obtained here from the power-law slopes of energetic flare parameters. The dissipated magnetic energies E_mag are found to have power-law indices that are close to the values of flare areas, i.e., ${\alpha }_{\mathrm{mag}}=2.8\pm 0.2$ versus ${\alpha }_{A}\approx 2.3$, which could be explained by the fact that the dissipation of free magnetic energies is concentrated in the chromospheric flare ribbons, which have a horizontal area-like surface geometry. This is also true for the bolometric flare energies observed in stellar flares, i.e., ${\alpha }_{\mathrm{bol}}=2.6\pm 1.0$ versus ${\alpha }_{A}\approx 2.3$, since the bolometric brightness is mostly irradiated in the photosphere of stellar surfaces. Also for the thermal EM, we find similarity with area-like features, i.e., ${\alpha }_{\mathrm{EM}}=2.8\pm 0.2$ versus ${\alpha }_{A}\approx 2.3$, which suggests that the soft X-ray EM in solar flares has an area-like extension with relatively low altitudes. For thermal flare energies, for which we expect mostly a volume dependence,

$$\begin{eqnarray}&&{E}_{\mathrm{th}}=2{n}_{e}{k}_{B}{T}_{e}V\propto V\propto {L}^{{D}_{3}},\end{eqnarray} \tag{ 30 }$$

as the electron density and the electron temperature vary much less than the flare volume. Indeed, we find values of ${\alpha }_{{E}_{\mathrm{th}}}=2.2\pm 0.1$ that are similar to the volume power-law index, ${\alpha }_{V}\approx 2$. For nonthermal flare energies we find a close coincidence of the power-law slope of ${\alpha }_{E,\mathrm{nth}}=1.57\pm 0.11$, which is close to the power-law slopes of fluences, ${\alpha }_{C}\approx 1.50$, which is expected because both the fluences (E = F T) and the nonthermal energies, ${E}_{\mathrm{nth}}=\int {e}_{\mathrm{nth}}(t)\ {dt}$, are time-integrated flare quantities.

These comparisons of power-law indices thus reveal interesting relationships that are relevant for determining the physical scaling laws that govern solar flares.

From the 23 solar (Figures 3, 4) and 14 stellar data sets (Figures 4, 5) we found little evidence for outliers among the most extreme events in each data set. The largest outliers for the most extreme events are found for a small data set (N = 172) of flare volumes with an excessive extreme event factor of q_ee = 2.75 (Figure 3(c)), for a small stellar flare data set (N = 209) with q_ee = 3.55 (Figure 4(g)), and for the two stars HD 2726 and CN Leo with even less statistics (N ≈ 15; Figure 5). The fact that all these most extreme events occur in small data sets, while larger data sets with N ≳ 10³ − 10⁴ events (Figure 2) reveal no excesses larger than q_ee ≲ 2.0, provides little evidence for the existence of super-extreme events that are not part of the “normal” population of small and large flare events. Hence, we conclude that extreme events of solar and stellar flares are largely scale-free (over an inertial range of ≈3–4 decades) in their volumes, fluxes, or energies.

What does it mean that we find no significant outliers in the cumulative size distributions of solar and stellar flares? As the size distribution of length scales reveals in Figure 3(a), solar flares have length scales of $L\approx 10\mbox{--}200$ Mm, which at the upper limit matches the maximum size of active regions on the solar surface. The largest possible flare that extends over the entire solar surface would have a circumference of $2\pi \times 696$ Mm, which is apparently impossible since we never observed a flare with a size larger than L ≲ 200 Mm. One effect that could help to increase the maximum flare size is the case of sympathetic flaring, which involves a magnetic coupling between adjacent (or even remote) active regions. However, the larger an unstable magnetic field configuration is, the smaller the filling factor of magnetic energy and heated plasma is with respect to the size of a flaring region. Consequently, the envisioned process of amplification and synchronization promoted in the Dragon-King hypothesis (Sornette 2009; Sornette & Ouillon 2012) seems to be questionable for solar and stellar flares.

It has been argued that real data deviate from ideal power-law size distributions, as they have been reproduced by cellular automaton simulations (e.g., Pruessner 2012). Some skeptics went even as far as to deny the existence of power laws at all (Stumpf & Porter 2012). Traditional least-square fitting techniques have been criticized to be inaccurate, while maximum-likelihood fitting methods with goodness-of-fit tests based on the Kolmogorov–Smirnov statistic were preferred (Clauset et al. 2009). However, we have demonstrated here and earlier (Aschwanden 2015) that least-square fitting methods of both differential and cumulative size distributions yield self-consistent fits of the power-law slope, as long as data truncation, undersampling of data in small events below some threshold, subtraction of event-unrelated background noise, and the steep fall-off at the upper bound of the cumulative size distribution is properly modeled. In other words, a straight power-law function represents an oversimplified model for most cases, while inclusion of the effects enumerated above is essential in defining adequate models for fitting statistical distributions.

Now, with the availability of large data sets (with $N\approx {10}^{5}\mbox{--}{10}^{8}$ events), we are becoming much more sensitive to the smallest deviations from ideal power laws. So it is no surprise that deviations from power laws become the rule rather than the exception. The deeper reason for these power-law deviations is the inhomogeneity of the medium and the physical conditions. It is instructive to consider the sandpile analogy of SOC systems. Although the cone shape of a sandpile appears to have a constant slope or angle of repose, when viewed from distance, closer inspection reveals macroscopic channels, erosions, bumps, and dents, on top of the microscopic fine structure. As an illustration we mention some empirical data sets (chosen by Clauset et al. 2009) that are shown in Figure 6 and described in the Appendix. Some data sets clearly reveal a complete inadequacy to fit any power-law function to data with a broken power law (Figure 6(b): wildfires; Figure 6(c): cities). The case of weblinks has such a tremendous statistic ($N=1.35\times {10}^{8}$ events) that it is highly sensitive to the smallest deformations from an ideal power law, i.e., amounting to χ_all = 223 for a slight deficit near $x\approx {10}^{5}\mbox{--}{10}^{6}$ (Figure 6(h)). A Dragon-King event with an excess of ${q}_{\mathrm{ee}}\approx 5$ is indicated in the data set of terrorism (Figure 6(e)).

Significant deviations from ideal power-law size distributions have also been detected in stellar flares, e.g., in the Kepler-based study of Sheikh et al. (2016), and most of the analyzed stars were found to exhibit no evidence for SOC avalanching. It is argued that these stellar flares are subject to “tuned criticality” rather than “SOC” (Sheikh et al. 2016).

The observations of SOC avalanches exhibit power-law–like size distributions, which constrain (scale-free) scaling laws of physical parameters (see Section 3.1), but do not directly reveal the physical mechanism of SOC processes as such. The power-law slopes of the size distributions of solar and flare data have been found to be consistent with the scale-free probability conjecture (Equation (15)), fractal-diffusive transport (Equation (18)), and with EUV and soft X-ray emission (flux or intensity) that is quasiproportional to the fractal volume of an avalanche event (Equation (20)). These scaling laws are consistent with standard solar flare models, which are driven by a magnetic reconnection process in the solar corona (Lu & Hamilton 1991). However, the observed scaling laws are also consistent with other physical mechanisms, such as turbulent photospheric convection, which may be (stochastically) coupled with coronal reconnection events (Uritsky & Davila 2012; Uritsky et al. 2013; Knizhnik et al. 2018). Additional interpretations of physical mechanisms operating in SOC systems can be found in Aschwanden (2011) and Aschwanden et al. (2016).

The interpretation of SOC systems in terms of physical mechanisms is beyond the scope of this study, but the fact that we find virtually no outliers or Dragon-King events implies that only one single physical mechanism with appropriate properties is needed to explain the observed power-law size distributions. A dual energy dissipation system is likely to produce broken power-law distributions. However, the same physical process, such as a magnetic reconnection process, can have multiple secondary energy dissipation processes, including nonthermal particle acceleration (hard X-ray fluxes), thermal flare plasma heating (soft X-ray fluxes), launching of a coronal mass ejection (white light emission), or electron beam formation (radio emission). Quantitative energy budget calculations have revealed energy closure for the primary flare-dissipated magnetic energy and secondary energy dissipation mechanisms (Aschwanden et al. 2017). All energy that is dissipated during a solar flare is supplied by the free (magnetic) energy in an active region. Interestingly, multiscale intermittent dissipation with power-law–like SOC properties has also been detected in quiet sun regions at spatial scales of ≳3 Mm, controlled by turbulent photospheric convection (Uritsky et al. 2013). It appears that we observe a coexistence of coronal SOC systems (in active regions) and photospheric SOC systems (in quiet sun regions).