A Statistical Search for a Uniform Trigger Threshold in Solar Flares from Individual Active Regions

The SDP process is a doubly stochastic renewal process that models the "stress" in the system as a function of time, X(t), as

$$\begin{eqnarray}&&X(t)=X(0)+t-\displaystyle \sum _{i=0}^{N(t)}{\rm{\Delta }}{X}_{i},\end{eqnarray} \tag{ 1 }$$

where X and t are expressed respectively in dimensionless units of X_cr (the critical stress in the system at which a stress-release event becomes certain) and τ (the time taken for the system to accumulate the critical stress X_cr, in the absence of any stress-release events). We attach a physical interpretation to X(t), in the solar flare context, in Section 2.3 and Appendix A. The amount of stress released at the ith event, ΔX_i , is a random variable, drawn from a user-specified PDF, $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$, where $X({t}_{i}^{-})$ is the stress in the system immediately prior to the ith event. In the standard configuration, $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$ is fixed as a power law, but other options exist (Carlin & Melatos 2019a, 2021). Making η conditional on $X({t}_{i}^{-})$ is necessary in order to ensure that the stress remains positive-definite and is plausible physically. The second random variable in Equation (1) is N(t), a stochastic function that counts the number of events up to time t. It is determined iteratively via the waiting time between each event.

We assume the instantaneous event rate is a monotonically increasing function of the stress in the system, viz.

$$\begin{eqnarray}&&\lambda [X(t)]=\displaystyle \frac{\alpha }{1-X(t)},\end{eqnarray} \tag{ 2 }$$

where α = λ₀ τ is a dimensionless control parameter and λ₀ = λ(X = 1/2)/2 is a reference rate. The long-term statistical output of the SDP process does not depend strongly on the functional form of λ[X(t)], so long as it diverges in the limit X → 1 as in Equation (2) (Fulgenzi et al. 2017; Carlin & Melatos 2019a). Equation (2) implies that the probability of an event occurring approaches one as the stress approaches the critical stress. As the stress increases deterministically between events, the PDF of waiting times Δt following the ith event is that of a variable-rate Poisson process (Cox 1955; Fulgenzi et al. 2017),

$$\begin{eqnarray}\begin{array}{rcl}p[{\rm{\Delta }}t\,| X({t}_{i}^{+})] & = & \lambda [X({t}_{i}^{+})+{\rm{\Delta }}t]\\ & & \times \,\exp \left\{-{\displaystyle \int }_{{t}_{i}^{+}}^{{t}_{i}^{+}+{\rm{\Delta }}t}{\rm{d}}t^{\prime} \lambda \left[X(t^{\prime} )\right]\right\},\end{array}\end{eqnarray} \tag{ 3 }$$

where $X({t}_{i}^{+})$ is the stress immediately following the event at time t_i .

Analytically solving the coupled Equations (1)–(3) to calculate the PDFs of event waiting times or sizes is usually not feasible, except for particular choices of $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right];$ see Section 6 of Fulgenzi et al. (2017) for an example. Instead, it is simple to run the following automaton to generate numerical solutions:

• 1.  
  Pick Δt from Equation (3), given the current stress X.
• 2.  
  Update the stress to X + Δt to account for the deterministic evolution.
• 3.  
  Pick ΔX from $\eta \left[{\rm{\Delta }}X\,| \,X+{\rm{\Delta }}t\right]$, and subtract it from the stress.
• 4.  
  Repeat from step 1.

Given α and $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$ the automaton generates a time-ordered sequence of waiting times and sizes. From the sequence, we can calculate the long-term PDFs for waiting times and sizes (Fulgenzi et al. 2017; Carlin & Melatos 2019a), as well as the cross-correlation (Melatos et al. 2018; Carlin & Melatos 2019a), autocorrelations (Carlin & Melatos 2019b), and other observables.

We identify the stress that accumulates between flares and relaxes at a flare with the spatially averaged magnetic energy density in a given active region. One could equally choose a different physical quantity, such as the magnetic shear, depending on the particular microphysics of the flare trigger. We assume that active regions have independent stress reservoirs, i.e., the coronal magnetic fields of different active regions do not interact strongly, and all flares from the same active region extract energy from one reservoir. An idealized toy model relating the foregoing definition of stress to magnetic energy input from the photosphere is outlined in Appendix A, as a simplified but concrete illustration of the physical picture under consideration.

Two major simplifying assumptions are that the rate at which energy is fed into the reservoir, τ⁻¹, is constant in time, and so is the critical, spatially averaged magnetic energy density, X_cr, for a given active region. These assumptions are motivated in Appendix A, but a key goal of the paper is to test their veracity. We also assume that flares reduce X(t) instantaneously. This assumption is defensible, as the typical time between flares (typically hours) is much larger than the typical duration of flaring events (typically minutes) (Fletcher et al. 2011). In what follows, we do not prescribe a particular functional form for $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$, to keep the SDP process as flexible as possible.

