* Paxton et al. (2019) Paxton B., et al., 2019, ApJS, 243, 10
  * Penny et al. (2008) Penny L. R., Ouzts C., Gies D. R., 2008, ApJ, 681, 554
  * Petrovic et al. (2005) Petrovic J., Langer N., Yoon S. C., Heger A., 2005, A&A, 435, 247
  * Plewa & Wlodarczyk (1993) Plewa T., Wlodarczyk K. J., 1993, Acta Astron., 43, 249
  * Podsiadlowski et al. (1992) Podsiadlowski P., Joss P. C., Hsu J. J. L., 1992, ApJ, 391, 246
  * Podsiadlowski et al. (2006) Podsiadlowski P., Morris T. S., Ivanova N., 2006, in Kraus M., Miroshnichenko A. S., eds, Astronomical Society of the Pacific Conference Series Vol. 355, Stars with the B[e] Phenomenon. p. 259
  * Pols et al. (1998) Pols O. R., Schröder K.-P., Hurley J. R., Tout C. A., Eggleton P. P., 1998, MNRAS, 298, 525
  * Qian et al. (2006) Qian S., Yang Y., Zhu L., He J., Yuan J., 2006, Ap&SS, 304, 25
  * Qian et al. (2007) Qian S. B., Yuan J. Z., Liu L., He J. J., Fernández Lajús E., Kreiner J. Z., 2007, MNRAS, 380, 1599
  * Raghavan et al. (2010) Raghavan D., et al., 2010, ApJS, 190, 1
  * Ramírez-Agudelo et al. (2017) Ramírez-Agudelo O. H., et al., 2017, A&A, 600, A81
  * Ramírez-Tannus et al. (2021) Ramírez-Tannus M. C., et al., 2021, A&A, 645, L10
  * Raucq et al. (2017) Raucq F., Gosset E., Rauw G., Manfroid J., Mahy L., Mennekens N., Vanbeveren D., 2017, A&A, 601, A133
  * Rauw et al. (2019) Rauw G., Nazé Y., Campos F., 2019, A&A, 627, A2
  * Riley et al. (2020) Riley J., Mandel I., Marchant P., Butler E., Nathaniel K., Neijssel C., Shortt S., Vigna-Gomez A., 2020, arXiv e-prints, p. arXiv:2010.00002
  * Ritter (1988) Ritter H., 1988, A&A, 202, 93
  * Rucinski (1999) Rucinski S. M., 1999, Acta Astron., 49, 341
  * Rucinski et al. (2007) Rucinski S. M., Pribulla T., van Kerkwijk M. H., 2007, AJ, 134, 2353
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Sana et al. (2013) Sana H., et al., 2013, A&A, 550, A107
  * Sana et al. (2014) Sana H., et al., 2014, ApJS, 215, 15
  * Sana et al. (2017) Sana H., Ramírez-Tannus M. C., de Koter A., Kaper L., Tramper F., Bik A., 2017, A&A, 599, L9
  * Sanyal et al. (2015) Sanyal D., Grassitelli L., Langer N., Bestenlehner J. M., 2015, A&A, 580, A20
  * Sanyal et al. (2017) Sanyal D., Langer N., Szécsi D., -C Yoon S., Grassitelli L., 2017, A&A, 597, A71
  * Schneider et al. (2016) Schneider F. R. N., Podsiadlowski P., Langer N., Castro N., Fossati L., 2016, MNRAS, 457, 2355
  * Schneider et al. (2018a) Schneider F. R. N., et al., 2018a, Science, 359, 69
  * Schneider et al. (2018b) Schneider F. R. N., et al., 2018b, A&A, 618, A73
  * Schneider et al. (2019) Schneider F. R. N., Ohlmann S. T., Podsiadlowski P., Röpke F. K., Balbus S. A., Pakmor R., Springel V., 2019, Nature, 574, 211
  * Schootemeijer et al. (2019) Schootemeijer A., Langer N., Grin N. J., Wang C., 2019, A&A, 625, A132
  * Selam (2004) Selam S. O., 2004, A&A, 416, 1097
  * Sen et al. (in prep) Sen K., et al., in prep., A&A
  * Shematovich et al. (2017) Shematovich V. I., Kaigorodov P. V., Bisikalo D. V., Fateeva A. M., 2017, Astronomy Reports, 61, 1038
  * Soberman et al. (1997) Soberman G. E., Phinney E. S., van den Heuvel E. P. J., 1997, A&A, 327, 620
  * Song et al. (2013) Song H. F., Maeder A., Meynet G., Huang R. Q., Ekström S., Granada A., 2013, A&A, 556, A100
  * Song et al. (2016) Song H. F., Meynet G., Maeder A., Ekström S., Eggenberger P., 2016, A&A, 585, A120
  * Stȩpień (2009) Stȩpień K., 2009, MNRAS, 397, 857
  * Stepien (2006) Stepien K., 2006, Acta Astron., 56, 199
  * Stickland et al. (1994) Stickland D. J., Koch R. H., Pachoulakis I., Pfeiffer R. J., 1994, The Observatory, 114, 107
  * Stępień & Gazeas (2012) Stępień K., Gazeas K., 2012, Acta Astron., 62, 153
  * Suijs et al. (2008) Suijs M. P. L., Langer N., Poelarends A. J., Yoon S. C., Heger A., Herwig F., 2008, A&A, 481, L87
  * Sukhbold et al. (2016) Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H.-T., 2016, ApJ, 821, 38
  * Sybesma (1985) Sybesma C. H. B., 1985, A&A, 142, 171
  * Sybesma (1986a) Sybesma C. H. B., 1986a, A&A, 159, 108
  * Sybesma (1986b) Sybesma C. H. B., 1986b, A&A, 168, 147
  * Szécsi (2017) Szécsi D., 2017, Proceedings of Science, 2017 [arXiv:1710.05655], PoS(MULTIF2017)065
  * Szécsi et al. (2015) Szécsi D., Langer N., Yoon S.-C., Sanyal D., de Mink S., Evans C. J., Dermine T., 2015, A&A, 581, A15
  * Szymanski et al. (2001) Szymanski M., Kubiak M., Udalski A., 2001, Acta Astron., 51, 259
  * Terrell et al. (2003) Terrell D., Munari U., Zwitter T., Nelson R. H., 2003, AJ, 126, 2988
  * Ulrich & Burger (1976) Ulrich R. K., Burger H. L., 1976, ApJ, 206, 509
  * Urushibata et al. (2018) Urushibata T., Takahashi K., Umeda H., Yoshida T., 2018, MNRAS, 473, L101
  * Vanbeveren et al. (1998) Vanbeveren D., De Loore C., Van Rensbergen W., 1998, A&ARv, 9, 63
  * Vanbeveren et al. (2013) Vanbeveren D., Mennekens N., Van Rensbergen W., De Loore C., 2013, A&A, 552, A105
  * Vaz et al. (1995) Vaz L. P. R., Andersen J., Rabello Soares M. C. A., 1995, A&A, 301, 693
  * Vigna-Gómez et al. (2019) Vigna-Gómez A., Justham S., Mandel I., de Mink S. E., Podsiadlowski P., 2019, ApJ, 876, L29
  * Vink et al. (2001) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, A&A, 369, 574
  * Webbink (1976) Webbink R. F., 1976, ApJ, 209, 829
  * Webbink (1977) Webbink R. F., 1977, ApJ, 215, 851
  * Wellstein et al. (2001) Wellstein S., Langer N., Braun H., 2001, A&A, 369, 939
  * Wilson (2001) Wilson R. E., 2001, Information Bulletin on Variable Stars, 5076, 1
  * Woosley & Bloom (2006) Woosley S. E., Bloom J. S., 2006, ARA&A, 44, 507
  * Woosley et al. (2007) Woosley S. E., Blinnikov S., Heger A., 2007, Nature, 450, 390
  * Wu et al. (2020) Wu S., Everson R. W., Schneider F. R. N., Podsiadlowski P., Ramirez-Ruiz E., 2020, ApJ, 901, 44
  * Yakut & Eggleton (2005) Yakut K., Eggleton P. P., 2005, ApJ, 629, 1055
  * Yang et al. (2019) Yang Y., Yuan H., Dai H., 2019, AJ, 157, 111
  * Yasarsoy & Yakut (2014) Yasarsoy B., Yakut K., 2014, New Astron., 31, 32
  * Yoon & Langer (2005) Yoon S. C., Langer N., 2005, A&A, 443, 643
  * Yoon et al. (2006) Yoon S. C., Langer N., Norman C., 2006, A&A, 460, 199
  * Zhang et al. (2020) Zhang X.-D., Qian S.-B., Liao W.-P., 2020, MNRAS, 492, 4112
  * Zhao et al. (2014) Zhao E., Qian S., Li L., He J., Liu L., Wang J., Zhang J., 2014, New Astron., 26, 112
  * Zhao et al. (2018) Zhao E.-G., Qian S.-B., Zejda M., Zhang B., Zhang J., 2018, Research in Astronomy and Astrophysics, 18, 059
  * de Mink & Mandel (2016) de Mink S. E., Mandel I., 2016, MNRAS, 460, 3545
  * de Mink et al. (2007) de Mink S. E., Pols O. R., Hilditch R. W., 2007, A&A, 467, 1181
  * de Mink et al. (2009) de Mink S. E., Cantiello M., Langer N., Pols O. R., Brott I., Yoon S. C., 2009, A&A, 497, 243
  * du Buisson et al. (2020) du Buisson L., et al., 2020, MNRAS, 499, 5941
  * van Marle, A. J. et al. (2008) van Marle, A. J. Langer, N. Yoon, S.-C. García-Segura, G. 2008, A&A, 478, 769


# Theoretical and Observational Evidence for Coriolis Effects in Coronal
Magnetic Fields Via Direct Current Driven Flaring Events

Darryl Z. Seligman Department of the Geophysical Sciences, University of
Chicago, Chicago, IL 60637, USA  Leslie A. Rogers Department of Astronomy and
Astrophysics, University of Chicago, Chicago, IL 60637, USA  Adina D.
Feinstein NSF Graduate Research Fellow  Department of Astronomy and
Astrophysics, University of Chicago, Chicago, IL 60637, USA  Mark R. Krumholz
Research School of Astronomy and Astrophysics, Australian National University,
Canberra, ACT 2611, Australia  ARC Centre of Excellence for Astronomy in Three
Dimensions (ASTRO-3D), Canberra, ACT 2611, Australia  James R. Beattie
Research School of Astronomy and Astrophysics, Australian National University,
Canberra, ACT 2611, Australia  Christoph Federrath Research School of
Astronomy and Astrophysics, Australian National University, Canberra, ACT
2611, Australia  ARC Centre of Excellence for Astronomy in Three Dimensions
(ASTRO-3D), Canberra, ACT 2611, Australia  Fred C. Adams Physics Department,
University of Michigan, Ann Arbor, MI 48109  Astronomy Department, University
of Michigan, Ann Arbor, MI 48109  Marco Fatuzzo  Physics Department, Xavier
University, Cincinnati, OH 45207  Maximilian N. Günther ESA Research Fellow
European Space Agency (ESA), European Space Research and Technology Centre
(ESTEC), Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

###### Abstract

All stars produce explosive surface events such as flares and coronal mass
ejections. These events are driven by the release of energy stored in coronal
magnetic fields, generated by the stellar dynamo. However, it remains unclear
if the energy deposition in the magnetic fields is driven by direct or
alternating currents. Recently, we presented observational measurements of the
flare intensity distributions for a sample of $\sim 10^{5}$ stars across the
main sequence observed by TESS, all of which exhibited power-law distributions
similar to those observed in the Sun, albeit with varying slopes. Here we
investigate the mechanisms required to produce such a distribution of flaring
events via direct current energy deposition, in which coronal magnetic fields
braid, reconnect, and produce flares. We adopt a topological model for this
process which produces a power-law distribution of energetic flaring events.
We expand this model to include the Coriolis effect, which we demonstrate
produces a shallower distribution of flare energies in stars that rotate more
rapidly (corresponding to a weaker decline in occurrence rates toward
increasing flare energies). We present tentative evidence for the predicted
rotation-power-law index correlation in the observations. We advocate for
future observations of stellar flares that would improve our measurements of
the power-law exponents, and yield key insights into the underlying dynamo
mechanisms that underpin the self-similar flare intensity distributions.

emerging flux tubes – solar dynamo – solar magnetic reconnection – optical
flares

††facilities: TESS (Ricker et al., 2015), Kepler (Borucki et al., 2010; Howell
et al., 2014)††software: lightkurve (Lightkurve Collaboration et al., 2018),
matplotlib (Hunter, 2007), astropy (Astropy Collaboration et al., 2013; Price-
Whelan et al., 2018), astroquery(Ginsburg et al., 2019), numpy (Van Der Walt
et al., 2011; Harris et al., 2020), stella (Feinstein et al., 2020a), emcee
(Foreman-Mackey et al., 2013)

