© L. Linan et al. 2020

1. Introduction

Magnetic helicity is a volume-integrated ideal magnetohydrodynamic (MHD) invariant describing the level of twist and entanglement of the magnetic field lines. Initially introduced by Elsasser (1956), magnetic helicity is a conserved quantity within the ideal MHD paradigm (Woltjer 1958). However, the strict definition of Elsasser is gauge invariant only for magnetically bounded system, a condition that is not satisfied in most cases, such as the solar atmosphere. This led to the introduction by Berger & Field (1984) and Finn & Antonsen (1985) of the relative magnetic helicity, a gauge-invariant quantity suitable for use in solar physics and more generally for natural plasmas.

Using a numerical simulation, Pariat et al. (2015) confirmed the hypothesis introduced by Taylor (1974) that even in presence of nonideal processes, the dissipation of magnetic helicity is negligible. Relative magnetic helicity cannot be dissipated or created within the corona, therefore it can only be transported or annihilated. This conservation property has several major consequences, one of which might be that coronal mass ejections (CMEs) are the consequence of the evacuation of an excess of helicity (Rust 1994; Low 1996).

In recent years, magnetic helicity has been at the heart of many studies dealing with various topics such as the generation of solar eruptions (e.g., Kusano et al. 2004; Longcope & Beveridge 2007; Priest et al. 2016), magnetic reconnection (e.g., Linton et al. 2001; Linton & Antiochos 2002; Del Sordo et al. 2010), solar filaments (e.g., Antiochos 2013; Knizhnik et al. 2015; Zhao et al. 2015), and solar and stellar dynamos (e.g., Brandenburg & Subramanian 2005; Simon 2012).

Magnetic energy is another relevant quantity in MHD with which eruptivity in the solar corona is studied because most of solar events are driven magnetically (e.g., Schrijver & Zwaan 2008). From the point of view of the energetic budget, magnetic energy is the only source of energy that can generate the powerful events that are observed in the solar atmosphere, such as coronal mass ejections, flares, and solar jets (Forbes 2000). Magnetic energy can be decomposed into a current-carrying energy, known as free energy, and a potential energy (cf. Sect. 3.1). Solar flares and CMEs are characterized by a rapid change of the coronal magnetic field that does not change the radial component of the photospheric field. Because the potential field is determined by the radial magnetic field at the photosphere, only the free energy can therefore be converted into kinetic and thermal energies during fast coronal events (Aulanier et al. 2009; Karpen et al. 2012). The potential energy thus represents the lowest energy state of the magnetic field in the solar corona (e.g., Priest 2014).

An analysis of magnetic energies combined with a study of magnetic helicities appears a powerful tool for characterizing active regions and their evolutions toward eruptive events. However, measuring these quantities from observational data remains challenging. One possibility is to estimate the accumulation of magnetic helicity and energy in the solar corona by integrating their fluxes across the solar photosphere over time (Kusano et al. 2002; Nindos et al. 2003; Yamamoto et al. 2005; Yamamoto & Sakurai 2009). This method cannot trace the coronal evolution and requires high-cadence time-series magnetograms as well as the velocity fields on the photosphere. Because no direct observation of the photospheric velocity is available, it is obtained by inferring the magnetic field on the solar surface. Despite the progress made in deducing the velocity field (Kusano et al. 2002; Welsch et al. 2004; Longcope 2004) as well as further improvement on flux estimations (Pariat et al. 2005; Chae 2007; Liu & Schuck 2012, 2013; Dalmasse et al. 2014, 2018), the computation of magnetic energy and helicity fluxes remains very sensitive to the method that is used and to the quality of the observations. A different approach is to compute energy and helicity in coronal volumes. Because magnetic energy and helicity are volume integrals, properly computing them with this method requires the full 3D knowledge of the magnetic field in the volume that is studied. Currently, only 2D measurement on the solar surface are provided, therefore a 3D extrapolation of the magnetic field is a necessary step. The diverse methods based on the volume-integration approach for estimating the magnetic relative helicity were benchmarked in Valori et al. (2016). Different solar active regions have previously been investigated (Valori et al. 2013; Moraitis et al. 2014, 2019; Guo et al. 2017; Polito et al. 2017; Temmer et al. 2017; James et al. 2018; Thalmann et al. 2019).

In parallel, the properties of both helicity and energy are still being studied in solar-like parametric simulations. Berger (2003) introduced the decomposition of the relative magnetic helicity into two gauge-invariant components: a current-carrying helicity related to the current-carrying magnetic field, and a complementary volume-threading helicity. Pariat et al. (2017) followed and estimated these quantities in a set of seven simulations of the formation of solar active regions (Leake et al. 2013, 2014). The different simulations led to either stable or eruptive configurations. The authors found that it is possible to distinguish the two configurations by studying the ratio of the current-carrying helicity to the relative helicity. The ratio before the eruption indeed presents high values only in the eruptive case. To better understand the properties of the relative helicity decomposition, Linan et al. (2018) provided the first analytical formulae of the time-variation of nonpotential and volume-threading helicity. They also computed and followed them in two simulations of Leake et al. (2013, 2014) and in a simulation of the generation of a coronal jet (Pariat et al. 2005). They found that the current-carrying helicity does indeed not evolve as a result of boundary fluxes, but builds up through its exchange with the volume-threading helicity. The evolutions of the current-carrying and the volume-threading helicities are partially controlled by a transfer term that reflects the exchange between these two types of helicity. This exchange term dominates the dynamics of the current-carrying helicity at different instants of the simulations. This means that this helicity does not only evolve as a result of boundary fluxes. The eruption phases of these simulation follow the same dynamics: the current-carrying helicity is first transformed into the volume-threading helicity, and then the latter is ejected from the domain by boundary fluxes. Moreover, the transfer term is expressed as a volume integral: consequently, these two helicities are not classically conserved quantities in the sense that they cannot be independently expressed as a flux through the boundaries, even in ideal MHD, unlike the relative magnetic helicity. This finding strengthens the knowledge of the properties of nonpotential and volume-threading helicity that was first studied by Moraitis et al. (2014).

Zuccarello et al. (2015) presented 3D parametric resistive MHD simulations of solar coronal eruptions. Simulations are distinguished by the different motions (line-tied) applied on the photosphere with similar but distinct flux cancelation drivers. Their eruptions were driven by the torus instability (Aulanier et al. 2009; Démoulin & Aulanier 2010) and occurred at a precise time identified by a series a relaxation runs. Recently, these models were used to investigate the increase in the downward component of the Lorentz force density around an polarity-inversion line in comparison with the photospheric observation (Barczynski et al. 2019). From these simulation, Zuccarello et al. (2018) studied the impact of the different boundary driving flows on the helicity and energy injection. They found that the helicity ratio of the current-carrying helicity to the relative helicity is clearly associated with the eruption trigger because the eruptions within the different runs took place exactly when the ratio reached the very same threshold value.

Recently, the first preliminary observational tests confirmed the idea that the helicity ratio is a good predictor of eruptivity. Based on 3D extrapolation and using different nonlinear force-free models, the time evolution of the helicity ratio has been investigated in three active regions: AR 12673 in Moraitis et al. (2019), the most active of the cycle 24; and AR 11158 and AR 12192 in Thalmann et al. (2019), two extensively studied active regions that generated eruptive and confined flares. However, complementary studies are still needed to understand how the different magnetic topologies observed in the solar corona are linked to the dynamics of the helicity ratio.

In the present study, we apply the helicity decomposition to the analysis of the simulations of Zuccarello et al. (2015, 2018) to investigate the time-variations of the different types of magnetic energy and helicity. In particular, we are interested in the way that the different boundary motion influence the helicity and energy dynamics. Zuccarello et al. (2018) showed that the different boundary motions lead to different efficiency in injecting helicity and energy in the domain. In the present work, we aim to explain the physical processes that are responsible for these differences: are they related to boundary fluxes, dissipation, or volume evolution? We also examine whether the transfer term between the two helicity components plays a major role in the helicity budgets, as has been observed in Linan et al. (2018).

