A&A 468, 313-321 (2007)
DOI: 10.1051/0004-6361:20054495

Force-free magnetic field extrapolation for MHD boundary conditions in simulations of the solar atmosphere

A. Otto¹ - J. Büchner² - B. Nikutowski²

1 - Geophysical Institute, University of Alaska, Fairbanks, AK 99775, USA
2 - Max-Planck-Institut für Sonnensystemforschung, Max-Planck-Str. 2, 37191, Katlenburg-Lindau, Germany

Received 8 November 2005 / Accepted 14 February 2007

Abstract
Context. In recent years the accuracy of magnetic field observations in the solar atmosphere has made considerable progress. Similar progress is being made in computer modeling of complex plasma systems and computer capabilities.
Aims. To consider observed solar magnetic field structure in numerical simulations a new extrapolation method for solar magnetic fields is used to incorporate such fields into three-dimensional MHD simulations
Methods. The simulation employs a new modified linear magnetic field extrapolation which is specifically designed to satisfy symmetry conditions which are generic to MHD models. The model domain includes photosphere, chromosphere, and corona.
Results. The new model is applied to several solar field configurations and results are compared to three-dimensional field structure from observations and another extrapolation method. The new model provides a simple and efficient method for the simulation of observed solar magnetic field structures by constructing a three-dimensional initial field that is consistent with symmetry boundary conditions of MHD simulations.

Key words: magnetohydrodynamics (MHD) - Sun: magnetic fields - Sun: corona - methods: numerical - plasmas

1 Introduction

Modern investigations of the solar atmosphere provide a wealth of new information on the magnetic and plasma structure. A valuable tool to interpret these data is numerical simulation of the magnetic coupling between the different regions of the solar atmosphere starting from the photosphere, where the magnetic field can be determined with the highest spatial resolution, through the chromosphere and the corona.

Existing solar simulation models usually use model magnetic fields on meso-scales (e.g., Choe & Cheng 2002; Cheng et al. 2003; Galsgaard et al. 2003; Antiochos et al. 2002; Priest & Forbes 2000) or address global solar phenomena (e.g., Linker et al. 2001; Mikic et al. 1999; Linker et al. 2003) to investigate the influence of the photospheric boundary conditions such as magnetic shear, plasma motion, emerging flux, etc. A first meso-scale simulation model based on observed line of sight magnetic fields has been used by Gudiksen & Nordlund (2002) to study coronal heating. The further development of such models is highly desirable to interpret the properties and evolution of structures in the solar atmosphere as observed in various wavelengths of the solar radiation. Here we describe a new method to incorporate observed photospheric magnetic fields in magnetohydrodynamic simulations to study the evolution of solar magnetic fields. For applications we refer to other publications (Büchner et al. 2005b,2004,2005a)

The method consists of two parts. First, the observed magnetic field along the direction of the line of sight has to be converted into three-dimensional magnetic fields. Our approach employs a solution of the force free magnetic field condition

$$\nabla\times{\vec B}=\kappa{\vec B}$$ (1)

with a constant coefficient $\kappa$ (we use $\kappa$ instead of the common $\alpha$ the latter is used as an angular variable in the Fourier expansion in this paper). Although there are various linear and nonlinear methods for the extrapolation of solar magnetic field most of these are not suitable to be used as initial conditions for MHD simulations because of the required MHD boundary conditions.

The second step consists of the use of the extrapolated magnetic fields as an initial condition in the MHD simulation. The simulation code uses a finite difference approximation for the full set of MHD equations where the plasma is coupled to a neutral fluid. The latter is important in the photosphere and chromosphere.

The following section focuses on the extrapolation of the magnetic field consistent with MHD boundary conditions. Properties of the initial state and the extrapolated magnetic field for a case of a coronal bright point analyzed by Brown et al. (2001) are discussed in Sect. 3. Section 4 addresses the simulation method and the last section presents a summary and discussion of the presented results.

  
2 Magnetic field extrapolation

Various methods have been suggested to derive three-dimensional magnetic field structures from the line-of-sight magnetic field component. Several of these methods use potential (Semel 1968; Schmidt 1964), linear (Semel 1988; Seehafer 1978; Chiu & Hilton 1977), or nonlinear force free fields (Cuperman et al. 1990; Amari et al. 1997; Aly 1989) with a non-constant value of the force-free parameter $\kappa$. The use of the extrapolated fields for MHD simulations, however, requires that the boundary conditions of the extrapolation are consistent with boundary conditions that can be applied to the MHD simulation. For the purpose of constructing a method which can easily be employed and automated to construct three-dimensional magnetic equilibria a constant value $\kappa$ is the most straightforward choice. The particular numerical value of $\kappa$ can always be chosen to represent characteristics of the observed magnetic field such as the stretching with altitude over the photosphere or a particular (inferred) magnetic helicity.

Our approach uses a method similar to the one outlined by Seehafer (1978). This solution is constructed through a Fourier series expansion of the magnetic field components where the coefficients of the horizontal components (B_x and B_y) are related to the vertical component z component B_z through the force free condition (1) and $\nabla\cdot{\vec B}=0$. This procedure implies periodic solutions in the x and y directions. The solution of Seehafer (1978) has the particular advantage of using a symmetry condition for B_zat the x and y boundaries of the system which allow a solution in a quarter of the full periodic system, i.e., the solution is constructed in the domain $0\leq x\leq L_{x}$, $0\leq y\leq L_{y}$ while the full period covers the domain $-L_{x}\leq x\leq L_{x}$, $-L_{y}\leq y\leq L_{y}$. This avoids the need for periodicity in the magnetic field.

