MODELING THE DISPERSAL OF AN ACTIVE REGION: QUANTIFYING ENERGY INPUT INTO THE CORONA

Duncan H. Mackay; L. M. Green; Aad van Ballegooijen

doi:10.1088/0004-637X/729/2/97

1. INTRODUCTION

The solar corona is a complex environment, where much of its complexity is due to magnetic fields. These magnetic fields are produced through a dynamo action near the base of the convection zone (Charbonneau 2005). Once formed they may become buoyantly unstable, rise through the convective zone, and break through the photosphere (Archontis et al. 2004; Magara 2004; Galsgaard et al. 2007; Murray & Hood 2008; Fan 2009), where they are observed as either sunspots in white light or active regions in magnetograms. Such active regions structure the Sun's atmosphere and provide energy for eruptive phenomena such as solar flares (Benz 2008) and coronal mass ejections (Cremades et al. 2006).

Once active regions form, small-scale motions such as granular and supergranular flows result in the decay of the active region and the dispersal of the magnetic flux across the solar surface in a random walk (Leighton 1964). This random walk leads to the convergence and cancellation of the magnetic flux along polarity inversion lines (PILs), in addition to the spreading of magnetic fields. Such evolution of magnetic elements in the photosphere acts as a driver for the buildup of free magnetic energy in the solar corona. The subsequent evolution of the coronal field through quasi-static equilibria, in response to these motions, may be mathematically modeled through force-free magnetic fields, magnetic fields which satisfy j × B = 0 where j = αB (Priest 1982).

In recent years, nonlinear force-free field (NLFFF) modeling has received much attention (Schrijver et al. 2006; Metcalf et al. 2008). An NLFFF is a special class of force-free fields, where the scalar function α(r) is a function of position, but must be constant along any field line. NLFFFs are of particular interest as they may contain free magnetic energy (Woltjer 1958). The construction of NLFFFs may be broadly split into two groups: static and time-dependent models. Examples of static modeling are the "extrapolation" of photospheric vector magnetic fields into the corona (Schrijver et al. 2006; Regnier & Priest 2007; De Rosa et al. 2009; Wheatland & Régnier 2009; Jing et al. 2010) and direct fitting of NLFFF models to observed filaments and coronal structures (van Ballegooijen 2004; Bobra et al. 2008; Su et al. 2009; Savcheva & van Ballegooijen 2009). Static modeling involves constructing either individual NLFFFs or independent sequences of magnetic configurations, where there is no direct correlation or evolution between the different fields.

In contrast, time-dependent quasi-static modeling evolves the coronal magnetic field through continuous sequences of related NLFFFs based on the evolution of a continuous time-dependent lower boundary condition. This boundary condition may be specified through either idealized magnetic field configurations (van Ballegooijen et al. 2000; Mackay & Gaizauskas 2003; Mackay & van Ballegooijen 2005, 2006a, 2009) or from observations (Mackay et al. 2000; Yeates et al. 2007, 2008). Only through using dynamic models can the buildup of free magnetic energy and magnetic helicity in the corona be studied, along with the evolution of isolated flux and current systems. Such a technique has been successfully applied in the past to model the evolution of the global corona (Mackay & van Ballegooijen 2006a, 2006b; Yeates & Mackay 2009), determine the origin and evolution of solar filaments (Yeates et al. 2008; Mackay & van Ballegooijen 2009), study the formation and lift off of magnetic flux ropes (Mackay & van Ballegooijen 2006b; Yeates & Mackay 2009; Yeates et al. 2010a), and finally the origin of Sun's open magnetic flux (Mackay & van Ballegooijen 2006b; Yeates et al. 2010b).

In this paper, we apply time-dependent quasi-static nonlinear force-free modeling to model the effect that the dispersal of an active region has on the overlying coronal magnetic field. A key difference between this paper and previous studies is the manner in which the time-dependent lower boundary condition is applied. In previous studies, either idealized magnetic field distributions or simplified surface configurations derived from observations were prescribed. In the present study, a new technique of prescribing the time-dependent lower boundary condition is presented. The technique allows the use of observed line-of-sight (LOS) component magnetograms directly as lower boundary conditions, reproducing as observed, the dispersal of the photospheric flux of the active region. From this, the effect of the surface motions on the coronal magnetic field and the subsequent energy input into the corona is determined.

The active region chosen to illustrate this new technique is NOAA AR 8005. It is chosen as it was an isolated region and during the period of observations there appears to be no significant flux emergence and the flux is well balanced throughout. The structure of the paper is as follows. In Section 2, the main properties of the decaying active region are discussed. In Section 3, the modeling technique for both the photospheric and coronal fields is described. The results of the simulation are given in Section 4, while in Section 5 some simple calculations are compared to the simulations to verify the results. Finally, a discussion and conclusions are given in Section 6.

2. OBSERVATIONS

The decaying active region considered in this study is NOAA AR 8005. It emerged on the far side of the Sun and as it rotated onto the limb only dispersed magnetic polarities without sunspots were present. In Figure 1, a full disk SOHO:MDI (Scherrer et al. 1995) 96 minute line-of-sight component magnetogram from 1996 December 18 can be seen. The active region of interest lies at central meridian just below the equator. It has a simple magnetic morphology with a single positive and negative polarity.

