In the following we attempt to model event statistics with a combination of (i) a Gaussian distribution (originating from incoherent random processes), and (ii) a power law distribution, e.g., created by spicular activity in the transition region, (Figure 1). The Gaussian noise is defined in the standard way, N F d F = N 0 ⁡ exp − F − F 0 2 2 σ F 2 d F , (1) where F is the flux averaged over the duration of an event (measured here at a wavelength of 1,400 Å), N(F) is the histogram of observed structures, F ₀ is the mean value, σ _F is one standard deviation, and N ₀ is the normalized number of events.

The second distribution we use in our analysis is a power law distribution function, which is defined in the simplest way by, N F d F = N 0 F F 0 − α F d F , (2) where α _F is the power law slope of the relevant part of the distribution function.

The flux F _IRIS of an IRIS pixel is defined by, F IRIS = 4 π f E λ k A Ω , e r g c m − 2 s − 1 (3) where f is the observed flux in [DN] (data number per second), E _λ is the energy of the photon, k is the factor that converts the DN to the number of photons, Ω is the pixel size in units of steradians, A [cm²] is the effective area of IRIS, and the unrelated background is subtracted (Vilangot Nhalil et al., 2020).

In this study we use a “pixelized” size distribution method that is more efficient and easier to calculate than standard size distributions. The standard method to sample size distributions N(F) of SOC avalanches is generally carried out by an algorithm that detects fluxes of an avalanche event above some given threshold F > F _thr , traces its spatial A(t) and temporal evolution F(t), and infers the size of an avalanche from the spatio-temporal evolution after saturation. Such avalanche detections have been accomplished for 12 IRIS datasets in the study of Vilangot Nhalil et al., 2020. Because the development of an automated feature recognition code is a complex and a time-consuming task, which needs extensive testing, we explore here a new method that is much simpler to apply and requires much less data to determine the underlying power law slope α _F .

We can parameterize a pixelized IRIS image with a Cartesian grid, i.e., F i , j = F ( x i , y j ) , i = 0 , … , n x , j = 0 , … , n y , where n _x and n _y are the dimensions of the image, and Δx = Δy is the pixel size. We can model a 2-D image with a superposition of n _k spatial structures with avalanche areas A _k and average fluxes F k , where the size distributions follow a power law distribution, i.e., N ( F ) ∝ F − α F (Eq. 2). The total flux F tot of such a 2-D distribution, which serves here as an analytical model of a 2-D (IRIS) image, can then be written as, F tot = Σ k = 1 n k F k − α F A k , (4) where the avalanche areas A _k are required to be non-overlapping, but area-filling. Areas without significant avalanche structures, ( F < F thr ) , can be included, in order to fulfill flux conservation, or can be neglected if the flux maximum is much larger than the threshold value, i.e., F max ≫ F thr .

In our new method we decompose the flux F k and area A _k of all avalanche components down to the pixel size level, Δx. The two requirements of non-overlapping and area-filling topology yield a unique mapping of the avalanche number k to the pixel ranges i = [i ₁(k), i ₂(k)] and j = [j ₁(k), j ₂(k)], i.e., k↦i ₁(k), …, i, … i ₂(k) and j ₁(k), …j, … , j ₂(k). For instance, in the case of a rectangular area A _k , the avalanche area A _k is then defined by, A k = i 2 k − i 1 k * j 2 k − j 1 k Δ x 2 (5) Adding the areas A and fluxes F of all k avalanche components, we obtain then the following total flux F tot,pix , F tot,pix = Σ i = 1 n x Σ j = 1 n y F i , j − α F Δ x 2 , (6) which can be set equal to the value of F tot of the standard method (Eq. 4) and proves this way that the power law slopes α F of the two methods are identical. Thus, our new method is parameterized just by a different decomposition of elementary components than in the standard size distribution sampling.

As a caveat, we have to be aware that the method determines size distribution from a single image. If the used 2-D image is not representative, additional 2-D images need to be included.

The new pixelation method is used in the calculations of the values α F 1 listed in Table 1 and Figure 4.