However, this solution is not suitable as an initial condition for an MHD simulation which requires well defined symmetry conditions at the boundaries. For B_z the solution by Seehafer (1978) employs

$$\begin{eqnarray*}B_{z}(-x,y,z)=-B_{z}(x,y,z), & & {\rm for}~~{\rm the}~~{\rm bou... ...(x,y,z) & & {\rm for}~~{\rm the}~~{\rm boundary}~~{\rm at}~~ y=0 \end{eqnarray*}$$

in any plane $z={\rm const.}$ such that B_z(0,y,z)=B_z(x,0,z)=0.

Using proper MHD symmetry boundary conditions (a detailed discussion will be given in the following section) at x=0 implies either B_x(-x,y,z)=-B_x(x,y,z) or B_x(-x,y,z)=B_x(x,y,z). Allowing for a non vanishing normal component at this boundary therefore implies the second choice. Furthermore $\nabla\cdot{\vec B}=0$requires $\partial B_{y}/\partial y=0$ at the x=0 boundary ( $\widetilde{x}=-x$, and $\partial/\partial\widetilde{x}=-\partial/\partial x$) which can only be satisfied in general if B_y(-x,y,z)=-B_y(x,y,z)or B_y(0,y,z)=0. Thus the normal component B_x is symmetric and nonzero while the tangential components are zero and reverse sign. However, using these symmetries and calculating the current at the boundary x=0 yields $j_{x}(-x,y,z)=\partial_{y}B_{z}(-x,y,z)-\partial_{z}B_{y}(-x,y,z)=-j_{x}(x,y,z)$. Thus the normal component of the current must reverse sign at the boundary.

The other boundaries at x and y have the same symmetry conditions. While there is no problem in the interior domain the conditions at the x and y boundaries are in contradiction to the force free condition (1) because the magnetic field normal to the boundary is continuous and the current normal to the boundary reverses sign. The reason for this mismatch is that the symmetry boundary condition is not the proper continuation of the magnetic field into the region x<0. The proper continuation into this region is determined by the fully periodic solution from -L_x to +L_x. For this reason the solution by Seehafer (1978) is not suitable as initial condition for an MHD simulation because an MHD simulation requires boundary conditions which are local symmetry conditions or requires to consider the full period of the system as considered by Gudiksen & Nordlund (2002). However, as will be illustrated the approach by Seehafer (1978) can be modified to accommodate well defined MHD symmetry conditions.

2.1 MHD symmetries for system boundaries

The MHD simulation will use the full set of MHD equations

    

$$\frac{\partial\rho}{\partial t} -\nabla\cdot\rho{\vec u}$$ = (2)

$$\frac{\partial\rho{\vec u}}{\partial t} -\nabla\cdot\rho{\vec u}{\vec u}-\nabla p+{\vec j}\times{\vec B}$$ = (3)

$$\frac{\partial{\vec B}}{\partial t} \nabla\times({\vec u}\times{\vec B}-\eta{\vec j})$$ = (4)

$$\frac{\partial p}{\partial t} -\nabla\cdot p{\vec u}-(\gamma-1)p\nabla\cdot{\vec u}+(\gamma-1)\eta{\vec j}^{2}$$ = (5)

$${\rm with} {\vec E}=-{\vec u}\times{\vec B}+\eta{\vec j}$$

$$\nabla\times{\vec B}=\mu_{0}{\vec j}.$$

The system boundaries are considered to at $x=\left[0,L_{x}\right]$and y=[-L_y/2,L_y/2]. We will now consider symmetry conditions at the boundary x=0 for a plane $z={\rm const.}$and distinguish two basic cases:

i)

Local mirroring at the plane x=0: $\widetilde{x}=-x$, $\widetilde{y}=y$, and $\widetilde{z}=z$, which implies $\partial_{\widetilde{x}}=-\partial_{x}$, $\partial_{\widetilde{y}}=\partial_{y}$, and $\partial_{\widetilde{z}}=-\partial_{z}$where $\partial_{x}\equiv\partial/\partial_{x}$.

The MHD variables must transform such that the MHD equations remain invariant for the transformation implied by the boundary condition. There are two possibilities for the individual variables. They can transform either as f(-x,y,z)=f(x,y,z) which we call symmetric (abbreviated "s'') or they may transform as g(-x,y,z)=-g(x,y,z)which we call antisymmetric (or "a''). Positive definite variables like $\rho$, p, and $\eta$ must transform symmetric for instance $\rho(-x,y,z)=\rho(x,y,z)$. Examining the full set of MHD equation with these rules yields two possible sets of transformations which maintain the invariance of the equations. The result for the individual variables is summarized in Table 1.

Table 1: Transformation properties for the transformation  $\widetilde{x}=-x$, $\widetilde{y}=y$, and $\widetilde{z}=z$.

The result demonstrates that the transformation of the magnetic field components is always different than the transformation of the current density components. Set "a'' is not suitable because of the reasons outlined above. Set "b'' implies a zero normal magnetic field (at x=0) but has a nonzero normal current density which is also inconsistent with the force free condition (1).

ii)

