# The Sun’s Non-Potential Corona over Solar Cycle 24

Anthony R. Yeates  Department of Mathematical Sciences, Durham University,
Durham, DH1 3LE, UK

###### keywords:

Corona, Models; Helicity, Magnetic; Magnetic fields, Corona; Solar Cycle,
Models

##  1 Introduction

There is little doubt that the Sun’s coronal magnetic field does not simply
evolve through a sequence of current-free equilibria as assumed by the
traditional Potential Field Source Surface, or PFSS model (Altschuler and
Newkirk, 1969; Schatten, Wilcox, and Ness, 1969). Indirect evidence for low
coronal currents includes the presence of sheared filament channels (e.g.,
Martin, 1998; Mackay et al., 2010; Sheeley et al., 2013), the twisted shape of
extreme-ultraviolet or X-ray loops (e.g., Pevtsov, 2002; Malanushenko et al.,
2014), and the presence of coronal flux rope cavities (e.g., Rumińska et al.,
2022).

On theoretical grounds one expects free magnetic energy – above that of a
current-free configuration – to build up only in regions of closed magnetic
field. This is because when an open field region is sheared by moving the
field line footpoints on the solar surface, the lack of line-tying at the
outer end means that this open field is free to relax back to a minimum-energy
(current-free) configuration. But in a closed-field region the field lines are
line-tied on the photosphere at both ends and will retain the shear.

What prevents free energy building up indefinitely in the corona? There are
three factors to consider (Mikić and Linker, 1994):

  1. i)

As coronal arcades are sheared, they expand and more of the field becomes
open, subsequently relaxing to potential and thus releasing energy.

  2. ii)

If strong enough current sheets form in the corona then these can support
magnetic reconnection allowing reconfiguration of the magnetic field to a
lower energy state (and possibly also the loss of magnetic energy through
ohmic heating). This can be gradual or sudden.

  3. iii)

If the magnetic field becomes too sheared or twisted, it can lose equilibrium
leading to an eruption that removes energy through the outer boundary. This is
thought to be the origin of coronal mass ejections.

One can obtain a feeling for how the free energy varies over the solar cycle
by looking at the observed rates of solar flares (e.g., Hathaway, 2015) or
coronal mass ejections (Webb and Howard, 2012; Lamy et al., 2019). These
indices of solar activity tend to follow (roughly) the observed sunspot
number, because most of these events originate from active regions. For active
regions, non-potential magnetic field extrapolations are now routinely made
(Wiegelmann and Sakurai, 2021), giving local estimates of free energy albeit
with some uncertainty (e.g., Régnier and Priest, 2007; DeRosa et al., 2015;
Thalmann, Gupta, and Veronig, 2022). But the magnetic flux from active regions
spreads out gradually as they decay and interact, and the effect of this
process on free energy is harder to pin down. Indirect observations show that
there is indeed free energy stored outside of active regions, in concentrated
filament channels or even large-scale sheared arcades. But to date there have
been relatively few quantitative studies of how the free energy in the global
solar corona varies over the solar cycle, two exceptions being Yeates,
Constable, and Martens (2010) and Chifu, Inhester, and Wiegelmann (2022).

Chifu, Inhester, and Wiegelmann (2022) used static nonlinear force-free field
(NLFFF) extrapolations from synoptic vector magnetograms, once per 27-day
Carrington rotation, to follow the evolution of the coronal magnetic field
over Solar Cycle 24 (from 2010 to 2019). They found a peak of free energy of
about $2.5\times 10^{33}\,\mathrm{erg}$ – occurring in late 2014 which is the
time when observed magnetic flux and flare activity peaked for that cycle.
These authors did not indicate the relative magnitude of this free energy
compared to the current-free field; however, for a single similar
extrapolation, Tadesse et al. (2014) previously obtained a figure of $15\%$.
On the other hand, this earlier work illustrated – as confirmed in the model
comparison study by Yeates et al. (2018) – that these global NLFFF
extrapolations from present-day vector magnetogram observations are unable to
accurately recover currents outside active regions, such as in filament
channels.

The study by Yeates, Constable, and Martens (2010) used an alternative,
magneto-frictional model with surface driving. This gives less accurate
results for the detailed structure within active regions but is better able to
reproduce the formation of filament channels and other non-potential
structures outside them. A similar model will be used in the present paper and
discussed further in Section 2. Unlike Chifu, Inhester, and Wiegelmann (2022),
the study of Yeates, Constable, and Martens (2010) covered Solar Cycle 23,
looking at six different periods over the cycle. They found that the free
magnetic energy increased by a factor of 15 between solar minimum and solar
maximum, while the total magnetic energy increased by a factor of 8. The
actual value of free energy during Maximum was around $2\times
10^{33}\,\mathrm{erg}$, a little below that found by Chifu, Inhester, and
Wiegelmann (2022) for Solar Cycle 24.

The time-dependent magneto-frictional models also have the advantage of
including the formation and eruption of magnetic flux ropes (Mackay and van
Ballegooijen, 2006; Yeates and Mackay, 2009a). Yeates, Constable, and Martens
(2010) defined flux ropes using the magnetic pressure and tension forces, then
identified ejections by looking at the outward velocity. They found that the
number of ejections increased from Minimum to Maximum by a factor of 8,
proportional to the sunspot number as in the observations of coronal mass
ejections. Lowder and Yeates (2017) analysed a continuous magneto-frictional
simulation from 1996 to 2014, with a new methodology for flux rope
identification based on thresholding of field line helicity. Again, they found
that the eruption rate followed solar activity. In addition, the model
predicted a non-zero number of eruptions during the very weak Minimum between
Cycles 23 and 24, which also seems to be consistent with observations. Both of
these previous magneto-frictional models were driven by idealised bipolar
active regions derived from US National Solar Observatory synoptic
magnetograms.

Whilst the general trend of increasing non-potentiality from Solar Minimum to
Solar Maximum is well established, the aim of this paper is to be more
quantitative. We use a model that is able to approximate the evolution of non-
potential structure on a global scale, rather than treating individual active
regions in isolation. In particular, we update our earlier findings on solar
cycle variation of non-potentiality during Cycle 23 to new magneto-frictional
simulations of Cycle 24 based on HMI data. After describing the model and
implementation (Section 2), we will consider in Section 3 a reference
simulation where all active regions emerge untwisted. Given the presence of
magnetic reconnection and outflow in our model, we would expect this to be a
practical lower bound for the amount of free energy in the corona. In Section
4, we will then discuss the effect of emerging twisted active regions, before
giving our conclusions in Section 5.

This study benefits from two technical advances. The first is the emergence of
active regions with observed shapes rather than idealised magnetic bipoles
(Yeates and Bhowmik, 2022). The surface flux transport study by Yeates (2020b)
showed that in Cycle 24 the idealised bipoles lead to a 24% overestimate of
the overall axial dipole moment, but the effect on the non-potential coronal
magnetic field has not yet been investigated. The second advance is the
development of field line helicity as an additional tool for probing the non-
potential structure of the coronal magnetic field (Berger, 1988; Yeates and
Hornig, 2016; Lowder and Yeates, 2017; Yeates, 2020c; Moraitis, Patsourakos,
and Nindos, 2021).

##  2 Model setup

We simulate the corona using a magneto-frictional model for the mean (large-
scale) magnetic field ${\bm{B}}$, coupled to an evolving surface flux
transport model on the solar surface. This approach was introduced by van
Ballegooijen, Priest, and Mackay (2000) and extended to the global corona by
Yeates, Mackay, and van Ballegooijen (2008). It has since been applied in
numerous studies, including investigations of filament formation (e.g, Mackay,
Gaizauskas, and van Ballegooijen, 2000; Yeates and Mackay, 2009b; Mackay et
al., 2018), coronal mass ejections (e.g., Mackay and van Ballegooijen, 2006;
Yeates and Mackay, 2009a; Lowder and Yeates, 2017; Bhowmik and Yeates, 2021),
sympathetic flares (Schrijver et al., 2013), the middle corona (Meyer et al.,
2020), or the interplanetary magnetic field/solar wind (Yeates et al., 2010;
Edwards et al., 2015). It has also been applied to other stars (e.g., Lehmann
et al., 2017), by other researchers (Hoeksema et al., 2020), and coupled to
full magnetohydrodynamic simulations (e.g., Yardley et al., 2021; Hayashi et
al., 2021). This experience has guided our choice of reference parameter
values; the effect of varying these parameters is considered in Appendix B.

###  2.1 Magneto-frictional model

Our magneto-frictional model for the mean/large-scale magnetic field is based
on the following main assumptions.

  1. i)

The mean magnetic field ${\bm{B}}$ in the corona originates from the emergence
of macroscopic active regions. In this paper we neglect small-scale ephemeral
regions although they do contribute to the total magnetic energy in the real
corona and potentially even to the magnetic helicity on large-scales, through
an inverse cascade (Mackay, DeVore, and Antiochos, 2014). Our treatment of
emerging regions is described in Section 2.4.

  2. ii)

The high magnetic Reynolds number implies a near-ideal evolution of
${\bm{B}}$, so it is important to retain the mean-field induction equation

$\frac{\partial{\bm{B}}}{\partial
t}=-\nabla\times{\bm{E}},\qquad\bm{E}=-{\bm{v}}\times{\bm{B}}+{\bm{N}},$ (1)

where the non-ideal term ${\bm{N}}$ is small. This means in particular that
${\bm{B}}$ holds a memory of previous interactions, allowing topological
structure to build over time.

  3. iii)

The high Alfvén speed in the corona makes ${\bm{B}}$ respond rapidly to
boundary motions on the solar surface (Section 2.3), and flux emergence from
the solar interior (Section 2.4). Accordingly, we neglect plasma forces and
magnetohydrodynamic waves, and approximate the fluid momentum equation with an
artificial velocity

${\bm{v}}=\frac{{\bm{J}}\times{\bm{B}}}{\nu B^{2}}+v_{\rm
out}(r)\hat{\bm{r}},\qquad{\bm{J}}=\nabla\times{\bm{B}}.$ (2)

