11institutetext: Institute of Physics and Astronomy, University of Potsdam,
Potsdam 14476, Germany

# Helicity shedding by flux rope ejection

B. Kliem  N. Seehafer

(Received 12 October 2021 / Accepted 13 November 2021)

###### Abstract

Aims. We quantitatively address the conjecture that magnetic helicity must be
shed from the Sun by eruptions launching coronal mass ejections in order to
limit its accumulation in each hemisphere.

Methods. By varying the ratio of guide and strapping field and the flux rope
twist in a parametric simulation study of flux rope ejection from
approximately marginally stable force-free equilibria, different ratios of
self- and mutual helicity are set and the onset of the torus or helical kink
instability is obtained.

Results. The helicity shed is found to vary over a broad range from a minor to
a major part of the initial helicity, with self helicity being largely or
completely shed and mutual helicity, which makes up the larger part of the
initial helicity, being shed only partly. Torus-unstable configurations with
subcritical twist and without a guide field shed up to about two-thirds of the
initial helicity, while a highly twisted, kink-unstable configuration sheds
only about one-quarter. The parametric study also yields stable force-free
flux rope equilibria up to a total flux-normalized helicity of 0.25, with a
ratio of self- to total helicity of 0.32 and a ratio of flux rope to external
poloidal flux of 0.94.

Conclusions. These results numerically demonstrate the conjecture of helicity
shedding by coronal mass ejections and provide a first account of its
parametric dependence. Both self- and mutual helicity are shed significantly;
this reduces the total initial helicity by a fraction of
$\sim\\!0.4\mbox{--}0.65$ for typical source region parameters.

###### Key Words.:

Instabilities – Magnetic fields – MHD – Sun: corona – Sun: coronal mass
ejections (CMEs) – Sun: flares