The SDP process generates sequences of dimensionless waiting times $\left\{{\rm{\Delta }}{t}_{i}^{\mathrm{SDP}}\right\}$ and sizes $\left\{{\rm{\Delta }}{X}_{i}^{\mathrm{SDP}}\right\}$. It is natural to ask how they correspond to directly observable sequences of waiting times $\left\{{\rm{\Delta }}{t}_{i}^{\mathrm{obs}}\right\}$ and sizes $\left\{{\rm{\Delta }}{s}_{i}\right\}$. The waiting times satisfy

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{i}^{\mathrm{obs}}={\rm{\Delta }}{t}_{i}^{\mathrm{SDP}}\tau ,\end{eqnarray} \tag{ 4 }$$

where τ is unknown a priori but can be estimated as discussed in Appendix B. The sizes are related in a more complicated way, because some energy is released through channels other than the soft X-ray flux. We make the simplifying assumption that the peak soft X-ray flux multiplied by the duration of the flare (henceforth Δs_i for the ith flare in a given active region) is proportional to the stress released from the reservoir. This implies

$$\begin{eqnarray}&&{\rm{\Delta }}{s}_{i}\propto {\rm{\Delta }}{X}_{i}^{\mathrm{SDP}}{X}_{\mathrm{cr}}.\end{eqnarray} \tag{ 5 }$$

The unknown constant of proportionality in Equation (5) hampers direct parameter estimation of X_cr from an observed size PDF.

Full parameter estimation for the SDP process given a set of observed waiting times and sizes is discussed in Melatos & Drummond (2019). It lies outside the scope of this paper. For completeness, in Appendix B we also describe an alternative parameter estimation procedure using a hierarchical Bayesian scheme. The advantage of the hierarchical scheme is that it clarifies, as a matter of principle, the information content and flow in the estimation problem, i.e., the parameter combinations that can be estimated uniquely, and the data components that inform each parameter estimate. The disadvantage is that it can be implemented in practice only if $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$ is drawn from a certain class of mathematical functions. The favored class produces waiting-time and size PDFs that do not match solar flare observations (see Section 5), so the hierarchical scheme is not applied to real data in this paper. Nonetheless, it is included for the benefit of the reader, who may wish to develop it further for solar flare analysis in the future.

When the system is slowly driven, we have X(t) ≪ 1, as events are usually triggered before the stress accumulates to near the threshold. When $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$ is a power law, the predicted waiting-time PDF p(Δt) is an exponential, and the predicted size PDF p(ΔX) is a power law over multiple decades. We show the qualitative behavior of the stress in the slowly driven regime in the top left panel of Figure 1, where for the sequence of 20 stress-release events shown, we have 0.1 < X(t) < 0.6. In the bottom left panels of the same figure, we show the predicted waiting-time and size PDFs. The size PDF p(ΔX) has a turnover in logarithmic slope at ΔX ≈ 10⁻² due to the particular choice of $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$ (see caption for the functional form).

Zoom In Zoom Out Reset image size

Depending on the choice of $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$, there may be detectable backward cross-correlations (i.e., a correlation between the size of an event and the waiting time since the preceding event). This is because the size of an event cannot exceed the amount of stress in the system, so longer periods of stress accumulation allow for the possibility of larger events. This backward cross-correlation is especially pronounced with the choice $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]\propto \delta \left[{\rm{\Delta }}X-X({t}_{i}^{-})\right]$, where δ(...) is the Dirac-delta function. This choice forces the stress reservoir to empty at each event, and it collapses the SDP process to other stochastic processes in the literature, e.g., forest fire models (Daly & Porporato 2006). In the solar flare context, a backward cross-correlation is a long-standing prediction of the Rosner & Vaiana (1978) stress buildup model, as well as the "reset" model of Hudson et al. (1998), Hudson (2019). If instead $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$ prefers small stress-release events, e.g., if it is a power law (Fulgenzi et al. 2017), the backward cross-correlation is small.

When the system is rapidly driven, we have $X({t}_{i}^{-})\lesssim 1$, i.e., the stress is driven close to the threshold before each event. As the dynamics of the system are strongly influenced by the presence of the threshold, we sometimes call this regime "threshold-driven." A prediction of the SDP process in this regime is that the size and waiting-time PDFs should "match," i.e., they should be the same distribution, up to a linear scaling (Fulgenzi et al. 2017). For example, if $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$ is a power law, both the waiting-time and size PDFs are power laws. We show this in the bottom right panels of Figure 1. We also show the qualitative behavior of the stress versus time in this regime in the top right panel of the same figure, where we see X(t) → 1 (i.e., unit critical stress) before every stress-release event.

A strong forward cross-correlation (i.e., a correlation between the size of an event and the waiting time to the next event) is observed in this regime. This is because a large stress-release event results in a longer delay before the system accumulates enough stress to approach the threshold. In the solar flare context, this cross-correlation is predicted by the "saturation" model of Hudson et al. (1998) and Hudson (2019). The smooth transition between the fast-driven regime (α ≪ 1), with high forward cross-correlations, and the slow-driven regime (α ≫ 1), is shown in Figure 2.

Zoom In Zoom Out Reset image size