The first term represents the magneto-frictional assumption (Chodura and
Schlüter, 1981; Yang, Sturrock, and Antiochos, 1986) and ensures relaxation
towards a force-free equilibrium ${\bm{J}}\times{\bm{B}}={\bm{0}}$, albeit
modified here by the second term. We also assume low plasma-beta and neglect
plasma pressure and gravity, although this is a simplification compared to
reality (cf. Chen, Rempel, and Fan, 2022).

  4. iv)

In the outer corona, the solar wind starts to influence ${\bm{B}}$, preventing
it from being force-free. This is accounted for by the radial outflow $v_{\rm
out}(r)=v_{w}(r/R_{1})^{11.5}$ in Equation (2), which models the radially
distending effect of the solar wind on closed-field arcades (Mackay and van
Ballegooijen, 2006), and reduces the solution’s sensitivity to the particular
choice of outer boundary height, $R_{1}$ (Rice and Yeates, 2021). We set
$R_{1}=2.5R_{\odot}$ and fix the outflow speed
$v_{w}=100\,\mathrm{km}\,\mathrm{s}^{-1}$ in our reference model.

  5. v)

The coronal ${\bm{B}}$ relaxes with respect to a rest frame rotating with the
Sun. Our simulations use the Carrington frame.

The non-ideal electric field, ${\bm{N}}$, in Equation (1) represents turbulent
diffusion – the net effect of unresolved small-scale fluctuations on the mean
magnetic field. Here we assume the hyperdiffusive form

${\bm{N}}=-\eta_{h}\frac{{\bm{B}}}{B^{2}}\nabla\cdot\left[B^{2}\nabla\left(\frac{{\bm{J}}\cdot\bm{B}}{B^{2}}\right)\right],$
(3)

which prevents volume dissipation of relative helicity (with respect to a
potential field) since the corresponding volume dissipation term
$\int_{V}{\bm{E}}\cdot{\bm{B}}\,dV$ (Berger and Field, 1984) becomes a
boundary integral (Boozer, 1986; van Ballegooijen and Cranmer, 2008). We set
the reference value $\eta_{h}=10^{11}\,\mathrm{km}^{4}\,\mathrm{s}^{-1}$
following van Ballegooijen and Mackay (2007). The recent study of Mackay and
Upton (2022) suggests that hyperdiffusion gives quite a conservative estimate
of the amount of free energy in the corona, compared to other possible forms
for ${\bm{N}}$, although further observational comparison is needed to be sure
of the most realistic approach.

The “friction” coefficient $\nu$ in equation (2) controls the speed of
relaxation relative to hyperdiffusion and boundary driving. To reduce
computational cost, we follow Mackay, Gaizauskas, and van Ballegooijen (2000)
and reduce the relaxation rate near the poles, setting
$\nu=\nu_{0}(r\cos\lambda)^{-2}$, where $\lambda$ denotes latitude. In our
reference model the amplitude is set to $\nu_{0}=2.8\times
10^{5}\,\mathrm{s}$.

In our numerical implementation, we write ${\bm{B}}=\nabla\times{\bm{A}}$ and
solve the uncurled version of Equation (1) for the vector potential
${\bm{A}}$, using a finite volume method with a staggered mesh so as to
conserve magnetic flux. The mesh has equally spaced points in
$\log(r/R_{\odot})$, sine latitude and longitude, with a resolution of
$(60,180,360)$ cells. The equations are discretised using second-order
“mimetic” spatial differences (discretising Stokes’ Theorem), with a modified
stencil at the polar grid points. Upwinding is used for the radial outflow
term. Time-stepping uses a second-order Runge-Kutta method. The simulations
run from 12 June 2010 to 31 December 2019.

###  2.2 Initial conditions

To initialise the simulation we use a potential field source surface (PFSS)
extrapolation taken from an imposed $B_{r}$ distribution on $r=R_{\odot}$,
assuming ${\bm{J}}={\bm{0}}$ in the corona with
${\bm{B}}\times\hat{\bm{r}}={\bm{0}}$ on the outer source surface boundary,
$r=R_{1}$. For the imposed $B_{r}$ distribution, we use the radial component,
pole-filled HMI map for Carrington rotation CR 2097 from the
hmi.synoptic_mr_polfil_720s series (Sun, 2018). To reduce its resolution to be
comparable to our simulations of the mean field, this map is first smoothed by
multiplying the spherical harmonic coefficients by a Gaussian filter
$\mathrm{e}^{-c_{w}l(l+1)}$, with $c_{w}=5\times 10^{-4}$.

The initial potential field extrapolation is computed using the author’s
finite difference code (Yeates, 2018; Stansby, Yeates, and Badman, 2020),
which ensures that ${\bm{J}}={\bm{0}}$ to machine precision for our particular
discretisation scheme.

###  2.3 Boundary conditions

At the inner boundary the horizontal electric field from Equation (1) is
replaced by

${\bm{E}}_{\perp}\big{|}_{r=R_{\odot}}={\bm{E}}_{\perp}^{\rm
em}-{\bm{v}}_{s}\times\big{(}B_{r}\hat{\bm{r}}\big{)}+\eta_{0}\nabla\times\big{(}B_{r}\hat{\bm{r}}\big{)}.$
(4)

The first term, ${\bm{E}}_{\perp}^{\rm em}$, represents emergence of new
active regions and will be described in Section 2.4. The final two terms
represent, respectively, advection by large-scale surface flows and
supergranular diffusion (the net effect of small-scale surface motions). The
large-scale flow velocity is assumed to take the steady form
${\bm{v}}_{s}=v_{\theta}(\theta)\hat{\bm{\theta}}+R_{\odot}\sin\theta\,\Omega(\theta)\hat{\bm{\phi}}$,
where $v_{\theta}$ denotes the meridional (poleward) flow and $\Omega(\theta)$
is the angular velocity of differential rotation. In this study we neglect any
additional injection of magnetic helicity into the mean field from small-scale
photospheric flows, although this could be included in parameterised form (see
Mackay, DeVore, and Antiochos, 2014).

The surface flux transport parameters ${\bm{v}}_{s}$ and $\eta_{0}$ were
previously calibrated for the same emerging regions by Yeates (2020b) to
ensure a reasonable match to the observed time series of axial dipole
strength. This led to $v_{\theta}(\theta)=D_{u}\cos\theta\sin^{p_{0}}\theta$,
with $p_{0}=2.33$ and $D_{u}=0.041\,\mathrm{km}\,\mathrm{s}^{-1}$ (giving a
peak flow speed of $15\,\mathrm{m}\,\mathrm{s}^{-1}$), and
$\eta_{0}=350\,\mathrm{km}^{2}\,\mathrm{s}^{-1}$. The calibration did not
constrain the differential rotation, so in this paper we use the established
angular velocity profile from Snodgrass and Ulrich (1990) of
$\Omega(\theta)=0.18-2.396\cos^{2}\theta-1.787\cos^{4}\theta$ (degrees per day
in the Carrington frame).

The resulting evolution of $B_{r}$ on the solar surface was illustrated in
Figure 5 of Yeates (2020b). In that paper, it was also determined that if
emerging regions are replaced by idealised BMRs, the meridional flow amplitude
should be increased to $D_{u}=0.055\,\mathrm{km}\,\mathrm{s}^{-1}$ (giving a
peak flow speed of $20\,\mathrm{m}\,\mathrm{s}^{-1}$), to counter the spurious
enhancement in end-of-cycle axial dipole strength resulting from the BMR
approximation. This is adopted in our simulation run BMR0 below.

At both inner and outer boundaries, we impose
${\bm{J}}\times\hat{\bm{r}}={\bm{0}}$, so as to avoid any flux of magnetic
energy into or out of the domain due to the magneto-friction term in Equation
(2). To compute ${\bm{N}}$ on these boundaries, we impose zero-gradient
conditions on the radial component of the hyperdiffusive flux
$B^{2}\nabla({\bm{J}}\cdot{\bm{B}}/B^{2})$.

###  2.4 Emerging regions

We use the dataset of active regions extracted by Yeates (2020b) from
Spaceweather HMI Active Region Patch (SHARP) data (Bobra et al., 2014).
Briefly, a single line-of-sight magnetogram for each region is selected at the
time closest to central meridian, and remapped to the simulation grid. The
dataset is filtered to remove regions that are (i) too unbalanced in flux;
(ii) too small to resolve; or (iii) repeat observations of the same region on
an earlier disk passage. This leaves surface $B_{r}$ maps for 1072 regions
during the simulation period. Figure 1(a) shows the time-latitude distribution
of regions, coloured by their unsigned magnetic flux, $\Phi_{0}$ (including
both polarities). The distribution of region fluxes is shown in Figure 1(b).

Figure 1.: Emerging regions determined from the HMI/SHARPs data, shown in a
time-latitude “butterfly” plot and coloured by either (a) unsigned flux,
$\Phi_{0}$ or (c) observed twist, $\alpha_{0}$. Panel (b) shows a histogram of
the $\Phi_{0}$ values, while (d) shows the mean and standard deviation of the
$\alpha_{0}$ values in equal sine-latitude bins. There is a hemispheric
tendency in $\alpha_{0}$ with
$\langle\alpha_{0}\rangle_{\lambda_{0}>0}=-0.0147\,\mathrm{Mm}^{-1}$ and
$\langle\alpha_{0}\rangle_{\lambda_{0}<0}=0.00510\,\mathrm{Mm}^{-1}$. The
percentage of regions obeying the hemispheric rule is $75.9\%$ in the North
and $60.6\%$ in the South.

Figure 1(c) shows the same SHARP regions but coloured instead by the observed
“twist” parameter $\alpha_{0}=\sum J_{z}B_{z}/(\sum B_{z}^{2})$ from the SHARP
metadata (MEANALP in Table 3 of Bobra et al., 2014). Here the sum is taken
over all observed pixels within the SHARP mask. This incorporates vector
magnetogram information through $J_{z}$, and will be used here to set the
helicity of emerging regions in two simulation runs. There is considerable
scatter in the observed $\alpha_{0}$ values. Some of this is likely real –
resulting from the variation in helicity between different active regions
(e.g., Georgoulis et al., 2009) – but some is likely due to measurement error.
In addition, it must be remembered that the $\alpha_{0}$ parameter does not
correspond directly with the full current helicity, nor the magnetic helicity
which is the underlying conserved quantity (cf. Russell et al., 2019).
Nevertheless, the latitude distribution in Figure 1(d) shows that $\alpha_{0}$
recovers the expected hemispheric helicity trend of negative in the North and
positive in the South (Pevtsov and Balasubramaniam, 2003; Liu, Hoeksema, and
Sun, 2014).

