© ESO 2022

1. Introduction

The solar dynamo generates flux with predominantly positive (right-handed) helicity in the southern hemisphere and predominantly negative (left-handed) helicity in the northern hemisphere (Hale 1925; Seehafer 1990; Pevtsov et al. 1995). New flux from the dynamo region emerges primarily at large, active-region scales into the atmosphere, where it undergoes a cascade to small scales driven by the motions in the convection zone (e.g., Parnell et al. 2009). However, the simultaneous cascade of the helicity is thought to proceed primarily toward large scales, although a weaker cascade toward small scales has also been found (Alexakis et al. 2006). While small-scale flux elements can submerge and re-emerge with the convection cells in the upper convection zone, the steep density stratification in the photosphere presents a barrier to the submergence of flux at the larger scales, with the dividing line estimated to lie at roughly 1 Mm (van Ballegooijen & Martens 1989). As a result, positive (negative) helicity is thought to accumulate in the large-scale coronal field in the southern (northern) hemisphere.

The Sun does have options to prevent the helicity from accumulating indefinitely, but two of them proceed very slowly, on the timescale of the magnetic cycle, because they are driven by the meridional flow. One is the annihilation of helicity when flux reconnects across the equatorial plane, which is facilitated by the equatorward migration of the active-region belts in the course of the solar cycle. The other is the transport of flux toward the adjacent pole, there joining the open flux, which can shed the helicity through torsional Alfvén waves propagating away from the Sun. As these processes operate on such a long timescale, the helicity would accumulate in each hemisphere to very high levels in the course of a cycle if there were not an option to shed the helicity much faster.

For the corona, whose flux is rooted in the photosphere and partly extends into the interplanetary space, the relative magnetic helicity is the relevant, gauge-invariant quantity to be considered (Berger & Field 1984). As the relative magnetic helicity is carried by the electric currents in the volume, which also carry the free magnetic energy, the accumulation of helicity is a plausible driver of the coronal field toward instability. Therefore, it has been conjectured that coronal mass ejections (CMEs) are the relevant process of helicity shedding from the corona (Rust 1994; Low 1996). CMEs presumably result from an instability of the coronal field on active-region scales (van Tend & Kuperus 1978; Forbes & Isenberg 1991; Kliem & Török 2006); they occur at least once in every major active region, as well as in the quiet Sun, and carry away a significant fraction of the source region flux in the form of a large-scale flux rope (Qiu et al. 2007).

The conjecture of helicity shedding by CMEs is observationally supported (Démoulin et al. 2002; Green et al. 2002) and widely accepted, although, to our knowledge, the effectiveness of erupting flux ropes at shedding helicity has not yet been investigated quantitatively. A preliminary simulation study of helicity shedding by a highly twisted, kink-unstable flux rope yielded a rather low effectiveness: only approximately one-fifth of the initial helicity was shed (Kliem et al. 2011). As many active regions experience only one CME in their lifetime, this effectiveness might be too small to prevent an unrealistic accumulation of helicity.

Here, we perform the first systematic investigation of helicity shedding by simulating the eruption and ejection of flux ropes susceptible to the far more relevant torus instability, and comparing this with the helicity shedding by a kink-unstable flux rope and by the flux rope in a data-constrained model of a specific solar event. Specifically, we focus on the dependence of the shedding on the initial partition into self- and mutual helicity.

2. Numerical methods

We model CMEs as ideal magnetohydrodynamic (MHD) instability of a force-free flux rope (Sakurai 1976; van Tend & Kuperus 1978; Kliem & Török 2006). The one-fluid ideal MHD equations are integrated numerically in the zero-beta limit, neglecting any thermodynamic effects, which are of only secondary importance for the helicity shedding, and neglecting gravity. Ohmic resistivity is also neglected to approach the almost ideal MHD conditions in the corona as closely as possible. Magnetic reconnection, required for a CME (Lin & Forbes 2000), is enabled by numerical diffusion in the vertical current sheet that steepens as the unstable flux rope rises. The resulting reconnection rate has been found to lie quite close (within a factor 2) to the observationally estimated one for two carefully modeled solar eruptions (Kliem et al. 2013; Hassanin & Kliem 2016), although the magnetic diffusion is only numerical and exceeds the coronal magnetic diffusion by many orders of magnitude. The equations are given in Török & Kliem (2003) for example, where the adopted numerical scheme is also detailed.

Our modified Lax-Wendroff scheme replaces the numerically stabilizing but highly diffusive Lax step by artificial diffusion (Sato & Hayashi 1979) which is of similar structure. The coefficient of the artificial diffusion is tailored for each simulation as a function of time, and if necessary also as a function of space, resulting in a reduced numerical diffusion of the magnetic field compared to the Lax term by a factor 10³. For numerical stability, the scheme also includes a small kinematic viscosity.

As the initial condition, we employ the analytical, approximate force-free equilibrium of a toroidal current channel and flux rope (the tokamak equilibrium) in the specific form suggested by Titov & Démoulin (1999) as a basic model for solar active regions. The current channel of major radius R and minor radius a is situated in the vertical (y − z) plane, with the center of the torus at r = (0, 0, −d). The stabilizing external poloidal (“strapping”) field is provided by a pair of magnetic point sources, ±q, located at r = ( ± L, 0, −d), respectively, whose resulting flux concentrations in the photosphere, {z = 0}, at x ≈ ±L model a bipolar solar active region. An external toroidal (“guide” or “shear”) field component is introduced by adding an antiparallel pair of dipoles under the footpoints of the flux rope. Their depth is ≈2.7R, determined such that the field line passing through the footpoints of the flux rope axis closely follows the axis in the corona. The model for the initial density is derived from the initial field through ρ₀(x) = B₀(x)^1.5, a choice implying a slow upward and sideward decrease of the initial Alfvén velocity from the core of the model active region, as on the Sun, and found to enable the modeling of CMEs in close quantitative agreement with observations (e.g., Török & Kliem 2005; Kliem et al. 2012, 2013). The initial velocity is set to zero, u₀ = 0.