Ideally, we fit the parameters of the SDP process (τ, X_cr, λ₀, and $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$) to a paired sequence of solar flare sizes and waiting times, and then infer what regime applies. In practice, however, at least some of these parameters vary between active regions, which limits us to sequences of length N ≲ 50 (typical active regions have at most dozens of events; see Section 4), which are insufficient to infer four or more parameters. Instead, we combine the results in Sections 3.1 and 3.2 to deduce qualitative features, which must be present in multiple observables simultaneously (e.g., PDFs and cross-correlations), if the SDP process operates in a particular α regime—or indeed operates at all. Specifically, we summarize the behavior in Sections 3.1 and 3.2 into the following observationally testable prediction. If solar flares are well-modeled with a rapidly driven process (e.g., the SDP process with α ≪ 1), then we should see large forward cross-correlations, accompanied by waiting-time and size PDFs with the same shape in individual active regions. This coincident signature should become more prominent as α (or a proxy thereof) decreases. We quantify this prediction in Section 5.2.

With the test above, we also ameliorate the impact of a potentially mis-specified model. Although it is agnostic about the microphysics, the SDP process is not the only possible prescription of stress accumulation and release. If some of the assumptions outlined in Section 2.3 do not hold, a different underlying model may underpin solar flares. For example, one may build a model in which stress accumulates according to a Brownian random walk with some underlying drift, until a threshold is breached, at which point a stress-release event is triggered (Carlin & Melatos 2020). In the Brownian stress-accumulation model, X(t) is always driven to the threshold before each event. If drift occurs faster than diffusion, we would predict the same coincident signature as for the SDP process: large forward cross-correlations are accompanied by matching waiting-time and size PDFs.

We use the publicly available flare summary data hosted by the National Geophysical Data Center (NGDC) ⁴ for flares before 2015 June 29, and by the Space Weather Prediction Center (SWPC) ⁵ for flares after 2015 June 28. We combine these two data sources into a homogeneous, cleaned database (henceforth "catalog"), with some anomalies corrected as described in Appendix C. These data are collated with the flare start epochs, t^s, peak epochs, t^p, and end epochs, t^e. The peak epoch is defined as the epoch that contains the highest peak flux, after the start epoch. The end epoch is defined as the epoch at which the flux returns to half of the difference between the peak flux and the background flux. The peak flux (irradiance) of each flare, f^p, is calculated with reference to the recorded flare class, in units of W m⁻². The longitude and latitude of the flare are also often available, if the flare is associated with an active region. For clarity, we henceforth denote as x_i,k an arbitrary measurement of the variable x in the ith flare from the kth active region.

We define the waiting time between two flare epochs as ${\rm{\Delta }}{t}_{i,k}={t}_{i+1,k}^{{\rm{p}}}-{t}_{i,k}^{{\rm{p}}}$. One could equally use the flare start or end epochs to define the waiting time. We use the flare peak epoch because it is less influenced by the background level, as discussed further in Section 4.2. We define the flare size as ${\rm{\Delta }}{s}_{i,k}={f}_{i,k}^{{\rm{p}}}\left({t}_{i,k}^{{\rm{e}}}-{t}_{i,k}^{{\rm{s}}}\right);$ that is, we multiply the irradiance by the duration of the flare to obtain a flare size with dimensions of energy per unit area. The duration of an active region is ${\rm{\Delta }}{T}_{k}={t}_{{N}_{k},k}^{{\rm{e}}}-{t}_{1,k}^{{\rm{s}}}$, where there are N_k flares in the region. The average flare rate in an active region is λ_k = N_k /ΔT_k . The average λ_k systematically overestimates the true flare rate, as the duration ΔT_k is between two flare epochs, both of which are counted in N_k . An unbiased estimate would take ΔT_k as the difference between the epochs when the active region appears and disappears, but neither epoch is recorded in the GOES flare summary data. ⁶ A concern we touch on further in Section 4.3 is that λ_k = N_k /ΔT_k assumes that the flare catalog is complete, i.e., all flares from a given region are detected and correctly attributed to that region.

As of 2022 September 1, there are 83,283 flares in the catalog, 49,328 of which are associated with an active region. There are 2429 active regions with N_k ≥ 5, accounting for 42982 flares. The median number of flares in an active region (when considering only regions with N_k ≥ 5) is 12, but the average is 17.7, as the distribution is peaked at N_k =5 and monotonically decreases with N_k . We show the probability mass function of the number of flares in an active region, Pr(N_k ), in the top panel of Figure 3 for regions with N_k ≥ 5. The PDF of observed flare rates, p(λ_k ), across all active regions with N_k ≥ 5, is displayed in the bottom panel of the same figure. The distribution of λ_k is well-described by a log-normal, with mean 0.08 hr⁻¹, and standard deviation 0.7 hr⁻¹. Henceforth, we typically only include regions with N_k ≥ 5 in our analysis, unless stated otherwise, as many of the subsequent statistical tests, such as the cross-correlation(s), have low statistical power with smaller sample sizes.

Zoom In Zoom Out Reset image size

As described in Section 2.2 of Wheatland (2001), the GOES flare detection algorithm involves a selection effect that obscures the detection of flares following a large flare, due to the enhanced background soft X-ray flux. To detect a flare, the flux must monotonically increase for four consecutive minutes, with the last value being 1.4 times the value from three minutes earlier. Hence, a flare must produce a 40% increase above the background flux. However, while flares typically rise rapidly to peak flux, the soft X-ray emission is observed to decay on a timescale of hours (Benz 2016). Thus, even large flares may be obscured by the enhanced background flux following, say, an X1 (f^p = 10⁻⁴ Wm⁻²) class flare.