**Figure 1.** *SOHO* MDI 96 minute line-of-sight component magnetogram for 1996 December 18 where white areas represent positive magnetic polarities and black negative magnetic polarities. The active region of interest, NOAA 8055, lies in the southern hemisphere at central meridian and is enclosed by the black box.
Download figure:
Standard image High-resolution image

To study the evolution of the active region, SOHO:MDI 96 minute magnetograms⁴ are obtained from 19.12.05 UT on 1996 December 16 to 20.47.05 UT on 1996 December 20, spanning a period of four days. The four-day period was chosen to include two days before and two days after central meridian passage, which occurs late on December 18. During this period as LOS projection effects are minimized, the measurements of the magnetic flux are known to be the most reliable. In total 61 96 minute magnetograms cover the period of interest. The data are corrected for the area foreshortening that occurs away from central meridian using the IDL Solar Software routine drot_map. An area of 181 × 126 pixels is cut out of the rotated magnetograms centered on NOAA 8055 where each pixel is 1 farcs 977. A time sequence of these can be seen in Figure 2. The area was chosen to be large enough to encompass the whole active region, but small enough that approximate flux balance is achieved. The active regions that lie to the north and south of NOAA 8005 have no contribution to the flux in the area considered.

**Figure 2.** Time sequence of derotated *SOHO* MDI 96 minute line-of-sight component magnetograms of NOAA 8005. The time sequence covers the evolution of the active region two days before and two days after central meridian passage (1996 December 18) where white areas represent positive magnetic polarities and black negative magnetic polarities. The images show the active region after it has been corrected for flux balance.(Animations (a and b) of this figure are available in the online journal.)
Download figure:
Standard image High-resolution image

On studying the magnetograms over the period of interest, the evolution of NOAA 8005 is dominated by the dispersal of the magnetic flux from the main flux concentrations. There is no significant emergence of new magnetic bipoles. On studying Figure 2 or the movie of the evolution (see animation a, associated with Figure 2, in the online version of the journal), the motion of the magnetic elements is clearly dominated by a random walk due to supergranular flows and no strong shear or vortical flows are present. In Figure 3, the separation distance (S(t)) between the centers of flux in both the positive and negative polarities can be seen, given as a function of time (t). This separation is defined as S = ∣S(t)∣, where S(t) is the vector pointing between the centers of flux in each polarity,

$\begin{equation} {\bf {S}}(t) = \frac{ \sum _{B_z > 0} B_z(i,j) {{\bf R}}_{i,j}}{\sum _{B_z > 0} B_z(i,j)} - \frac{ \sum _{B_z < 0} B_z(i,j) {{\bf R}}_{i,j}}{\sum _{B_z < 0} B_z(i,j)}, \end{equation} \tag{ 1 }$

where B_z(i, j) is the LOS component of the field at the ith, jth pixel and R_i,j is the position vector of the ith, jth pixel from the origin. The origin is defined to be the lower left corner of the magnetogram. Throughout the four-day period, as the active region crosses central meridian, the separation of the polarities is seen to slightly increase from 68, 000 km to 80, 000 km. This increase is only 18% that of the initial separation. While in general, the region is diverging, there are several instances where patches of strong magnetic flux separate from the main polarities and cancel along the PIL. A key point to note is that during the evolution no systematic shear motions or vortical motions of the magnetic flux are seen. The motions of magnetic elements are mainly due to the buffeting as a result of supergranular motions. Although supergranular motions provide the only obvious flow pattern, due to the large latitudinal extent of the active region differential rotation may have a small non-negligible effect. The evidence and consequences of this will be discussed in Section 6.

**Figure 3.** Separation between the center of flux of the positive and negative polarities as a function of time in days from 1996 December 16 19:12 UT.
Download figure:
Standard image High-resolution image

In Figure 4(a), the variation of the total unsigned magnetic flux within the selected active region can be seen as a function of time. The flux values are given in units of 10²² Mx and time in days from 1996 December 16 19.12.05 UT. The solid line in Figure 4(a) shows the total unsigned flux as derived from the individual magnetograms. The flux values are consistent with a medium-sized active region. Initially over the first day the flux values are seen to increase slightly, possibly due to LOS changes as the region rotates toward central meridian and the emergence of a small bipole on the southwestern edge of the active region. However, between days 1–4 the flux decays away as cancellation becomes significant along the internal PIL (Green & Kliem 2009). By the end of day 4 the flux has dropped by 20% of its original value. This behavior fits well with the observation that it is a decaying active region. The peak flux densities found within both the positive and negative polarities are around 1 kG. In Figure 4(b), the flux difference (flux imbalance) calculated for each magnetogram can be seen as a function of time. Over the four-day period the imbalance is systematically negative, possibly due to dominantly negative surrounding small-scale fields. At worst it is 10% that of the total flux of the active region. As this value is small, the active region may be regarded to be in flux balance to a high degree of accuracy.