This configuration is discretized on a Cartesian grid with a closed bottom boundary at z = 0, closed side boundaries at x = ±l_x and y = l_y, and an open upper boundary at z = l_z. Point-symmetric mirroring about the z axis of the variables in the boundary at y = 0 takes advantage of the symmetry in the configuration which is preserved throughout the evolution. Photospheric motions are negligible during the short timescales of solar eruptions, and so we set u(x, y, 0, t) = 0, which implies that the vertical photospheric field, B_z(x, y, 0), is conserved. The side boundaries are chosen to be closed for numerical convenience, as open boundaries can introduce numerical artifacts, especially if inflows develop; this requires the boundaries to be placed sufficiently far away such that the erupting flux never reaches them. Here u = 0 in the boundary, which preserves the normal field component and the density in the boundary, while the tangential field components are allowed to vary to minimize the influence of the boundary on the rising flux rope. The fluid advects freely across the upper boundary; this is implemented by using extrapolation of the velocity onto the ghost layers of the grid.

The MHD variables are normalized in the natural way (e.g., Török & Kliem 2003) by using the apex height of the geometric flux rope axis, h_a = R − d, as the length unit, the initial Alfvén velocity (and field strength; density) at the apex point of the flux rope axis, V_A0 = B₀(0, 0, h_a)/(μ₀ρ₀(0, 0, h_a))^1/2, as the velocity (and field; density) unit, and the resulting Alfvén time, τ_A = h_a/V_A0, as the unit of time. A core set of configurations uses a uniform geometry of the current channel, given by d = 0.1 and a = 0.6. To evaluate the influence of the geometry, two further sets, one with varying d, the other with varying a, are also included. For each configuration, L is numerically determined such that the initial configuration is close to marginal stability (slightly supercritical) with respect to the torus instability, with the lower and upper bounds on the critical value, L_cr, differing by less than 5%. For the kink-unstable comparison run, a smaller a yields a supercritical twist and L is chosen to be moderately subcritical. For each set of these parameters, the source strength, q, of the strapping field is given by the analytical equilibrium condition. The guide field strength is a free parameter and varied from run to run, as detailed in Sect. 4.2. The grid of 251 × 180 × 240 points for the half cube is stretched outward from the origin to permit a large box size, l_x = l_y = 32 and l_z = 40, while the central volume, which fully includes the initial flux rope, is well resolved, Δ_x, y, z = 0.04. The stretching also very efficiently damps any outward traveling perturbation (e.g., that from the initial relaxation of the analytical equilibrium), so that we do not observe any reflection of perturbations at the side boundaries.

3. Computation of helicities

We calculate the relative magnetic helicity (Berger & Field 1984; Finn & Antonsen 1985) in the rectangular box V = {x, y, z  :   − l_x ≤ x ≤ l_x,   − l_y ≤ y ≤ l_y,  0 ≤ z ≤ l_z} according to

$$\begin{array}{cc}\hfill H=& {\displaystyle {\int }_{\mathrm{Box}\mathrm{volume}}\mathit{A}·\mathit{B}{\mathrm{d}}^3\mathit{r}}\hfill \\ \hfill & +{\displaystyle {\int }_{\mathrm{Side}\mathrm{faces}}{\varphi }_{\mathrm{C}}{\mathit{A}}_{\mathrm{C}}\mathrm{d}\mathit{S}+{\int }_{\mathrm{Top}\mathrm{face}}(\mathit{A}\times {\mathit{A}}_{\mathrm{C}})\mathrm{d}\mathit{S},}\hfill \end{array}$$(1)

where A and A_C are special vector potentials for B and for the current-free field with the same normal component at the box boundary, B_C, and ϕ_C is a scalar potential for B_C, that is

$$\begin{array}{cc}& \mathit{B}=\mathrm{\nabla }\times \mathit{A},{\mathit{B}}_{\mathrm{C}}=\mathrm{\nabla }\times {\mathit{A}}_{\mathrm{C}}=-\mathrm{\nabla }{\varphi }_{\mathrm{C}},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \mathrm{\Delta }{\varphi }_{\mathrm{C}}=0,{\mathrm{\nabla }{\varphi }_{\mathrm{C}}·\widehat{\mathit{n}}|}_{\partial V}={-B_n|}_{\partial V},\hfill \end{array}$$(3)

where $\widehat{\mathit{n}}$ is the exterior unit normal at the boundary ∂V of V. The two vector potentials are given by

$$\begin{array}{c}\hfill \mathit{A}={\mathit{A}}_{\mathrm{C}}(x,y,0)-\widehat{\mathit{z}}\times {\displaystyle {\int }_0^z\mathit{B}(x,y,z^{\prime })\mathrm{d}{z}^{\prime }}\end{array}$$(4)

and

$$\begin{array}{c}\hfill {\mathit{A}}_{\mathrm{C}}=\widehat{\mathit{z}}\times \mathrm{\nabla }{\displaystyle ({\int }_0^z{\varphi }_{\mathrm{C}}(x,y,z^{\prime })\mathrm{d}{z}^{\prime }+g_{\mathrm{C}}(x,y)),}\end{array}$$(5)

where $\widehat{\mathit{z}}$ is the unit vector in the vertical (z) direction and g_C(x, y) is a solution of the two-dimensional Poisson equation

$$\begin{array}{c}\hfill ({\partial }_x^2+{\partial }_y^2)g_{\mathrm{C}}(x,y)=B_z(x,y,0).\end{array}$$(6)

