© ESO 2018

1 Introduction

Magnetic helicity is a key geometrical parameter used to describe the structure of solar coronal magnetic fields (e.g. Berger 1999). Magnetic helicity in a volume V can be determined as $$H_{\text{M}}={\displaystyle {\int }_{V}A}\cdot B\text{d}V,$$(1)

where A is a vector potential and B is magnetic field in this volume. Magnetic helicity H_M is conserved in an ideal magnetoplasma (Woltjer 1958). As long as the overall magnetic Reynolds number is large, it is still approximately conserved, even in the case of relatively slow magnetic reconnection taking place (Taylor 1974; Berger & Field 1984). Since the vector potential A is not uniquely defined, H_M is not gauge-invariant. The magnetic helicity in a volume has a well-determined value only when the magnetic field at the boundary is exclusively tangential, i.e. if $B\cdot \widehat{n}{|}_S=0$. On the other hand, Berger & Field (1984) have shown that in the case of boundaries open to magnetic flux penetration, relative magnetic helicity (H_R), given by the Finn-Antonsen (1985) formula $$H_{\text{R}}={\displaystyle \int \left(A+A_{\text{p}}\right)\cdot \left(B-P\right)}\text{d}V,$$(2)

is gauge-invariant if the magnetic field is chosen as a reference field such that $P\cdot \widehat{\text{n}}{|}_S=B\cdot \widehat{\text{n}}{|}_S$. It is customary to choose the potential field (∇×P = 0) as the reference field. The concept of magnetic helicity has been successfully applied to characterize the solar coronal process and the solar dynamo to interpret solar observations and corona simulations. Chae (2001) calculate the relative helicity by applying the Fourier transform (FFT) approach to obtain vector potential and the local correlation tracking (LCT) to obtain the flow velocity to the Michelson Doppler Imager (MDI) data. Since then the helicity calculation method has been developed further (Démoulin & Berger 2003; Pariat et al. 2005; Schuck 2008) by improving the helicity density map and velocity tracking techniques. The analysis of magnetic helicity based on observations and simulations has also been developed in, for example the analysis of magnetic helicity injection in the course of flux emergence (Yang et al. 2009a; Liu et al. 2014),the correlation between helicity change and solar eruption (Zhang et al. 2008; Park et al. 2008; Yang et al. 2015), magnetic energy-helicity relation analysis in solar eruptions (Tziotziou et al. 2012), magnetic helicity distribution in different scales in the solar dynamo (Seehafer et al. 2003; Pipin & Pevtsov 2014), magnetic helicity estimation in the interplanetary magnetic cloud (Janvier et al. 2015), and testing the Taylor conjecture in quasi-ideal magnetohydrodynamics (MHD) simulations (Pariat et al. 2015). For a recent review of modelling and observations of magnetic helicity, see e.g. Démoulin & Pariat (2009).

Although the magnetic helicity is conserved in the fast reconnection process in the close volume, the redistribution of magnetic helicity, i.e. helicity transport could still happen as it is strongly coupled with the magnetic energy release process. For example, the magnetic helicity exchange process has been found between neighbouring emerging active regions (Yang et al. 2009a). Thus, the correlation study between magnetic helicity distribution and evolution and solar eruption become important. Zhang et al. (2006) noted that the accumulation of magnetic helicity in the corona plays a significant role in storing magnetic energy. They propose that there is an upper bound on the total magnetic helicity that a force-free field can contain. Nindos & Andrews (2005) found that magnetic helicity of coronal mass ejection (CME) productive active regions (ARs) is higher than other ARs. The survey of LaBonte et al. (2007) of helicity accumulation in 393 ARs revealed that a necessary condition for the occurrence of an X-flare is that the peak helicity flux has a magnitude larger than 3 × 10³⁶Mx²∕s. In the simulations the debates also exist for the possible upper-bound helicity before the solar eruption happens (Amari et al. 2003; Jacobs et al. 2006). In these papers, the corresponding normal relative magnetic helicity (H_∕Φ²) when the eruption happens reaches approximately −0.16 for the case of Kink instability, and approximately −0.18 for the case of torus instability in the simulations of Fan & Gibson (2007). However, from the observations the range of normalized helicity for the active region ranges from 0.02 to 0.08 (LaBonte et al. 2007; Yang et al. 2009a; Tian & Alexander 2008). This value is one order of magnitude smaller than the above simulations. Yang et al. (2013) investigate the value of normalized helicity, which reached 0.0298 just prior to a drastic energy release by magnetic reconnection, by using the magnetic field above active region NOAA 8210 obtained by a photospheric-data-driven MHD model (Santos et al. 2011). Yang et al. (2015) studied an emerging and quickly decaying active region (NOAA 9729) in detail as it passed across the solar disk. There was only one CME associated with that active region. This provided a good opportunity to determine the consequences of a single CME after the injection of magnetic helicity to the solar corona. The absolute value of normalized magnetic helicity was 0.036 just before solar eruption happened (Yang et al. 2015), which is close to the simulation of NOAA 8210 and the theoretical prediction value of Zhang & Flyer (2008) only in a multipolar force-free magnetic field structure. However, whether the relative magnetic helicity has a upper-bound and whether it plays an important role in solar eruptions are still open questions.