##  1 Introduction

Since the first recorded observation of a Solar superflare in 1859 by
Carrington (1859), the mechanisms driving explosive flaring in stars have been
the subject of detailed inquiry from plasma physicists and astronomers alike.
The origin of flares is closely related to the still-open problem of what
mechanism is responsible for heating the Solar corona (Withbroe & Noyes,
1977). In the working picture that has emerged since 1859, it has become clear
that magnetic fields are an important component of the heating (Golub &
Pasachoff, 1997; Parker, 1989; Schrijver & Zwaan, 2000). The concentration and
twisting of magnetic field lines can lead to the release of energy via the
process of magnetic reconnection in flaring events (Sturrock et al., 1984;
Kulsrud, 1998; Priest & Forbes, 2000). In this study, we are interested in
understanding the role of the Coriolis force in the dynamo of other stars, and
to accomplish this we apply theoretical insights about the heating of the
Solar corona and the energy mechanisms driving Solar flares to other stars.

The role of the Coriolis force in the Solar dynamo produces observational
signatures in the inclination, twist and location of coronal loops and active
regions. Loops tend to be inclined with respect to the latitudinal line (Hale
et al., 1919; Howard, 1993; Choudhuri & Gilman, 1987; Choudhuri & D’Silva,
1990; Howard, 1993; D’Silva & Choudhuri, 1993; Longcope & Fisher, 1996;
Longcope et al., 1998; Weber et al., 2011, 2013). Magnetic fields in active
regions tend to twist in the opposite direction in the North and South
hemisphere (Seehafer, 1990; Pevtsov et al., 1995; Abramenko et al., 1997; Bao
& Zhang, 1998; Longcope & Welsch, 2000; Longcope & Pevtsov, 2003; Seligman et
al., 2014; Manek et al., 2018; Manek & Brummell, 2021). Faster rotating
sunspots are more likely to be towards the equator and slower rotating spots
towards the pole (D’Silva & Howard, 1995). Stars spin down as they age (Noyes
et al., 1984; Soderblom, 2010; Notsu et al., 2013; Candelaresi et al., 2014;
Doyle et al., 2018; Ilin et al., 2019; Doyle et al., 2019, 2020). Recent work
has considered the role of the coriolis force in the dynamo of other stars
with different rotation rates (Holzwarth, 2007; Kitchatinov & Olemskoy, 2015;
Weber & Browning, 2016; Browning et al., 2016; Solanki et al., 2006; Parker,
2009; Rempel, 2011). In this paper we consider an additional signature of the
Coriolis force: its impact on the distribution of flares.

Observations at different wavelengths have demonstrated that the distribution
of peak intensity of Solar flares follow power-laws (Drake, 1971; Datlowe et
al., 1974; Dennis, 1985; Lin et al., 1984), where

$\frac{dN}{dE}\sim E^{-\alpha}\,,$ (1)

and $N$ is the number of flares, $E$ is the total energy released in the flare
(where the flare magnitude $m\propto\log E$), and the power-law exponent falls
roughly in the range of $\alpha\sim 1.25-1.5$. Wheatland (2000) found that
typical flare frequency distribution (FFD) indices did not vary for individual
active regions on the Sun. However, it is important to note that converting
from the intensity to energy requires assumptions about the geometry and
physical conditions of the flaring region, the mechanisms producing extremely
energetic photons (Lee et al., 1993; Brown et al., 1998), and the flare height
as a function of amplitude (Mitra-Kraev & Benz, 2001).

A useful parameter to describe the potential mechanisms driving coronal
heating and the coupling of the convective interior and the corona is the
plasma $\beta$ parameter,

$\beta=2\frac{c_{s}^{2}}{v_{A}^{2}}\,,$ (2)

where $c_{s}$ is the sound speed and $v_{A}$ is the Alfvén speed. When
$\beta>1$ acoustic modes dominate the dynamics and transport of energy, and
when $\beta<1$ Alfvénic modes dominate. This is useful for identifying the
importance of electrodynamic coupling relative to mechanical coupling between
the inner and outer atmosphere of the Sun. The electrodynamic coupling will
dominate if $\beta<1$ in the outer atmosphere and $\beta\sim 1$ in the inner
atmosphere, where mechanical dynamics such as convection or differential
rotation can couple to Alfvénic perturbations in the corona and drive the
heating (Ionson, 1985). This is generally believed to be the case for the Sun,
and similarly for other stars.

Two broad mechanisms that have been proposed to explain the heating of the
Solar corona (Withbroe & Noyes, 1977) are direct current (DC) and alternating
current (AC). Magnetic stress or DC heating dominates when large scale
subsurface fluid motions have timescales that are much longer than the
Alfvénic coronal crossing timescale, $\tau_{A}=L_{c}/v_{A}$ (where $L_{c}$ is
the radial extent of the corona) (Ionson, 1982), while AC heating dominates in
the opposite limit.

The fundamental idea behind DC models originates in the work of Parker (1972),
who demonstrated that a large scale magnetic field could not be in hydrostatic
equilibrium if field lines were topologically braided or knotted. The
nonexistence of an equilibrium is robust to any pressure perturbations applied
along the field line. Parker concluded that braided or knotted flux tubes
produce rapid dissipation and merging of field lines into a one-dimensional
topology. Parker (1983) connected this simple model to heating of the Solar
corona via DC currents, whereby slow random walks in the footpoint positions
generate magnetic braids that must reconnect and release energy. This braiding
mechanism for heating of the Solar corona was extended by many authors (Mikic
et al., 1989; Berger, 1993; Parker, 1994; Longcope & Sudan, 1994; Galsgaard &
Nordlund, 1996; Berger & Asgari-Targhi, 2009; Berger et al., 2015).

In the DC regime, a simplified view is that reconnection is triggered when the
angle, $\theta$, between neighboring magnetic field vectors is greater than a
critical angle $\theta_{c}$, which may trigger smaller events in neighboring
field lines (Sturrock et al., 1984; Porter et al., 1987; Parker, 1988;
Sturrock et al., 1990; Berger, 1993; Krucker & Benz, 2000), sometimes referred
to as nanoflares. These were recently observed by Antolin et al. (2021), who
verified that these could be explained as reconnection events at small angles.
Parker (1988) conjectured that $\theta_{c}\sim 30^{\circ}$, although Dahlburg
et al. (2005) argued that $\theta_{c}\sim 45^{\circ}$ with a more detailed
analysis that included secondary instability. It is important to note,
however, that models of coronal heating such as that presented by van
Ballegooijen (1986) required reconnection events even at small $\theta_{c}$ in
a cascade of magnetic energy transport. This model was furthered by Cargill
(1994), who presented a model of an active region as hundreds of small
elemental flux loops randomly heated by nanoflares. In any case, when
$\theta<\theta_{c}$, magnetic reconnection proceeds slowly, and magnetic
energy is deposited in the form of braided fields quadratically in time
(Parker, 1983; Moffatt & Tsinober, 1990; Berger, 1993). The combination of
twisting and braiding of field lines could explain the flare frequency
spectrum observed in active regions (Zirker & Cleveland, 1993a, b).

For AC models, the timescale of the footpoint motions is of order or shorter
than the coronal crossing timescale, so footpoint motions excite Alfvénic
modes that travel back and forth along coronal magnetic field lines (van
Ballegooijen et al., 2011; Asgari-Targhi & van Ballegooijen, 2012; Asgari-
Targhi et al., 2013; van Ballegooijen et al., 2014). This regime is driven by
smaller scale surface convection mechanisms (Ionson, 1985). Asgari-Targhi et
al. (2014) demonstrated that the AC model was consistent with observations of
non-thermal widths of coronal emission of Fe xii, Fe xiii, Fe xv, and Fe xvi
from the Extreme-ultraviolet Imaging Spectrometer on the Hinode spacecraft.
However, it is not clear that AC models can explain the existence of
exceptionally hot $>5$ MK coronal loops (Asgari-Targhi et al., 2015), and it
has been speculated that DC events are important there.

The heating of the Solar corona has been linked to the theory of self-
organized criticality (Bak et al., 1987, 1988), which describes dissipative
dynamical systems that remain in a critical state with no intrinsic length or
time scale. The theory requires a local instability mechanism that can trigger
neighboring instabilities. This avalanching mechanism produces energetic
events at all length scales (Kadanoff et al., 1989; Babcock & Westervelt,
1990). Applications of the theory have been hypothesized in turbulence,
percolation systems (Turcotte, 1999), neuroscience (Ribeiro et al., 2010;
Hesse & Gross, 2014), landslides (Bak et al., 1990; Turcotte et al., 2002),
atmospheric dynamics (Grieger, 1992; Andrade et al., 1998), astrophysical
accretion disks (Dendy et al., 1998), traffic patterns (Nagel & Herrmann,
1993), evolution (Bak & Sneppen, 1993), extinction events (Newman, 1996),
financial markets (Bak et al., 1997), earthquakes (Gutenberg & Richter, 1956;
Bak & Tang, 1989; Sornette & Sornette, 1989; Olami et al., 1992; Carlson &
Langer, 1989), and Conway’s game of Life (Bak et al., 1989). Lu & Hamilton
(1991) proposed that the Solar coronal magnetic field is also in a self-
organized critical state to explain the power-law observed in the magnitude of
Solar flares (Lu et al., 1993; Crosby et al., 1993; Aschwanden et al., 1998;
Charbonneau et al., 2001; de Arcangelis et al., 2006). This approach has been
powerful in prediction of extreme flares (Morales & Santos, 2020).

In the DC picture of coronal heating, the slow buildup of braided fields
provides all of the requisites for a self-organized critical system. Twisted
coronal fields are generated via dynamo mechanisms in the fluid-dominated
interior (Charbonneau, 2010), via convective and Coriolis driven vortical
subsurface plasma flows (Parker, 1955; Moffatt, 1978; Longcope et al., 1998;
Seligman et al., 2014). Prior & MacTaggart (2016) suggested that active Solar
regions could be formed via the injection of rising magnetic field topologies
that were only braided, and not twisted. MacTaggart et al. (2021) demonstrated
numerically that active regions are created by the emergence of a large flux
tube of pre-twisted magnetic fields by examining the evolution of the
topological quantity magnetic winding in the emergence of active regions.
Berger & Asgari-Targhi (2009) and Berger et al. (2015) demonstrated that a
topological model that included braiding and reconnection of coronal fields
exhibited power-law distributions of energetic events.

It is important to note that while there exist physical connections between DC
heating models and self-organized criticality, the theory is not the only way
to create a power-law distribution of energetic events (Rosner & Vaiana, 1978;
Litvinenko, 1996). For example, Newman & Sneppen (1996) demonstrated a
dynamical system driven by coherent noise could arrive at a similar stationary
state characterized by power-law distributions of avalanches, but without
maintaining the “critical” state. MHD turbulence invoked in the AC regime can
produce flare-like energy occurrence distributions without relying on the
theory of self-organized criticality (Longcope & Sudan, 1994; Einaudi et al.,
1996; Galsgaard & Nordlund, 1996; Dmitruk & Gómez, 1997; Galtier & Pouquet,
1998; Georgoulis et al., 1998; Einaudi & Velli, 1999; Galtier, 1999).

To date, most studies of heating and flaring mechanisms have focused on the
Sun. While this is obviously the system for which we can obtain the richest
and highest-quality data, it remains a sample of one. Studies of flaring in
other stars can therefore offer a unique and valuable perspective. Flare-like
X-ray emissions from other stars exhibit power-law energy distributions with
similar indices to those seen in the Sun (Shakhovskaia, 1989; Osten & Brown,
1999; Audard et al., 2000). Aschwanden & Güdel (2021) found a power-law
dependence of energies for optical flares observed with Kepler (Davenport,
2016). The TESS mission provided 2-minute cadence light curves for
$\sim$200,000 stars, which allowed for the identification of a statistically
significant sample of flaring events (Günther et al., 2020). Feinstein et al.
(2021) demonstrated that the these flares follow a power-law distribution of
intensity (with slopes $\alpha^{\prime}\approx 0.9-1.5$) for all main sequence
stars.111$\alpha^{\prime}$ indicates the slope analogous to $\alpha$ in
Equation 1, for the normalized intensity distribution. These newly measured
values are close to the median of previously measured slopes (Figure 3 in
Feinstein et al. (2021); Shibayama et al., 2013; Günther et al., 2020; Ilin et
al., 2019; Lin et al., 2019; Howard et al., 2019; Yang & Liu, 2019; Feinstein
et al., 2020b; Raetz et al., 2020; Ilin et al., 2021; Aschwanden & Güdel,
2021). In this paper, we attempt to interpret these power-law distributions
under the framework of the DC heating mechanism.