Our procedure for calculating H is an extension of that derived by DeVore (2000) for the half space z ≥ 0, where H is given by the first line in Eq. (1). The adaption of the method of DeVore (2000), characterized in particular by the gauge $\mathit{A}·\widehat{\mathit{z}}={\mathit{A}}_{\mathrm{C}}·\widehat{\mathit{z}}=0$ of the vector potentials, to the case of a finite rectangular box was first proposed by Valori et al. (2012); Eqs. (1)–(6) are easily obtained from their expressions (see Appendix A). Other algorithms for calculating the relative magnetic helicity in a rectangular box as a functional of the field B in the box were suggested, for instance, by Rudenko & Myshyakov (2011), Thalmann et al. (2011), and Yang et al. (2013, 2018). A review of helicity calculation methods for finite volumes was given by Valori et al. (2016).

Defining the closed field in V,

$$\begin{array}{c}\hfill {\mathit{B}}_{\mathrm{cl}}=\mathit{B}-{\mathit{B}}_{\mathrm{C}},\end{array}$$(7)

and its vector potential, A_cl = A − A_C, the helicity can be decomposed into self helicity of the closed field,

$$\begin{array}{c}\hfill H_{\mathrm{self}}={\displaystyle {\int }_V{\mathit{A}}_{\mathrm{cl}}·{\mathit{B}}_{\mathrm{cl}}{\mathrm{d}}^3\mathit{r},}\end{array}$$(8)

and mutual helicity between the closed and open potential field,

$$\begin{array}{c}\hfill H_{\mathrm{mutual}}=2{\displaystyle {\int }_V{\mathit{A}}_{\mathrm{C}}·{\mathit{B}}_{\mathrm{cl}}{\mathrm{d}}^3\mathit{r};}\end{array}$$(9)

see Berger (1999). The relation H = H_self + H_mutual is satisfied by our numerically computed helicities to better than 0.3%. The closed field B_cl is identical to the current-carrying part of the field, and therefore the self helicity is also referred to as the current-carrying helicity.

4. Helicity shedding

4.1. Reference run: no guide field

First, we describe the evolution of a torus-unstable configuration with vanishing guide field, which sheds helicity quite effectively and is taken to be our reference (Run T1). For d = 0.1 and a = 0.6, marginal stability with respect to the torus instability is given for a “sunspot distance” slightly above the value of L = 1.3, which is chosen for this run to initiate the eruption. The threshold of the torus instability is given in terms of the decay index of the external poloidal field, n(z) = − dlnB_ep(z)/dlnz, at the apex point of the flux rope’s magnetic axis, which is slightly offset from the geometric axis, h_m = 1.1. Experience with many different systems has shown that the threshold value lies in the range n_cr ∼ 1 − 2 and depends on various parameters of the equilibrium (in ways yet to be determined). The canonical value of n_cr ∼ 3/2 (Bateman 1978) appears to provide a reasonable representative value for approximately semicircular flux ropes in the absence of a guide field (Török & Kliem 2007; Aulanier et al. 2010; Fan 2010; Zuccarello et al. 2015). For L = 1.3, we have n(z = h_m) = 1.29, which is slightly above the numerically determined threshold for the given configuration and is also close to the value of the full analytical expression for n(R) in Kliem & Török (2006). The twist, averaged over the cross-section of the current channel, is |Φ| = 2.45π, and sufficiently small for stability against the helical kink mode (Hood & Priest 1981; Török et al. 2004).

Figure 1 shows the rise profile of the fluid element at the apex point of the flux rope axis. The run starts with a phase of numerical relaxation that attempts to reach a numerical equilibrium from the initial approximate analytical equilibrium. The torus instability develops out of this equilibrium from t ≈ 18, initiated by the velocities that develop during the relaxation and act as a small perturbation. The occurrence of instability is clear from the initially nearly perfect exponential rise up to t ≈ 50. The flux rope reaches a rise velocity of u_z ≈ 0.6, and its apex point leaves the box at t ≈ 133. Figures 2 and 3 show snapshots of the flux rope, the initial current channel, the vertical current sheet spawned by the eruption, and the remaining volume currents. The latter are concentrated in the area of the two prominent flux bundles in the re-closed field of the model active region. These remain from the initial flux rope after the rope has reconnected with ambient flux in the vertical current sheet and connect the flux rope footpoints with the opposite-polarity sunspots (see the snapshots at t ≥ 82 in Fig. 2), which is close to the potential field (see, e.g., Fig. 6e in Hassanin & Kliem 2016). Such reconnection has been observed in many eruptions of a Titov-Démoulin flux rope (e.g., Török & Kliem 2005) and other models of solar eruptions (Aulanier & Dudík 2019) and certainly occurs on the Sun, although the frequency of its occurrence is yet to be determined.

Fig. 1. Rise profile of the flux rope apex in the reference run T1, showing height, h(t), and velocity, log uz(t).

Fig. 2. Snapshots of the unstable flux rope in Run T1 visualized by rainbow-colored field lines. Green field lines are started at the z-axis (with uniform spacing) to visualize the overlying flux and the vertical current sheet that forms under the rising flux rope. Pink field lines visualize the ambient flux in the vicinity of the flux rope, which recloses after the eruption. The magnetogram, Bz(x, y, 0), is shown in gray scale; a logarithmic scaling is used because the sunspot flux is rather concentrated due to the small value of d.

Fig. 3. Isosurfaces of current density, J, for Run T1 showing, from left to right: the initial flux rope; the vertical current sheet during the evolution of the instability (also wrapping around the bottom part of the flux rope legs in the form of a sigmoid); and the currents remaining after the eruption (which are already visible at t = 82).

Figure 4 shows the relative magnetic helicity in the box (Eq. (1)), normalized by the unsigned flux in the magnetogram area, F = (F₊ − F₋)/2, versus time for this run. The helicity begins to decrease rapidly when the upper edge of the current-carrying flux rope reaches the upper boundary at t ≈ 127, and levels off after the top part of the flux rope has nearly completely left the box at t ≈ 200. A major part (64%) of the initial helicity is shed by t ≈ 300 (the very slow subsequent increase in the helicity is considered to be due to the development of weak current density ripples in the central-difference numerical integration and is disregarded).