Despite its important role in the dynamical evolution of solar plasmas, so far only a few attempts have been made to estimate the helicity of coronal magnetic fields based on observations and numerical simulations (see e.g. Thalmann et al. 2011; Rudenko & Myshyakov 2011; Valori et al. 2012).Yang et al. (2013); developed a method for calculating the relative magnetic helicity in a finite 3D volume which was applied to a simulated flaring AR 8210 (Santos et al. 2011). However, this method required the magnetic flux to be balanced on each of the side boundaries in the considered volume. In this paper, the corresponding scheme is presented which does not require this restriction on the flux balance. The methodhas already been applied in Valori et al. (2016), and it is the most accurate of the finite volume methods employing the Coulomb gauge.

We applied our new method to the magnetic field above active region NOAA 11429 as it was obtained by a photospheric-data-driven MHD model GOEMHD3 (Skála et al. 2015). In Sect. 2, we describe the restriction of vector potential in the previous paper. In Sect. 3, we present the new scheme to calculate the vector potentials on the six boundaries. In Sect. 4, we use a non-linear force-free magnetic field model (Low & Lou 1990) and a MHD simulation model to check our scheme. The summary and discussions are given in Sect. 5.

2 Restricted method to obtain A_p and A at the boundaries

The computation of H_R in Eq. (2) requires the knowledge of A, A_p, and P from the known B in the volume V. We adopt the Coulomb gauge for the vector potentials (see Sect. 2.2 in Yang et al. 2013 for details). Let us define a finite three-dimensional (3D) rectangular volume in Cartesian coordinates. The magnetic field B (x, y, z) is given in this volume. The volume is restricted to x = [0, l_x], y = [0, l_y], and z = [0, l_z].

In order to solve for the vector potentials, the boundary conditions must be specified. First, we obtain the values of A _p and A on all six boundaries (x = 0, l_x;y = 0, l_y;z = 0, l_z). Taking the bottom boundary (z = 0) as an example, we define a new scalar function φ(x, y) that determines the vector potential A_p of potential magnetic field P corresponding to B on this boundary as follows: $$A_{\text{p}x}=-\frac{\partial \phi }{\partial y},A_{\text{p}y}=\frac{\partial \phi }{\partial x},{A_{\text{p}z}|}_{z=0}=0.$$(3)

According to the definition of the vector potential, the scalar function φ(x, y) satisfies the Poisson equation: $$\Delta \phi (x,y)=B_z(x,y,z=0).$$(4)

In our previous work (Yang et al. 2013) we set the values of ∂φ∕∂n at the four edges of the plane z = 0 to zero in Eq. (4). As a consequence, the corresponding magnetic flux at the boundary also vanish according to Ampère’s law. The values of A_p on the otherfive boundaries can be obtained in the same way. For the vector potential A at all boundaries, the same values are taken as for A_p. When the magnetic fluxes through the six boundaries are finite, we should calculate the values of vector potentials at the 12 edges of the 3D volume to obtain a Neumann boundary condition for the Poisson Eq. (4) at each side boundary. In the next section, we introduce a scheme to calculate the vector potentials at the 12 edges.

3 General method to obtain A_p and A at the boundaries