Additionally, we study the dynamics of the helicities in comparison with the dynamics of their energy counterparts, for instance, current-carrying helicity and free energy, and volume-threading helicity and potential energy. Our goal is to highlight the differences and the similarities in the helicity and energy buildup. This study aims to improve our knowledge on magnetic helicities and energies, and it is a necessary step to better understand the full topological and energetic complexity of solar active regions.

Our paper is divided into different sections that are organized as follows. First, we present the time-variation of nonpotential and volume-threading helicities (see Sect. 2). In the same way, we then introduce the different components of the magnetic energy and also their time-variation written for the specific case of resistive MHD (see Sect. 3). After presenting the simulations (see Sect. 4), we present the time evolution of the different quantities (see Sect. 5). Using our set of simulations, we investigate the role of the transfer term between the helicity components in their evolutions (see Sect. 6). While we investigate the difference in helicity dynamics in Sect. 7, we focus on the similarities between magnetic energy and helicity fluxes in Sect. 8. In the conclusion, we discuss the effect of the different boundary-driving motions on the energy and helicity injection and on the eruptivity helicity ratio.

2. Nonpotential and volume-threading helicities

In the fixed volume V bounded by the surface S, the magnetic helicity is defined as

(1)

with A the vector potential of the studied magnetic field B, i.e ∇ × A = B. In practice, this scalar description of the geometrical properties of magnetic field lines is general only if the magnetic field is tangential to the surface, that is, if V is a magnetically bounded volume. The magnetic helicity is gauge invariant if and only if this condition is respected. For the study of natural plasmas, especially in solar physics, the magnetic field does not satisfy this condition, the solar photosphere being subject to significant flux. Berger & Field (1984) introduced the relative magnetic helicity, a gauge-invariant quantity, based on a reference field. Throughout the paper, we use the potential reference field B_p that has the same normal distribution of B throughout the surface S and satisfies

(2)

where n is the outward-pointing unit vector normal on S. The potential field can thus be defined by a scalar function, such as ∇ϕ = B_p, and ϕ is the solution of the Laplace equation,

(3)

When A_p is the vector potential of the potential field B_p  =  ∇×A_p, the relative magnetic helicity provided by Finn & Antonsen (1985) is defined as

(4)

In this form, the relative magnetic helicity is independently invariant to gauge transformation of both A and A_p. The difference between the potential field and the magnetic field can be written as a nonpotential magnetic field, B_j = B − B_p, associated with the vector A_j, defined as A_j = A − A_p, such as ∇ × A_j = B_j. When we use this unique decomposition of B and following the work of Berger (2003), Eq. (4) can be divided into two gauge-invariant quantities:

(5)

(6)

(7)

where H_j is the current-carrying magnetic helicity associated with only the current-carrying component of the magnetic field B_j, and H_pj is the volume-threading helicity involving both B and B_p. By construction, both H_j and H_pj are gauge invariant because by virtue of Eq. (3), B_j has a vanishing normal component on the surface.

In resistive MHD, where E = −v × B + η∇×B (η being the resistivity, which is here assumed to be constant), Linan et al. (2018) established the following equation for the time evolution of the current-carrying magnetic helicity H_j:

(8)

with

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

From this decomposition, Linan et al. (2018) obtained an equation for the time-variation that is composed only of gauge-invariant terms:

(17)

with

(18)

All these terms initially provided by Linan et al. (2018) are recalled here because we analyze them and comment on them in the next sections.

Similarly, the time evolution of the volume-threading magnetic helicity H_pj can be decomposed as

(19)

with

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

The time-variation of H_pj can also be constructed with gauge-invariant terms only:

(28)

with

(29)

These decompositions were obtained without any particular hypothesis on the gauge that is used. In particular, we are free to use the Coulomb gauge for A and A_p. With this choice, the volume terms dH_j/dt|_Bp,  var and dH_pj/dt|_Bp,  var both vanish. Thus dH_j/dt|_Own and dH_pj/dt|_Own only contain boundary-flux contributions. The transfer term dH_j/dt|_Transf expresses the exchange between the helicities H_j and H_pj without any consequence on the evolution of the total relative helicity H_V. Furthermore, this quantity being a volume term, the time-variations of H_j and H_pj cannot be expressed solely through boundary fluxes. Therefore H_j and H_pj are not conserved quantities in resistive or ideal MHD.

3. Magnetic energy

3.1. Free and potential energies

With the decomposition of the magnetic field into current-carrying and potential components, B = B_p + B_j, the total magnetic energy E_v can be also decomposed as

(30)

with

(31)

In this decomposition, E_j is the energy of the current-carrying magnetic field, also known as free energy, and E_p the energy of the solenoidal magnetic field. Because the potential field shares the same surface distribution as the total magnetic field, the surface integral vanishes in Eq. (30). Numerically, the discretization of the mesh grid unavoidably induces a finite level of non-solenoidality (∇ ⋅ B ≠ 0), and consequently, the last term in Eq. (30) is not exactly null (cf. Valori et al. 2013). However, considering a solenoidal field, Eq. (30) can be simplified into

(32)

This decomposition is similar to the decomposition of the helicity obtained in Eq. (5). However, here, the potential energy E_p, unlike the volume-threading helicity, only depends on the potential field without a dependence on the nonpotential field.

3.2. Time-variation of the total magnetic energy

We aim to determine the time-variation of the total magnetic energy E_v in a fixed volume V,

(33)

In the resistive MHD, we use the Faraday law, ∂B/∂t = −∇×E, and we obtain

(34)

Using the Gauss divergence theorem, we can decompose the first term of Eq. (34):

(35)

Here, the surface term corresponds to the surface integral of the poynting vector and can be divided into two terms:

(36)

Finally, assuming for simplicity that the resistivity is constant in space, the variation of the total magnetic energy can be decomposed as

(37)

with

(38)

(39)

(40)

(41)

We performed a similar time variation for E_v as Linan et al. (2018) did for helicities as we summarized in Sect. 2. We find two fluxes: F_Vn,  E a shearing term associated with horizontal motion, and an emerging term F_Bn,  E that is related to the emergence. Longcope et al. (2007) decomposed the shearing term at the photospheric level into two contributions by differentiating the motion between the different flux patches and the spin motion of isolated flux patches. As with H_j and H_pj, the time-variation of E_v cannot be expressed through boundary fluxes, and thus E_v is not a conserved quantity. Even in ideal MHD, when dE_v/dt|_Diss is null, a volume term dE_v/dt|_var survives. Unlike dH_j/dt and dH_pj/dt, the time-variation of the total energy E_v depends on neither A nor A_p: each term of Eq. (37) is gauge invariant.

3.3. Time-variation of the potential and free magnetic energies

Similarly to the analysis of dE_v/dt in the previous section, it is possible to obtain the time-variation of E_p. Using the scalar potential ϕ of B_p such as ∇ϕ = B_p and the Gauss divergence theorem, we write

(42)

with

(43)

(44)

For a purely solenoidal potential field, dE_p/dt|_ns is null and thus the time-variation of E_p is written in a simple way as a single surface term depending on the variation of the potential field. A decrease in potential energy can therefore be associated in particular with a cancelation at the polarity-inversion line (PIL) or a dispersion of the potential magnetic field. Unlike E_v, E_p is a conserved quantity in resistive MHD.

For the study of the evolution of nonpotential energy E_j, the easiest way is to consider only the difference between the decompositions obtained from Eqs. (37) and (42):

(45)

Still, this identity is truly accurate only in the case of purely solenoidal magnetic fields, that is, ∇ ⋅ B_j = 0 and when B_j|_S = 0 (cf. Sect. 3.1). Hereafter, we therefore introduce a generic formulation where terms that explicitly account for nonsolenoidal errors are retained,

(46)

The last two terms on the right are very similar to dE_p/dt and dE_v/dt. They can thus be decomposed in the same way,

(47)

and

