Eruptivity Criteria for Two-Dimensional Magnetic Flux Ropes in the Solar Corona

1 Introduction

Flux ropes are twisted bundles of magnetic flux in the solar corona (Liu, 2020). Accurately predicting their behaviour is essential for reliable space weather predictions, as unstable flux ropes can erupt and lead to large coronal mass ejections (Forbes et al., 2006; Chen, 2011). The causes of such eruptions are not yet indisputably understood and a variety of mechanisms have been proposed. These have been explored with approaches ranging from analytic two-dimensional models to three-dimensional full magnetohydrodynamic (MHD) simulations.

Analytic two-dimensional models of flux rope behaviour date back to Kuperus and Raadu (1974) and van Tend and Kuperus (1978), who modelled a horizontal line current and its interaction with a specified background magnetic field. They established conditions on the current that allow it to be stable, and showed that, for an eruption to occur, the background magnetic field strength must decrease rapidly with increasing height.

Further analytic approaches have introduced the torus instability (e.g., Kliem and Török, 2006) where the flux rope is modelled as a current ring rather than a line current. Rather than a conditionon the decay of the background field with altitude, it is instead proposed that for instability to occur the background field component orthogonal to the torus must decrease sufficiently quickly. It has since been shown that the conditions for such an instability are essentially the same as for instability of a line current, and they are just two special cases of a continuous theory of more general current paths (Kliem et al., 2014). Although these works provide useful theoretical background as to the nature of flux rope behaviour, they are difficult to relate to physical ropes as these tend to be significantly more complex and much of the physics is by necessity not taken into account. There are indeed other conflicting explanations for instability, such as that of Ishiguro and Kusano (2017), who propose that fast magnetic reconnection can occur below a flux rope even if the overlying magnetic field does not decay with altitude.

Full 3D MHD simulations of flux rope formation and behaviour are also now possible. Leake et al. (2013), Leake et al. (2014) presented a series of such simulations, based on magnetic flux emergence into a preexisting field. Some ropes were found to be eruptive and others stable, depending on the configuration of the background magnetic field. Pariat et al. (2017) then analysed various properties of the system for these simulations, looking for properties that could predict whether or not the flux rope would erupt. Most diagnostics, such as the magnetic free energy and relative magnetic helicity, were shown not to have any strong correlation to eruptivity. However, they identified an “eruptivity index”—the fraction of the relative magnetic helicity that is in the current-carrying component—which became large only in erupting cases. In a different set of numerical simulations, where coronal flux ropes are created by a variety of solar surface footpoint motions, this eruptivity index again increased prior to eruption and the eruption was found to occur for a fixed value of the index, irrespective of the pattern of footpoint motions that injected the energy (Zuccarello et al., 2018). However, while promising, a complete theoretical understanding of this eruptivity index and its generality remains lacking.

Our modeling approach is chosen as a compromise between the simple analytic models of the past and the expensive full MHD simulations now possible. We use the magneto-frictional model, pioneered by Yang et al. (1986), whereby a fictional velocity field is determined explicitly from the magnetic field as opposed to using the fluid equations. For tracking the quasi-static injection of magnetic energy into the corona through solar surface motions, the model provides a viable alternative to full MHD simulations, at a fraction of the computational cost. Magneto-frictional models have previously been coupled with time-dependent lower boundary conditions based on imposed surface flux transport (e.g., Yeates et al., 2008) and realistic flux rope formation, as well as loss of equilibrium, is observed in such simulations (Mackay and van Ballegooijen, 2006; Yeates and Mackay, 2009; Lowder and Yeates, 2017; Hoeksema et al., 2020). Flux ropes are also formed on a smaller scale when the model is driven by high-resolution observations within active regions (e.g., Yardley et al., 2018). Since the magneto-frictional model neglects the full equation of motion, it cannot accurately follow the dynamics once a flux rope loses equilibrium and erupts. Nevertheless, if ropes produced by magneto-friction are used to initialise full MHD simulations, it is found that the magnetic flux ropes do indeed lose equilibrium at the same point, and can lead to CME eruptions (Kliem et al., 2013; Pagano et al., 2013). Thus magneto-friction can act as an accurate model for the pre-eruption evolution.

In our study, we use the magneto-frictional approach but simplify it to a 2.5-dimensional (translationally invariant) cartesian domain, which significantly increases computational speed while still exhibiting the fundamental features of fully 3D models. In particular, both the flux rope and the (sheared) overlying arcade are formed self-consistently by flux cancellation and shearing motions on the photosphere. In our model these shearing motions are purely large-scale, essentially modelling the differential rotation. As such, we focus on the more gradual evolution in decaying active regions—a potentially different scenario to the active region eruptions where the eruptivity index was previously studied, but nevertheless an important source of solar eruptions.

In a similar manner to the work of Pariat et al. (2017), we attempt to find a scalar quantity that acts as a proxy for the eruptivity of our 2.5D flux rope. The large number simulations we are able to run (in the order of 500) allows for an extensive parameter study based on the variation of the magneto-friction coefficient and the rate of photospheric flux cancellation. Using a probabilistic approach, a large number of eruptive and non-eruptive ropes will be compared against diagnostic measurements of the system at various points of their evolution. We note that, in addition to the formation of flux ropes, our model exhibits “arcade eruptions” whereby we observe periodic reconnection at the top of a sheared magnetic arcade. It is possible that these represent streamer blowouts or even “stealth CMEs,” as coined by Webb and Howard (2012), since they are characterised by the lack of a detectable signature on the solar surface. The nature and period of these eruptions (around every 25–30 days) matches well between our 2.5D model, global magneto-frictional simulations and observations (Bhowmik et al., 2021). However, the main focus of this paper is not on arcade eruptions but “flux rope eruptions,” when we observe magnetic reconnection below a flux rope and the rope itself moves upwards out of the domain. Such eruptions are larger than arcade eruptions and it is likely that in reality a significant number of CMEs occur as a result of such a mechanism.

