# Avalanches in magnetohydrodynamical simulations

Henri Lamarre, Paul Charbonneau, Quentin Noraz, Antoine Strugarek, Alexis Blaise, Allan Sacha Brun, Mats Carlsson, Boris Vilhelm Gudiksen

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Henri Lamarre, Paul Charbonneau, Quentin Noraz, Antoine Strugarek, Alexis Blaise, et al. Avalanches in magnetohydrodynamical simulations. Physical Review E, 2025, 112 (6), pp.L063202. 10.1103/28ws-d57z. cea-05496232

HAL Id: cea-05496232

https://cea.hal.science/cea-05496232v1

Submitted on 5 Feb 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Avalanches in Magnetohydrodynamical simulations

Henri Lamarre,<fnref n="1" /> Paul Charbonneau,<fnref n="1" /> Quentin Noraz,<fnref n="2" /> Antoine Strugarek,³ Alexis

Blaise,³ Allan Sacha Brun,³ Mats Carlsson,<fnref n="2" />,⁴ and Boris Vilhelm Gudiksen,<fnref n="2" />,⁴

¹Département de Physique, Université de Montréal, Montréal, Québec H2V 0B3, Canada

²Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029 Blindern, 0315 Oslo, Norway

³Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France

⁴Rosseland Centre for Solar Physics, University of Oslo, P.O. Box 1029 Blindern, 0315 Oslo, Norway

(Dated: September 30, 2025)

Scale invariance is a hallmark of many natural systems, including solar flares, where energy release spans a vast range of scales. Recent computational advances, at the level of both algorithmics and hardware, have enabled high-resolution magnetohydrodynamical (MHD) simulations to span multiple scales, offering new insights into magnetic energy dissipation processes. Here, we study scale invariance of magnetic energy dissipation in two distinct MHD simulations. Current sheets are identified and analyzed over time. Results demonstrate that dissipative events exhibit scale invariance, with power-law distributions characterizing their energy dissipation and lifetimes. Remarkably, these distributions are consistent across the two simulations, despite differing numerical and physical setups, suggesting universality in the process of magnetic energy dissipation. Comparisons between the evolution of dissipation regions reveals distinct growth behaviors in high plasma-$\beta$ regions (convective zone) and low plasma-$\beta$ regions (atmosphere). The latter display spatiotemporal dynamics similar to those of avalanche models, suggesting self-organized criticality and a common universality class.

Scale invariance is found in many natural systems, from the evolution of forest fires to the distribution of earthquake magnitudes [1]. In the Sun, this property emerges for solar flares, where energy is released across a vast range of scales [2–4]. The most energetic of these events can have significant impacts on Earth [5] and as such, understanding them is key to produce predictions. However, solar flares have been hard to model from first principles due to the wide disparity of spatial and temporal scales at play. Furthermore, the small-scale nature of magnetic reconnection dynamics, and the complex behaviour of the resulting magnetic energy dissipation [6, 7] makes the modelling even more arduous.

Avalanche models have been successful in reproducing the statistical properties of flares [8–11], and other scale invariant processes, from earthquakes to forest fires [12, 13]. In the case of solar flares, these models might be capturing the intricate collective dynamics of reconnection events [14–16]. In fact, a nanoflare [17] can be conceptually linked to the idea of a local instability.

Bridging the gap between scale invariant models for solar flares and comprehensive magnetohydrodynamics (MHD) simulations is nevertheless lacking. Despite capturing the global statistics [8], avalanche models rely on ad hoc evolutionary rules, in contrast to the detailed physics described by MHD; reconciling them into a common physical framework is far from trivial [15, 18–21]. However, MHD simulations can provide physical grounding to these ad hoc rules. State-of-the-art MHD simulations can now resolve finer spatial scales, providing insights into energy dissipation mechanisms that were previously unattainable [6, 7]. It is now possible to investigate whether energy dissipation in simulated low plasma-$\beta$ environments ($\beta$ is the ratio of gas pressure to magnetic pressure [22]), is scale invariant, and if it occurs through

avalanches.