By definition, the curl of the potential field is null, and thus the second volume integral on the right-hand side formally vanishes. Finally, by grouping all the terms, the time-variation of the nonpotential magnetic energy can be written as

(48)

with

(49)

(50)

(51)

(52)

(53)

(54)

As expected, the decomposition we obtain is similar to the one obtained with dE_v/dt (see Eq. (37)). The dissipation term dE_j/dt|_ns is null if the magnetic field is solenoidal. However, we compute it in order to quantify the effect of purely numerical errors.

4. Line-tied eruptive simulations

In order to comparatively analyze the evolution of helicities and energies and to study the time-variation of the energy, we used magnetic field data produced by parametric 3D MHD simulations of eruptive events of the solar corona that were initially presented in Zuccarello et al. (2015).

For this set of simulations, the OHM-MPI code (Aulanier et al. 2005) solves MHD equations in the system’s nondimensional units for a volume covering the domain x ∈ [ − 10, 10], y ∈ [ − 10, 10] and z ∈ [0, 30]. The employed mesh is nonuniform and the smallest cell is centered at x = y = z = 0. In order to facilitate the computation of the energies and helicities, our study was performed on a subdomain excluding the z = 0 plane, interpolated into a uniform Cartesian grid composed of 333 cells in the x and y direction, and 500 cells in the z directions. The resulting analyzed volume is x ∈ [ − 10, 10], y ∈ [ − 10, 10] and z ∈ [0.006, 30]. As a result of the interpolation on a uniform grid (whose cell sizes are 0.06, to be compared to the original grid, whose cell sizes range from 0.006 to 0.32), the magnetic field we obtained has a lower solenoidality than the initial grid. This reduces the accuracy of the magnetic helicity and energy computations (Valori et al. 2016).

The system is delimited by a set of boundaries subject to “open” boundary conditions (except at z = 0). In the analyzed datasets, the bottom boundary is at z = 0.006, one mesh point above the surface corresponding to the photospheric level where line-tied boundary conditions were imposed in the original numerical experiments. All the physical MHD quantities, such as the magnetic field, can leave the simulation domain through lateral and top boundaries during the evolution of the system.

Four parametric simulations were performed, all starting with a common phase that is referred to as the “shearing phase” in Zuccarello et al. (2015, 2018). Initially, the magnetic field is potential and generated by two unbalanced and asymmetric subphotospheric polarities. During the shearing phase, asymmetric vortices centered around the local maxima of |B_z(z = 0)| slowly evolve the initial potential magnetic field into a current-carrying magnetic field. This shearing flow motion induces a magnetic shear close to the PIL and at the end of this phase, creates a current-carrying magnetic field arcade surrounded by a quasi-potential background field. During this entire phase, the distribution of B_z at the bottom boundary remains unchanged. This phase lasts from the time t = 10t_A until t = 100t_A in the system time coordinate, where t_A is the reference Alfvén time used in Zuccarello et al. (2015).

Then the four parametric simulations differ by the motion pattern that is imposed at the bottom boundary (cf. Fig. 1). We refer to this phase as the “pre-eruption” phase. In a first motion profile, labeled “convergence” (Fig. 1, left panel), the velocity flows only follow the horizontal direction x and are only applied close to the PIL. This creates a cancellation of the magnetic flux around the PIL but only slightly affects the periphery of the active regions. Unlike the previous case, for the run labeled “stretching” (Fig. 1, middle left panel), these horizontal motions are also applied at the periphery of the active region. For the other two runs, labeled “dispersion central” and “dispersion peripheral” (Fig. 1, middle right and right panels), the motions spread in all directions from the center. The only difference is in the portion of the active region that is subjected to these motions. In the dispersion peripheral run, only the periphery of the active region is concerned, while in dispersion central, the dispersion also occurs in the center of the polarity where the magnetic field is strongest.

Fig. 1. Applied boundary-driving motions for the four different numerical experiments. White represents the positive polarity (Bz(z = 0.006) > 0) and black the negative polarity (Bz(z = 0.006) < 0). Orange and cyan arrows indicate the distribution of the velocity flows we applied to the negative and positive polarity, respectively.

The four runs all present a cancellation of magnetic flux at the PIL that is permitted by a finite photospheric diffusion. The sheared-arcade configuration at the end of the shearing phase evolves into a bald-patch topology, and the magnetic reconnection process leads to the formation of a flux rope. The system then evolves until it reaches the instant where it becomes unstable and erupts (cf. Aulanier et al. 2009, for a description of the eruption process). The onset of the eruption, that is, the time t₁ of the onset of the instability, is accurately determined by a series of relaxation runs for each simulation (Zuccarello et al. 2015). It occurs at t₁ = 196, 214, 220 and 164t_A for the convergence, stretching, dispersion peripheral, and dispersion central runs, respectively.

To ensure the numerical stability of the code, a finite resistivity η and a pseudo-viscosity ν are necessary (Zuccarello et al. 2015). During the common shearing phase until t <  100t_A and during the pre-eruption phase, the coronal diffusivity is η_cor = 4, 8 × 10⁻⁴ and the pseudo-viscosity is fixed to ν′ = 25. After the eruption, during a phase referred to as eruption phase, η_cor is 2, 1 × 10⁻³ and ν′ is 41, 7. To allow later flux cancellation at the PIL, a photospheric resistivity η_phot = η_cor = 4, 8 × 10⁻⁴ is imposed only during the pre-eruption phase. The photospheric resistivity is set to zero in the shearing phase and before the eruption. The change in resistivity follows a ramp-down time profile during the time t₁ − 5t_A <  t <  t₁ + 5t_A. This transitional period is called “eruption onset phase” and is represented as the yellow band in all the figures. We note that the time t₁ corresponds to the middle of the ramp-down time profile, therefore the boundary flows are null only at t >  t₁ + 5t_A. In this paper, we removed the first mesh point in the z-direction, which corresponds to the bottom boundary level. We therefore have a uniform resistivity throughout the domain in order to facilitate the calculation of the so-called “nonideal” terms.

5. Energy and helicity evolutions

In this section we first introduce the method for numerically computing the energies and helicities in our set of simulations. Then, we discuss the computation of the time-variations.

5.1. Energy and helicity estimations

In order to compute the different helicities and energies at each time t_A, we used the method of Valori et al. (2012). The datacubes of the magnetic field B, of the plasma-velocity field v, and of the plasma thermodynamic quantities allowed us to compute all the quantities that appear in Eqs. (8), (19), (48), and (42).

First the scalar potential ϕ(t) of the potential magnetic field B_p(t) was obtained from a numerical solution of the Laplace Eq. (3). The numerical methods we employed to solve this equation required an uniform grid and thus led to the interpolation of the initial grid, as mentioned in Sect. 4. The potential vectors A(t) and A_p(t) were computed according to Eq. (14) in Valori et al. (2012) and follow the DeVore-Coulomb gauge defined in Pariat et al. (2015):

(55)

(56)

This choice of gauge was complemented by the following relationship inherent in the integration method:

(57)

where ⊥ means the normal component. It corresponds to the 1D integration of magnetic fields starting at the top boundary of the domain at height z_top. Finally, we obtained the helicities and energies from Eqs. (5) and (30). In Fig. 2 we plot the time evolution of these quantities for the different runs. In order to facilitate the comparison between the different runs, the time is plotted with a modified time variable t − t₁ in each figure, where t₁ is the onset time defined in Sect. (4) and is different for each of the four simulations.

Fig. 2. Evolution of the different magnetic helicities (top panels), from left to right: relative magnetic helicity (Hv, Eq. (4)), volume-threading helicity (Hpj, Eq. (7)), and current-carrying helicity (Hj, Eq. (6)). Time evolution of the different magnetic energies (bottom panel), from left to right: total magnetic energy (Ev, Eq. (32)), potential energy (Ep, Eq. (31)), and free energy (Ej, Eq. (31)). The different simulations are dispersion central (red line), dispersion peripheral (green line), stretching (yellow line) and convergence (blue line). The yellow vertical band corresponds to the onset phase of the eruption.

