Evolution of Relative Magnetic Helicity: New Boundary Conditions for the Vector Potential
Shangbin Yang Joerg Buechner Jan Skala Hongqi Zhang
astro-ph.SR

 

   

   Key Laboratory of Solar Activity, National
Astronomical Observatories, Chinese Academy of Sciences, 100012
Beijing, China 

             yangshb@bao.ac.cn
         Max-Planck Institute for Solar System Research, 37077 Göttingen, Germany
     University of Chinese Academy of Sciences, 100049 Beijing,China
     Astronomical Institute of Czech Academy of Sciences, 25165 Ondrejov, Czech Republic
     University J. E. Purkinje, 40096 Usti nad Labem, Czech Republic





For a better understanding of the dynamics of the solar corona it is important to analyze the evolution of the helicity of the magnetic field. Since the helicity cannot be directly determined by observations  we recently proposed a method to calculate the relative magnetic helicity in a finite volume for a given magnetic field in <cit.>, which however required the flux to be balanced separately on all the sides of the considered volume.


Development of a scheme to obtain the vector potential in a volume without the above restriction at the boundary. Study of the dissipation and escape of relative magnetic helicity from an active region.


In order to allow finite magnetic fluxes through the boundaries, a Coulomb gauge is constructed that allows for global magnetic flux balance. The property of sinusoidal function is used to obtain the vector potentials at the twelve edges of the considered rectangular volume extending above an active region. We tested and verified our method in a theoretical fore-free magnetic field model.


We apply the new method to the former calculation data and found a difference of less than 1.2%. We also applied our method to the magnetic field above active region NOAA 11429 obtained by a new photospheric-data-driven MHD model code GOEMHD3.  We analyzed the magnetic helicity evolution in the solar corona using our new method. It was found that the normalized magnetic helicity (H_/Φ^2) is equal to -0.038 when fast magnetic reconnection is triggered.  This value is comparable to the previous value (-0.029) in the MHD simulations when magnetic reconnection happened and the observed normalized magnetic helicity (-0.036) from the eruption of newly emerging active regions. We found that only 8% of the accumulated magnetic helicity is dissipated after it is injected through the bottom boundary. This is in accordance with the Woltjer conjecture. Only 2% of magnetic helicity injected from the bottom boundary escapes through the corona. This is consistent with the observation of magnetic clouds, which could take away  magnetic helicity into the interplanetary space, in the case considered here, several halo CMEs and two X-class solar flares origin from this active region.






Yang  et al.
Evolution of RMH: New Boundary Conditions for the Vector Potential

Evolution of Relative Magnetic Helicity: New Boundary Conditions for the Vector Potential
    
Shangbin Yang 1,2,3
Jörg Büchner 2
Jan Skála 2,4
Hongqi Zhang   1 
    December 30, 2023
=========================================================================================






§ INTRODUCTION

Magnetic helicity is a key geometrical parameter to describe the
structure of solar coronal magnetic fields <cit.>. Magnetic helicity in a volume V can be determined as

    H_ M=∫_V𝐀·𝐁dV,

where A is a vector potential and B is  magnetic field
in this volume. Magnetic helicity H_ M is conserved in an ideal
magneto-plasma <cit.>. 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 <cit.>. 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·n̂|_S = 0. On the other hand, Berger and Field (1984) have
shown that in the case of boundaries open to magnetic flux
penetration, instead, relative magnetic helicity (H_ R),
given by the Finn-Antonsen (1985) formula


    H_ R = ∫( 𝐀 + 𝐀_ p)
    ·( 𝐁 - 𝐏) dV ,

is gauge-invariant if the magnetic field is chosen as a
reference field such that 𝐏· n|_S=𝐁·
n|_S. It is customary to choose potential field (∇×𝐏=0) as the reference
field.