Fig. 4. Flux-normalized relative magnetic helicity in the box vs. time for the key runs analyzed in this paper (absolute values are shown). The initial helicity of Run S is scaled, placing it in the middle between H0/F2 for Runs T2 and T3 based on the similarity of the parameters. The length unit and resulting Alfvén time of Run S are rescaled from the values in Kliem et al. (2013) to conform to the choice for the other runs here (the apex height of 0.1 R⊙ at the onset of instability is used as length unit). The curves are labeled with the respective values of the ratio of external toroidal (guide) and external poloidal (strapping) field, Bet/Bep, decay index, n, at the initial flux rope apex, and twist averaged over the cross-section of the current channel, Φ. The range of the conjectured base level of the mutual helicity in the configurations T1–T3, ≈0.06 − 0.08, is shaded in light blue.

The initial helicity, H₀, in this configuration is partitioned into self- and mutual helicity as H_self, 0 = 0.33 H₀ and H_mutual, 0 = 0.67 H₀. The initial self helicity due to the twist in the initial flux rope completely leaves the system with the flux of the CME bubble (Fig. 2); there is even a slight overshoot (buildup of a small opposite self helicity) in this run, which is not essential in the total helicity budget. Of the initial mutual helicity, 41% is shed by the eruption, and the rest is carried by the volume currents in and around the two flux bundles that remain from the flux rope (Fig. 3 at t = 299).

The computation of the helicity was checked by extending the integrals in Eq. (1) to a lower part of the box only, z = 0 − 10. The initial helicity was found to be identical to the one in Fig. 4 (to within 1 × 10⁻³); that is, in the top part of the box, the computed potential field is nearly identical to the ambient field of the Titov-Démoulin flux rope, so that the contribution to the total relative helicity is negligible. The helicity versus time was found to follow the curve in Fig. 4 up to the arrival of the top part of the flux rope at the virtual boundary at t ≈ 73 and to show a slightly more rapid decrease by the same amount, ΔH/F² = 0.13, by t ≈ 155, i.e., the same shedding effectiveness.

4.2. Parametric study

4.2.1. Effect of the guide field

The strong difference in the effectiveness of helicity shedding between the torus-unstable configuration of Sect. 4.1 and the kink-unstable configuration in Kliem et al. (2011) motivates us to study the dependence of the shedding on the partition into mutual- and self helicity of the initial configuration. The twist and corresponding self helicity of the flux rope play at most a secondary role in the dynamics of the former case, but are the driver of the initial instability in the latter. First we study the effect of the partition on purely torus-unstable configurations, which are most likely relevant for the majority of CMEs. To this end, we take advantage of the fact that the strength of the external toroidal field component in the initial equilibrium can be varied freely without significantly changing the residual forces. In two further runs, we keep the geometry from Run T1 and set the strength of the external toroidal field such that its ratio to the equilibrium value of the external poloidal field at the apex of the flux rope axis is B_et/B_ep = 0.2 (Run T2) and B_et/B_ep = 0.5 (Run T3). As the external toroidal field has a stabilizing effect on the torus and helical kink instabilities, we must simultaneously reduce the sunspot distance L in order to raise the decay index at the position of the flux rope apex. The critical decay index for the two configurations is numerically found to lie slightly below n = 1.40 and 1.80, which respectively correspond to L = 1.2 and 0.9 chosen here. It should be noted that the stabilizing effect of the guide field is presumably somewhat overemphasized in the analytical equilibrium compared to the conditions in the solar corona because the observations of the shear in flare loops indicate that the guide field decreases strongly already above a few initial flux rope heights (e.g., Su et al. 2006; Warren et al. 2011). The parameters defining this core set of equilibria and of several configurations considered for comparison purposes are compiled in Table 1, together with the (absolute values of the) resulting helicities.

Figure 4 shows that the initial helicity decreases for increasing B_et/B_ep in Runs T1–T3. The mutual helicity decreases as the flux rope is more aligned with the ambient flux, and the self helicity decreases as the twist in the flux rope is reduced by the higher external toroidal field (Table 1). Here the mutual helicity contributes only slightly more to the overall decrease than the self helicity.

The fraction of helicity shed by the eruption decreases as well, which is due to two effects. First, as in Run T1, almost all of the initial self helicity is shed in Runs T2 and T3, but self helicity makes up a progressively smaller fraction of the initial helicity in these runs. Second, mutual helicity is shed less effectively, which appears to result from an inability to shed mutual helicity beyond a certain base level (see the rather similar final helicities for T1–T3 in Fig. 4). In other words, the system cannot fully approach the potential field. This base level of normalized mutual helicity in the configurations T1–T3 is ≈0.06 − 0.08, but its value differs for other configurations. It is smaller in Run K presented below, but is expected to be higher if the legs of the flux rope do not reconnect in the vertical current sheet (as in the 2D standard flare model). In this case, the strongest flux bundles corresponding to the potential field, that is, those connecting the flux rope footpoints with the sources of the strapping field, do not form.

The complete or nearly complete shedding of the initial self helicity from our configurations does not mean that the post-eruption state does not possess any remaining current density, helicity, or free energy. Currents remain in the two prominent flux bundles that connect the main polarities similarly to the potential field (Fig. 3 at t = 299). These flux bundles have the same handedness as the original flux rope. However, their ambient flux has built up a weak opposite handedness in a large volume (different from the purely potential ambient field of the original configuration), and so contributes self helicity of opposite sign. Opposing contributions to the total remaining self helicity are found in all of our runs. These cancel out nearly perfectly in Runs T3 and S and yield a weak self helicity opposite to the initial one in Runs T1, T2, and K. Their underlying currents permit a significant mutual helicity.

