1 Introduction

Magnetic helicity is a key geometrical parameter to describe the structure and evolution of solar coronal magnetic fields (e.g. Berger 1999). The concept of helicity has successfully been applied to basic solar processes such as the solar dynamo (e.g. Seehafer et al. 2003; Kuzanyan et al. 2003), flares (e.g. Zhang 2006), filaments (Devore 2000), coronal mass ejections (CMEs; Zhang, Flyer, and Low 2006; Kliem, Rust, and Seehafer 2011), interplanetary coronal mass ejections (ICMEs; Nakwacki et al. 2011), and also to the interaction between active regions (Yang, Büchner, and Zhang 2009). A recent review about modeling and observations of photospheric magnetic helicity can be found in Démoulin and Pariat (2009).

Magnetic helicity in a volume V can be determined as

$$ H_\mathrm{M}=\int_{V}\mathbf{A}\cdot\mathbf{B}\,\mathrm{d}V, $$
(1)

where A is the vector potential for the magnetic field B in this volume. Magnetic helicity is conserved in an ideal magneto-plasma (Chae et al. 2001). As long as the overall magnetic Reynolds number is high, however, it is still approximately conserved, even in the course of relatively slow magnetic reconnection (Berger and Field 1984). 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, Inhester, and Wiegelmann 2011; Rudenko and Myshyakov 2011).

The problem with the helicity defined by Equation (1) is that the vector potential A in this formula is not directly measurable, nor is it obtainable from general numerical simulations unless they are carried out in terms of vector potential instead of electromagnetic fields. Even more – since the vector potential A is not gauge-invariant, the magnetic helicity in a closed volume has a well-determined value only when the magnetic field at the boundary have only tangential components, i.e., if \({\mathbf{B}} \cdot \hat{\mathbf{n}}|_{S} = 0\). On the other hand, Berger and Field (1984) has shown that in the case of boundaries that are instead open to magnetic flux penetration, the relative magnetic helicity (H R), given by the Finn–Antonsen (1985) formula

$$ H_\mathrm{R} = \int {( {\mathbf{A} + \mathbf{A}_\mathrm{p} } ) \cdot ( {\mathbf{B} - \mathbf{P}} )} \,\mathrm{d}V , $$
(2)

is gauge-invariant if the potential magnetic field is chosen as a reference field. In Equation (2) the quantities A and B are the actual vector potential and the original magnetic field, respectively, while P is the corresponding potential (i.e., current-free) magnetic field (∇×P=0) and A p is its vector potential. The vector potentials A and A p can be chosen to have the same gauge. For the correct calculation of H R the corresponding potential magnetic field P has to satisfy the conditions

$$ {\nabla\times\mathbf{P}=0 }, $$
(3)
$$ \mathbf{P}\cdot\hat{\mathrm{n}}|_{S}=\mathbf{B}\cdot\hat{\mathrm{n}}|_{S}, $$
(4)

i.e., the normal component of the potential field P is equal to that of B at the boundaries. Here \(\hat{\mathbf{n}}|_{S}\) denotes the unit vector normal to the boundary S.

In ideal magnetohydrodynamics (MHD) the rate of change of relative magnetic helicity in a given finite volume is given by the flux of relative helicity across the boundary S of this volume (e.g. Berger and Field 1984; Berger 1999):

$$ {\frac{\mathrm{d}H_\mathrm{R}}{\mathrm{d}t}= -2\int_{S}\bigl( (\mathbf{A}_\mathrm{p}\cdot\mathbf{V})\mathbf{B} -(\mathbf{A}_\mathrm{p}\cdot\mathbf{B})\mathbf{V}\bigr) \cdot \hat{\mathrm{n}}|_{S}}\, \mathrm{d}S, $$
(5)

where V is the plasma velocity at the boundary and \(\mathbf{A}_{\mathrm{p}}\cdot\widehat{n}|_{S}=0\).

In a non-ideal plasma system, in contrast, relative magnetic helicity changes not only by the helicity flow through the boundary, but also by internal dissipation. The latter can be determined from E⋅B (Berger and Field 1984), which adds a term to the evolution equation of relative helicity