We begin in Section 2 by outlining the mathematical basis of our model, including the mechanism by which we represent surface shearing (influenced by the differential rotation of the Sun at different latitudes), photospheric diffusion and the effect of the solar wind. We identify the variable parameters in the model and discuss which of these we can use to produce an array of differing flux rope behaviour. We then discuss the system diagnostics we have chosen to focus on, including a newly-defined measure of the relative helicity in two dimensions (Section 2.4.8). This measure allows us to compare our results against 3D equivalents and calculate the eruptivity index. The results of our parameter study are presented in Section 3, and the findings discussed in Section 4.

2 Methods

2.1 Magneto-Frictional Model

Rather than using full magnetohydrodynamics (MHD), we adopt a simplified approach that has been developed for global modelling of the solar coronal magnetic field: the magneto-frictional model (e.g., Mackay and Yeates, 2012). On a global scale, this technique has been applied with time-dependent lower boundary conditions from surface flux transport models (e.g., Yeates et al., 2008; Yeates, 2014). We adopt a 2.5-dimensional (translationally invariant) version of this approach.

In full MHD models, the velocity field, v, is determined by the momentum equation

where ρ is the fluid density, p is the fluid pressure, Ψ is the gravitational potential, and $\mathbf{j}={\mu }_0^{-1}\nabla \times \mathbf{B}$ is the current density. This is coupled to the induction equation

along with additional fluid equations to close the system. In the magneto-frictional method, Eq. 2 is retained, but the inertial terms, pressure gradients and gravity are neglected and instead a frictional velocity is imposed as

Combined with the induction equation, this leads to monotonic relaxation towards a stationary force-free field with j ×B = 0. The friction coefficient ν is typically given the form ν = ν₀|B|² (with some minimum value imposed) so that the overall evolution is independent of the magnitude of B and relaxation is not unduly slow near to magnetic null points (Yang et al., 1986).

In the outer corona, the solar wind outflow prevents the magnetic field from being force-free, but this effect can be approximated in the magneto-frictional model by relaxing towards an equilibrium with a specified outflow v_out, thus choosing v according to

This ad hoc approach was introduced by Mackay and van Ballegooijen (2006), and has subsequently been used in global magneto-frictional models of solar and stellar coronae (e.g., Yeates, 2014; Gibb et al., 2016; Mackay et al., 2018; Meyer et al., 2020).

2.2 Our 2.5-Dimensional Model

In our model, we simplify the problem by using a Cartesian coordinate system, and removing any dependence of the solution on the y coordinate. The domain ratio is taken to be (x:z) = (2:1). The resulting 2.5-dimensional system captures many of the essential features of the evolution in the lower solar corona, while affording a vast reduction in computational expense. This has allowed us to run an extensive parameter study, comprising hundreds of simulations.

The state of the system can be described in terms of a vector potential A(x, z, t), such that

The vector field A is then evolved according to

where E is the electric field, satisfying Ohm’s Law

in terms of the magnetic field B, current density j and frictional velocity v as given in Eq. 4. Here η is a constant representing coronal turbulent diffusivity, which is assumed to be much more significant than ohmic diffusivity in the highly-conducting corona, though still smaller than the other effects in the model. The outflow velocity used to represent the solar wind is taken to be

so the effect is minimal near the lower boundary and increases rapidly near the top boundary at z = z₁ at which the parameter v₁ gives the maximum speed.

The simulations are initialised with an “outflow field” (Rice and Yeates, 2021), a variation on a potential arcade that takes into account the effect of v_out, such that the system is initially in equilibrium. As seen in Figure 1 these fields are similar to potential arcades, except in the upper half of the domain where the field lines open out to become more vertical.

FIGURE 1. Examples of initial conditions for the 2.5D magneto-frictional code in a cartesian domain. Here (A) is a potential field and (B) is an equivalent outflow field, with the same lower boundary condition. The black lines represent magnetic field lines.

The boundary conditions determine the overall behaviour of the system. On the top boundary (z = z₁) we have the condition that B ⊥n = 0, ensuring that the magnetic field lines are vertical/radial here. This condition is consistent with observations that the field lines in the real corona become almost radial above a certain altitude. The notable exception to this is during eruptions, but these last a relatively short time and occur in our model even with the radial condition imposed. On the sides of the domain we set B ⋅n = 0, i.e., there is no magnetic flux through the sides. The lower boundary condition is more complex, and incorporates two effects: photospheric shearing and photospheric diffusion.

In our model, the photospheric shearing is assumed to originate from the differential rotation of the solar surface, namely that the Sun rotates more quickly at the equator than the poles. A magnetic arcade that has footpoints at different latitudes will thus be sheared, stretching the field lines along the polarity inversion line (PIL). We represent this shearing by imposing a velocity in the out-of-plane (y) direction at the lower boundary, following the profile

This profile is chosen as it ensures symmetry and simplicity while approximating the qualitative effect of differential rotation. The photospheric diffusion, η₀ represents the large-scale effect of supergranular flows, and is also imposed on the lower boundary (Mackay and van Ballegooijen, 2006). The combined effect of both shearing and diffusion results in the lower boundary condition

The effect of this is to bring the arcade footpoints closer together, eventually reconnecting to form a twisted flux rope (van Ballegooijen and Martens, 1989). Thus the ropes form self-consistently, without the need for any imposed flux emergence.

Numerically, Eq. 6 is solved on a finite-difference staggered grid (Yee, 1966) similarly to Cheung and DeRosa, 2012, with typical resolutions of 256 × 128 cells. The code is initialised with Python but runs using Fortran 90, making use of MPI parallelisation.