Our technique for emerging each region in the three-dimensional simulation was
motivated, described, and illustrated in detail by Yeates and Bhowmik (2022).
Briefly, for a fixed time interval $T_{\rm em}=24\,\mathrm{hr}$ we impose a
steady (surface) electric field

${\bm{E}}^{\rm
em}_{\perp}(\theta,\phi)={\bm{E}}^{0}_{\perp}(\theta,\phi)-\nabla\frac{\Psi(\theta,\phi)}{T_{\rm
em}}.$ (5)

The first term ${\bm{E}}^{0}_{\perp}$ is the local inductive electric field,
which generates the given $B_{r}$ distribution while producing a coronal
magnetic field with minimum possible complexity (close to potential). This is
termed “local” because we fix ${\bm{E}}^{0}_{\perp}={\bm{0}}$ outside of the
emerging region.

The second term in (5) does not change the final $B_{r}$ distribution but
generates additional helicity in the three-dimensional field, allowing us to
account for active regions that emerge in a non-potential state. In the
reference run T0, this additional twisting is omitted, generating coronal
active regions that are close to potential with low internal helicity. An
example is shown in Figure 2(a).

In runs where additional twisting is included, we follow Yeates and Bhowmik
(2022) and set the twisting potential to

$\Psi(\theta,\phi)=\tau b_{0}\langle f_{\rm pil}\rangle\overline{B}_{r}.$ (6)

Here $\overline{B}_{r}$ is a smoothed magnetogram used to identify the
polarity inversion lines, and $\langle f_{\rm pil}\rangle$ is a function which
localises the twist near to them. The normalisation factor $b_{0}$ is chosen
to ensure that the horizontal magnetic field generated by twisting is
$\approx|\tau B_{r}|$. The single dimensionless “twist” parameter $\tau$ then
controls the helicity injection within the emerging region, with $\tau>0$
leading to positive helicity and $\tau<0$ to negative helicity. Figure 2(b)
shows the same region as Figure 2(a), but emerged with $\tau=-0.1$ giving it
negative helicity – a reverse-S shaped twist.

Since the optimum value of $\tau$ to model each active region is unknown, we
follow Yeates and Bhowmik (2022) and vary $\tau$ in a parameter study. The
simulation runs are summarised in Table 1. In TU0.05 and TU0.1, we simply fix
a uniform $\tau$ value for all regions in each hemisphere. In some sense this
maximises the energisation since it minimises cancellation between positive
and negative helicity during the relaxation. The largest value of $|\tau|=0.1$
is chosen because larger values lead to unrealistically twisted coronal field
lines (Yeates and Bhowmik, 2022). In order to consider the effect of varying
helicity between active regions, runs TOb5 and TOb10 use different values of
$\tau$ for each emerging region, chosen proportional to the observed twist
values $\alpha_{0}$. Since the dimensions of $\tau$ and $\alpha_{0}$ are
different, the proportionality constant $\ell$ in Table 1 has dimensions of
length. The values of $\ell$ are chosen so that the two runs have a comparable
overall helicity to the two runs with uniform $\tau$.

Figure 2.: Three-dimensional renderings of an emerging region (SHARP 187/NOAA
11109 on 2010-09-28), immediately after completion of emergence in three of
the simulation runs. The twist parameter is $\tau=0$ in (a,c) and $\tau=-0.1$
in (b). Grayscale pixels show $B_{r}$ on $r=R_{\odot}$ (white positive/black
negative, saturated at $\pm 25\,\mathrm{G}$), while magnetic field lines are
coloured by field line helicity (red positive/blue negative, saturated at $\pm
10^{22}\,\mathrm{Mx}$).

Finally, for most runs, we use the observed SHARP active region shapes. For
run BMR0, the regions are replaced by equivalent, idealised BMRs (bipolar
magnetic regions), according to the prescription in Yeates (2020b). We take
these to be untwisted, as illustrated in Figure 2(c).

Table 1.:  Model runs with different emerging regions. Coronal parameters take values in all of these runs, namely $\nu_{0}=2.8\times 10^{5}\,\mathrm{s}$, $\eta_{h}=10^{11}\,\mathrm{km}^{4}\,\mathrm{s}^{-1}$, $v_{w}=100\,\mathrm{km}\,\mathrm{s}^{-1}$.  Run | Emerging Regions | Twist | SFT Parameters |
---|---|---|---|---
|  | $\tau$ |  $\eta_{0}$ [$\mathrm{km}^{2}\,\mathrm{s}^{-1}$] |  $D_{u}$ [$\mathrm{km}\,\mathrm{s}^{-1}$] | $p_{0}$ |
T0 | Observed | 0.00 | 350 | 0.041 | 2.33 |
TU0.05 | Observed | $\pm 0.05$ | 350 | 0.041 | 2.33 |
TU0.1 | Observed | $\pm 0.10$ | 350 | 0.041 | 2.33 |
TOb5 | Observed |  $\ell\alpha_{0}$, $\ell=5\,\mathrm{Mm}$ | 350 | 0.041 | 2.33 |
TOb10 | Observed |  $\ell\alpha_{0}$, $\ell=10\,\mathrm{Mm}$ | 350 | 0.041 | 2.33 |
BMR0 | Idealised BMRs | 0.00 | 350 | 0.055 | 2.33 |

###  2.5 Cycle variation of the photospheric magnetic field

The distribution of $B_{r}$ on the solar surface is unaffected by the choice
of coronal parameters or emerging region twist. The longitude-averaged $B_{r}$
was shown as a function of time and latitude in Figure 5 of Yeates (2020b).
Here, in Figures 3(a-c), we show the emergence rate (number of emerging
regions in 27-day bins), total unsigned magnetic flux,

$\Phi_{0}=\oint_{r=R_{\odot}}|B_{r}|\,\mathrm{d}S,$ (7)

and total polarity inversion line (PIL) length. The darker curves show the
reference run T0 (or equivalently TU0.05, TU0.1, TOb5, or TOb10). The feinter
curves show run BMR0. This has the same emergence rate. However, whilst the
BMR approximation does not change $\Phi_{0}$ for individual regions, it does
slightly increase $\Phi_{0}$ around Solar Maximum because of differences in
flux cancellation between regions, and the modified $B_{r}$ pattern also
slightly increases the PIL length.

Figure 3.: Cycle variation of the solar surface magnetic field in runs T0
(original region shapes) and BMR0 (idealised bipoles). Panels (a-c) show the
quantities with their original units, while panel (d) shows all curves
normalised to their mean value during 2018 and 2019 (to right of the vertical
dashed line). Grey shading shows the 13-month smoothed monthly International
Sunspot Number, normalised in the same way. (Source: WDC-SILSO, Royal
Observatory of Belgium, Brussels.)

Figure 3(d) shows the same three quantities but normalised to their mean
values during 2018 and 2019 (to the right of the vertical dashed line). This
period is chosen to represent Solar Minimum rather than the start of the
simulation because the HMI dataset driving the simulation does not extend all
the way back to Solar Minimum. Similar to Yeates, Constable, and Martens
(2010), we observe that (i) $\Phi_{0}$ is about 4-5 times higher at Maximum
than Minimum, and (ii) the emergence rate in the simulation is about 7-8 times
higher at Maximum than Minimum. Note that the increase in observed sunspot
number is rather higher than in the previous study, perhaps 20 times more at
Maximum than Minimum. This is likely because the SHARP regions can contain
multiple spots.

At the resolution considered, the PIL length shows a much more modest increase
of less than a factor 2 from Minimum to Maximum, although the length
fluctuates substantially over 2018 and 2019. This modest increase is again in
line with Yeates, Constable, and Martens (2010).

##  3 Reference simulation

In this section we analyse the solar cycle variation of the reference
simulation: run T0 in Table 1. This represents a conservative model where all
active regions emerge with no internal twist/helicity. As evidenced by the
selected snapshots in Figure 4, this still leads to a corona that is non-
potential with concentrations of magnetic helicity in sheared arcades or
twisted flux ropes (cf. Yeates and Hornig, 2016; Bhowmik and Yeates, 2021).
The associated free energy can build up at active latitudes, as in Figure
4(c), but it can also be stored in closed-field arcades at higher latitudes,
as seen on the Northern polar crown in Figure 4(d). In the following
subsections we will consider how quantitative properties of the simulated
magnetic field vary between Solar Minimum and Solar Maximum.

Figure 4.: Four snapshots from the reference simulation, run T0, viewed from
the direction of Earth at the corresponding time. Magnetic field lines are
coloured by their field line helicity (Section 3.2): red positive/blue
negative, saturated at $\pm 10^{22}\,\mathrm{Mx}$. Gray shading shows $B_{r}$
on the solar surface, black negative/white positive, saturated at $\pm
15\,\mathrm{G}$.

###  3.1 Energy

Figure 5(a) shows total magnetic energy,

$E=\int_{V}\frac{B^{2}}{8\pi}\,\mathrm{d}V,$ (8)

along with the energy of a corresponding PFSS extrapolation. The cyan curve
shows free energy,

$E_{\rm free}=E-E_{\rm p},\qquad E_{\rm p}=\int_{V}\frac{(B_{\rm
p})^{2}}{8\pi}\,\mathrm{d}V,$ (9)

where ${\bm{B}}_{\rm p}$ is the current-free magnetic field in the coronal
volume matching $B_{{\rm p}r}=B_{r}$ on both $r=R_{\odot}$ and $r=R_{1}$.
Whilst this reference potential field technically differs from the PFSS field
shown in Figure 5(a), their energies are indistinguishable on the scale of the
plot. The peak around late 2014/early 2015 arises from a particularly strong
activity complex – the “Great Solar Active Region” NOAA 12192 (Sun et al.,
2015) which was the largest since 1990. This already causes a substantial peak
of magnetic energy in the PFSS model, and coincides with the peak of solar
activity as noted by Chifu, Inhester, and Wiegelmann (2022).