To use the magnetograms as a lower boundary condition in the numerical simulations complete flux balance is required. Therefore in each magnetogram the imbalance per pixel is determined and subtracted from each pixel. This correction varies from one magnetogram to the next. In general, it is less than ∣6∣ G and so is comparable to the noise level within the magnetograms and significantly less than the peak flux density. After correction, the variation of the total unsigned flux is given by the dotted line in Figure 4(a). It can be seen that this line follows the uncorrected flux very closely and applying the correction does not significantly affect the net flux of the active region.

3. THE MODEL

3.1. Coronal Model

To simulate the evolution of the coronal magnetic field above the active region, a magnetofrictional relaxation technique is applied. This technique evolves the coronal magnetic field through a sequence of related nonlinear force-free equilibria. The technique is similar to that described in the papers of van Ballegooijen et al. (2000) and Mackay & van Ballegooijen (2006a) where the code described in these papers is adapted to a local Cartesian frame of reference. Previously the technique has been successfully applied to consider the coronal field based on evolving photospheric boundary conditions in either idealized setups (Mackay & Gaizauskas 2003; Mackay & van Ballegooijen 2005, 2009) or approximations to observed magnetograms (Mackay et al. 2000; Yeates et al. 2007, 2008). A key element of these simulations is that through considering a time sequence of evolved equilibria, this allows the memory of previous flux connectivities and currents to be maintained from one time to the next as energy and magnetic helicity is injected into the corona. Such a feature is significantly different compared to independent extrapolations which do not maintain such memory.

The magnetic field, B = ∇ × A, is evolved by the induction equation,

$\begin{eqnarray} &&\frac{ \partial {\bf {A}}}{ \partial t} = {\bf {v}} \times {\bf {B}}, \end{eqnarray} \tag{ 2 }$

where v(r, t) is the plasma velocity and for the present study non-ideal terms are not included. To ensure that the coronal field evolves through a series of force-free states, a magnetofrictional method is employed (Yang et al. 1986). We assume that the plasma velocity within the interior of the box (representing the solar corona) is given by

$\begin{eqnarray} &&{\bf { v}} = \frac{1}{\nu } \frac{ {\bf {j}} \times {\bf {B}}}{ B^2}, \end{eqnarray} \tag{ 3 }$

where j = ∇ × B and ν is the coefficient of friction. Equation (3) represents, in an approximate manner, the fact that in the corona the Lorentz force is dominant and simulates the plasma experiencing a frictional force as it moves with respect to the reference frame. This ensures that as the coronal field is perturbed via boundary motions it evolves through a series of quasi-static NLFFF distributions, satisfying j × B = 0. To carry out the computations, a staggered grid is used to obtain second-order accuracy for all derivatives and closed boundary conditions are applied on the side and top boundaries to match those given in Section 3.2.

A key new feature within this study is that the evolution of the coronal field is driven directly through the insertion of a time sequence of observed LOS component magnetograms. This means that the evolution of the photospheric field occurs in exactly the same way as is seen in the observations. No simplifications, idealization, or smoothing is applied to the observed data and the data are inserted as observed (with only a small correction for flux balance applied). The technique of specifying the sequences of the lower boundary conditions is discussed next.

3.2. Lower Boundary Condition and Initial Condition

To model the evolution of the active region, 61 magnetograms taken at discrete intervals of 96 minutes are available, covering the four-day period around central meridian passage. From these magnetograms, a continuous time sequence of lower boundary conditions is produced. This sequence of lower boundary conditions is designed to match each observed magnetogram, pixel by pixel, every 96 minutes. However, to model the evolution of the coronal field, A_xb and A_yb the horizontal components of the vector potential A on the base that correspond to this magnetogram must be determined. To determine these and produce a continuous time sequence, the following process is applied.

1.
Each of the observed magnetograms, B_z(x, y, k) for k = 1 → 61 are taken, where k represents the discrete 96 minute time index.
2.
Next, the horizontal components of the vector potential at the base, z = 0, are written in the form
$\begin{eqnarray*} &&A_{xb} (x,y,k) = \frac{\partial \Phi }{\partial y}, \nonumber \\ &&A_{yb} (x,y,k) = - \frac{\partial \Phi }{\partial x} \nonumber, \end{eqnarray*}$
where Φ is a scalar potential.
3.
For each discrete time index k, the equation
$\begin{eqnarray*} &&B_{z} = \frac{ \partial A_{yb}}{\partial x} - \frac{ \partial A_{xb}}{\partial y} \nonumber \end{eqnarray*}$
then becomes
$\begin{eqnarray} &&\frac{ \partial ^{2} \Phi }{\partial x^{2}} + \frac{ \partial ^{2} \Phi }{\partial y^{2}}= - B_{z}, \end{eqnarray} \tag{ 4 }$
which is solved using a multigrid numerical method. Details of this method can be found in the papers by Finn et al. (1994) and Longbottom (1998) and references therein. For a full description of the boundary conditions applied see Mackay & van Ballegooijen (2009).