2.3 Parameters and Units

There are a number of free parameters in our model that can be easily varied. The aim is to produce a set of simulations that model realistic behaviour well, but provide enough variation to be able to make predictions of future behaviour. The parameters in the model are:

• η: The coronal diffusion.

• ν₀: The magneto-frictional friction coefficient.

• v₁: The imposed solar wind outflow speed.

• η₀: The photospheric diffusion.

• V_y0: The photospheric shearing velocity.

However, in our analysis we fix some of these parameters. Firstly, we take the coronal diffusion to be very small, η ≈ 1 × 10⁻⁶, so that its effect on short-term flux rope behaviour is negligible. (The corresponding diffusion time across the height of the domain would be 10⁶ time units.) Secondly, we effectively set the time unit by fixing the maximum shearing velocity V_y0 to unity. Finally, we fix the outflow speed v₁ = 50. This is reasonable since the outflow velocity has little effect on flux rope behaviour; the solar wind mainly serves to induce currents in the upper corona, at the top of the domain and far from the flux rope formation region. The only notable effect of an increased outflow velocity is a slight increase in the frequency of arcade eruptions. We choose v₁ = 50 to reflect the fact that the solar wind outflow is faster than the shearing velocity from differential rotation.

This leaves two remaining parameters: ν₀ and η₀. These can be varied considerably (between certain bounds) and still exhibit realistic behaviour. By running simulations with different combinations of these parameters we observe different flux rope behaviours, with a good mix of eruptive and non-eruptive simulations. The results from this parameter study are described in Section 3.

Note that we adopt dimensionless units throughout this paper, with a maximum photospheric shearing velocity of unity according to Eq. 9. For comparison to the real corona, one would choose the length unit—equivalently the height of the domain, z₁—and specify an observed shearing velocity caused by the differential rotation of the coronal arcade. These would then fix the time that corresponds to one dimensionless time unit in our paper. For example, if we take the angular velocity of differential rotation at latitude λ on the Sun to be Ω(λ) = 0.18–2.396 sin²λ − 1.787 sin⁴λ degrees per day (Snodgrass, 1983) and choose the latitudinal limits of our domain to be 10°–40°, this results in a maximum shearing velocity |v_ϕ| ≈ 0.086 km s⁻¹. Taking z₁ = 1.8 × 10⁵ km (half of the latitudinal extent) would then imply that a code time unit is of the order $\approx 25\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$.

2.4 System Diagnostics

In this section we briefly describe the diagnostic measurements of the system used to identify events and ultimately try to make predictions of future behaviour such as eruptions. There are innumerable measurements that could theoretically be taken of the state of the magnetic field, but for practical purposes we have selected nine, which are as follows:

2.4.1 Open Flux

The open flux is the sum of the unsigned magnetic flux through the top boundary of the domain. (Due to the symmetry of the system, the sum of the signed flux is zero.) The open flux increases due to the effect of the solar wind, as the field lines become more vertical near the top of the domain and fewer of them loop back down in the arcade.

2.4.2 Maximum Current

For practical computation we set μ₀ = 1 and define the current density as the curl of the magnetic field: j = ∇ ×B. As a diagnostic, we take the maximum value that this attains in the domain. A potential field has zero current and an outflow field only has current concentrated near the top boundary, due to the effect of the solar wind.

2.4.3 Rope Current

We also measure the integral of the current within the rope, in the direction of the rope axis (the y direction). In our two-dimensional model the rope is easily identifiable as it consists of the infinitely-long field lines that never reach either the photospheric or outer boundary. As the rope is in general in the lower half of the domain, this current is unaffected by the behaviour at the top boundary.

2.4.4 Magnetic Energy

The magnetic energy is defined in our unitless system as $E_{\mathrm{M}}={\int }_0^1{\int }_{-1}^1\frac{1}{2}|\mathbf{B}{|}^2\mathrm{d}x\mathrm{d}z$. A potential arcade corresponds to a minimum energy solution for given normal-field boundary conditions, and an outflow field has only a slightly increased energy (less than 1%). There is a large increase in magnetic energy during flux rope formation as the system evolves far from potentiality.

2.4.5 Free Magnetic Energy

In addition to the overall magnetic energy, we can calculate the “free magnetic energy,” defined as $E_{\mathrm{F}}=E_{\mathrm{M}}-{\int }_0^1{\int }_{-1}^1\frac{1}{2}|{\mathbf{B}}_P{|}^2\mathrm{d}x\mathrm{d}z$, where B_P is a potential magnetic field as calculated in Section 2.4.8 below. This quantity has the advantage that the contribution to the overall energy from the background magnetic field is lessened, such that the contribution from the non-potential flux rope is more significant.

2.4.6 Poloidal Rope Flux

This is a measure of the poloidal (in-plane) magnetic flux contained within the flux rope (the region with infinitely-long field lines), defined as the flux intersecting a chord between the centre of the rope and the edge of the rope (usually the lower edge of the domain).

2.4.7 Axial Rope Flux

This is defined as the integral of the magnetic flux in the rope (the region with infinitely-long field lines) in the y direction, parallel to the axis of the rope itself. This appears to be roughly correlated to the rope current.

2.4.8 Relative Helicity

The classical helicity within a volume V would be defined as h(V) = ∫_VA ⋅B dV, where A is the vector potential of B. This quantity is dependent in general on the gauge of A, and so we use the alternative relative helicity instead (Berger and Field, 1984). In a 3D domain, this would be calculated by finding a potential field B_P matching the original magnetic field on the boundary, and a corresponding vector potential, A_P. The relative helicity would then be