To address this question, we analyze magnetic energy dissipation in two distinct MHD simulations. More precisely, we identify regions of high magnetic gradients where ohmic dissipation is expected to occur and analyze their evolution. The two simulation setups we use in this study are displayed in figure 1. The first (top panel) is a state of the art Bifrost simulation [23], which is a 3D MHD solar-like atmosphere extending from the top of the convective zone to the lower-corona <fnref n="1" />. It is a quiet-Sun simulation (ch012023_by800) experiencing a network-like flux emergence episode. This simulation is designed to capture two dynamically different regimes. First, the convective zone along with the photosphere and lower chromosphere, hereafter CONV, spans approximately the lower third of the simulation and are in a high plasma-$\beta$ regime. The upper-chromosphere, the transition region and the beginning of the corona, hereafter ATMO, spans the upper two thirds of the simulation and are mostly in a low plasma-$\beta$ regime. We also study a completely different setup in a PLUTO simulation <fnref n="2" />, shown in figure 1 (bottom panel). This simulation, designed as a test case to study dissipative events, is a bipolar arcade embedded in a stratified low-$\beta$ atmosphere and slowly twisted at its base. Saying that these two simulations are different is an understatement. Table I lists various simulation characteristics and parameter values, to highlight the fact they differ markedly from one another, in terms of both numerical and physical setups.

We focus on locations with strong magnetic dissipation,

<footnote n="1">¹ Q. Noraz, et al. A&A (submitted)</footnote>

<footnote n="2">² A. Blaise, et al. in Prep</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 1. [Top] The ch012023_by800 Bifrost simulation. Magnetic field lines in the atmosphere are color coded with the Poynting flux. The photosphere is color coded with $B_z$ and the convective zone with vertical convective velocity. [Bottom] PLUTO simulation of a magnetic dipole initially in equilibrium with a stratified atmosphere. The plasma is then rotated at the foot of the dipole. In both panels, green isosurfaces enclose active zones (see text).</figcaption></figure>

which can be interpreted as current sheets. They reveal valuable insights about transfers between the magnetic energy and the internal energy. These locations are identified via the criterion:

$$ \frac{\|\nabla \times \mathbf{B}\|}{\|\mathbf{B}\|} > \frac{1}{\Delta}, \qquad (1) $$

where $\Delta$ corresponds to a length scale<fnref n="3" />. A low threshold

<table><caption>TABLE I. Characteristics of the MHD simulations studied in this work.</caption><thead><tr><th></th><th>PLUTO</th><th>Bifrost</th></tr></thead><tbody><tr><td>Magnetic Topology</td><td>Twisted bipolar arcade</td><td>Network-like flux emergence in the Quiet-Sun</td></tr><tr><td>Dissipative Scheme</td><td>Global constant resistivity and augmented local resistivity (YS-94 [7])</td><td>Explicit operator split into a weak global term and a localized hyperdif- fusive term [23]</td></tr><tr><td>Plasma-β range</td><td>10<sup>-4</sup> - 10<sup>2</sup></td><td>ATMO: 10<sup>-3</sup> - 10<sup>5</sup><br>CONV: 10<sup>-2</sup> - 10<sup>7</sup></td></tr><tr><td>Average Plasma-β</td><td>10<sup>-1</sup></td><td>ATMO: 10<sup>-1</sup><br>CONV: 10<sup>3</sup></td></tr><tr><td>Resolution</td><td>256x256x512</td><td>512x512x512</td></tr><tr><td>Cell dimension</td><td>70km</td><td>Horizontal: 20km<br>Vertical: 12-70km</td></tr></tbody></table>

old would yield only few very large zones, whereas a high threshold produces multiple small disjointed zones. The value of $\Delta = 0.01$Mm (Bifrost) and $\Delta = 0.035$Mm (PLUTO) which corresponds to $\Delta = \frac{1}{2}$ computational cell width is chosen to maximize simultaneously the number of identified zones and the range of their volumes. However, our results are robust to variations of $\Delta$, as discussed in Appendix A.

All cells for which the inequality 1 holds are defined as active cells. Active zones (hereafter zones) are defined using the Hoshen-Kopelman algorithm [24, 25] which identifies clusters defined as diagonally and horizontally connected active cells, displayed as green isosurfaces in figure 1. Each zone is assigned a label and an ancestry algorithm<fnref n="4" /> is applied to track their temporal evolution.

A first step in characterizing these dissipation zones is to calculate the joule heating