To measure the overall energisation, Figure 5(b) shows relative free energy,
$E_{\rm free}/E_{\rm p}$. This fluctuates from one snapshot to the next, but
on average is about $13.8\%$ over much of the cycle. This is rather lower than
many of the non-potential models compared by Yeates et al. (2018), where the
ratio $E_{\rm free}/E_{\rm p}$ was in the range $40-50\%$. However, we will
see in Section 4 that $E_{\rm free}$ is substantially increased if the active
regions emerge twisted.

Figure 5.: Cycle variation of the coronal magnetic energy in run T0, computed
at 27-day intervals. Dark blue curves show the total energy, $E$, blue-grey
the energy of a corresponding PFSS extrapolation, $E_{\rm p}$, and cyan the
free energy, $E_{\rm free}$. Panel (a) shows the original values while (b)
shows the relative free energy, $E_{\rm free}/E_{\rm p}$. Panel (c) shows the
same curves as (a) but normalised as in Figure 3 and compared to the observed
sunspot number.

Figure 5(c) shows how the three energies vary over the solar cycle, compared
to the 2018-2019 period. Since energy is quadratic in ${\bm{B}}$, by comparing
to the variation of photospheric flux in Figure 3(d) one would expect roughly
a factor $16-25$ increase from Minimum to Maximum, and this is indeed the
order of magnitude found. Note that $E$ and $E_{\rm p}$ show a similar
relative increase, as does $E_{\rm free}$ albeit a little higher. However, one
has to bear in mind that these are snapshots rather than running means; in
particular, since $E_{\rm free}$ tends to lag $E_{\rm p}$ (because of the time
taken for differential rotation to shear the field), the two energies – which
are fluctuating on a faster timescale – are not necessarily sampled in phase.
Thus we should not read too much into the precise heights of the peaks.

We note that Yeates, Constable, and Martens (2010) found $E_{\rm free}$ to
increase rather more from Minimum to Maximum compared to $E$, but here we find
it much closer. However, this is again an effect of assuming the emerging
regions to be untwisted, as we will see in Section 4.

###  3.2 Current and helicity

Free energy is not the only way to characterise the complexity of the coronal
magnetic field. Figure 6 shows two further measures: mean current density and
unsigned helicity.

Figure 6.: Cycle variation of (a) the mean current density $\langle J\rangle$
and (b) the unsigned helicity $\overline{H}$, in run T0, computed at 27-day
intervals. Panel (c) shows the same curves as (a) and (b) but normalised as in
Figure 3 and compared to the observed sunspot number. The feint line in (b)
shows $\overline{H}$ for corresponding PFSS extrapolations from the same
surface $B_{r}$ distribution.

Mean current density – shown in Figure 6(a) is defined as

$\langle{J\rangle}=\frac{\int_{V}|\nabla\times{\bm{B}}|\,\mathrm{d}V}{\int_{V}\,\mathrm{d}V},$
(10)

where $V$ is the full simulation volume. Figure 6(c) shows that it follows a
similar 4-5 fold increase to the magnetic flux from Figure 3.

Unsigned helicity – shown in Figure 6(b) – is a measure of the topological
complexity of the coronal magnetic field. It is defined as

$\overline{H}=\frac{1}{2}\int_{\partial
V}\big{|}\mathcal{A}B_{r}\big{|}\,\mathrm{d}S,$ (11)

where $\mathcal{A}$ is the field line helicity

$\mathcal{A}(L)=\int_{L}{\bm{A}}^{*}\cdot\,\mathrm{d}{\bm{l}},$ (12)

for a magnetic field line $L$ traced from a given point on $\partial V$, and
${\bm{A}}^{*}$ is an appropriate vector potential. Figure 6(c) shows a similar
rate of increase in $\overline{H}$ from Minimum to Maximum compared to
$\langle J\rangle$, albeit a rather different time profile, relatively lower
in the increasing phase of the cycle. Like magnetic energy, $\overline{H}$
also shows a peak in late 2014. This peak is present even in the PFSS model
(as explained by Yeates, 2020c), but is significantly enhanced in the non-
potential model – even in run T0 – because the strong active region acts a
seed for shearing by differential rotation.

Our rationale for using $\overline{H}$ rather than the simpler-to-compute
$\int_{V}|{\bm{A}}\cdot{\bm{B}}|\,\mathrm{d}V$ is that $\overline{H}$ is a
topological invariant: it is unchanged under ideal deformations within the
domain that fix the field line footpoints on the boundary, whereas the latter
is not (and neither is $\langle J\rangle$). Thus it is a more robust measure
than $\int_{V}|{\bm{A}}\cdot{\bm{B}}|\,\mathrm{d}V$. Our rationale for using
$\overline{H}$ rather than the commonly-used relative helicity invariant
(Berger and Field, 1984) is the fact that relative helicity is a signed
quantity. It would be a poor measure of topology for our global model because
it would incorporate cancellation between regions of positive and negative
helicity, even if these regions are remotely located on the Sun and do not
interact in reality. This departure from what has become routine practice in
solar physics is further discussed in Appendix A.

The specific values of $\mathcal{A}(L)$ and hence of $\overline{H}$ depend on
the chosen gauge of ${\bm{A}}^{*}$. However, recent work suggests that there
is a meaningful “canonical” choice of gauge (Berger and Hornig, 2018; Yeates,
2020c; Xiao, Prior, and Yeates, 2023), having the form
${\bm{A}}^{*}=\nabla\times(P\hat{\bm{r}})+T\hat{\bm{r}}$ for suitable $P$ and
$T$. Since ${\bm{A}}^{*}$ differs from the gauge of ${\bm{A}}$ used to solve
Equation (1), we compute ${\bm{A}}^{*}$ numerically to find $\mathcal{A}$ and
$\overline{H}$. Whilst not unique, using this gauge consistently for all of
our calculations allows us to compare different time snapshots and different
runs. Moreover, as we will see in Section 4, this gauge choice leads to a
meaningful correlation between $\overline{H}$ and free energy.

Note that even potential fields can have non-zero $\overline{H}$ if their
boundary $B_{r}$ distribution lacks axisymmetry – this behaviour was
demonstrated for typical solar-like configurations in Yeates (2020c). For run
T0, the feint curve in Figure 6(b) shows that $\overline{H}$ in the PFSS model
is approximately $33\%$ of $\overline{H}$, following a similar cycle trend.
Correlation between them arises because sub-structure in the PFSS field acts
as a seed for subsequent shearing by differential rotation, which amplifies
$\overline{H}$ in run T0. However, $\overline{H}$ also fluctuates more rapidly
in the magneto-frictional simulation due to sudden reconfigurations such as
flux rope ejections.

###  3.3 Eruptivity

In this paper we avoid identifying individual magnetic flux ropes, as this is
a significant task in itself and the statistics would depend on the method
used and definition of a flux rope (but for two possible approaches see Yeates
and Mackay, 2009a; Lowder and Yeates, 2017). Instead, we have identified a
proxy for the eruption rate that can be computed without the need to identify
individual structures or track them over time. Importantly, this proxy can be
computed from global diagnostics which can be saved at high time cadence,
rather than requiring three-dimensional magnetic snapshots.

Our chosen eruptivity proxy is the second time derivative of unsigned open
magnetic flux,

$\Phi_{1}=\oint_{r=R_{1}}|B_{r}|\,\mathrm{d}S.$ (13)

As illustrated by Bhowmik and Yeates (2021), eruptions of individual
structures in the model lead to temporary enhancements in $\Phi_{1}$. Taking
the second time derivative, $\ddot{\Phi}_{1}$, removes the underlying secular
variation and extracts these peaks. The evolution of this quantity is
illustrated for run T0 in Figures 7(a-d).

Figure 7.: Time series used to compute the eruptivity proxy,
$\langle\ddot{\Phi}_{1}\rangle$, shown for run T0. Panel (a) shows the open
magnetic flux, $\Phi_{1}$ over the full simulation, with enlargement of a
selected period in (b). Panels (c-d) show the second time derivative,
$\ddot{\Phi}_{1}$. For comparison, panels (e-f) show the Poynting flux,
$S_{1}$, through the outer boundary, with dots in (f) indicating peaks. Panel
(g) shows the 27-day mean $\langle\ddot{\Phi}_{1}\rangle$, overlayed with the
number of $S_{1}$ peaks in 27-day bins, both normalised to their maximum
values. Panel (h) shows $\langle\ddot{\Phi}_{1}\rangle$ normalised as in
Figure 3 and compared to the observed sunspot number.

Another signature of eruptions are short “bursts” in the Poynting flux of
magnetic energy through the outer boundary. From Equations (1)-(3), and our
boundary conditions, this takes the form

$S_{1}=-\frac{1}{4\pi}\oint_{r=R_{1}}{\bm{E}}\times{\bm{B}}\cdot\hat{\bm{r}}\,dS=-\frac{v_{\rm
out}(R_{1})}{4\pi}\oint_{r=R_{1}}|{\bm{B}}\times\hat{\bm{r}}|^{2}\,dS.$ (14)

Due to the radial outflow, the energy flux is purely outward (negative).
Moreover, it depends only on the horizontal magnetic field, so is negligible
in equilibrium, spiking only during eruptive events. It is shown for run T0 in
Figures 7(e-f). Comparing the enlarged Figures 7(d) and (e) shows that the
timings of peaks in $S_{1}$ (labelled with small circles) correspond well to
peaks in $\ddot{\Phi}_{1}$.

For clearer comparison over the full solar cycle, we take our eruptivity proxy
to be the 27-day average, $\langle\ddot{\Phi}_{1}\rangle$. This is shown by
the green curve in Figure 7(g). The purple curve shows the number of peaks in
$S_{1}$ in 27-day bins, which is seen to follow a similar variation over the
cycle. The peaks were counted requiring a mininium-to-maximum height of at
least the mean value $2\times 10^{24}\,\mathrm{erg}\,\mathrm{s}^{-1}$, but
varying this threshold changes only the overall normalisation in Figure 7(g)
and not the relative variation over the cycle. Thus
$\langle\ddot{\Phi}_{1}\rangle$ gives a meaningful proxy for the normalised
eruption rate. Note that $S_{1}$ itself would peak more strongly at Maximum
because it is quadratic in ${\bf B}$. Our chosen measure instead counts the
rate of eruptions, with less regard to size. This is simpler to compare with
observations.