The above detection algorithm introduces a secondary bias in the catalog. A flare, say with recorded peak flux ${f}_{i}^{{\rm{p}}}$, that occurs during the decay of another flare, say with recorded peak flux ${f}_{i-1}^{{\rm{p}}}$, has an enhanced peak flux compared to the case of the same flare occurring during a period of low background flux, viz.

$$\begin{eqnarray}&&{f}_{i}^{{\rm{p}}}\approx {f}_{i,{\rm{true}}}^{{\rm{p}}}+{f}_{i-1}^{{\rm{p}}}\exp \left(\displaystyle \frac{-{\rm{\Delta }}{t}_{i-1}}{\tau }\right),\end{eqnarray} \tag{ 6 }$$

where τ is the timescale on which the peak flux from the (i − 1)th flare decays, and ${f}_{i,{\rm{true}}}^{{\rm{p}}}$ is the "true" peak flux of the ith flare. A secondary effect is that the duration of the ith flare will have a reduced duration, due to the decaying contribution of the previous flare to the background flux. That is, during the ith flare, the background flux decays by a factor of $\exp \left[-\left({t}_{i}^{{\rm{e}}}-{t}_{i}^{{\rm{s}}}\right)/\tau \right]$, which reduces the time taken for the flux to decay back to half the difference between the peak and pre-flare fluxes. These two effects run counter to one another, as Δs_i is defined as the product of peak flux and duration. They are likely present in a small subset of the catalog, as only 1.2% of all flares from regions with N_k ≥ 5 have start times before the end time of the previous flare.

A comprehensive reanalysis of the soft X-ray flux time series, with appropriate background subtraction, may reveal flares that were not detected with the original flare detection algorithm and correct the biases in flare size outlined above. Such a reanalysis lies outside the scope of this paper. Acknowledging that the GOES catalog is incomplete, we mitigate the impact on our analysis by creating a secondary masked catalog that only includes flares with Δs ≥ 10⁻³ Jm⁻² (i.e., class C1 multiplied by the median duration of 10³ s and higher), where we know the fraction of missed flares is smaller. This masked catalog is used in Sections 4.3 and 5.1 when we perform comparative studies and parametric fits for the waiting-time and size PDFs.

When we aggregate flare waiting times and sizes across all active regions with N_k ≥ 5, we obtain the PDFs shown in the top and bottom panels of Figure 4, respectively. The data are shown as the black histograms. The gray region in the bottom panel shows the flare sizes that are not included in the masked catalog, as described in Section 4.2. The difference between the masked and full catalog is not visible in the top panel, as the effect on p(Δt) is minimal.

Zoom In Zoom Out Reset image size

Three empirical trial distributions with common analytic forms are overlaid on each PDF for the masked catalog. The overlaid distributions have parameters fixed to their maximum likelihood values, given the data. By eye, it is clear that the log-normal distribution best describes the aggregated waiting-time PDF ⁷ , p(Δt), while a power-law distribution best describes the aggregated size PDF, p(Δs). We formalize this comparison using the corrected Akaike Information Criterion (AICc; Akaike 1974; Hurvich & Tsai 1989). The AICc calculates the model, among a set of possible models, that minimizes the information loss, while accounting for potential bias due to the number of model parameters and the sample size. When calculated for p(Δt), the relative probability for a log-normal describing the data over an exponential or a power law is ${e}^{3\times {10}^{3}}$ or ${e}^{3\times {10}^{4}}$, respectively. When calculated for p(Δs), the relative probability for a power law describing the data over an exponential or a log-normal is ${e}^{3\times {10}^{4}}$ or ${e}^{4\times {10}^{3}}$, respectively. If we instead only mask flares with peak flux less than 10⁻⁶ Wm⁻² (i.e., class C1), we find that a log-normal best describes the data. This preference is also noted in Verbeeck et al. (2019). It arises because of the small number of flares in this alternatively masked catalog with 10⁻⁴ < Δs/Jm⁻² ≲ 10⁻³. If we fit the full, unmasked catalog (i.e., include the flares shaded in gray in the bottom panel of Figure 4), we find that p(Δs) is again fitted best with a log-normal. We remind the reader that p(Δs) is an observed distribution; it is a function of both the underlying generative physics (i.e., how much energy is released in each flare) and the systematic observational biases (i.e., how many and what flares are detected or not). The masking described in Section 4.2 is a first pass at accounting for some of these biases.

If we assume that the same α and $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$ apply to all regions, Figure 4 is inconsistent with the SDP framework for any choice of α and $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$. However, as we discuss in Section 2.3, there is good reason to believe that α (and perhaps even $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$) may vary from region to region.