Solving for the scalar potential Φ determines the horizontal components of the vector potential on the base (A_xb, A_yb) for each discrete time interval, 96 minutes apart. To produce a continuous time sequence between each of the observed distributions, a linear interpolation of A_xb and A_yb between each time interval k and k + 1 is carried out. To interpolate the fields 500 interpolation steps are used. By linearly interpolating the horizontal components of the vector potential on the base, this effectively evolves the magnetic field from one state to the other. Therefore, techniques such as local correlation tracking are not required to determine the horizontal velocity. Numerically it also means that undesirable effects such as numerical overshoot or flux pileup at cancellation sites do not occur and no additional numerical techniques to remove these problems have to be applied.

The technique described above means that there are two timescales involved in the evolution of the lower boundary condition. The first, which is 96 minutes, is the timescale between observations, the second, which is 11.52 s, is the timescale introduced to produce the advection of the magnetic polarities between observed states by interpolation along with the relaxation of the coronal field. The process described above reproduces the observed magnetograms at each 96 minute discrete time interval and therefore produces a highly accurate description of the magnetogram observations.

In the animation (b) associated with Figure 2, a comparison of the actual observations (left) and normal magnetic field component produced by the above technique (right) can be seen. A good agreement between the two is found. The new technique means that neither feature tracking nor local correlation tracking is required to derive horizontal velocity profiles, in order to simulate the evolution of magnetic fields between observed distributions, therefore removing a source of uncertainty. At this point, it should be noted that the above technique only specifies A_xb, A_yb on z = 0. It however does not specify A_z which lies 1/2 a grid point in the box and is determined by Equation (2), the coronal evolution equation. Non-potential effects near the base, as a result of the evolving lower boundary condition, may be contained within this term. These are systematically produced as the horizontal field components on the lower boundary evolve in time.

The choice and construction of the initial condition for the coronal field is now discussed. For the initial condition, a wide variety of choices may be made. These range from potential to NLFFFs. One method to constrain this is to use coronal images of the active region to determine the field geometry. The most realistic form of the initial condition is the NLFFF. However, since no vector magnetic field data are available to us, such a field cannot be constructed. On comparing various linear force-free field solutions with coronal images we find that no single value of the force-free parameter α fits the whole corona structure. Therefore, as vector magnetic field data are not available, we choose the initial condition to be a potential field, as the solution is both unique and stable. By choosing this initial condition, when analyzing the evolution of the active region we restrict ourselves to studying trends in quantities, rather than absolute values. In the future, with SDO:HMI vector field observations this restriction of only considering trends will be removed. In Section 6, the consequences of changing the initial condition from that of a potential field are discussed.

To construct a potential field, the equation

$\begin{eqnarray} &&\nabla ^{2}{\bf A} = 0 \end{eqnarray} \tag{ 5 }$

is solved in a cube with sides ranging from 0 < x, y, z < 6 on a 256³ grid where 1 unit = 60, 733 km. The cube represents an isolated region of the solar surface, where flux may only enter or leave through the lower boundary (z = 0). As flux may only enter/leave through the lower boundary, all field lines must start and end in the z = 0 plane, and no field lines may leave through the side boundaries. This means that B_n, the normal component of B, vanishes on each of the faces except z = 0 (which represents the photosphere). To satisfy this, it is assumed that the tangential components of A are zero on each of the faces, except z = 0. In addition, the normal derivative of the normal component of A is set to zero on each of these faces. If Equation (5) is solved subject to the above boundary conditions, it is straight forward to show that the solution will have ∇.A = 0 everywhere within the domain (Finn et al. 1994). Thus, the initial condition is a potential field associated with the imposed normal field on the boundary with the choice of the Coulomb gauge (∇.A = 0). The initial potential field is constructed from the first observed magnetogram at 19.12.05 UT on 1996 December 16.

4. SIMULATION RESULTS

4.1. Magnetic Field Lines

In Figure 5, an illustration of the field lines of the potential field used as the initial condition can be seen. From these field lines, it is clear that the decaying active region has a simple bipolar form. The field lines connecting between the positive flux (solid contours) and the negative flux (dashed contours) form semi-circular loops.

In Figure 6(a), the illustration shows the field lines for the NLFFF after four days of evolution. This figure can be compared to Figure 6(b) where the corresponding field lines for a potential field deduced from the same normal field component on the base are plotted. In each case, the starting point for the field lines is taken on the base within the positive flux region. On comparing the two images, it is clear that there are many differences in the connectivity and structure of the field. The main differences occur low down along the PIL between the two main polarities. Here the field lines of the non-potential field have a much more sheared structure as a result of energy and helicity being injected along the field lines by the small-scale convective motions.

Another major difference is that for the field lines lying at the southern end, the connectivity of the field is very different from that of the potential field. This is because the connectivity within the nonlinear force-free simulation is initially defined at the start of the simulation and is preserved throughout the simulation (except where numerical diffusion becomes large). The largest field lines within the simulation are only slightly different as the small-scale convective motions have not been able to inject helicity along the full length of these field lines during the time period of the simulation.

4.2. Magnetic Energy

In Figure 7(a), the graph of total magnetic energy stored within the coronal field can be seen as a function of time. The dotted line is for the NLFFF simulation, while the solid line is for a potential field deduced from the same normal field component on the lower boundary as that of the NLFFF (see Section 3.2). Both values are initially equal to one another as the initial condition is a potential field. For both cases, the total energy is of the order of 10³² erg which is typical for a small active region.