Point (line) mirroring (at x=0) with respect to x=y=0: $\widetilde{x}=-x$, $\widetilde{y}=-y$, and $\widetilde{z}=z$, which implies $\partial_{\widetilde{x}}=-\partial_{x}$, $\partial_{\widetilde{y}}=-\partial_{y}$, and $\partial_{\widetilde{z}}=\partial_{z}$. This resulting set of boundary conditions implies a line symmetry along the line x=y=0.

Similar to the prior symmetry all positive definite MHD variables must be symmetric, e.g., $\rho(-x,-y,z)=\rho(x,y,z)$. The result for all variables is given in Table 2.

Table 2: Transformation properties for the transformation $\widetilde{x}=-x$, $\widetilde{y}=-y$, and $\widetilde{z}=z$.

Both set "a'' and set "b'' have the same transformation properties for the magnetic field and the current density and are thus suitable for the MHD equations satisfying the force-free condition. In the following we will employ the transformation set "b'' for the magnetic field extrapolation. Note that while for MHD models the symmetry set (a) is equally suitable but is lost if the Hall term is included in Ohm's law. The boundary condition has some implications regarding a preconditioning of the initial magnetic field. The symmetry also applies to the plane x=0, i.e., B_z(0,-y,z)=B_z(0,y,z) such that the magnetic field along the z axis must be symmetric with respect to y=0 which can be obtained through

$$B_{z}^{\rm precond}(0,y,z)=\frac{1}{2}\left(B_{z}(0,y,z)+B_{z}(0,-y,z)\right).$$ (6)

Applying this boundary condition to each boundary along x and y also implies that all corner point must have the same value for B_z

$$B_{z}^{\rm corner}=\frac{1}{4}\sum_{i={\rm corners}}B_{z}^{(i)}.$$ (7)

In addition the solution will produce a magnetic field where the z axis (and the corresponding midpoint lines for the other boundaries are magnetic field lines because the symmetry along these lines implies B_x,B_y=0. However, this preconditioning appears acceptable compared to method by Seehafer (1978) which requires B_z=0 on all x and y boundary planes.

2.2 Line symmetry and fourier expansion

In the x, y plane we assume an expansion for B_z with base functions of the form

 

$$\Psi_{mn} c_{1}\sin\frac{\pi m}{L_{x}}x\sin\frac{\pi n}{L_{y}}y+c_{2}\sin\frac{\pi m}{L_{x}}x\cos\frac{\pi n}{L_{y}}y$$ =

$$+c_{3}\cos\frac{\pi m}{L_{x}}x\sin\frac{\pi n}{L_{y}}y+c_{4}\cos\frac{\pi m}{L_{x}}x\cos\frac{\pi n}{L_{y}}y$$ (8)

in a system 0<x<L_x, -L_y/2<y<L_y/2 as illustrated in Fig. 1. The line symmetry conditions are

$$\begin{eqnarray*}\Psi_{mn}\left(-x,-y\right)=\Psi_{mn}\left(x,y\right) & & {\rm ... ...}}{2}+y\right)\\ & & {\rm at}~~ y=y_{\max}=\frac{L_{y}}{2}\cdot\end{eqnarray*}$$

Applying the symmetry conditions to the expansion yields the base functions

$$\Psi_{mn}=c_{1}\sin\frac{\pi m}{L_{x}}x\sin\frac{\pi n}{L_{y}}y+c_{2}\cos\frac{\pi m}{L_{x}}x\cos\frac{\pi n}{L_{y}}y$$ (9)

with m, n both even or m, n both odd for the expansion in x and y.

	Figure 1: Illustration of the geometry of line symmetry.
Open with DEXTER

  
2.3 Solution to the force free condition

The solution is found similar to the work by Seehafer using an expansion in the form:

 

$$B_{z,\alpha\beta} = c_{1}\exp\left(-\lambda_{\alpha\beta}z\right)... ...n\beta y +c_{2}\exp\left(-\lambda_{\alpha\beta}z\right)\cos\alpha x~\cos\beta y$$ (10)

and

  

$$B_{x\alpha\beta} = \exp\left(-\lambda_{\alpha\beta}z\right)\left(c_{1x}\sin\alpha x~\cos\beta y+c_{2x}\cos\alpha x~\sin\beta y\right)$$ (11)

$$B_{y\alpha\beta} = \exp\left(-\lambda_{\alpha\beta}z\right)\left(c_{1y}\sin\alpha x~\cos\beta y+c_{2y}\cos\alpha x~\sin\beta y\right).$$ (12)

Applying the force free condition allows to formulate the general solution which satisfies the line symmetric boundary conditions.

   

$$\sum_{m,n~{\rm odd}}\frac{c_{1mn}}{\kappa^{2}+\lambda_{mn}^{2}}\e... ...-\lambda_{mn}z\right)\left(\beta_{n}\kappa\sin\alpha_{m}x~\cos\beta_{n}y\right.$$                             B_x =

$$\left.-\alpha_{m}\lambda_{mn}\cos\alpha_{m}x~\sin\beta_{n}y\right)$$

$$+\sum_{m,n~{\rm odd}}\frac{c_{2mn}}{\kappa^{2}+\lambda_{mn}^{2}}\... ...a_{mn}z\right)\left(\alpha_{m}\lambda_{mn}\sin\alpha_{m}x~\cos\beta_{n}y\right.$$

$$\left.-\beta_{n}\kappa\cos\alpha_{m}x~\sin\beta_{n}y\right)$$ (13)