##  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^{3}$. 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 $\bm{r}=(0,0,-d)$. The stabilizing external poloidal
(“strapping”) field is provided by a pair of magnetic point sources, $\pm q$,
located at $\bm{r}=(\pm L,0,-d)$, respectively, whose resulting flux
concentrations in the photosphere, $\\{z\\!=\\!0\\}$, at $x\approx\pm 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 $\approx 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 $\rho_{0}(\bm{x})=B_{0}(\bm{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, $\bm{u}_{0}=0$.

This configuration is discretized on a Cartesian grid with a closed bottom
boundary at $z=0$, closed side boundaries at $x=\pm 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 $\bm{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 $\bm{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_{\mathrm{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_{\mathrm{A0}}=B_{0}(0,0,h_{\mathrm{a}})/(\mu_{0}\rho_{0}(0,0,h_{\mathrm{a}}))^{1/2}$,
as the velocity (and field; density) unit, and the resulting Alfvén time,
$\tau_{\mathrm{A}}=h_{\mathrm{a}}/V_{\mathrm{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_{\mathrm{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 Section 4.2. The grid of $251\times 180\times
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,
$\Delta_{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\,\mbox{:}\,-l_{x}\leq x\leq
l_{x},\,-l_{y}\leq y\leq l_{y},\,0\leq z\leq l_{z}\\}$ according to

$\displaystyle H=$ $\displaystyle\int_{\text{Box
volume}}\@vec{A}\cdot\@vec{B}\,\mathrm{d}^{3}\@vec{r}$
$\displaystyle+\int_{\text{Side
faces}}\phi_{\mathrm{C}}\@vec{A}_{\mathrm{C}}\,\mathrm{d}\@vec{S}+\int_{\text{Top
face}}(\@vec{A}\times\@vec{A}_{\mathrm{C}})\,\mathrm{d}\@vec{S}\,,$ (1)

where $\@vec{A}$ and $\@vec{A}_{\mathrm{C}}$ are special vector potentials for
$\@vec{B}$ and for the current-free field with the same normal component at
the box boundary, $\@vec{B}_{\mathrm{C}}$, and $\phi_{\mathrm{C}}$ is a scalar
potential for $\@vec{B}_{\mathrm{C}}$, that is

$\@vec{B}=\nabla\times\@vec{A}\,,\;\@vec{B}_{\mathrm{C}}=\nabla\times\@vec{A}_{\mathrm{C}}=-\nabla\phi_{\mathrm{C}}\,,$
(2)
$\Delta\phi_{\mathrm{C}}=0\,,\;\left.\nabla\phi_{\mathrm{C}}\cdot\@vec{\hat{n}}\right|_{\partial
V}=\left.-B_{n}\right|_{\partial V}\,,$ (3)

where $\@vec{\hat{n}}$ is the exterior unit normal at the boundary $\partial
V$ of $V$. The two vector potentials are given by

$\@vec{A}=\@vec{A}_{\mathrm{C}}(x,y,0)-\@vec{\hat{z}}\times\int_{0}^{z}\@vec{B}(x,y,z^{\prime})\,\mathrm{d}z^{\prime}$
(4)

and

$\@vec{A}_{\mathrm{C}}=\@vec{\hat{z}}\times\nabla\left(\int_{0}^{z}\phi_{\mathrm{C}}(x,y,z^{\prime})\,\mathrm{d}z^{\prime}+g_{\mathrm{C}}(x,y)\right)\,,$
(5)

where $\@vec{\hat{z}}$ is the unit vector in the vertical ($z$) direction and
$g_{\mathrm{C}}(x,y)$ is a solution of the two-dimensional Poisson equation

$\left(\partial_{x}^{2}+\partial_{y}^{2}\right)g_{\mathrm{C}}(x,y)=B_{z}(x,y,0)\,.$
(6)

Our procedure for calculating $H$ is an extension of that derived by DeVore
(2000) for the half space $z\geq 0$, where $H$ is given by the first line in
Eq. (1). The adaption of the method of DeVore, characterized in particular by
the gauge
$\@vec{A}\cdot\@vec{\hat{z}}=\@vec{A}_{\mathrm{C}}\cdot\@vec{\hat{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 $\@vec{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$,

$\@vec{B}_{\mathrm{cl}}=\@vec{B}-\@vec{B}_{\mathrm{C}},$ (7)

and its vector potential,
$\@vec{A}_{\mathrm{cl}}=\@vec{A}-\@vec{A}_{\mathrm{C}}$, the helicity can be
decomposed into self helicity of the closed field,

$H_{\mathrm{self}}=\int_{V}\@vec{A}_{\mathrm{cl}}\cdot\@vec{B}_{\mathrm{cl}}\,\mathrm{d}^{3}\@vec{r,}$
(8)

and mutual helicity between the closed and open potential field,

$H_{\mathrm{mutual}}=2\int_{V}\@vec{A}_{\mathrm{C}}\cdot\@vec{B}_{\mathrm{cl}}\,\mathrm{d}^{3}\@vec{r;}$
(9)

see Berger (1999). The relation $H=H_{\mathrm{self}}+H_{\mathrm{mutual}}$ is
satisfied by our numerically computed helicities to better than 0.3%. The
closed field $\@vec{B}_{\mathrm{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)=-d\ln
B_{\mathrm{ep}}(z)/d\ln z$, at the apex point of the flux rope’s magnetic
axis, which is slightly offset from the geometric axis, $h_{\mathrm{m}}=1.1$.
Experience with many different systems has shown that the threshold value lies
in the range $n_{\mathrm{cr}}\sim 1\dots 2$ and depends on various parameters
of the equilibrium (in ways yet to be determined). The canonical value of
$n_{\mathrm{cr}}\sim 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_{\mathrm{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
$|\Phi|=2.45\pi$, and sufficiently small for stability against the helical
kink mode (Hood & Priest, 1981; Török et al., 2004).

Figure 1: Rise profile of the flux rope apex in the reference run T1, showing
height, $h(t)$, and velocity, $\log u_{z}(t)$.

Figure 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, $B_{z}(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$.

Figure 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 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\approx 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\approx 50$.
The flux rope reaches a rise velocity of $u_{z}\approx 0.6$, and its apex
point leaves the box at $t\approx 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\geq 82$ in Figure 2), which is close to the potential field
(see, e.g., Figure 6(e) 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.

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\approx 127$,
and levels off after the top part of the flux rope has nearly completely left
the box at $t\approx 200$. A major part (64%) of the initial helicity is shed
by $t\approx 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).

The initial helicity, $H_{0}$, in this configuration is partitioned into self-
and mutual helicity as $H_{\mathrm{self,0}}=0.33H_{0}$ and
$H_{\mathrm{mutual,0}}=0.67H_{0}$. The initial self helicity due to the twist
in the initial flux rope completely leaves the system with the flux of the CME
bubble (Figure 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 (Figure 3 at $t=299$).

The computation of the helicity was checked by extending the integrals in
Equation (1) to a lower part of the box only, $z=0\mbox{--}10$. The initial
helicity was found to be identical to the one in Figure 4 (to within $1\times
10^{-3}$); 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 Figure 4 up to the
arrival of the top part of the flux rope at the virtual boundary at $t\approx
73$ and to show a slightly more rapid decrease by the same amount, $\Delta
H/F^{2}=0.13$, by $t\approx 155$, i.e., the same shedding effectiveness.

###  4.2 Parametric study

####  4.2.1 Effect of the guide field

Figure 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 $H_{0}/F^{2}$
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.1R_{\odot}$ 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,
$B_{\mathrm{et}}/B_{\mathrm{ep}}$, decay index, $n$, at the initial flux rope
apex, and twist averaged over the cross-section of the current channel,
$\Phi$. The range of the conjectured base level of the mutual helicity in the
configurations T1–T3, $\approx\\!0.06\mbox{--}0.08$, is shaded in light blue.
Table 1: Parameter values of the Runs T0–T4, K, and S.

Run | T0 | T1 | T2 | T3 | T4 | K | S
---|---|---|---|---|---|---|---
$B_{\mathrm{et}}/B_{\mathrm{ep}}$ | 0 | 0 | 0.2 | 0.5 | 0.5 | 0.5 | $\sim\\!0.4$
$d$ | 0.1 | 0.1 | 0.1 | 0.1 | 0.375 | 0.1 | $\dots$
$a$ | 0.9 | 0.6 | 0.6 | 0.6 | 0.6 | 0.27 | $\dots$
$L$ | 1.3 | 1.3 | 1.2 | 0.9 | 0.65 | 1.0 | $\dots$
$n(z\\!=\\!h_{\mathrm{m}})$ | 1.45 | 1.29 | 1.40 | 1.80 | 1.73 | 1.53 | 1.5–1.75
$\Phi/\pi$ | 1.7 | 2.5 | 2.3 | 2.1 | 2.4 | 4.9 | $\sim\\!0.5$
$H_{0}/F^{2}$ | 0.29 | 0.21 | 0.15 | 0.11 | 0.10 | 0.063 | 0.028
$H_{\mathrm{self,0}}/F^{2}$ | 0.084 ( 29%) | 0.068 ( 33%) | 0.044 ( 29%) | 0.022 ( 20%) | 0.020 ( 19%) | 0.010 ( 16%) | 0.008 (28%) aa$a$Percentage of $H_{0}/F^{2}$ is given in parentheses.
$H_{\mathrm{mutual,0}}/F^{2}$ | 0.21 ( 71%) | 0.14 ( 67%) | 0.11 ( 71%) | 0.086 ( 80%) | 0.084 ( 81%) | 0.053 ( 84%) | 0.021 (72%) aa$a$Percentage of $H_{0}/F^{2}$ is given in parentheses.
$\Delta H/F^{2}$ | 0.20 ( 70%) | 0.13 ( 64%) | 0.081 ( 54%) | 0.043 ( 39%) | 0.040 ( 38%) | 0.015 ( 24%) | 0.012 (42%) bb$b$Shedding effectiveness is given in parentheses as a percentage of the initial value.
$\Delta H_{\mathrm{self}}/F^{2}$ | 0.084 (100%) | 0.075 (100%) | 0.044 (100%) | 0.022 (100%) | 0.021 (100%) | 0.012 (100%) | 0.007 (87%) bb$b$Shedding effectiveness is given in parentheses as a percentage of the initial value.
$\Delta H_{\mathrm{mutual}}/F^{2}$ | 0.118 ( 57%) | 0.057 ( 41%) | 0.037 ( 34%) | 0.021 ( 24%) | 0.019 ( 23%) | 0.003 ( 6%) | 0.005 (25%) bb$b$Shedding effectiveness is given in parentheses as a percentage of the initial value.

111

The strong difference in the effectiveness of helicity shedding between the
torus-unstable configuration of Section 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_{\mathrm{et}}/B_{\mathrm{ep}}=0.2$ (Run T2) and
$B_{\mathrm{et}}/B_{\mathrm{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_{\mathrm{et}}/B_{\mathrm{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 Figure 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 $\approx\\!0.06\mbox{--}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 (Figure 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 $\Phi\approx 4\pi$ 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_{\mathrm{et}}/B_{\mathrm{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

Figure 5: Same as Figure 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.

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_{\mathrm{cc}}$, and external poloidal flux, $F_{\mathrm{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\mbox{--}1$. This
corresponds to a flux ratio
$\mathcal{R}=F_{\mathrm{cc}}/F_{\mathrm{ep}}=0.26\mbox{--}1.06$
($\mathcal{R}=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_{\mathrm{f}}=1.1$ in Run T1 to $D_{\mathrm{f}}=1.3$ and
1.6, respectively. Both sets of runs use $B_{\mathrm{et}}/B_{\mathrm{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_{0}/F^{2}=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 $\mathcal{R}$. Consequently, the
effectiveness of helicity shedding increases from $\Delta H/F^{2}=0.084$ (54%)
for $a=0.45$ to $\Delta H/F^{2}=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_{\mathrm{cr}}<1.35$ for $a=0.45\mbox{--}0.9$ and
within $1.35<L_{\mathrm{cr}}<1.4$ for $a=1$. The decay index values
$n(z=h_{\mathrm{m}})$ given in Figure 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).

The highest normalized helicity of a stable configuration in our parametric
study is $H_{0}/F^{2}=0.27$, found for $a=0.9$ with $d=0.1$ and $L=1.35$. When
analyzed with higher spatial resolutions of $\Delta_{x,y,z}=0.02$ and 0.01,
this value reduces slightly to $H_{0}/F^{2}=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 $\Delta_{x,y,z}=0.04$,
does not change the computed helicities.) The corresponding values of other
relevant parameters are $H_{\mathrm{self,0}}/H_{0}=0.315$ and
$\mathcal{R}=F_{\mathrm{cc}}/F_{\mathrm{ep}}=0.94$. To our knowledge, these
are the highest values of $H/F^{2}$ and $\mathcal{R}$ 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\sim 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_{\mathrm{ep}}$
links with the flux rope, increases with $D_{\mathrm{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_{\mathrm{et}}/B_{\mathrm{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
$\sim\\!0.06\mbox{--}0.08$ obtained in Runs T1–T3 (see Table 1 and Figure 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 $\approx\\!520$ km s-1, 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\approx
1.5\mbox{--}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 $\sim\\!0.5\pi$. Approximating the external field by the
potential field, $B_{\mathrm{et}}/B_{\mathrm{ep}}\sim 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.42H_{0}$, similar to Runs T2–T3
as well.

The brief increase of the helicity in this run during
$t\approx\\!100\mbox{--}130$ is a signature of the steepening of the current
layer in front of the erupting flux rope (see Figure 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 $\approx\\!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_{\mathrm{et}}/B_{\mathrm{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_{\mathrm{et}}/B_{\mathrm{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, $\sim\\!(0.4\mbox{--}0.65)H_{0}$.

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_{\mathrm{mutual}}/F^{2}\sim 0.06\mbox{--}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_{0}/2$. Therefore, further mechanisms
may be relevant in regulating the helicity budget of the corona (Section 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_{\mathrm{mutual,0}}/H_{\mathrm{self,0}}=2.1\dots 3.9$ (corresponding to
ratios of the equilibrium guide and strapping field components in the range
$B_{\mathrm{et}}/B_{\mathrm{ep}}=0\dots 0.5$), we find that a fraction of
$\sim\\!(0.4\mbox{--}0.65)H_{0}$ of the initial helicity $H_{0}$ is shed. A
data-constrained simulation of a solar eruption (from Active Region NOAA 11060
on 2010 April 8) yields a shedding of $0.42H_{0}$, 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 $\sim\\!0.06\mbox{--}0.08$ for our torus-unstable model configurations.

The parametric study yields stable equilibria up to a normalized helicity of
$H/F^{2}=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,
$\mathcal{R}=F_{\mathrm{cc}}/F_{\mathrm{ep}}=0.94$, is significantly higher
than previous estimates of the limiting flux ratio, which lie in the range of
$\sim\\!0.1\mbox{--}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_{\mathrm{self}}/H=0.315$, which is only slightly higher
than a previously suggested limit of $\simeq 0.29\pm 0.01$ for torus-unstable
flux ropes (Zuccarello et al., 2018), but lower than a limit of $\sim\\!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^{2}$ .

###### Acknowledgements.

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

  * Alexakis et al. (2006) Alexakis, A., Mininni, P. D., & Pouquet, A. 2006, ApJ, 640, 335
  * Aulanier & Dudík (2019) Aulanier, G. & Dudík, J. 2019, A&A, 621, A72
  * Aulanier et al. (2010) Aulanier, G., Török, T., Démoulin, P., & DeLuca, E. E. 2010, ApJ, 708, 314
  * Bateman (1978) Bateman, G. 1978, MHD Instabilities (MIT Press, Cambridge, MA)
  * Berger (1999) Berger, M. A. 1999, in Geophys. Monogr., Vol. 111, Magnetic Helicity in Space and Laboratory Plasmas, ed. M. R. Brown, R. C. Canfield, & A. A. Pevtsov, 1
  * Berger & Field (1984) Berger, M. A. & Field, G. B. 1984, J. Fluid Mechanics, 147, 133
  * Bobra et al. (2008) Bobra, M. G., van Ballegooijen, A. A., & DeLuca, E. E. 2008, ApJ, 672, 1209
  * Cheng et al. (2013) Cheng, X., Zhang, J., Ding, M. D., Liu, Y., & Poomvises, W. 2013, ApJ, 763, 43
  * Cheng et al. (2020) Cheng, X., Zhang, J., Kliem, B., et al. 2020, ApJ, 894, 85
  * Démoulin et al. (2002) Démoulin, P., Mandrini, C. H., van Driel-Gesztelyi, L., et al. 2002, A&A, 382, 650
  * DeVore (2000) DeVore, C. R. 2000, ApJ, 539, 944
  * Fan (2010) Fan, Y. 2010, ApJ, 719, 728
  * Finn & Antonsen (1985) Finn, J. M. & Antonsen, T. M. 1985, Comm. Plasma Phys. Contr. Fusion, 9, 111
  * Forbes & Isenberg (1991) Forbes, T. G. & Isenberg, P. A. 1991, ApJ, 373, 294
  * Green et al. (2002) Green, L. M., López Fuentes, M. C., Mandrini, C. H., et al. 2002, Sol. Phys., 208, 43
  * Hale (1925) Hale, G. E. 1925, PASP, 37, 268
  * Hassanin & Kliem (2016) Hassanin, A. & Kliem, B. 2016, ApJ, 832, 106
  * Hassanin et al. (2016) Hassanin, A., Kliem, B., & Seehafer, N. 2016, AN, 337, 1082
  * Hood & Priest (1981) Hood, A. W. & Priest, E. R. 1981, GAFD, 17, 297
  * Ji et al. (2003) Ji, H., Wang, H., Schmahl, E. J., Moon, Y.-J., & Jiang, Y. 2003, ApJ, 595, L135
  * Kliem et al. (2011) Kliem, B., Rust, S., & Seehafer, N. 2011, in Advances in Plasma Astrophysics, ed. A. Bonanno, E. de Gouveia Dal Pino, & A. G. Kosovichev, Vol. 274, 125
  * Kliem et al. (2013) Kliem, B., Su, Y. N., van Ballegooijen, A. A., & DeLuca, E. E. 2013, ApJ, 779, 129
  * Kliem & Török (2006) Kliem, B. & Török, T. 2006, Phys. Rev. Lett., 96, 255002
  * Kliem et al. (2012) Kliem, B., Török, T., & Thompson, W. T. 2012, Sol. Phys., 281, 137
  * Kusano et al. (2020) Kusano, K., Iju, T., Bamba, Y., & Inoue, S. 2020, Science, 369, 587
  * LaBonte et al. (2007) LaBonte, B. J., Georgoulis, M. K., & Rust, D. M. 2007, ApJ, 671, 955
  * Lin & Forbes (2000) Lin, J. & Forbes, T. G. 2000, J. Geophys. Res., 105, 2375
  * Low (1996) Low, B. C. 1996, Sol. Phys., 167, 217
  * Mackay & Yeates (2012) Mackay, D. H. & Yeates, A. R. 2012, Living Reviews in Solar Physics, 9, 6
  * Nindos et al. (2015) Nindos, A., Patsourakos, S., Vourlidas, A., & Tagikas, C. 2015, ApJ, 808, 117
  * Pariat et al. (2017) Pariat, E., Leake, J. E., Valori, G., et al. 2017, A&A, 601, A125
  * Parnell et al. (2009) Parnell, C. E., DeForest, C. E., Hagenaar, H. J., et al. 2009, ApJ, 698, 75
  * Patsourakos et al. (2013) Patsourakos, S., Vourlidas, A., & Stenborg, G. 2013, ApJ, 764, 125
  * Pevtsov et al. (1995) Pevtsov, A. A., Canfield, R. C., & Metcalf, T. R. 1995, ApJ, 440, L109
  * Qiu et al. (2007) Qiu, J., Hu, Q., Howard, T. A., & Yurchyshyn, V. B. 2007, ApJ, 659, 758
  * Reeves & Golub (2011) Reeves, K. K. & Golub, L. 2011, ApJ, 727, L52
  * Rudenko & Myshyakov (2011) Rudenko, G. V. & Myshyakov, I. I. 2011, Sol. Phys., 270, 165
  * Rust (1994) Rust, D. M. 1994, Geochim. Res. Lett., 21, 241
  * Sakurai (1976) Sakurai, T. 1976, PASJ, 28, 177
  * Sato & Hayashi (1979) Sato, T. & Hayashi, T. 1979, Physics of Fluids, 22, 1189
  * Savcheva & van Ballegooijen (2009) Savcheva, A. & van Ballegooijen, A. 2009, ApJ, 703, 1766
  * Savcheva et al. (2012) Savcheva, A. S., Green, L. M., van Ballegooijen, A. A., & DeLuca, E. E. 2012, ApJ, 759, 105
  * Schrijver et al. (2008) Schrijver, C. J., Elmore, C., Kliem, B., Török, T., & Title, A. M. 2008, ApJ, 674, 586
  * Seehafer (1990) Seehafer, N. 1990, Sol. Phys., 125, 219
  * Su et al. (2007) Su, Y., Golub, L., & Van Ballegooijen, A. A. 2007, ApJ, 655, 606
  * Su et al. (2011) Su, Y., Surges, V., van Ballegooijen, A., DeLuca, E., & Golub, L. 2011, ApJ, 734, 53
  * Su et al. (2006) Su, Y. N., Golub, L., van Ballegooijen, A. A., & Gros, M. 2006, Sol. Phys., 236, 325
  * Thalmann et al. (2011) Thalmann, J. K., Inhester, B., & Wiegelmann, T. 2011, Sol. Phys., 272, 243
  * Titov & Démoulin (1999) Titov, V. S. & Démoulin, P. 1999, A&A, 351, 707
  * Török & Kliem (2003) Török, T. & Kliem, B. 2003, A&A, 406, 1043
  * Török & Kliem (2005) Török, T. & Kliem, B. 2005, ApJ, 630, L97
  * Török & Kliem (2007) Török, T. & Kliem, B. 2007, AN, 328, 743
  * Török et al. (2004) Török, T., Kliem, B., & Titov, V. S. 2004, A&A, 413, L27
  * Valori et al. (2012) Valori, G., Démoulin, P., & Pariat, E. 2012, Sol. Phys., 278, 347
  * Valori et al. (2016) Valori, G., Pariat, E., Anfinogentov, S., et al. 2016, Space Sci. Rev., 201, 147
  * van Ballegooijen & Martens (1989) van Ballegooijen, A. A. & Martens, P. C. H. 1989, ApJ, 343, 971
  * van Tend & Kuperus (1978) van Tend, W. & Kuperus, M. 1978, Sol. Phys., 59, 115
  * Warren et al. (2011) Warren, H. P., O’Brien, C. M., & Sheeley, Neil R., J. 2011, ApJ, 742, 92
  * Yang et al. (2013) Yang, S., Büchner, J., Santos, J. C., & Zhang, H. 2013, Sol. Phys., 283, 369
  * Yang et al. (2018) Yang, S., Büchner, J., Skála, J., & Zhang, H. 2018, A&A, 613, A27
  * Zuccarello et al. (2015) Zuccarello, F. P., Aulanier, G., & Gilchrist, S. A. 2015, ApJ, 814, 126
  * Zuccarello et al. (2018) Zuccarello, F. P., Pariat, E., Valori, G., & Linan, L. 2018, ApJ, 863, 41

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

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

$\@vec{A}=\@vec{b}+\@vec{\hat{z}}\times\int_{z}^{z_{2}}\@vec{B}\,\mathrm{d}z^{\prime}\,,\quad\@vec{A}_{\mathrm{C}}=\@vec{b}_{\mathrm{C}}+\@vec{\hat{z}}\times\int_{z}^{z_{2}}\@vec{B}_{\mathrm{C}}\,\mathrm{d}z^{\prime}\,,$
(10)

where

$\@vec{b}=\begin{pmatrix}b_{x}(x,y)\\\ b_{y}(x,y)\\\
0\end{pmatrix}\,,\quad\@vec{b}_{\mathrm{C}}=\begin{pmatrix}b_{\mathrm{C},x}(x,y)\\\
b_{\mathrm{C},y}(x,y\\\ 0\end{pmatrix}$ (11)

satisfy

$\nabla\times\@vec{b}=\nabla\times\@vec{b}_{\mathrm{C}}=B_{z}(x,y,z_{2})\,\@vec{\hat{z}}\,.$
(12)

With the choice $z_{2}=0$ and $\@vec{b}=\@vec{b}_{\mathrm{C}}$, Eqs. (10)
become

$\@vec{A}=\@vec{b}_{\mathrm{C}}-\@vec{\hat{z}}\times\int_{0}^{z}\@vec{B}\,\mathrm{d}z^{\prime}\,,\quad\@vec{A}_{\mathrm{C}}=\@vec{b}_{\mathrm{C}}-\@vec{\hat{z}}\times\int_{0}^{z}\@vec{B}_{\mathrm{C}}\,\mathrm{d}z^{\prime}\,,$
(13)

where $\@vec{b}_{\mathrm{C}}$ obeys

$\nabla\times\@vec{b}_{\mathrm{C}}=B_{z}(x,y,0)\,\@vec{\hat{z}}\,.$ (14)

As

$\@vec{b}_{\mathrm{C}}(x,y)=\@vec{A}_{\mathrm{C}}(x,y,0)\,,$ (15)

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

As suggested by Valori et al., we choose $\@vec{b}_{\mathrm{C}}$ solenoidal,
which is possible because both the curl (given by Eq. (14)) and the divergence
of $\@vec{b}_{\mathrm{C}}$ can be prescribed. Therefore,
$\@vec{b}_{\mathrm{C}}$ can be represented in the form

$\@vec{b}_{\mathrm{C}}=\begin{pmatrix}-\partial_{y}g_{\mathrm{C}}\\\
\partial_{x}g_{\mathrm{C}}\\\ 0\end{pmatrix}=\@vec{\hat{z}}\times\nabla
g_{\mathrm{C}}$ (16)

by means of a stream function $g_{\mathrm{C}}(x,y)$. The second equation in
Eqs. (13) then becomes

$\@vec{A}_{\mathrm{C}}=\@vec{\hat{z}}\times\nabla
g_{\mathrm{C}}+\int_{0}^{z}\nabla\phi_{\mathrm{C}}\,\mathrm{d}z^{\prime}=\@vec{\hat{z}}\times\nabla\left(g_{\mathrm{C}}+\int_{0}^{z}\phi_{\mathrm{C}}\,\mathrm{d}z^{\prime}\right)\,,$
(17)

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

$\left(\partial_{x}^{2}+\partial_{y}^{2}\right)g_{\mathrm{C}}(x,y)=B_{z}(x,y,0)\,,$
(18)

which is identical to Eq. (6).

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

$\displaystyle H=$
$\displaystyle\int_{V}\left(\@vec{A}+\@vec{A}_{\mathrm{C}}\right)\cdot\left(\@vec{B}-\@vec{B}_{\mathrm{C}}\right)\,\mathrm{d}^{3}\@vec{r}=\int_{V}\@vec{A}\cdot\@vec{B}\,\mathrm{d}^{3}\@vec{r}$
$\displaystyle-\int_{V}\@vec{A}_{\mathrm{C}}\cdot\@vec{B}_{\mathrm{C}}\,\mathrm{d}^{3}\@vec{r}+\int_{V}\left(\@vec{A}_{\mathrm{C}}\cdot\@vec{B}-\@vec{A}\cdot\<EMAIL_ADDRESS>(19)

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

$\displaystyle\int_{V}\@vec{A}_{\mathrm{C}}\cdot\@vec{B}_{\mathrm{C}}\,\mathrm{d}^{3}\@vec{r}$
$\displaystyle=-\int_{V}\nabla\phi_{\mathrm{C}}\cdot\@vec{A}_{\mathrm{C}}\,\mathrm{d}^{3}\@vec{r}$
$\displaystyle=-\int_{V}\nabla\cdot\left(\phi_{\mathrm{C}}\@vec{A}_{\mathrm{C}}\right)\,\mathrm{d}^{3}\@vec{r}$
$\displaystyle=-\int_{\partial
V}\phi_{\mathrm{C}}\@vec{A}_{\mathrm{C}}\cdot\mathrm{d}\@vec{S}$
$\displaystyle=-\int_{\text{Side
faces}}\phi_{\mathrm{C}}\@vec{A}_{\mathrm{C}}\cdot\mathrm{d}\@vec{S}\,,$ (20)

where we make use of $\nabla\cdot\@vec{A}_{\mathrm{C}}=0$ and
$\@vec{A}_{\mathrm{C}}\cdot\@vec{\hat{z}}=0$. The second integral in the
second line of Eq. (19) can finally be rewritten as

$\displaystyle\int_{V}\left(\@vec{A}_{\mathrm{C}}\cdot\@vec{B}-\@vec{A}\cdot\@vec{B}_{\mathrm{C}}\right)\,\mathrm{d}^{3}\@vec{r}$
$\displaystyle=\int_{V}\nabla\cdot\left(\@vec{A}\times\@vec{A}_{\mathrm{C}}\right)\,\mathrm{d}^{3}\@vec{r}$
$\displaystyle=\int_{\partial
V}\left(\@vec{A}\times\@vec{A}_{\mathrm{C}}\right)\cdot\mathrm{d}\@vec{S}$
$\displaystyle=\int_{\text{Top
face}}\left(\@vec{A}\times\@vec{A}_{\mathrm{C}}\right)\cdot\mathrm{d}\@vec{S}\,,$
(21)

where we take advantage of the fact that $\@vec{A}\times\@vec{A}_{\mathrm{C}}$
is purely vertical (parallel to the $z$ axis) and
$\@vec{A}(x,y,0)=\@vec{A}_{\mathrm{C}}(x,y,0)$. Inserting Eqs. (20) and (21)
into Eq. (19) leads to Eq. (1).