The 12 analyzed 1,400 Å SJI images F(x, y) of IRIS are shown in Figure 2, which are identical in time and FOV (field-of-view) with those of Vilangot Nhalil et al., 2020, and are also identical with those used in the study on fractal dimension measurements (Aschwanden and Vilangot Nhalil, 2022). Note that events #6 and #7; are almost identical, except for a time difference of 20 min, which can be used for stability tests.

The 12 IRIS maps shown in Figure 2 have the following color code: The Gaussian distribution with values F(x, y) < F _thr below a threshold of F _thr is rendered with orange-to-red colors, while the power law function with the fat-tail F(x, y) > F _thr is masked out with white color. In other words, all the orange-to-red regions in the IRIS maps visualize the locations of incoherent random noise while the white regions mark the location of SOC-driven coherent avalanches (probably produced by spicular dynamics in the transition region). An even crispier representation of the spicular component F(x, y) > F _thr , is displayed with a black-and-white rendering (Figure 3), where black depicts locations with power law distributions, and white demarcates locations with Gaussian distributions.

The information content of an IRIS image can be described with a 2-D array of flux values F(x, y) at a given time t, or alternatively with a 1-D histogram N(F). Since we want to fit a two-component distribution function (i.e., with a Gaussian and a power law), we need to introduce a separator between the two distributions, which we derive from the full width at half maximum (see F ₂ in Figure 1). We fit then both distribution functions (Eqs 1, 2) separately, the Gaussian function in the range of [F ₁, F ₂], and the power law function in the range of [F ₂, F ₃], as depicted in Figure 1. The minimum flux (F ₁) and maximum flux (F ₃) are determined from the minimum and maximum flux value in the image. We are fitting the distribution functions with a standard Gaussian fit method, and with a standard linear regression fit for the logarithmic flux function. Note that the power law function N(S) appears to be a straight line in a logarithmic display only (Figure 1 bottom panel), i.e., log(N)-log(S), but not in a linear representation (Figure 1 top panel), i.e., lin(N)-lin(S), as used here.