Initially over the first day the energy of both the NLFFF and the potential field increases. This increase is due to the increasing flux values; however, after one day the magnetic energy in both peaks and then starts to decrease as the level of flux decreases due to cancellation. While both the NLFFF and the potential field have a similar behavior, it can be seen that the NLFFF has a systematically higher energy compared to the potential field. This increased energy is due to the small-scale convective motions which inject a Poynting flux into the corona and evolve the initial coronal field away from potential.

In Figure 7(b), a graph of the total free magnetic energy within the coronal field is plotted as a function of time. The total free energy is defined as the difference between that of the non-potential and potential fields when integrated over the entire volume. This quantity is always positive and by the end of the simulation the free magnetic energy is approximately 8 × 10³⁰ erg which is 10% that of the potential field. An interesting feature is that throughout the simulation the free magnetic energy as a result of the convective motions increases steadily. The rate of increase is around 2.5 × 10²⁵ erg s⁻¹, indicating that the small-scale motions deduced from the magnetograms may inject large amounts of energy into the corona.

In Figure 8, the images illustrate the locations of free magnetic energy storage for days 1–4. The plots are in the x–z plane and give the free magnetic energy summed along the LOS,

$\begin{equation} E(x,z)= A \int \frac{ \big({\bf B}^2 - {\bf B_p^{\rm 2}}\big)}{8 \pi } dy, \end{equation} \tag{ 6 }$

where B is the magnetic field of the NLFFF and ${\bf {B_p}}$ is the magnetic field of the potential field satisfying the same normal field components on the boundaries. The factor A represents the area of the column being summed over (A = 2.02 × 10¹⁶ cm²). The area factor is included so that the free energy along the LOS is computed in units of erg. In each plot, the x-direction represents the full length of the computational box, but the z-direction is only half the height. White locations denote where the NLFFF has a higher energy than the potential field when integrated along the LOS (i.e., locations where free magnetic energy is stored). Black denotes where the NLFFF has a lower energy when integrated along the LOS (i.e., locations where there is no free energy). We note that while the NLFFF must and does have a volume-integrated energy that is greater than that of the potential field (see Figure 7), there is no restriction that in any subvolume this must always be true. Hence the negative values of the line integral quantity in Equation (6) are physically valid.

As the small-scale convective motions advect the field between days 1 and 4, the locations of free magnetic energy storage expand up into the coronal volume and mainly lie in the center of the box. From these images, it is clear that the main locations of free magnetic energy are low down in the corona (below 30 Mm) and between the two main flux concentrations. It is built up at this location, as here the magnetic fields are the strongest and the energy injection through Poynting flux is the greatest. In addition, these field lines are the shortest and occupy a smaller volume.

The images shown in Figure 8 are useful for illustrating the general spatial locations of free magnetic energy. Often such scaled images lead to misleading interpretation on the exact values. To quantify this, in Figure 9, graphs of the variation of free magnetic energy can be seen for day 4, along the paths given by the dashed lines in Figure 8. In Figure 9(a), the free magnetic energy can be seen as a function of x at a height of 7300 km (z = 0.12 units). Positive values denote locations where the nonlinear force-free field has a higher energy, negative values where it is lower (when integrated along the LOS). From this graph it is clear that in the coronal volume the free energy is mainly positive. Near the center of the domain there is a small region where the NLFFF has a lower energy than the potential field. However, the negative value of the LOS-integrated free energy is small compared to the locations where it is positive. From this it can be seen that the scaling and color shown in Figure 8 exaggerates the level where the free magnetic energy is negative. In Figure 9(b), the variation of free magnetic energy can be seen as a function of z, taken at x = 3 units. It is clear that along this cut the free energy is always positive. In general, the amount of free energy available decreases with height in the box.

**Figure 9.** Graphs showing the free magnetic energy in units of 10²⁷ erg as a function of (a) x at z = 0.12 units and (b) z at x = 3 units along the dashed lines shown in Figure 8. Positive values denote locations where the NLFFF energy has a higher energy compared to that of a potential field.
Download figure:
Standard image High-resolution image

Alternative views of the free magnetic energy storage can be seen in Figure 10. In Figure 10(a), it is shown in the y–z plane summed along the x-direction (Equation (7)) and in Figure 10(b) for the x–y plane summed along the z-direction (Equation (8)). In each case results are shown for day 4 and the same scaling is used as in Figure 8:

$\begin{eqnarray} &&E(y,z) = A \int \frac{ \big({\bf B}^2 - {\bf B_p^{\rm 2}}\big)}{8 \pi } dx, \end{eqnarray} \tag{ 7 }$

$\begin{eqnarray} &&E(x,y) = A \int \frac{ \big({\bf B}^2 - {\bf B_p^{\rm 2}}\big)}{8 \pi } dz. \end{eqnarray} \tag{ 8 }$

**Figure 10.** Images showing the locations of free magnetic energy storage on day 4 in (a) the y–z plane summed along the x-direction and (b) the x–y plane summed along the z-direction for the NLFFF. White areas denote locations where there is excess free magnetic energy compared to that of a potential field. The same scaling is used as in Figure 8.
Download figure:
Standard image High-resolution image