$$ Q_J = 4\pi\eta\frac{\mathbf{J}^2}{c^2} \qquad (2) $$

where $\eta$ is the magnetic diffusivity and $\mathbf{J}$ is the electric current density. It quantifies the transfer of energy from the magnetic field to internal energy. We compute the amount of magnetic energy dissipated by a zone by integrating $Q_J$ over its time-varying volume, over its lifetime. From this, we construct probability density functions (PDFs) for the dissipated magnetic energy and lifetime of the zones, as shown in figure 2. These frequency distributions are characterized by extended power-law tails, spanning more than eight orders of magnitudes in energy and two orders of magnitude in lifetimes. The best-fit exponents are listed in the figure. These values are roughly consistent with previous analyzes carried out on a similar

<footnote n="3">³ Note that equation 1 allows for cases where null points or points where $\|\mathbf{B}\|$ is small are tagged as active cells. These points, however, only represent a small (< 0.7%) subset of the active cells. In addition, those cells don't dissipate enough magnetic energy to impact the foregoing analysis.</footnote>

<footnote n="4">⁴ Our ancestry algorithm is analogous to the one described in [26]. However, in mergers and splitting, which account for 10% of events, the largest zone keeps its label whereas the smallest is considered ended, or created.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 2. Probability density function of the dissipated magnetic energy (A) and lifetime (B) of the active zones in the Bifrost simulation (ATMO and CONV) and the PLUTO simulation. The cited exponents and uncertainties are recovered with a random sample consensus method [28].</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 3. Evolution of the average number of active cells in dissipative zones in the Bifrost simulation, (ATMO and CONV), the PLUTO simulation, and in a Lu and Hamilton 3D model (AVALANCHE). We only retain averages with more than 50 zone samples for statistical significance ($\log_{10}(Lifetime) < 2$, for PLUTO for example).</figcaption></figure>

Bifrost simulation, using a completely distinct zone identification scheme [27], which strengthens our confidence in our identification algorithm. For all distributions, the low energies generally associated to zones comprised of a few cells are discarded from the fit performed with a random sample consensus method [28]. It is noteworthy that the low plasma-$\beta$ regimes of the two simulations, with completely different setups, are described by similar power laws in energies and lifetimes. This result not only confirms that dissipative events are scale invariant in MHD simulations, but also hints at the universality of this process.

Many physical mechanisms can lead to the existence of power-laws for the distribution of energy dissipated in localized regions. A paramount example can be found in MHD turbulence (akin to the CONV region in our Bifrost simulation), which is well-known for producing such power-laws [26, 29]. However, other generators of power laws exist in physical systems [30, 31], such as for instance avalanching processes and criticality. Asserting which frameworks can be used to describe this dissipation thus requires going beyond power-laws.

Avalanching systems have a characteristic evolution [32] controlled by a set of local rules of interaction between the smallest components of the system. Scale invariance and universality emerge naturally from such systems. A large event is simply a collection of smaller ones. Under slow driving, many such systems evolve autonomously towards a self-organized critical state, which naturally leads to scale invariance.

The Lu and Hamilton [33] avalanche model for solar flares (hereafter LH93) is an archetypal example. It is akin to a classical sandpile model [34], in which sand grains are added until the local slope exceeds a threshold and the pile becomes unstable. At this point, an unstable site will redistribute some sand to its neighbour to restore stability, but doing so may render its neighbours

unstable and trigger another redistribution, which can lead to an avalanche. In this model, scale invariance of avalanches emerges naturally from the stability and redistribution rules, which are purely local, and remains robust to changes in model parameters such as the magnitude of the driving [9], provided a good separation of timescales exists between energy loading and avalanching.

The local interactions will then determine the way unstable regions grow in these systems. This is typically assessed by critical exponents, which quantify the spatiotemporal evolution of these systems and define universality classes [32]. In this study, we focus on a spreading exponent [24, 35, 36] which characterizes the rate at which unstable zones grow over the lattice. It is recovered from the relation

$$N(t) = t^n, \qquad (3)$$

where $N$ is the number of new active cells at lifetime $t$, averaged over all dissipation events occurring over the duration of a given simulation, and $n$ is the spreading exponent.