In order to determine A_p, we define the magnetic flux Φ_i (i = 1, …, 6) respectively at each side boundary (z = 0; z = l_z; x = 0; x = l_x; y = 0; y = l_y). The integrals of ∫ A_p ⋅dl at the 12 edges are defined as a_i (i = 1, …, 12). The 12 integrals and the corresponding directions are depicted in Fig. 1. Applying Ampère’s law to the six boundaries independently, we obtain a linear system of equations for obtaining the integral value a_i $$TX=F,$$(5)

where $F={({\Phi }_1,{\Phi }_2,{\Phi }_3,{\Phi }_4,{\Phi }_5,{\Phi }_6)}^T$, $X={(a_1,a_2,a_3,\mathrm{...},a_{12})}^T$, and $\widehat{T}$ is a 6 × 12 matrix given by $$\widehat{T}=\left[\begin{array}{cccccccccccc}1& 1& 1& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 1& 1& 1& 0& 0& 0& 0\\ 0& 0& 0& -1& 0& 0& 0& -1& 1& 0& 0& 1\\ 0& -1& 0& 0& 0& -1& 0& 0& 0& 1& 1& 0\\ -1& 0& 0& 0& -1& 0& 0& 0& -1& -1& 0& 0\\ 0& 0& -1& 0& 0& 0& -1& 0& 0& 0& -1& -1\end{array}\right].$$(6)

The six rows in this matrix are not linearly independent because the magnetic field is divergence-free and the sum of Φ_i through all six boundaries vanishes. Therefore, there are no unique solutions for the 12 integrals a_i . We use the freedom in the gauge to construct 12 independent equations to obtain unique solutions for a_i . We remove the last row of Eq. (6) and add another seven rows to define a new matrix $\widehat{T^{\prime }}$ to make sure the determinant of $\widehat{T^{\prime }}$ is not zero without limiting the validity of the solution for a_i. We choose the following seven rows to construct the new matrix $\widehat{T^{\prime }}$: $$\widehat{T^{\prime }}=\left[\begin{array}{cccccccccccc}1& 1& 1& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 1& 1& 1& 0& 0& 0& 0\\ 0& 0& 0& -1& 0& 0& 0& -1& 1& 0& 0& 1\\ 0& -1& 0& 0& 0& -1& 0& 0& 0& 1& 1& 0\\ -1& 0& 0& 0& -1& 0& 0& 0& -1& -1& 0& 0\\ 1& 0& -1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& -1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& -1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& -1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& -1& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1& 0& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& -1& 0& 0& 0\end{array}\right].$$(7)

This choice for the lower seven rows of $\widehat{T^{\prime }}$ is not unique, but any influence on the gauge will be removed later. We can calculate that the determinant of $\widehat{T^{\prime }}$ does not vanish. According to Cramer’s rule, a unique solution of the new linear equation $$\widehat{T^{\prime }}X=F^{\prime },$$(8)

can exists if $X={({\text{a}}_1,{\text{a}}_2,{\text{a}}_3,\mathrm{...},{\text{a}}_{12})}^{\text{T}}$ and $F^{\prime }={({\Phi }_1,{\Phi }_2,{\Phi }_3,{\Phi }_4,{\Phi }_5,0,0,0,0,0,0,0)}^T$. Then we can obtain the integralsof a_i at the 12 edges such that Gauss’s theorem applied to the Coulomb vector potentials in the volume V is satisfied. Since only the integral values a_i are relevant in Eq. (3), the vector potentials at the 12 edges can be obtained representing the components of A _p in the following way: $$\begin{array}{l}{\text{A}}_{\text{px}}\left(a_i\right)=\frac{\pi a_i}{2L_x}\sin (\pi x/L_x),i=1,3,5,7\\ {\text{A}}_{\text{py}}\left(a_i\right)=\frac{\pi a_i}{2L_y}\sin (\pi y/L_y),i=2,4,6,8\\ {\text{A}}_{\text{pz}}\left(a_i\right)=\frac{\pi a_i}{2L_z}\sin (\pi z/L_z),i=9,10,11,12\end{array}$$(9)