$$-\sum_{m,n~{\rm odd}}\frac{c_{1mn}}{\kappa^{2}+\lambda_{mn}^{2}}\... ...da_{mn}z\right)\left(\beta_{n}\lambda_{mn}\sin\alpha_{m}x~\cos\beta_{n}y\right.$$ B_y =

$$\left.+~\alpha_{m}\kappa\cos\alpha_{m}x~\sin\beta_{n}y\right)$$

$$+\sum_{m,n~{\rm odd}}\frac{c_{2mn}}{\kappa^{2}+\lambda_{mn}^{2}}\... ...\lambda_{mn}z\right)\left(\alpha_{m}\kappa\sin\alpha_{m}x~\cos\beta_{n}y\right.$$

$$\left.+\beta_{n}\lambda_{mn}\cos\alpha_{m}x~\sin\beta_{n}y\right)$$ (14)

$$\sum_{m,n~{\rm odd}}c_{1mn}\exp\left(-\lambda_{mn}z\right)\sin\alpha_{m}x~\sin\beta_{n}y$$ B_z =

$$+\sum_{m,n~{\rm odd}}c_{2mn}\exp\left(-\lambda_{mn}z\right)\cos\alpha_{m}x~\cos\beta_{n}y$$ (15)

with $\kappa^{2}+\lambda_{mn}^{2}=\alpha_{m}^{2}+\beta_{n}^{2}$, $\alpha_{m}=\pi m/L_{x}$, and $\beta_{n}=\pi n/L_{y}$. A given polarity image in B_z is expanded in the set of these base functions at the solar surface z=0 subject to the preconditioning of the image data according to (6) and (7). The extrapolated solution is then used as an initial condition in the full three-dimensional MHD simulation model.

It is important to note that neither the Seehafer (1978) nor the presented method require a total magnetic flux balance different from a fully periodic solution as used by Gudiksen & Nordlund (2002). The Seehafer (1978) solution is restricted in that it requires that the magnetic field has only a normal component at the boundaries in x and y while our extrapolation allows that all components of the magnetic field are non-zero at these boundaries except for the lines of symmetry.

3 Properties of the magnetic field extrapolation

The goal of the magnetic field extrapolation is to construct a suitable initial condition for an MHD simulation of the corresponding configuration. This does not necessarily imply that the initial configuration has to reflect the exact state of the magnetic field at the start of the simulation because this configuration undergoes continuous changes caused by the photospheric plasma motion. For this reason the magnetic and plasma structure can never be expected to be in an exact equilibrium configuration. However, after sufficient time has passed it can reasonably be expected that the precise initial state has only a minor influence on the dynamical evolution. In other words the system has a finite memory.

	Figure 2: Magnetogram of the line of sight components at the beginning of the $\Pi~$phase of an observed EUV bright point (Brown et al. 2001).
Open with DEXTER

Nevertheless, it is desirable that basic geometrical magnetic field properties are accurately reflected in the initial magnetic field. Boundary conditions always influence this geometry and the magnetic field evolution. Therefore it is important to choose a simulation domain of sufficient size such that the main magnetic field concentrations and their inter-connections are contained within the system. As an example for the model properties we used a magnetogram related to an EUV bright point described by Brown et al. (2001). The original magnetogram is shown in Fig. 2 with a pixel size of about 500 km. We note that a decomposition by 32 modes provides a highly accurate representation of the original magnetogram. However, for a system size of 28 Mm as in the case of the Brown et al. (2001) observation, the scale height of harmonics larger than 16 is less than 600 km, decreasing with increasing mode number. Thus the higher order harmonics do not contribute significantly to the magnetic field structure at transition region heights and in the solar corona.

	Figure 3: Representation of the magnetogram in Fig. 2 by 16 ( top) and by 8 Fourier modes ( bottom) with an overlay of magnetic field lines representing the potential field. Field lines are started from the photospheric boundary and projected onto this boundary.
Open with DEXTER

The Fourier decomposed B_z magnetic field component based on the solution in described in Sect. 2.3 for $\kappa=0$ is shown in Fig. 3. Here the magnetogram is resolved by 8 and by 16 Fourier modes for each direction. The plots also include the projection of magnetic field lines starting from the photospheric boundary where a reddish color indicates field lines extending higher into the photosphere. Both decompositions show the main magnetic polarities of the original data albeit the size of the major polarities increases with decreasing mode number and the magnitude of the magnetic field decreases with decreasing mode number. Although magnetic field lines have some similarity Fig. 3 demonstrates that the inclusion of higher harmonics describes more complex details of the field structure. Also more of the uniformly space field lines are short and close within the photosphere consistent with the small scale height of higher harmonics.

In comparison to our 8 mode solution Fig. 4 shows the Seehafer (1978) extrapolation with 8 modes. The two results show very similar distributions of the polarities and a similar shape of field lines close to the center of the system. However, field lines close to the boundary are rather different. The figure also demonstrates that all boundary field lines in the Seehafer solution are normal to the boundaries where they leave the domain as expected by the symmetry conditions discussed in Sect. 2. In fact the extrapolation after Seehafer (1978) allows only a normal component of the magnetic field at the boundary.

	Figure 4: Same as Fig. 3 for the Seehafer (1978) 8 mode extrapolation.
Open with DEXTER