For regions with N_k ≥ 5, we disaggregate the data and ask what PDF shape best characterizes each region's waiting-time and size PDFs, rather than considering only the aggregated data set as in Section 4.3. Using the AICc, we find that, for waiting-time PDFs, 63% of regions are fitted best with an exponential distribution, 20% with a log-normal, and the remainder with a power law. For size PDFs, 63% of regions are fitted best with a power-law distribution, 24% with an exponential, and the remainder with a log-normal. These proportions are broadly consistent with previously published results, i.e., while the aggregated p(Δt) is fitted best with a log-normal, individual regions are often fitted best with an exponential (with varying rates) (Wheatland 2000b). Individual regions have p(Δs) that is usually fitted best with a power law, but can occasionally be better represented with log-normal or exponential distributions, especially in regions with lower N_k . As an example of the different shapes that distributions of waiting times and sizes can have in different regions, we display four arbitrary but representative per-active-region fits in Appendix E.

One may reasonably ask whether there are clear trends, or predictors, for the waiting-time and size PDFs in any given region. An exhaustive search for predictors, e.g., with a multiple regression analysis, lies outside the scope of this paper, but one sensible first step is to see if the shape that fits best evolves with λ_k . This test is performed for the waiting-time and size PDFs in the top and bottom panels of Figure 5, respectively. This figure is constructed by binning all active regions with N_k ≥ 5 into 20 uniformly spaced bins between λ_k = 0.025 hr⁻¹ and λ_k = 0.3 hr⁻¹, then calculating which distribution fits best the waiting times and sizes using the AICc. For the waiting-time PDFs, the proportion of regions fitted best with a power law stays roughly constant as λ_k increases, at around 15%, while the proportion of regions fitted best with a log-normal grows with λ_k , at the expense of the exponential. For the size PDFs, we see a similar trend, with the proportion of regions fitted best with a power law staying roughly constant as λ_k increases, this time at around 65%, while the proportion of regions fitted best with a log-normal grows with λ_k , at the expense of the exponential. In both panels, we plot the cumulative distribution function (CDF) of λ_k to remind the reader that each λ_k bin does not contain the same number of regions; 90% of regions with N_k ≥ 5 have 0.028 < λ_k /(1 hr⁻¹) < 0.27.

Zoom In Zoom Out Reset image size

Can the evolution of the proportion of each shape versus λ_k be explained with the SDP model? Under the assumption that λ_k is a tracer of the driving rate, i.e., λ_k ∝ α⁻¹, we expect to see a greater proportion of exponentially distributed waiting times at low λ_k (high α). This is the regime in which the stress does not approach the threshold at X = X_cr before each event, so waiting times are uncorrelated with sizes and are (broadly) Poissonian, i.e., the waiting times are exponentially distributed. This expectation broadly conforms with what we see in the top panel of Figure 5.

The evolution versus λ_k in the bottom panel of Figure 5 is harder to explain with the SDP model, if each region has the same $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$. We expect p(Δs) ∝ η, when α is low, but we see all of the three possible shapes represented at the highest values of λ_k . This implies at least one of the following: (i) η truly varies from one region to the next, which implies different stress-release mechanisms are at play in different regions; or (ii) λ_k ≳ 0.2 hr⁻¹ does not correspond to α ≪ 1, and hence p(Δs) ∝ η; or (iii) the small sample size of events in each region results in the AICc not favoring the "true" size distribution; or (iv) the obscuration effects described in Section 4.2 is stronger with higher λ_k , due to the enhanced background flux in regions that have many flares in a short period of time. To test (iii), we generate Figure 5 again for regions with N_k ≥ 10, instead of N_k ≥ 5. For both the waiting times and sizes, the proportion fitted best by an exponential drops by ∼10% in each λ_k bin, while the proportion fitted best by a log-normal increases. Qualitatively, however, the evolution with λ_k remains consistent with what is seen in Figure 5.