For run T0, Figure 7(h) shows the cycle variation of
$\langle\ddot{\Phi}_{1}\rangle$ when normalised by its values in 2018 and
2019, as for previous quantities. We observe an increase from Minimum to
Maximum by a factor comparable with $\Phi_{0}$, $\langle J\rangle$, or
$\overline{H}$, but rather less than $E_{\rm free}$. This contrasts with
Yeates, Constable, and Martens (2010) where the rate of eruptions was found to
increase more substantially over the cycle, more akin to the sunspot number or
total energy. Furthermore, $\langle\ddot{\Phi}_{1}\rangle$ does not follow the
shape of the sunspot curve, peaking in late 2014 but weaker before 2014. In
particular, the eruptions themselves do not correlate with the timing of
individual active region emergences – some are caused more-or-less directly by
emergence, but others arise from the longer term shearing of remnant flux, as
in previous magneto-frictional models (e.g. Mackay and van Ballegooijen,
2006). However, we will see in the next section that the eruptivity is also
affected by the helicity of emerging regions.

##  4 Effect of emerging region properties

We have repeated the analysis of coronal properties for each simulation run in
Table 1.

###  4.1 Emerging region twist – cycle variation

To investigate the effect of injecting helicity in the emerging regions, we
compare runs TU0.05, TU0.1, TOb5 and TOb10 to the reference run T0. All of
these runs share the same evolution of $B_{r}$ on the solar surface as run T0,
with the only difference being the non-zero twist parameters $\tau$.

Figure 8 compares the time evolution of $E_{\rm free}$ and $\overline{H}$ for
the different runs. Both quantities are enhanced in all of the runs with
$\tau\neq 0$, compared to T0. Again, this is clearly related to periods of
active region emergence, as one would expect given that the additional free
energy is coming from the twisting of emerging regions. Figure 8(a) for
$E_{\rm free}$ can be compared with the NLFFF extrapolations in Figure 8 of
Chifu, Inhester, and Wiegelmann (2022). Even in run TOb10, our $E_{\rm free}$
looks to be perhaps $50\%$ lower on average, likely due at least in part to
the lower resolution – Chifu, Inhester, and Wiegelmann (2022) use
(180,280,720) cells for analysis, and their $E_{\rm free}$ peaks at about
$2.5\times 10^{33}\,\mathrm{erg}$ in late 2014. On the other hand we do see a
similar overall pattern of short-term enhancements, with most of our models
peaking similarly in late 2014. Our $E_{\rm free}$ shows sharper fluctuations,
during which the magnitude can be comparable to Chifu, Inhester, and
Wiegelmann (2022). On the other hand, our $E_{\rm free}$ values around March
2015 look broadly consistent with those obtained by Mackay and Upton (2022) in
their magneto-frictional runs with hyperdiffusion.

Figure 8.: Time evolution of (a) $E_{\rm free}$ and (b) $\overline{H}$, in the
coronal simulations with different emerging region properties. Here colours
distinguish the different runs, with symbols indicating the peaks. The curves
for run T0 are repeated from Figures 5(a) and 6(b).

Figures 9(b,c,e,f) show the normalised cycle variation of $E_{\rm free}$ and
$\overline{H}$ for each run, as well as the eruptivity
$\langle\ddot{\Phi}_{1}\rangle$. Again we see that, as we increase $|\tau|$ –
either the uniform value in runs TU0.05 and TU0.1 or the factor multiplying
$\alpha_{0}$ for each region in runs TOb5 and TOb10 – there is a greater
contrast in non-potentiality of the corona from Minimum to Maximum. This is
particularly seen in $E_{\rm free}$, which now shows a greater increase than
the sunspot number. The pattern of variation in
$\langle\ddot{\Phi}_{1}\rangle$ is also brought more into line with that of
the sunspot number, in particular being enhanced during the rising phase of
the cycle. Both of these findings bring the results into line with Yeates,
Constable, and Martens (2010), removing the discrepancies seen in run T0.

Figure 9.: Cycle variation of parameters in the coronal simulations with
different emerging region properties. Each panel corresponds to a different
simulation run (Table 1), with all curves normalised as in Figure 3 and
compared to the observed sunspot number.

###  4.2 Emerging region twist – cycle averages

Figure 10 shows some interesting scatter plots for the simulations in Table 1,
namely cycle-averaged free energy and eruptivity against cycle-averaged
$\overline{H}$. In Figure 10(a), we see a clear power-law relation between the
cycle-averages of $E_{\rm free}$ and $\overline{H}$, which holds for all of
the runs (including BMR0). It is also clear how both quantities increase as
the emerging region twist $\tau$ is increased. The intermediate runs TU0.05
and TOb5 are roughly comparable to one another in both $\overline{H}$ and
$E_{\rm free}$. Compared to the reference run T0, these two runs have
approximately $20\%$ more unsigned helicity and $55\%$ more free energy.

Figure 10.: Relation between cycle-averaged unsigned helicity, $\overline{H}$
and (a) free energy, $E_{\rm free}$; or (b) eruptivity,
$\langle\ddot{\Phi}_{1}\rangle$. Each symbol denotes a run with different
emerging region properties. In (a) the dashed line shows a power-law fit with
the given coefficients (in the indicated units), while in (b) separate dashed
lines guide the eye between runs with uniform $\tau$ (red) and
$\tau\propto\alpha_{0}$ (blue).

A similar relation between free energy and magnetic helicity has been found
for data-driven force-free models of solar active regions (Tziotziou,
Georgoulis, and Raouafi, 2012). Those authors found a fit of
$|H_{R}|=1.37\times 10^{14}\,E_{\rm free}^{0.897}$, where $E_{\rm free}$ is in
ergs and $H_{R}$ is the usual relative magnetic helicity in $\mathrm{Mx}^{2}$.
If this is extrapolated to a global value of $H_{R}=1.5\times
10^{44}\,\mathrm{Mx}^{2}$ – approximately the cycle mean of $\overline{H}$ for
run TU0.05 – it predicts $E_{\rm free}\approx 3\times 10^{33}\,\mathrm{erg}$,
which is an order of magnitude greater than our simulations. One possible
explanation for the difference is the fact that $|H_{R}|\leq\overline{H}$,
with equality only if $\mathcal{A}$ is uniform in sign throughout the region.
Another possibility is that active regions have free energy further above the
helicity threshold than the corona does as a whole.

Figure 11.: Relation between unsigned helicity, $\overline{H}$ and free
energy, $E_{\rm free}$, for individual time snapshots with no cycle averaging.
Each symbol denotes a run with different emerging region properties,
separately for (a) runs where all regions have the same twist and (b) runs
with varying amplitudes and signs of twist.

Figure 11 shows that when the time averaging is removed from Figure 10(a),
there is more scatter in the relation between unsigned helicity and free
energy. For run T0 the relation is still approximately linear for all time
snapshots, but when emerging regions are twisted, there is more variation. In
this case, the value for run T0 seems to act as an approximate lower bound on
the free energy, at a given level of unsigned helicity.

Figure 10(b) shows that the eruptivity seems to depend not only on
$\overline{H}$ but also on whether all regions have the same twist (runs
TU0.05 and TU0.1) or have varying amplitudes and signs of twist (runs TOb5 and
TOb10). Specifically, for a given overall unsigned helicity, there are more
eruptions when $\tau$ varies between regions rather than being uniform. This
effect is also visible by comparing Figures 9(b) and 9(e). We speculate that
the cause could be the distribution of many different magnitudes $|\tau|$ (of
both signs), with a sizeable tail of regions having $|\tau|>0.1$, namely 714
out of 1072 regions in run TOb10, even though the overall average helicity is
less than run TU0.1 where all regions have $|\tau|=0.1$. The most twisted of
these regions are likely to erupt rapidly, before the magnetic field has a
chance to relax or dissipate strong currents.

###  4.3 Idealised bipole approximation

Figures 9 and 10 also show the run BMR0 where the emerging regions were
replaced with idealised bipolar magnetic regions (BMRs). The corresponding
values of $\overline{H}$ and $E_{\rm free}$ are approximately $8\%$ and $18\%$
higher than run T0, and this increase is consistent with the increase in
surface flux $\Phi_{0}$ seen in Figure 3(b). In Figure 10(b), we see that
despite the increase $\overline{H}$ the eruptivity
$\langle\ddot{\Phi}_{1}\rangle$ in run BMR0 is comparable to run T0. This
suggests that the additional helicity is well distributed and not leading to
the more concentrated flux rope structures that would generate more eruptions.

##  5 Conclusions

In this paper we have presented a picture of how the global coronal magnetic
field may have varied over Solar Cycle 24, based on a simplified model that
nevertheless retains some of the key physics discarded by the PFSS model. Our
model confirms the finding of Chifu, Inhester, and Wiegelmann (2022) that free
energy peaked in late 2014, and we have shown that this is also true of
(unsigned) magnetic helicity. It also confirms the result of Yeates,
Constable, and Martens (2010) that the magneto-frictional model can reproduce
the observed pattern of eruptivity varying roughly like the sunspot number. On
the other hand, we have shown that this pattern of eruptivity relies on active
regions emerging twisted (i.e., with magnetic helicity). This is consistent
with previous modelling of active region formation by emerging magnetic flux
tubes, which are found to emerge more readily if they are twisted (Cheung and
Isobe, 2014). Interestingly, Mackay and Upton (2022) show that one can obtain
similar coronal energisation if there is sufficient small-scale helicity
injection which accumulates along PILs in the mean field.