We compute the spreading exponent $n$ for CONV and ATMO in the Bifrost simulation, the PLUTO simulation, and a Lu and Hamilton avalanche model (LH93). As shown in figure 3 in all cases, the average number of new active cells $N(t)$ grows as a power-law of time. These are fitted in logspace to measure their slope, directly yielding the spreading exponent $n$. Zones in the CONV region start on average with $\sim 3 (10^{0.5})$ active cells and the number of newly added cells per successive snapshot never exceeds 10. In contrast, zones in the ATMO region start on average with $\sim 25 (10^{1.4})$ active cells and the number of newly added active cells exhibits marked power-law growth. Similar evolution is observed
in the PLUTO simulation and the avalanche model<fnref n="5" />. The spreading exponent that quantifies these evolutions are listed in figure 3 as well as in table II, appendix B. Videos of the evolution of sample zones in ATMO and CONV are provided in supplementary material. In those videos, orange cubes are new active cells and green cubes are persisting active cells. In CONV, the new active cells show marked spatiotemporal variability and are vertically aligned, which is to be expected from a high plasma-β turbulent regime dominated by convective downflows. In ATMO, the new active cells are typically contiguous to active neighbours and the zone forms a sheet-like structure.

To summarize, externally driven homogeneous isotropic turbulence and SOC both generate scale invariance and both dissipate energy at the smallest length scales. However, in turbulence, energy is transfered from large scales to small scales, whereas the opposite can be argued for SOC where small-scale instabilities can destabilize larger areas through avalanches. In that regard, distinguishing between generators requires characterizing the spatiotemporal evolution of energy dissipation regions. In light of this, we assert that the evolution of zones in low plasma-β regimes (ATMO and PLUTO), are compatible with avalanche models, and confirm that zones in high plasma-β environments, such as CONV evolve differently than in low plasma-β environments. The much smaller value of *n* in CONV further underlines the distinct nature of the dissipative dynamics in the MHD-convective turbulent regime, where we do not expect avalanches to occur.

Systems near criticality have characteristic relations between their critical exponents [35]. These relations should hold for any system at or near criticality [35, 36]. Thus, we verify if one such relation holds for our models. This relation, detailed in appendix B, is based on another critical exponent, $\tau$, describing the probability distribution for zone sizes. Figure 4 shows that the low-$\beta$ regimes, ATMO and PLUTO, along with the LH93 model (AVALANCHE) all yield measured values for $\tau$ which are in acceptable agreement with the theoretical value. On the other hand, the measured value of $\tau$ for high plasma-$\beta$ (CONV) is in marked disagreement with the theoretical value $\tau^*$, showing that the CONV region is likely far from a critical regime. These reinforce the stark contrast between the high and low $\beta$ regimes in the Bifrost simulation and supporting the conjecture that the associated power-law reflect the turbulent state of the plasma, and does not indicate criticality.

The good agreement between the measured and computed critical exponents $\tau$ and $\tau^*$ (viz. Fig 4) in the low-$\beta$ region of the Bifrost and PLUTO simulations is

<figure><img src="image_7.png" /><figcaption>FIG. 4. Consistency test for the measured value of $\tau$ (Eq. B1) and its theoretical value $\tau^*$ (Eq. B2) where a $\pm 2\sigma$ interval is used for the error bars. The larger range on $\tau^*$ (gray box) for PLUTO is caused by smaller sample size for zones in that simulation.</figcaption></figure>

indicative of these simulations operating at or near a critical point. Since neither simulation is subjected to fine tuning via an external control parameter, this criticality must be of the self-organized variety. Self-organized criticality (SOC)[37, 38] is understood to materialize in slowly-driven open dissipative systems subjected to a local threshold instability [37, 39]. In the Bifrost simulation, the evolution of the ATMO region is driven by the magnetic footpoints being anchored in a turbulent convective layer with typical turnover timescale $\sim 10^3$ s; whereas in PLUTO a slow global twist is imposed at the footpoint of the pre-existing magnetic arcade, with timescale $\simeq 10^5$ s. This is much longer than the lifetime of all but the longest dissipation events observed in both simulations (viz. Fig. 2B), so that both systems can be legitimately considered slowly-driven. Activation of enhanced dissipation represents a form of threshold instability [6, 7], which affect plasma conditions locally and eventually brings the system back below threshold. The basic requirements for SOC are thus satisfied. Our analyzes thus support the conjecture that the low-$\beta$ regions of the solar chromosphere and lower corona are in a self-organized critical state, with conversion of magnetic energy to heat taking place through avalanches of small dissipative events [17], and scale-free avalanches mapping naturally on scale-free energy release in flares [39].

