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

The global magnetic field in the solar corona is known to contain free magnetic energy and magnetic helicity above that of a current-free (potential) state. But the strength of this non-potentiality and its evolution over the solar cycle remain uncertain. Here we model the corona over Solar Cycle 24 using a simplified magneto-frictional model that retains the magnetohydrodynamic induction equation but assumes relaxation towards force-free equilibrium, driven by solar surface motions and flux emergence. The model is relatively conservative compared to some others in the literature, with free energy approximately 20 − 25% of the potential field energy. We find that unsigned helicity is about a factor 10 higher at Maximum than Minimum, while free magnetic energy shows an even greater increase. The cycle averages of these two quantities are linearly correlated, extending a result found previously for active regions. Also, we propose a practical measure of eruptivity for these simulations, and show that this increases concurrently with the sunspot number, in accordance with observed coronal mass ejection rates. Whilst shearing by surface motions generates 50% or more of the free energy and helicity in the corona, we show that active regions must emerge with their own internal helicity otherwise the eruptivity is substantially reduced and follows the wrong pattern over time.


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 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 × 10 33 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 magnetofrictional 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).

Model setup
We simulate the corona using a magneto-frictional model for the mean (largescale) magnetic field 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.

Magneto-frictional model
Our magneto-frictional model for the mean/large-scale magnetic field is based on the following main assumptions.
i) The mean magnetic field 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.ii) The high magnetic Reynolds number implies a near-ideal evolution of B, so it is important to retain the mean-field induction equation where the non-ideal term N is small.This means in particular that B holds a memory of previous interactions, allowing topological structure to build over time.iii) The high Alfvén speed in the corona makes 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 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 J × B = 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).iv) In the outer corona, the solar wind starts to influence B, preventing it from being force-free.This is accounted for by the radial outflow v 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 ⊙ and fix the outflow speed v w = 100 km s −1 in our reference model.v) The coronal B relaxes with respect to a rest frame rotating with the Sun.
Our simulations use the Carrington frame.
The non-ideal electric field, 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 which prevents volume dissipation of relative helicity (with respect to a potential field) since the corresponding volume dissipation term V E • B dV (Berger and Field, 1984) becomes a boundary integral (Boozer, 1986;van Ballegooijen and Cranmer, 2008).We set the reference value η h = 10 11 km 4 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 N , although further observational comparison is needed to be sure of the most realistic approach.
The "friction" coefficient ν 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 ν = ν 0 (r cos λ) −2 , where λ denotes latitude.In our reference model the amplitude is set to ν 0 = 2.8 × 10 5 s.
In our numerical implementation, we write B = ∇×A and solve the uncurled version of Equation (1) for the vector potential 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 ⊙ ), 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.

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 ⊙ , assuming J = 0 in the corona with B × r = 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 e −cwl(l+1) , with c w = 5 × 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 J = 0 to machine precision for our particular discretisation scheme.