$$ \frac{{\mathrm{d}H_\mathrm{R} }}{{\mathrm{d}t}} = - 2\int_{V} ( {\mathbf{E} \cdot \mathbf{B} } )\mathrm{d}V -2\int_{S} \bigl( (\mathbf{A}_\mathrm{p}\cdot\mathbf{V})\mathbf{B} -(\mathbf{A}_\mathrm{p}\cdot\mathbf{B})\mathbf{V}\bigr) \cdot \mathbf{n} \,\mathrm{d}S. $$
(6)

To calculate the evolution of the relative helicity H R of a magnetic field B in a volume, one first has to calculate the appropriate vector potential A, the corresponding potential field P and its vector potential A p and, to account for internal dissipation processes, the electric field E. Note that Equations (5) and (6) can be derived from Equation (2) if the vector potential satisfies the Coulomb gauge and the normal component of A p vanishes at the boundary (Berger and Field 1984).

A rigorous theory of vector potential in a half-space was developed by Boulmezaoud (1999). Devore (2000) suggested a simplified method for calculating the vector potential to obtain the relative magnetic helicity in a half-space. In most practical cases, however, the volume under consideration is limited in all three spatial directions. Recently, Thalmann, Inhester, and Wiegelmann (2011) presented a method for calculating the vector potential of a solenoidal magnetic field by decomposing it to the sum of a Laplacian part and a current-carrying part according to the method of Kusano, Suzuki, and Nishikawa (1995). Rudenko and Myshyakov (2011) introduced an algorithm for calculating the gauge-invariant helicity of a magnetic field in a rectangular box. Valori, Démoulin, and Pariat (2012) employed a gauge different from that of Thalmann, Inhester, and Wiegelmann (2011) and Rudenko and Myshyakov (2011), which greatly simplified the computation of helicity in a finite volume. However, all these studies did not consider the evolution of relative magnetic helicity by comparing the accumulated relative magnetic helicity and the helicity flux through the boundary. In fact, the dissipation of relative magnetic helicity inside the considered coronal volume becomes important if non-ideal plasma processes take place, such as fast magnetic reconnection. The existing methods for calculating the relative magnetic helicity in a finite volume were tested on certain magnetic field structures such as an analytical magnetic field (Rudenko and Myshyakov 2011; Thalmann, Inhester, and Wiegelmann 2011) and theoretical models (Valori, Démoulin, and Pariat 2012).

Our goal is to present a method that allows calculating the relative magnetic helicity according to Equation (2). We apply our method to an analytical field and to simulated data including both ideal and non-ideal coronal plasma processes above an active region. We also investigate the inward/outward flux of the relative magnetic helicity through the boundaries, according to Equation (5) for ideal plasmas and the loss of helicity according to Equation (6) for non-ideal plasma processes to check whether the calculation of relative magnetic helicity in a finite volume coincides with the helicity flux across the boundary.

In Section 2 we describe our new method for calculating vector potentials needed to determine relative magnetic helicity and its flux through the boundaries. In Section 3 we use the theoretical model of Low and Lou (1990) to check our method. Then, in Section 4, we apply this method to data obtained by the numerical simulation of the field and plasma dynamics above an observed active region on the Sun. The summary and discussion of our results are given in Section 5.

2 Calculating the Relative Magnetic Helicity in a Finite Volume and Its Flux through the Boundaries

Let us define a finite three-dimensional (3D) volume (“box”) in Cartesian coordinates with a magnetic field B(x,y,z) given in this volume. Let the volume be bounded by x=[0,l x ], y=[0,l y ], and z=[0,l z ].

2.1 Determining the Boundary Values of A p and A

First one has to provide the values of A p and A on all six boundaries (x=0,x=l x ;y=0,y=l y ;z=0,z=l z ). To take the bottom boundary (z=0), for example, we define a new scalar function φ(x,y) that determines the vector potential A p of the potential magnetic field P on this boundary as follows:

$$ A_{\mathrm{p} x} = - \frac{{\partial \varphi }}{{\partial y}}, \qquad A_{\mathrm{p} y} = \frac{{\partial \varphi }}{{\partial x}}, \qquad {A_{\mathrm{p} z} } |_{z = 0} = 0. $$
(7)