## ACKNOWLEDGMENTS

This research has been supported by the NSERC Discovery grant RGPIN-2024-04050, the Centre de Recherche en Astrophysique du Québec (CRAQ)S, the European Research Council through the Synergy Grant number 810218 (“The Whole Sun”, ERC-2018-SyG), the Centre National d'Etudes Spatiales (CNES) Solar Orbiter, the Institut National des Sciences de l'Univers (INSU) via the Programme National Soleil-Terre (PNST), the French Agence Nationale de la Recherche (ANR) project STORMGENESIS #ANR-22-CE31-0013-01, the Institut d'Astrophysique Spatiale (IAS) by the Research Council of Norway through its Centres

<footnote n="5">⁵ Time in the LH93 is measured in discrete iterations, but was rescaled to second here for comparison purposes in figure 3, which does not affect the value of the exponent</footnote>
<figure><img src="image_8.png" /><figcaption>FIG. 5. Variation of the distribution of the magnetic energy dissipation as a function of the threshold length scale $\Delta$ used to identify the studied dissipation zones. The results are shown for $\Delta = 0.01$ Mm (0.5 cells) in green, 0.004 Mm (0.2 cells) in red and 0.04 Mm (2 cells) in black. The power-law behaviour is essentially unaffected by this choice and $\Delta = 0.01$ Mm (0.5 cells) is used in this study.</figcaption></figure>

of Excellence, project number 262622, and Programme for Supercomputing, and by computing HPC and storage resources by GENCI thanks to the grant 2024-A0160410133.

## Appendix A: Threshold Variations

We have verified that the choice of the activation threshold length scale $\Delta$ in Eq. B3 does not affect significantly the power-law form and logarithmic slopes of the distribution of dissipation event sizes. As can be seen on Figure 5, only the small event end of the distribution is affected by variations of $\Delta$, while the power-law tail of large events, of primary interest in our study, remains unaffected over more than six orders of magnitude in dissipated energy.

## Appendix B: Critical Exponents

Critical exponents characterize the spatiotemporal evolution of systems near their critical point, and can be used collectively to define universality classes [32]. These exponents can be recovered by fitting associated distributions of event size measures and are related to one another via theoretically-derived numerical relationships [35]. One such exponent, $\tau$, quantifies the growth of an unstable region via the power-law relationship:

$$P(S) \sim S^{-\tau}, \qquad (B1)$$

where $P(S)$ is the probability that an active zone reaches a size $S$ in its lifetime. This exponent is directly accessible from simulation output, and in critical systems can be related to other critical exponents [35, 36]:

$$\tau^* = \frac{1 + n + 2\delta}{1 + n + \delta}, \qquad (B2)$$

where $n$ is the spreading exponent defined in equation 3, and $\delta$ characterizes the probability distribution $P(t)$ that a zone is still active at time $t$:

$$P(t) \sim t^{-\delta}. \qquad (B3)$$

The three critical exponents $n$, $\delta$ and $\tau$ are fitted on distributions constructed from simulation output, using the random consensus method [28]. Their numerical values are listed in table II along with the value of $\tau^*$ computed from Eq. B2. These are the data used to generate Fig. 4.

<table><caption>TABLE II. Measured and theoretical values for critical exponents in MHD simulations and an avalanche model.</caption><thead><tr><th>Critical Exponent</th><th>n - Eq.3</th><th>δ - Eq.B3</th><th>τ - Eq.B1</th><th>τ* - Eq. B2</th></tr></thead><tbody><tr><td>ATMO (Bifrost)</td><td>1.08 ± 0.05</td><td>1.8 ± 0.1</td><td>1.56 ± 0.03</td><td>1.46 ± 0.05</td></tr><tr><td>CONV (Bifrost)</td><td>0.22 ± 0.02</td><td>2.09 ± 0.08</td><td>2.24 ± 0.05</td><td>1.63 ± 0.04</td></tr><tr><td>PLUTO</td><td>0.9 ± 0.1</td><td>1.1 ± 0.2</td><td>1.16 ± 0.03</td><td>1.4 ± 0.1</td></tr><tr><td>Avalanche Model (LH93)</td><td>0.83 ± 0.02</td><td>1.04 ± 0.05</td><td>1.49 ± 0.03</td><td>1.36 ± 0.03</td></tr></tbody></table>