Once again it is clear that the locations of free magnetic energy storage exist in the center of the computational box where the main flux concentrations are at their strongest. In Figure 10(b), the main spatial location where the LOS integrated free energy is negative is surrounded by positive free-energy locations. To determine why such a region exists, Figure 11 compares the photospheric magnetic field distribution of the active region to the location where the LOS free magnetic energy is negative. In this plot, thin solid lines denote the positive flux, dashed lines the negative flux, and the thick solid line where the LOS-integrated free energy is −6 × 10²⁶ erg or less. It can be clearly seen that the negative region is spatially correlated with the cancellation site between the positive and negative polarities. Therefore, it is due to the removal of photospheric and coronal flux at this location as a result of flux cancellation (Figure 4).

**Figure 11.** Comparison of the underlying magnetic field distribution with locations where the NLFFF has less energy than the potential field on day 4. In the plot, the positive and the negative photospheric magnetic flux are given by the thin solid and dashed contours. The thick solid contour denotes the area where the NLFFF has an energy less than 6 × 10²⁶erg of the corresponding potential field when summed in the z-direction.
Download figure:
Standard image High-resolution image

4.3. Magnetic Helicity

As the coronal magnetic field is subjected to the small-scale random motions, these motions inject not only magnetic energy into the field but also magnetic helicity (Démoulin & Pariat 2009). This is seen by the increasing complexity of the field. To quantify this injection, a calculation of relative magnetic helicity (H_r) is carried out. This quantity is defined as

$\begin{equation} H_r = \int _v {\bf {A.B}} d\tau - \int _v {\bf {A_p.B_p}}d\tau, \end{equation} \tag{ 9 }$

where A and B are the vector potential and magnetic field of the NLFFF. Similarly, $\bf A_p$ and $\bf B_p$ are the vector potential and magnetic field of a potential field which satisfies the same normal field components as that of the NLFFF on all boundaries (namely, B_z on z = 0 and B_n = 0 on all other faces). The definition of relative magnetic helicity was first introduced by Berger & Field (1984) as an invariant and therefore meaningful measure of magnetic helicity. It should be noted that in constructing the potential field for the relative helicity calculation the technique described in Section 3.2 is used. The horizontal components of the vector potential on the base, A_xb and A_yb, are identical for both the nonlinear force-free and potential magnetic fields. This means that the expression for the relative helicity does not include an addition surface integral term.

In Figure 12, the graph shows the variation of the relative magnetic helicity as a function of time (solid line). From the behavior of the curve it can be seen that over the first day there is no definite trend of helicity injection by the surface motions. The relative helicity oscillates between positive and negative values. In contrast, between days 1 and 4 there is a definite trend of positive helicity injection, with the relative helicity showing an increasing positive value. Positive helicity injection is consistent with observations that show that active regions in the southern hemisphere have a dominant positive helicity (Pevtsov et al. 1995). The dash-dotted line in Figure 12 is a linear best fit to the relative helicity curve between days 1 and 4. The gradient of the line gives a helicity injection rate of 1.218 × 10³⁴ Mx² s⁻¹. Such an increase of helicity within the coronal field can only be the result of helicity injection through the lower boundary as a result of the relative motion of the magnetic fragments (Démoulin & Berger 2003; Démoulin & Pariat 2009) as during the period of observations there is no increase in the magnetic flux. However, through considering the evolution of the magnetic elements in the magnetogram, there is no clear indication as to the origin of the helicity injection, as the dominant motion acting on the active region are small-scale random motions.

In contrast, on comparing Figure 12 to the graphs of flux variation (Figure 4) and separation between the centers of flux elements (Figure 3), it can be seen that the positive helicity injection may possibly be related to large-scale properties of the evolution of the active region. To begin with on day 1 the flux peaks, after which it systematically starts to decrease. This decrease is the result of first convergence and then cancellation of the polarities along the internal PIL of the bipole. Second, around the same time the separation of the sources increases, so there is an overall divergence of the active region polarities. Either of these two processes (convergence/divergence) may result in the net positive helicity injection. However, within the present simulation we are unable to determine which is the cause. To compute this, techniques such as local correlation tracking would have to be applied (Chae 2001; Pariat et al. 2005; Jeong & Chae 2007). While there is no clear indication of its origin, a third possibility is that such a positive injection of helicity is consistent with the effect of differential rotation as will be discussed in Section 6.

While the relative helicity has a positive value, the amount of helicity is small. Active regions containing fluxes of 1 × 10²² Mx may well contain helicity on the order of 10⁴²–10⁴⁴ Mx² (Tian & Alexander 2008). Upon studying the coronal field in more detail, it is found that the random motions inject large amounts of both positive and negative helicity, but in near equal amounts over the four days of the simulation. This is what is expected from such motions, therefore resulting in a low net helicity.