4.2.2. Kink-unstable case

For comparison, we include a kink-unstable run (Run K in Table 1), whose parameters are similar to the run in Kliem et al. (2011); except for a flatter flux rope shape (d = 0.83) and a somewhat larger minor radius, a = 0.33, in that run. A minor fraction of solar eruptions appear to be initiated by the onset of the helical kink instability (see, e.g., the event in Ji et al. 2003 and its modeling in Hassanin & Kliem 2016 and Hassanin et al. 2016, which yielded an estimate of Φ ≈ 4π for the average flux rope twist). As the helical kink saturates already for moderate displacements of the current channel, the subsequent onset of the torus instability is required for evolution into a CME; otherwise a confined eruption is produced. The helical kink lifts the flux rope into the torus-unstable height range in such CMEs; this is also the case in our simulation. For Run K, we keep the ratio B_et/B_ep = 0.5 from Run T3 and reduce the minor radius a to raise the twist. The sunspot distance L is set to a moderately subcritical value in order to initially stabilize the torus mode but allow its occurrence after the development of the helical kink. Both parts of the initial helicity are further reduced due to the smaller flux content of the thinner flux rope. This run also continues the trends in the effectiveness of the shedding found for Runs T1–T3: (1) self helicity, although shed completely, contributes less because it makes up an even smaller fraction of the initial helicity, and (2) the effectiveness of shedding mutual helicity is further reduced. As a result, the fraction of total helicity shed is only about one-quarter, which is comparable to the value of about one-fifth or less estimated for the configuration in Kliem et al. (2011).

4.2.3. Effect of flux rope geometry

To cover a broader range of the conditions on the Sun, we also consider two sets of very weakly torus-unstable equilibria with varying geometric parameters of the flux rope. The kink-unstable configuration of Run K indicates that the ratio of the flux in the current channel, F_cc, and external poloidal flux, F_ep, which is small in this run due to the small minor radius, may considerably influence both the initial helicity and the conjectured base level of the mutual helicity, and therefore the effectiveness of the shedding. To study these effects, the minor radius is varied in the range a = 0.45 − 1. This corresponds to a flux ratio ℛ = F_cc/F_ep = 0.26 − 1.06 (ℛ = 0.42 in Run T1). In the second set, the rope is made flatter by setting d = 0.375 and d = 0.833, which raises the half-separation of its footpoints from D_f = 1.1 in Run T1 to D_f = 1.3 and 1.6, respectively. Both sets of runs use B_et/B_ep = 0, and are to be compared to Run T1.

Figure 5 shows a significant increase in the normalized helicity with increasing flux ratio, up to a value of H₀/F² = 0.29 for a = 0.9, which we refer to as Run T0. The helicity remaining after the ejection of the flux rope shows only a weak increase with ℛ. Consequently, the effectiveness of helicity shedding increases from ΔH/F² = 0.084 (54%) for a = 0.45 to ΔH/F² = 0.20 (70%) for a = 0.9. As in Run T1, all of the self helicity but only part of the mutual helicity (17−57%) is shed. These configurations possess a very weakly varying critical “sunspot distance”, lying in the range 1.3 < L_cr < 1.35 for a = 0.45 − 0.9 and within 1.35 < L_cr < 1.4 for a = 1. The decay index values n(z = h_m) given in Fig. 5 correspond to the lower edges of these ranges and vary slightly more than L because the position of the magnetic flux rope axis increases slightly with increasing flux rope thickness (the large-aspect-ratio approximation underlying the expression for the tokamak equilibrium degrades as a increases).

Fig. 5. Same as Fig. 4 for the additional sets of weakly torus-unstable runs with varying minor radius (blue curves) and flatness (red curves) of the initial current channel.

The highest normalized helicity of a stable configuration in our parametric study is H₀/F² = 0.27, found for a = 0.9 with d = 0.1 and L = 1.35. When analyzed with higher spatial resolutions of Δ_x, y, z = 0.02 and 0.01, this value reduces slightly to H₀/F² = 0.25 and appears to converge at this value. (It is worth noting that a decrease in the box size by factors of 2 and 4 in each direction, keeping the resolution at Δ_x, y, z = 0.04, does not change the computed helicities). The corresponding values of other relevant parameters are H_self, 0/H₀ = 0.315 and ℛ = F_cc/F_ep = 0.94. To our knowledge, these are the highest values of H/F² and ℛ found so far for stable force-free flux rope equilibria. However, such high-arching and simultaneously very thick flux ropes (filling all space between the flux rope axis and the photosphere) may be rare on the Sun. The observations of extreme ultraviolet (EUV) hot channels indicate smaller cross-sections, typically of the order of a ∼ 0.5, when the channel arches high up in the corona before erupting (e.g., Reeves & Golub 2011; Cheng et al. 2013; Patsourakos et al. 2013).

Both parts of the normalized helicity decrease as flatter flux ropes are considered. At first sight, this is counter-intuitive because the longer flux rope has a higher twist combined with the same toroidal flux and because the area under the flux rope, through which the external flux F_ep links with the flux rope, increases with D_f. Both self and mutual helicity indeed increase with increasing d, but the flux also increases (because the larger distance to the sources of the strapping field requires a higher source strength in the adopted normalization) and both normalized helicities decrease. This can perhaps be best understood by envisioning a normalization to R, that is, a fixed torus with a fixed helicity in the whole space. Increasing d raises the photospheric level toward the upper edge of the torus, so that a progressively smaller fraction of the helicity remains in the coronal half space. The initial partition into self and mutual helicity is the same as in Run T1, but the critical sunspot distance decreases, and the decay index increases (both do so weakly in the considered range), because the larger footpoint area of the flatter flux ropes increases the stabilizing effect of the line-tying.