This paper is organized as follows. In §2, we review the DC braiding model of
reconnection events presented by Berger & Asgari-Targhi (2009) and Berger et
al. (2015). We expand the braiding model to include a Coriolis-driven bias in
handedness of braids injected, and calculate the resulting distributions of
energetic events. In §3, we present observations of flare frequency
distributions (FFDs) for slow and fast rotating stars, and show that they are
consistent with the analytic predictions. In §4 we conclude and outline future
observational and theoretical work.

##  2 Generalised Topological Braiding Model

Reconnection is at its base a phenomenon whereby (anomalous) resistivity
dissipates the currents that sustain magnetic fields in a plasma, leading to a
violation of flux freezing and allowing rearrangement of the magnetic
topology. Regions where reconnection occurs must have a current flowing
through them (since otherwise there would be nothing to dissipate), and thus
necessarily have

$|\,\vec{\nabla}\times\vec{B}\,|\neq 0\,.$ (3)

This result follows from consideration of Ampére’s Law,

$\vec{\nabla}\times\vec{B}=\frac{1}{c}\left(4\pi\vec{J}+\frac{\partial\vec{E}}{\partial
t}\right)\,.$ (4)

In a plasma that is overall electrically neutral, the electric field in the
rest frame vanishes, $|\vec{E}|=0$, so that $\partial\vec{E}/\partial t=0$. It
is worth noting that in non-inertial frames, $\vec{E}=\vec{B}\times\vec{V}/c$,
where $\vec{V}$ is the plasma velocity and c is the speed of light. However,
$\partial\vec{E}/\partial t<<1$, because in a non-relativistic plasma,
$|\vec{V}|/c<<1$. The presence of a current, $|\vec{J}|\neq 0$, thus implies a
non-vanishing curl in the magnetic field, $|\vec{\nabla}\times\vec{B}|\neq 0$.
This condition, in turn, implies that the reconnection rate — and the flaring
rate — must be connected to the distribution of the curl of the magnetic
field. This insight motivates the analysis that follows.

In this section, we review the braiding model for DC coronal heating via
reconnection events presented in Berger & Asgari-Targhi (2009) and further
developed in Berger et al. (2015). We then generalise the braiding model to
include a Coriolis-induced handedness bias. Before we begin, it is worth
noting that other authors have investigated the role of helicity and braiding
in active regions and found that both are important for supplying energy to
the corona (Longcope et al., 2007; Liu et al., 2014).

###  2.1 Estimate of Free Magnetic Energy

We consider a small patch of a stellar atmosphere. We define height along the
$\hat{z}$ direction, defined such that $z=0$ is the base of the fields and
$z=L$ represents some radial extent into the corona. At $z=0$, which is
nominally at the photosphere, although the exact position is irrelevant,
subsurface convection provides stochastic forcing of the position for the
magnetic flux tubes.

As was introduced in Section 2 of Berger & Asgari-Targhi (2009), we consider a
braid consisting of $N$ flux tubes or magnetic field lines that permeate the
volume. The set of flux tubes is characterized by the number of crossing
points, $C$, at which magnetic flux tubes cross. Each flux tube has a number
of crossings, $n_{\rm C}$, on average, given by

$n_{\rm C}=\bigg{(}\frac{2C}{N}\,\bigg{)}\,,$ (5)

so the typical height between crossing junctures is

$\delta z=\frac{L}{n_{\rm C}}\,,$ (6)

as was demonstrated in Equation (1) of Berger & Asgari-Targhi (2009). Let $D$
be the typical diameter of a flux tube. Note that this quantity is related to
typical separations between magnetic field lines, and is given by the
resistivity of the ambient medium. If two of the tubes wrap around each other
between heights of $z$ and $z+\delta z$, the distance that they travel in the
horizontal direction is roughly $\delta\ell=\pi D/2$. Therefore the ratio
between the perpendicular and parallel components of the magnetic field with
respect to the radial direction, $B_{\bot}$ and $B_{\parallel}$, is
approximately given by

$\bigg{(}\,\frac{\delta\ell}{\delta
z}\bigg{)}\,\sim\bigg{(}\,\frac{B_{\bot}}{B_{\parallel}}\bigg{)}\,\sim\bigg{(}\,\frac{\pi
D}{NL}\,\bigg{)}\,C\,.$ (7)

This step is important, because it demonstrates that the perpendicular
magnetic field is directly proportional to the number of crossings in a flux
tube. Berger (1993) demonstrated that the free magnetic energy density is
proportional to the square of the perpendicular or transverse magnetic field,
so that the free energy per unit volume of the braided field obeys the scaling

$E_{\rm Free}\propto C^{2}\,,$ (8)

as given by Equation (5) of Berger & Asgari-Targhi (2009).

Figure 1: A schematic diagram of three braided flux tubes that cause an
incoherent reconnection event - resulting in the merging of opposite signed
sequences - that releases energy. In this setup, the two sequences on the left
hand side are braided in opposing directions. After the reconnection event,
which removes the interchange, the two sequences unwind and release the
magnetic energy density that was stored in transverse components of the
fields. This figure is similar to Figure 4 in Berger & Asgari-Targhi (2009).

###  2.2 Evolution to Equilibrium

Berger & Asgari-Targhi (2009) considered a topological model of the braiding
of flux tubes to explain the distribution of energy in reconnection events
observed in the Solar corona. In the model, there are two defining topological
quantities for a nest of individual magnetic flux tubes. The first is the
winding number $w$ of an individual sequence describing the braiding of two
flux tubes around each other. This quantity encodes both the number of
crossings, where $|w|=C$, and the handedness of the braiding itself. An
interchange represents a region where a third flux tube crosses one of the two
original flux tubes, and is inserted into the original braiding strand. These
interchanges are the sites of reconnection events.

In Figure 1, we present a schematic diagram of the topological model, similar
to Figure 4 in Berger & Asgari-Targhi (2009). The left hand side shows a braid
consisting of two sequences and one interchange. The primary flux tube is
shown in red, and takes part in both of the sequences. The blue flux tube is
braided around the red flux tube from the base where $z=0$, and these two
create a sequence with $w_{1}=-4$, of negative handedness. Above the sequence,
the blue flux tube is removed from the braid in an interchange with the yellow
flux tube. The yellow flux tube then creates a sequence where $w_{2}=+4$ with
the red flux tube above the interchange, of positive handedness.

Berger & Asgari-Targhi (2009) considered a model in which the coronal magnetic
field consisted of many of these types of braided flux tubes, nested together.
The defining feature of the model is that a reconnection event is modeled as
the removal of an interchange. This process is demonstrated in the middle
portion of Figure 1, where the interchange is removed and the blue and yellow
flux tubes merge on the right hand side. After the removal of the interchange,
the two sequences above and below the interchange are able to merge. The final
winding, $w_{f}$ of the resulting single sequences is given by
$w_{f}=w_{1}+w_{2}$. In the configuration shown in the schematic figure, the
equal and opposite windings above and below the interchange cancel each other
out completely. In other words, the entire braid completely unwinds itself.
This behavior is shown on the right hand side of the figure, where the three
flux tubes are now relaxed to an equilibrium where there is no braiding.
Because of the results in the previous subsection, and namely Equation (8),
the free magnetic energy that is released in the reconnection event encoded by
the removal of the interchange is given by the number of winds that are
undone.

Remarkably, Berger & Asgari-Targhi (2009) presented a method to evolve this
simplistic model through time to an equilibrium state, which we summarize
here. At some arbitrary initial time $t=0$, consider a neighborhood in the
corona containing an integer number of braiding sequences, $\Lambda$, which by
construction contains $\Lambda-1$ interchanges. The system then evolves
through a series of idealised discrete time steps. During any given time step,

  * •

1 sequence is added,

  * •

1 interchange is added,

  * •

1 reconnection event occurs, removing 1 pre-existing interchange. Sequences on
either side of the interchange merge.

Berger & Asgari-Targhi then investigate what steady-state distribution of
windings $w$ such an evolutionary sequence produces. At any point, the number
of sequences where $C=w$ is given by $\xi(w)$. Therefore,

$\Lambda=\sum_{w=-\infty}^{\infty}\xi(w)\,,$ (9)

as defined with different notation in Equation (7) of Berger & Asgari-Targhi
(2009). The probability that there will be a sequence of length $w$ is given
by the normalized probability function, $f(w)$, which is defined in Section
5.2 of Berger & Asgari-Targhi (2009) as

$f(w)=\frac{\xi(w)}{\Lambda}\,.$ (10)

Next, the probability of adding a new sequence of $w$ at each time step is
given by a function defined to be $p(w)$. At each time step in the algorithm,
the probability function evolves to a new value through the mapping

$f(w)\rightarrow f(w)+\delta f(w)\,,$ (11)

defined by the three events that are allowed to happen in each time step
itemized above. Therefore, this process can be described at each time step by
the following equation,

$\begin{split}\Lambda\delta f(w)&=\delta\xi(w)=p(w)-2f(w)+\\\
&\int_{-\infty}^{\infty}f(\omega_{1})d\omega_{1}\int_{-\infty}^{\infty}f(\omega_{2})\delta\big{(}w-(\omega_{1}+\omega_{2})\big{)}d\omega_{2}\,.\end{split}$
(12)

The three terms on the right hand side denote (in order from left to right)
the addition of a new sequence, the removal of the surrounding sequences, and
the addition of the resulting merged sequence. In order to evolve this model
to an equilibrium, it is sufficient to set Equation (12) equal to zero to
indicate that a steady state has been reached. When the equilibrium is
reached, the mapping no longer changes the probability density function and
$\delta\xi(w)=0$. At this point, the equation reduces to

$p(w)=2f(w)-(f\star f)(w)\,,$ (13)

where the $\star$ indicates a convolution, as in Equation (11) of Berger &
Asgari-Targhi (2009). It is straightforward to demonstrate, as in Section 5.2
of Berger & Asgari-Targhi (2009), that this can be solved in Fourier space via
Fourier transforms and inverse Fourier transforms. Assuming that the input of
interchanges is a Poisson process such that the windings inserted in the
sequences between interchanges follow an exponential distribution,

$p(w)=\frac{\lambda}{2}e^{-\lambda|w|}\,,$ (14)

for some constant $\lambda$, then the equilibrium solution exhibits a
probability distribution of windings characterized by the equation

$f(w)=\frac{\lambda}{2}\bigg{[}L_{-1}(\lambda w)-I_{1}(w\lambda)\bigg{]}\,,$
(15)

where $I_{1}$ is the Bessel-$I$ function and $L_{-1}$ is a Struve-$L$
function. We note that this function is slightly different than Equation (17)
in Berger & Asgari-Targhi (2009), although we have verified that our solution
has the same form as Figure 5 in their paper. Importantly, for winding numbers
with magnitude greater than unity, $|w|>1$, the probability distribution has
the form of a power-law where $f(w)\sim w^{-2}$.

###  2.3 Energy Distribution of Flares

Consider the distribution of reconnection events in the previous subsection.
Since we know that the magnetic free energy scales with the square of the
crossing number of the set of flux tubes, this must also correspond to a
distribution or power law of energetic flaring events once the equilibrium
state has been reached, since $E_{\rm Free}\sim C^{2}$ as in Equation (8).
Therefore, the difference in energy before and after a reconnection event,
$\Delta E$, is given by the difference in the square of the initial and final
values of $C$, i.e.,

$\Delta E\sim C_{0}^{2}-C_{F}^{2}\,.$ (16)

Berger & Asgari-Targhi (2009) assumed that the reconnection occurs after the
winding is greater than some critical value, $C_{\rm crit}$. They assumed the
interchange removal merged two sequences with twist numbers $\omega_{1}$ and
$\omega_{2}$, where by convention $|\omega_{2}|>|\omega_{1}|$. Critically,
they assume that the merged sequences have opposite sign, on the basis that
mergers of sequences with the same winding direction will not create a release
of energy in a flaring event. The merging after interchange produces a single
sequence with $\omega_{f}=|\omega_{2}|-|\omega_{1}|$, and energy

$\Delta E\sim C_{\rm crit}^{2}-\bigg{(}C_{\rm
crit}-2|\omega_{1}|\bigg{)}^{2}\,.$ (17)

If $\omega_{1}\ll C_{\rm crit}$, then to first order we have

$E\sim 4C_{\rm crit}|\omega_{1}|\,,$ (18)

and we can deduce the distribution of flare energies, $F(E)$ from the
relationship between energies before and after reconnection and the
distribution of winding numbers $w$:

$F(E)\sim E^{-\alpha}\,,$ (19)

where $\alpha=2\gamma-1$, and $\gamma$ is the exponent of the winding
distribution $f(w)\sim|w|^{-\gamma}$, as demonstrated in Section 5.3 of Berger
& Asgari-Targhi (2009).