The results of the fitting of the observed histograms are shown for all 12 datasets in Figure 4, where the Gaussian fit is rendered with a blue color, and the power law fit with a red color. We see that our two-component model for the distribution function produces accurate fits to the analyzed IRIS data (histograms in Figure 4) for seven datasets (# 4–9, 11), while it fails in 5 cases (# 1–3, 10, 12). On the other hand, 4 cases contain sunspots (# 1–3, 10) and coincide with the cases with power law fit failures.

If we would assume that all fluxes are generated by incoherent random noise, we would not be able to fit the histogrammed data at all. Obviously, we would under-predict most of the fluxes substantially (blue dashed curves in Figure 4), which underscores that the “fat-tail” power law function, a hallmark of SOC processes, is highly relevant for fitting the observed IRIS 1400 Å data here.

In a next step we investigate the numerical values of the power law slopes α _F of the flux distribution parameters F, which are listed in the third column of Table 1. At a first glance, it appears that these values vary wildly in a range of α _F = 0.94 to 2.13. However, Vilangot Nhalil et al., 2020 classified the 12 analyzed datasets into 4 cases containing sunspots, and 8 cases containing plages in the transition region without sunspots. From this bimodal behavior it was concluded that the power law index of the energy distribution is larger in plages (α _E > 2), compared with sunspot-dominated active regions (α _E < 2), (Vilangot Nhalil et al., 2020). In our investigation here, the 4 cases with sunspots exhibit substantially flatter power law slopes α _F (except #3), which indicates that sunspot-dominant distributions are indeed significantly different from those without sunspots (Table 1). Actually, we find an even better predictor of this bimodal behavior, by using the maximum flux F _max (Column 6 in Table 1). We find that flux distributions N ( F ) ∝ F − α F with maximum fluxes less than F _max ≲ 50 [DN] exhibit a power law value of α F obs ≈ 1.70 ± 0.15 , F max < 50 D N , (7) which includes the five datasets #6–9, 11. In contrast, the seven other datasets #1–5, 10, 12 have consistently higher maximum values, F _max ≳ 50 DN. Instead of using the maximum values F _max , we can also use the average fluxes and find the same bimodal behavior.

Even more significant is that this power law value is consistent with the theoretical prediction of the power law slopes (Aschwanden, 2012; Aschwanden, 2022a; Aschwanden et al., 2016), α F,SOC = 9 5 = 1.80 . (8) Thus we conclude that flux distributions have a power law slope that agree with the theoretial prediction under special conditions, such as for small maximum fluxes. Moreover we find that magnetic flux distributions with sunspots and large magnetic flux imbalances produce flatter slopes and failed power law fits (Tables 1 and 2), see Section 2.4.

In order to test the universality of the results we repeat the same analysis for 12 coincident HMI magnetograms onboard the Solar Dynamics Observatory (SDO), which have simultaneous times and identical spatial field-of-views. The 12 analyzed HMI images are shown in Figure 5, where black features indicate negative magnetic polarity, and white features indicate positive magnetic polarity. We see sunspots in at least four magnetograms (#1–3, 10), with two sunspots having a negative magnetic polarity (#1, 2), and two cases with positive magnetic polarity (#3, 10). All 12 magnetograms show mixed polarities, but some are heavily unbalanced (#1–5, 10–12).

We quantify the magnetic flux balance with the ratio q _pos , q pos = ∑ pos F i j ∑ pos F i j + | ∑ neg F i j | . (9) If the magnetic flux (line-of-sight) component is well-balanced, we would expect a value of q _pos = 0.5, assuming ∑ _pos = |∑ _neg |. Only 4 cases have approximately balanced fluxes (#6–9), namely, q _pos = [0.44, 0.43, 0.38, 0.44], while the other 6 cases have large flux imbalances, from q _pos = 0.04 to 0.99 (Table 2; Figure 6). The associated power law slopes of the four well-balanced cases are α _F = [1.67, 1.64, 1.79, 1.78] = 1.72 ± 0.07, which closely coincide with the theoretical SOC-prediction of α _F ≈ 1.80 (Aschwanden, 2012; Aschwanden, 2022a; Aschwanden et al., 2016).

We analyze the HMI data in the same way as the IRIS data, by fitting Gaussian distributions (blue curves in Figure 6) and power law distribution functions (red curves in Figure 6), which clearly show a “fat-tail” feature that is far in excess of the Gaussian function (blue dashed curves in Figure 6). We compare the power law slopes α _F obtained with the two completely different datasets from IRIS and HMI in Figure 7, using the “pixelation” method. The two datasets are found to be highly correlated (with CCC = 0.90, if we ignore the outlier #3). Nevertheless, the power law slopes α _F shown in Figure 7 are concentrated in two regimes, one that is consistent with our theoretical SOC prediction of α _F,IRIS = α _F,HMI ≈ 1.80, while a second cluster is centered around α _F,IRIS ≈ 1.0–1.5 and α _F,HMI ≈ 1.0–1.5 (Figure 7). In essence, we find four datasets (# 6–9) that are consistent with the SOC prediction for events with well-balanced flux q _pos ≈ 0.5, while a second group cannot reproduce the SOC model, but can be characterized with large unbalanced magnetic fluxes (# 1–5, 10–12). The flux imbalance, however, is not always decisive. Tests with variations of the FOV reveal that the arbitrary choice of the FOV (in HMI data) can be more important in deciding whether the calculated power law slope is universally consistent with SOC models, e.g., see event #11.

The physical interpretation of the HMI data is, of course, different for the IRIS data. In the previous analysis of IRIS data we interpreted the coherent statistics (in terms of SOC-controlled power law functions) due to spicular activity in the transition region. In contrast, using the HMI data, which provides the magnetic field line-of-sight component B _z , we can interpret the statistics of incoherent random distributions in terms of “salt-and-pepper” small-scale magnetic fields in the photosphere, and the coherent avalanche statistics in terms of SOC-controlled magnetic reconnection processes in nanoflares and larger flares (Table 3). Note that the two parameters α _F,IRIS amd α _F,HMI are observed independently from different spacecraft, as well as in markedly different wavelength bands, i.e., λ ≈ 1,400 Å for IRIS, and λ = 6,173 Å for HMI/SDO magnetograms, which measures the mean flux F from the line-of-sight magnetic field component B _z (x, y). Despite of the very different instruments and wavelengths, the power law slope α _F of the mean flux appears to be universally valid and consistent with the theoretical SOC prediction for datasets with approximate magnetic flux balance (Figure 7). However we learned that the magnetic flux balance and the absence of sunspots represent additional requirements to warrant the universality of the SOC slopes. This yields a testable prediction: If the field-of-view of each HMI magnetogram is readjusted so that the enclosed magnetic flux becomes more balanced and no sunspot appears in the FOV, the power law slope is expected to approach the theoretical universal value of α _F,IRIS ≈ α _F,HMI ≈ 1.80.

The most extensive statistical study on the size distribution of magnetic field features on the solar surface has been undertaken by Parnell et al. (2009). Combining magnetic field data from three instruments (SOT/Hinode, MDI/NFI, and MDI/FD on SOHO, a combined occurrence frequency size distribution was synthesized that extends over 5 decades, in the range of Φ = 2 × 10¹⁷–10²³ Mx (Parnell et al., 2009), N Φ ∝ Φ 0 − 1.85 ± 0.14 M x − 1 c m − 2 , (10) where the magnetic flux Φ is obtained from integration of the magnetic field B(x, y) over a thresholded area A = ∫ dx dy, Φ = ∫ B x , y d x d y M x . (11) If we equate the magnetic flux Φ with the mean flux F of an event in standard SOC models, we predict a power law slope of (Aschwanden, 2012; Aschwanden, 2022a; Aschwanden et al., 2016), using d = 3, D _V = 5/2, and γ = 1, α F,SOC = 1 + d − 1 γ D V = 9 5 = 1.80 , (12) which agree well with the result (Eq. 10) observed by Parnell et al. (2009). A lower value was found from cellular automaton simulations, N(Φ) ≈ Φ^−1.5±0.05 (Fragos et al., 2004), where flux emergence is driven by a percolation rule, similar to the percolation model of Seiden and Wentzel (1996), or Balke et al. (1993). Mathematical models have been developed to model the percolation phenomenon, based on combinatorial and statistical concepts of connectedness that exhibit universality in form of powerlaw distributions.

Power law-like size distributions are the hallmark of self-organized criticality systems. Statistical studies in the past have collected SOC parameters such as length scales L, time scales T, peak flux rates P, mean fluxes F, fluences and energies E = F × T, mono-fractal and multi-fractal dimensions (Mandelbrot, 1977), in order to test whether the theoretically expected power law size distributions, or the power law slopes of waiting times, agree with the observed distributions (mostly observed in astrophysical systems). The universality of SOC models (Aschwanden, 2012; Aschwanden, 2022a; Aschwanden et al., 2016) is based on four scaling laws: the scale-free probability conjecture N(L) ∝ L ^−d , classical diffusion L ∝ T ^β/2, the flux-volume relationship F ∝ V ^γ , and the Euclidean scaling law, P ∝ L ^γd , where d = 3 is the Euclidean dimension, β ≈ 1 is the classical diffusion coefficient, γ ≈ 1 the flux-volume proportionality, while D _A = 3/2 and D _V = 5/2 are the mean fractal dimensions in 2-D and 3-D Euclidean space. The standard SOC model is expressed in terms of these universal constants: d = 3, γ = 1, β = 1. Consequently, the four basic scaling laws reduce to N(L) ∝ L ⁻³, L ∝ T ^1/2, F ∝ L ^2.5, and P ∝ L ³. Since we measure the mean flux F in this study, our main test of the universality of SOC models if formulated in terms of the flux-volume relationship F ∝ V ^γ , leading to the power law slope α _F,SOC = 1.80 (Eq. 12).

The SOC-inferred scaling laws hold for a large number of phenomena. This implies that our SOC formalism is universal in the sense that the statistical size distributions are identical for each phenomenon, displaying a universal power law slope of α _F,SOC = 1.80. When we conclude that the power law slope α _F is universal, the SOC model implies that the flux-volume proportionality (γ ≈ 1) as well as the mean fractal dimension (d = 3, D _V ≈ 2.5) are universal too.

Once we establish the self-consistency of power law slopes between theoretical (SOC) and observed size distributions, the next question is what physical processes are at work. We envision four different types of phenomena (Table 3): (i) Gaussian random noise in IRIS data); (ii) spicular plage events in the transition region (described by the power law size distribution in IRIS data); (iii) salt-and-pepper small-scale magnetic structures (described by the random noise distributions in HMI magnetograms); and (iv) magnetic reconnection processes in flares and nanoflares (described by the power law size distribution in HMI data). However, there are deviations from these rules. We found that the power law distributions are modified in the presence of sunspots, when the magnetic flux is unbalanced, or when the FOV is arbitrarily chosen. Under ideal conditions, the SOC scaling laws are fulfilled universally, independent of the wavelength or plasma temperature. Magnetic field data (from HMI/SDO) or λ ≈ 1,400 Å (from IRIS) appear to produce emission in volumes that are proportional in the photosphere or transition zone, even when they are formed at quite different temperatures, i.e., T _phot ≈ 5,800 K in the photosphere and T _TR ≈ 10⁴–10⁵ K in the transition region.

Another ingredient of the SOC model is the scale-free probability conjecture, i.e., N(L) ∝ L ^−d = L ⁻³, which cannot be uniquely linked to a particular physical process. Parnell et al. (2009) conclude that a combination of emergence, coalescence, cancellation, and fragmentation may possibly produce power law size distributions of spatial scales L. Alternative models include the turbulence and the Weibull distributions (Parnell, 2002). Munoz-Jaramillo et al. (2015) study the best-fitting distribution functions for 11 different databases of sunspot areas, sunspot group areas, sunspot umbral areas, and magnetic fluxes. They find that a linear combination of Weibull and log-normal distributions fit the data best (Munoz-Jaramillo et al., 2015). Weibull and log-normal distributions combine two distribution functions, similar to our synthesis of a Gaussian-plus-power-law distribution.

A general physical scenario of a power law size distribution is the evolution of avalanches by exponential growth (Rosner and Vaiana, 1978), with subsequent saturation (logistic growth) after a random time interval, which produces an exact power law function (Aschwanden et al., 1998). Our approach to model the size distribution of solar phenomena with two different functions, employing a Gaussian noise and a power law tail, reflects the duality of incoherent and coherent random components, in both the data from IRIS and HMI (Table 3). In summary, incoherent random components include salt-and-pepper small-scale magnetic features, while coherent components include spicular avalanches, and magnetic reconnection avalanches from nanoflares to large flares.

The physical understanding of solar (or stellar) granulation has been advanced by numerical magneto-convection models and N-body dynamic simulations, which predict the evolution of small-scale (granules) into large-scale features (meso- or super-granulation), organized by surface flows that sweep up small-scale structures and form clusters of recurrent and stable granular features (Berrilli et al., 1998; Berrilli et al., 2005; Hathaway et al., 2000; Martinez-Sykora et al., 2008; Rieutord et al., 2008; Rieutord et al., 2010; Cheung and Isobe, 2014). An analytical model of convection-driven generation of ubiquitous coronal waves is considered in Aschwanden et al. (2018b). The fractal multi-scale dynamics has been found to be operational in the Quiet-Sun photosphere, in quiescent non-flaring states, as well as during flares (Uritsky et al., 2007; Uritsky et al., 2013; Uritsky and Davila, 2012). The fractal structure of the solar granulation is obviously a self-organizing pattern that is created by a combination of subphotospheric magneto-convection and surface flows, which are turbulence-type phenomena.

The interpretation of granulation as the cause of the Gaussian “noise” in IRIS data is controversial for two reasons: (i) The intensity measured by IRIS 1400 in non-magnetic areas has densities that originate from the middle chromosphere, rather than from the underlying photosphere. (ii) No convective signal propagates to these heights and densities, and thus the scale of granulation cannot be probed at these heights (Martinez-Sykora et al., 2015).

One prominent feature in the transition region is the phenomenon of “moss”, which appears as a bright dynamic pattern with dark inclusions, on spatial scales of L ≈ 1–3 Mm, which has been interpreted as the upper transition region above active region plages, and below relatively hot loops (De Pontieu et al., 1999; De Pontieu et al., 2014). Our measurement of structures in the IRIS 1,400 Å channel is sensitive to a temperature range of T _e ≈ 10⁴–10⁵ K, and thus is likely to include chromospheric and transition region phenomena such as: spicules II (De Pontieu et al., 2007), macro-spicules, dark mottles, dynamic fibrils, surges, miniature filament eruptions, etc. Theoretical models include the rebound shock model (Sterling and Hollweg, 1988), pressure-pulses in the high atmosphere (Singh et al., 2019), Alfvénic resonances (Sterling, 1998), magnetic reconnection models for type II spicules (De Pontieu et al., 2007), ion-neutral collisional damping (De Pontieu, 1999), leakage of global p-mode oscillations (De Pontieu et al., 2004), MHD kink waves (Zaqarashvili and Erdelyi, 2009), vortical flow models (Kitiashvili et al., 2013), and magneto-convective driving by shock waves (De Pontieu et al., 2007).

The fact that we obtain a power law size distribution (α _F = 1.70 ± 0.15, Table 1), which is very similar to solar flares in general, α _F,SOC = 1.80, implies the universality of the SOC framework. Furthermore we find power law-like size distributions for spicular events, rather than a Gaussian distribution, which tells us that spicule events need to be modeled in terms of SOC-driven avalanches, instead of Gaussian random distributions. As mentioned above, the difference between incoherent and coherent random processes is the following: Gaussian statistics reflect the operation of a memoryless stationary random process; while avalanche processes such as occurring in SOC systems are characterized by extended spatial and temporal correlations, i.e., the unfolding of an avalanche is influenced by the imprint of earlier avalanches on the system.

We propose that spicules around magnetic elements are responsible for the power law slope α _F observed in those areas. This appears to be a plausible interpretation since these dynamical phenomena are very relevant in plage and network areas. For instance, event #11 shows a fully unbalanced magnetic configuration, which supports the idea that strong magnetic fields, fragmented in small-scale elements in plage and/or network seems to be the relevant characteristics, rather than flux balance over an arbitrary FOV.

We interepret the random noise Gaussian distribution of magnetic fluxes in Quiet-Sun regions as small-scale magnetic field “pepper-and-salt” structures, also called magnetic carpet (Priest et al., 2002), where the black and white color in magnetograms (Figure 5) corresponds to negative and positive polarity. The fact that we obtain two distinctly different size distributions (Gaussian vs power law) indicates at least two different physical mechanisms, one being an incoherent random (Gaussian) process, the other one being a coherent (power law-like) avalanche process. The salt-and-pepper structure is generated apparently by an incoherent random process, rather than by a coherent avalanching process, according to our fits. This may constrain the origin of the solar magnetic field, being created by emergence, submergence, coalescence, cancellation, fragmentation, and/or small-scale dynamos, etc. Not all would be expected to yield Gaussian statistics (e.g., fragmentation processses often yield log-normal distributions; Verbeeck et al., 2019).

The re-arrangement of the stress-induced solar magnetic field requires ubiquitous and permanent (but intermittent) magnetic reconnection processes on all spatial and temporal scales. Our study finds power law size distributions, with a slope of α _F = 1.72 ± 0.07 from HMI magnetograms, which is similar to flares in general (see Aschwanden et al., 2016 for a review of all wavelengths (e.g., gamma-rays, hard X-rays, soft X-rays, UV, EUV, FUV, etc.). This tells us that there is a strong correlation between the photospheric field (in HMI images) and the transition region (in IRIS images), as evident from the cross-correlation coefficient of CCC = 0.90 shown in Figure 7. The fractal multi-scale dynamics apparently operates in the quiet photosphere, in the quiescent non-flaring state, as well as during flares in active regions (Uritsky and Davila, 2012).