According to the condition imposed on the reference potential field (Equation (4)), the scalar function φ(x,y) should satisfy the Poisson equation:

$$ \Delta \varphi(x,y) = B_z (x,y,z = 0). $$
(8)

The value of ∂ϕ/∂n on the four sides of the plane z=0 is set to zero in Equation (8). According to Equation (7), A px and A py will vanish at y=0,l y and at x=0,l x , respectively, on the z=0 plane. Therefore, the corresponding magnetic flux at the boundary should also vanish because of Ampère’s law. In the examples given in Sections 3 and 4, the magnetic fluxes at the six boundaries are all close to zero and satisfy the above requirement for magnetic flux. We use a fourth-order accurate fast direct scheme according to Boisvert (1984), which is included, e.g., in the IMSL Numerical Libraries, to solve the Poisson equation (Equation (8)). By solving Equation (8) on the boundary (z=0), one can obtain the values of A p. The values of A p on the other five boundaries can be obtained in a similar way. For the vector potential A at all boundaries the same values are taken as for A p.

2.2 Determining A p and A in the Volume

According to the definition of the potential field (Equation (3)), the vector potential A p of P satisfies the Laplace equation

$$ \Delta \mathbf{A}_\mathrm{p} = 0 $$
(9)

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 \mathbf{A}=-\mathbf{J}, $$
(10)

where J=μ 0 j denotes the current density in the plasma system if A satisfies the Coulomb gauge. Again using the method of Boisvert (1984), we solve the Laplace and Poisson equations (Equations (9) and (10)) with the boundary conditions given in Section 2.1. The potential field P in this volume is defined as ∇ψ, which should satisfy the following Poisson equation

(11)

which can be solved using the same method.

2.3 Clean the Divergence of A p and A

In calculating A p, however, finite divergence values may arise; the method cannot ensure the divergence of A p to vanish under the boundary condition (Equation (8)). To clean the residual divergence and obtain the correct solution of Equation (9), we first introduce a solenoidal modification vector ∇×M that satisfies the condition

$$ \nabla \times ( {\nabla \times \mathbf{M}} ) = \mathbf{P} - \nabla \times \mathbf{A}_\mathrm{p}, $$
(12)

in which the potential field P and the corresponding vector potential A p are obtained as described in Section 2.2. The components of M satisfy the three Poisson equations