###  2.4 Expansion of Braiding Model

Figure 2: A schematic diagram of a coherent reconnection event. In this setup,
the two sequences on the left hand side both have positive winding. After the
reconnection event that removes the interchange, the resulting braid on the
left hand side has a winding of $w=+8$.

The Berger & Asgari-Targhi (2009) model assumes that injection of braids with
positive and negative winding are equally probable, which is a reasonable
approach for a slowly-rotating star like the Sun, where the Rossby number in
the sub-surface convection zone is relatively large (Greer et al., 2016).
However, this assumption must begin to fail in more rapidly rotating stars
with smaller Rossby numbers, and we must therefore generalize the model to
include a bias in the sign of $w$.

To generalize this (but keeping the same distribution on the overall amplitude
of $w$ in the injected sequences) we just adjust the prefactor. We define a
parameter $\eta\in[0,1]$ to be the probability that a sequence has a positive
handed winding number. The probability distribution function for the input of
sequences can be written as

$p(w)=\begin{dcases}\eta\lambda e^{-\lambda w}&{\rm if}\,\,w>0\\\
\left(1-\eta\right)\lambda e^{+\lambda w}\,&{\rm if}\,\,w<0\,.\end{dcases}$
(20)

In Equation (20), the value of $\eta$ determines the percentage of sequences
that are locally injected into the corona with the same sign or handedness in
the twist. If there is a bias towards more sequences with the same sign, then
one would expect that there should be a bias towards more coherent interchange
removals, where the surrounding sequences add constructively.

In Figure 2, we show a schematic in which the two sequences have the same sign
and number of windings. After the interchange is removed, only one braid
remains, but it has twice the number of windings as the original two braids,
so no energy is released from the unwinding of braids. It is natural to expect
that as $\eta$ increases or decreases away from $1/2$ (which represents an
equal injection of positive and negative sequences), the steady state winding
distribution will be skewed towards higher numbers. We can demonstrate this
trend analytically by taking the Fourier transform of Equation (20), which
takes the form

$\tilde{p}(k)=\left(\frac{\lambda^{2}-i\lambda
k(2\eta-1)}{\lambda^{2}+k^{2}}\right)\,.$ (21)

Using this result in Equation (13), the solution for the steady state winding
distribution function in Fourier space is

$\tilde{f}(k)=1-\sqrt{\frac{k^{2}+i\lambda k(2\eta-1)}{\lambda^{2}+k^{2}}}\,,$
(22)

and the corresponding real-space distribution is

$\begin{split}f(w)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\left[1-\sqrt{\frac{k^{2}+i\lambda
k(2\eta-1)}{\lambda^{2}+k^{2}}}\right]e^{ikw}dk\,.\end{split}$ (23)

We verified analytically that Equation (23) is symmetric for positive and
negative $w$, such that $f_{\eta}(w)=f_{1-\eta}(-w)$, and that it is
independant of rescaling of $\lambda$, such that for $x=\lambda w$,
$f(x)=f(w)/\lambda$. In order to compute this numerically, we decompose
Equation (22) into real and imaginary components,

$\mathfrak{Re}[\tilde{f}(k)]=1-\sqrt{\frac{\left(\phi+k^{2}\right)}{2\left(k^{2}+\lambda^{2}\right)}}$
(24)

and

$\begin{split}\mathfrak{Im}[\tilde{f}(k)]=-\frac{\lambda
k(2\eta-1)}{\sqrt{2\left(k^{2}+\lambda^{2}\right)\left(\phi+k^{2}\right)}}\,,\end{split}$
(25)

where

$\phi=\bigg{(}k^{4}+k^{2}\lambda^{2}(2\eta-1)^{2}\bigg{)}^{1/2}.$ (26)

It is important to note that $\phi$ in Equation (26) is an even function of
$k$. Therefore $\mathfrak{Re}[\tilde{f}(k)]$ (Equation [24]) is an even
function of $k$ and $\mathfrak{Im}[\tilde{f}(k)]$ (Equation [25]) is an odd
function of $k$. Equation (23) then becomes,

$\begin{split}f(w)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\,\bigg{(}\,\mathfrak{Re}[\tilde{f}(k)]\,\cos(kw)\,\bigg{)}\,dk\\\
-\frac{1}{2\pi}\int_{-\infty}^{\infty}\,\bigg{(}\,\mathfrak{Im}[\tilde{f}(k)]\,\sin(kw)\,\bigg{)}\,dk\,,\end{split}$
(27)

since the imaginary component is an odd function of $k$ and integrates to 0.

Figure 3: The distribution function for winding of sequences in an equilibrium
state. We calculate these distributions by numerically performing the inverse
Fourier transform defined by Equation 23 for $\eta=1/2,3/4$ and 1. The purple
line shows the solution when $\eta=1/2$ for an equal probability of positive
and negative twists, which has a power-law slope of $\gamma=2$, consistent
with the findings of Berger & Asgari-Targhi (2009). As $\eta$ increases, and
more coherent and same signed braids are injected into the corona, the slopes
of the distribution functions get shallower and approach $\gamma=3/2$. This is
plotted using $\lambda=1$, where $|w|=1/\lambda$ is the mean winding magnitude
per injected sequence. Changing $\lambda$ simply re-scales $w$, leaving the
power-law dependencies unchanged.  Figure 4: The distribution function for
winding of sequences in equilibrium. This is calculated in the same way as in
Figure 3, and we show the distribution for positive and negative winding
numbers with $\eta=0,1/4,4/10,1/2,6/10,3/4$ and 1. The distributions for
$\eta$ and $1-\eta$ are symmetric about $w=0$. The equilibrium state reflects
the same sign bias as the input of sequences in Equation 20.

We evaluate Equation (27) numerically and show the resulting distributions in
Figure 3 for $\eta=1/2,3/4$ and 1. For the case of $\eta=1/2$, the problem
reduces to the analytic solution presented in the previous section and in
Berger & Asgari-Targhi (2009). As $\eta$ approaches 1 (or 0), and more
coherent braids are injected into the corona, the FFD becomes shallower. This
result is consistent with the findings of Berger & Asgari-Targhi (2009), who
performed Monte Carlo realizations of the problem with an asymmetric injection
of positive to negative handed sequences.

The distributions of $f(w)$ in Figure 3 are symmetric about the $y$-axis for
positive and negative winding numbers when $\eta=1/2$. In Figure 4, we show
the equilibrium distributions for positive and negative winding numbers for a
range of $\eta\in[0,1]$. When $\eta\neq 1/2$, the resulting equilibrium
distribution is asymmetric about $w=0$ and contains more sequences with the
same sign that is preferentially input. As expected, the equilibrium
distributions for $\eta$ and $1-\eta$ (shown with solid and dashed lines) are
symmetric about $w=0$ $\forall\eta\in[0,1]$. We note that the end-member case,
where $\eta=1$ or $0$, is formally not a physically plausible scenario. This
is because, as shown in Figure 4, the distribution functions for these cases
only have sequences with a single sign, which means the star would not be able
to produce any energetic flares.

###  2.5 Generalized Energy Distribution

It is important to note that the transformation from the distribution of
winding to energy given by Equation 19 is no longer valid for asymmetric
distributions where $\eta\neq 1/2$. Following the notation in Subsection 2.3,
assuming that $|w_{2}|>|w_{1}|$ and allowing for two cases where $w_{2}<0$,
$w_{1}>0$ and $w_{2}>0$, $w_{1}<0$, the energy distribution, denoted by
capital $F(E)$, presented in Equation 20 in Berger & Asgari-Targhi (2009) can
be generalized,

$\begin{split}F(E)\sim\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\\\
\int_{0}^{\infty}\int_{-\infty}^{-w_{1}}\,\bigg{(}\,f(w_{1})f(w_{2})\,\delta\big{(}E-4|w_{1}|C_{\rm
crit}\big{)}\,\bigg{)}\,dw_{2}dw_{1}\\\
+\int_{-\infty}^{0}\int_{|w_{1}|}^{\infty}\,\bigg{(}\,f(w_{1})f(w_{2})\,\delta\big{(}E-4|w_{1}|C_{\rm
crit}\big{)}\,\bigg{)}\,dw_{2}dw_{1}\,.\end{split}$ (28)

In this equation, the first term on the right hand side corresponds to the
reconnection events where $w_{1}>0$. Similarly, the second term on the right
corresponds to reconnection events where $w_{1}<0$. Note that the $\delta$
function in both integrands has an absolute value sign because of the form of
Equation 18, where $E\sim|w_{1}|$. From Figures 3 and 4, the distribution for
$w$ can be approximated as a power law times a negative exponential. The
approximate distribution function is given by

$f(w)\sim\begin{dcases}e^{-(2\eta-1)w/w_{+}}|w|^{-\beta}&{\rm if}\,\,w>0\\\
e^{(2\eta-1)w/w_{-}}|w|^{-\psi}\,&{\rm if}\,\,w<0\,,\end{dcases}$ (29)

for some truncation winding numbers and power law exponents for the positive
and negative cases, $w_{+}$ and $w_{-}$, and power law exponents $\beta$ and
$\psi$. All four of these parameters depend on the value of $\eta$, and
generally $w_{+}>>w_{-}$ and $\beta>\psi$ for $\eta\in[1/2,1]$. Equation 29
was constructed to be symmetric for $\eta$ and $1-\eta$. For the remainder of
this subsection, we solve for the energy distribution for the case of
$\eta\in[1/2,1]$, since the resulting distribution is the same for $\eta$ and
$1-\eta$. By substituting Equation 29 into Equation 28, the distribution is
given by,

$\begin{split}F(E)\sim\int_{0}^{\infty}\int_{-\infty}^{-w_{1}}\,\bigg{(}\,\delta\big{(}E-4|w_{1}|C_{\rm
crit}\big{)}\\\
\,f(w_{1})(e^{(2\eta-1)w_{2}/w_{-}}|w_{2}|^{-\psi})\,\bigg{)}\,dw_{2}dw_{1}+\\\
\int_{-\infty}^{0}\int_{|w_{1}|}^{\infty}\,\bigg{(}\,\delta\big{(}E-4|w_{1}|C_{\rm
crit}\big{)}\,\\\
f(w_{1})(e^{-(2\eta-1)w_{2}/w_{+}}|w_{2}|^{-\beta})\,\bigg{)}\,dw_{2}dw_{1}\,.\end{split}$
(30)

This can be solved analytically, by recalling the definition of the incomplete
gamma function,

$\Gamma(a,x)=\int_{x}^{\infty}t^{a-1}e^{-t}dt\,.$ (31)

We perform two change of variables to integrate Equation (30). For the first
term on the right hand side, we define the variable $v=-(2\eta-1)w_{2}/w_{-}$,
such that $dw_{2}=-w_{-}/(2\eta-1)dv$. Similarly, for the second term on the
right hand side, we define the variable $u=(2\eta-1)w_{2}/w_{+}$, such that
$dw_{2}=w_{+}/(2\eta-1)du$. After evaluating the first integral of the double
integral in each term, Equation (30) reduces to

Figure 5: The energy distribution of flaring events for a range of $\eta$.
These distributions were computed by performing numerical integrations of
Equation 28 for the distributions presented in Figure 4. The energy
distribution becomes shallower as $\eta$ increases.

$\begin{split}F(E)\sim\int_{0}^{\infty}\bigg{[}\bigg{(}\frac{w_{-}}{2\eta-1}\bigg{)}^{1-\psi}f(w_{1})\\\
\Gamma\bigg{(}1-\psi,\frac{2\eta-1}{w_{-}}w_{1}\bigg{)}\delta\big{(}E-4|w_{1}|C_{\rm
crit}\big{)}\bigg{]}dw_{1}+\\\
\int_{-\infty}^{0}\bigg{[}\bigg{(}\frac{w_{+}}{2\eta-1}\bigg{)}^{1-\beta}f(w_{1})\\\
\Gamma\bigg{(}1-\beta,\frac{2\eta-1}{w_{+}}|w_{1}|\bigg{)}\delta\big{(}E-4|w_{1}|C_{\rm
crit}\big{)}\bigg{]}dw_{1}\end{split}$ (32)

Using the composition properties of the Dirac delta function, the coefficient
$\delta(E-4|w_{1}|C_{\rm crit})$ can be written as,

$\begin{split}\delta(E-4|w_{1}|C_{\rm crit})=\\\ \frac{1}{4C_{\rm
crit}}\bigg{(}\delta(w_{1}-\frac{E}{4C_{\rm
crit}})+\delta(w_{1}+\frac{E}{4C_{\rm crit}})\bigg{)}\,.\end{split}$ (33)

By substituting Equation (33) into Equation (32) and performing the integral,
Equation (21) in Berger & Asgari-Targhi (2009) generalizes to the following
energy distribution function,