We caution that the actual numbers obtained for many of the quantities in our
study are still likely to be model dependent. The comparison study by Yeates
et al. (2018) makes clear that different non-potential models have widely
differing $E_{\rm free}/E_{\rm p}$, due to the inclusion of different forms of
electric current. Even within the scope of the magneto-frictional model,
Mackay and Upton (2022) have shown how different formulations of the non-ideal
term ${\bm{N}}$ can have a substantial effect on the free energy, with fully
ideal simulations (${\bm{N}}={\bm{0}}$) having $E_{\rm free}/E_{\rm p}$ in
excess of $60\%$. Compared to this, our simulations are rather conservative in
the amount of free energy (the cycle-average ratio $E_{\rm free}/E_{\rm p}$ is
about $20\%$ in TU0.05 or TOb5, and about $25\%$ in TU0.1 or TOb10). Mackay
and Upton (2022) also showed that assimilating active regions at their time of
maximum flux rather than central meridian crossing leads to substantial
increases in $\Phi_{0}$ and $E$, though not always in $E_{\rm free}/E_{\rm
p}$.

Another problem is numerical resolution. Appendix B.1 shows that halving the
resolution of our simulations still gives reasonable estimates of helicity and
free energy in 2014, though it underestimates eruptivity. But in general, the
NLFFF modelling of active regions shows that extrapolations at higher
resolution typically contain more free energy (DeRosa et al., 2015; Thalmann,
Gupta, and Veronig, 2022).

In future, more stringent calibration against observations will likely help to
choose between the different models. Other than the surface magnetic field
(which was calibrated by Yeates, 2020b), this was beyond the scope of this
study because (a) the comparison is necessarily indirect owing to the lack of
coronal magnetic observations, so it is difficult to be quantitative; and (b)
optimising the coronal parameters in the model is computationally very
expensive. Therefore, whilst our model is based on physical principles, it is
still only indicative of likely trends. Possible contraints that have
previously been used – albeit not for a full solar cycle – include sheared or
twisted fields in filament channels (Yeates, Mackay, and van Ballegooijen,
2008; Mackay et al., 2018; Mackay and Upton, 2022), the shapes of coronal
streamers in extreme-ultraviolet (Meyer et al., 2020; Mackay and Upton, 2022;
Wagner et al., 2022), the locations of coronal holes (Yeates et al., 2018),
the total open/heliospheric flux (Linker et al., 2017), or even the timing of
specific eruptions (Yardley, Mackay, and Green, 2021). We hope to pursue this
in future.

#### Acknowledgements

This work was supported by UK STFC grants ST/W00108X/1 and ST/V00235X/1. The
SDO data are courtesy of NASA and the SDO/HMI science team. Figures 2 and 4
used PyVista (Sullivan and Kaszynski, 2019). The author thanks the anonymous
reviewer for useful suggestions.

## Data Availability

The database of emerging regions is available as a pre-generated ascii file on
the Harvard Dataverse (Yeates, 2020a). This may be reproduced with either the
same or different parameters using our open-source Python code at
https://github.com/antyeates1983/sharps-bmrs. The author will be happy to
share data from the magneto-frictional simulations upon reasonable request.

##  Appendix A Unsuitability of relative helicity

In solar physics, the relative helicity, $H_{R}$, of Berger and Field (1984)
is frequently used as a measure of global magnetic topology. Like the unsigned
helicity $\overline{H}$ used in this paper, it is unchanged under ideal
deformations within the domain that fix the field line footpoints on the
boundary. Unlike $\overline{H}$, it has the advantage of being a gauge
invariant quantity, but it replaces gauge dependence by dependence on a choice
of reference field. It is conventional, however, to choose a potential field
as the reference.

Unfortunately, $H_{R}$ is not a useful measure for global simulations of the
solar corona, because it is a signed quantity that suffers from cancellation
between positive and negative values of the integrand. These can come from
regions that are remotely located on the Sun and do not interact in reality.
To illustrate this, Figure 12 repeats Figure 11, but using $|H_{R}|$ on the
horizontal axis instead of $\overline{H}$. There is greatly increased scatter
compared to Figure 11, and, in particular, $E_{\rm free}$ does not go to zero
as $|H_{R}|$ goes to zero, unlike the behaviour with $\overline{H}$. This
highlights how a low value of $|H_{R}|$ can result from cancellation of
opposite signs and need not indicate a corona near to potential.

Figure 12.: Relation between (unsigned) relative helicity, $H_{R}$ and free
energy, $E_{\rm free}$, for individual time snapshots with no cycle averaging.
As in Figure 11, Each symbol denotes a run with different emerging region
properties, separately for (a) runs where all regions have the same twist and
(b) runs with varying amplitudes and signs of twist.

Relative helicity is most commonly computed using the Finn and Antonsen (1985)
formula,

$\displaystyle H_{R}=\int_{V}({\bf A}+{\bf A}_{\rm ref})\cdot({\bf B}-{\bf
B}_{\rm ref})\,\mathrm{d}V,$ (15)

where ${\bf B}_{\rm ref}=\nabla\times{\bf A}_{\rm ref}$ is the reference
field. Since $H_{R}$ is gauge invariant for both ${\bf A}$ and ${\bf A}_{\rm
ref}$, if we use the canonical gauge ${\bf A}={\bf A}^{*}$ from Section 3.2
and choose the gauge of ${\bf A}_{\rm ref}$ such that ${\bf A}^{*}\times{\bf
A}_{\rm ref}\cdot{\bf n}=0$ on $\partial V$, then we can write

$H_{R}=\int_{V}{\bf A}^{*}\cdot{\bf B}\,\mathrm{d}V=\frac{1}{2}\int_{\partial
V}\mathcal{A}|B_{r}|\,\mathrm{d}S,$ (16)

where the second equality requires all field lines to intersect the boundary
so as to be counted in the integral (Berger, 1988). Note that if ${\bf
B}\equiv{\bf B}_{\rm ref}$ then the integrand in (15) vanishes everywhere,
while in (16) the integral vanishes only after integrating over all field
lines.

Figure 13.: Scatter plots between unsigned helicity, $\overline{H}$, and
(unsigned) relative helicity, $H_{R}$, for (a) individual time snapshots, and
(b) cycle-averaged values in each simulation run. The dashed line in (b) shows
the indicated linear fit, using all six points.

Comparing Equation (16) with Equation (11) shows that
$|H_{R}|\leq\overline{H}$, and this is clearly seen in Figure 13(a), which
plots $\overline{H}$ against $|H_{R}|$ for individual time snapshots. It is
very clear from this figure that the two measures are not correlated at
individual times, and also the fluctuation of $|H_{R}|$ makes its change from
Minimum to Maximum less meaningful. On the other hand, there is a reasonable
correlation between $|H_{R}|$ and $\overline{H}$ when averaged over all
snapshots during the cycle (Figure 13b), albeit with some scatter – for
example, run BMR0 has lower than expected $\overline{H}$ given $|H_{R}|$,
indicating a more uniform sign of $\mathcal{A}$ owing to the simpler surface
$B_{r}$ pattern. This correlation persists in Figure 10 if $\overline{H}$ is
replaced by $|H_{R}|$ (not illustrated), albeit weaker. Taking the linear fit
for $\overline{H}$ against $|H_{R}|$ from Figure 13(b) – with its significant
non-zero offset – and inserting this into the power law from Figure 10(a),
suggests an offset power law $|H_{R}|\,[\mathrm{Mx}^{2}]=3.58\times
10^{31}E_{\rm free}^{0.38}-0.17\times 10^{44}\,[\mathrm{erg}]$.

##  Appendix B Parameter dependence

The parameters in the magneto-frictional model are not directly constrained by
observations. Therefore, to check the robustness of our results, we have
carried out further 12-month runs, varying one coronal parameter at a time
compared to run T0. Each of these runs was initialised using the same three-
dimensional snapshot of the non-potential ${\bm{B}}$ from run T0 on 2013
December 29. Figures 14(a,b) show an analogous plot to Figure 10, but with the
quantities $\overline{H}$, $E_{\rm free}$ and $\langle\ddot{\Phi}_{1}\rangle$
averaged only over the 12 months of 2014. The parameters which have been
varied are indicated in the legend: we try both reducing and increasing by
$20\%$ the friction coefficient $\nu_{0}$, the hyperdiffusivity $\eta_{h}$,
and the outflow speed $v_{w}$. In all of these cases, the other parameters
keep their default values and the solar surface $B_{r}$ and emerging region
properties remain identical. The figure also shows runs T0, TU0.1 and TU0.5,
and notice that the average values are all increased compared to Figure 10
because the averages are taken over 2014 only (an active period), rather than
the full simulation. Note that the exponent of the power-law fit in Figure
14(a) is comparable to that for the whole cycle in Figure 10.

Figure 14.: Parameter dependence of the relations between averaged unsigned
helicity, $\overline{H}$ and (a) free energy, $E_{\rm free}$; (b) eruptivity,
$\langle\ddot{\Phi}_{1}\rangle$, computed over 2014 (see text). Panels (c) and
(d) show enlargements of the regions around run T0. The dashed line in (a)
shows a power-law fit between T0, TU0.05 and TU0.1 only, for context, while in
(b) the dashed line guides the eye between these three runs.

The most significant conclusion is that the $20\%$ parameter variations do not
substantially affect the mean quantities in Figure 14, leading to standard
deviations (relative to run T0) of only $1\%$ in $\overline{H}$, and $5\%$ in
$E_{\rm free}$ or $\langle\ddot{\Phi}_{1}\rangle$. This is much smaller than
the effect of twisting the emerging regions, as seen by comparing runs TU0.05
or TU0.1.

Looking more closely, Figures 14(c,d) show enlargements of the region around
T0 in Figures 14(a,b), respectively. Firstly, note that increasing/decreasing
the friction coefficient $\nu_{0}$ has the effect of increasing/decreasing
$E_{\rm free}$, because higher $\nu_{0}$ means slower relaxation. On the other
hand, increasing/decreasing $\eta_{h}$ reduces/increases $E_{\rm free}$,
because the hyperdiffusion acts to dissipate magnetic energy. Varying $v_{w}$
has only a weak effect on $E_{\rm free}$, because most of the energy is stored
at lower heights which are weakly affected by the outflow. On the other hand,
$v_{w}$ has a relatively stronger effect on the eruptivity in Figure 14(d),
with higher $v_{w}$ increasing $\langle\ddot{\Phi}_{1}\rangle$ and vice versa.
Increasing/decreasing $\eta_{h}$ reduces/increases
$\langle\ddot{\Phi}_{1}\rangle$, which is consistent with Yeates and Mackay
(2009a) who found that decreased ohmic diffusion leads to more flux rope
eruptions because ropes that form are larger and more twisted. Finally
increasing/decreasing $\nu_{0}$ leads to a small reduction/increase in
$\langle\ddot{\Phi}_{1}\rangle$ consistent with the change in $E_{\rm free}$.