The difference in the B_z component between our extrapolation and the Seehafer extrapolation is shown in Fig. 5 for the 16 mode expansion. The largest differences occur close to the boundaries as can be expected because of the different boundary conditions. The differences are also most pronounced in the highest order modes indicating that the difference decreases fast with height.

	Figure 5: Difference of the line of sight magnetic field between the Seehafer and our extrapolation for 16 modes.
Open with DEXTER

However, the comparison of the extrapolated magnetic fields does not necessarily imply that the magnetic structure is similar. A quantitative measure of this structure is difficult and the comparison of individual field lines is a poor indicator of the basic field structure because magnetic flux boundaries can be shifted by a small amount such that individual magnetic flux tubes started at the same physical location can have entirely different shape. Here we employ the differential magnetic flux tube volume

$$V_{B}=\int_{B}\frac{{\rm d}s}{\left\vert{\vec B}\right\vert}$$ (16)

as a measure of the magnetic configuration. The differential magnetic flux tube volume has the properties that it increases with the length of individual field lines. It also increases with decreasing magnitude of the field along a field line, such that a corresponding plot shows large flux tube volume in regions where field lines stretch to higher altitudes and where field lines enter regions of particularly weak magnetic field.

	Figure 6: Flux tube volume for our magnetic field extrapolation ( top) and Seehafer's model.
Open with DEXTER

The differential magnetic flux tube volume for our field model and for the Seehafer extrapolation is shown in Fig. 6 for the cases of 8 modes and 16 modes for our model and the 8 mode expansion for the Seehafer model. Comparison with Fig. 3 demonstrates that small flux tube volume is indicative for short flux tubes which close below the transition region. These are typically arranged along polarity inversion boundaries. Large flux tube volume indicates foot point locations of magnetic structures which rise to considerable height. Large flux tube volume can also be indicative for complex magnetic structure because the differential flux tube volume becomes very large (or diverges) close to magnetic neutral points or lines.

The comparison of the 8 mode expansions show good agreement in the central region of the system whereas there are considerable differences closer to the boundaries. Noticeable are very large flux tube volumes close to the boundaries in the Seehafer extrapolation. It turns out that these are actually an artifact of the Seehafer solution. Since the potential field in this expansion has only a normal component at the boundaries, any change of the polarity of the normal component implies a zero value of the magnetic field. Thus a reversal of the normal component implies the presence of artificial neutral lines at these boundaries. Figure 7 shows the normal magnetic field and the neutral lines for the Seehafer extrapolation at the x=0 boundary.

	Figure 7: Normal magnetic field for the Seehafer solution at x=0.
Open with DEXTER

A generic topological property of a three-dimensional field are isolated neutral points rather than neutral lines. Thus our model contains only isolated neutral points. The area of photospheric regions which map into the vicinity of neutral points is small. In contrast the presence of neutral lines (rather than neutral points) at the boundaries of the Seehafer extrapolation leads to the extended photospheric regions of large magnetic flux tube volumes illustrated in Fig. 6.

The 8 mode plot in Fig. 6 shows a relatively simple structure where areas with bluish color indicate longer field lines extending into the corona or mapping to the vicinity of neutral points. The case of 16 modes shows the same basic structure, however, with smaller dark regions which are cut out of the regions of higher flux tube volume. To understand this one has to remember that a larger number of modes implies the presence of small scale structure which also has a small scale height, i.e., which is not reaching into the corona. The strong magnetic field regions are more confined for cases with larger mode numbers. The dark cutouts in the flux tube volume reflect these small loops that do not reach into the corona and therefore have small flux tube volume. The same tendency, i.e., additional small size cutouts from the larger flux tube volume regions continuous in the case of 32 or more modes (not shown).

Another important question is the deviation of the tangential (B_xand B_y) components of the magnetic field, which arise in both, the Seehafer and our extrapolation methods although they start with the normal (B_z) component of the photospheric magnetic field alone. The reason is that the force-free condition $\kappa=\rm const.$requires for selfconsistency finite B_x and B_y components as well. In order to test our new extrapolation method we applied it to a case, where the full vector magnetogram information was available. As such we choose the active region NOAA AR8210 which was previously discussed and investigated in detail by Longcope (2004) and Welsch et al. (2004). We choose the moment 17:30 UT, just before a major eruption took place. As a quantitative measure for the deviation we calculated the relative tangential magnetic field deviation