The concept of magnetic helicity has successfully been applied to characterize solar coronal process as well as the solar dynamo to interpret solar observations and corona simulations. <cit.> obtained the relative helicity by applying Fourier transform (FFT) method and local correlation tracking (LCT) to MDI data. Since then the helicity calculation method has been developed further <cit.> by improving the helicity density map and velocity tracking techniques. The analysis of magnetic helicity based on observations and simulations has also been developed such as in the analysis of magnetic
helicity injection in the course of flux emergence <cit.>,the correlation between helicity
change and solar eruption <cit.>, magnetic energy-helicity relation analysis in solar eruptions by <cit.>, magnetic helicity distribution in different scales in the solar dynamo <cit.>, magnetic helicity estimation in the interplanetary magnetic cloud
<cit.>, testing the conjecture of Taylor conjecture in quasi-ideal MHD simulations <cit.>. For a recent review of modeling and observations of  magnetic helicity see, e.g., <cit.>.








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 neighboring emerging active regions <cit.>. Thus the correlation study between magnetic helicity distribution and evolution and solar eruption become important.

<cit.> 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. <cit.> found that magnetic helicity of CME productive ARs is higher than other ARs. The survey of <cit.> to helicity accumulation in 393 ARs a 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^36Mx^2/s.










In the simulations the debates also exist for the possible upper-bound helicity before the solar eruption happens <cit.>. In these papers, the corresponding normal relative magnetic helicity (H_/Φ^2) when the eruption happens reaches approximately -0.16 for case of Kink instability and approximately -0.18 for case torus instability in the simulations of <cit.>. However, from the observations the range of normalized helicity for the active region is from 0.02 to 0.08 <cit.>. This is  one order of magnitude smaller value than the above simulations. <cit.> investigate the value of normalized helicity reached 0.0298 just prior to drastic energy release by magnetic reconnection by using the magnetic field above active region NOAA 8210 obtained by a photospheric-data-driven MHD model <cit.>. <cit.> 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 find that the consequences of 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 <cit.>, which is close to the simulation of NOAA 8210 and the theoretical prediction value of <cit.> only in a multipolar force-free magnetic field structure. However, whether the relative magnetic helicity  has a upper-bound and how it plays a important role in solar eruptions is still an open question.


Despite of 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, Inhester, and Wiegelmann,
2011; Rudenko and Myshyakov, 2011). <cit.> developed a method for calculating the relative magnetic helicity in a finite 3D volume which was applied to a simulated flaring AR 8210 <cit.>. However, this method required that
the magnetic flux is balanced  on each of the side boundaries
in the considered volume. In this paper, the corresponding scheme is presented which does not require such
restriction on the flux balance. Such method has already been applied in <cit.> and it is the most accurate one among 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 <cit.>.
In Sec. <ref>, we describe
the restriction of vector potential in the previous paper. In
Sec. <ref>, we present the  new scheme to
calculate the vector potentials on the six boundaries. In
sec. <ref>, we use a Non-Linear Force-Free magnetic field model <cit.>  and MHD simulation model to check our scheme. The  summary and some discussions are given in
Sec. <ref>.








§ RESTRICTED METHOD TO OBTAIN A_ P AND A AT THE BOUNDARIES


The computation of H_R of 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 Section 2.2 of <cit.> for details.
Let us define a finite three-dimensional (3-D) 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, one has to provide 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_ p x  =  - ∂φ/∂ y,     
    A_ p  y  = ∂φ/∂ x,     .A_ p  z|_z = 0  = 0.



According to the definition of the vector potential, the scalar
function φ (x, y) satisfies the Poisson equation:

    Δφ(x,y)  = B_z (x,y,z = 0).

In our previous work <cit.> we set the values of
∂φ/∂ n at the four edges of the
plane z=0 to zero in Equation (<ref>).
 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 other five 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, one should
calculate the values of vector potentials at the twelve edges of the
three-dimensional volume to obtain a Neumann boundary condition for
the Poisson equation (<ref>) at each side boundary.


In next section, we will introduce a scheme to calculate the vector
potentials at the twelve edges.


§ GENERAL METHOD TO OBTAIN 𝐀_P AND 𝐀 AT THE
BOUNDARIES
 
In order to determine 𝐀_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· d l at the twelve edges are defined as
a_i (i=1,...,12). The twelve integrals and the corresponding
directions are depicted in Fig. <ref>.

Applying Ampère's law to the six boundaries independently, one obtain a linear system of equations for obtaining the integral value a_i :

    T X =  F,