$\begin{split}F(E)\sim\left(\frac{1}{4C_{\rm
crit}}\right)\bigg{(}\frac{w_{-}}{2\eta-1}\bigg{)}^{1-\psi}\\\
\Gamma\bigg{(}1-\psi,\frac{2\eta-1}{w_{-}}\frac{E}{4C_{\rm
crit}}\bigg{)}f\,\bigg{(}\,\frac{E}{4C_{\rm crit}}\,\bigg{)}\\\
+\left(\frac{1}{4C_{\rm
crit}}\right)\bigg{(}\frac{w_{+}}{2\eta-1}\bigg{)}^{1-\beta}\\\
\Gamma\bigg{(}1-\beta,\frac{2\eta-1}{w_{+}}\frac{E}{4C_{\rm
crit}}\bigg{)}f\,\bigg{(}\,\frac{-E}{4C_{\rm crit}}\,\bigg{)}\,.\end{split}$
(34)

By substituting Equation (29) into Equation (34), the final energy
distribution is given by,

$\begin{split}F(E)\sim\left(\frac{1}{4C_{\rm
crit}}\right)\bigg{(}\frac{w_{-}}{2\eta-1}\bigg{)}^{1-\psi}\Gamma\bigg{(}1-\psi,\frac{2\eta-1}{w_{-}}\frac{E}{4C_{\rm
crit}}\bigg{)}\\\ \,\bigg{(}\,e^{-(2\eta-1)E/(w_{+}4C_{\rm
crit})}|\frac{E}{4C_{\rm crit}}|^{-\beta}\,\bigg{)}\\\ +\left(\frac{1}{4C_{\rm
crit}}\right)\bigg{(}\frac{w_{+}}{2\eta-1}\bigg{)}^{1-\beta}\Gamma\bigg{(}1-\beta,\frac{2\eta-1}{w_{+}}\frac{E}{4C_{\rm
crit}}\bigg{)}\\\ \,\bigg{(}\,e^{-(2\eta-1)E/(w_{-}4C_{\rm
crit})}|\frac{E}{4C_{\rm crit}}|^{-\psi}\,\bigg{)}\\\ \,.\end{split}$ (35)

For large $x$, the incomplete Gamma function can be approximated as the
asymptotic series (Temme, 1975),

$\Gamma(a,x)\approx
x^{a-1}e^{-x}\bigg{(}1+(a-1)x^{-1}+(a-1)(a-2)x^{-2}+...\bigg{)}\,,$ (36)

for $x>>|a-n|$, where $n$ is the index of the term number in the series, or
for a series truncated at the point where $(a-n)/x\sim 1$. Therefore, Equation
(35) has the form of a sum of power laws times a negative exponential for
large flare energies. We verified numerically that Equation (35) has the form
of a power law with an exponent $\alpha\approx\beta+\psi-1$ and an exponential
cutoff set by $w_{-}$ (where $w_{-}<w_{+}$), analagous to the simpler case
found in Berger & Asgari-Targhi (2009), where $\beta=\psi$ and
$\alpha=2\beta-1$. For example, for $w_{-}=10$, the inertial range extends
from $0<w<100$, while for $w_{-}=1$, it extends from $0<w<10$. As can be seen
in Figure 4, although the negative exponential dominates for most of the
domain for negative winding numbers, $\psi$ slightly decreases as $\eta$
increases. Therefore, as $\eta$ changes, $\alpha$ scales primarily with
$\beta$.

To validate this analytic calculation, we computed the energy distribution
$F(\tilde{E})$ for a proxy for the flare energy, $\tilde{E}\sim w/4C_{\rm
crit}$. We computed this by performing a numerical integration of Equation 28
directly, for the distribution functions presented in Figure 4. We show the
energy distributions in Figure 5, for a range of $\eta$. As expected, the
energy distribution becomes shallower as winding parity violation is
increased.

Our main conclusion is that the slope decreases as $\eta\rightarrow 1$. This
result can be interpreted physically as follows. As more coherent sequences
are injected into the corona, the removal of interchanges yields a larger
number of coherent additions of sequences, as shown in Figure 2. Compared to
the symmetric case ($\eta=1/2$), this has the effect of allowing larger values
of $|w|$ to build up, since the injection of a new braid is now more likely to
increase $|w|$ than to decrease it. The winding distribution $f(w)$ thus
becomes flatter. Since the flare energy index $\alpha$ is directly related to
the winding number index, a flatter distribution of winding numbers
immediately produces a flatter energy distribution. This change in slope is
compensated for by a reduction in the overall number of flaring events at
$\eta\to 1$, so that the total amount of energy dissipated remains the same.

Figure 6: The Gaia DR2 $B_{p}-R_{p}$ color of every star in our sample vs. the
measured rotation periods from the literature (Feinstein et al., 2020b; Howard
et al., 2020). Each point is colored by the flare rate, taken from Günther et
al. (in prep). All of these stars were observed by TESS at 2-minute cadence.
Our sample is biased towards stars with rotation periods of $<13$ days due to
the observing strategy of TESS.

##  3 Connection with TESS Observations

Studying stellar activity on a statistical level requires the observation and
identification of a large number of flaring events. Historically, such a
catalogue has been difficult to compile due to the long observational
baselines with high temporal cadences required for each star to capture short-
lived events. It is also a non-trivial task to identify and characterize
flares. Recently, surveys designed to discover extrasolar planets such as
Kepler/K2 (Borucki et al., 2010; Howell et al., 2014) and the Transiting
Exoplanet Survey Satellite (TESS; Ricker et al., 2015) have provided a wealth
of observations that can be used for compiling a catalogue of flares (see also
Section 1).

###  3.1 Sample & Flare Identification

Our calculations suggest that stars with strong Coriolis forces should exhibit
a higher proportion of high-energy flares relative to low-energy ones. To
verify this empirically, we consider the subset of stars that were observed by
TESS at 2-minute cadence, are included in our sample from Feinstein et al.
(2021), and which have rotation periods available in the literature. The
rotation periods used here were measured from TESS and Evryscope light curves,
and were taken directly from Feinstein et al. (2020b) and Howard et al.
(2020), which yields a sample of 1380 stars ranging from
$T_{\textrm{eff}}=2300-9300$ K. A summary of this sub-sample is presented in
Figure 6. Although requiring a measured period significantly decreases the
number of stars in our sample, we still have enough stars to check for
signatures of the Coriolis force.

The flares were detected using the convolutional neural network models (CNNs)
presented in Feinstein et al. (2020b). These models were specifically designed
to find flares in TESS 2-minute cadence light curves and assign a
“probability” to each flare of being real (1) or not (0). We run the 10 CNNs
recommended by Feinstein et al. (2020b) and average the prediction outputs for
each light curve. The average prediction is used as our probability that each
flare is real.

Feinstein et al. (2020b) included all flares with a threshold of $\geq 0.5$ to
be potential true flares, meaning the flares in their sample have a 50%
probability of being a true flare. However, we adopt a more conservative
threshold of 0.9 to ensure that at least 90% of the samples are classified as
true flares.

Given the limited size of the training set for the Feinstein et al. (2020b)
CNNs, we applied additional false-positive filters from Günther et al. (in
prep), which are summarized in Feinstein et al. (2021). The filters remove
flares with peak intensities that are within $3\times$root-mean-square of the
light curve, flares with equivalent durations of $\leq 4$ minutes, and flares
that are associated with periodic events (e.g., eclipsing binaries or variable
stars). Once these filters are applied, we are left with a sample of 16,184
flares with true-flare probabilities $\geq 0.9$ from 869 stars.

###  3.2 Observational Results

Figure 7: Flare frequency distributions (FFDs) for the flares in our sample,
binned by the intensity of the flare in percentage of the star’s normalized
flux. The flare amplitude, $A$, is divided into bins of equal width in log-
space (with 5 bins per dex) ranging from $-2.0\leq log_{10}(A)\leq 1$, and the
vertical axis plots the number of flares observed per star per day in each
amplitude bin. The top panel shows the distribution of flares on stars with
rotation periods, $P_{\rm rot}$, of $<3$ days (522 stars with 11,614 flares
total), and the bottom shows the distribution of flares on stars with $P_{\rm
rot}\geq 3$ days (347 stars with 4,570 flares total). While our sample of
slower rotators is incomplete, there is evidence that the more rapidly
rotating objects have shallower slopes characterized by more energetic events,
which implies stronger winding parity violation.  Figure 8: Five of the
highest amplitude flares (highlighted in blue) from five different stars in
our sample. The inset panels are of the full light curve, displaying various
forms of starspot driven modulation, used to measure the rotation periods of
these stars. Each of these stars are fast rotators, and contribute to the top
panel of Figure 7.  Figure 9: The distributions of the best fitting power-law
exponent and amplitude cutoff parameter for the slow and fast rotators, from
the MCMC fitting process described in §3.2. This shows the full 2D posterior
PDF, including both slopes and cutoffs, with histograms showing the marginal
1D PDFs. In the central panel showing the 2D PDF, the shaded regions are the
68%, 90%, and 95% CIs, which visually demonstrate that the CIs are overlapping
at 68% confidence. This suggests that the distributions of flares are
marginally different for fast and slow rotators.

In Figure 7, we show the intensity distributions of flaring events for stars
that we identify as fast ($P_{\rm rot}<3$ days) and slow rotators ($P_{\rm
rot}\geq 3$ days). In these plots, the flare amplitude, $A$ (expressed as a
percentage of the star’s normalized flux), is divided into bins of equal width
in log-space from $-2.0\leq log_{10}(A)\leq 1$ (with 5 bins per dex, and 30
bins in total), and the vertical axis shows the number of flares, weighted by
their true-flare probability, observed per star per day within each bin. The
bins represent the flare rate for all flares with a threshold of $\geq 0.9$.
The error bars on the flare rate in each bin are defined as follows. Our upper
estimate for the flare rate includes all flares with a threshold $\geq 0.5$.
Our lower estimate is for all flares with a threshold $\geq 0.99$. We have
visually verified that even the very large amplitude flares we measure are
real, and not the result of a failure in the detection algorithm. In Figure 8,
we show the light curves of five of the highest amplitude flares in the fast
rotators to reinforce this conclusion.

In order to evaluate the slope of the distribution of flare amplitudes, and
test whether there is a significant difference between these two samples, we
fit a truncated power-law distribution $dp/dA\propto
A^{-\alpha^{\prime}}e^{-A/A_{*}}$ for amplitudes $A>A_{\rm min}$, where
$A_{\rm min}=10^{-2.0}$ is the smallest amplitude flare for which our data are
complete – there is a visible turn-down in the flare frequency distribution
below this limit. We perform this fit separately to each of the two sub-
samples (fast- and slow-rotators). Note that the form of this equation is
similar to Equation 1, although we include the possibility of a truncation at
high amplitude, and denote the slope here $\alpha^{\prime}$ rather than
$\alpha$ to indicate that we are fitting the index for flare amplitude, not
energy. We fit the slopes using a Markov chain Monte Carlo (MCMC) method
implemented with the emcee package (Goodman & Weare, 2010; Foreman-Mackey et
al., 2013), using the log-likelihood function

$\log\mathcal{L}=\sum_{i}\log\left(\mathcal{N}A_{i}^{-\alpha^{\prime}}e^{-A/A_{*}}\right),$
(37)

where $A_{i}$ is the amplitude of the $i$ th flare in the sample, $A_{*}$ is a
flare amplitude cutoff parameter to be fit, and $\mathcal{N}$ is a
normalization factor chosen to ensure that $\int_{A_{\rm
min}}^{\infty}(dp/dA)\,dA=1$. The MCMC fit has two free parameters – $A_{*}$
and $\alpha^{\prime}$ – and we adopt priors that are flat in $\log A_{*}$ and
$\alpha^{\prime}$.

We initialized our MCMC fit with 100 walkers and iterated for 2500 steps; we
discard the first 1000 steps for burn-in, and verify visually that the chains
are well-converged. We show the median fits obtained by this method, as well
as the confidence intervals around them, overlaid on the data in Figures 7; we
show the joint and marginal posterior probability distribution functions
(PDFs) for the slope $\alpha^{\prime}$ and cutoff parameter $A_{*}$ in Figure
9. As the plots show, there is a marginally significant difference in the
posterior PDFs for the short- and long-period samples — the $2\sigma$ regions
do not overlap, but are close, and the marginal parameter values do overlap at
the $2\sigma$ level. Repeating the fits using different breaks between fast
and slow rotators of 1 day or two days yields qualitatively similar results.
Formally, we find $\alpha^{\prime}=1.866\pm 0.035$ and $\log
A_{*}=0.004_{-0.106}^{+0.131}$ for the short-period sample, and
$\alpha^{\prime}=1.967_{-0.070}^{+0.068}$ and $\log
A_{*}=-0.017_{-0.223}^{+0.336}$ for the long-period sample; the quoted
uncertainties here indicate the 5th to 95th percentile range. To the extent
that the differences are real, and not simply statistical noise, we find the
fast rotators showing a flatter slope but with a cutoff at somewhat smaller
amplitude, and the slow rotators showing a steeper slope but more gradual
cutoff. A shallower slope for fast rotators would be consistent with that seen
in the literature (Hartmann & Noyes, 1987; Maggio et al., 1987; Doyle et al.,
2020), but we emphasise that our present sample does not provide a
statistically significant detection.