[1] L. de Arcangelis, C. Godano, E. Lippiello, and M. Nicodemi, Universality in solar flare and earthquake occurrence, Physical review letters 96, 051102 (2006).

[2] M. J. Aschwanden, T. D. Tarbell, R. W. Nightingale, C. J. Schrijver, C. C. Kankelborg, P. Martens, H. P. Warren, et al., Time variability of the “quiet” sun observed

with trace. ii. physical parameters, temperature evolution, and energetics of extreme-ultraviolet nanoflares, The Astrophysical Journal 535, 1047 (2000).

[3] T. Shimizu, Energetics and occurrence rate of active-region transient brightenings and implications for the heating of the active-region corona, Publications of the
Astronomical Society of Japan, v. 47, p. 251-263. 47, 251 (1995).

[4] N. B. Crosby, M. J. Aschwanden, and B. R. Dennis, Frequency distributions and correlations of solar x-ray flare parameters, Solar Physics 143, 275 (1993).

[5] Z. M. K. Abda, N. F. Ab Aziz, M. Z. A. Ab Kadir, and Z. A. Rhazali, A review of geomagnetically induced current effects on electrical power system: Principles and theory, IEEE Access 8, 200237 (2020).

[6] Ø. H. Færder, D. Nóbrega-Siverio, and M. Carlsson, A comparative study of resistivity models for simulations of magnetic reconnection in the solar atmosphere, Astronomy & Astrophysics 675, A97 (2023).

[7] Ø. H. Færder, D. Nóbrega-Siverio, and M. Carlsson, A comparative study of resistivity models for simulations of magnetic reconnection in the solar atmosphere-ii. plasmoid formation, Astronomy & Astrophysics 683, A95 (2024).

[8] P. Charbonneau, S. W. McIntosh, H.-L. Liu, and T. J. Bogdan, Avalanche models for solar flares (invited review), Solar Physics 203, 321 (2001).

[9] A. Strugarek, P. Charbonneau, R. Joseph, and D. Pirot, Deterministically driven avalanche models of solar flares, Coronal Magnetometry, 371 (2014).

[10] C. Thibeault, A. Strugarek, P. Charbonneau, and B. Tremblay, Forecasting solar flares by data assimilation in sandpile models. solar physics, arXiv preprint arXiv:2206.13583 (2022).

[11] M. L. Fuentes and J. A. Klimchuk, A nanoflare-based cellular automaton model and the observed properties of the coronal plasma, The Astrophysical Journal 828, 86 (2016).

[12] M. Paczuski and S. Boettcher, Universality in sandpiles, interface depinning, and earthquake models, Physical review letters 77, 111 (1996).

[13] D. L. Turcotte and B. D. Malamud, Landslides, forest fires, and earthquakes: examples of self-organized critical behavior, Physica A: Statistical Mechanics and its Applications 340, 580 (2004).

[14] L. Chitta, D. Pontin, E. Priest, D. Berghmans, E. Kraaikamp, L. Rodriguez, C. Verbeeck, A. Zhukov, S. Krucker, R. A. Cuadrado, et al., A magnetic avalanche as the central engine powering a solar flare, arXiv preprint arXiv:2503.12235 (2025).

[15] J. Reid, J. Threlfall, and A. W. Hood, Self-consistent nanoflare heating in model active regions: MHD avalanches, MNRAS 518, 1584 (2023).

[16] G. Cozzo, J. Reid, P. Pagano, F. Reale, and A. W. Hood, Coronal energy release by mhd avalanches-effects on a structured, active region, multi-threaded coronal loop, Astronomy & Astrophysics 678, A40 (2023).

[17] E. N. Parker, Nanoflares and the solar x-ray corona, Astrophysical Journal, Part 1 (ISSN 0004-637X), vol. 330, July 1, 1988, p. 474-479. 330, 474 (1988).