Again, the self helicity is shed completely, indicating the existence of a base level of mutual helicity. This level also decreases weakly from the one in Run T1, meaning that the effectiveness of helicity shedding decreases moderately from 64% in Run T1 to 53% for the flattest flux rope in our set of equilibria. Within the range of parameters considered, changing the flatness of the flux rope has a weaker influence on the effectiveness of helicity shedding than changing the guide field or the flux content of the rope.

Finally, in Run T4, we include a guide field in one of the weakly twisted, flatter flux rope equilibria in order to obtain a more reliable estimate of the minimum helicity shedding by torus-unstable flux ropes. Here we choose B_et/B_ep = 0.5 and d = 0.375 because the equilibrium with d = 0.833 is completely stable for this guide field strength. The shedding effectiveness of 38% is found to be only slightly lower than in Run T3, and the remaining mutual helicity also stays in the range of ∼0.06 − 0.08 obtained in Runs T1–T3 (see Table 1 and Fig. 5).

4.3. Model of a solar event

To judge the parametric study against the conditions on the Sun, we repeat the modeling of the CME on 2010 April 8 from Active Region NOAA 11060 (Kliem et al. 2013) with an open upper boundary. This was a weak event (a relatively slow CME with a projected velocity of ≈520 km s⁻¹, associated with a weak flare of magnitude B3.7) from a decaying, approximately bipolar active region. Due to the simplicity of the source region, the event should be comparable to our model runs above. It was estimated that the eruption was driven by the torus instability from a threshold value of n ≈ 1.5 − 1.75, similar to the range of thresholds in Runs T2 and T3. The modeling of the source region magnetic field indicated a flux rope with a weak to moderate guide field, also similar to runs T2 and T3, and a very small average twist of ∼0.5π. Approximating the external field by the potential field, B_et/B_ep ∼ 0.4 at the axis of the flux rope in the center of the active region. For this simulation (Run S in Table 1), we obtain a helicity shedding of 0.42 H₀, similar to Runs T2–T3 as well.

The brief increase of the helicity in this run during t ≈ 100 − 130 is a signature of the steepening of the current layer in front of the erupting flux rope (see Fig. 5 in Hassanin & Kliem 2016 for the visualization of a similar configuration). By Lenz’s law, the toroidal direction of the current in this layer is opposite to the toroidal flux rope current. However, the toroidal field direction is the same as inside the rope. Therefore, the current layer locally changes the handedness of the field to the opposite of the flux rope’s handedness. When this flux leaves the box, the total helicity increases weakly and briefly until the flux rope begins to leave the box.

It is noteworthy that the normalized initial helicity lies far below the values in our model configurations. This results from the fact that the solar active region contains a large amount of dispersed and rather remote flux that does not contribute to the equilibrium of the flux rope. Such “exterior flux” does not exist in the Titov-Démoulin equilibrium, where all external flux passes over the polarity inversion line that hosts the flux rope. Exterior flux exists typically in the source regions of solar eruptions, except during the early stage of active-region emergence, which on the other hand does not (yet) launch CMEs. The typical existence of exterior flux can also be seen from the small maximum value of ≈0.022 for the normalized helicity injected into solar active regions through the photosphere prior to eruptions (LaBonte et al. 2007). Our configuration in T1 contains an approximately ten times higher normalized helicity very near the marginal stability state.

5. Discussion

The ratio B_et/B_ep of the external toroidal (guide) to the external poloidal (strapping) field is a key parameter in determining self- and mutual helicity of force-free flux rope equilibria and their stability against the torus mode. Here, both B_et/B_ep and n are varied in ranges that appear to be representative of a large fraction of solar eruptions (e.g., Cheng et al. 2020), and so may be the range of helicity shedding found, ∼(0.4 − 0.65)H₀.

The geometric range covered should be representative as well. Flux ropes forming in filament channels tend to be thick (e.g., Savcheva & van Ballegooijen 2009; Su et al. 2011), as in our studied range, and flat, flatter than can reasonably be studied using the Titov-Démoulin equilibrium. However, the slow rise of filaments or EUV hot channels up to the onset height of eruption always yields a considerably arched shape (Nindos et al. 2015; Cheng et al. 2020), with a semicircular shape being not uncommon (e.g., Schrijver et al. 2008; Cheng et al. 2013).

The essentially complete shedding of the initial self helicity in our simulation runs is facilitated by the open top boundary. When the top part of the flux rope has left the box completely, the legs can untwist freely. On the Sun, such untwisting is possible if the flux rope reconnects with open ambient flux. However, even if the flux rope remains intact, nearly all of its twist propagates into the interplanetary space. A significant but still minor part of the initial twist and self helicity remain in the corona only if the flux rope legs reconnect with each other (Hassanin & Kliem 2016; Hassanin et al. 2016) or with closed ambient flux.

The incomplete shedding of mutual helicity in our runs suggests that a base level of mutual helicity may exist that limits the fraction of helicity that can be shed by erupting flux ropes. This level appears to be weakly dependent on the parameters of the initial equilibrium, falling in the range H_mutual/F² ∼ 0.06 − 0.08 for typical source region parameters of solar eruptions. An incomplete shedding of helicity is indeed typically indicated by the considerable shear of the first forming flare loops (Su et al. 2007) and has also been found in simulations of long-term coronal evolution (Mackay & Yeates 2012). Further study is needed to understand its origin.

Overall, we find that a significant part of the initial relative helicity is shed by a typical flux rope eruption. However, the remaining part is significant as well, on the order of H₀/2. Therefore, further mechanisms may be relevant in regulating the helicity budget of the corona (Sect. 1) and merit a quantitative study.

6. Conclusions