Unlike Zuccarello et al. (2018), we are interested here in the evolution of the quantities after the common shearing phase, including the eruption phase (which was not studied by Zuccarello et al. 2018). We also recall that the domain studied here is slightly different from the one studied in Zuccarello et al. (2018) (cf. Sect. 4). Figure 2 shows that the dynamics of energies and helicities are qualitatively similar from one simulation to the other during the three different phases. The different boundary-forcing mainly affects the magnitude of the different quantities, but not the quality of their dynamical behaviors.

We also note that the dynamic of helicities and energies changes during the eruption. For the current-carrying helicity, H_j, and the free energy, E_j, we observe an overall increase during the pre-eruption phase, followed by a decrease during the eruption phase. The quantities related to the potential magnetic fields, E_p and H_pj, both decrease in the pre-eruption phase. In the eruption phase E_p remains constant while H_pj continues to decrease, although at a lower rate than in the pre-eruption phase.

Because both H_j and H_pj decrease, the relative helicity H_v also decreases in the eruption phase. The total energy, E_v, has a dynamics similar to H_v: the system loses energy throughout the simulation, but at a different rate before and after the eruption.

Overall, the behavior of the helicities is similar to their energy counterparts, for example, H_v to E_v, H_j to E_j, H_pj to E_p. This is particularly visible for the current-carrying helicity and the free energy: when H_j increases (decreases), E_j also increases (decreases). In the next sections we focus more on the reasons of these trends by studying the fluxes of the different quantities.

5.2. Time-variation estimation

After all the vectors were calculated, we computed the instantaneous time-variations of energies and helicities obtained from Eqs. (8), (19), (48), and (42). The surface integrals were calculated systematically as the sum of the contributions from the six boundaries. A study focusing on the contribution of the lower boundary alone is conducted in Sect. 9.3.

Linan et al. (2018) validated the time-variations equations established for the volume-threading helicity, H_pj, and the current-carrying helicity, H_j. The accuracy of the computation is related to difference factors such as the discretization and the remapping of the data, spatial and temporal, as well as to nonexplicit numerical diffusive terms that are not accounted for in our analytical resistive MHD model. In our study, we present a complementary test by comparing the time-variation of the relative helicity, dH_v/dt computed from Eq. (23) in Pariat et al. (2015), with the sum of the time-variations of nonpotential and volume-threading helicities, dH_j/dt + dH_pj/dt from (8) to (19). In this way, both sides are computed with the same temporal accuracy.

The result of this comparison is presented in the right panel of Fig. 3. In this figure we plot dH_v/dt and dH_j/dt + dH_pj/dt for the dispersion peripheral simulation. For the other runs, the difference is on the same order of magnitude and varies in a similar way. We therefore do not plot this here. The difference is very low, with an average deviation smaller than 0.1%. This confirms the robustness of our calculation method as well as the validity of our analytical equations.