In Figure 13, a comparison of (a) the normal magnetic field component on the lower boundary and (b) the distribution of α on the lower boundary can be seen for the final time snapshot of the simulation. Similar plots were found at other times. In Figure 13(a) white/black represents positive/negative flux, while in Figure 13(b) red/blue denotes positive/negative values of α. The magnetic field values are set to saturate at ±100 G while the α values saturate at ±2.73 × 10⁻⁹ m⁻¹. On comparing the distributions it can be seen that the spatial distribution of α is well correlated to that of the magnetic field. The distribution of α takes two forms. First, in locations of strong field concentrations the α distribution takes the form of large coherent patterns. Second, it takes the form of small isolated patches of intermingled signs which lie in small flux concentrations. It is clear from comparing the two images that within the large concentrations of flux a single magnetic element of one polarity may contain both positive and negative values of α. Such a mixture of positive and negative values of α, as a result of small-scale convective motions, supports the results that the coronal field of the active region cannot be modeled by a linear force-free field as found in Section 2. The positive and negative α distributions are similar to the divided distribution of current used in the idealized model of Régnier (2009). An interesting feature of the comparison is that within the large flux concentrations the northern portions appear to contain more positive values of α, while the southern portions have more negative values. In addition to the α values that lie within the center of the simulation, values may be seen along the edge of the magnetogram area. Such edge values are a pure boundary effect and should not be considered further.

For the active region considered here no vector magnetogram data are available for comparing observed α distributions or trends of α in the observations with those obtained within the simulation. In future with the launch of SDO and the SDO:HMI instrument, regular vector magnetic field data will be available which can be compared to future simulations applying this technique. Through this we will be able to determine how much helicity is injected due to small-scale motions, and in addition how much must be injected by other effects (e.g., torsional Alfven waves). This provides a new tool for the quantification of such effects which may be applied directly to observations.

5. SIMPLE CALCULATIONS

To determine if the rate of energy buildup within the present simulation is realistic, a simple order of magnitude calculation is carried out to determine the rate of inflow of energy into the lower boundary via the Poynting vector, S = (c/4π)E × B. The Poynting flux (P) through the lower boundary is

$\begin{equation} P = - \frac{c}{4 \pi } \int _s ({\bf E} \times {\bf B}).{\bf {da}}, \end{equation} \tag{ 10 }$

where E is the electric field and ${\bf {da}}= da {\hat{\bf n}}$ . In order of magnitude terms this equation may be simplified to

$\begin{equation} P \approx \frac{v_h B_z B_h A_m}{4 \pi }, \end{equation} \tag{ 11 }$

where v_h is the typical horizontal velocity, B_z is a typical normal field component, B_h is a typical horizontal field strength, and A_m is the area of energy input. The area of input is taken as the area of the observed magnetogram A_m = 4.7 × 10²⁰ cm²; v_h = 0.5 km s⁻¹ is the peak flow rate of a supergranular cell. The values of the average field components are then taken to be B_z = 12 G and B_h = 2 G, where these values are determined from the initial potential field. Through this an upper estimate of the Poynting flux becomes 4.5 × 10²⁵ erg s⁻¹. This value is around a factor of two greater than the value deduced within the simulations of 2.5 × 10²⁵ erg s⁻¹. Although the order of magnitude calculation is larger by a factor of two, this value is an overestimate as it assumes (1) that the flow is always at the peak rate and (2) flows that are steady and systematic. In reality, all elements in the magnetogram do not move at all times and in the simulations the flow is irregular. Taking this into account it can be seen that the simulation produces a realistic energy input into the coronal volume.

In the second calculation, we consider how important the energy input from the small convective motions is relative to an estimate of the radiative losses for the active region. The radiative loss in an optically thin plasma at coronal temperatures is given by

$\begin{equation} E_R = n_e^2 R(T) \, \rm {erg\, cm^{-3}\, s^{-1}}, \end{equation} \tag{ 12 }$

where n_e is the electron number density, T is the temperature, and R(T) is the radiative loss function which is given by Rosner et al. (1978) as

$\begin{equation} R(T) \approx 10^{-17.73} T^{-2/3}, \quad \rm {6.3 < log \, {\it T} < 7.0}. \end{equation} \tag{ 13 }$

The average temperature of the active region is determined to be log T = 6.41 by using the soft x-ray telescope filter ratio technique on the Al.1 and AlMg filters. We take the volume of the coronal emitting plasma to be 10²⁹ cm³. Note that this is not the volume of the computation box but rather an upper estimate of the volume of emitting plasma in the active region. Taking also an upper estimate of the electron number density to be 10⁹ cm⁻³ provides an upper estimate of the coronal radiative losses as ∼1 × 10²⁵ erg s⁻¹. This value can be seen to be of the same order as the energy buildup due to the small-scale convective motions.

While the two values are consistent with one another, it should be noted that the free magnetic energy in the simulation is held within the coronal field and therefore not immediately accessible to heat the corona and account for the radiative losses. For such a process to take place, an additional energy release mechanism would have to act. Possible mechanisms include nanoflares (Parker 1988; Browning et al. 2008; Wilmot-Smith et al. 2010) or turbulent reconnection (Heyvaerts & Priest 1984). At this point, it is unclear if such a mechanism may release the energy. Also, the simple calculation does not take into account the effect of thermal conduction which is an important factor in the energy balance of hot coronal loops. However, taking these limitations into account, the calculation does show that convective motions may play a key role in the energy buildup of the corona.