Care is required to define the relative helicity for our 2.5D field. A two-dimensional helicity measure for h(V) has been proposed before (Hu et al., 1997), but we are not aware of a previously published two-dimensional analogue for the relative helicity.

We start by considering the 3D formula (Eq. 11) on a finite volume $V_{y_1}$, where −1 < x < 1, − y₁ < y < y₁ and 0 < z < 1. Although B(x, z) in our 2.5D field is independent of y, the corresponding potential reference field B_P will, in general, vary in the y direction. This arises from the fact that it is potential, coupled with the need to match B_Py(x, ±y₁, z) = B_y(x, z) on y = ±y₁. We define the relative helicity per unit length to be

where B_P and its vector potential A_P are calculated on $V_{y_1}$. However, a physically meaningful helicity measure for our 2.5D field cannot possibly require integration in y. We will show that $H_R^{(2.5D)}$ not only converges as y₁ → ∞ but can indeed be calculated by a two-dimensional integral in x and z.

To do this, we decompose B_P into three components,

where the first component is a uniform field accounting for the net flux in the y direction,

and the other two components are both potential fields satisfying Δϕ₁ = Δϕ₂ = 0, with corresponding boundary conditions

Notice that ∇ϕ₀ and ∇ϕ₁ are independent of both y and y₁. The y dependence is concentrated only in ∇ϕ₂.

The potential ϕ₂ has the important property that, as y₁ increases, it becomes more and more concentrated near to the end boundaries y = ±y₁, irrespective of B_y(x, z). To see this, note that in the Cartesian domain $V_{y_1}$, the solution for ϕ₂ may be written as a Fourier series

where m, n are integers, k² = (m/2)² + n², and the sum includes all terms except m = n = 0 (which has been separated as the ϕ₀ component). The coefficients are then determined by the boundary condition on ∂ϕ₂/∂y, which gives

Now consider the value of ϕ₂(x, ay₁, z) for some fixed fraction |a| < 1. Then

where F(x, z) contains the x and z dependence from (Eqs 17, 18). Noting that

we see that the non-zero part of ∇ϕ₂ becomes an increasingly smaller fraction of the domain length as y₁ → ∞. This is illustrated in Figure 2. It follows that the contribution from ∇ϕ₂ to $H_R^{(2.5D)}$ in Eq. 12 vanishes in the limit, so that

where ${\mathbf{B}}_P^{(2.5D)}(x,z)=\nabla {\varphi }_1(x,z)+{\mathrm{\Phi }}_0{\mathbf{e}}_y$. The gauge of ${\mathbf{A}}_P^{(2.5D)}$ does not affect the integral (as usual for H_R), so one is free to choose an ${\mathbf{A}}_P^{(2.5D)}$ that is independent of y and hence evaluate $H_R^{(2.5D)}$ with a purely 2D integral. This is our approach in this paper.

FIGURE 2. Plot of the function $f(a)=\frac{\cosh (k\pi a{y}_1)}{\cosh (k\pi y_1)}$ for k = 1 and various values of y₁.

2.4.9 Eruptivity Index

The 3D expression that has been proposed (Pariat et al., 2017) as an eruptivity index is |H_J/H_R|, where

is the helicity of the current-carrying part of the field (Berger, 1999) and H_R is the relative helicity as above. We define the 2.5D version of H_J analogously to $H_R^{(2.5D)}$, thus we consider the ratio $|H_J^{(2.5D)}/H_R^{(2.5D)}|$ with

2.4.10 Ratios

We also consider the ratios between the above quantities, generally defined such that the ratios are independent of the magnetic field strength.

3 Results

3.1 Qualitative Behaviour

We first illustrate the two types of eruption that can occur in the system.

3.1.1 Zero Photospheric Diffusion

When we have no photospheric diffusion (η₀ = 0), the effect of the shearing causes the arcade footpoints to move only in the y direction. The magnetic energy increases as the system is no longer in a relaxed state, but it is not possible for flux ropes to form in the low corona as the footpoints are not brought closer together and there is no flux cancellation. In this case, the only free parameter is the friction coefficient ν₀.

We choose ν₀ of the order unity, in our code units. In general, the equivalent friction coefficents used in 3D global magneto-frictional simulations (e.g., Yeates and Hornig, 2016) are slightly smaller than this if directly compared, but this is by no means a precise measurement. There is thus a compromise between these more realistic values and the significant improvements in computational speed gained from increasing ν₀. Altering the domain shape also has a significant effect on the ideal ν₀, but for our (x, z) = (2:1) ratio setting ν₀ ≈ 0.5 produces realistic behaviour. For ν₀ significantly larger than this, the frictional relaxation is unrealistically slow relative to the footpoint shearing.

Consequently, our parameter study is confined to 0.5 < ν₀ < 2. For these ν₀ values, the magnetic arcade becomes more sheared in the y direction, as expected. After a certain period during which the magnetic energy and open flux increase gradually, we observe fast magnetic reconnection at the top of the domain accompanied by a sharp decrease in these diagnostics. This is what we refer to as an “arcade eruption” (Linker and Mikic, 1995). A sequence of snapshots of one of these eruptions is shown in Figure 3. The first eruption is shown at t = 2.92.

FIGURE 3. Snapshots showing the shearing of the potential arcade resulting in an arcade eruption and subsequent reformation preceding another eruption at t = 4. The black lines are magnetic field lines projected onto the (x, z) plane, and the heatmap represents the magnetic field strength into the page (in the y direction). In this simulation, η₀ = 0 and ν₀ = 0.5.

In general this process then repeats, leading to periodic eruptions occuring every time unit or so, as shown clearly by the oscillation of the diagnostic measures in Figure 4. We observe that the arcade becomes more sheared up to time t = 2.76, at which point there is an eruption, denoted by a blue circle. The process then repeats, building up to another eruption at t = 4 and so on, with the times of eruptions (as determined by the mid-point of the drop in open flux) shown by blue circles.