Fig. 3. Left panel: time evolution of the instantaneous time-variation, dEv/dt (dashed black curves, Eq. (37)), and of the sum, dEj/dt + dEp/dt (continuous blue curves, Eqs. (8) and (19)) for the dispersion peripheral run. Right panel: time evolution of the instantaneous time-variation, dHv/dt (dashed black curves, Eq. (23) in Pariat et al. 2015, and of the sum, dHj/dt + dHpj/dt (continuous blue curves, Eqs. (8) and (19)) for the dispersion peripheral run. The yellow bands correspond to the eruption onset phase.

In the same way, we compare in Fig. 3 (left panel) the time-variation of the total magnetic energy, dE_v/dt, computed from Eq. (37), with the sum of the time-variations of free and potential energies, dE_j/dt + dE_p/dt, from Eqs. (42) to (48). Here, the difference is not negligible, with an average relative difference of 27%. The cause of this difference is likely mainly the non-solenoidality of the magnetic field. As mentioned Sect. 3, with a finite level of non-solenoidality, the equality of Eq. (45) is not fully accurate because E_v ≠ E_j,s + E_p (see Sect. 3.1).

In order to estimate the artificial non-solenoidal energy contributions, Valori et al. (2013) introduced the following decomposition:

(58)

where E_j,s and E_p, s are the energies of the current-carrying and potential solenoidal magnetic field. E_j,ns and E_p, ns are the nonsolenoidal components, whereas E_mix corresponds to all the remaining cross terms. For a solenoidal field we have E_j,ns = E_p, ns = E_mix = 0, E_j,ns = E_j, E_p, ns = E_p, and therefore dE_j/dt + dE_p/dt = dE_v/dt.

The finite non-solenoidality (∇ ⋅ B ≠ 0) affects the precision of the helicity and energy computations, as studied by Valori et al. (2016). Following Valori et al. (2013, 2016), Thalmann et al. (2019), we used the energy criteria E_div/E_v to quantify the non-solenoidality effect. The divergence-based energy is defined as

(59)

For the four simulation analyzed in this paper we find an average of E_div/E_v ≃ 2%. According to Valori et al. (2016), this value for the average of non-solenoidality leads to a precision of ≤6% for our helicity computations, which is much lower than the 27% discrepancy found in Fig. 3. One possible cause of non-solenoidality is our interpolation of the original data from a highly nonuniform mesh onto a uniform grid, which can increase the nondivergence of the magnetic field. Different tests have been made to degrade and also improve the interpolation to give a rough estimate of the effect. Finer interpolations, however, have required considering only a fraction of the whole numerical domain to keep the number of grid points manageable. The outcome of these tests is that neither presented results that differed significantly from our baseline interpolation. In particular, E_div/E_v always remained on the order of 2%, which means that this error therefore seems intrinsic to the numerical models. The level of interpolation chosen in our study therefore is a good compromise between the required computed power and the quality of our data. In addition, it is worth noting that various terms in our equations depend on time variations, so that the accuracy of our flux computations can also be limited by a relative coarseness of the time outputs of the available data. Testing for this would require recalculating the simulations with a higher cadence for its outputs, which is beyond the scope of this paper.

6. Helicity transfer

As mentioned in the introduction (see Sect. 1), Linan et al. (2018) showed for two different eruptive simulations that the exchange between H_pj and H_j is controlled by the gauge-invariant transfer term dH_j/dt|_Transf (see Eq. (10)), which therefore plays a key role in evolving these helicities. In order to confirm these different results in the particular case of our line-tied simulations, we plot in Fig. 4 the gauge-invariant terms of dH_j/dt (see Eq. (17)), and dH_pj/dt (see Eq. (28)) for the convergence and the dispersion peripheral runs.

Fig. 4. Time evolution of the helicity variation rates, dHj/dt and dHpj/dt (dashed black curves; Eqs. (8) and (19)), of the helicity transfer term, dHj/dt|Transf and dHpj/dt|Transf (solid red curves; Eqs. (10) and (21)), of the own terms, dHj/dt|Own and dHpj/dt|Own (solid blue curves; Eqs. (18) and (29)), and of the dissipation terms, dHj/dt|Diss and dHpj/dt|Diss (solid green curves; Eqs. (9) and (20)) for the convergence simulation (top panels) and for the dispersion peripheral simulation (bottom panels). Left and right panels: evolution of the current carrying helicity, Hj, and volume-threading helicity Hpj, respectively.

In the convergence case (see Fig. 4, top left panel), the conversion of helicity from H_pj to H_j and the boundary flux have similar amplitudes during the pre-eruption phase. Both contributions are positive, thus H_j increases (e.g., dH_j/dt is positive). The transfer term, dH_j/dt|_Transf, dominates the evolution of H_j when it is close in time to the eruption. Meanwhile, H_pj decreases mainly because of the strong helicity transfer (cf. Fig. 4, top right panel), that is, dH_j/dt|_Transf is the dominating term of dH_pj/dt. In comparison, the flux of H_pj related to the own term, dH_pj/dt|_Own, is very low. This means that the injection of H_pj is not enough to compensate for its conversion to H_j.

Moreover, as with the flux emergence simulations presented in Linan et al. (2018), the eruption is accompanied by a sharp decrease of the transfer term. However, unlike in Linan et al. (2018), the sign of dH_j/dt|_Transf does not change here during the onset phase of the eruption. During the eruption phase, the dissipation terms, dH_j/dt|_Diss and dH_pj/dt|_Diss, dominate the variation of the helicities. This is related to the increase in resistivity η imposed in the numerical experiment during that period.

In the dispersion peripheral simulation (see Fig. 4, bottom panels), the variations of H_j and H_pj during the pre-eruption phase are noticeably dominated by the boundary fluxes, dH_j/dt|_Own and dH_pj/dt|_Own. The transfer terms are significantly less intense than the injection of H_j and H_pj, except close to the eruption. During the eruption phase, the dynamics is mostly dominated by the resistive dissipation terms, dH_j/dt|_Diss and dH_pj/dt|_Diss.

We thus observe that while the trends of H_j and H_pj are similar for the convergence and the dispersion peripheral simulations, as discussed in Sect. 5, these variations are in fact due to noticeably different dynamics of the helicity fluxes. For instance, the decrease of H_pj in the pre-eruption phase is due to an intense conversion of helicity (high values of dH_pj/dt|_Transf) for the dispersion peripheral simulation, whereas in the convergence run, a similar evolution of H_pj is explained by an intense negative boundary flux, dH_pj/dt|_Own during that period.

Finally, as has been noted in Linan et al. (2018), for both the simulations analyzed here but also for the two others, the transfer terms cannot be neglected. We confirm that the estimations of boundary fluxes are not sufficient to follow the dynamics of H_j and H_pj. However, unlike with the flux emergence and solar jet simulations studied in Linan et al. (2018), the precise mechanism of the buildup of H_j and H_pj does depend on the simulations during the pre-eruption phase. Studying this dependence is the goal of the next section.

7. Distinguishing between simulations in terms of helicity dynamics

In the previous section, we described that the magnitude of the different gauge-invariant helicity variation terms could be significantly different in the dispersion peripheral and in the convergence simulations. This demonstrates that even if the general trends of H_pj and H_j are similar (see Sect. 5), their dynamics can be significantly different.

In order to estimate the effect of the different boundary driver, Fig. 5 displays the different gauge-invariant variation terms for the four different numerical simulations: the boundary fluxes dH_j/dt|_Own and dH_pj/dt|_Own; the transfer term, dH_j/dt|_Transf; and the dissipations terms dH_j/dt|_Diss and dH_pj/dt|_Diss.

Fig. 5. Time evolution of the different gauge invariant terms of dHj/dt (top panels), from left to right: dissipation term (dHj/dt|Diss, Eq. (9)), own term (dHj/dt|Own, Eq. (18)), and helicity transfer term (dHj/dt|Transf, Eq. (10)). Time evolution of the different gauge invariant terms of dHpj/dt (bottom panels), from left to right: dissipation term (dHpj/dt|Diss, Eq. (20)), and own term (dHpj/dt|Own, Eq. (29)). Each curve color corresponds to a particular simulation: dispersion central (red line), dispersion peripheral (green line), stretching (yellow line), and convergence (blue line). The yellow band corresponds to the onset phase of the eruption.

The dissipations terms (cf. Fig. 5, left panels) are not significantly different from one simulation to the other. The sudden increase in absolute values of the dissipations terms, observed during the eruption onset phase, is related to the imposed numerical increase in resistivity. The variations in magnitude, particularly in the pre-eruption phase, are minor compared to the variations in other gauge-invariant terms.

The boundary flux term dH_j/dt|_Own is also very similar from one simulation to another, except for the dispersion central run, for which more nonpotential helicity is markedly injected during the pre-eruption phase (see Fig. 5, middle top panel). Unlike dH_j/dt|_Own, the boundary flux of H_pj, dH_pj/dt|_Own, is strongly sensitive to the boundary-driving pattern (see Fig. 5, bottom left panel). Tthe sign and magnitude of dH_pj/dt|_Own depend on the simulation. For the dispersion simulations, there is a significant injection of negative H_pj, whereas in the convergence and stretching case, the flux is significantly weaker, if not of the opposite sign.

Finally, Fig. 5 (top right panel) shows that the helicity transfer rate, dH_j/dt|_Transf, is higher for the convergence and stretching simulations than for the dispersion cases in the pre-eruption phase. Unlike with the other cases where the transfer term is almost constant during an early period, in the dispersion central run dH_j/dt|_Transf increases from the first moments of the simulation.

In summary, we observed that during the pre-eruption phase, the increase of H_j and reciprocally the decrease of H_pj (cf. Fig. 5) are not explained by the same physical process in the different simulations. We observe three significantly different dynamics:

• The convergence and stretching simulations present a similar dynamics for their fluxes of helicity. They are characterized by a relatively weak boundary flux of H_pj counterbalanced by a strong transfer from H_pj to H_j. The own term of H_j is positive, involving an injection of current-carrying helicity. Its magnitude is almost identical in these two runs.

• For the dispersion peripheral run, H_pj decreases mostly because of the boundary flux, unlike with the previous cases. In comparison to the boundary flux, the transfer from H_pj to H_j is less important. The flux of H_j is similar to the convergence and stretching runs.

• The dispersion central shares some similarities with the dispersion peripheral run regarding the variations of H_pj. However, this simulation is characterized by a high boundary flux of H_j that is distinct and significantly higher than the three other cases.

Finally, simulations with the largest injection of helicities due to their own terms (whether H_j or H_pj) have the lowest magnitude of the transfer term. Inversely, a strong exchange between H_j and H_pj is accompanied by lower fluxes through the surfaces. Both lead to a similar trend for the relative helicity H_v. This shows that the boundary fluxes of H_j or H_pj as well as the volume term, dH_j/dt|_Transf, are directly related to the morphology and the evolution of the magnetic field at the bottom boundary. In Sect. 9.3 we discuss that a specific boundary-driven pattern may influence the different physical mechanisms of the evolution of the magnetic helicities.

8. Distinguishing between simulations in terms of energy dynamics

As shown in Sect. 5, the evolutions of H_j and E_j are very similar. Likewise, H_pj and E_pj evolve in the same way during the pre-eruption phase. The main difference appears after the eruption, where H_pj still decreases while E_pj remains constant. However, the similar overall behaviors of H_pj and H_j hide very different physical mechanisms, depending on the simulation, as shown in the previous section. We identified three types of evolution for the dynamics of the helicities. In this section we focus on the following questions: how do E_pj and E_j evolve? Does the dynamics of the energy fluxes also distinguish between the different simulations, as the helicity dynamics do?

For this purpose, we present in Fig. 6 all the flux that appear in the decomposition of dE_p/dt (see Eq. (42)) for the convergence (Fig. 6, left panel) and the dispersion central simulations (Fig. 6, middle panel). In both simulations, the nonideal term is almost null because of a very low non-solenoidality of the potential field, for instance, ∇ ⋅ B_p ≃ 0. The evolution of the potential energy therefore depends only on F_ϕ, Bz, which results from the evolution of the magnetic field at the boundaries. During the pre-eruption phase, the magnitude of F_ϕ, Bz as the relative change of the magnetic field at the boundary becomes weaker. Then, during the eruption onset phase, F_ϕ, Bz decreases strongly before becoming null during the eruption phase. The magnetic field is indeed kept fixed at the bottom boundary during that period.

Fig. 6. Left and middle panel: time evolution of the potential energy variation term (dashed black line; dEp/dt; Eq. (42)) and the different terms constituting the instantaneous time-variation of Ep (Eq. (42)): Fϕ, Bz (blue line; Eq. (43)), and dEp/dt|ns (orange line; Eq. (44)). Left and middle panels: evolution for the convergence and dispersion central simulation, respectively. Right panel: time evolution of Fϕ, Bz for the four simulations dispersion central (red line), dispersion peripheral (green line), stretching (yellow line), and convergence (blue line). The yellow band corresponds to the onset phase of the eruption.

From comparing F_ϕ, Bz in the different simulations, we note that F_ϕ, Bz presents the same evolution for all simulations except for the dispersion central (see Fig. 6, right panel). This indicates that except for the dispersion central simulation, the differences in boundary-driven motions do not affect the injection of the potential energy (cf. the discussion in Sect. 9.3). However, this run has the same functional form as the other, is only more efficient, and therefore quicker, in achieving the eruption.

Regarding E_p, only the dispersion central run presents a different behavior. The same conclusion was obtained for the fluxes of H_j, for instance, dH_j/dt|_Own (see Fig. 5, middle top panel), but not for dH_pj/dt|_Own (cf. Fig. 5l). First, F_ϕ, Bz is negative for the entire simulation, while the sign dH_pj/dt|_Own depends on the simulation. This confirms that there is no direct link between the dynamics of E_p and the injection of H_pj.

In Fig. 7 we observe the different terms of dE_j/dt (Eq. (48)) for the four simulations. Unlike dH_j/dt and dH_pj/dt, the trends and dominant terms of dE_j/dt are similar in the simulations. Only the magnitude of each term may differ. Before the eruption, dE_j/dt is dominated by the emergence term F_Bn, Ej despite a significant magnitude of the dissipation term, dE_j/dt|_Diss. Then, during the eruption, F_Bn, Ej becomes null as a result of the interruption of the boundary-driving flows. During the eruption phase, dE_j/dt is negative and dominated by dE_j/dt|_Diss. The free energy mainly decreases because it is dissipated and not ejected out of the volume.

Fig. 7. Time evolution of the free-energy variation rate (dashed black line; dEj/dt; Eq. (48)) and the different terms constituting the instantaneous time-variation of Ej (Eq. (48)): dEj/dt|Diss (blue line; Eq. (49)), dEj/dt|Var (orange line; Eq. (50)), FBn, Ej (green line; Eq. (51)), FVn, Ej (red line; Eq. (52)), and dEj/dt|ns (purple line; Eq. (54)). Each panel corresponds to a different simulation: convergence (top left), stretching (top right), dispersion peripheral (bottom left), and dispersion central (bottom right). The yellow band corresponds to the onset phase of the eruption.

The dissipation term, dE_j/dt|_Diss does not vary much between the simulations (see Fig. 8, bottom right panel). Similarly, the differences of F_Vn, Ej and dE_j/dt|_Var are small between the runs during the pre-eruption phase. Only F_Bn, Ej (see Fig. 7, top left panel) presents significant differences in the simulations that affect the evolution of the free energy, E_j. The dispersion central simulation presents a distinctive trend. The magnitude of F_Bn, Ej starts very high and then decreases to values similar to the other runs during the onset phase of the eruption.

Fig. 8. Time evolution of the different gauge-invariant terms of dEj/dt: FBn, Ej (top right panel, Eq. (51)), FVn, Ej (top left panel, Eq. (52)), dEj/dt|Var (bottom left panel, Eq. (50)), and dEj/dt|Diss (bottom right panel, Eq. (49)). The different colors present one simulation: dispersion central (red line), dispersion peripheral (green line), stretching (yellow line), and convergence (blue line). The yellow band corresponds to the onset phase of the eruption.

Unlike the evolution of the helicities, H_j and H_pj, only the dispersion central simulation stands out from the other runs. This simulation is characterized by a higher decrease of the potential field (see Fig. 6, right panel) and by a higher initial injection of E_j caused by F_Bn, Ej (see Fig. 8, top left panel). Before the eruption, another difference with the helicities is that the variations of the trend of E_j and E_p are purely related to the boundary fluxes. Finally, one key outcome of our study is that the dynamics of the energy fluxes do not distinguish between the simulations, unlike the helicity fluxes.

This shows that even if volume helicities and energies follow similar trends (cf. Fig. 2), the physical mechanisms that drive their dynamics are very different. First, the evolution of free energy, E_j, and potential energy, E_p, are independent, while the current-carrying helicity, H_j, evolves in a correlated way with the volume-threading helicity, H_pj. Additionally, different boundary forcing only affects the magnitude of the energy fluxes. The dynamics of the helicity is more complex and varies drastically from one simulation to another. One group of simulations (dispersion central and peripheral) is dominated by the flux through the surfaces, while a second group (convergence and stretching) is controlled by volume exchange within the domain. We conclude that energy, helicity, and their decompositions have distinct properties whose analysis should be complementary for the study of the eruptivity of active regions.

9. Discussion

9.1. Summary

In Sect. 2 we introduced the formulation of the magnetic relative helicity, H_v, as well as the formulation of its decomposition into the current-carrying helicity, H_j, and the volume-threading helicity, H_pj. We also recalled the analytical equations of their time-variations obtained in Linan et al. (2018). Similarly, we introduced the decomposition of magnetic energy, E_v, into the sum of the potential energy, E_p, and the free energy, E_j. Then, we obtained the time-variation of E_v (see Eq. (37)), E_p (see Eq. (42)), and E_j (see Eq. (48)) by analytically deriving their time derivative (see Sect. 3). These formulae are valid for any gauge choices and in the presence of finite level of non-solenoidality for the magnetic field.

Our numerical study of time-variations of energies and helicities is based on a series of four eruptive numerical MHD simulations of solar active regions (see Sect. 4) that have been investigated in Zuccarello et al. (2015). The evolution of each simulation is characterized by different boundary forcing (line-tied) until the eruption (see Sect. 4). After the same shearing phase, four driving photospheric flows were considered: convergence, stretching, and peripheral and central dispersion flows. In this study we were particularly interested in the fluxes of energies and helicities during the flux rope formation, during the eruption onset phase, when the torus instability occurs, and during a short time interval after the eruption onset, called the eruption phase.

Initially, the magnetic energy, E_v, decreases as a result of the decrease in the potential magnetic field, E_p, despite the increase in free energy, E_j. At the same time, the decrease in volume-threading helicity, H_pj, compensates for the injection of the current-carrying helicity, H_j, which leads to a quasi-constant evolution of the relative helicity H_v (see Sect. 5). The relative helicity was mostly injected during the earlier shearing phase.

The fluxes of free and potential energies, E_j and E_p (see Sect. 8) showed that the effectiveness of the buildup of free energy within the domain is purely related to the magnitude of one surface term. Similarly, the evolution of potential energy is fully linked with its boundary fluxes.

We also used these simulations to investigate the importance of the exchange between H_j and H_pj, which was previously highlighted by Linan et al. (2018). The exchange of helicity between H_j and H_pj is controlled by a gauge-invariant term, dH_j/dt|_Transf (cf. Eq. (10)). As in Linan et al. (2018), we observed that this term plays a key role in the dynamics of both H_j and H_pj, in particular during the buildup phases where the transfer terms, dH_pj/dt|_Transf, dominate the evolution of dH_pj/dt for the runs convergence and stretching. In the dispersion simulations, the evolutions of dH_j/dt and dH_pj/dt are dominated by their boundary fluxes, dH_pj/dt|_Own and dH_pj/dt|_Own (cf. Eqs. (18) and (29)), even though the magnitude of the transfer term remains significant. This means that neither H_j nor H_pj evolve as a result of boundary fluxes alone. This conclusion is consistent with the results of Linan et al. (2018) that were obtained for different numerical experiments. Specifically, H_j and H_pj cannot be estimated in observed active solar active regions by time integration of its flux through the solar photosphere, but rather with a volume-integration method (Valori et al. 2016). This approach requires a 3D reconstruction of the coronal magnetic field from the 2D photospheric measurement with coronal field extrapolation techniques (Wiegelmann & Sakurai 2012; Wiegelmann et al. 2014). A more detailed discussion of the effect of the transfer term on the estimation of H_j and H_pj can be found in the conclusion of Linan et al. (2018).

The key outcome of the study is the observation that the dynamics of the transfer and fluxes of H_j and H_pj depend on the simulation and thus on the imposed driving motions, even though the variations in H_j and H_pj seemed relatively independent of the simulation setup (see Sect. 6). Despite the four boundary forcings, the different simulations remain very similar in terms of the magnetic field topology. Nonetheless, the dominant terms of dH_j/dt and dH_j/dt are not the same from one simulation to another.

We highlighted three distinct types of dynamics of the evolution of the helicities in the simulations. In the convergence and stretching simulations, H_pj does not evolve as a result of boundary fluxes, but because of its conversion into H_j. The opposite is observed during the two dispersion simulations, for which the evolution of H_pj is mainly related to its fluxes through the boundary, with a weak transfer term. The dispersion central simulation stands out from the others because its boundary flux of the current-carrying helicity, H_j, is significantly higher.

Thus we were able to process several photospheric forcings to approach the diversity of active regions that are observed at the solar photosphere. We have come to the conclusion that the evolution of helicity during the formation of a flux rope is a complex process whose origin can be related to fluxes through the surface as well as to volume contributions.

9.2. Buildup of the helicity ratio

Zuccarello et al. (2018) have shown that the trigger of the eruptions is related to a threshold in the helicity ratio H_j/H_v: this ratio reaches the same value, |H_j/H_v|_thresh, for all simulations at the onset of the eruptions. In our runs, this threshold is 0.29 ± 0.01. However, as discussed in Sect. 7 of Zuccarello et al. (2018), the exact value of this threshold needs to be taken with care because relative helicity is not simply an additive quantity. We investigated how this helicity ratio is built up and eventually reached by studying the specific dynamics of H_j and H_pj (see Sects. 6 and 7). Despite the different boundary forcings, the simulations are very similar, so that it might have been thought that the fluxes of H_j and H_pj would also be similar. However, the key outcome of our study is that the terms that dominate the evolution of dH_j/dt or dH_pj/dt sensitively depend on the simulation even if the overall trends are the same (H_j increases and H_pj decreases, see Sect. 5).

Three very distinct ways to reach the helicity eruptivity threshold were found. We observed that the eruption was triggered at a specific value of H_j/H_v independently of the dynamics of H_j and H_pj to reach this threshold. Our work suggests that active regions could reach an eruptive state, either through strong increases of helicity fluxes or through magnetic configurations that induce strong helicity transfer.

The different ways to reach the helicity eruptivity threshold are not all equally effective. The eruption occurs more or less quickly after the end of the shearing phase. The dispersion central run is the most rapid simulation. Then we find the stretching and convergence runs, and last the dispersion peripheral case. The dispersion central simulation stands out from the others because it presents higher helicity fluxes (due to dH_j/dt|_Own and dH_pj/dt|_Own) and energy fluxes (due to F_ϕ, Bz and F_Bn, Ej) than the other cases.

Finally, using a set of resistive MHD simulations, we provided new knowledge of the energy and helicity properties. In particular, our analytical and numerical work emphasizes recent studies (Pariat et al. 2017; Zuccarello et al. 2018; Moraitis et al. 2019; Thalmann et al. 2019) that demonstrated how promising the helicity ratio is as a marker of eruptivity. Further studies are still needed, whether to analytically establish the link between the helicity ratio and torus instability or to properly estimate it from data that are measured in the solar atmosphere.

From direct observational data, the evolution of the ratio H_j/H_v was also analyzed in three active regions with different eruptive profiles. Moraitis et al. (2019) investigated the most active region of cycle 24 (AR 12673), while Thalmann et al. (2019) focused on an eruptive and a confined flare (AR 11158 and AR 12192). These recent observational analyses seem to qualitatively confirm that the H_j/H_v ratio is tightly related to the eruptivity of solar active regions.

9.3. Effect of the different flows on the helicity and energy injections

A key result of our study is that the specific driving flows that are applied at the bottom boundary are of significant importance on the dynamics of magnetic helicities and energies. They influence the magnitude and the sign of the own terms for the helicity as well as those of the main fluxes of dE_j/dt and dE_p/dt. Moreover, as mentioned in the previous section, even if the way to reach the helicity eruptivity threshold matters less than reaching the threshold, the spatial velocity and magnetic distributions at the boundary affect the time that is required to reach the threshold.

To reach an understanding of the effect of the line-tied forcing, we here briefly discuss the distribution of different quantities at the bottom boundary. A more complete study is beyond the scope of the present paper but will be performed, however.

Our goal is to present quantities that might eventually be obtained from observed photospheric magnetograms. Figure 9 shows four quantities that are related to the energy and helicity fluxes in a (x − y) view at z = 0.006 for all simulations at the same modified time to the eruption (t₁ − t_A = −58): (v ⋅ A_j)B_z, the integrand of F_Vn, Aj, for the injection of H_j (see the first row in Fig. 9); (v ⋅ A_j)B_z, the integrand of F_Vn, Ap, for the injection of H_pj (see the second row in Fig. 9); (∂ϕ/∂t)B_p, the integrand of F_ϕ, Bp (see the third row in Fig. 9); and (v ⋅ B_j)B_z, the integrand of F_Bn, Ej (see the fourth row in Fig. 9). The modified time was taken arbitrarily in the pre-eruption phase.

Fig. 9. From top to bottom: dimensionless magnitude of (v ⋅ Aj)Bz, (v ⋅ Ap)Bz, (∂ϕ/∂t)(Bp), and (v ⋅ Bj)Bz viewed in the (x, y) plane at z = 0.006, at the relative time of t − t1 = −58. Isocontours of |Bz| (dashed line for negative values, solid line for positive values) correspond to values of |Bz| = − 4.5, −2.0, 0, 2.0, and 4.5. Each column in the panels presents one simulation, from left to right: convergence, stretching, dispersion peripheral, and dispersion central.

We recall that the injections of H_j or H_pj related to their own terms cannot be simply reduced to F_Vn, Aj and F_Vn, Ap. Other terms such as dH_pj/dt|_Transf must be considered. The second issue is that neither F_Vn, Aj nor F_Vn, Ap are gauge-invariant quantities. They are still good indicators of the spatial distribution of the helicity injection, however, because the same gauge was used for all four runs. Moreover, the emergence term of the relative helicity is one quantity that has traditionally been investigated at the solar photosphere (Liu & Schuck 2012; Bi et al. 2018). The distributions of F_Vn, Aj and F_Bn, Ej are very similar (see the first and last row in Fig. 9). This reveals that the injection of free energy appears to be directly connected to the injection of current-carrying helicity.

The area where F_Vn, Aj and F_Bn, Ej are intense is along the PIL, except for the dispersion central simulation. In the convergence and stretching runs, the shear is favored because the angle between the velocity and the magnetic field at the PIL is high. For the dispersion peripheral case, the velocity is perpendicular to the PIL. Thus the free energy increases more slowly than in the other cases, and more time is needed to build up the flux rope. We also note a low contribution of F_Vn, Aj related to the dispersion of the magnetic field at the periphery of the polarities in the stretching and dispersion peripheral runs.

As highlighted in Sect. 7, the dynamics in the dispersion central simulation stands out from the others. Unlike the other cases, the intense regions of F_Vn, Aj and F_Bn, Ej are located close to the center of the polarities. The difference between this simulation and the other three appears clearly for F_ϕ, Bp (see the third row in Fig. 9). In Sect. 8 we observed that F_ϕ, Bp was distinctly higher for the dispersion central than in the other cases. The magnitude of F_ϕ, Bp is higher for the dispersion central simulation because the center of the polarity, where the magnetic field is the most intense, is displaced. The positive polarity, which possess a higher magnetic flux than the negative polarity, contributes more to F_ϕ, Bp. The three other cases are very similar because the regions subjected to the flow do not influence F_ϕ, Bp.

The distribution of the F_Vn, Ap term (see the second row in Fig. 9) is markedly different from the other quantities presented in Fig. 9. Two different behaviors appear. For the convergence and stretching runs, like F_Vn, Aj, the shear along the PIL is the main contribution of (v ⋅ A_j)B_z. For the dispersion runs, F_Vn, Ap is distributed in a symmetric quadrupole at the location where the flow is applied. This leads to an almost null total flux, F_Vn, Ap.

Finally, the best way to reach the helicity threshold for these parametric simulations is to facilitate the dispersion of the intense magnetic field. For instance, in the dispersion central case, this corresponds to distributing the flow from the center of each polarity. For the other cases with lower dispersion, the efficient ways to build up the energies and helicities are mostly linked with flux-cancellation converging motions that contain a shearing-flow component parallel to the PIL. The convergence and stretching runs have more intense contributions there.

The observations made above offer only a limited insight on the effect of the different flows on the evolution of energy and helicity. We saw that even if the four simulations look similar overall, some differences related to the applied motions clearly appear in the dynamic of the helicity and energy. New investigations are required to provide a better understanding of the energy and helicity dynamics. Our study especially provides new information for interpreting the injection of helicity and energy in observed active regions.

With HMI vector magnetic field data, the evolution of energy and helicity flux has previously been studied in different active regions (Liu & Schuck 2012). This is commonly based on the decomposition of the flux into two components: a shear component provided by photospheric tangential flow, and a vertical component linked with normal flows due to emergence (Bi et al. 2018). At the same time, new methods for properly measuring helicities in the solar corona are still being developed (Dalmasse et al. 2013, 2014, 2018; Valori et al. 2016; Guo et al. 2017; Moraitis et al. 2018; Gosain & Brandenburg 2019). For instance, new analytic expressions for the helicity transport allow us to estimate the injection of relative magnetic helicity into the solar atmosphere over an entire solar cycle (Hawkes & Yeates 2019; Pipin et al. 2019). However, the quality of the helicity and energy estimations from observational data greatly depends on the accuracy of the magnetic field measurements. In particular, it is still difficult to measure the tangential component of the magnetic field. The instrumentation on board the new Solar Orbiter mission, for example, the PHI magnetogram, will help us to make significant progress on the magnetic field measurement and consequently on the helicity and energy computations. Future results that will benefit from our study will provide new insight for a better understanding of the solar coronal activity.

Acknowledgments

The authors thank the anonymous reviewer for careful reading of the manuscript and helpful comments. L.L., E.P., K.M. acknowledge support of the French Agence Nationale pour la Recherche through the HELISOL project ANR-15-CE31-0001. L.L., E.P. acknowledge the support of the french Programme National Soleil-Terre. GV acknowledges the support from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 824135 and of the STFC grant number ST/T000317/1. The numerical simulation used in this work was performed on the HPC center MesoPSL financed by the project Equip@Meso (reference ANR-10-EQPX-29-01) and the Region Ile-de-France. This article profited from discussions during the meetings of the ISSI International Team Magnetic Helicity in Astrophysical Plasmas.

References

All Figures

	Fig. 1. Applied boundary-driving motions for the four different numerical experiments. White represents the positive polarity (Bz(z = 0.006) > 0) and black the negative polarity (Bz(z = 0.006) < 0). Orange and cyan arrows indicate the distribution of the velocity flows we applied to the negative and positive polarity, respectively.
In the text

Fig. 2. Evolution of the different magnetic helicities (top panels), from left to right: relative magnetic helicity (Hv, Eq. (4)), volume-threading helicity (Hpj, Eq. (7)), and current-carrying helicity (Hj, Eq. (6)). Time evolution of the different magnetic energies (bottom panel), from left to right: total magnetic energy (Ev, Eq. (32)), potential energy (Ep, Eq. (31)), and free energy (Ej, Eq. (31)). The different simulations are dispersion central (red line), dispersion peripheral (green line), stretching (yellow line) and convergence (blue line). The yellow vertical band corresponds to the onset phase of the eruption.

In the text

Fig. 3. Left panel: time evolution of the instantaneous time-variation, dEv/dt (dashed black curves, Eq. (37)), and of the sum, dEj/dt + dEp/dt (continuous blue curves, Eqs. (8) and (19)) for the dispersion peripheral run. Right panel: time evolution of the instantaneous time-variation, dHv/dt (dashed black curves, Eq. (23) in Pariat et al. 2015, and of the sum, dHj/dt + dHpj/dt (continuous blue curves, Eqs. (8) and (19)) for the dispersion peripheral run. The yellow bands correspond to the eruption onset phase.

In the text

Fig. 4. Time evolution of the helicity variation rates, dHj/dt and dHpj/dt (dashed black curves; Eqs. (8) and (19)), of the helicity transfer term, dHj/dt|Transf and dHpj/dt|Transf (solid red curves; Eqs. (10) and (21)), of the own terms, dHj/dt|Own and dHpj/dt|Own (solid blue curves; Eqs. (18) and (29)), and of the dissipation terms, dHj/dt|Diss and dHpj/dt|Diss (solid green curves; Eqs. (9) and (20)) for the convergence simulation (top panels) and for the dispersion peripheral simulation (bottom panels). Left and right panels: evolution of the current carrying helicity, Hj, and volume-threading helicity Hpj, respectively.

In the text

Fig. 5. Time evolution of the different gauge invariant terms of dHj/dt (top panels), from left to right: dissipation term (dHj/dt|Diss, Eq. (9)), own term (dHj/dt|Own, Eq. (18)), and helicity transfer term (dHj/dt|Transf, Eq. (10)). Time evolution of the different gauge invariant terms of dHpj/dt (bottom panels), from left to right: dissipation term (dHpj/dt|Diss, Eq. (20)), and own term (dHpj/dt|Own, Eq. (29)). Each curve color corresponds to a particular simulation: dispersion central (red line), dispersion peripheral (green line), stretching (yellow line), and convergence (blue line). The yellow band corresponds to the onset phase of the eruption.

In the text

Fig. 6. Left and middle panel: time evolution of the potential energy variation term (dashed black line; dEp/dt; Eq. (42)) and the different terms constituting the instantaneous time-variation of Ep (Eq. (42)): Fϕ, Bz (blue line; Eq. (43)), and dEp/dt|ns (orange line; Eq. (44)). Left and middle panels: evolution for the convergence and dispersion central simulation, respectively. Right panel: time evolution of Fϕ, Bz for the four simulations dispersion central (red line), dispersion peripheral (green line), stretching (yellow line), and convergence (blue line). The yellow band corresponds to the onset phase of the eruption.

In the text

Fig. 7. Time evolution of the free-energy variation rate (dashed black line; dEj/dt; Eq. (48)) and the different terms constituting the instantaneous time-variation of Ej (Eq. (48)): dEj/dt|Diss (blue line; Eq. (49)), dEj/dt|Var (orange line; Eq. (50)), FBn, Ej (green line; Eq. (51)), FVn, Ej (red line; Eq. (52)), and dEj/dt|ns (purple line; Eq. (54)). Each panel corresponds to a different simulation: convergence (top left), stretching (top right), dispersion peripheral (bottom left), and dispersion central (bottom right). The yellow band corresponds to the onset phase of the eruption.

In the text

Fig. 8. Time evolution of the different gauge-invariant terms of dEj/dt: FBn, Ej (top right panel, Eq. (51)), FVn, Ej (top left panel, Eq. (52)), dEj/dt|Var (bottom left panel, Eq. (50)), and dEj/dt|Diss (bottom right panel, Eq. (49)). The different colors present one simulation: dispersion central (red line), dispersion peripheral (green line), stretching (yellow line), and convergence (blue line). The yellow band corresponds to the onset phase of the eruption.

In the text

Fig. 9. From top to bottom: dimensionless magnitude of (v ⋅ Aj)Bz, (v ⋅ Ap)Bz, (∂ϕ/∂t)(Bp), and (v ⋅ Bj)Bz viewed in the (x, y) plane at z = 0.006, at the relative time of t − t1 = −58. Isocontours of |Bz| (dashed line for negative values, solid line for positive values) correspond to values of |Bz| = − 4.5, −2.0, 0, 2.0, and 4.5. Each column in the panels presents one simulation, from left to right: convergence, stretching, dispersion peripheral, and dispersion central.

In the text