Figure 10: Same as Figure 7, with the two sub-samples grouped by Rossby
number. The top panel shows the distribution of flares on stars with Rossby
number, Ro, of $<0.13$ (458 stars with 11,1148 flares total), and the bottom
panel shows the distribution of flares on stars with $\textrm{Ro}\geq 0.13$
(349 stars with 3,507 flares total).  Figure 11: Same as Figure 9, with the
two sub-samples grouped by Rossby number. For $\textrm{Ro}<0.13$, the fits
produced a slope $\alpha^{\prime}=1.847_{-0.033}^{+0.032}$ and log($A_{*}$) =
$-0.012_{-0.102}^{+0.117}$. For $\textrm{Ro}>0.13$ the fits produced
$\alpha^{\prime}=2.155_{-0.115}^{+0.112}$ and log($A_{*}$) =
$-0.111_{-0.348}^{+1.996}$.

It has been established that stellar properties other than rotation period,
such as mass and surface temperature, affect the dynamo. The Rossby number,
Ro, incorporates the rotational period, mass and temperature, and is defined
as $\textrm{Ro}=P/\tau$, where $\tau$ is the convective turnover time (Noyes
et al., 1984). For a given stellar type and $\tau$, the Rossby number is a
direct proxy for rotational period. X-ray luminosity and overall activity
increases with rotational period for a range of stellar types, including Sun-
like stars (Pallavicini et al., 1981; Wright et al., 2011; Candelaresi et al.,
2014) and low mass, fully convective stars (Wright & Drake, 2016). However,
this relationship does not apply for very fast rotators, and the X-ray
luminosity saturates for Ro$<0.13$, independant of spectral type (Vilhu, 1984;
Wright & Drake, 2016).

Motivated by this work, we repeat the same analysis but group both sub-samples
based on Rossby number instead of just the rotational period. We convert
stellar mass to the convective turnover time, $\tau$, using Equation 11 in
Wright et al. (2011),

${\rm log}_{10}(\tau)=1.16–1.49\,{\rm
log}_{10}\bigg{(}\frac{M}{M_{\odot}}\bigg{)}-0.54\,{\rm
log}_{10}^{2}\bigg{(}\frac{M}{M_{\odot}}\bigg{)}\,.$ (38)

By binning the stars based on the saturation break of Rossby number, the two
resulting sub-samples include 458 stars with 11,1148 flares total where
$\textrm{Ro}<0.13$, and 349 stars with 3,507 flares total with
$\textrm{Ro}\geq 0.13$. The resulting distributions are shown in Figure 10,
with the best fit distributions overplotted. The fits were calculated in the
same manner as the fits for sub-samples grouped by rotational period. In
Figure 11, we show the joint and marginal posterior probability distributions
functions for slope $\alpha^{\prime}$ and cutoff parameter $A_{*}$, with the
two sub-samples binned by Rossby number. We find that for $\textrm{Ro}<0.13$,
$\alpha^{\prime}=1.847_{-0.033}^{+0.032}$ and log($A_{*}$) =
$-0.012_{-0.102}^{+0.117}$. For $\textrm{Ro}>0.13$,
$\alpha^{\prime}=2.155_{-0.115}^{+0.112}$ and log($A_{*}$) =
$-0.111_{-0.348}^{+1.996}$. The results of all fits are presented in Table 1.
It appears that the differences in the distributions are more pronounced when
grouped by Rossby number. Since the Rossby number scales with the rotational
period, it is plausible that the coriolis effect investigated in this paper is
responsible for the differences in the slopes of these FFDs.

There are four potential confounding factors that might influence our results.
First and most obviously, the sample is relatively small: only $\approx 3000$
flares in total with amplitude $A>10^{-0.5}$ where we are reasonably complete.
It may simply be that the sample is too small to yield statistically-
significant results. Second, the measured rotation periods are biased towards
$<12.5$ days due to the observing strategy of TESS. Thus our sample of slow
rotators is perhaps better described as a sample of moderator rotators; truly
slow rotators like the Sun are largely excluded. Third, our sample has not
been corrected for detection biases. Methods such as injection-recovery tests
have been used to adjust observed flare rates in previous studies, as in
Figure 11 in Kővári et al. (2020). We do not perform injection-recovery tests
on this sample for two reasons. For one, Feinstein et al. (2020b) demonstrated
that injection recovery tests did not produce accurate results for their CNN,
due to the differences in the shape of the flares themselves. Additionally, we
only start at 0.32% amplitudes here and therefore ignore the smallest flares
(which are most affected by this bias). A final confounding factor is that,
while rotation period is a reasonable proxy for the importance of Coriolis
forces in the dynamo, it is certainly not the only factor that might affect
flare distributions. In particular, stars of different ages and effective
temperature almost certainly have different surface magnetic field strengths,
and our fast and slow rotator samples mix together a wide range of stellar age
and effective temperature. Similar to our marginal agreement here, Feinstein
et al. (2020b) measured the FFD slope for stars with age $\leq 50$ Myr and
$>50$ Myr and found they were in a 1-$\sigma$ agreement with each other. For
the specific details of this analysis, we refer the reader to Figures 13-14
and Section 4.1 in Feinstein et al. (2020b). While the results from our
analysis incorporating the Rossby number are highly suggestive, it is possible
that an analysis that separates all of these contributing factors – age,
rotation period, and effective temperature – might yield stronger results.

Table 1: Summary of the parameter fits for the flare frequency distributions binned by rotation period and Rossby number.  Sample | slope, $\alpha^{\prime}$ | log($A_{*}$)
---|---|---
$P_{\textrm{rot}}<3$ days | $1.866\pm 0.035$ | $0.004_{-0.106}^{+0.131}$
$P_{\textrm{rot}}\geq 3$ days | $1.967_{-0.070}^{+0.068}$ | $-0.017_{-0.223}^{+0.336}$
$\textrm{Ro}<0.13$ | $1.847_{-0.033}^{+0.032}$ | $-0.012_{-0.102}^{+0.117}$
$\textrm{Ro}\geq 0.13$ | $2.155_{-0.115}^{+0.112}$ | $-0.111_{-0.348}^{+1.996}$

##  4 Conclusions

In this paper, we have investigated the role of the Coriolis force in stellar
dynamos. More specifically, we considered the direct current model of coronal
heating (Parker, 1972) along with the model of magnetic braiding and
reconnection (proposed by Berger & Asgari-Targhi 2009; see also Berger et al.
2015). We expanded the braiding model to incorporate the effects of the
Coriolis force, which should be more dominant in faster rotating stars.
Specifically, we incorporated a bias in handedness of injected braids, which
should be present in stars with dynamos that are strongly affected by the
Coriolis force. As increasingly coherent braids are injected into the corona,
the slopes of the resulting power-law distributions of energetic flaring
events decreases in magnitude (corresponding to a weaker decline in occurrence
rates toward increasing flare energies). We search for this effect in the
flare frequency distributions for stars observed by TESS that have measured
rotation periods; however, while the results are suggestively consistent with
the theoretical prediction, the sample is too small to yield a definitive
statistical conclusion.

While the results presented in this study are suggestive, there are many
opportunities for future theoretical work. It would be informative to perform
numerical experiments of reconnection events for both the DC and AC regimes of
coronal heating. Specifically, numerical simulations of avalanching
reconnection events similar to those presented by Lu & Hamilton (1991) may
also exhibit shallower slopes with driving that mimics effects of the Coriolis
force. Being able to directly probe the relation between the magnitude of the
Coriolis force and the $\eta$ (handedness) parameter would be an invaluable
measurement to test the role of rotation on the topology of the flux tubes and
hence the energy of the flaring events. This type of calculation would also
yield insights into the connection between the theory of self-organized
criticality and the DC braiding and reconnection picture of coronal heating.

Alternatively, it is possible that the power-law distributions of stellar
flares are simply realizations of forced MHD turbulence, along the framework
of the AC regime. Therefore, numerical simulations of forced coronal plasmas
that track reconnection events may reveal an alternative explanation for the
differences in the observed slopes.

Several possible adjustments to the observational data could further validate
the hypothesis presented here. We have included a relatively small ($\approx
3,000$ stars) sample of stars with measured rotation periods. However, there
are $10^{5}$ stars with high-cadence observations from TESS. By measuring the
rotation periods for more stars with known flaring events, we would be able to
vastly increase the statistical sample in this study. This analysis would also
be improved by the proposed increase from a 10- to 3-minute cadence of TESS
Full-Frame Images for the 2nd extended mission. This improvement would
increase our sample of stars by an order of magnitude, which would improve the
slope measurements.

It will also be necessary to develop a reliable method for measuring long
rotation periods ($P>13$ days) in TESS (e.g. Lu et al., 2020; Breton et al.,
2021; Claytor et al., 2021). In the current state, we are limited both in the
baseline for the TESS sector observing strategy ($\sim 27$ days) and the
orbital gap halfway through each sector. Traditional methods of measuring
rotation periods often identify the systematics that are associated with the
beginning or end of each orbit. We are thus limited to only reliably being
able to measure short rotation periods. Alternatively, indirect estimates of
stellar rotation periods may be useful for expanding our sample (Mamajek &
Hillenbrand, 2008).

It would also be useful to improve the convolutional neural networks presented
in Feinstein et al. (2020b) to detect flares on rapid rotating stars ($P<1$
day). When these methods were developed, there were relatively few examples of
these light curves, so that the original models were trained on a limited
number of high energy events on rapidly rotating stars. As a result, the sharp
rotational features are often confused for flare events. Now that more high
energy events have been observed, we have the opportunity to improve our
method for flare identification using machine learning.

Finally, it would be of significant interest to test the impact of stellar age
versus rotation period on the flare rates. We were not able to do this
analysis with our relatively small sample. However, populations within young
stellar clusters, such as Pleiades, Hyades, or Praesepe, have well constrained
masses and ages. The activity in young clusters has already been examined
using K2 long-cadence data (Ilin et al., 2019, 2021), although this cadence is
inadequate for detecting low-amplitude flares. Examining activity in young
stellar clusters using new high-cadence TESS data may yield key insights into
the relationship between rotation, dynamo, and spin-down.

##  5 Acknowledgements

We thank Mitch Berger, Mahboubeh Asgari-Targhi, Eugene Chiang, Amir Siraj,
Jared Siegel, Dimitar Sasselov, Konstantin Batygin, Jacob Bean, Juliette
Becker, Fred Ciesla, and Harriet Lau for useful conversations. We thank the
two anonymous reviewers for insightful comments and constructive suggestions
that strengthened the scientific content of this manuscript, especially for
suggesting searching for differences in the power law distribution based on
the Rossby number.

ADF acknowledges support from the National Science Foundation Graduate
Research Fellowship Program under Grant No. DGE-1746045. JRB acknowledges
financial support from the Australian National University, via the Deakin PhD
and Dean’s Higher Degree Research (theoretical physics) Scholarships and the
Australian Government via the Australian Government Research Training Program
Fee-Offset Scholarship. MRK acknowledges support from the Australian Research
Council through its Future Fellowships scheme, award FT180100375. LAR
gratefully acknowledges support from the Research Corporation for Science
Advancement through a Cottrell Scholar Award. C. F. acknowledges funding
provided by the Australian Research Council (Future Fellowship FT180100495),
and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). MF is
grateful for the support he received from the Hauck Foundation.