FIGURE 4. A selection of diagnostics for a simulation exhibiting repeated arcade eruptions, which are represented by blue circles. In this simulation, η₀ = 0 and ν₀ = 0.5. The time of an eruption is taken to be the midpoint in time between the maximum and minimum open flux values either side of the eruption.

During an arcade eruption there is a rapid decrease in open flux as the reconnection at the top of the domain results in the closing down of field lines at the top of the arcade. The corresponding decrease in magnetic energy occurs as these newly closed field lines are in a more potential state immediately after the eruption, although the free energy has a non-zero “floor” because the background field is non-potential due to the outflow velocity. The peaks in current during an eruption occur at the top boundary at the current sheet that temporarily forms above the PIL. The relative helicity does not follow the same pattern as the other diagnostics and instead increases during an eruption. This increase results from the fact that the newly-closed potential field lines at the top of the arcade now have a mutual “linkage” with the still sheared field lines beneath.

3.1.2 Nonzero Photospheric Diffusion

When the photospheric diffusion η₀ is nonzero, we observe the formation of flux ropes. These are characterised by twisted bundles of magnetic field lines that do not meet the boundary of the domain. As a result, the flux ropes are effectively infinitely long as there is no variation in the y direction. At the same time, we observe similar periodic buildups and releases in magnetic energy to the case with no photospheric diffusion. These overlying arcade eruptions do not affect the amount of poloidal magnetic flux within the flux rope, but they do cause the rope to oscillate vertically.

However, in addition to these eruptions, we also observe more significant events whereby the flux rope itself erupts, which in a 2D setting causes it to move rapidly upward out of the domain. Unlike with arcade eruptions, we observe magnetic reconnection below the flux rope, between it and the polarity inversion line on the lower boundary. There is a significant decrease in free magnetic energy, open flux and relative helicity, and in most cases there is no longer any flux contained within the rope after an eruption (although this is not always the case).

Snapshots showing a flux rope eruption are shown in Figure 5, and diagnostics corresponding to the same simulation in Figure 6. Here the red circles denote the flux rope eruptions, shown at the time when the poloidal rope flux drops. In most cases, after a flux rope eruption a new rope will form and the process will repeat, although depending on the specific parameters the system may just relax into a steady state. Compared to Figure 4, Figure 6 includes new plots of the rope flux and current as well as the current-carrying helicity $H_J^{(2.5D)}$ and eruptivity index $|H_J^{(2.5D)}/H_R^{(2.5D)}|$. The figure shows both the poloidal/in-plane flux and the axial/out-of-plane flux, contained within the region of infinitely-long magnetic field lines. In our model, prior to eruption, there is no null point beneath the rope. So the poloidal rope flux is simply the difference in the vector potential, max(A_y(x, z)) − A_y(0, 0).

FIGURE 5. Snapshots showing the formation of a flux rope and its subsequent eruption at time t = 8. This is followed by the formation of a second flux rope, which experiences an arcade eruption at t = 15.8, before the process repeats. The black lines are magnetic field lines projected onto the (x, z) plane, and the heatmap represents the magnetic field strength into the page (in the y direction). In this simulation, η₀ = 7 × 10⁻³ and ν₀ = 0.6.

FIGURE 6. A selection of diagnostics for a simulation exhibiting flux rope eruptions (red circles) interspersed with multiple arcade eruptions (blue circles). In this simulation, η₀ = 7 × 10⁻³ and ν₀ = 0.6. The time of a flux rope eruption is taken to be the time of the maximum poloidal rope flux before this rapidly decreases.

The effect of a flux rope eruption on the diagnostic measurements is more significant than an arcade eruption. We still observe rapid decreases in open flux and magnetic energy, but unlike with arcade eruptions there is now a significant decrease in the relative helicity during the eruption, rather than an increase. Arcade eruptions do not affect the poloidal flux in the rope at all as can be seen in Figure 6, but naturally as the rope no longer exists after a flux rope eruption these fluxes—as well as the axial rope current—immediately fall to zero, before increasing again if flux cancellation is ongoing. The axial rope current and axial rope flux are affected by arcade eruptions, but not significantly. The eruptivity index will be discussed below.

3.2 Dependence on Photospheric Diffusion

In the parameter study we primarily focus on the effect of changing the photospheric diffusion parameter η₀. For very small η₀ the system exhibits periodic arcade eruptions as described above, and the flux ropes form too slowly to erupt within the timeframe of the simulation, which we set at 50 time units (corresponding to several years on the real Sun). It is likely, however, that these ropes will eventually erupt.

When the photospheric diffusion is in the range 2 × 10⁻³ < η₀ < 1 × 10⁻¹ we observe both arcade and flux rope eruptions before t = 50. For runs where η₀ is not too large, there are several arcade eruptions between each flux rope eruption that cause the rope to oscillate. At some point, seemingly arbitrarily, one of these arcade eruptions coincides with reconnection below the rope and causes the flux rope to move rapidly upwards and out of the top of the domain. Usually after one of these eruptions the system returns to a potential-like state and the process repeats, although the out-of plane magnetic strength usually diminishes after each eruption. This repeats until the system relaxes to a steady state, with no rope present.

Varying the value of the photospheric diffusion (η₀) produces a wide range of flux rope behaviours. Figure 7 shows 100 simulations with ν₀ = 0.5 and varying η₀. For low η₀ (the left side of the figure), we observe regular arcade eruptions (indicated by the blue circles) similarly to the cases with no flux ropes. Small ropes do form in this region (indicated by the red lines), and it is likely that they would erupt completely given a sufficiently long integration time.