###  B.1 Mesh resolution

There are two additional points in Figure 14 labelled $0.5n_{i}$. These
represent completely new simulation runs (from 12 June 2010) at half of the
original mesh resolution, namely $(30,90,180)$ cells, firstly with $\tau=0$
and secondly with uniform twist $|\tau|=0.1$. In these runs, the emerging
regions were re-extracted from HMI/SHARPs, and there are fewer of them because
some of the original regions fell below the resolution limit. Nevertheless,
Figure 14(a) shows that these much less computationally intensive runs were
still able to predict the 2014 averages of $E_{\rm free}$ to within $5\%$ for
both T0 and TU0.1, and of $\overline{H}$ to within $\pm 7\%$. On the other
hand, Figure 14(b) shows that the eruptivity $\langle\ddot{\Phi}_{1}\rangle$
is underestimated by $20\%$ for run T0 and almost $40\%$ for run TU0.1.

## References

  * Altschuler and Newkirk (1969) Altschuler, M.D., Newkirk, G.: 1969, Magnetic Fields and the Structure of the Solar Corona. I: Methods of Calculating Coronal Fields. Solar Phys. 9, 131\. DOI. ADS.
  * Berger (1988) Berger, M.A.: 1988, An energy formula for nonlinear force-free magnetic fields. Astron. Astrophys. 201, 355\. ADS.
  * Berger and Field (1984) Berger, M.A., Field, G.B.: 1984, The topological properties of magnetic helicity. J. Fluid Mech. 147, 133\. DOI. ADS.
  * Berger and Hornig (2018) Berger, M.A., Hornig, G.: 2018, A generalized poloidal-toroidal decomposition and an absolute measure of helicity. J. Phys. A Math. Gen. 51, 495501\. DOI. ADS.
  * Bhowmik and Yeates (2021) Bhowmik, P., Yeates, A.R.: 2021, Two Classes of Eruptive Events During Solar Minimum. Solar Phys. 296, 109\. DOI. ADS.
  * Bobra et al. (2014) Bobra, M.G., Sun, X., Hoeksema, J.T., Turmon, M., Liu, Y., Hayashi, K., Barnes, G., Leka, K.D.: 2014, The Helioseismic and Magnetic Imager (HMI) Vector Magnetic Field Pipeline: SHARPs - Space-Weather HMI Active Region Patches. Solar Phys. 289, 3549\. DOI. ADS.
  * Boozer (1986) Boozer, A.H.: 1986, Ohm’s law for mean magnetic fields. J. Plasma Phys. 35, 133\. DOI. ADS.
  * Chen, Rempel, and Fan (2022) Chen, F., Rempel, M., Fan, Y.: 2022, A Comprehensive Radiative Magnetohydrodynamics Simulation of Active Region Scale Flux Emergence from the Convection Zone to the Corona. Astrophys. J. 937, 91\. DOI. ADS.
  * Cheung and Isobe (2014) Cheung, M.C.M., Isobe, H.: 2014, Flux Emergence (Theory). Living Rev. Solar Phys. 11, 3\. DOI. ADS.
  * Chifu, Inhester, and Wiegelmann (2022) Chifu, I., Inhester, B., Wiegelmann, T.: 2022, Coronal magnetic field evolution over cycle 24. Astron. Astrophys. 659, A174. DOI. ADS.
  * Chodura and Schlüter (1981) Chodura, R., Schlüter, A.: 1981, A 3D Code for MHD Equilibrium and Stability. J. Comp. Phys. 41, 68\. DOI. ADS.
  * DeRosa et al. (2015) DeRosa, M.L., Wheatland, M.S., Leka, K.D., Barnes, G., Amari, T., Canou, A., Gilchrist, S.A., Thalmann, J.K., Valori, G., Wiegelmann, T., Schrijver, C.J., Malanushenko, A., Sun, X., Régnier, S.: 2015, The Influence of Spatial resolution on Nonlinear Force-free Modeling. Astrophys. J. 811, 107\. DOI. ADS.
  * Edwards et al. (2015) Edwards, S.J., Yeates, A.R., Bocquet, F.-X., Mackay, D.H.: 2015, Influence of Non-Potential Coronal Magnetic Topology on Solar-Wind Models. Solar Phys. 290, 2791\. DOI. ADS.
  * Finn and Antonsen (1985) Finn, J.M., Antonsen, J. Thomas M.: 1985, Magnetic helicity: What is it and what is it good for? Comm. Plasma Phys. Controlled Fusion 9, 111\. ADS.
  * Georgoulis et al. (2009) Georgoulis, M.K., Rust, D.M., Pevtsov, A.A., Bernasconi, P.N., Kuzanyan, K.M.: 2009, Solar Magnetic Helicity Injected into the Heliosphere: Magnitude, Balance, and Periodicities Over Solar Cycle 23. Astrophys. J. Lett. 705, L48. DOI. ADS.
  * Hathaway (2015) Hathaway, D.H.: 2015, The Solar Cycle. Living Rev. Solar Phys. 12, 4\. DOI. ADS.
  * Hayashi et al. (2021) Hayashi, K., Abbett, W.P., Cheung, M.C.M., Fisher, G.H.: 2021, Coupling a Global Heliospheric Magnetohydrodynamic Model to a Magnetofrictional Model of the Low Corona. Astrophys. J. Supp. Ser. 254, 1\. DOI. ADS.
  * Hoeksema et al. (2020) Hoeksema, J.T., Abbett, W.P., Bercik, D.J., Cheung, M.C.M., DeRosa, M.L., Fisher, G.H., Hayashi, K., Kazachenko, M.D., Liu, Y., Lumme, E., Lynch, B.J., Sun, X., Welsch, B.T.: 2020, The Coronal Global Evolutionary Model: Using HMI Vector Magnetogram and Doppler Data to Determine Coronal Magnetic Field Evolution. Astrophys. J. Supp. Ser. 250, 28\. DOI. ADS.
  * Lamy et al. (2019) Lamy, P.L., Floyd, O., Boclet, B., Wojak, J., Gilardy, H., Barlyaeva, T.: 2019, Coronal Mass Ejections over Solar Cycles 23 and 24. Space Sci. Rev. 215, 39\. DOI. ADS.
  * Lehmann et al. (2017) Lehmann, L.T., Jardine, M.M., Vidotto, A.A., Mackay, D.H., See, V., Donati, J.-F., Folsom, C.P., Jeffers, S.V., Marsden, S.C., Morin, J., Petit, P.: 2017, The energy budget of stellar magnetic fields: comparing non-potential simulations and observations. Mon. Not. R. Astron. Soc. 466, L24. DOI. ADS.
  * Linker et al. (2017) Linker, J.A., Caplan, R.M., Downs, C., Riley, P., Mikic, Z., Lionello, R., Henney, C.J., Arge, C.N., Liu, Y., Derosa, M.L., Yeates, A., Owens, M.J.: 2017, The Open Flux Problem. Astrophys. J. 848, 70\. DOI. ADS.
  * Liu, Hoeksema, and Sun (2014) Liu, Y., Hoeksema, J.T., Sun, X.: 2014, Test of the Hemispheric Rule of Magnetic Helicity in the Sun Using the Helioseismic and Magnetic Imager (HMI) Data. Astrophys. J. Lett. 783, L1. DOI. ADS.
  * Lowder and Yeates (2017) Lowder, C., Yeates, A.: 2017, Magnetic Flux Rope Identification and Characterization from Observationally Driven Solar Coronal Models. Astrophys. J. 846, 106\. DOI. ADS.
  * Mackay and Upton (2022) Mackay, D.H., Upton, L.A.: 2022, A Comparison of Global Magnetofrictional Simulations of the 2015 March 20 Solar Eclipse. Astrophys. J. 939, 9\. DOI. ADS.
  * Mackay and van Ballegooijen (2006) Mackay, D.H., van Ballegooijen, A.A.: 2006, Models of the Large-Scale Corona. I. Formation, Evolution, and Liftoff of Magnetic Flux Ropes. Astrophys. J. 641, 577\. DOI. ADS.
  * Mackay, DeVore, and Antiochos (2014) Mackay, D.H., DeVore, C.R., Antiochos, S.K.: 2014, Global-scale Consequences of Magnetic-helicity Injection and Condensation on the Sun. Astrophys. J. 784, 164\. DOI. ADS.
  * Mackay, Gaizauskas, and van Ballegooijen (2000) Mackay, D.H., Gaizauskas, V., van Ballegooijen, A.A.: 2000,. Astrophys. J. 544, 1122\. DOI. ADS.
  * Mackay et al. (2010) Mackay, D.H., Karpen, J.T., Ballester, J.L., Schmieder, B., Aulanier, G.: 2010, Physics of Solar Prominences: II—Magnetic Structure and Dynamics. Space Sci. Rev. 151, 333\. DOI. ADS.
  * Mackay et al. (2018) Mackay, D.H., DeVore, C.R., Antiochos, S.K., Yeates, A.R.: 2018, Magnetic Helicity Condensation and the Solar Cycle. Astrophys. J. 869, 62\. DOI. ADS.
  * Malanushenko et al. (2014) Malanushenko, A., Schrijver, C.J., DeRosa, M.L., Wheatland, M.S.: 2014, Using Coronal Loops to Reconstruct the Magnetic Field of an Active Region before and after a Major Flare. Astrophys. J. 783, 102\. DOI. ADS.
  * Martin (1998) Martin, S.F.: 1998, Conditions for the Formation and Maintenance of Filaments (Invited Review). Solar Phys. 182, 107\. DOI. ADS.
  * Meyer et al. (2020) Meyer, K.A., Mackay, D.H., Talpeanu, D.-C., Upton, L.A., West, M.J.: 2020, Investigation of the Middle Corona with SWAP and a Data-Driven Non-Potential Coronal Magnetic Field Model. Solar Phys. 295, 101\. DOI. ADS.
  * Mikić and Linker (1994) Mikić, Z., Linker, J.A.: 1994, Disruption of Coronal Magnetic Field Arcades. Astrophys. J. 430, 898\. DOI. ADS.
  * Moraitis, Patsourakos, and Nindos (2021) Moraitis, K., Patsourakos, S., Nindos, A.: 2021, Relative field line helicity of a large eruptive solar active region. Astron. Astrophys. 649, A107. DOI. ADS.
  * Pevtsov (2002) Pevtsov, A.A.: 2002, Active-Region Filaments and X-ray Sigmoids. Solar Phys. 207, 111\. DOI. ADS.
  * Pevtsov and Balasubramaniam (2003) Pevtsov, A.A., Balasubramaniam, K.S.: 2003, Helicity patterns on the sun. Adv. Space Res. 32, 1867\. DOI. ADS.
  * Régnier and Priest (2007) Régnier, S., Priest, E.R.: 2007, Free Magnetic Energy in Solar Active Regions above the Minimum-Energy Relaxed State. Astrophys. J. Lett. 669, L53. DOI. ADS.
  * Rice and Yeates (2021) Rice, O.E.K., Yeates, A.R.: 2021, Global Coronal Equilibria with Solar Wind Outflow. Astrophys. J. 923, 57\. DOI. ADS.
  * Rumińska et al. (2022) Rumińska, A., Ba̧k-Stȩślicka, U., Gibson, S.E., Fan, Y.: 2022, Coronal Cavities in CoMP Observations. Astrophys. J. 926, 146\. DOI. ADS.
  * Russell et al. (2019) Russell, A.J.B., Demoulin, P., Hornig, G., Pontin, D.I., Candelaresi, S.: 2019, Do Current and Magnetic Helicities Have the Same Sign? Astrophys. J. 884, 55\. DOI. ADS.
  * Schatten, Wilcox, and Ness (1969) Schatten, K.H., Wilcox, J.M., Ness, N.F.: 1969, A model of interplanetary and coronal magnetic fields. Solar Phys. 6, 442\. DOI. ADS.
  * Schrijver et al. (2013) Schrijver, C.J., Title, A.M., Yeates, A.R., DeRosa, M.L.: 2013, Pathways of Large-scale Magnetic Couplings between Solar Coronal Events. Astrophys. J. 773, 93\. DOI. ADS.
  * Sheeley et al. (2013) Sheeley, J. N. R., Martin, S.F., Panasenco, O., Warren, H.P.: 2013, Using Coronal Cells to Infer the Magnetic Field Structure and Chirality of Filament Channels. Astrophys. J. 772, 88\. DOI. ADS.
  * Snodgrass and Ulrich (1990) Snodgrass, H.B., Ulrich, R.K.: 1990, Rotation of Doppler Features in the Solar Photosphere. Astrophys. J. 351, 309\. DOI. ADS.
  * Stansby, Yeates, and Badman (2020) Stansby, D., Yeates, A., Badman, S.: 2020, pfsspy: A Python package for potential field source surface modelling. J. Open Source Software 5, 2732\. DOI. ADS.
  * Sullivan and Kaszynski (2019) Sullivan, B., Kaszynski, A.: 2019, PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK). Journal of Open Source Software 4, 1450\. DOI. https://doi.org/10.21105/joss.01450.
  * Sun (2018) Sun, X.: 2018, Polar Field Correction for HMI Line-of-Sight Synoptic Data. arXiv e-prints, arXiv:1801.04265. DOI. ADS.
  * Sun et al. (2015) Sun, X., Bobra, M.G., Hoeksema, J.T., Liu, Y., Li, Y., Shen, C., Couvidat, S., Norton, A.A., Fisher, G.H.: 2015, Why Is the Great Solar Active Region 12192 Flare-rich but CME-poor? Astrophys. J. Lett. 804, L28. DOI. ADS.
  * Tadesse et al. (2014) Tadesse, T., Pevtsov, A.A., Wiegelmann, T., MacNeice, P.J., Gosain, S.: 2014, Global Solar Free Magnetic Energy and Electric Current Density Distribution of Carrington Rotation 2124. Solar Phys. 289, 4031\. DOI. ADS.
  * Thalmann, Gupta, and Veronig (2022) Thalmann, J.K., Gupta, M., Veronig, A.M.: 2022, The effect of spatial sampling on magnetic field modeling and helicity computation. Astron. Astrophys. 662, A3. DOI. ADS.
  * Tziotziou, Georgoulis, and Raouafi (2012) Tziotziou, K., Georgoulis, M.K., Raouafi, N.-E.: 2012, The Magnetic Energy-Helicity Diagram of Solar Active Regions. Astrophys. J. Lett. 759, L4. DOI. ADS.
  * van Ballegooijen and Cranmer (2008) van Ballegooijen, A.A., Cranmer, S.R.: 2008, Hyperdiffusion as a Mechanism for Solar Coronal Heating. Astrophys. J. 682, 644\. DOI. ADS.
  * van Ballegooijen and Mackay (2007) van Ballegooijen, A.A., Mackay, D.H.: 2007, Model for the Coupled Evolution of Subsurface and Coronal Magnetic Fields in Solar Active Regions. Astrophys. J. 659, 1713\. DOI. ADS.
  * van Ballegooijen, Priest, and Mackay (2000) van Ballegooijen, A.A., Priest, E.R., Mackay, D.H.: 2000, Mean Field Model for the Formation of Filament Channels on the Sun. Astrophys. J. 539, 983\. DOI. ADS.
  * Wagner et al. (2022) Wagner, A., Asvestari, E., Temmer, M., Heinemann, S.G., Pomoell, J.: 2022, Validation scheme for solar coronal models: Constraints from multi-perspective observations in EUV and white light. Astron. Astrophys. 657, A117. DOI. ADS.
  * Webb and Howard (2012) Webb, D.F., Howard, T.A.: 2012, Coronal Mass Ejections: Observations. Living Rev. Solar Phys. 9, 3\. DOI. ADS.
  * Wiegelmann and Sakurai (2021) Wiegelmann, T., Sakurai, T.: 2021, Solar force-free magnetic fields. Living Rev. Solar Phys. 18, 1\. DOI. ADS.
  * Xiao, Prior, and Yeates (2023) Xiao, D., Prior, C.B., Yeates, A.R.: 2023, Spherical winding and helicity. J. Phys. A Math. Gen. 56, 205201\. DOI. ADS.
  * Yang, Sturrock, and Antiochos (1986) Yang, W.H., Sturrock, P.A., Antiochos, S.K.: 1986, Force-free Magnetic Fields: The Magneto-frictional Method. Astrophys. J. 309, 383\. DOI. ADS.
  * Yardley, Mackay, and Green (2021) Yardley, S.L., Mackay, D.H., Green, L.M.: 2021, Simulating the Coronal Evolution of Bipolar Active Regions to Investigate the Formation of Flux Ropes. Solar Phys. 296, 10\. DOI. ADS.
  * Yardley et al. (2021) Yardley, S.L., Pagano, P., Mackay, D.H., Upton, L.A.: 2021, Determining the source and eruption dynamics of a stealth CME using NLFFF modelling and MHD simulations. Astron. Astrophys. 652, A160. DOI. ADS.
  * Yeates (2018) Yeates, A.: 2018, antyeates1983/pfss: First release of PFSS Code, v1.0, Zenodo. DOI. https://zenodo.org/record/1472183.
  * Yeates (2020a) Yeates, A.R.: 2020a, Bipolar magnetic regions determined from HMI SHARPs data. DOI. https://doi.org/10.7910/DVN/1Z7YMT.
  * Yeates (2020b) Yeates, A.R.: 2020b, How Good Is the Bipolar Approximation of Active Regions for Surface Flux Transport? Solar Phys. 295, 119\. DOI. ADS.
  * Yeates (2020c) Yeates, A.R.: 2020c, The Minimal Helicity of Solar Coronal Magnetic Fields. Astrophys. J. Lett. 898, L49. DOI. ADS.
  * Yeates and Bhowmik (2022) Yeates, A.R., Bhowmik, P.: 2022, Automated Driving for Global Nonpotential Simulations of the Solar Corona. Astrophys. J. 935, 13\. DOI. ADS.
  * Yeates and Hornig (2016) Yeates, A.R., Hornig, G.: 2016, The global distribution of magnetic helicity in the solar corona. Astron. Astrophys. 594, A98. DOI. ADS.
  * Yeates and Mackay (2009a) Yeates, A.R., Mackay, D.H.: 2009a, Initiation of Coronal Mass Ejections in a Global Evolution Model. Astrophys. J. 699, 1024\. DOI. ADS.
  * Yeates and Mackay (2009b) Yeates, A.R., Mackay, D.H.: 2009b, Modelling the Global Solar Corona: III. Origin of the Hemispheric Pattern of Filaments. Solar Phys. 254, 77\. DOI. ADS.
  * Yeates, Constable, and Martens (2010) Yeates, A.R., Constable, J.A., Martens, P.C.H.: 2010, Solar Cycle Variation of Magnetic Flux Ropes in a Quasi-Static Coronal Evolution Model. Solar Phys. 263, 121\. DOI. ADS.
  * Yeates, Mackay, and van Ballegooijen (2008) Yeates, A.R., Mackay, D.H., van Ballegooijen, A.A.: 2008, Modelling the Global Solar Corona II: Coronal Evolution and Filament Chirality Comparison. Solar Phys. 247, 103\. DOI. ADS.
  * Yeates et al. (2010) Yeates, A.R., Mackay, D.H., van Ballegooijen, A.A., Constable, J.A.: 2010, A nonpotential model for the Sun’s open magnetic flux. J. Geophys. Res. (Space Phys.) 115, A09112. DOI. ADS.
  * Yeates et al. (2018) Yeates, A.R., Amari, T., Contopoulos, I., Feng, X., Mackay, D.H., Mikić, Z., Wiegelmann, T., Hutton, J., Lowder, C.A., Morgan, H., Petrie, G., Rachmeler, L.A., Upton, L.A., Canou, A., Chopin, P., Downs, C., Druckmüller, M., Linker, J.A., Seaton, D.B., Török, T.: 2018, Global Non-Potential Magnetic Models of the Solar Corona During the March 2015 Eclipse. Space Sci. Rev. 214, 99\. DOI. ADS.