## References

  * Abramenko et al. (1997) Abramenko, V. I., Wang, T., & Yurchishin, V. B. 1997, Sol. Phys., 174, 291, doi: 10.1023/A:1004957515498
  * Andrade et al. (1998) Andrade, R. F. S., Schellnhuber, H. J., & Claussen, M. 1998, Physica A Statistical Mechanics and its Applications, 254, 557, doi: 10.1016/S0378-4371(98)00057-0
  * Antolin et al. (2021) Antolin, P., Pagano, P., Testa, P., Petralia, A., & Reale, F. 2021, Nature Astronomy, 5, 54, doi: 10.1038/s41550-020-1199-8
  * Aschwanden et al. (1998) Aschwanden, M. J., Dennis, B. R., & Benz, A. O. 1998, ApJ, 497, 972, doi: 10.1086/305484
  * Aschwanden & Güdel (2021) Aschwanden, M. J., & Güdel, M. 2021, The Astrophysical Journal, 910, 41, doi: 10.3847/1538-4357/abdec7
  * Asgari-Targhi et al. (2015) Asgari-Targhi, M., Schmelz, J. T., Imada, S., Pathak, S., & Christian, G. M. 2015, ApJ, 807, 146, doi: 10.1088/0004-637X/807/2/146
  * Asgari-Targhi & van Ballegooijen (2012) Asgari-Targhi, M., & van Ballegooijen, A. A. 2012, ApJ, 746, 81, doi: 10.1088/0004-637X/746/1/81
  * Asgari-Targhi et al. (2013) Asgari-Targhi, M., van Ballegooijen, A. A., Cranmer, S. R., & DeLuca, E. E. 2013, ApJ, 773, 111, doi: 10.1088/0004-637X/773/2/111
  * Asgari-Targhi et al. (2014) Asgari-Targhi, M., van Ballegooijen, A. A., & Imada, S. 2014, ApJ, 786, 28, doi: 10.1088/0004-637X/786/1/28
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Audard et al. (2000) Audard, M., Güdel, M., Drake, J. J., & Kashyap, V. L. 2000, ApJ, 541, 396, doi: 10.1086/309426
  * Babcock & Westervelt (1990) Babcock, K. L., & Westervelt, R. M. 1990, Phys. Rev. Lett., 64, 2168, doi: 10.1103/PhysRevLett.64.2168
  * Bak et al. (1989) Bak, P., Chen, K., & Creutz, M. 1989, Nature, 342, 780, doi: 10.1038/342780a0
  * Bak et al. (1990) Bak, P., Chen, K., & Tang, C. 1990, Physics Letters A, 147, 297, doi: 10.1016/0375-9601(90)90451-S
  * Bak et al. (1997) Bak, P., Paczuski, M., & Shubik, M. 1997, Physica A Statistical Mechanics and its Applications, 246, 430, doi: 10.1016/S0378-4371(97)00401-9
  * Bak & Sneppen (1993) Bak, P., & Sneppen, K. 1993, Phys. Rev. Lett., 71, 4083, doi: 10.1103/PhysRevLett.71.4083
  * Bak & Tang (1989) Bak, P., & Tang, C. 1989, J. Geophys. Res., 94, 15,635, doi: 10.1029/JB094iB11p15635
  * Bak et al. (1987) Bak, P., Tang, C., & Wiesenfeld, K. 1987, Phys. Rev. Lett., 59, 381, doi: 10.1103/PhysRevLett.59.381
  * Bak et al. (1988) —. 1988, Phys. Rev. A, 38, 364, doi: 10.1103/PhysRevA.38.364
  * Bao & Zhang (1998) Bao, S., & Zhang, H. 1998, ApJ, 496, L43, doi: 10.1086/311232
  * Berger (1993) Berger, M. A. 1993, Phys. Rev. Lett., 70, 705, doi: 10.1103/PhysRevLett.70.705
  * Berger & Asgari-Targhi (2009) Berger, M. A., & Asgari-Targhi, M. 2009, ApJ, 705, 347, doi: 10.1088/0004-637X/705/1/347
  * Berger et al. (2015) Berger, M. A., Asgari-Targhi, M., & Deluca, E. E. 2015, Journal of Plasma Physics, 81, 395810404, doi: 10.1017/S0022377815000483
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977, doi: 10.1126/science.1185402
  * Breton et al. (2021) Breton, S. N., Santos, A. R. G., Bugnet, L., et al. 2021, A&A, 647, A125, doi: 10.1051/0004-6361/202039947
  * Brown et al. (1998) Brown, J. C., McArthur, G. K., Barrett, R. K., McIntosh, S. W., & Emslie, A. G. 1998, Sol. Phys., 179, 379, doi: 10.1023/A:1005011107402
  * Browning et al. (2016) Browning, M. K., Weber, M. A., Chabrier, G., & Massey, A. P. 2016, ApJ, 818, 189, doi: 10.3847/0004-637X/818/2/189
  * Candelaresi et al. (2014) Candelaresi, S., Hillier, A., Maehara, H., Brandenburg, A., & Shibata, K. 2014, ApJ, 792, 67, doi: 10.1088/0004-637X/792/1/67
  * Cargill (1994) Cargill, P. J. 1994, ApJ, 422, 381, doi: 10.1086/173733
  * Carlson & Langer (1989) Carlson, J. M., & Langer, J. S. 1989, Phys. Rev. Lett., 62, 2632, doi: 10.1103/PhysRevLett.62.2632
  * Carrington (1859) Carrington, R. C. 1859, MNRAS, 20, 13, doi: 10.1093/mnras/20.1.13
  * Charbonneau (2010) Charbonneau, P. 2010, Living Reviews in Solar Physics, 7, 3, doi: 10.12942/lrsp-2010-3
  * Charbonneau et al. (2001) Charbonneau, P., McIntosh, S. W., Liu, H.-L., & Bogdan, T. J. 2001, Sol. Phys., 203, 321, doi: 10.1023/A:1013301521745
  * Choudhuri & D’Silva (1990) Choudhuri, A. R., & D’Silva, S. 1990, A&A, 239, 326
  * Choudhuri & Gilman (1987) Choudhuri, A. R., & Gilman, P. A. 1987, ApJ, 316, 788, doi: 10.1086/165243
  * Claytor et al. (2021) Claytor, Z. R., van Saders, J. L., Llama, J., et al. 2021, arXiv e-prints, arXiv:2104.14566.  https://arxiv.org/abs/2104.14566
  * Crosby et al. (1993) Crosby, N. B., Aschwanden, M. J., & Dennis, B. R. 1993, Sol. Phys., 143, 275, doi: 10.1007/BF00646488
  * Dahlburg et al. (2005) Dahlburg, R. B., Klimchuk, J. A., & Antiochos, S. K. 2005, ApJ, 622, 1191, doi: 10.1086/425645
  * Datlowe et al. (1974) Datlowe, D. W., Elcan, M. J., & Hudson, H. S. 1974, Sol. Phys., 39, 155, doi: 10.1007/BF00154978
  * Davenport (2016) Davenport, J. R. A. 2016, The Astrophysical Journal, 829, 23, doi: 10.3847/0004-637X/829/1/23
  * de Arcangelis et al. (2006) de Arcangelis, L., Godano, C., Lippiello, E., & Nicodemi, M. 2006, Phys. Rev. Lett., 96, 051102, doi: 10.1103/PhysRevLett.96.051102
  * Dendy et al. (1998) Dendy, R. O., Helander, P., & Tagger, M. 1998, A&A, 337, 962.  https://arxiv.org/abs/astro-ph/9907055
  * Dennis (1985) Dennis, B. R. 1985, Sol. Phys., 100, 465, doi: 10.1007/BF00158441
  * Dmitruk & Gómez (1997) Dmitruk, P., & Gómez, D. O. 1997, ApJ, 484, L83, doi: 10.1086/310760
  * Doyle et al. (2020) Doyle, L., Ramsay, G., & Doyle, J. G. 2020, MNRAS, 494, 3596, doi: 10.1093/mnras/staa923
  * Doyle et al. (2019) Doyle, L., Ramsay, G., Doyle, J. G., & Wu, K. 2019, MNRAS, 489, 437, doi: 10.1093/mnras/stz2205
  * Doyle et al. (2018) Doyle, L., Ramsay, G., Doyle, J. G., Wu, K., & Scullion, E. 2018, MNRAS, 480, 2153, doi: 10.1093/mnras/sty1963
  * Drake (1971) Drake, J. F. 1971, Sol. Phys., 16, 152, doi: 10.1007/BF00154510
  * D’Silva & Choudhuri (1993) D’Silva, S., & Choudhuri, A. R. 1993, A&A, 272, 621
  * D’Silva & Howard (1995) D’Silva, S., & Howard, R. F. 1995, Sol. Phys., 159, 63, doi: 10.1007/BF00733032
  * Einaudi & Velli (1999) Einaudi, G., & Velli, M. 1999, Physics of Plasmas, 6, 4146, doi: 10.1063/1.873679
  * Einaudi et al. (1996) Einaudi, G., Velli, M., Politano, H., & Pouquet, A. 1996, ApJ, 457, L113, doi: 10.1086/309893
  * Feinstein et al. (2020a) Feinstein, A., Montet, B., & Ansdell, M. 2020a, The Journal of Open Source Software, 5, 2347, doi: 10.21105/joss.02347
  * Feinstein et al. (2020b) Feinstein, A. D., Montet, B. T., Ansdell, M., et al. 2020b, AJ, 160, 219, doi: 10.3847/1538-3881/abac0a
  * Feinstein et al. (2021) Feinstein, A. D., Seligman, D. Z., Günther, M. N., & Adams, F. C. 2021, arXiv e-prints, arXiv:2109.07011.  https://arxiv.org/abs/2109.07011
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Galsgaard & Nordlund (1996) Galsgaard, K., & Nordlund, Å. 1996, J. Geophys. Res., 101, 13445, doi: 10.1029/96JA00428
  * Galtier (1999) Galtier, S. 1999, ApJ, 521, 483, doi: 10.1086/307537
  * Galtier & Pouquet (1998) Galtier, S., & Pouquet, A. 1998, Sol. Phys., 179, 141, doi: 10.1023/A:1005056102064
  * Georgoulis et al. (1998) Georgoulis, M. K., Velli, M., & Einaudi, G. 1998, ApJ, 497, 957, doi: 10.1086/305486
  * Ginsburg et al. (2019) Ginsburg, A., Sipőcz, B. M., Brasseur, C. E., et al. 2019, AJ, 157, 98, doi: 10.3847/1538-3881/aafc33
  * Golub & Pasachoff (1997) Golub, L., & Pasachoff, J. M. 1997, The Solar Corona
  * Goodman & Weare (2010) Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65, doi: 10.2140/camcos.2010.5.65
  * Greer et al. (2016) Greer, B. J., Hindman, B. W., & Toomre, J. 2016, ApJ, 824, 128, doi: 10.3847/0004-637X/824/2/128
  * Grieger (1992) Grieger, B. 1992, Physica A Statistical Mechanics and its Applications, 191, 51, doi: 10.1016/0378-4371(92)90505-K
  * Günther et al. (2020) Günther, M. N., Zhan, Z., Seager, S., et al. 2020, The Astronomical Journal, 159, 60, doi: 10.3847/1538-3881/ab5d3a
  * Gutenberg & Richter (1956) Gutenberg, B., & Richter, C. F. 1956, Annals of Geophysics, 53, doi: 10.4401/ag-4588
  * Hale et al. (1919) Hale, G. E., Ellerman, F., Nicholson, S. B., & Joy, A. H. 1919, ApJ, 49, 153, doi: 10.1086/142452
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hartmann & Noyes (1987) Hartmann, L. W., & Noyes, R. W. 1987, ARA&A, 25, 271, doi: 10.1146/annurev.aa.25.090187.001415
  * Hesse & Gross (2014) Hesse, J., & Gross, T. 2014, Frontiers in Systems Neuroscience, 8, 166, doi: 10.3389/fnsys.2014.00166
  * Holzwarth (2007) Holzwarth, V. 2007, Mem. Soc. Astron. Italiana, 78, 271.  https://arxiv.org/abs/0709.3008
  * Howard (1993) Howard, R. F. 1993, Sol. Phys., 145, 105, doi: 10.1007/BF00627986
  * Howard et al. (2019) Howard, W. S., Corbett, H., Law, N. M., et al. 2019, The Astrophysical Journal, 881, 9, doi: 10.3847/1538-4357/ab2767
  * Howard et al. (2020) —. 2020, ApJ, 895, 140, doi: 10.3847/1538-4357/ab9081
  * Howell et al. (2014) Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398, doi: 10.1086/676406
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ilin et al. (2019) Ilin, E., Schmidt, S. J., Davenport, J. R. A., & Strassmeier, K. G. 2019, A&A, 622, A133, doi: 10.1051/0004-6361/201834400
  * Ilin et al. (2021) Ilin, E., Schmidt, S. J., Poppenhäger, K., et al. 2021, A&A, 645, A42, doi: 10.1051/0004-6361/202039198
  * Ionson (1982) Ionson, J. A. 1982, ApJ, 254, 318, doi: 10.1086/159736
  * Ionson (1985) —. 1985, Sol. Phys., 100, 289, doi: 10.1007/BF00158433
  * Kadanoff et al. (1989) Kadanoff, L. P., Nagel, S. R., Wu, L., & Zhou, S.-M. 1989, Phys. Rev. A, 39, 6524, doi: 10.1103/PhysRevA.39.6524
  * Kővári et al. (2020) Kővári, Z., Oláh, K., Günther, M. N., et al. 2020, A&A, 641, A83, doi: 10.1051/0004-6361/202038397
  * Kitchatinov & Olemskoy (2015) Kitchatinov, L. L., & Olemskoy, S. V. 2015, Research in Astronomy and Astrophysics, 15, 1801, doi: 10.1088/1674-4527/15/11/003
  * Krucker & Benz (2000) Krucker, S., & Benz, A. O. 2000, Sol. Phys., 191, 341, doi: 10.1023/A:1005255608792
  * Kulsrud (1998) Kulsrud, R. M. 1998, Physics of Plasmas, 5, 1599, doi: 10.1063/1.872827
  * Lee et al. (1993) Lee, T. T., Petrosian, V., & McTiernan, J. M. 1993, ApJ, 412, 401, doi: 10.1086/172929
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration, Cardoso, J. V. d. M., Hedges, C., et al. 2018, Lightkurve: Kepler and TESS time series analysis in Python, Astrophysics Source Code Library.  http://ascl.net/1812.013
  * Lin et al. (2019) Lin, C. L., Ip, W. H., Hou, W. C., Huang, L. C., & Chang, H. Y. 2019, The Astrophysical Journal, 873, 97, doi: 10.3847/1538-4357/ab041c
  * Lin et al. (1984) Lin, R. P., Schwartz, R. A., Kane, S. R., Pelling, R. M., & Hurley, K. C. 1984, ApJ, 283, 421, doi: 10.1086/162321
  * Litvinenko (1996) Litvinenko, Y. E. 1996, Sol. Phys., 167, 321, doi: 10.1007/BF00146342
  * Liu et al. (2014) Liu, Y., Hoeksema, J. T., Bobra, M., et al. 2014, The Astrophysical Journal, 785, 13, doi: 10.1088/0004-637x/785/1/13
  * Longcope & Fisher (1996) Longcope, D. W., & Fisher, G. H. 1996, ApJ, 458, 380, doi: 10.1086/176821
  * Longcope et al. (1998) Longcope, D. W., Fisher, G. H., & Pevtsov, A. A. 1998, ApJ, 507, 417, doi: 10.1086/306312
  * Longcope & Pevtsov (2003) Longcope, D. W., & Pevtsov, A. A. 2003, Advances in Space Research, 32, 1845, doi: 10.1016/S0273-1177(03)90618-1
  * Longcope et al. (2007) Longcope, D. W., Ravindra, B., & Barnes, G. 2007, ApJ, 668, 571, doi: 10.1086/521095
  * Longcope & Sudan (1994) Longcope, D. W., & Sudan, R. N. 1994, ApJ, 437, 491, doi: 10.1086/175013
  * Longcope & Welsch (2000) Longcope, D. W., & Welsch, B. T. 2000, ApJ, 545, 1089, doi: 10.1086/317846
  * Lu & Hamilton (1991) Lu, E. T., & Hamilton, R. J. 1991, ApJ, 380, L89, doi: 10.1086/186180
  * Lu et al. (1993) Lu, E. T., Hamilton, R. J., McTiernan, J. M., & Bromund, K. R. 1993, ApJ, 412, 841, doi: 10.1086/172966
  * Lu et al. (2020) Lu, Y., Angus, R., Agüeros, M. A., et al. 2020, AJ, 160, 168, doi: 10.3847/1538-3881/abada4
  * MacTaggart et al. (2021) MacTaggart, D., Prior, C., Raphaldini, B., Romano, P., & Guglielmino, S. 2021, arXiv e-prints, arXiv:2106.11638.  https://arxiv.org/abs/2106.11638
  * Maggio et al. (1987) Maggio, A., Sciortino, S., Vaiana, G. S., et al. 1987, ApJ, 315, 687, doi: 10.1086/165170
  * Mamajek & Hillenbrand (2008) Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264, doi: 10.1086/591785
  * Manek & Brummell (2021) Manek, B., & Brummell, N. 2021, ApJ, 909, 72, doi: 10.3847/1538-4357/abd859
  * Manek et al. (2018) Manek, B., Brummell, N., & Lee, D. 2018, ApJ, 859, L27, doi: 10.3847/2041-8213/aac723
  * Mikic et al. (1989) Mikic, Z., Schnack, D. D., & van Hoven, G. 1989, ApJ, 338, 1148, doi: 10.1086/167265
  * Mitra-Kraev & Benz (2001) Mitra-Kraev, U., & Benz, A. O. 2001, A&A, 373, 318, doi: 10.1051/0004-6361:20010524
  * Moffatt (1978) Moffatt, H. K. 1978, Magnetic field generation in electrically conducting fluids
  * Moffatt & Tsinober (1990) Moffatt, H. K., & Tsinober, A. 1990, Topological Fluid Mechanics
  * Morales & Santos (2020) Morales, L. F., & Santos, N. A. 2020, Sol. Phys., 295, 155, doi: 10.1007/s11207-020-01713-0
  * Nagel & Herrmann (1993) Nagel, K., & Herrmann, H. J. 1993, Physica A Statistical Mechanics and its Applications, 199, 254, doi: 10.1016/0378-4371(93)90006-P
  * Newman (1996) Newman, M. E. J. 1996, Proceedings of the Royal Society of London Series B, 263, 1605.  https://arxiv.org/abs/adap-org/9607002
  * Newman & Sneppen (1996) Newman, M. E. J., & Sneppen, K. 1996, Phys. Rev. E, 54, 6226, doi: 10.1103/PhysRevE.54.6226
  * Notsu et al. (2013) Notsu, Y., Shibayama, T., Maehara, H., et al. 2013, ApJ, 771, 127, doi: 10.1088/0004-637X/771/2/127
  * Noyes et al. (1984) Noyes, R. W., Hartmann, L. W., Baliunas, S. L., Duncan, D. K., & Vaughan, A. H. 1984, ApJ, 279, 763, doi: 10.1086/161945
  * Olami et al. (1992) Olami, Z., Feder, H. J. S., & Christensen, K. 1992, Phys. Rev. Lett., 68, 1244, doi: 10.1103/PhysRevLett.68.1244
  * Osten & Brown (1999) Osten, R. A., & Brown, A. 1999, ApJ, 515, 746, doi: 10.1086/307034
  * Pallavicini et al. (1981) Pallavicini, R., Golub, L., Rosner, R., et al. 1981, ApJ, 248, 279, doi: 10.1086/159152
  * Parker (1955) Parker, E. N. 1955, ApJ, 122, 293, doi: 10.1086/146087
  * Parker (1972) —. 1972, ApJ, 174, 499, doi: 10.1086/151512
  * Parker (1983) —. 1983, ApJ, 264, 642, doi: 10.1086/160637
  * Parker (1988) —. 1988, ApJ, 330, 474, doi: 10.1086/166485
  * Parker (1989) —. 1989, Sol. Phys., 121, 271, doi: 10.1007/BF00161700
  * Parker (1994) —. 1994, Spontaneous current sheets in magnetic fields : with applications to stellar x-rays. International Series in Astronomy and Astrophysics, 1
  * Parker (2009) —. 2009, Space Sci. Rev., 144, 15, doi: 10.1007/s11214-008-9445-x
  * Pevtsov et al. (1995) Pevtsov, A. A., Canfield, R. C., & Metcalf, T. R. 1995, ApJ, 440, L109, doi: 10.1086/187773
  * Porter et al. (1987) Porter, J. G., Moore, R. L., Reichmann, E. J., Engvold, O., & Harvey, K. L. 1987, ApJ, 323, 380, doi: 10.1086/165835
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Priest & Forbes (2000) Priest, E., & Forbes, T. 2000, Magnetic Reconnection
  * Prior & MacTaggart (2016) Prior, C., & MacTaggart, D. 2016, Geophysical and Astrophysical Fluid Dynamics, 110, 432, doi: 10.1080/03091929.2016.1216552
  * Raetz et al. (2020) Raetz, S., Stelzer, B., Damasso, M., & Scholz, A. 2020, Astronomy and Astrophysics, 637, A22, doi: 10.1051/0004-6361/201937350
  * Rempel (2011) Rempel, M. 2011, Solar Convection Zone Dynamics, ed. M. P. Miralles & J. Sánchez Almeida, Vol. 4, 23
  * Ribeiro et al. (2010) Ribeiro, T. L., Copelli, M., Caixeta, F., et al. 2010, PLoS ONE, 5, e14129, doi: 10.1371/journal.pone.0014129
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003, doi: 10.1117/1.JATIS.1.1.014003
  * Rosner & Vaiana (1978) Rosner, R., & Vaiana, G. S. 1978, ApJ, 222, 1104, doi: 10.1086/156227
  * Schrijver & Zwaan (2000) Schrijver, C. J., & Zwaan, C. 2000, Solar and Stellar Magnetic Activity
  * Seehafer (1990) Seehafer, N. 1990, Sol. Phys., 125, 219, doi: 10.1007/BF00158402
  * Seligman et al. (2014) Seligman, D., Petrie, G. J. D., & Komm, R. 2014, ApJ, 795, 113, doi: 10.1088/0004-637X/795/2/113
  * Shakhovskaia (1989) Shakhovskaia, N. I. 1989, Sol. Phys., 121, 375, doi: 10.1007/BF00161707
  * Shibayama et al. (2013) Shibayama, T., Maehara, H., Notsu, S., et al. 2013, The Astrophysical Journals, 209, 5, doi: 10.1088/0067-0049/209/1/5
  * Soderblom (2010) Soderblom, D. R. 2010, ARA&A, 48, 581, doi: 10.1146/annurev-astro-081309-130806
  * Solanki et al. (2006) Solanki, S. K., Inhester, B., & Schüssler, M. 2006, Reports on Progress in Physics, 69, 563, doi: 10.1088/0034-4885/69/3/R02
  * Sornette & Sornette (1989) Sornette, A., & Sornette, D. 1989, EPL (Europhysics Letters), 9, 197, doi: 10.1209/0295-5075/9/3/002
  * Sturrock et al. (1990) Sturrock, P. A., Dixon, W. W., Klimchuk, J. A., & Antiochos, S. K. 1990, ApJ, 356, L31, doi: 10.1086/185743
  * Sturrock et al. (1984) Sturrock, P. A., Kaufman, P., Moore, R. L., & Smith, D. F. 1984, Sol. Phys., 94, 341, doi: 10.1007/BF00151322
  * Temme (1975) Temme, N. M. 1975, Mathematics of Computation, 29, 1109, doi: 10.2307/2005750
  * Turcotte (1999) Turcotte, D. L. 1999, Reports on Progress in Physics, 62, 1377, doi: 10.1088/0034-4885/62/10/201
  * Turcotte et al. (2002) Turcotte, D. L., Malamud, B. D., Guzzetti, F., & Reichenbach, P. 2002, Proceedings of the National Academy of Sciences, 99, 2530, doi: 10.1073/pnas.012582199
  * van Ballegooijen (1986) van Ballegooijen, A. A. 1986, ApJ, 311, 1001, doi: 10.1086/164837
  * van Ballegooijen et al. (2014) van Ballegooijen, A. A., Asgari-Targhi, M., & Berger, M. A. 2014, ApJ, 787, 87, doi: 10.1088/0004-637X/787/1/87
  * van Ballegooijen et al. (2011) van Ballegooijen, A. A., Asgari-Targhi, M., Cranmer, S. R., & DeLuca, E. E. 2011, ApJ, 736, 3, doi: 10.1088/0004-637X/736/1/3
  * Van Der Walt et al. (2011) Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Vilhu (1984) Vilhu, O. 1984, A&A, 133, 117
  * Weber & Browning (2016) Weber, M. A., & Browning, M. K. 2016, ApJ, 827, 95, doi: 10.3847/0004-637X/827/2/95
  * Weber et al. (2011) Weber, M. A., Fan, Y., & Miesch, M. S. 2011, ApJ, 741, 11, doi: 10.1088/0004-637X/741/1/11
  * Weber et al. (2013) —. 2013, Sol. Phys., 287, 239, doi: 10.1007/s11207-012-0093-7
  * Wheatland (2000) Wheatland, M. S. 2000, ApJ, 532, 1209, doi: 10.1086/308605
  * Withbroe & Noyes (1977) Withbroe, G. L., & Noyes, R. W. 1977, ARA&A, 15, 363, doi: 10.1146/annurev.aa.15.090177.002051
  * Wright & Drake (2016) Wright, N. J., & Drake, J. J. 2016, Nature, 535, 526, doi: 10.1038/nature18638
  * Wright et al. (2011) Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48, doi: 10.1088/0004-637X/743/1/48
  * Yang & Liu (2019) Yang, H., & Liu, J. 2019, The Astrophysical Journals, 241, 29, doi: 10.3847/1538-4365/ab0d28
  * Zirker & Cleveland (1993a) Zirker, J. B., & Cleveland, F. M. 1993a, Sol. Phys., 144, 341, doi: 10.1007/BF00627598
  * Zirker & Cleveland (1993b) —. 1993b, Sol. Phys., 145, 119, doi: 10.1007/BF00627988