FIGURE 7. Overview of 100 flux rope simulations, for ν₀ = 0.5 and varying η₀. Each simulation (for a given η₀) is represented by a vertical red line, and the thickness of this line is proportional to the poloidal flux in the rope at that time. Arcade eruptions are represented by blue circles and flux rope eruptions are represented by red squares. The size of these points is proportional to the decrease in open flux and poloidal rope flux during the eruption, respectively.

In the region η₀ > 2 × 10⁻³ we observe flux rope eruptions before t = 50 (the end of the simulations). In general, these eruptions occur sooner as η₀ increases, although there is not a simple, continuous dependence of the eruption time on η₀ for all η₀. They are usually preceded by a number of arcade eruptions, varying in magnitude. In most cases, after a flux rope eruption a second flux rope forms, and we observe further eruptions of both types.

We observe that the time of the first flux rope eruption, if any, is negatively correlated to the diffusion η₀. This is visible in the pattern of red squares in Figure 7, which roughly follow a curve from the top-left of the diagram, flattening out as η₀ increases. We are unsure whether this curve tends to an asymptote at η₀ = 0 or at a finite, minimum value of η₀ below which flux rope eruptions do not occur. The time period between successive flux rope eruptions is also roughly negatively correlated to η₀.

For η₀⪆0.1 the flux ropes form and erupt very quickly, and there is no time for arcade eruptions (which have a frequency roughly independent of η₀). The flux ropes are small but have a high poloidal flux (resulting in the thicker vertical lines on Figure 7).

3.3 Diagnostic Measurements as Predictors of Eruptivity

Here we aim to quantify which, if any, of the considered diagnostic measurements in Section 2.4 can be used as a proxy for the eruptivity of a flux rope. In this section we only consider flux rope eruptions and not arcade eruptions, which are present at almost all times. We calculate the usefulness of each of the diagnostic measures using a probabilistic approach.

The data used here are from 500 simulation runs up to time t = 50, covering the parameter space 5 × 10⁻⁴ < η₀ < 0.5 and 0.5 < ν₀ < 2.0. For reference, a realistic value for ν₀ is likely of order unity. For ν₀ > 2, the relaxation is unrealistically slow relative to the driving, and eruptions become very frequent and less well-defined, with flux ropes not having time to form properly. Each simulation logs 5000 data points, for a total of 2.5 × 10⁶ measurements for each diagnostic. We separate each point according to whether it precedes a flux rope eruption, within a certain specified time cutoff.

With multiple flux rope eruptions observed within each simulation, we are able to observe a large distribution of eruption magnitudes and rope sizes/strengths. None of the diagnostics alone are enough to reliably predict imminent eruptions. The magnetic energy, helicity, and open flux all decrease significantly during flux rope eruptions, but in the preceding time interval there is no significant indication of an imminent eruption, and in fact most of the variation in these quantities is due to arcade eruptions. In contrast, the rope fluxes and currents build up from zero with each new rope, until they reach a threshold at which there is an eruption, making them more promising predictors. However, this threshold does not have a constant value, as we can see by observing in Figure 7 that there is large variation in the size of flux rope eruptions. Figure 6 also shows that the peak levels of the rope fluxes and currents tend to reduce for each successive eruption, even in a single simulation.

However, by calculating the ratios of (squared) rope fluxes and current to the other diagnostics, we observe that some of these ratios do indeed appear to have an eruptivity threshold with a roughly constant value. To motivate this observation, Figure 8 shows scatter plots of the points immediately preceding eruptions, with the appropriate diagnostics plotted against one other. The diagnostics on each axis are weighted to be proportional to the square of the magnetic field strength.

FIGURE 8. Pairwise scatter plots of the diagnostic values for all snapshots with later flux rope eruptions, in order to establish whether the ratios of the diagnostics are good predictors. The sizes of the points are weighted based on proximity to the eruption, such that the larger points are close to eruptions and vice versa.

In particular, we observe good linear correlations between the square of the axial rope current immediately prior to eruption, and both the magnetic energy and relative helicity of the system. This suggests that the ratios $I_y^2/H_R$, $I_y^2/H_J$, $I_y^2/E_M$ and $I_y^2/E_F$, where I_y is the axial rope current and E_M, E_F are the total and free magnetic energies respectively, should be good indicators of an imminent flux rope eruption. These particular ratios seem to work because the denominators are relatively steady between each successive flux rope eruption, but have a different level each time (this is visible in Figure 6). Since $I_y^2$ increases prior to a flux rope eruption, the ratios also increase but reach a common, normalized threshold between different eruptions. Note that H_R, H_J, E_M and E_F all include important contributions from the sheared background field surrounding the rope. Quantities derived only from the rope itself are not found to be successful denominators.

The distributions of eruptive and non-eruptive points for these ratios and a selection of the diagnostics are shown in Figure 9. At any point in the simulations when a flux rope is present, we check whether the rope will erupt within a time cutoff t = 0.3 (roughly equivalent to 8 days). If so, the diagnostics at that time are added to the histograms coloured red. The points corresponding to ropes that will not erupt are coloured blue. The areas under the histograms are corrected for an equal weighting of eruptive and non-eruptive ropes, but in principle this could be altered if the overall prior/uninformed probability of eruption were known. Thus the diagnostics with least red/blue overlap are better predictors of eruptivity than those with a large overlap.

FIGURE 9. Histograms of five of the diagnostic parameters, and five of the ratios between them. The points that precede an eruption within t = 0.3 time units are shown in red, and those that do not are shown in blue. The diagnostics with less red/blue overlap are better predictors of eruptions and vice versa. The eruptivity skill score E is given for each diagnostic. The curves are normalised to have an area of unity.