The correlation between flare sizes and subsequent waiting times, i.e., the "forward" cross-correlation, is denoted as ρ_+,k , while the correlation between flare sizes and the preceding waiting times, i.e., the "backward" cross-correlation, is denoted as ρ_−,k . We calculate these correlations using the Spearman correlation coefficient (Lehmann & D'Abrera 2006).

The PDFs of forward and backward cross-correlations measured in all active regions with N_k ≥ 5 are displayed in the top panel of Figure 6. The PDFs are broad, mostly because the median N_k is 12 (i.e., small). However, p(ρ_+,k ) and p(ρ_−,k ) are different; the latter has a median of ρ_−,k = 0.0, while the former has a median of ρ_+,k = 0.08. A Kolmogorov–Smirnov (KS) two-sample test provides quantitative evidence of the difference, returning a p-value of 10⁻¹⁴.

Zoom In Zoom Out Reset image size

The SDP process predicts a high forward cross-correlation for α ≪ 1, as discussed in Section 3.2. The data do not show evidence in favor of such a correlation in the majority of active regions. In the bottom panel of Figure 6, we do not see a clear visual trend between ρ_+,k and λ_k (for regions with N_k ≥ 10), although a Spearman correlation test returns a small but nonzero correlation of 5 × 10⁻² (p-value of 0.06). We remind the reader that both ρ_+,k and λ_k are empirical estimates for each active region, and the number of events per region is small. We do not plot uncertainties in the bottom panel of Figure 6, but they are typically large, on the order of the magnitude of the central estimate. These results imply that either (i) all active regions have α ≳ 1, where no forward cross-correlation is expected, i.e., λ_k ≳ 0.2 hr⁻¹ does not correspond to the threshold-limited regime; or (ii) some active regions have α ≪ 1, but either the trigger threshold or the driving rate is not constant with time, i.e., the random process triggering flares does not conform to the assumptions made in the SDP framework.

As described in Section 4.2, the GOES catalog is not complete, i.e., while a given active region is visible, not all flares that occur are recorded. Wheatland (2000a) noted that the obscuration in Section 4.2 creates an artificial forward cross-correlation, which explains why the median ρ_+,k is positive.

Another way to probe these data is to ask whether regions with relatively high ρ_+,k exhibit "matching" functional forms for p(Δt) and p(Δs), as the SDP model predicts for α ≲ 1. We quantify the degree to which the PDFs match via the p-value ${{ \mathcal M }}_{k}$ of the KS two-sample test applied to the sampled PDFs, {Δt_i,k } and {Δs_i,k }. To calibrate, we perform a Monte Carlo simulation using sequences of events drawn from the SDP model. For each value of α, we simulate 10⁵ "regions". For each region, we generate M events, where M is a random number drawn from p(N_k ), restricted to values N_k ≥ 10, as empirically measured for the GOES regions (the distribution for N_k ≥ 5 shown in Figure 3). From these M events, we calculate summary statistics, e.g., ${ \mathcal M }$ and ρ₊. This results in 10⁵ pairs of $({ \mathcal M },{\rho }_{+})$, from which we construct a two-dimensional kernel density estimate (KDE; Wand & Jones 1995). The KDE estimates the true joint PDF $p({ \mathcal M },{\rho }_{+})$, and it is qualitatively equivalent to a smoothed, two-dimensional histogram.

The result of this procedure for four values of α is shown in the top panel of Figure 7. For α = 0.01 (dark orange lines), the 10% and 50% credible interval contours are invisible in the top right corner of the panel, as the vast majority of simulated regions have ρ₊ ≈ 1 and ${ \mathcal M }\approx 1$. For α = 0.1, the KDE spreads out slightly, with the 50% credible interval reaching ρ₊ ≈ 0.8 and ${ \mathcal M }\approx 0.9$. For α = 1, the PDF shifts such that the 50% credible interval contour stretches from ${ \mathcal M }=1$ to ${ \mathcal M }\approx 0.8$, while we have 0 ≲ ρ₊ ≲ 0.8. For α = 10, we find ρ₊ centered around zero, while the 50% credible interval of ${ \mathcal M }$ extends to ${ \mathcal M }\approx 0.3$. The bottom panel of Figure 7 repeats the exercise for the backward cross-correlation. It confirms that ρ₋ and ${ \mathcal M }$ are uncorrelated, as predicted elsewhere (Fulgenzi et al.2017; Carlin & Melatos 2019a, 2019b). Integrating over ρ_±, the reader would find that the marginal distribution of ${ \mathcal M }$ is broad in both panels, for α ≳ 1, as it is difficult to confidently reject the null hypothesis that two sets of samples are from the same distribution, when dealing with small sample sizes.

Zoom In Zoom Out Reset image size

Monte Carlo calibration in hand, we present the equivalent data for all regions in the GOES catalog with N_k ≥ 10 in Figure 8. To compute ${{ \mathcal M }}_{k}$ for each region, we first normalize both the waiting times and the sizes by their respective means for that region, before applying the KS two-sample test. In the top (bottom) panel, we see no clear relationship between ρ_+,k (ρ_−,k ) and ${{ \mathcal M }}_{k}$. Marginalizing over ${{ \mathcal M }}_{k}$, we recover the PDFs of ρ_+,k and ρ_−,k , shown in the top panel of Figure 6. If, for the sake of argument, we assume that all regions have the same value of α, we can compare the KDEs in Figure 8 to those in Figure 7. Under this assumption, we are pushed into the regime 1 ≲ α ≲ 10, as the $p({\rho }_{+,k},{{ \mathcal M }}_{k})$ KDE is slightly offset from the horizontal axis (i.e., median ρ_+,k > 0), but the 50% credible interval for ${{ \mathcal M }}_{k}$ only extends to ${{ \mathcal M }}_{k}\lesssim 0.6$. Even if α is not exactly the same in each region, we can interpret the above result as ruling out the notion that a large proportion of regions have α ≲ 1, as if that were the case we would see higher values of ${{ \mathcal M }}_{k}$ associated with higher ρ_+,k more often than in Figure 8.

Zoom In Zoom Out Reset image size

While there is a wealth of information available in flare waiting-time and size PDFs, as in Section 5.1, and cross-correlations, as in Sections 5.2 and 5.3, one may also consider statistics that quantify the longer-term memory of the stress in the system. For example, as in the context of neutron star glitches (Carlin & Melatos 2019b), studying the autocorrelation between consecutive waiting times, ρ_Δt , or between consecutive sizes, ρ_Δs , allows one to place constraints on applicable SDP model parameters. In the fast-driven regime, α ≪ 1, we predict ρ_Δt = 0 and ρ_Δs = 0, as we have X(t) → X_cr before every stress-release event, resulting in the system resetting at every event. If we have $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]\propto \delta [X({t}_{i}^{-})]$, i.e., all stress in the system is released at each event, the same prediction of no autocorrelations holds—again the system is reset at every event. Observations of nonzero autocorrelations therefore rule out certain regimes in the SDP framework.