Our parametric simulation study of flux rope eruption by ideal MHD instability of a line-tied tokamak equilibrium (Titov & Démoulin 1999) demonstrates the helicity shedding conjectured for CMEs. For the most relevant case of torus-unstable flux ropes in bipolar source regions and expected typical partitions of the initial relative helicity into mutual and self helicity in the range H_mutual, 0/H_self, 0 = 2.1 − 3.9 (corresponding to ratios of the equilibrium guide and strapping field components in the range B_et/B_ep = 0 − 0.5), we find that a fraction of ∼(0.4 − 0.65)H₀ of the initial helicity H₀ is shed. A data-constrained simulation of a solar eruption (from Active Region NOAA 11060 on 2010 April 8) yields a shedding of 0.42 H₀, which is consistent with the model results.

While the initial self helicity (flux rope twist) is completely or nearly completely shed, the effectiveness of shedding mutual helicity decreases with a decreasing initial flux-normalized value of this component. This appears to result from an inability to shed mutual helicity beyond a certain base level, which depends on the parameters of the initial equilibrium and is in the range of ∼0.06 − 0.08 for our torus-unstable model configurations.

The parametric study yields stable equilibria up to a normalized helicity of H/F² = 0.25, which, to our knowledge, is higher than what has been found so far (Fan 2010). For this equilibrium, the ratio of the flux in the current channel to the external poloidal flux, ℛ = F_cc/F_ep = 0.94, is significantly higher than previous estimates of the limiting flux ratio, which lie in the range of ∼0.1 − 0.5 (Bobra et al. 2008; Savcheva et al. 2012; Kusano et al. 2020). The corresponding ratio of self (current-carrying) helicity to total helicity is H_self/H = 0.315, which is only slightly higher than a previously suggested limit of ≃0.29 ± 0.01 for torus-unstable flux ropes (Zuccarello et al. 2018), but lower than a limit of ∼0.45 for eruptions from newly emerging flux (Pariat et al. 2017).

The flux-normalized helicity in solar data depends strongly on the amount of exterior flux in the magnetogram, that is, flux that is not part of the flux rope equilibrium. Such flux must be separated in order to find a reliable estimate of H/F².

Acknowledgments

We acknowledge the careful reading and constructive comments by an anonymous referee. This work was supported by the DFG and by NASA grants 80NSSC17K0016, 80NSSC18K1705, 80NSSC19K0860, 80NSSC19K0082, and 80NSSC20K1274.

References

Appendix A: Derivation of Eqs. (1)–(6)

We start from Eqs. (19–22) in Valori et al. (2012), which can be written as

$$\begin{array}{c}\hfill \mathit{A}=\mathit{b}+\widehat{\mathit{z}}\times {\displaystyle {\int }_z^{z_2}\mathit{B}\mathrm{d}{z}^{\prime },{\mathit{A}}_{\mathrm{C}}={\mathit{b}}_{\mathrm{C}}+\widehat{\mathit{z}}\times {\int }_z^{z_2}{\mathit{B}}_{\mathrm{C}}\mathrm{d}{z}^{\prime },}\end{array}$$(A.1)

where