Using the rates of energy injection determined in the simulation and assuming that it remains constant over the lifetime of the active region, it is possible to estimate how much free energy the active region could have on December 16, the start day of the present simulations, due to prior evolution as it rotated onto the disk and toward central meridian. Assuming that it emerged just before disk passage the active region would have evolved for five days. As a result, the small-scale convective motions would inject a total of 1.1 × 10³¹erg of free magnetic energy into the coronal field. This value is 10% of the active region energy on December 16. While the energy may have been 10% higher, we cannot use this to constrain the initial condition without having vector magnetic field data to constrain the distribution of α on the lower boundary.

6. DISCUSSION AND CONCLUSIONS

In this paper, a new technique for modeling NLFFFs directly from LOS component magnetogram observations has been presented. The key new feature of this technique is that sequences of LOS magnetograms can be directly used as lower boundary conditions to drive the evolution of coronal magnetic fields between successive force-free equilibria. As a result, the lower boundary condition in the model closely follows the observed magnetic field evolution as given by the observed LOS magnetograms over the period of observations.

The technique is illustrated by applying it to MDI observations of a decaying active region, NOAA AR 8005. The active region was close to central meridian and in flux balance to a high degree. It was modeled over a four-day period, where 61 individual SOHO:MDI 96 minute magnetograms were available to use as lower boundary conditions. During this time, the dispersal of the active region was mainly dominated by random motions due to small-scale convective cells and flux cancellation along the internal PIL.

Through applying the 61 magnetograms as lower boundary conditions and simulating a continuous evolution of the surface and coronal field between them, the buildup of free magnetic energy into the corona was studied. It should be noted here that the initial condition for the coronal field in the simulation was taken to be a potential field. Due to the availability of only LOS component magnetograms this was the only unique choice. However, at the time corresponding to the start of the simulation, the actual coronal field of the active region need not have been potential. Therefore our aim was not to reproduce the evolution of the coronal loops as seen in coronal observations, but rather to study the trends of energy input over the period of four days, where this is deduced from using observed magnetograms.

When the initial coronal field is taken to be potential, the small-scale motions inject around 2.5 × 10²⁵ erg s⁻¹ of free magnetic energy, this energy is mainly stored in the low corona below 30 Mm. However, even if the initial state on the Sun was different, this would not significantly alter the trends of energy injection found within the simulation. To illustrate this, we have repeated the simulation using an initial condition of a linear force-free field, with α values of 1.55 × 10⁻⁸ m⁻¹ and 1.86 × 10⁻⁸ m⁻¹. These values are 80% and 88% of the first resonant value of α. Only positive values of α are chosen as the active region lies in the southern hemisphere. With these alternative initial conditions similar results and trends are found to those shown in Figure 7. The only difference is that the magnetic field has a higher initial energy. The subsequent rates of energy injection are 2.65 × 10²⁵ erg s⁻¹ and 2.8 × 10²⁵ erg s⁻¹, respectively. These values are consistent with those found when using the potential field. The slightly higher values occur due to the increased horizontal component of field in the linear force-free field leading to a higher Poynting flux. For all cases, the energy injection through this new technique is consistent with estimates derived from Poynting flux calculations. After just four days of evolution, the free energy is 10% that of the potential field. Such an energy buildup is sufficient to explain the radiative losses at coronal temperatures within the active region (1 × 10²⁵ erg s⁻¹). Therefore small-scale convective motions play an integral part in the energy balance of the corona.

Through studying the evolution of the relative magnetic helicity in the coronal field, it was found that near equal amounts of positive and negative helicity are injected into the coronal field, but with a slight preference for positive helicity injection in the later stages. This is consistent with the active region being buffeted by small-scale random motions and observations that show that active regions in the southern hemisphere have a preference for positive helicity (Pevtsov et al. 1995). One possible cause of this slight positive helicity injection is the effect of differential rotation. In the southern/northern hemisphere differential rotation will inject positive/negative helicity into the corona (DeVore 2000). Démoulin et al. (2002) showed that helicity injection by differential rotation is composed of two terms, which for a decayed and dispersed active region are of opposite signs. The injection rate also increases with active region latitude and varies with active region tilt angle. During the short time period of our study, where the active region lies near the equator, we expect differential rotation to play only a minor role in ejecting positive helicity into the active region.

The new technique presented here has wide ranging applications for modeling the evolution of photospheric magnetic fields observed on the Sun and the subsequent effect these have on the coronal field over long periods of time. In the future, we shall use new high-resolution, high-cadence full disk vector magnetic field data from SDO:HMI in combination with coronal observations to constrain the initial condition. With this we may then compare the evolution of the coronal magnetic field with that seen in coronal observations. Under such constraints it will become a powerful tool for studying the nature and properties of coronal magnetic fields.

D.H.M. thanks UK STFC for their financial support and the Royal Society for funding an 81 processor supercomputer under their Research Grants Scheme. L.M.G. thanks the Royal Society and Leverhulme Trust for financial support. We acknowledge the use of data provided by the SOHO/MDI instrument. The authors thank the referee for constructive comments which have improved this paper.

MODELING THE DISPERSAL OF AN ACTIVE REGION: QUANTIFYING ENERGY INPUT INTO THE CORONA

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. OBSERVATIONS