Boundary conditions
At the inner boundary the horizontal electric field from Equation ( 1) is replaced by The first term, E 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 smallscale surface motions).The large-scale flow velocity is assumed to take the steady form v s = v θ (θ) θ + R ⊙ sin θ Ω(θ) φ, where v θ denotes the meridional (poleward) flow and Ω(θ) 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 v s and η 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 θ (θ) = D u cos θ sin p0 θ, with p 0 = 2.33 and D u = 0.041 km s −1 (giving a peak flow speed of 15 m s −1 ), and η 0 = 350 km 2 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 Ω(θ) = 0.18 − 2.396 cos 2 θ − 1.787 cos 4 θ (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 km s −1 (giving a peak flow speed of 20 m s −1 ), to SOLA: ms.tex; 24 May 2024; 2:44; p. 6 There is a hemispheric tendency in α 0 with ⟨α 0 ⟩ λ 0 >0 = −0.0147Mm −1 and ⟨α 0 ⟩ λ 0 <0 = 0.00510 Mm −1 .The percentage of regions obeying the hemispheric rule is 75.9% in the North and 60.6% in the South.
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 J × r = 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 N on these boundaries, we impose zerogradient conditions on the radial component of the hyperdiffusive flux B 2 ∇(J • B/B 2 ).

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 timelatitude distribution of regions, coloured by their unsigned magnetic flux, Φ 0 (including both polarities).The distribution of region fluxes is shown in Figure 1(b).
Figure 1(c) shows the same SHARP regions but coloured instead by the observed "twist" parameter α 0 = J z B z /( B 2 z ) 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 α 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 α 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 α 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 em = 24 hr we impose a steady (surface) electric field The first term E 0 ⊥ 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 E 0 ⊥ = 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 Here B r is a smoothed magnetogram used to identify the polarity inversion lines, and ⟨f pil ⟩ 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 ≈ |τ B r |.The single dimensionless "twist" parameter τ then controls the helicity injection within the emerging region, with τ > 0 leading to positive helicity and τ < 0 to negative helicity.Figure 2(b) shows the same region as Figure 2(a), but emerged with τ = −0.1 giving it negative helicity -a reverse-S shaped twist.
Since the optimum value of τ to model each active region is unknown, we follow Yeates and Bhowmik (2022) and vary τ in a parameter study.The simulation runs are summarised in Table 1.In TU0.05 and TU0.1, we simply fix a uniform τ 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 |τ | = 0.1 is chosen because larger SOLA: ms.tex; 24 May 2024; 2:44; p. 8  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 τ for each emerging region, chosen proportional to the observed twist values α 0 .Since the dimensions of τ and α 0 are different, the proportionality constant ℓ in Table 1 has dimensions of length.
The values of ℓ are chosen so that the two runs have a comparable overall helicity to the two runs with uniform τ .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).

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), SOLA: ms.tex; 24 May 2024; 2:44; p. 9 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 Φ 0 for individual regions, it does slightly increase Φ 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(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) Φ 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).

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.

Energy
Figure 5(a) shows total magnetic energy, along with the energy of a corresponding PFSS extrapolation.The cyan curve shows free energy, where B p is the current-free magnetic field in the coronal volume matching B pr = B r on both r = R ⊙ 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 free /E 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 free /E p was in the range 40 − 50%.However, we will see in Section 4 that E free is substantially increased if the active regions emerge twisted.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 p show a similar relative increase, as does E free albeit a little higher.However, one has to bear in mind that these are snapshots rather than running means; in particular, since E free tends to lag E 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 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.

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.
Mean current density -shown in Figure 6(a) is defined as 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 where A is the field line helicity for a magnetic field line L traced from a given point on ∂V , and A * is an appropriate vector potential.Figure 6(c) shows a similar rate of increase in H from Minimum to Maximum compared to ⟨J⟩, albeit a rather different time profile, relatively lower in the increasing phase of the cycle.Like magnetic energy, 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 H rather than the simpler-to-compute V |A•B| dV is that 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 ⟨J⟩).Thus it is a more robust measure than V |A•B| dV .Our rationale for using 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 A(L) and hence of H depend on the chosen gauge of 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 A * = ∇ × (P r) + T r for suitable P and T .Since A * differs from the gauge of A used to solve Equation (1), we compute A * numerically to find A and 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 H and free energy.
Note that even potential fields can have non-zero 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 H in the PFSS model is approximately 33% of 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 H in run T0.However, H also fluctuates more rapidly in the magneto-frictional simulation due to sudden reconfigurations such as flux rope ejections.

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, As illustrated by Bhowmik and Yeates (2021), eruptions of individual structures in the model lead to temporary enhancements in Φ 1 .Taking the second time derivative, Φ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).
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 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 Φ1 .
For clearer comparison over the full solar cycle, we take our eruptivity proxy to be the 27-day average, ⟨ Φ1 ⟩.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 × 10 24 erg s −1 , but varying this threshold changes only the overall normalisation in Figure 7(g) and not the relative variation over the cycle.Thus ⟨ Φ1 ⟩ 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 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 ⟨ Φ1 ⟩ 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 Φ 0 , ⟨J⟩, or H, but rather less than E 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, ⟨ Φ1 ⟩ 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.