$$\begin{array}{c}\hfill \mathit{b}=\left(\begin{array}{c}{b}_x(x,y)\\ b_y(x,y)\\ 0\end{array}\right),{\mathit{b}}_{\mathrm{C}}=\left(\begin{array}{c}{b}_{\mathrm{C},x}(x,y)\\ b_{\mathrm{C},y}(x,y\\ 0\end{array}\right)\end{array}$$(A.2)

satisfy

$$\begin{array}{c}\hfill \mathrm{\nabla }\times \mathit{b}=\mathrm{\nabla }\times {\mathit{b}}_{\mathrm{C}}=B_z(x,y,z_2)\widehat{\mathit{z}}.\end{array}$$(A.3)

With the choice z₂ = 0 and b = b_C, Eq. (A.1) becomes

$$\begin{array}{c}\hfill \mathit{A}={\mathit{b}}_{\mathrm{C}}-\widehat{\mathit{z}}\times {\displaystyle {\int }_0^z\mathit{B}\mathrm{d}{z}^{\prime },{\mathit{A}}_{\mathrm{C}}={\mathit{b}}_{\mathrm{C}}-\widehat{\mathit{z}}\times {\int }_0^z{\mathit{B}}_{\mathrm{C}}\mathrm{d}{z}^{\prime },}\end{array}$$(A.4)

where b_C obeys

$$\begin{array}{c}\hfill \mathrm{\nabla }\times {\mathit{b}}_{\mathrm{C}}=B_z(x,y,0)\widehat{\mathit{z}}.\end{array}$$(A.5)

As

$$\begin{array}{c}\hfill {\mathit{b}}_{\mathrm{C}}(x,y)={\mathit{A}}_{\mathrm{C}}(x,y,0),\end{array}$$(A.6)

the first equation in Eq. (A.4) is identical to Eq. (4).

As suggested by Valori et al. (2012), we choose b_C solenoidal, which is possible because both the curl (given by Eq. (A.5)) and the divergence of b_C can be prescribed. Therefore, b_C can be represented in the form

$$\begin{array}{c}\hfill {\mathit{b}}_{\mathrm{C}}=\left(\begin{array}{c}-{\partial }_y{g}_{\mathrm{C}}\\ {\partial }_x{g}_{\mathrm{C}}\\ 0\end{array}\right)=\widehat{\mathit{z}}\times \mathrm{\nabla }{g}_{\mathrm{C}}\end{array}$$(A.7)

by means of a stream function g_C(x, y). The second equation in Eq. (A.4) then becomes

$$\begin{array}{c}\hfill {\mathit{A}}_{\mathrm{C}}=\widehat{\mathit{z}}\times \mathrm{\nabla }{g}_{\mathrm{C}}+{\displaystyle {\int }_0^z\mathrm{\nabla }{\varphi }_{\mathrm{C}}\mathrm{d}{z}^{\prime }=\widehat{\mathit{z}}\times \mathrm{\nabla }(g_{\mathrm{C}}+{\int }_0^z{\varphi }_{\mathrm{C}}\mathrm{d}{z}^{\prime }),}\end{array}$$(A.8)

which is identical to Eq. (5), while Eq. (A.5) takes the form of the Poisson equation

$$\begin{array}{c}\hfill ({\partial }_x^2+{\partial }_y^2)g_{\mathrm{C}}(x,y)=B_z(x,y,0),\end{array}$$(A.9)

which is identical to Eq. (6).

The relative magnetic helicity in the box volume V is calculated using the formula of Finn & Antonsen (1985):

$$\begin{array}{cc}\hfill H=& {\displaystyle {\int }_V(\mathit{A}+{\mathit{A}}_{\mathrm{C}})·(\mathit{B}-{\mathit{B}}_{\mathrm{C}}){\mathrm{d}}^3\mathit{r}={\int }_V\mathit{A}·\mathit{B}{\mathrm{d}}^3\mathit{r}}\hfill \\ \hfill & -{\displaystyle {\int }_V{\mathit{A}}_{\mathrm{C}}·{\mathit{B}}_{\mathrm{C}}{\mathrm{d}}^3\mathit{r}+{\int }_V({\mathit{A}}_{\mathrm{C}}·\mathit{B}-\mathit{A}·{\mathit{B}}_{\mathrm{C}}){\mathrm{d}}^3\mathit{r}.}\hfill \end{array}$$(A.10)

For the first integral in the second line of Eq. (A.10) we have

$$\begin{array}{cc}\hfill {\displaystyle {\int }_V{\mathit{A}}_{\mathrm{C}}·{\mathit{B}}_{\mathrm{C}}{\mathrm{d}}^3\mathit{r}}& =-{\displaystyle {\int }_V\mathrm{\nabla }{\varphi }_{\mathrm{C}}·{\mathit{A}}_{\mathrm{C}}{\mathrm{d}}^3\mathit{r}}\hfill \\ \hfill & =-{\displaystyle {\int }_V\mathrm{\nabla }·\left({\varphi }_{\mathrm{C}}{\mathit{A}}_{\mathrm{C}}\right){\mathrm{d}}^3\mathit{r}}\hfill \\ \hfill & =-{\displaystyle {\int }_{\partial V}{\varphi }_{\mathrm{C}}{\mathit{A}}_{\mathrm{C}}·\mathrm{d}\mathit{S}}\hfill \\ \hfill & =-{\displaystyle {\int }_{\mathrm{Side}\mathrm{faces}}{\varphi }_{\mathrm{C}}{\mathit{A}}_{\mathrm{C}}·\mathrm{d}\mathit{S},}\hfill \end{array}$$(A.11)

where we make use of ∇ ⋅ A_C = 0 and ${\mathit{A}}_{\mathrm{C}}·\widehat{\mathit{z}}=0$. The second integral in the second line of Eq. (A.10) can finally be rewritten as

$$\begin{array}{cc}\hfill {\displaystyle {\int }_V({\mathit{A}}_{\mathrm{C}}·\mathit{B}-\mathit{A}·{\mathit{B}}_{\mathrm{C}}){\mathrm{d}}^3\mathit{r}}& ={\displaystyle {\int }_V\mathrm{\nabla }·(\mathit{A}\times {\mathit{A}}_{\mathrm{C}}){\mathrm{d}}^3\mathit{r}}\hfill \\ \hfill & ={\displaystyle {\int }_{\partial V}(\mathit{A}\times {\mathit{A}}_{\mathrm{C}})·\mathrm{d}\mathit{S}}\hfill \\ \hfill & ={\displaystyle {\int }_{\mathrm{Top}\mathrm{face}}(\mathit{A}\times {\mathit{A}}_{\mathrm{C}})·\mathrm{d}\mathit{S},}\hfill \end{array}$$(A.12)

where we take advantage of the fact that A × A_C is purely vertical (parallel to the z axis) and A(x, y, 0) = A_C(x, y, 0). Inserting Eqs. (A.11) and (A.12) into Eq. (A.10) leads to Eq. (1).

All Tables

All Figures

	Fig. 1. Rise profile of the flux rope apex in the reference run T1, showing height, h(t), and velocity, log uz(t).
In the text

Fig. 2. Snapshots of the unstable flux rope in Run T1 visualized by rainbow-colored field lines. Green field lines are started at the z-axis (with uniform spacing) to visualize the overlying flux and the vertical current sheet that forms under the rising flux rope. Pink field lines visualize the ambient flux in the vicinity of the flux rope, which recloses after the eruption. The magnetogram, Bz(x, y, 0), is shown in gray scale; a logarithmic scaling is used because the sunspot flux is rather concentrated due to the small value of d.

In the text

	Fig. 3. Isosurfaces of current density, J, for Run T1 showing, from left to right: the initial flux rope; the vertical current sheet during the evolution of the instability (also wrapping around the bottom part of the flux rope legs in the form of a sigmoid); and the currents remaining after the eruption (which are already visible at t = 82).
In the text

Fig. 4. Flux-normalized relative magnetic helicity in the box vs. time for the key runs analyzed in this paper (absolute values are shown). The initial helicity of Run S is scaled, placing it in the middle between H0/F2 for Runs T2 and T3 based on the similarity of the parameters. The length unit and resulting Alfvén time of Run S are rescaled from the values in Kliem et al. (2013) to conform to the choice for the other runs here (the apex height of 0.1 R⊙ at the onset of instability is used as length unit). The curves are labeled with the respective values of the ratio of external toroidal (guide) and external poloidal (strapping) field, Bet/Bep, decay index, n, at the initial flux rope apex, and twist averaged over the cross-section of the current channel, Φ. The range of the conjectured base level of the mutual helicity in the configurations T1–T3, ≈0.06 − 0.08, is shaded in light blue.

In the text

	Fig. 5. Same as Fig. 4 for the additional sets of weakly torus-unstable runs with varying minor radius (blue curves) and flatness (red curves) of the initial current channel.
In the text