$$ \begin{array}{l} \left\{ \begin{array}{l} \Delta M_z = ( {\nabla \times \mathbf{A}_\mathrm{p} } )_z - P_z , \\ \noalign {\vspace {3pt}} M_z ( {z = 0,l_z } ) = 0, \\ \noalign {\vspace {3pt}} \frac{{\partial M_z }}{{\partial x}}( {x = 0,l_x } ) = 0, \\ \noalign {\vspace {3pt}} \frac{{\partial M_z }}{{\partial y}}( {y = 0,l_y } ) = 0, \\ \end{array} \right.\qquad \left\{ \begin{array}{l} \Delta M_y = ( {\nabla \times \mathbf{A}_\mathrm{p} } )_y - P_y , \\ \noalign {\vspace {3pt}} M_y ( {y = 0,l_y } ) = 0, \\ \noalign {\vspace {3pt}} \frac{{\partial M_y }}{{\partial x}}( {x = 0,l_x } ) = 0, \\ \noalign {\vspace {3pt}} \frac{{\partial M_y }}{{\partial z}}( {z = 0,l_z } ) = 0, \\ \end{array} \right.\quad \\ \noalign {\vspace {5pt}} \left\{ \begin{array}{l} \Delta M_x = ( {\nabla \times \mathbf{A}_\mathrm{p} } )_x - P_x , \\ \noalign {\vspace {3pt}} M_x ( {x = 0,l_x } ) = 0, \\ \noalign {\vspace {3pt}} \frac{{\partial M_x }}{{\partial y}}( {y = 0,l_y } ) = 0, \\ \noalign {\vspace {3pt}} \frac{{\partial M_x }}{{\partial z}}( {z = 0,l_z } ) = 0. \\ \end{array} \right. \end{array} $$
(13)

These Poisson equations can be solved with the same method as described in Section 2.2. The magnetic field from the modification vector ∇×(∇×M) will vanish at the six boundaries, so that it does not change the normal component of the magnetic field \({\mathbf{B}}\cdot\hat{\mathbf{n}}|_{S}\).

Let us consider the z-component of ∇×M on the boundary of z=0 or z=l z as an example. We know that

$$ \nabla \times \mathbf{M} = \hat{\mathbf{i}} \biggl( {\frac{{\partial M_z }}{{\partial y}} -\frac{{\partial M_y }}{{\partial z}}} \biggr) + \hat{\mathbf{j}} \biggl( {\frac{{\partial M_x }}{{\partial z}} -\frac{{\partial M_z }}{{\partial x}}} \biggr) + \hat{\mathbf{k}} \biggl( {\frac{{\partial M_y }}{{\partial x}} -\frac{{\partial M_x }}{{\partial y}}} \biggr) . $$
(14)

Because M z =0 at z=0 or z=l z , the z-component of ∇×M on the same boundary can be simplified as

$$ \nabla \times {( {\nabla \times \mathbf{M}} )_z} |_{z = 0,l_z } = {\frac{\partial }{{\partial x}} \biggl({\frac{{\partial M_x }}{{\partial z}}} \biggr)} \bigg|_{z = 0,l_z } + {\frac{\partial }{{\partial y}}\biggl( {\frac{{\partial M_y}}{{\partial z}}} \biggr)} \bigg|_{z = 0,lz}. $$
(15)

The quantities \(\frac{\partial M_{y}}{\partial z}\) and \(\frac{\partial M_{x}}{\partial z}\) also vanish on the boundary z=0 or z=l z as required by Equation (13). Hence, the z-component of ∇×M vanishes on the boundary z=0 or z=l z . In the same way, the x- and y-components of ∇×M also vanish on the boundary y=0,l y and x=0,l x , respectively. That is to say, the normal component of the magnetic field does not change on the boundary after adding ∇×M to the vector potential A p.

Note that ∇×M is divergence-free, although A p is not. However, the vector potential A p cannot vanish on the six boundaries after adding ∇×M to it. Therefore, we introduce another scalar field ϕ(x,y,z) which satisfies the Poisson equation

(16)

This Poisson equation can be solved with the same method as described in Section 2.2. The final form for the modified vector potential \({\mathbf{A}}_{\mathrm{p}}'\) is given as

$$ \mathbf{A}'_\mathrm{p} = \mathbf{A}_\mathrm{p} +\nabla \times \mathbf{M} + \nabla \phi. $$
(17)

In the same way, one can obtain the corrected vector potential A′ by replacing the right-hand-side term of Equation (13) with ∇×A−B. One can confirm that the final modified vector potential \({\mathbf{A}}'_{\mathrm{p}}\) and A′ satisfy the Coulomb gauge; they are also the solutions to Equations (9) and (10) given in Section 2.2. The normal component of the fields \(\nabla\times \mathbf{A}'_{\mathrm{p}}\) and ∇×A′ is equal to that of B on the six boundaries, as required by Equation (4). Now one can calculate the relative helicity in the volume according to Equation (2) as well as the change rates of relative magnetic helicity for an ideal plasma evolution according to Equation (5) and for a non-ideal evolution according to Equation (6), if additionally the electric fields and the velocity on S are known.

3 Testing the Method

To test the method we use the axisymmetric nonlinear force-free fields of Low and Lou (1990). We used the model labeled P 1,1 with l=0.3 and Φ=π/4 in the notation of their paper. We calculated the magnetic field on a grid of 64×64×64. Figure 1 depicts the original magnetic field B, the reconstructed magnetic field ∇×A calculated according to our method, and the corresponding errors at the bottom boundary. We found that the reconstructed magnetic field is very similar to the original magnetic field. The ratio in magnetic energy between the two magnetic fields (∑ i |B|2/∑ i |∇×A|2) is 0.998, and the corresponding vector correlation (Schrijver et al. 2006, Equation (28)) for them is 0.980.

Figure 1
figure 1

Comparison between B (first row) and ∇×A (second row) at the bottom boundary z=0. The third row shows the difference between them. The left, middle, and right columns show the x-, y-, and z-components of the magnetic field, respectively.

Full size image

We also used the fractional flux error f i in a small discrete volume ΔV i (with surface ΔS i ) over all grid points to check the divergence of A and A p (Wheatland, Sturrock, and Roumeliotis 2000),

$$ f_i = \frac{{\int_{\Delta S_i } {\mathbf{A} \cdot \mathrm{d}\mathbf{S}} }}{{\int_{\Delta S_i } {| {\mathbf{A}} |\mathrm{d}S} }} \approx \frac{{\biggl( {\nabla \cdot \mathbf{A}} \biggr)_i \Delta V_i }}{{A_i S_i }}, $$
(18)

where S i is the surface area of a small volume V i . The average values of the magnitude of f i (〈|f i |〉) for A and A p are 6.69×10−5 and 6.12×10−5, respectively. Before cleaning the divergence of the vector potential, the average value for A is 0.0055 and for A p it is 0.0038. Figure 2 shows the distribution of the fractional flux error over all grid points after cleaning the divergence. It indicates that A p and A are close to divergence-free, meaning that our method for cleaning the divergence of vector potentials is effective.

Figure 2
figure 2

Distributions of the fractional flux error f i over all grid points in A (left) and A p (right).

Full size image

4 Relative Magnetic Helicity Flux in the Simulated Flaring Active Region 8210

We used a version of the photospheric-data-driven 3D MHD code LINMOD3D (Büchner, 2006a, 2006b; Santos et al. 2008; Santos, Büchner, and Otto 2011) to obtain the magnetic and electric field values necessary to calculate the evolution of relative magnetic helicity in a finite volume above a specific active region (AR 8210) and to compare it with the relative magnetic helicity flux through the boundaries obtained according to the method described in Section 2.

For this particular case the simulation code solved the set of MHD equations

$$ \frac{{\partial \rho }}{{\partial t}} = - \nabla \cdot \rho \mathbf{u} ,\\ $$
(19)
$$ \frac{{\partial \rho \mathbf{u}}}{{\partial t}} = - \nabla \cdot \rho (\mathbf{u}\mathbf{u}) - \nabla p + \mathbf{j} \times \mathbf{B} - \nu \rho (\mathbf{u} - \mathbf{u}_\mathrm{n}) ,\\ $$
(20)
$$ \frac{{\partial \mathbf{B}}}{{\partial t}} = - \nabla \times ({\mathbf{u} \times \mathbf{B} - \eta \mathbf{j}} ), $$
(21)
$$ \frac{{\partial p}}{{\partial t}} = - \nabla \cdot p\mathbf{u} - ( {\gamma - 1} )p\nabla \cdot \mathbf{u} + 2\biggl( { {{\gamma - 1}} \over \beta_{0} } \biggr)\eta j^2, $$
(22)

together with Ohm’s law

$$ \mathbf{E} = \mathbf{u} \times \mathbf{B} - \eta \mathbf{j}, $$
(23)

and Ampère’s law

$$ \nabla \times \mathbf{B} = \mu _0 \mathbf{j}. $$
(24)

Here ρ denotes the plasma density, u the plasma velocity, B the magnetic field, p the thermal pressure. The set of equations is closed by assuming the adiabatic condition, i.e., p∝ρ γ with γ=5/3.

To numerically solve the above set of equations a discretization scheme was chosen that suppresses numerical dissipation by using second-order accurate leap-frog finite differences to solve the conservative, non-dissipative hyperbolic partial differential equations on a staggered grid and Dufort–Frankel differences to cope with the parabolic type equations that arise from the dissipative Ohm law (Equation (24)). The variables are made dimensionless by normalizing them to the following quantities: B norm=100 G (gauss), L norm=107 m, \(p_{\mathrm{norm}} = \beta_{0} B_{\mathrm{norm}}^{2} /(2\mu_{0})\) with β 0=0.1, ρ norm=p norm M p/2k B T norm, T norm=106 K, \(v_{\mathrm{norm}} = B_{\mathrm{norm}} /\sqrt{\mu_{0} \rho_{\mathrm{norm}}}\), t norm=L norm/v norm, and the ideal gas law defines the temperature by p=2ρk B T/M p, where k B is the Boltzmann constant and M p is the proton mass. For the normalized plasma beta value (the ratio of plasma over magnetic pressures) we used β 0=0.1 to represent the strong magnetic fields observed in AR 8210. The quantity \({\bf u}_{\mathrm{n}}\) in the momentum equation (Equation (20)) denotes the velocity of the neutral gas to which the plasma is coupled in the chromosphere. In our model the coupling via collisions and charge-exchange processes is quantified by a symbolic friction coefficient ν. In the chromosphere the plasma-neutral gas friction coefficient ν is large compared to the inverse Alfvén transit time that characterizes the propagation of magnetic perturbations. As a consequence, the plasma motion in the chromosphere is strongly coupled to the neutral gas. In contrast, in the upper solar atmosphere, the corona, the friction coefficient ν is much smaller than the inverse Alfvén transit time, i.e., plasma and neutral gas are dynamically decoupled from each other. Since we are interested mainly in calculating the evolution of relative magnetic helicity of the simulation data, radiative losses and heat conduction were neglected as well as the gravity of the Sun.

To start the simulation, first an initial potential magnetic field P is obtained by extrapolating the measured B z component of the magnetic field on the photosphere (z=0), which is compatible with the side boundary conditions of the simulation box according to Otto, Büchner, and Nikutowski (2007). To sufficiently resolve the magnetic field gradients by the 1313 simulation grid, all Fourier modes of the 2D photospheric magnetic field higher than eight are neglected. The resulting filtered photospheric magnetic field distribution is shown in the left panel of Figure 3. The resulting new smoothed magnetic field is relaxed to resume an equilibrium that is to be used as the initial condition for the simulation. Since solar gravity is neglected, the plasma is assumed to be initially in hydrostatic equilibrium (p=const). A horizontal plasma motion is then applied via the neutral gas coupling in the chromosphere. A rotational motion would inject helicity into the system without magnetic flux emergence through the photosphere if the velocity field satisfies ∇⋅u n=0, which inhibits the pile-up of plasma and magnetic flux in the photosphere.

Figure 3
figure 3

(a) The B z component of the magnetic field at z=0 and t=0 s. (b) The horizontal velocity field at z=0 and t=0 s. (c) The horizontal velocity field at z=0 and t=522 s. (d) The horizontal velocity field at z=0 and t=686 s.

Full size image

For definiteness we started the simulation with the line-of-sight component of the photospheric magnetic field in active region NOAA 8210, which is well observed by SOHO/MDI (Scherrer et al. 1995) and analyzed (see, e.g., Régnier and Canfield 2006). The photospheric plasma motion responsible for the evolution of the coronal magnetic field above NOAA AR8210 between 17:13 UT and 21:29 UT on 1 May 1998 was determined by Welsch et al. (2004), Longcope (2004), and Georgoulis and LaBonte (2006). Using the local correlation-tracking technique combined with the solution of the induction equation, Welsch et al. (2004) found for that period of time a clockwise rotation around the major negative magnetic polarity inside AR 8210. This motion produced a shear flow around the polarity inversion line, which was suggested to be responsible for the energy injection into the corona above the active region. Even before this time, a large amount of magnetic flux was injected into the solar atmosphere (see Santos, Büchner, and Otto 2011). Furthermore, the rotational motion itself will also inject magnetic helicity into the corona. In the following we determine the injection, accumulation, and dissipation of magnetic helicity by applying the method described in Section 2 to the results obtained from the numerical simulation of electric and magnetic fields above AR 8210. For this sake a simulation run was started with the photospheric magnetic field corresponding to the one observed at 10:30:04 UT on 1 May 1998, well before the eruption. Figure 3(a) depicts the Fourier-filtered normal magnetic field B z threading the lower boundary. The coronal response to this energy and helicity inflow was simulated in a volume extending 217.5 Mm in all three (X, Y, and Z) directions above the active region. The physical quantities in this volume were calculated on a grid containing 131×131×131 grid points. Figures 3(c) – 3(d) depict the velocity vortex applied to the plasma via its neutral gas interaction in the chromosphere at t=0 s, 522 s, and 686 s. (cf. Santos, Büchner, and Otto 2011, case 3).

While running the simulation, it appeared that for the low initial background magnetic diffusivity (η/ρ 0) of the code (on the order of 1 m2 s−1), which corresponds to a magnetic Reynolds number on the order of Rem≅107, the corona indeed evolves over 522 s like an ideal plasma system. Meanwhile, the Poynting flux through the lower boundary enhances both the magnetic energy content and, at a slower rate, also the kinetic energy generated by accelerated plasma flows in the corona. Only a very small portion of the currents that develop in the course of this almost ideal evolution of the corona is Joule-dissipated into heat (cf. Santos, Büchner, and Otto 2011, Figure 11). The almost ideal evolution continued for 522 s. Then the enhanced currents caused the formation of strong micro-turbulence and microscopic structures, which essentially increases the dissipation, according to Büchner and Elkina (2005, 2006). The anomalous resistivity is switched on when the current carrier velocity V ccv=j/n e reaches a certain threshold (Büchner and Elkina 2005, 2006). After t=522 s the background and anomalous resistivity are re-normalized based on the expected ion skin depth in the corona (Santos, Büchner, and Otto 2011) by an increase of the anomalous current-dependent resistivity over six orders of magnitude. This substantially increases the dissipation and lowers the magnetic Reynolds number locally to the order of unity. Although the dissipation is enhanced first only locally in quasi-separatrix layers below the grid-resolution of the MHD simulation, they will finally cause fast energy conversion via large-scale 3D magnetic reconnection (cf., Büchner 2007). In the actual simulation this transition takes place after 522 s. After this time, the magnetic energy is massively dissipated. This increases the thermal energy in the corona to more than z=1.5, which surpasses even the kinetic energy, to which an increasing portion of the magnetic energy is also transferred by magnetic reconnection (Santos, Büchner, and Otto 2011, Figure 11).

Using the method described in Section 2, we determined the flux of the relative magnetic helicity through the boundaries and the accumulated amount of the relative magnetic helicity during this period. Figure 4(a) shows the temporal evolution of the relative magnetic helicity in the simulation box (pluses) and the accumulated relative helicity (diamonds) based on the helicity change rate of Equation (6). The mean difference in relative magnetic helicity between the two calculations is within 4 %, which further confirms the rationality of our method. It was also found that the temporal evolution of the relative magnetic helicity is consistent with the evolution of the magnetic energy (Santos, Büchner, and Otto 2011, Figure 11). Before t=522 s, the relative magnetic helicity in the volume steadily continues to increase as a result of continuous driving on the boundary. The maximum accumulated relative magnetic helicity is 3.72×1042 Mx2 (Mx: maxwell) when the enhanced dissipation happened. At t=522 s the driving is stopped and the resistivity is re-normalized. As a consequence, helicity decreases in time due to the combined effect of dissipation and outflow through the boundaries. The accumulated relative magnetic helicity in the non-ideal phase is − 8.33×1041 Mx2, which has the opposite sign to the magnetic helicity flux across the boundaries in the ideal case before t=522 s.

Figure 4
figure 4

(a) Temporal evolution of relative magnetic helicity in the simulation box (pluses) and the accumulated relative helicity (diamonds) based on the helicity change rate of Equation (6). (b) Temporal evolution of the normalized relative helicity H m/Φ2. (c) Helicity flux dH m/dt across the boundaries based on Equation (5). (d) Helicity change rate dH m/dt including the helicity flux across the boundaries and dissipation in the simulation box based on Equation (6).

Full size image

Figure 4(b) shows the temporal evolution of normalized relative helicity H m/Φ2, which could indicate how much the magnetic configuration is stressed (Démoulin and Pariat 2009). Φ is the average of the absolute value of the magnetic flux at the bottom boundary in both magnetic polarities. The evolution of the normalized magnetic helicity is similar to the evolution of the relative magnetic helicity, and its maximum value is 0.0298 at t=522 s at the time of the enhanced dissipation of the relative magnetic helicity. Figure 4(c) shows the temporal evolution of the helicity flux dH m/dt across the six boundaries based on Equation (5). We found that at the beginning of the ideal evolution, the change rate of relative helicity is relative low. It starts to increase later and reaches an almost stationary level due to the continued driving through the photospheric motion after about 100 s. The helicity change rate slightly decreases due to helicity back-flow to the chromosphere below the transition region. After t=522 s, the helicity flux across the boundaries decreases significantly because the driving has been stopped; only 2.76×1041 Mx2 of the relative magnetic helicity across the boundaries. Note that the magnetic helicity flux is still positive although negative magnetic helicity is accumulated after t=522 s.

Figure 4(d) shows the temporal evolution of the helicity change rate dH m/dt including the helicity flux across the boundaries and dissipation in the simulation box based on Equation (6). We found that the difference in the profiles between Figures 4(c) and 4(d) is quite small before t=522 s (less than 3 %), which reflects the low initial background magnetic diffusivity and negligible dissipation of magnetic helicity. After t=522 s, significant negative magnetic helicity is accumulated in the simulation volume, i.e., positive magnetic helicity is lost from the simulation box. The first term of Equation (6) mainly contributed to the negative magnetic helicity in the simulation box. The dissipation of relative magnetic helicity, represented by the first term of Equation (6), is − 1.11×1042 Mx2.

5 Summary and Discussion

In this paper we introduced a new method for determining the vector potentials in a given volume (a simulation box) in order to be able to calculate the relative magnetic helicity in the volume for applications, e.g., to the solar corona. The method requires the magnetic flux to be balanced separately on all the sides of the considered volume. We verified the reliability of our approach by applying it to a theoretical nonlinear force-free field to check whether the computed vector potentials indeed reproduce the corresponding magnetic fields and satisfy the Coulomb gauge. Using our method, we also analyzed the evolution of the relative magnetic helicity in a photospheric-data-driven MHD simulation of the active region NOAA 8210.

In the course of an ideal plasma evolution in the corona before and after an eruption, our method appropriately described the evolution of the relative magnetic helicity in both cases. The temporal evolution of the relative magnetic helicity was consistent with the evolution of the magnetic energy. When comparing the accumulated amount of relative helicity inside the volume with the time-integrated flow of relative helicity through the lower boundary and the integrated rate of dissipation inside, we obtained a good agreement (the difference is less than 4 %). The maximum value of the accumulated relative magnetic helicity is 3.72×1042 Mx2 just before the enhanced dissipation. In particular, we found that the loss of relative helicity from the coronal volume at the time of enhanced dissipation (i.e., R m locally on the order of one) cannot be neglected, as had often been assumed. The dissipation of relative magnetic helicity was −1.11×1042 Mx2 in the phase of fast energy conversion via large-scale 3D magnetic reconnection.

The efficiency of our method is such that to simulate the 131×131×131 electromagnetic fields used in this paper, it takes only about 5 min on an outdated Pentium IV 3.0 GHz PC processor to obtain the relative magnetic helicity for one instant of time.

The maximum normalized helicity we found is 0.0298, compared to the similar value, 0.022, of Démoulin and Pariat (2009, Figure 5). Equation (10) of Yang, Zhang, and Büchner (2009) describes the relation between the accumulated helicity ΔH and magnetic flux Φm for active regions:

$$ \log\frac{\Delta H}{{H_0 }}=a\log \frac{{\Phi _m }}{{\Phi_0}}+b, $$
(25)

where a=1.85, b=− 0.41, H 0=1041 Mx2, and Φ0=1021 Mx. Note that the corresponding accumulated magnetic helicity is 3.38×1042 Mx2 if one substitutes the maximum magnetic flux, 1.12×1022 Mx, obtained in this study into Equation (25). This quantity differs by less than 9 % from the actual maximum accumulated magnetic helicity of 3.72×1042 Mx2 at the time of the magnetic reconnection. It is necessary to point out that the time of the reconnection is not controlled by the artificial resistivity; it is controlled by the current carrier velocity reaching a certain threshold that is physically determined by the plasma, at which strong micro-turbulence is excited and microscopic structures form that essentially enhance the dissipation. This implies that the normalized magnetic helicity (\(H_{\mathrm{m}}/\Phi_{\mathrm{m}}^{2}\)) may predict the occurrence of solar events in the solar corona. However, this still needs to be confirmed from observations and simulations in the future.