$$\Delta B_{\rm t}\equiv\sqrt{\frac{\left(B_{x,{\rm ex}}-B_{x,{... ...\right)^{2}+\left(B_{y,{\rm ex}}+B_{y,{\rm meas}}\right)^{2}}}$$

where the lower index ex indicates the field components from the extrapolation (ours or Seehafer's) and the index meas indicates the magnetic field components obtained by polarimetry methods in the photospheric plane. Figure 8 shows the relative tangential magnetic field deviation for our extrapolation method (top) and for the Seehafer (bottom) method. In both extrapolations we assumed $\kappa=0$, i.e., potential fields. Both methods reveal similar regions with deviations in the tangential magnetic field components as compared to the observed ones. Because of the similarity we interpret most regions with large deviations as regions where larger field-aligned (and possibly perpendicular ) currents are present in the observed fields. Both methods also indicate larger deviations close to the boundaries which is expected because of the influence of the boundary conditions on the solution.

	Figure 8: Deviation of the tangential magnetic field between the extrapolated (our model - top; Seehafer - bottom) field and the vector magnetogram for the active region NOAA AR8210.
Open with DEXTER

4 Simulation model

The MHD equations are solved with a Leapfrog scheme which is second order accurate and has very low dissipation. A small dissipation is switched on if oscillations develop on the grid scale similar to flux corrected transport (FCT) schemes. Terms involving second order derivatives are treated with the Dufort-Frankel method which allows to consider very small resistivity.

The grid is chosen nonuniform in the z direction with the highest resolution at $z_{\min}$ and $\Delta z=0.3$ (corresponding to about 160 km). At least 8 grid points should be considered for the smallest wave lengths of the Fourier series expansion to ensure sufficiently small discretization errors for the finite difference scheme, i.e., for an area of 64 pixels modeled with a maximum of 8 modes per direction at least 64 grid points should be used for the simulation. In fact it turned out that 16 grid points are desirable for the shortest wavelength of the 16 mode expansion examples presented here. This corresponds to a $256\times256$ grid in the x, y plane. Since magnetic flux and energy decrease fast with with altitude, only modes which contribute significantly at coronal heights (z>4) are important for reconnection and magnetic reconfiguration in the corona.

4.1 Equations and normalization

An important consideration for a model of the photosphere/chromosphere/corona region is the interaction with the neutral environment. In the collision dominated photosphere and chromosphere this implies ionization and recombination, friction, and thermal contact between the plasma and neutrals. Radiative processes are not included in the present model. The plasma is coupled to a neutral background fluid which is present only in the lower (z<4) portion of the system representing photosphere and chromosphere.

The focus of the current implementation of the model is to provide a framework to simulate solar processes based on the extrapolation (Sect. 2) of observed photospheric magnetic field and to examine basic effects caused by the collisional environment below the transition region.

    

$$\frac{\partial\rho}{\partial t} -\nabla\cdot\rho{\vec u}-\mu(\rho-\rho_{0})$$ = (17)

$$\frac{\partial\rho{\vec u}}{\partial t} -\nabla\cdot\left[\rho{\vec u}{\vec u}+\frac{1}{2}\left(p+B^{2}\r... ...underline{\underline{1}}-{\vec B}{\vec B}\right]-\mu\rho({\vec u}-{\vec u}_{0})$$ = (18)

$$\frac{\partial{\vec B}}{\partial t} \nabla\times({\vec u}\times{\vec B}-\eta{\vec j})$$ = (19)

$$\frac{\partial h}{\partial t} -\nabla\cdot h{\vec u}-\frac{(\gamma-1)}{\gamma}h^{1-\gamma}\left(2(\gamma-1)\eta{\vec j}^{2}-\mu(h-h_{0})\right).$$ = (20)

The pressure p is substituted by the variable $h=(p/2)^{1/\gamma}$because this yields a continuity equation in the absence of source terms for the internal energy.

The above set of equations is in normalized units. We normalize magnetic field to $B_{0}=1~{\rm G}=10^{-4}~{\rm T}$, plasma density to the density in the corona just above the transition region $n_{0}=2\times10^{15}$ m^-3for protons (the actual number density is smaller but a presence of heavier ions and convenience let us choose this number), and length scales close to the length of an MDI pixel size L₀=500 km. This yields a normalization of velocities to the Alfvén speed $v_{\rm A0}=50$ km s^-1 and of times to Alfvén times $\tau_{\rm A0}=10$ s. Pressure is normalized to typical magnetic pressure $p_{0}=B_{0}^{2}/\left(2\mu_{0}\right)=4.0\times10^{-3}$and temperature to $T_{0}=p_{0}/\left(2n_{0}k_{\rm B}\right)=7.2\times10^{4}~{\rm K}$(the factor 2 is because of electrons and ions). This yields a typical thermal velocity of $v_{\rm th0}=\sqrt{2k_{\rm B}T_{0}/m}=\sqrt{B_{0}^{2}/\left(2\mu_{0}mn\right)}=v_{\rm A0}/\sqrt{2}$. Note that the physical temperature $T_{\rm p}$ is determined by $T_{\rm p}=T_{0}\sqrt{p/\rho}$with p and $\rho$ being the pressure and plasma mass density (assuming protons) in simulation units.

The ion neutral collision frequency is $\nu_{\rm in}=n_{n}\sigma_{n}v_{\rm th}$with $\sigma_{n}\approx10^{-19}$ m². In normalized units this yields $\mu=70\rho_{n}\sqrt{p/\rho}$, where $\rho_{n}$ is the neutral mass density normalized to $\rho_{0}$ (assuming hydrogen). This illustrates that the normalized collision frequency is large anywhere below the transition region (increasing toward the solar surface).

Electron - ion collision frequencies are about $\nu_{\rm ei}\approx10^{7}$ s^-1in the photosphere/chromosphere and $\nu_{\rm ei}\approx10^{3}$ s^-1in the lower corona. The normalized resistivity is $\eta=\lambda_{\rm e}^{2}\nu_{\rm c}\tau_{\rm A}/L_{0}^{2}$with the electron inertia length $\lambda_{\rm e}=c/\omega_{\rm pe}$. Note that the diffusion time is $\tau_{\rm diff}=L_{0}^{2}/(\lambda_{\rm e}^{2}\nu_{\rm c})$which is 10⁷ s in the photosphere/chromosphere and 10¹⁰ s in the corona. The corresponding normalized resistivity is between $\eta\approx10^{-6}$ and $\eta\approx10^{-9}$. However, the classical resistivity is far too small to explain observed time-scales of reconnection in the solar corona. In the simulation model we use a current dependent resistivity based on the drift speed of the current carriers. Our model assumes particle scattering as a result of plasma turbulence when the drift velocity surpasses a threshold value (such as for ion-acoustic or other micro-instabilities). Since these are sub grid processes it is assumed that the actual current are filamentary and on a scale of a few km, i.e. about 1/100 of our horizontal grid resolution. The resistivity is switched on when the parallel plasma velocity is larger then some critical value ( $0.2~~ v_{\rm A}(z_{\rm trans})~~\sqrt{~~ T(z)/T(z_{\rm trans}}$) where $z_{\rm trans}$ is the height of the transition region.

Note that for strong magnetic fields the resistivity and Alfvén time decrease with the inverse magnitude of the magnetic field. To obtain fast reconnection the thickness of the diffusion region has to decrease to about 1 km below the transition region (if fast reconnection were possible with the strong neutral coupling) and to about 1 m in the corona (and below these values if the magnetic field is larger than B₀=1 G).

4.2 Initial conditions

The value for the ion neutral collision frequency is chosen in part for numerical reasons (very large collision frequencies may require to reduce the time step). The mass density below the transition region is assumed to increase to $\rho=100$ (in units of the corona density which is set to 1) and the normalized collision frequency is assumed with a maximum of $\mu=30$ in the photosphere. This underestimates the actual collisional effects such as friction, however, whether the plasma is forced with 99% of the neutral velocity or with 99.99% is not significant for the physics. On the considered time scales the forcing is always large.

In the simulation we do not distinguish between effective frequencies for ionization/recombination, friction, and energy exchange. The main effect of the collisions is to keep plasma densities, velocities, etc close to equilibrium or neutral values. For this purpose the exact values do not matter as long as the frequencies are $\gg$1. The pressure is assumed constant with a value of 1 (this is equal to the plasma beta $\beta=1$ based on B₀=1 G) which yields a plasma $\beta\gg1$for all regions with ${\vec B}\gg1$.

The dynamical evolution in the simulation is usually initiated through a photospheric neutral vortex motion. This neutral velocity is chosen to satisfy $\nabla\cdot{\vec u}_{n}=0$. Since the plasma motion is strongly coupled to the neutral motion, this condition implies ${\rm d}\rho/{\rm d}t=0$ and thus avoids potentially unrealistic density perturbations. The neutral velocity is contained in x, y planes and implemented through potential U such that ${\vec u}_{n}=\nabla\times(U{\vec e}_{z})$. The neutral motion is specified throughout the simulation domain but is effective only where the neutral density and collision frequency are large. Thus the neutral motion forces the plasma to move in the same manner below the transition region.

4.3 Boundary conditions

At the x and y boundaries line symmetry with respect to the center of each boundary plane is assumed as discussed in Sect. 2.

At the lower boundary the tangential velocity is specified as defined by the neutral motion. The normal velocity is set to 0. The normal magnetic field uses $\nabla\cdot{\vec B}=0$ and the horizontal components are computed from $\nabla\times{\vec B}=\alpha{\vec B}$. As an alternative one could just fix the initial magnetic field for cases without the neutral flow or one can assume symmetry, i.e., $\partial_{z}B_{x}=0$and $\partial_{z}B_{y}=0$. Density and pressure are assumed symmetric (zero normal derivative). At the top boundary symmetric conditions (zero normal derivatives) are assumed.

5 Summary and discussion

There are many applications of the proposed modeling method. For instance the method can be used to identify the stability of solar magnetic field structures or to identify specific modes of magnetic reconnection. The presented model can employ inferred convection patterns for the foot point motion of the magnetic field in the photosphere. In the case of potential field this motion provides a pointing flux into the corona. In cases of a force free initial condition the resulting Poynting flux can be up- or downward.

	Figure 9: Magnetic field and magnetogram for the 16 ( top) and 8 mode ( bottom) expansions 99 min into the simulation. Field lines are computed from the plane z=5 just above the transition region.
Open with DEXTER

As an example for the dynamic evolution Fig. 9 shows magnetic field lines and the magnetogram 99 min into simulations for the 16 and the 8 mode expansions for the configuration illustrated in Fig. 2. Here field lines are computed from a plane at z=5 corresponding to about 2500 km, i.e., just above the transition region. The photospheric flow boundary condition have been presented in Büchner et al. (2005b). While Fig. 3 illustrates large differences in the magnetic field, the field structure is very similar in the plots of Fig. 9. The reason for this is the selection of field lines starting from the plane z=5. Since the amplitude of higher modes decreases rapidly with altitude (13)-(15) the magnetic structure of the 8 and 16 mode expansions are converging with increasing height. The difference between these expansions is largest at the photospheric boundary and for the current example the 8 mode expansions provides a reasonable approximation for the dynamics above the transition region.

	Figure 10: Extrapolated potential field for the 16 ( top) and 8 mode ( bottom) based on the observed at 15:36 UT by MDI line-of-sight magnetic field components, i.e., the same time as for the simulation results shown in Fig. 9.
Open with DEXTER

Our simulations demonstrate that currents are generated in the corona in response to the energy input through the photosphere . Indeed, we found that both parallel as well as perpendicular currents are created, which deform the magnetic fields compared to the potential field extrapolation. The two plots in Fig. 10 depict the potential fields, extrapolated based on the observed MDI line-of-sight magnetic field components at 15:36 UT, 99 min after the moment, for which the simulation was started (cf. Fig. 3) and at the same time that shows the simulation result in 9. The results of Figs. 9 and 10 demonstrate that the simulated magnetic field corresponds to a large extend to the magnetic field, extrapolated for the observed photospheric line of sight field. However there are many differences if one examines the detailed field structure which reflect the fact that the field from the simulation is not anymore a potential field and contains field-aligned as well as perpendicular currents. For instance, the elongated structural element, extending from the left main polarity to the lower right region, exists only in the simulated magnetic field structure. This structural element is due to the local current dissipation and reconnection, which we are able to describe by our simulation approach in contrast to any force-free extrapolation approach, which cannot reveal regions of perpendicular currents and reconnection.

Additional results on this case are presented in Büchner et al. (2005b). Those results demonstrate a good agreement between the TRACE observations of a corresponding EUV bright point and the field-aligned parallel potential. Other results using the presented method can be found in Büchner et al. (2005b,a). It is worth mentioning that model results frequently show large parallel electric fields forming in and just above the transition region, which demonstrates the importance of including this region into the MHD simulation.

The proposed method contains various assumptions and limitations, however, several of which can be overcome or can be improved upon. The main limitations are due to (a) the limited spatial resolution of observations and simulation; (b) the assumption of a force free equilibrium with a constant coefficient; (c) some of the physics included in the basic dynamical equations; and (d) the absence of newly emerging magnetic flux which is not implemented yet. Limitations in the third category address for instance radiative energy loss and heat conduction. The model development is being continued with the goal to incorporate additional important physics.

In summary the presented solar magnetic field expansion and simulation model provide a straightforward method which can incorporate observed solar magnetic fields into MHD simulations of the dynamics of magnetic structure. The model uses an extrapolation of the solar fields which is consistent with generic MHD boundary conditions. The model includes the photosphere, chromosphere, and solar transition region and results indicate that this transition region may be important for the dynamics of solar magnetic structure.

References

• Aly, J. J. 1989, Sol. Phys., 120, 19 [NASA ADS] [CrossRef]
• Amari, T., Aly, J. J., Luciani, J. F., Boulmezaoud, T. Z., & Mikic, Z. 1997, Sol. Phys., 174, 129 [NASA ADS] [CrossRef]
• Antiochos, S. K., Karpen, J. T., & DeVore, C. R. 2002, ApJ, 575, 578 [NASA ADS] [CrossRef]
• Brown, D. S., Parnell, C. E., Deluca, E., Golub, L., & McMullen, R. A. 2001, Sol. Phys., 201, 305 [NASA ADS] [CrossRef] (In the text)
• Büchner, J., Nikutowski, B., & Otto, A. 2004, in SOHO 15 Coronal Heating, ed. J. Ireland, & R. W. Walsh (Noordwijk: ESA Publ. Div.), ESA SP-575, 23
• Büchner, J., Nikutowski, B., & Otto, A. 2005a, in Proc. IAU Symp. 223, ed. A. V. Stepanov, E. E. Benevolenskaya, & A. G. Kosovichev, 353
• Büchner, J., Nikutowski, B., & Otto, A. 2005b, in Astrophysical Particle Acceleration in Geospace and Beyond, Geophysical monograph, ed. J. L. Horwitz, J. D. Perez, R. D. Preece, & J. Queenby (American Geophys. Union), 156, 161
• Cheng, C. Z., Ren, Y., Choe, G. S., & Moon, Y.-J. 2003, ApJ, 596, 1341 [NASA ADS] [CrossRef]
• Chiu, Y. T., & Hilton, H. H. 1977, ApJ, 212, 873 [NASA ADS] [CrossRef]
• Choe, G. S., & Cheng, C. Z. 2002, Physics of Plasmas, 9, 2330 [NASA ADS] [CrossRef]
• Cuperman, S., Ofman, L., & Semel, M. 1990, A&A, 230, 193 [NASA ADS]
• Galsgaard, K., Titov, V. S., & Neukirch, T. 2003, ApJ, 595, 506 [NASA ADS] [CrossRef]
• Gudiksen, B. V., & Nordlund, Å. 2002, ApJ, 572, L113 [NASA ADS] [CrossRef] (In the text)
• Linker, J. A., Lionello, R., Mikic, Z., & Amari, T. 2001, J. Geophys. Res., 25165
• Linker, J. A., Mikic, Z., Lionello, R., et al. 2003, Phys. Plasmas, 10, 1971 [NASA ADS] [CrossRef]
• Longcope, D. W. 2004, ApJ, 612, 1181 [NASA ADS] [CrossRef] (In the text)
• Mikic, Z., Linker, J. A., Schnack, D. D., Lionello, R., & Tarditi, A. 1999, Physics of Plasmas, 6, 2217 [NASA ADS] [CrossRef]
• Priest, E. R., & Forbes, T. 2000, Magnetic Reconnection: MHD Theory and Applications (Cambridge Univ. Press)
• Schmidt, H. U. 1964, in The Physics of Solar Flares, 107
• Seehafer, N. 1978, Sol. Phys., 58, 215 [CrossRef]
• Semel, M. 1968, Ann. Astrophys., 30, 513 [NASA ADS]
• Semel, M. 1988, A&A, 198, 293 [NASA ADS]
• Welsch, B. T., Fisher, G. H., Abbett, W. P., & Regnier, S. 2004, ApJ, 610, 1148 [NASA ADS] [CrossRef] (In the text)