Effect of emerging region properties
We have repeated the analysis of coronal properties for each simulation run in Table 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 τ .
Figure 8 compares the time evolution of E free and H for the different runs.Both quantities are enhanced in all of the runs with τ ̸ = 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 free can be compared with the NLFFF extrapolations in Figure 8 of Chifu, Inhester, and Wiegelmann (2022).Even in run TOb10, our E 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 free peaks at about 2.5 × 10 33 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 free shows sharper fluctuations, during which the magnitude can be comparable to Chifu, Inhester, and Wiegelmann (2022).On the other hand, our E free values around March 2015 look broadly consistent with those obtained by Mackay and Upton (2022) in their magneto-frictional runs with hyperdiffusion.
Figures 9(b,c,e,f) show the normalised cycle variation of E free and H for each run, as well as the eruptivity ⟨ Φ1 ⟩.Again we see that, as we increase |τ | -either the uniform value in runs TU0.05 and TU0.1 or the factor multiplying α 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 free , which now shows a greater increase than the sunspot number.The pattern of variation in ⟨ Φ1 ⟩ 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.

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 H.In Figure 10(a), we see a clear power-law relation between the cycle-averages of E free and H, which holds for all of the runs (including BMR0).It is also clear how both quantities increase as the emerging region twist τ is increased.The intermediate runs TU0.05 and TOb5 are roughly comparable to one another in both H and E free .Compared to the reference run T0, these two runs have approximately 20% more unsigned helicity and 55% more free energy.
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 × 10 14 E 0.897 free , where E free is in ergs and H R is the usual relative magnetic helicity in Mx 2 .If this is extrapolated to a global value of H R = 1.5 × 10 44 Mx 2 -approximately the cycle mean of H for run TU0.05 -it predicts E free ≈ 3 × 10 33 erg, which is an order of magnitude greater than our simulations.One possible explanation for the difference is the fact that |H R | ≤ H, with equality only if 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 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 SOLA: ms.tex; 24 May 2024; 2:44; p. 18  value for run T0 seems to act as an approximate lower bound on the free energy, at a level of unsigned helicity.Figure 10(b) shows that the eruptivity seems to depend not only on 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 τ 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 |τ | (of both signs), with a sizeable tail of regions having |τ | > 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 |τ | = 0.1.The most SOLA: ms.tex;24 May 2024;2:44;p. 19 twisted of these regions are likely to erupt rapidly, before the magnetic field has a chance to relax or dissipate strong currents.

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 H and E free are approximately 8% and 18% higher than run T0, and this increase is consistent with the increase in surface flux Φ 0 seen in Figure 3(b).In Figure 10(b), we see that despite the increase H the eruptivity ⟨ Φ1 ⟩ 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.

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 free /E 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 N can have a substantial effect on the free energy, with fully ideal simulations (N = 0) having E free /E p in excess of 60%.Compared to this, our simulations are rather conservative in the amount of free energy (the cycle-average ratio E free /E 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 Φ 0 and E, though not always in E free /E 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 timing of specific eruptions (Yardley, Mackay, and Green, 2021).We hope to pursue this in future.(a) where 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 B ≡ B ref then the integrand in (15) vanishes everywhere, while in ( 16) the integral vanishes only after integrating over all field lines.Comparing Equation ( 16) with Equation (11) shows that |H R | ≤ H, and this is clearly seen in Figure 13

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 B from run T0 on 2013 December 29.Figures 14(a,b) show an analogous plot to Figure 10, but with the quantities H, E free and ⟨ Φ1 ⟩ 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 ν 0 , the hyperdiffusivity η 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.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 H, and 5% in E free or ⟨ Φ1 ⟩.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 ν 0 has the effect of increasing/decreasing E free , because higher ν 0 means slower relaxation.On the other hand, increasing/decreasing η h reduces/increases E free , because the hyperdiffusion acts to dissipate magnetic energy.Varying v w has only a weak effect on E 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 ⟨ Φ1 ⟩ and vice versa.Increasing/decreasing η h reduces/increases ⟨ Φ1 ⟩, 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 ν 0 leads to a small reduction/increase in ⟨ Φ1 ⟩ consistent with the change in E 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 τ = 0 and secondly with uniform twist |τ | = 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 free to within 5% for both T0 and TU0.1, and of H to within ±7%.On the other hand, Figure 14(b) shows that the eruptivity ⟨ Φ1 ⟩ is underestimated by 20% for run T0 and almost 40% for run TU0.1.

Figure 1 .
Figure 1.Emerging regions determined from the HMI/SHARPs data, shown in a time-latitude "butterfly" plot and coloured by either (a) unsigned flux, Φ 0 or (c) observed twist, α 0 .Panel (b) shows a histogram of the Φ 0 values, while (d) shows the mean and standard deviation of the α 0 values in equal sine-latitude bins.There is a hemispheric tendency in α 0 with ⟨α 0 ⟩ λ 0 >0 = −0.0147Mm −1 and ⟨α 0 ⟩ λ 0 <0 = 0.00510 Mm −1 .The percentage of regions obeying the hemispheric rule is 75.9% in the North and 60.6% in the South.

Figure 2 .
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 τ = 0 in (a,c) and τ = −0.1 in (b).Grayscale pixels show Br on r = R ⊙ (white positive/black negative, saturated at ±25 G), while magnetic field lines are coloured by field line helicity (red positive/blue negative, saturated at ±10 22 Mx).

Figure 3 .
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 4 .
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 ±10 22 Mx.Gray shading shows Br on the solar surface, black negative/white positive, saturated at ±15 G.

Figure 5 .
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, Ep, and cyan the free energy, E free .Panel (a) shows the original values while (b) shows the relative free energy, E free /Ep.Panel (c) shows the same curves as (a) but normalised as in Figure 3 and compared to the observed sunspot number.

Figure 5 Figure 6 .
Figure 6.Cycle variation of (a) the mean current density ⟨J⟩ and (b) the unsigned helicity 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 H for corresponding PFSS extrapolations from the same surface Br distribution.

Figure 7 .
Figure 7. Time series used to compute the eruptivity proxy, ⟨ Φ1 ⟩, shown for run T0.Panel (a) shows the open magnetic flux, Φ 1 over the full simulation, with enlargement of a selected period in (b).Panels (c-d) show the second time derivative, Φ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 ⟨ Φ1 ⟩, overlayed with the number of S 1 peaks in 27-day bins, both normalised to their maximum values.Panel (h) shows ⟨ Φ1 ⟩ normalised as in Figure 3 and compared to the observed sunspot number.

Figure 8 .
Figure 8.Time evolution of (a) E free and (b) 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).

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