where F=(Φ_1,Φ_2,Φ_3,Φ_4,Φ_5,Φ_6)^T,
X=(a_1,a_2,a_3,...,a_12)^T and T̂ is a 6×12 matrix given by

    T̂=
    [
    [  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;    ]].

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
twelve integrals a_i.  We use the freedom in the gauge to construct twelve 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 T̂'̂ to make sure the determinant of T̂'̂ is not zero without limiting
the validity of the solution for a_i. we choose the following seven rows  to construct the new matrix  T̂'̂:

    T̂'̂=
    [
    [  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;    ]].

This choice for the lower seven rows of T̂'̂ is not unique, but any influence on the gauge will be removed later.
One can calculate that  the determinant of T̂'̂ does not
vanish. According to Cramer's rule, an unique solution of
the new linear equation

    T̂'̂ X = F',

can exists if X=(a_1,a_2,a_3,...,a_12)^T and F'=(Φ_1,Φ_2,Φ_3,Φ_4,Φ_5,0,0,0,0,0,0,0)^T.
Then we can obtain the integrals of a_i at the twelve edges such that Gauss 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 twelve edges can be obtained representing the components of A_ p in the following way:

    A_px( a_i ) = π a_i /2L_x sin (π x/
     .
     -L_x ),i = 1,3,5,7 
    A_py( a_i ) = π a_i /2L_y sin (π y/
     .
     -L_y ),i = 2,4,6,8 
    A_pz( a_i ) = π a_i /2L_z sin (π z/
     .
     -L_z ),i = 9,10,11,12


Note 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. (<ref>).  Now we solve the Poisson equations (<ref>) 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 Sec. 2.2 and Sec 2.3 of <cit.>.  Using the boundary conditions above, 𝐀_p and 𝐀 in the volume are obtained by solving the following Poisson and Laplace equations:

    Δ𝐀_ p = 0

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

    Δ𝐀=-𝐉,