This can be quantified by assigning a “probability” of eruption to each value of each diagnostic. This is essentially the height of the red histogram curve on Figure 9 divided by the combined height of both curves, and is output as a number 0 < P_e < 1 for each data point. The probability P_e can then be compared against whether the rope will actually erupt or not to produce a skill score, defined as

If a diagnostic is a perfect indicator of eruptivity then it would have a skill score E = 1. If the diagnostic is a no better predictor than random chance then it would have skill score E = 0.5. The values of E for each diagnostic are given in the headings on Figure 9. We observe that the free energy and open flux are not suitable predictors, as they are little better than chance. The relative helicity and eruptivity index are slightly better, correctly predicting around 60% of eruptions within the time frame. We note, however, that these diagnostics are inversely correlated to the probability of eruption. This is likely because eruptions are more frequent for lower values of the helicity, open flux and energy, as smaller ropes tend to erupt more frequently in our model.

In contrast, the ratios of the squared axial current to the helicities and energies are excellent predictors in this time frame, all with skill scores above E = 0.86 and significantly better than the constituent diagnostics on their own. The correlations are also all positive, indicating that these ratios increase until they reach a given threshold—represented by the peaks in the red curves on Figure 9. The best predictor for t = 0.3 before an eruption (the time at which the predictors are best) is $I_y^2/H_R$, with a skill score of E = 0.93, or a predictive accuracy of 93.0%.

If we disregard the effect of arcade eruptions, the free energy and helicity are roughly constant during the periods between successive flux rope eruptions. However, the ratio between these quantities has also been identified as a suitable predictor, as seen in Figure 9. Immediately before an eruption the relative helicity decreases significantly, before the corresponding decrease in free energy, and thus there is a brief peak in the ratio between the two quantities.

The chosen time cutoff naturally has an effect on the skill scores of each of these quantities. This is illustrated in Figure 10, where the skill scores for each of the quantities are plotted against the time cutoff, up to t = 1.0 before the eruption. Of the raw diagnostics, the relative helicity is consistently the best predictor, with skill scores higher than E = 0.65 at all times. At almost all times the eruptivity index does not perform as well, and indeed becomes little better than chance for time cutoffs approaching t = 1.0. The free energy and open flux fare little better, with skill scores around E = 0.6.

FIGURE 10. Variation of the skill scores for each of the measured quantities, depending on the time cutoff within which an eruption must occur. The raw diagnostics and eruptivity index are plotted as solid lines, and the ratios are plotted as dashed lines.

In contrast, the ratios of the axial current to the chosen diagnostics perform consistently well, although their skill score decreases as the time cutoff increases. The quantities E_F/H_R and $I_y^2/H_R$ perform best at small time cutoffs, with maximum skill scores of E = 0.938 and E = 0.937 respectively. The ratios of the axial current to magnetic energy also perform well, and indeed for time cutoffs greater than t = 0.5 these ratios perform equally well or better than the other ratios. The maximum skill scores for most of the ratios occurs at a time cutoff of around t = 0.3, which corresponds to several days in reality. As such these predictors have potential for predicting solar eruptions a useful time before they occur.

4 Discussion

In an extensive parameter study, we have analysed the behaviour of two-dimensional magnetic flux ropes to examine which, if any, properties of the system can be used to predict whether or not the rope will erupt, an event which in reality is likely to cause a CME. For our simulations we used the magneto-frictional model rather than full MHD, reducing computational complexity further and allowing us to run thousands of simulations over a wide parameter space. We observe repeated arcade eruptions as in equivalent 3D simulations, as well as the formation of flux ropes. The main focus of our results is the prediction of the eruption of the flux ropes, rather than eruptions within the overlying arcade.

We have observed that the behaviour depends greatly on the values of the magneto-friction coefficient ν₀ and the photospheric diffusion η₀. In particular, flux ropes form and erupt more quickly for higher η₀, whereas the frequency of arcade eruptions within the overlying field is roughly independent of η₀. By comparing the probability of a rope erupting within a certain time to a number of diagnostic measurements, we have assigned a “skill-score” to each of them and the ratios between them. Of particular interest are the relative helicity (which we have newly defined for 2.5D systems) and its current-carrying component, H_J. The ratio of H_J to the total relative helicity constitutes the eruptivity index of Pariat et al. (2017).

We found that none of the diagnostics considered were by themselves good predictors of eruptivity, with skill scores not significantly greater than random chance. Of note, the eruptivity index has similar predictive skill to the other diagnostics, and is in fact negatively correlated to the likelihood of a flux rope eruption, unlike in the previous MHD simulations (Pariat et al., 2017; Zuccarello et al., 2018). However, when we consider the ratios between the diagnostics, weighted so as to be independent of the overall magnetic field strength, we find that certain of these can be good predictors.

Of particular note are the ratios with the axial current as the numerator, and either a helicity or energy measure as a denominator. The axial current indicates the “strength” of the rope, and the helicity/energy is in effect a measure of the strength of the background field. These denominators are roughly constant in between eruptions, whereas the axial current steadily increases as the rope becomes larger. Upon the ratio reaching a certain threshold, the rope will erupt. Notably, this threshold appears to be independent of both η₀ and ν₀. The ratio of free energy to relative helicity is also a very good predictor for eruptions in the immediate future, as a rapid decrease in relative helicity is often followed by an eruption.

4.1 Effect of Background Magnetic Field on the Eruptivity Index

In this section we propose a straightforward explanation for why Pariat et al. (2017) observed a high eruptivity index prior to the flux rope eruption, whereas we do not. We propose that this difference is due to the direction of the overlying magnetic field. We will consider a simple analytical model, which shows that the eruptivity index will naturally be higher when the background/overlying horizontal field direction is opposite to that of the arcade.