When we calculate ρ_Δt,k and ρ_Δs,k for all regions in the GOES catalog with N_k ≥ 10, we find that the PDFs p(ρ_Δt,k ) and p(ρ_Δs,k ) are broad, akin to the p(ρ_+,k ) shown in Figure 6. The medians are 0.11 and 0.10, respectively. This result is incongruent with the assumption α ≪ 1 in all regions, as α ≪ 1 implies both PDFs should have a median of zero. These results likely stem from the obscuration effects described in Section 4.2. One tentative physical interpretation of a positive ρ_Δt and ρ_Δs is via analogy to terrestrial earthquakes, which exhibit aftershocks, i.e., large events are often followed by larger-than-average events, with smaller-than-average waiting times (Utsu & Ogata 1995).

Solving Equations (1)–(3) for the long-term observable PDFs, p(Δt) and p(ΔX), is analytically intractable for most choices of $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$. However, for the special case

$$\begin{eqnarray}&&\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]\propto {\left[X({t}_{i}^{-})-{\rm{\Delta }}{X}_{i}\right]}^{\delta },\end{eqnarray} \tag{ B1 }$$

we obtain the analytic result

$$\begin{eqnarray}&&p{(z)=(\alpha +\delta +1)(1-z)}^{\alpha +\delta },\end{eqnarray} \tag{ B2 }$$

where the variable z is either Δt or ΔX, i.e., the waiting-time and size PDFs are identical (Fulgenzi et al. 2017). In Equation (B1), we have δ > 0 and the constant of proportionality is set by the condition that η must integrate to unity between ΔX_i = 0 and ${\rm{\Delta }}{X}_{i}=X({t}_{i}^{-});$ see Section 6 and Appendix D of Fulgenzi et al. (2017) for a full derivation. Equation (B1) is a monotonically decreasing function of ΔX_i , i.e., small stress-release events are preferred over large events. Its specific functional form is reasonable but arbitrary; it is not inferred from solar flare data. We adopt it here as a pedagogical device to illustrate the advantages and disadvantages of a hierarchical Bayesian approach to analyzing flare data.

When we restore the dimensions to Equation (B2), and consider the set of observations in the kth active region, D_k = {Δt_i,k , Δs_i,k }, with 1 ≤ i ≤ N_k , we can write the likelihood

$$\begin{eqnarray}&&{ \mathcal L }({D}_{k}\,| \,{{\boldsymbol{\theta }}}_{k})=\displaystyle \prod _{i=1}^{{N}_{k}-1}p({\rm{\Delta }}{t}_{i,k}\,| \,{{\boldsymbol{\theta }}}_{k})\displaystyle \prod _{i=1}^{{N}_{k}}p({\rm{\Delta }}{s}_{i,k}\,| \,{{\boldsymbol{\theta }}}_{k})\end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray}\begin{array}{rcl} & = & {\tau }_{k}^{1-{N}_{k}}\,{\xi }_{k}^{-{N}_{k}}{\left({\beta }_{k}+1\right)}^{2{N}_{k}-1}\displaystyle \prod _{i=1}^{{N}_{k}-1}\\ & & \times \,{\left(1-\displaystyle \frac{{\rm{\Delta }}{t}_{i,k}}{{\tau }_{k}}\right)}^{{\beta }_{k}}\displaystyle \prod _{i=1}^{{N}_{k}}{\left(1-\displaystyle \frac{{\rm{\Delta }}{s}_{i,k}}{{\xi }_{k}}\right)}^{{\beta }_{k}},\end{array}\end{eqnarray} \tag{ B4 }$$

where ξ_k is the constant of proportionality in Equation (5), which translates the drop in magnetic energy density ΔX to the observed flare size, Δs, and one has β_k = λ_0,k τ_k + δ_k . While there are N_k flare sizes in the region, there are only N_k − 1 waiting times. The three model parameters θ _k = {τ_k , ξ_k , β_k } are assumed constant in time for the lifetime of the active region. Bayes' theorem calculates the posterior probability distribution (henceforth "posterior"), p( θ _k ∣ D_k ), i.e., the probability density of parameter vector θ _k given the data D_k , using the prior probability distribution (henceforth "prior"), π( θ _k ), viz.

$$\begin{eqnarray}&&p({{\boldsymbol{\theta }}}_{k}\,| \,{D}_{k})\propto { \mathcal L }({D}_{k}\,| \,{{\boldsymbol{\theta }}}_{k})\pi ({{\boldsymbol{\theta }}}_{k}).\end{eqnarray} \tag{ B5 }$$

The normalizing constant of proportionality is often called the evidence, and while essential for model comparison, it is not relevant for the parameter estimation exercise below (Gelman et al. 2013).

Suppose, for the sake of illustration, that X_cr, λ₀, and every component of θ _k are approximately the same in all active regions. In the language of hierarchical Bayesian inference, this corresponds to assuming that the value of θ _k for each region is a random number drawn from a population-level distribution known as a "hyperprior," with narrow extent. For concreteness, we assume the hyperprior for each model parameter is a Gaussian with mean and standard deviation μ_a and σ_a respectively, with a ∈ {τ, ξ, β}. While β_k does depend on τ_k , the inference remains accurate so long as the covariance between τ_k , λ_0,k , and δ_k is minimal. The marginal posterior distribution for the parameters describing the hyperpriors is calculated as