where J = μ_0  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:

    ∇×( ∇×𝐌) = 𝐏
    - ∇×𝐀_ p.

 The components of M satisfy the three Poisson equations:

    [ {[ ΔM_z = ( ∇×A⃗_p)_z - P_z;      M_z( z = 0,l_z) = 0; ∂M_z/∂ x( x = 0,l_x) = 0; ∂M_z/∂ y( y = 0,l_y) = 0 ]., {[ ΔM_y = ( ∇×A⃗_p)_y - P_y;      M_y( y = 0,l_y) = 0; ∂M_y/∂ x( x = 0,l_x) = 0; ∂M_y/∂ z( z = 0,l_z) = 0 ].,; {[ ΔM_x = ( ∇×A⃗_p)_x - P_x;      M_x( x = 0,l_x) = 0; ∂M_x/∂ y( y = 0,l_y) = 0; ∂M_x/∂ z( z = 0,l_z) = 0 ]..                                                                                                               ]

Finally, in order to remove the residual errors in the solenoidal property of A_
p, we introduce  scalar
field ϕ (x, y, z) which satisfies the following Poisson
equation:

    {[                      Δϕ  =  - ∇·𝐀_ p; . ∂ϕ/∂n|_s  = .  - ( ∇×𝐌) ·3pt⌢ n|_s;                                      ]..

After solving Eq. <ref>, one obtains a new modified vector potential A_
p' is represented as

    𝐀'_ p  = 𝐀_ p  +∇×𝐌 +
    ∇ϕ.

and inn the same way, one can obtain the corrected vector potential A' just by replacing the right-hand-side term of Eq. (<ref>) as
∇× A -  B. Now one can calculate the relative helicity in the
volume according to Eq. (<ref>). As already mentioned when introducing T̂'̂, the choice of the lower seven rows of T̂'̂ and the final seven elements of F' is not unique. 


§ TESTING THE SCHEME

 For testing the new scheme to obtain the vector
potentials at the boundaries, we use analytical nonlinear
force-free fields of <cit.>. We utilized the model labeled
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 calculate the magnetic fluxes Φ_0 through
the six boundaries and substitute it to equation (<ref>) to obtain the integrals
a_i along the twelve edges of the 3D volume. Then we substitute a_i
into equation (<ref>) respectively to obtain the boundary value for
solving the Poisson equation (<ref>) at the
six boundaries. After we obtained A_ p on the six
boundaries, we calculate the magnetic fluxes
 Φ from the computed vector potential using
 the relation between the vector potential and the
magnetic field: B·n̂=∇× A_
p·n̂. Table. <ref> represents the
quantitative result after applying the above scheme. As one can see in the table the calculated magnetic fluxes through the six boundaries by using our new scheme is actually better flux-balanced than the theoretical model. Such small errors are partly due to the fact that the total magnetic flux of the
analytical model does not completely vanish. It is required that the total magnetic flux must be zero for resolving the linear equation (<ref>). Numerical errors in solving the Poisson equation on the other hand will also introduce finite total magnetic flux as well. The calculated relative magnetic helicity using the previous method is -14441.45 and the value is -13253.36 using the new method. There is 8% difference for the two methods.




We then applied our new scheme to re-calculate the relative magnetic helicity above active
region NOAA AR8210 in <cit.>. Fig. <ref> depicts the comparison of relative magnetic helicity evolution by using the previous and the new scheme.  The solid (dashed) line represents the new (previous) relative magnetic helicity evolution.  The
difference is very small (<1.2%). This is due to the magnetic field structure of the corona above AR8210 obtained by <cit.> in which the fluxes through the side boundaries were almost balanced.

Finally, we apply our new scheme to the 3D magnetic field data obtained from simulated evolution of the solar corona above active region 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 solving a second-order-accurate MHD equations <cit.>. It was successfully tested and applied to study the magnetic coupling between
the solar photosphere and corona based on multi-wavelength 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 fulfill 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 programing paradigm, adopting a standard two-dimensional domain-decomposition approach. This allows 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.








Fig.<ref> shows the evolution of relative magnetic helicity by using our
method to the simulation data from the application of GOEMHD3 to the evolution of the AR 11429 on March 07. 2012 <cit.>. The vertical  dash-dot line in the Fig.<ref> indicate the time (t=1035s) when fast reconnection started. Fig.<ref>a 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 <cit.> is extended to all six boundaries and the dissipative volume integral in Eq. (6) of <cit.> is also included in the computation.
Fig.<ref>b depicts the normalized magnetic helicity evolution (H_/Φ^2). We use two types of fluxes Φ as the normalized parameter. One is unsigned magnetic flux only through the bottom boundary (solid line in Fig. <ref>b), which is usually used in observations. The other is unsigned magnetic flux through the six side boundaries (dashed line in Fig. <ref>b). A value of -0.038 (-0.025)  is reached after 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 the magnetic flux from the observations of newly emerging active regions <cit.>:

    logΔ H/H_0 =alogΦ _m /Φ_0+b,

where a=1.85, b=-0.41, H_0 = 10^41 Mx^2, and Φ_0 =
10^21 Mx.
, which describes the relation between the accumulated helicity Δ H
and magnetic flux Φ _m for active regions. Fig.<ref>c depicts the magnetic helicity flux through the bottom boundary into the solar corona. The total magnetic helicity flux is shown in Fig.<ref>d. It is found that negative magnetic helicity is injected into the Corona due to the sub-photospheric plasma motion. The total injected helicity is
-8.5×10^42Mx^2. Fig.<ref>e 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.<ref>f. It is found 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^42Mx^2. 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
by <cit.>.
Fig.<ref>g depicts the magnetic helicity flux which escape through the upper five boundaries. The total escaped magnetic helicity evolution is shown in
Fig.<ref>h. It is found that only -0.16×10^42 Mx^2 of magnetic helicity flux escapes through the upper boundaries. This is
23% of the dissipated helicity flux (0.71×10^42Mx^2, see
Fig.<ref>f) and only 2% of the total injected magnetic helicity ( -8.5×10^42Mx^2, see Fig.<ref>d).  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 <cit.>. In the calculation of Fig.<ref> c-g at the separate boundary, we use the Coulomb gauge.



§ 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 <cit.>. We verified the new method using a analytical theoretical force-free model magnetic field  <cit.> and to the simulated data in <cit.>. We also apply the new method to simulated coronal
magnetic fields using the newly developed GOEMHD3 to investigate the magnetic helicity evolution of solar corona in the course of evolution of AR11429.

We found 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 <cit.>. Only 2% of magnetic helicity injected through the bottom boundary escapes to the solar wind. This shows that the magnetic helicity cannot efficiently escape.

This may help to understand the lack of magnetic clouds [Near-Earth Interplanetary Coronal Mass Ejections Since January 1996: www.srl.caltech.edu/ACE/ASC/DATA/level3/icmetable2.htm. It is described in <cit.>.] with considerable magnetic helicity in the interplanetary space <cit.>. 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_/Φ^2) 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 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 <cit.>. 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 <cit.> when eruption happened. This might be due to the anomalous resistivity caused by the strong micro-turbulence and microscopic structures <cit.> used in the GOEMHD3
model. The anomalous current dissipations allows essentially increase the dissipation rate of magnetic energy following a global MHD instability (Büechner et al. 2017, in preparation).

<cit.> as well as <cit.> 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 obtained value 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 <cit.>. The upper-bound normalized magnetic helicity value deviates evidently 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 <cit.>.  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 the future study.

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 three-dimensional magnetic structure to analyze the evolution of relative magnetic helicity.






We would like to thank the referee for carefully reading our
manuscript and for giving such constructive comments which
substantially helped improving the paper. This study is supported by grants 10733020, 10921303,
41174153,11173033 11178016 and 11573037 of National Natural Science
Foundation of China, 2011CB811400 of National Basic Research Program
of China, a sandwich-PhD grant of the Max-Planck Society and 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.The
authors also like to thank the Supercomputing Center of Chinese
Academy of Sciences (SCCAS) and Max Planck Computing and Data Facility (MPCDF)
for the allocation of computing time.





 [Berger & Field(1984)]BER84 Berger, M. A.,& Field, G. B. 1984, J. Fluid Mech., 147, 133

[Berger(1999)]ber99Berger, M. A.: 1999, Plasma Phys. Contr. Fusion, 41, 167

[Büchner et al.(2004)]bue04Büchner, J., Nikutowski, B., & Otto, A. 2004, in Multi-Wavelength
Investigations of Solar Activity, eds. A. V. Stepanov, E. E. Benevolenskaya,
& A. G. Kosovichev, IAU Symp., 223, 353

[Büchner & Elkina(2006)]bue06bBüchner, J., Elkina, N. 2006, Phys. Plasmas,13, 082304.

[Cane & Richardson(2003)]RC03Cane, H. V., & Richardson, I. G. 2003, JGR, 108, A4.

[Chae (2001)]Cha01 Chae, J. 2001, ApJ, 560, L95

[Démoulin & Berger(2003)]Dem03
Démoulin, P., & Berger, M. A. 2003, Solar Phys., 215, 203.

[Démoulin & Pariat(2009)]Dem09Démoulin, P., Pariat, E.: 2009, Adv. Space Res., 43, 1013

[Démoulin et al.(2015)]Dem15Demoulin, P.; Janvier, M.; Dasso, S. 2015, Solar physics, in press

[Fan (2005)]Fan05Fan, Y. 2005, ApJ, 630, 543.

[Fan & Gibson (2007)]Fan07
Fan, Y., & Gibson, S. E. 2007, ApJ, 668, 1232

[Jacobs et al.(2006)]Jac06Jacobs, C., Poedts, S., van der Holst, B. 2006, A&A, 450,793
[Amari et al.(2004)]Amari04Amari, T., Luciani, J. F., Aly, J. J., Mikic, Z., Linker, J. 2003,ApJ,585,1073

[Jing et al.(2012)]Jing12Jing, J., Park, S.,  Liu, C., Lee, J., Wiegelmann, T., Xu, Y., Deng, N.,& Wang, H. M. 2012, , 752, L9

[Labonte et al.(2007)]Lan07B. J. LaBonte,1 M. K. Georgoulis, & D. M. Rust 2007,
ApJ, 671,955

[Low & Lou(1990)]low90 Low, B. C.,& Lou, Y.Q. 1990, , 352, 343

[Liu et al.(2014)]liu14Liu, Y., Hoeksema, J. T., Bobra, M., Hayashi, K., Schuck, P. W., Sun, X. 2014, 785, 13L

[Nindos & Andrews (2005)]Nin05Nindos, A., Andrews, M. D. 2005, IAUS, 226, 194N

[Park et al.(2008)]Park08
Park, S., Lee, J., Choe, G. S., Chae, J., Jeong, H., Guo, Y., Jing,
J., & Wang, H. M. 2008, ApJ, 686, 1397

[Richardson & Cane(2010)]RC10Richardson, I. G., & Cane H. V., 2010, Solar physics 264, 189.

[Rudenko & Myshyakov(2011)]Ruden11Rudenko, G. V., & Myshyakov, I. I. 2011, Solar phys., 270, 165

[Thalmann, Inhester, Wiegelmann(2011)]TIW11Thalmann, J. K., Inhester, B., & Wiegelmann, T. 2011, Solar Phys., 272, 243

[Valori, Démoulin, Pariat(2012)]VDP12Valori, G., Démoulin, P., Pariat, E. 2012, Solar Phys., 278, 347
[Valori et al.(2016)]Val16Valori, G., Pariat, E., Anfinogentov, S., Chen,F., Georgoulis, M. K., Guo, Y., Liu, Y., Moraitis, K.,Thalmann, J. K., & Yang, S. 2016, Space Science Review, 201,147

[Woltjer(1958)]wol58Woltjer, L. 1958, Proc. Natl Acad. Sci. USA, 44, 480

[Pariat, Démoulin, & Berger(2005)]Par05
Pariat, E., Démoulin, P., & Berger, M. A. 2005, A&A, 439, 1191
[Pariat et al.(2015)]par15Pariat, E.,Valori, G., Démoulin, P., & Dalmasse, K. 2015, A&A, 580, A128
[Pipin & Pevtsov (2014)]Pipin14Pipin, V. V., & Pevtsov, A. A. 2014, ApJ, 789, 21
[Santos et al.(2011)]Santos11Santos, J. C., Büchner, J., & Otto, A. 2011,

[Schuck (2008)]Sch08Schuck, P. W.	2008, ApJ, 683, 1134

[Seehafer et al.(2003)]Seehafer03
Seehafer, N., Gellert, M., Kuzanyan, K. M., Pipin, V. V. 2003, AdSpR 32, 1819.


[Skala et al.(2015)]Skala15Skála, J., Baruffa, F. , J. Büchner, & RamppA M. 2015 A&A, 580, A48.
[Tian & Alexander(2008)]Tian08Tian, Lirong,& Alexander, David 2008 ApJ, 673, 532


[Tziotizou et al.(2012)]Tzi12Tziotziou, Kostas; Georgoulis, Manolis K.; Raouafi, Nour-Eddine 2012, ApJ, 759, L4

[Taylor (1974)]Tay74 Taylor, J. B. 1974, Phys. Rev. Lett., 33, 1139

[Yang et al.(2009)]Yang09aYang, S., Büchner, J., Zhang, H. 2009, ApJ, L25.

[Yang, Zhang, & Büchner(2009b)]Yang09b Yang, S., Zhang, H., & Büchner, J. 2009, A&A, 502, 333.

[Yang et al.(2013)]Yang13Yang, S., Büchner, J. , Santos, J. C., & Zhang, H. 2013, Solar Physics, 283, 369.
[Yang et al. (2015)]Yang15Yang, S., Xie, W., Liu, J. 2015, AdSpR, 55, 1553.

[Zhang & Flyer(2008)]Zhang08Zhang, M., & Flyer, N. 2008, ApJ, 683.

[Zhang et al.(2006)]Zhang06Zhang, M., Flyer, M., & Low, B. 2006,ApJ, 644,575-586

[Zhang et al.(2008)]ZhangY08
Zhang, Y., Tan, B. L., & Yan, Y. H. 2008, , 682, L133