Figure 10 .
Figure 10.Relation between cycle-averaged unsigned helicity, H and (a) free energy, E free ; or (b) eruptivity, ⟨ Φ1 ⟩.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 τ (red) and τ ∝ α 0 (blue).

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

Figure 13 .
Figure 13.Scatter plots between unsigned helicity, 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.
if H is replaced by |H R | (not illustrated), albeit weaker.Taking the linear fit for 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 | [Mx 2 ] = 3.58 × 10 31 E 0.38 free − 0.17 × 10 44 [erg].

Figure 14 .
Figure 14.Parameter dependence of the relations between averaged unsigned helicity, H and (a) free energy, E free ; (b) eruptivity, ⟨ Φ1 ⟩, 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.
Finn and Antonsen (1985)ed) relative helicity, H R and free energy, E free , for individual time snapshots with no cycle averaging.As in Figure11, 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.There is greatly increased scatter compared to Figure11, and, in particular, E free does not go to zero as |H R | goes to zero, unlike the behaviour with 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.Relative helicity is most commonly computed using theFinn and Antonsen (1985)formula, H.
the reference field.Since H R is gauge invariant for both A and A ref , if we use the canonical gauge A = A * from Section 3.2 and choose the gauge of A ref such that A * × A ref • n = 0 on ∂V , then we can write SOLA: ms.24 May 2024;024; 2:44; p. 23