$$\begin{eqnarray}&&p({\boldsymbol{\Lambda }}\,| \,{ \mathcal D })\propto \pi ({\boldsymbol{\Lambda }})\displaystyle \prod _{k}^{M}\int {\rm{d}}{{\boldsymbol{\theta }}}_{k}\,{ \mathcal L }({D}_{k}\,| \,{{\boldsymbol{\theta }}}_{k})\pi ({{\boldsymbol{\theta }}}_{k}\,| \,{\boldsymbol{\Lambda }}),\end{eqnarray} \tag{ B6 }$$

with Λ = {μ_a , σ_a }. In Equation (B6), π(Λ) is the prior on Λ, and ${ \mathcal D }=\{{D}_{1},\,\ldots ,\,{D}_{k},\,\ldots \,{D}_{M}\}$ is the set of all data, D_k , from M different active regions.

In practice, we can use a Monte Carlo sampler to estimate Equation (B6), and thus infer the posterior predictive distributions, p( θ ), where we drop the subscript k to signify that these distributions bring together information from all active regions to inform the inference. We opt to sample using a Hamiltonian Monte Carlo No U-Turn Sampler (Betancourt 2018), as implemented in the Stan programming language (Stan Development Team 2022).

To estimate the efficacy of the scheme in Appendix B.2 to infer population-level parameters, we first apply it to a set of synthetic data generated directly from the model. That is, we generate M = 50 fake "active regions," each with N_k = 100 events, with flare sizes and waiting times generated from the SDP framework with $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$ as in Equation (B1). In what follows, we use the compact notation $x\sim { \mathcal N }(\mu ,\sigma )$ to denote that the random variable x is drawn from a Gaussian distribution with mean μ and standard deviation σ. To generate the data for each region, we select a value ${\tau }_{k}\sim { \mathcal N }({\mu }_{\tau }=3\,{\rm{day}},{\sigma }_{\tau }=0.5\,{\rm{day}})$, a value ${\xi }_{k}\sim { \mathcal N }({\mu }_{\xi }\,=2\,{\rm{arb.}}\,{\rm{units}},{\sigma }_{\xi }=0.3\,{\rm{arb.}}\,{\rm{units}})$, and a value ${\beta }_{k}\sim { \mathcal N }({\mu }_{\beta }\,=5,{\sigma }_{\beta }=0.3)$. Each Gaussian is truncated such that all variates are positive.

After running the sampler with the synthetic data, we perform a posterior predictive check, i.e., we compare the samples from our posterior distributions p(τ), p(ξ), and p(β) with the injected distributions and the priors on the hyperparameters. We show the results in Figure 9. We see that the posterior samples (blue) appropriately overlap with the injected ground truth (red) for τ, ξ, and β.

Zoom In Zoom Out Reset image size

Despite the encouraging results in Figure 9, we do not apply this hierarchical Bayesian scheme to the GOES catalog in this paper. This is because, in most regions, neither the waiting-time nor the size PDFs follow the functional form in Equation (B2). When we include Equation (B2) in the set of options available for the AICc to select between, as in Section 5.1, fewer than 0.1% of active regions are fitted best with Equation (B2). Therefore, Equation (B4) is not an appropriate likelihood for the data. Other choices of the jump distribution in Equation (B1) may alleviate this problem, but they lie outside the scope of this paper.

An alternative approach is to build the likelihood out of a summary statistic, such as the average waiting time in a given region, 〈Δt〉_k , rather than each observed waiting time and/or size. In the SDP framework, we find (Fulgenzi et al. 2017; Millhouse et al. 2022)

$$\begin{eqnarray}&&\langle {\rm{\Delta }}t{\rangle }^{\mathrm{SDP}}=\displaystyle \frac{1}{\alpha +{\alpha }_{0}},\end{eqnarray} \tag{ B7 }$$

where α₀ ∼ 1 is a dimensionless constant whose exact value depends on the particular form of $\eta \left[{\rm{\Delta }}{X}_{i}\,| \,X({t}_{i}^{-})\right]$. Restoring the dimensions to Equation (B7), we can write the likelihood as

$$\begin{eqnarray}&&{ \mathcal L }\left(\langle {\rm{\Delta }}t{\rangle }_{k}\,| \,{{\boldsymbol{\theta }}}_{k}\right)=\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{k}}}\exp \left[\displaystyle \frac{-{\left(\langle {\rm{\Delta }}t{\rangle }_{k}-\tfrac{{\tau }_{k}}{{\lambda }_{0,k}{\tau }_{k}+{\alpha }_{0,k}}\right)}^{2}}{2{\sigma }_{k}^{2}}\right],\end{eqnarray} \tag{ B8 }$$

with θ _k = {τ_k , λ_0,k , α_0,k , σ_k }, where we assume the residual of 〈Δt〉_k with the model subtracted is normally distributed, with zero mean and standard deviation σ_k . With the likelihood in Equation (B8), one could perform population-level parameter estimation as described in Section B.2. We leave this, as well as an exploration of how θ _k depends on parameters intrinsic to each active region, to future work.