Two magnetic field configurations are presented in Figure 11. The left pane shows a configuration similar to the fields we generate naturally by shearing a potential field, where the background magnetic field is orientated in the same direction as the arcade. The right pane has the background field in the opposite direction, leading to a magnetic null point above the arcade. When flux ropes emerge into the second type of field, they are more likely to erupt. This was shown clearly by Pariat et al. (2017), who compared simulations with both orientations of the overlying field. We endeavour to show here that the right-hand field configuration fundamentally results in a higher eruptivity index, for a given sheared field component. By contrast, our simulations correspond to the left-hand field configuration, so even though they do erupt, this is not accompanied by a high eruptivity index.

FIGURE 11. Comparison between two magnetic arcades with overlying magnetic fields in opposite directions. The black lines represent the magnetic field projected into the (x, z) plane and the heatmap represents the magnetic field strength out of this plane.

The model magnetic field plotted in Figure 11 comes from the analytical expression $\mathbf{B}=\tilde{\mathbf{B}}+B_0{\mathbf{e}}_x$, where

with ξ = 4(x² + z²). The sheared, out-of-plane component, B_y, is fixed, and the only parameter is the strength of the background magnetic field, given by the (constant) parameter B₀. The left pane in Figure 11 has B₀ = 0.25, the right B₀ = −0.25. We proceed to observe the dependence on B₀ of $H_R^{(2.5D)}$ and $H_J^{(2.5D)}$, as defined in Section 2.4.8, and the resultant effect on the eruptivity index $|H_J^{(2.5D)}/H_R^{(2.5D)}|$.

Following Section 2.4.8, in order to calculate the relative helicity we choose a vector potential $\mathbf{A}=\tilde{\mathbf{A}}(x,z)-B_{0}z{\mathbf{e}}_y$, where $\nabla \times \tilde{\mathbf{A}}=\tilde{\mathbf{B}}$ and $\tilde{\mathbf{A}}$ is independent of B₀. Since the B₀e_x component of B is a potential field, we can similarly choose ${\mathbf{A}}_P^{(2.5D)}={\tilde{\mathbf{A}}}_P(x,z)-B_{0}z{\mathbf{e}}_y$. It follows that, for these fields, the “current-carrying” helicity $H_J^{(2.5D)}=\iint ({\mathbf{A}-\mathbf{A}}_P^{(2.5D)})\cdot ({\mathbf{B}-\mathbf{B}}_P^{(2.5D)})\mathrm{d}x\mathrm{d}z$ has no dependence on the background field B₀, because the B₀ terms from A and ${\mathbf{A}}_P^{(2.5D)}$ will cancel. For the relative helicity, however, $H_R^{(2.5D)}={\int }_V({\mathbf{A}+\mathbf{A}}_P^{(2.5D)})\cdot ({\mathbf{B}-\mathbf{B}}_P^{(2.5D)})\mathrm{d}V$, so the two terms add together and it does depend on B₀. Thus we can write

where $H_0^{(2.5D)}$ is the relative helicity with B₀ = 0. Since $H_0^{(2.5D)}$ and $H_J^{(2.5D)}$ have no dependence on B₀, the eruptivity index can simply be expressed as

For a sheared field with B_y non-zero (and non-uniform), we can clearly see that there will be a particular background field strength B₀ where the eruptivity index will become infinite as the denominator vanishes. For the magnetic field specified in Eqs 25–27 the constants take the values

which results in a peak in the eruptivity index at B₀ ≈ − 0.15, when in particular the overlying magnetic field is oppositely directed to the magnetic field in the upper part of the arcade (as in the right pane of Figure 11). By contrast, if the overlying magnetic field has the same direction as that in the arcade, so that B₀ > 0 (as in the left pane of Figure 11), then the denominator of (Eq. 29) will not become very small so the eruptivity index will not become large.

In all of our simulations—where the flux rope is formed by shearing of a pre-existing potential arcade—the background field has the same direction as that of the arcade, whether or not the flux rope erupts. Generalising from the analytical model with B₀ > 0, this explains why our eruptions are not preceded by a high eruptivity index.

The simulations of (Pariat et al., 2017), which were driven by flux emergence, included cases with both directions of background field. The eruptivity index behaved as predicted by the simple model in this section, but in that case only the cases with oppositely-directed field (and high eruptivity index) erupted. Our work shows that there is a whole class of eruptions that will not have high eruptivity index owing to the fact that they occur despite having the same direction of overlying field.

4.2 Implications for Space Weather Forecasting

The simplified nature of our 2.5D system means that any quantitative predictions will not be valid in 3D or indeed for any differing domain size or shape. However, the qualitative patterns of behaviour that we observe (such as those of the ratios of axial rope current to overall helicity/energy) should be equally valid in all systems, including global coronal models, where flux ropes are formed by footpoint shearing from differential rotation. Moreover, since these flux ropes are formed by gradual shearing over days to weeks, and are located in the magnetically-dominated low corona, we do not expect that the general conclusions would change significantly if we were to move from the magneto-frictional model to full MHD. For example, it has been shown that the linear stability criteria in magneto-frictional and MHD systems are the same (Craig and Sneyd, 1986). Nevertheless, to make quantitative predictions about specific 3D magnetic configurations on the Sun will require further work to understand how the behaviour of the diagnostics depends on the local coronal magnetic stucture.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author Contributions

OR did all of the calculations and worked out the proof in Section 2.4.8, under the guidance of AY who conceived the project. Both authors contributed to writing the paper.

Funding

OR was supported by a UKRI/STFC PhD studentship, and AY by UKRI/STFC research grant ST/S000321/1.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank the two referees for valuable comments that have greatly improved the paper.

References