[18] E. T. Lu, Avalanches in Continuum Driven Dissipative Systems, Phys. Rev. Lett. 74, 2511 (1995).

[19] H. Isliker, A. Anastasiadis, D. Vassiliadis, and L. Vlahos, Solar flare cellular automata interpreted as discretized mhd equations, Astronomy and Astrophysics, v. 335, p. 1085-1092 (1998) 335, 1085 (1998).

[20] H.-L. Liu, P. Charbonneau, A. Pouquet, T. Bogdan, and S. McIntosh, Continuum analysis of an avalanche model for solar flares, Physical Review E 66, 056111 (2002).

[21] H. Lamarre, P. Charbonneau, and A. Strugarek, Energy definition and minimization in avalanche models for solar flares, Solar Physics 299, 13 (2024).

[22] E. Priest, Magnetohydrodynamics of the Sun (Cambridge University Press, 2014).

[23] B. V. Gudiksen, M. Carlsson, V. H. Hansteen, W. Hayek, J. Leenaarts, and J. Martínez-Sykora, The stellar atmosphere simulation code bifrost-code description and validation, Astronomy & Astrophysics 531, A154 (2011).

[24] L. F. Morales, P. Dmitruk, and D. O. Gómez, Energy dissipation in coronal loops: Statistical analysis of intermittent structures in magnetohydrodynamic turbulence, The Astrophysical Journal 894, 90 (2020).

[25] J. Hoshen and R. Kopelman, Percolation and cluster distribution. i. cluster multiple labeling technique and critical concentration algorithm, Physical Review B 14, 3438 (1976).

[26] V. Zhdankin, D. A. Uzdensky, and S. Boldyrev, Temporal analysis of dissipative structures in magnetohydrodynamic turbulence, The Astrophysical Journal 811, 6 (2015).

[27] C. Kanella and B. V. Gudiksen, Investigating 4d coronal heating events in magnetohydrodynamic simulations, Astronomy & Astrophysics 617, A50 (2018).

[28] O. Chum, J. Matas, and J. Kittler, Locally optimized ransac, in Pattern Recognition: 25th DAGM Symposium, Magdeburg, Germany, September 10-12, 2003. Proceedings 25 (Springer, 2003) pp. 236-243.

[29] G. Boffetta, V. Carbone, P. Giuliani, P. Veltri, and A. Vulpiani, Power laws in solar flares: self-organized criticality or turbulence?, Physical review letters 83, 4662 (1999).

[30] M. Aschwanden, Power Laws in Astrophysics: Self-Organized Criticality Systems (Cambridge University Press, 2024).

[31] D. Sornette, Critical phenomena in natural sciences: chaos, fractals, selforganization and disorder: concepts and tools (Springer Science & Business Media, 2006).

[32] G. Ódor, Universality classes in nonequilibrium lattice systems, Reviews of modern physics 76, 663 (2004).

[33] E. T. Lu, R. J. Hamilton, J. McTiernan, and K. R. Bromund, Solar flares and avalanches in driven dissipative systems, Astrophysical Journal, Part 1 (ISSN 0004-637X), vol. 412, no. 2, p. 841-852. 412, 841 (1993).

[34] P. Bak, C. Tang, and K. Wiesenfeld, Self-organized criticality: An explanation of the 1/f noise, Physical review letters 59, 381 (1987).

[35] M. A. Munoz, R. Dickman, A. Vespignani, and S. Zapperi, Avalanche and spreading exponents in systems with absorbing states, Physical Review E 59, 6175 (1999).

[36] L. F. Morales and P. Charbonneau, Scaling laws and frequency distributions of avalanche areas in a self-organized criticality model of solar flares, Geophysical research letters 35 (2008).

[37] H. J. Jensen, Self-organized criticality: emergent complex behavior in physical and biological systems, Vol. 10 (Cambridge university press, 1998).

[38] M. Aschwanden, Self-organized criticality in astrophysics: The statistics of nonlinear processes in the universe (Springer Science & Business Media, 2011).

[39] E. T. Lu, The Statistical Physics of Solar Active Regions and the Fundamental Nature of Solar Flares, APJL 446, L109 (1995).