It should be noted that such A_p by construction vanishes at the ends of every edge, i.e. at every corner of the box as required by Eq. (3). Now we solve the Poisson Eq. (4) to obtain A_p at the six boundaries. The vector potential A at the six boundaries is equal to A _p . The following procedure corresponds to the Sects. 2.2 and 2.3 of Yang et al. (2013). Using the boundary conditions above, A _p and A in the volume are obtained by solving the Poisson and Laplace equations $$\Delta A_{\text{p}}=0$$(10)

everywhere in the volume if A_p satisfies the Coulomb gauge. The vector potential A of the original magnetic field B satisfies the Poisson equation $$\Delta A=-J,$$(11)

where J = μ₀j denotes the current density. In order to remove numerical errors and imperfections in the gauge, we implement a projection method that removes violations of the solenoidal constraint in the volume. In particular, we introduce a solenoidal modification vector ∇ × M satisfying the following condition: $$\nabla \times \left(\nabla \times M\right)=P-\nabla \times A_{\text{p}}.$$(12)

The components of M satisfy the three Poisson equations: $$\begin{array}{cc}\{\begin{array}{l}\Delta M_z={\left(\nabla \times A_{\text{p}}\right)}_z-P_z\hfill \\ M_z\left(z=0,l_z\right)=0\hfill \\ \frac{\partial M_z}{\partial x}\left(x=0,l_x\right)=0\hfill \\ \frac{\partial M_z}{\partial y}\left(y=0,l_y\right)=0\hfill \end{array},& \{\begin{array}{l}\Delta M_y={\left(\nabla \times A_{\text{p}}\right)}_y-P_y\hfill \\ M_y\left(y=0,l_y\right)=0\hfill \\ \frac{\partial M_y}{\partial x}\left(x=0,l_x\right)=0\hfill \\ \frac{\partial M_y}{\partial z}\left(z=0,l_z\right)=0\hfill \end{array},\\ \{\begin{array}{l}\Delta M_x={\left(\nabla \times A_{\text{p}}\right)}_x-P_x\hfill \\ M_x\left(x=0,l_x\right)=0\hfill \\ \frac{\partial M_x}{\partial y}\left(y=0,l_y\right)=0\hfill \\ \frac{\partial M_x}{\partial z}\left(z=0,l_z\right)=0\hfill \end{array}.& \end{array}$$(13)

Finally, in order to remove the residual errors in the solenoidal property of A_p, we introduce a scalar field ϕ(x, y, z) which satisfies the following Poisson equation: $$\{\begin{array}{l}\Delta \varphi =-\nabla \cdot A_{\text{p}}\hfill \\ {\frac{\partial \varphi }{\partial n}|}_s={-\left(\nabla \times M\right)\cdot \stackrel{⌢}{n}|}_s\hfill \end{array}.$$(14)

After solving Eq. (14), we obtain a new modified vector potential A′_p, which is represented as $${A^{\prime }}_{\text{p}}=A_{\text{p}}+\nabla \times M+\nabla \varphi$$(15)

and in the same way, we can obtain the corrected vector potential A ′ just by replacing the right-hand side term of Eq. (13) as ∇×A −B. Now we can calculate the relative helicity in the volume according to Eq. (1). As already mentioned when introducing $\widehat{T^{\prime }}$, the choice of the lower seven rows of $\widehat{T^{\prime }}$ and the final seven elements of F′ is not unique.

Fig. 1 Magnetic flux Φi (i = 1, …, 6) through the six boundaries and integration paths ai (i = 1, …, 12) along the 12edges of the volume.

4 Testing the scheme

To test the new scheme for obtaining the vector potentials at the boundaries, we used the analytical non-linear force-free fields of Low & Lou (1990). We utilized the model labelled P_1,1 with l = 0.3 and Φ = π∕2 in the notation of their paper. We calculated the magnetic field on a uniform grid of 64 × 64 × 64.

We first calculated the magnetic fluxes Φ₀ through thesix boundaries and substituted it into Eq. (8) to obtain the integrals a_i along the 12 edges of the 3D volume. Then we substituted a_i into Eq. (9) to obtain the boundary value to solve the Poisson Eq. (4) at the six boundaries. After we obtained A _p on the six boundaries, we calculated the magnetic fluxes Φ from the computed vector potential using the relation between the vector potential and the magnetic field: $B\cdot \widehat{\text{n}}=\nabla \times A_{\text{p}}\cdot \widehat{\text{n}}$. Table. 1 represents the quantitative result after applying the above scheme. As can be seen in the table, the magnetic fluxes through the six boundaries calculated via our new scheme is actually better flux-balanced than the theoretical model. Suchsmall errors are partly due to the fact that the total magnetic flux of the analytical model does not completely vanish. The total magnetic flux must be zero in order to resole the linear Eq. (8). Numerical errors in solving the Poisson equation, on the other hand, will also introduce finite total magnetic flux. The calculated relative magnetic helicity using the previous method is −14 441.45, while the value is −13 253.36 using the new method. There is an 8% difference for the two methods.

We then applied our new scheme to recalculate the relative magnetic helicity above active region NOAA AR8210 in Yang et al. (2013). Figure 2 depicts the comparison of relative magnetic helicity evolution by using the previous and the newscheme. The solid (dashed) line represents the new (previous) relative magnetic helicity evolution. The very small difference (<1.2%) is due to the magnetic field structure of the corona above AR8210 obtained by Yang et al. (2013) where the fluxes through the side boundaries were almost balanced.

Finally, we applied our new scheme to the 3D magnetic field data obtained from the simulated evolution of the solar corona above activeregion NOAA AR11429. The new GOEMHD3 code was used to reveal the magnetic field evolution in the solar atmosphere in response to the energy influx from the chromosphere through the transition region. The weak Joule current dissipation and a finite viscosity is taken into account in the almost dissipationless solar corona. The GOEMHD3 code is a massively parallel code for solving second-order accurate MHD equations (Skála et al. 2015). It was successfully tested and applied to study the magnetic coupling between the solar photosphere and corona based on multiwavelength observations. GOEMHD3 discretizes the ideal part of the MHD equations using a fast and efficient leap-frog scheme that is second-order accurate in space and time and whose initial and boundary conditions can easily be modified. For the investigation of diffusive and dissipative processes the corresponding terms are discretized by a DuFort-Frankel scheme. To always fulfil the Courant-Friedrichs-Lewy stability criterion, the time step of the code is adapted dynamically. Non-equidistant grids enhance the spatial resolution near the transition region. GOEMHD3 is parallelized based on a hybrid MPI-OpenMP programming paradigm, adopting a standard 2D domain-decomposition approach. This allows us to investigate the long time evolution of the relative magnetic helicity in the solar corona both in ideal and non-ideal magnetohydrodynamics by the non-ideal magnetohydrodynamical phase triggered by the switching on of resistivity.

Figure 3 shows the evolution of the relative magnetic helicity, using our method to simulate the data from the application of GOEMHD3, to the evolution of AR 11429 on March 07, 2012 (Skála et al. 2015). The vertical dash-dotted line in Fig. 3 indicates the time (t = 1035s) when fast reconnection started. Figure 3a depicts the evolution of relative magnetic helicity in the simulation box (solid line) and the accumulated relative helicity (dashed line) calculated using the helicity change rate Eq. (6) in the Yang et al. (2013). The surface integral in Eq. (6) of Yang et al. (2013) is extended to all six boundaries and the dissipative volume integral in Eq. (6) of Yang et al. (2013) is also included in the computation.

Figure 3b depicts the normalized magnetic helicity evolution (H_∕Φ²). We used two types of fluxes Φ as the normalized parameter. One is unsigned magnetic flux only through the bottom boundary (solid line in Fig. 3b), which is usually used in observations. The other is unsigned magnetic flux through the six side boundaries (dashed line in Fig. 3b). A value of −0.038 (−0.025) was reachedafter fast reconnection started while choosing unsigned magnetic flux through the bottom boundary (side boundaries) as the normalized parameter. The calculated normalized value −0.039 (−0.026) is consistent with the above simulation results if we substitute the two types of fluxes into the formula between the magnetic helicity and themagnetic flux from the observations of newly emerging active regions (Yang et al. 2009b) $\log \frac{\Delta H}{H_0}=a\log \frac{{\Phi }_m}{{\Phi }_0}+b,$(16)

where a = 1.85, b = −0.41, H₀ = 10⁴¹Mx², and Φ₀ = 10²¹Mx., which describes the relation between the accumulated helicity ΔH and magnetic flux Φ_m for active regions. Figure 3c depicts the magnetic helicity flux through the bottom boundary into the solar corona. The total magnetic helicity flux is shown in Fig. 3d. We find that the negative magnetic helicity is injected into the Corona due to the subphotospheric plasma motion. The total injected helicity is − 8.5 × 10⁴²Mx². Figure 3e depicts the dissipation of magnetic helicity flux calculated as − 2 ∫ E ⋅ BdV . The helicity dissipation is enhanced due to fast magnetic reconnection, i.e. non-ideal magnetohydrodynamics. The integrated dissipated magnetic helicity flux is shown in Fig. 3f. We find that the sign of dissipation is positive, i.e. opposite to the helicity injected through the bottom boundary. The total magnetic helicity dissipation in the reconnection is 0.71 × 10⁴²Mx². We note that only 8% of the injected magnetic is dissipated. This reflects the fact that magnetic helicity is approximately conserved during fast magnetic reconnection as pointed out by (Taylor 1974). Figure 3g depicts the magnetic helicity flux which escapes through the upper five boundaries. The total escaped magnetic helicity evolution is shown in Fig. 3h. We find that only − 0.16 × 10⁴²Mx² of the magnetic helicity flux escapes through the upper boundaries. This is 23% of the dissipated helicity flux (0.71 × 10⁴²Mx², see Fig. 3f) and only 2% of the total injected magnetic helicity (− 8.5 × 10⁴²Mx², see Fig. 3d). We note that only the helicity flux through all six boundaries is gauge invariant, and that fluxes though individual boundaries have values that may change by changing gauge (Pariat et al. 2015). In the calculation of Fig. 3c-g at the separate boundary, we use the Coulomb gauge.

Fig. 2 Comparison of the magnetic helicity calculation by applying the previous and new schemes to the magnetic field data of Yang et al. (2013). The solid (dashed) line represents the new (previous) relative magnetic helicity evolution.

Fig. 3 Relative magnetic helicity evolution calculation using the result of the GOEMHD3 simulation. The dash-dotted line represents the time when the fast reconnection started (t = 1035s). (a) Relative magnetic helicity in the simulation box (dashedline) and the accumulated relative helicity (solid line) sum of the helicity injection through the boundary and the dissipated in the volume helicity according to Eq. (6) in Yang et al. (2013). (b) Evolution of the normalized relative helicity, where the magnetic helicity is normalized by unsigned flux at the bottom boundary (solid line) and the six side boundaries (dashed line). (c) Injected magnetic helicity flux through the bottom boundary. (d) Injected magnetic helicity from the bottom boundary. (e) Magnetic helicity dissipation rate in the simulation box. (f) Magnetic helicity dissipation in the simulation box. (g) Magnetic helicity change rate escaped from the other five boundaries except the bottom boundary. (h) Magnetic helicity escaped from the upper five boundaries.

5 Summary and discussion

We propose a new generalized scheme to calculate the vector potential at the boundaries of a closed volume without the restrictions of the method applied in the previous paper (Yang et al. 2013). We verified the new method using an analytical theoretical force-free model magnetic field (Low & Lou 1990) and the simulated data in Yang et al. (2013). We also applied the new methodto simulated coronal magnetic fields using the newly developed GOEMHD3 to investigate the magnetic helicity evolution of solar corona in the course of the evolution of AR11429.

We find that only 8% of the accumulated injected magnetic helicity is dissipated. This is consistent with the Taylor conjecture and also with the recent simulation results of Pariat et al. (2015). Only 2% of the magnetic helicity injected through the bottom boundary escapes to the solar wind. This shows that the magnetic helicity cannot efficiently escape. This may help us to understand the lack of magnetic clouds¹ with considerable magnetic helicity in the interplanetary space (Janvier et al. 2015). This is true even though several halo CMEs originated from the simulated active region AR11492.

Our simulation results further confirm that the absolute normalized magnetic helicity (H_∕Φ²) reaches 0.038 after the reconnection starts. It is interesting to note that recently an isolated and quickly decaying active region (NOAA 9729) was studied in detail as it passed across the solar disk. There was only one CME associated with that active region. This provided a good opportunity to investigate the consequences of a single CME after the injection of magnetic helicity to the solar corona. The absolute value of normalized magnetic helicity was 0.036 just after solar eruption happened (Yang et al. 2015). This is also similar to the obtained simulation results. A normalized helicity is one order of magnitude smaller than the value −0.16 (−0.18) obtained by the MHD simulation results of kink or torus instability of ideal MHD (Fan & Gibson 2007) when eruption happened. This might be due to the anomalous resistivity caused by the strong micro-turbulence and microscopic structures (Büchner & Elkina 2006) used in the GOEMHD3 model. The anomalous current dissipations allows the dissipation rate of magnetic energy to essentially increase following a global MHD instability (Büechner et al., in preparation).

Zhang et al. (2006) and Zhang & Flyer (2008) proposed that there is an upper bound of the total magnetic helicity that a force-free field can contain in a multipolar force-free magnetic field structure before eruption, which is very close to the value obtained above. However, the allowable level of helicity in a force-free field is only a sufficient condition for an eruption. A CME expulsion may still occur even before this helicity limit is reached (Zhang et al. 2006). The upper-bound normalized magnetic helicity value clearly deviates for different magnetic structures in the theoretical force-free field model. For multipolar fields, the helicity upper bound can be 10 times smaller than that of a dipolar field (Zhang & Flyer 2008), but in our observations and simulation result, the normalized helicity reaches the same order when eruption happens. It is essential to investigate the evolution of magnetic helicity and energy using more observations and simulations in future studies.

Above all, after introducing a new scheme removing the former restriction on the magnetic flux through the boundaries, we can calculate the relative magnetic helicity of any magnetic field structure in Cartesian coordinates. In the observations, we could use a force-free extrapolation to obtain the 3D magnetic structure to analyse the evolution of relative magnetic helicity.

Acknowledgements

We would like to thank the referee for carefully reading our manuscript and for giving such constructive comments which substantially helped to improve the paper. This study is supported by grants 11427901, 10921303, 11673033, U1731113,11611530679, and 11573037 of the National Natural Science Foundation of China; by grant No. XDB09040200, XDA15010700 of the Strategic Priority Research Program of Chinese Academy of Sciences; by the Max-Planck Society Interinstitutional Research Initiative Turbulent transport and ion heating, reconnection, and electron acceleration in solar and fusion plasmas of Project No. MIF-IF-A-AERO8047; by Max-Planck-Princeton Center for Plasma Physics PS AERO 8003; by ISSI International Team on Magnetic Helicity estimations in models and observations of the solar magnetic field. The authors would also like to thank the Supercomputing Center of the Chinese Academy of Sciences (SCCAS) and the Max-Planck Computing and Data Facility (MPCDF) for the allocation of computing time.

References

All Tables

All Figures

	Fig. 1 Magnetic flux Φi (i = 1, …, 6) through the six boundaries and integration paths ai (i = 1, …, 12) along the 12edges of the volume.
In the text

	Fig. 2 Comparison of the magnetic helicity calculation by applying the previous and new schemes to the magnetic field data of Yang et al. (2013). The solid (dashed) line represents the new (previous) relative magnetic helicity evolution.
In the text

Fig. 3 Relative magnetic helicity evolution calculation using the result of the GOEMHD3 simulation. The dash-dotted line represents the time when the fast reconnection started (t = 1035s). (a) Relative magnetic helicity in the simulation box (dashedline) and the accumulated relative helicity (solid line) sum of the helicity injection through the boundary and the dissipated in the volume helicity according to Eq. (6) in Yang et al. (2013). (b) Evolution of the normalized relative helicity, where the magnetic helicity is normalized by unsigned flux at the bottom boundary (solid line) and the six side boundaries (dashed line). (c) Injected magnetic helicity flux through the bottom boundary. (d) Injected magnetic helicity from the bottom boundary. (e) Magnetic helicity dissipation rate in the simulation box. (f) Magnetic helicity dissipation in the simulation box. (g) Magnetic helicity change rate escaped from the other five boundaries except the bottom boundary. (h) Magnetic helicity escaped from the upper five boundaries.

In the text