A&A 429, 335-351 (2005)
DOI: 10.1051/0004-6361:20041507

Simulations of magneto-convection in the solar photosphere

Equations, methods, and results of the MURaM code

A. Vögler¹ - S. Shelyag¹ - M. Schüssler¹ - F. Cattaneo² - T. Emonet^3, - T. Linde³

1 - Max-Planck-Institut für Sonnensystemforschung, Max-Planck-Strasse 2, 37191 Katlenburg-Lindau, Germany
2 - Department of Mathematics, University of Chicago, 5734 South University Avenue, Chicago, IL 60637, USA
3 - Department of Astronomy and Astrophysics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA

Received 22 June 2004 / Accepted 23 August 2004

Abstract
We have developed a 3D magnetohydrodynamics simulation code for applications in the solar convection zone and photosphere. The code includes a non-local and non-grey radiative transfer module and takes into account the effects of partial ionization. Its parallel design is based on domain decomposition, which makes it suited for use on parallel computers with distributed memory architecture. We give a description of the equations and numerical methods and present the results of the simulation of a solar plage region. Starting with a uniform vertical field of 200 G, the processes of flux expulsion and convective field amplification lead to a dichotomy of strong, mainly vertical fields embedded in the granular downflow network and weak, randomly oriented fields filling the hot granular upflows. The strong fields form a magnetic network with thin, sheet-like structures extending along downflow lanes and micropores with diameters of up to 1000 km which form occasionally at vertices where several downflow lanes merge. At the visible surface around optical depth unity, the strong field concentrations are in pressure balance with their weakly magnetized surroundings and reach field strengths of up to 2 kG, strongly exceeding the values corresponding to equipartition with the kinetic energy density of the convective motions. As a result of the channelling of radiation, small flux concentrations stand out as bright features, while the larger micropores appear dark in brightness maps owing to the suppression of the convective energy transport. The overall shape of the magnetic network changes slowly on a timescale much larger than the convective turnover time, while the magnetic flux is constantly redistributed within the network leading to continuous formation and dissolution of flux concentrations.

Key words: magnetohydrodynamics (MHD) - Sun: magnetic fields - Sun: photosphere - Sun: granulation - Sun: faculae, plages

1 Introduction

The interaction between convective flows and magnetic fields in the solar photosphere and the uppermost layers of the convection zone gives rise to a number of processes which are of interest both in a solar and in a general astrophysical context. These include flux expulsion and convective field amplification, which lead to highly evacuated flux concentrations with superequipartition field strength, the modification or suppression of convective energy transport by strong fields, and the dynamics resulting from wave excitation and magnetic reconnection. In the context of the magnetic connectivity of the entire solar atmosphere, information on the structure and dynamics of the photospheric fields represents a crucial input.

While the progress in observational techniques (e.g. infrared spectroscopy, speckle polarimetry, adaptive optics) has greatly improved our knowledge of photospheric magnetic fields, observations provide information about the variation of physical quantities along the line of sight in a structured atmosphere only in highly convoluted form, which does usually not allow an unambiguous interpretation, while subphotospheric layers are entirely inaccessible to direct observations. Numerical simulations of MHD processes in the photosphere, on the other hand, provide a way to obtain information on the full three-dimensional structure of magnetic field configurations, as well as on spatial scales which are not resolved by current observations, thus helping to develop consistent models of the physical processes underlying the observed phenomena. In order to approximate solar conditions, comprehensive simulations of magneto-convection in the photosphere and convection zone include elaborate physics: radiative transfer, which is the main driver of convection and has an important influence on the temperature structure and brightness of magnetic field concentrations, and partial ionization, which strongly affects the efficiency of convective energy transport. With continuum, spectral line and polarization diagnostics, realistic numerical models can be directly compared with observations.

Following the work by Nordlund (1983), fully three-dimensional simulations of photospheric magneto-convection (e.g. Nordlund & Stein 1990; Bercik 2002; Stein & Nordlund 2003b) have provided valuable insight into the spatial structure and time evolution of magnetic field concentrations. A different approach was taken by Grossmann-Doerth et al. (1998) and Steiner et al. (1998) who carried out two-dimensional simulations of magnetic structures in Cartesian geometry. While geometrically constrained, these simulations achieved relatively high spatial resolution at affordable computational cost, thus permitting a more realistic modelling of small-scale magnetic structures, boundary layers and shocks. Comprehensive reviews of realistic simulations of solar magneto-convection have been given by Schüssler (2001) and Schüssler & Knölker (2001). Following the lines of more idealized models, Cattaneo (1999), Emonet & Cattaneo (2001) and Cattaneo et al. (2003) have carried out simulations of thermal magneto-convection in the Boussinesq approximation which highlighted the possibility that a substantial fraction of the magnetic field in the quiet photosphere is generated locally by fast dynamo action associated with the granular and supergranular flows.

In this paper we present MURaM, a newly developed 3D MHD simulation code for applications in the solar photosphere and convection zone which meets the requirements of realistic simulations (Vögler 2003; Vögler et al. 2003). It includes non-local and non-grey radiative transfer and takes into account partial ionization effects. The code is parallelized by means of domain decomposition, which allows to make full use of the computational capabilities of large, massively parallel computers with distributed memory architecture. We give a detailed description of the code and discuss simulation results obtained with it. Early other results already obtained with the MURaM code have been presented by Vögler & Schüssler (2003), Schüssler (2003), Schüssler et al. (2003), and Keller et al. (2004). The paper is organized as follows. In Sect. 2 we describe the physical model, including the basic equations and the boundary conditions. The numerical methods are addressed in Sect. 3. In Sect. 4 we present the simulation of a typical solar plage region with an average field strength of 200 G. We discuss morphology, time evolution and statistical properties of the magnetic field as well as the relation between field strength and brightness of magnetic structures. Our conclusions are given in Sect. 5.

2 Equations and boundary conditions

2.1 The MHD equations

We solve the equations of (non-ideal) MHD in three spatial dimensions in an inertial frame of reference with constant gravitational acceleration, using Cartesian coordinates. Rotation of the Sun is not accounted for since the Coriolis force is negligible on typical length-scales of solar granulation of the order of 1000 km. An estimate for the Rossby number gives  $Ro \equiv u/2\Omega L \simeq 300$, whereas the Coriolis force becomes dominant for $Ro \ll 1$. We write the system of magnetohydrodynamic equations in conservative form in terms of density, $\varrho$, momentum density,   $\varrho{\vec u}$, total energy per volume, e, and magnetic field strength, ${\vec B}$, as the independent variables. The system consists of the continuity equation,

$$\frac{\partial\varrho}{\partial t} + \nabla\cdot(\varrho{\vec u}) = 0 ~$$ (1)

representing mass conservation, the equation of motion,

$$\frac{\partial\varrho{\vec u}}{\partial t} + \nabla\cdot \... ...i} \right]= \varrho {\vec g} + \nabla\cdot{\underline{\tau}},$$ (2)

the energy equation,

 

$$\frac{\partial e}{\partial t} + \nabla\cdot\left[ {\vec u}\left(e... ...}\vert^2}{8\pi}\right) -\frac{1}{4\pi}{\vec{B}}( \vec{u}\cdot \vec{B}) \right]=$$ (3)

$$\frac{1}{4\pi}\nabla\cdot ({\vec B}\times \eta\nabla\times {\vec ... ...line{\tau}}) + \nabla\cdot(K\nabla T) +\varrho ({\vec g\cdot u}) + Q_{\rm rad},$$

and the induction equation

$$\frac{\partial{\vec B}}{\partial t} + \nabla\cdot({\vec{uB}-\vec{Bu}}) = -\nabla\times (\eta\nabla\times {\vec B}).$$ (4)

In Eqs. (1)-(4), u is the flow velocity, p is the gas pressure and g is the gravitational acceleration. ${\vec{uu}}$,  ${\vec{uB}}$,  ${\vec{Bu}}$ and   ${\vec{BB}}$ are dyadic products and  $\underline{{\vec 1}}$is the $3\times 3$ unit matrix. ${\underline{\tau}}$ is the viscous stress tensor. The energy equation is written for the total energy density per volume, e, which is the sum of internal, kinetic and magnetic energy densities:  $e=e_{\rm int}+ \varrho \vert{\vec u}\vert^2/2+ \vert{\vec B}\vert^2/8\pi$. T is the temperature and K the thermal conductivity. $Q_{\rm rad}$ is the source term which accounts for radiative heating and cooling. This term is discussed in detail in Sect. 2.2.

In our code, the diffusive terms are implemented in two alternative ways. One version assumes constant scalar diffusion coefficients, $\mu$ (the dynamic viscosity), K and  $\eta$. In this case, the components of the vicous stress tensor  $\underline{\underline{\tau}}$are given by

$$\tau_{ij} = \mu\left(\frac{\partial u_{i} }{\partial x_j} +\... ...eft(\nabla\cdot{\vec u}\right) \right), \quad\quad i,j=1,2,3.$$ (5)

The alternative implementation uses artificial shock-resolving and hyperdiffusivities, which provide efficient numerical stabilization and high effective Reynolds numbers in large eddy simulations, which do not aim at resolving the actual diffusive length scales. The numerical implementation is described in Sect. 3.2. In the simulations discussed in Sect. 4, the viscous and thermal diffusion terms were replaced by their artificial counterparts.

The system of MHD Eqs. (1)-(4) is completed by the equation of state, which describes the relations between the thermodynamical quantities of the fluid. At temperatures typically encountered in the uppermost few Mm of the convection zone, the solar plasma is partly ionized and the simple thermodynamical relations for an ideal gas do not apply. As hydrodynamical calculations have shown, changes in the thermodynamical properties due to partial ionization have considerable consequences for the character of convection. In the hydrogen ionization zone, about 2/3 of the enthalpy flux is transported by latent heat. Buoyancy driving is strongly enhanced since partial ionization reduces the adiabatic gradient  $\nabla_{\rm ad}$, and the increase in specific heat tends to suppress the radiative damping of temperature fluctuations (Rast & Toomre 1993). For our equation of state, we take into account the first ionization of the eleven most abundant elements. Since we write the MHD equation with $\varrho$ and  $e_{\rm int}$ as the fundamental thermodynamical quantities whose time development is explicitly described, the system is closed by relations specifying the dependence of T and p on $\varrho$ and  $e_{\rm int}$:

$$T=T(\varrho,e_{\rm int}); \quad\quad p=p(\varrho,e_{\rm int}).$$ (6)

These relations are stored in tables from which the required quantities are interpolated during a simulation run. The procedure by which we obtain these relations is sketched in Appendix A.

2.2 The radiative source term

 

The photosphere is the region where radiation takes over from convection as the dominant energy transport mechanism and where the plasma becomes transparent for radiation in the visible wavelength range. The energy exchange between gas and radiation determines the outgoing radiation intensity as well as the temperature structure of the photosphere and is responsible for the strong entropy gradient which acts as the main driver of convection. Therefore, an accurate representation of the radiative energy exchange rate,  $Q_{\rm rad}$, is crucial both for the dynamics and temperature structure as well as for the diagnostics of magneto-convection. Since the mean free path of photons becomes large as the atmosphere becomes transparent in the photosphere, radiative transfer at this height is essentially non-local, rendering the diffusion approximation of radiative energy transport inadequate. Hence, an accurate calculation of   $Q_{\rm rad}$ requires the solution of the (time-independent) radiative transfer equation (RTE hereafter)

$$\frac{{\rm d}I_\nu}{{\rm d}\tau_\nu} = S_{\!\nu}-I_\nu \;\;$$ (7)

for a number of frequencies and ray directions. $I_\nu$ is the (spectral) specific intensity, $S_{\!\nu}$ is the source function and  ${\rm d}\tau_\nu = \kappa_\nu~\varrho ~{\rm d}s$ is the optical thickness of the path element ${\rm d}s$. $\kappa_\nu$ is the opacity of the plasma. We neglect departures from local thermal equilibrium (LTE) and assume that the source function is given by the Planck function,  $S_\nu = B_\nu$. The numerical treatment of radiative transfer is dicussed in Sect. 3.5. Once the radiation field is known, the radiative energy flux,

$${\vec F}_{\nu}=\int_{4\pi} I_\nu(\vec{\mu})~{\vec{\mu}}~ {\rm d}\omega,$$ (8)

and the average intensity,

$$J_\nu= \frac{1}{4\pi}\int_{4\pi} I_\nu (\vec{\mu})~{\rm d}\omega,$$ (9)

can be calculated. Here  $\vec{\mu}(\omega)$ denotes the unit vector in the ray direction. The radiative heating rate then follows from

$$Q_{\rm rad} = -\int_\nu(\nabla\cdot{\vec F}_{\nu})~{\rm d}\nu \;=\; 4\pi\varrho\int_\nu\kappa_\nu (J_\nu-B_\nu)~{\rm d}\nu.$$ (10)

   
2.3 Boundary conditions

2.3.1 Upper boundary

In the current vesion of the code, the upper boundary of the computational domain is assumed to be closed, with stress-free conditions for the horizontal velocity components. This assumption leads to unphysical reflection of waves and shocks. In our simulations, however, the upper boundary is located near the temperature minimum, where densities are very low. Therefore, possible inaccuracies in the modelling of the flow field in the uppermost layers of the computational domain due to the closed-box assumption have only a small influence on the deeper layers around  $\tau_{500} =1$ which we are primarily interested in. For a realistic modelling of wave propagation near the upper boundary, an open and transmitting boundary is required. The implementation of such a boundary condition is under development. Physically, a closed boundary implies that the vertical convective fluxes of mass, energy, and horizontal momentum vanish on the boundary itself. This is achieved by setting the vertical velocity and the vertical gradients of mass density and energy density to zero:

$$u_z\;\vert _{\rm top}\equiv 0, \qquad \left. \frac{\partial \... ...frac{\partial e}{\partial z}\;\right\vert _{\rm top} \equiv 0.$$ (11)

The stress-free boundary conditions read

$$\left. \frac{\partial u_{x}}{\partial z}\;\right\vert _{\rm t... ...{\partial u_{y}}{\partial z}\;\right\vert _{\rm top} \equiv 0.$$ (12)

The treatment of the pressure gradient at the upper boundary is described in Sect. 3.4. For the magnetic field we assume that the field lines are vertical at the boundary, i.e.

$$B_{x,y}\;\vert _{\rm top}\equiv 0.$$ (13)

By virtue of solenoidality of ${\vec B}$, this implies

$$\left. \frac{\partial B_{z}}{\partial z}\;\right\vert _{\rm top} \equiv 0.$$ (14)

With this boundary condition free footpoint movement of the field lines is allowed.

2.3.2 Lower boundary

The lower boundary of the computational domain is located in the convectively unstable upper layers of the convection zone which are characterized by bulk fluid motions. Ideally, an implementation of an open lower boundary should allow free motions across the boundary without interfering with the fluid. Some degree of interference is unavoidable, though, since in principle any formulation of an open boundary condition requires knowledge of the physical conditions outside the domain, which is not available. The missing information has to be replaced with reasonable assumptions regarding the physical properties of the incoming fluid at the boundary.

The first assumption we make is that the total pressure $p_{\rm tot} = p + \vert{\vec B}^2\vert/8\pi$ is uniform across the lower boundary:

$$\left.p_{\rm tot}~\right\vert _{\rm bot} \equiv p_{{\rm tot},0}.$$ (15)

This is justified as long as fluid motions are slow compared to the speed of magneto-acoustic waves, $c_{\rm mag} = \left( c_{\rm Alf\acute{e}n}^2 + c_{\rm sound}^2 \right)^{1/2}$, since in this case fluctuations in  $p_{\rm tot}$ are balanced on a timescale much smaller than the flow timescale. For the magnetic field, the same vertical conditions as for the upper boundary are specified. Forcing the magnetic field lines to be vertical on the boundary implies that the net vertical magnetic flux in the computational box remains height-independent and constant in time. In downflow regions (u_z < 0) we assume a smooth isentropic outflow by setting

$$\left. \frac{\partial u_x}{\partial z}\;\right\vert _{\rm bot... ...rac{\partial u_z}{\partial z}\;\right\vert _{\rm bot} \equiv 0$$ (16)

and

$$\left. \frac{\partial s}{\partial z} \;\right\vert _{\rm bot} \equiv 0$$ (17)

where s is the entropy density per unit mass.

Owing to the strong stratification of the convection zone, upflows undergo a strong adiabatic expansion which tends to smooth out any initial or subsequent fluctuation during their rise. Therefore, one expects upflows entering the computational domain to be quite uniform. In our implementation, the internal energy density per unit mass, $\varepsilon_{\rm int}$, in upflows is assumed to be uniform across the lower boundary and set to a global value,  $\varepsilon_0$. We assume that all inflows are vertical with a smooth vertical profile of u_z:

$$u_x~\vert _{\rm bot} =u_y~\vert _{\rm bot} \equiv 0,\qquad ... ...ac{\partial u_z}{\partial z}\;\right\vert _{\rm bot} \equiv 0.$$ (18)

The parameter  $\varepsilon_0$ determines the inflow of energy into the computational box and thus can be used to control the net vertical energy flux through the domain. We adjust the value of  $\varepsilon_0$ in time in order to maintain the time-averaged outgoing radiation flux density at the top of the box at the solar value  $F_\odot = 6.34 \times 10^{10}~{\rm erg~s^{-1}~cm^{-2}}$. Since the rate at which energy is radiated does not react instantaneously to changes in  $\varepsilon_0$, the adjustment is done on the Kelvin-Helmholtz timescale  $\tau_{\rm KH}$ of the system, which is defined as the time over which the net energy flux through the domain equals the total internal energy of the system:

$$\tau_{\rm KH} = \frac{\int_{\rm box}~ e_{\rm int} ~{\rm d}V} {\int_{\rm top}F_\odot ~ ~{\rm d}x ~{\rm d}y }\cdot$$ (19)

At each timestep,  $\varepsilon_0$ is corrected according to

$$\varepsilon_0 \rightarrow \varepsilon_0\cdot \left(1 + \frac... ...rm KH}}\cdot \frac{F_\odot - F_{\rm top} }{F_\odot} \right),$$ (20)

where  $F_{\rm top}$ is the current value of the simulated radiation flux at the top of the domain and $\Delta t$ is the time increment.

Since our numerical simulations are carried out in a rather small computational box containing a limited number of convective flow cells, statistical fluctuations of the total mass in the simulated volume can become significant if no measures are taken to control the mass flux across the lower boundary. We use the total pressure at the bottom,  $p_{{\rm tot},0}$, as control parameter in order to preserve the total mass of the system. Let $\delta M$ be the (relative) mass deficit inherited from the previous timestep,

$$\delta M = \frac{M - M_0}{M_0},$$ (21)

where M is the current total mass and M₀the initial value. In order for the adjustment of  $p_{{\rm tot},0}$ to be smooth, we require that the deficit $\delta M$ be balanced on a timescale $\tau_M = 30$ s. For each timestep the deficit to be corrected is

$$\delta M_{\rm corr} = \delta M \cdot\frac{\Delta t}{\tau_M}\cdot$$ (22)

Based on the value of  $p_{{\rm tot},0}$ from the previous timestep, this deficit is balanced by adjusting the gas pressure in the upflow regions:

$$p_{\rm up} = p_{{\rm tot},0} + \Delta p.$$ (23)

The resulting changes in density and, consequently, in the vertical mass flux  $f_{\varrho} = \varrho~u_z$ at the bottom allow to efficiently control the mass influx. Subsequently, the newly found  $p_{\rm up}$ is used as the starting value  $p_{{\rm tot},0}$ for the following timestep. As a consequence, the downflow regions always lag behind one timestep in the pressure adjustment. Since the pressure correction for a single timestep is very small, this does not lead to the build-up of significant horizontal pressure gradients. The correct value for the pressure adjustment $\Delta p$is found by means of an iteration.

2.3.3 Side boundaries

All quantities are taken to be periodic in both horizontal directions.

3 Numerical methods

3.1 Spatial and temporal discretization

 The MURaM code solves the system of MHD equations on a three-dimensional uniform Cartesian grid. For the spatial discretization of the partial derivatives, centered, fourth-order accurate, explicit finite differences are used for both first and second derivatives. The code is parallelized using a domain-decomposition scheme. The computational domain is divided into a three-dimensional array of rectangular subdomains, each of which is assigned to a separate process on the computer. Parallelization is done with MPI and follows the distributed-memory concept, i.e. each process only possesses knowledge of the variables in its own subdomain. In order to apply the 5-point stencil of the fourth-order scheme to grid cells bordering on subdomain boundaries, knowledge of two layers of cells outside the subdomain is required. Therefore each subdomain is endowed with two layers of "ghost'' cells on each subdomain boundary. The required information is communicated between the neighbouring subdomains and stored in the ghost layers before the numerical derivatives are evaluated. If the subdomain boundary coincides with a global boundary, the ghost cells are assigned values according to the boundary conditions. The numerical solution of the system is advanced in time with an explicit, fourth-order accurate Runge-Kutta scheme. The centered-difference operators and the time-stepping scheme are described in Appendices B and C, respectively.

   
3.2 Artificial diffusivities

The low viscosity and corresponding high Reynolds number in the solar convection zone lead to structure formation down to a diffusive length scale of the order of centimeters. Since the computing power necessary to resolve these scales in a numerical simulation exceeds the capacity of contemporary computers by many orders of magnitude, all simulations of the solar photosphere and convection zone are necessarily large-eddy simulations, which simulate flows on resolvable scales and cut off the part of the energy spectrum which lies below the grid scale. In order to balance the cascading of energy to small scales due to the nonlinearity of the momentum equation, a large-eddy simulation requires some kind of numerical viscosity to prevent the build-up of energy at the grid scale. One way of handling this problem is based on numerical schemes which are highly non-diffusive (like e.g. higher-order centered-difference schemes as used in our code) and use explicit viscous terms to dissipate energy at the grid scale. We follow this approach and use the methods described by Stein & Nordlund (1998) and Caunt & Korpi (2001) as a guideline. The diffusive terms in the momentum and energy equations are replaced by artificial equivalents; likewise, in the continuity equation an artificial diffusive term is introduced, which has no physical counterpart. In the induction equation, we retain the fourth-order centered-difference term with constant $\eta$ and add artificial diffusion only near the lower boundary of the computational domain where we found it necessary to stabilize the numerical scheme.

For each physical quantity subject to diffusion and for each coordinate direction, a separate diffusion coefficient, consisting of a shock-resolving and a hyperdiffusive part, is defined:

$$\nu_l(u) = \nu_l^{\rm shk} + \nu_l^{\rm hyp}(u).$$ (24)

Here u stands for the quantity to be diffused and the index l indicates the coordinate direction.

The shock-resolving part is designed to have significant values only in those regions where converging flows with strong cell-to-cell velocity jumps lead to the build-up of strong gradients in advected quantities. The rate at which gradients grow in converging flows is determined by the local value of the flow divergence. Following Caunt & Korpi (2001), we define

$$\nu_l^{\rm shk} = \left\{ \begin{array}{cl} c_{\rm shk}\cd... ...uad & \qquad\nabla\cdot {\vec u} \ge 0 \end{array} \right..$$ (25)

Here,  $c_{\rm shk}$ is a scaling factor of order unity and  $\Delta x_l$ is the grid spacing in direction l. With the shock-diffusivity defined this way, the timescale for the build-up of gradients in shocks is approximately equal to the timescale of diffusion across the mesh width  $\Delta x_l$, which ensures that the solution remains resolved at the shock. We included the shock diffusivity in the momentum and energy equations. In the other equations it was not found to be necessary for stability and therefore not included.

The hyperdiffusive part is defined on cell interfaces normal to the direction of diffusion. For a physical quantity u and direction l, we define, following Stein & Nordlund (1998):

$$\nu_l^{\rm hyp}(u) = c_{\rm hyp}\cdot c_{\rm tot}\cdot \Delta... ...ac{{\rm max}_3\;\Delta_l^3 u}{{\rm max}_3\;\Delta_l^1 u} \cdot$$ (26)

Here  $c_{\rm hyp}$ is a scaling factor of order unity,  $c_{\rm tot}$ is the sum of flow velocity, speed of sound, and Alfvén velocity. At the interface  $i + \frac{1}{2}$ between cells i and i+1 (i being the grid index of direction l), the expressions  $\Delta_l^3 u$ and  $\Delta_l^1 u$ are defined as the third and first differences of u, respectively:

$$(\Delta_l^3 u)_{i + \frac{1}{2}} = \vert 3(u_{i+1}-u_{i}) - (u_{i+2}-u_{i-1})\vert$$ (27)

and

$$(\Delta_l^1 u)_{i + \frac{1}{2}} = \vert u_{i+1}-u_{i}\vert.$$ (28)

The symbol  ${\rm max}_3$ indicates the maximum over three adjacent interfaces. The expression  $\Delta_l^3 u /\Delta_l^1 u$ in Eq. (26), which is proportional to the ratio of third and first derivatives, detects small-scale fluctuations and leads to significant values of the hyperdiffusion only where numerical noise on the grid level needs to be dissipated while resolved structures remain largely unaffected. The factor  $c_{\rm tot} \cdot \Delta x_l$ results in diffusive timescales across the mesh width,  $\Delta x_l$, which balance the timescale on which noise on the grid scale grows as the result of the information exchange between neighbouring grid cells. The explicit form of the diffusion terms is given in Appendix D.

   
3.3 Calculation of the timestep

The maximum allowed timestep,  $\Delta t_{\rm max}$, is determined by a criterion which ensures that in a single timestep flow velocities, wave speeds, and diffusion transport information across distances not larger than the mesh width. When the artificial diffusion terms are used, the timestep criterion in our code reads

$$\Delta t_{\rm max} = {\rm min}(\Delta t_{\rm ad}, \Delta t_{\nu}).$$ (29)

Here

$$\Delta t_{\rm ad}= c_{\rm ad}\cdot \frac{{\rm min}(\Delta x, \Delta y, \Delta z)} {c_{\rm tot}}$$ (30)

is the advective timestep, which is determined by the usual CFL criterion. $c_{\rm tot}$ is the total wave speed:

$$c_{\rm tot} = u + c_{\rm sound} + c_{\rm Alfv\acute{e}n}.$$ (31)

Each coefficient of artificial diffusion imposes a timestep limit of the form

$$\Delta t_\nu = c_\nu \cdot \frac{\Delta x^2}{\nu},$$ (32)

which is based on the diffusion time across the mesh width, $\Delta x$. $c_{\rm ad}$ and $c_\nu$ are safety factors of order unity (typically 0.5-0.7).

In the case of constant scalar diffusion coefficients, a more strigent timestep criterion, based on the analysis of linear advection-diffusion equations discretized with fourth-order centered differences (e.g. Hirsch 1988), is used:

$$\Delta t_{\rm max} = {\rm min}(\Delta t_{1}, \Delta t_{2}),$$ (33)

with

$$\Delta t_{1} = \frac{8}{3}\cdot\frac{{\rm min}(\nu,\eta,\kappa)}{c_{\rm tot}^2}$$ (34)

and

$$\Delta t_{2} = \frac{2}{3} \left[ {\rm max}(\nu,\eta,\kappa... ...\Delta y^2} + \frac{1}{\Delta z^2} \right) \right]^{-1}\cdot$$ (35)

Here,  $\nu=\mu/\varrho$ and  $\kappa=K/c_p\varrho$.

   
3.4 Implementation of the boundary conditions

In order to set a quantity u to value u₀ on the upper or lower boundary, the values in the two layers of ghost cells are fixed in such a way that the difference, u-u₀, is antisymmetric with respect to the boundary,

$$u^{\rm gc}_{1,2}-u_0 = u_0-u^{\rm dc}_{1,2}\cdot$$ (36)

The superscripts dc and gc refer to domain cells and ghost cells, respectively, and the subscripts 1 and 2 denote, respectively, the first and second grid plane on either side of the boundary. A vertical derivative  $\partial_z u$ is set to zero by assigning the ghost cell values for u symmetrically with respect to the boundary:

$$u^{\rm gc}_{1,2} = u^{\rm dc}_{1,2}.$$ (37)

Consequently, the closed upper boundary condition implies that the ghost-cell values for   $\varepsilon$ and $\varrho$ are to be filled symmetrically (cf. Sect. 2.3.1). While the symmetric condition for   $\varepsilon$ is consistent with the generally small vertical temperature gradients encountered in the upper photosphere near the temperature minimum, the symmetric conditons for $\varrho$, and hence, p, are certainly not realistic since one expects a roughly exponential density and pressure stratification in this region. In order to avoid unphysical behaviour as a result of these assumptions, the vertical pressure gradient is treated separately at the upper boundary. We replace the fourth-order centered differences with first- and second-order expressions that do not make use of the ghost cell pressure values:

$$\left[\frac{\partial p}{\partial z}\right]^{\rm dc}_1 = \fr... ...rm dc}_2 = \frac{p_1^{\rm dc}-p_3^{\rm dc}}{2\Delta z}\cdot$$ (38)

The lower formal order of the representation of the vertical pressure derivative affects only the uppermost two grid planes.

3.5 Numerical treatment of radiative transfer 

   
3.5.1 Treatment of frequency dependence

We use the multigroup method (Nordlund 1982; Ludwig 1992; see also Stein & Nordlund 2003a) in order to account for the frequency dependence of the radiative transfer. The underlying idea of this approach is to sort frequencies into a small number of groups according to the height at which they mainly contribute to the radiative heating rate. As the sorting criterion we use the geometrical depth in a 1D reference atmosphere at which optical depth  $\tau_\nu=1$ is reached. For each frequency group, an average opacity is defined and a transfer equation for the group-integrated intensity with a group-integrated source function is solved. For the simulation described in Sect. 4 we have used four opacity groups. We took the opacity distribution functions and continuum opacities from the ATLAS9 stellar atmosphere package of Kurucz (1993) as the basis for the multigroup sorting and averaging procedure. A comprehensive description of our multigroup implementation is given by Vögler et al. (2004).

   
3.5.2 Integration scheme

Since the computation of the radiative source term is usually the most time-consuming part of a realistic simulation, efficiency considerations have high priority in the choice of numerical methods. Furthermore, the domain decomposition scheme used in our code demands a solver which is local and reduces the communication overhead to a minimum. Given these requirements, the method of choice is the short-characteristics formal solver (Mihalas et al. 1978; Olson & Kunasz 1987; Kunasz & Olson 1988; Kunasz & Auer 1988). It is based on the discretized form of the formal solution of the radiative transfer equation. In order to calculate the intensity for a certain direction and frequency on a given grid point, the transfer equation is solved along the ray segment (short characteristic) between the grid point and the nearest upwind intersection of the ray with a cell boundary. The grid on which we solve the radiative transfer coincides with the cell corners of the grid used by the MHD solver. The values of T, p and $\varrho$ on the radiative grid are interpolated from the adjacent MHD grid cells. An example for a short characteristic within a grid cell is shown in Fig. 1. Applying the formal solution to the short characteristic between the points F and E, the intensity at point F is given by

$$I_{\rm F} = {I_{\rm E}~{\rm e}^{-\Delta\tau_{\rm EF}} } + {... ..._{\rm E}}B(\tau) ~{\rm e}^{\tau_{\rm F} - \tau}~ {\rm d}\tau}.$$ (39)

Here $\tau$ measures the optical depth along the ray, starting from point F, and

$$\Delta\tau_{\rm EF} = \tau_{\rm E} - \tau_{\rm F} =\int_{\rm F}^{\rm E} \kappa (s)\varrho (s) {\rm d}s.$$ (40)

Unless the ray direction is aligned with one of the coordinate directions (or one of the diagonals of the grid cell), the point of intersection E is located on a cell interface and the intensity $I_{\rm E}$ needs to be interpolated from the grid points at the interface corners (points A to D in Fig. 1). On the 3D grid used here, we calculate $I_{\rm E}$ by bilinear interpolation on the rectangle ABCD. As pointed out by Kunasz & Auer (1988), the short-characteristics scheme with linear interpolation of upwind intensities leads to an artificial broadening of rays. For the calculation of angle-integrated quantities, this behaviour is not entirely undesirable since, in this context, the intensity for a given ray direction is supposed to represent the radiation coming from a cone of directions. For the evaluation of the integral contribution, $\int_{\rm F}^{\rm E}B(\tau) ~{\rm e}^{\tau_{\rm F} - \tau}~ {\rm d}\tau$, and the optical depth interval,  $\Delta\tau_{\rm EF}$, we adopt the linear method described by Kunasz & Auer (1988): density, opacity and source function are approximated as linear functions along the ray segment  $\rm\overline{EF}$, which leads to analytical expressions for  $\Delta\tau_{\rm EF}$and the integral.

	Figure 1: The intensity at gridpoint F is obtained by solving the transfer equation along the short characteristic  $\overline {\rm EF}$. The intensity at the upwind point E is interpolated from the (already known) intensity values at the surrounding gridpoints, A to D.
Open with DEXTER

	Figure 2: The walking order of the short-characteristics method in a 2D grid for a ray direction pointing into the upper right quadrant. Black circles represent gridpoints on the upwind boundaries, where the intensity values are assumed to be given.
Open with DEXTER

Since the interpolation of the upwind intensity requires knowledge of the intensity at the four surrounding gridpoints (points A to D in Fig. 1), the grid must be traversed in a sequence which makes sure that the required upwind information is always available. To this end, for a given ray direction, the scheme starts in each subdomain at those boundaries through which the radiation enters (the upwind boundaries). The intensity values at these boundaries are assumed to be given. Then the traversal of the subdomain systematically proceeds in the downwind direction, propagating the boundary information across the grid (see Fig. 2 for a 2D example). Since the correct initial values on the upwind boundaries of a subdomain are a priori unknown unless these boundaries coincide with the top or bottom of the computational box, this procedure must be iterated until convergence on the boundaries is obtained. The intensities at a given upwind boundary are updated after each iteration with the new values provided by the neighbouring subdomain. Clearly, the number of iterations required depends on the accuracy of the initial guess. We use a linear extrapolation of the boundary values of the previous two timesteps. With this choice, on average 2-3 iteration steps per frequency and ray direction are sufficient to keep the relative error in intensity below 10^-3.

At the global top and bottom boundaries of the computational domain, the incoming intensity must be specified in order to start the short-characteristics scheme. As long as the medium is optically thin at the top of the box for a given frequency group $\nu$, the incoming radiation from outside the computational domain is negligible and the boundary condition

$$I_\nu({\vec{\mu}})\left.\right\vert _{\rm top} = 0 \qquad \forall~\mu_z < 0$$ (41)

can be used.

In the simulations discussed here, this assumption is valid for all frequency groups except the one representing the strongest line opacities. For this group, the  $\tau_\nu=1$ level is close to the top of the box and setting the incoming radiation to zero would lead to artificial cooling of the uppermost layers. In order to derive a more realistic boundary condition, we assume that above the computational domain there is a isothermal plane-parallel layer with temperature  $T_{\rm top}$, whose optical thickness in the most opaque group is  $\tau_{\rm top}$. According to the formal solution of the RTE, the radiation entering the box from such a layer is given by

$$I_\nu({\vec{\mu}}) \left.\right\vert _{\rm top} = B_\nu(T_{... ...\rm e}^{\tau_{\rm top}/\mu_z} \right) \qquad \forall~\mu_z <0.$$ (42)

In our simulations, we use this relation with  $\tau_{\rm top}=0.2$ and $T_{\rm top} = 4000 \;{\rm K}$ as boundary condition for the strongest opacity group. The choice  $\tau_{\rm top}=0.2$ corresponds to the optical depth $\tau_\nu$ of the strong-opacity group at similar geometrical height in 1D solar standard atmospheres.

The bottom of the simulation box is located in the optically thick regions where the diffusion approximation holds. At the bottom, incoming radiation is set to the local value of the source function:

$$I_\nu({\vec x},{\vec{\mu}})\left.\right\vert _{\rm bot} = B_\nu({\vec x}) \qquad \forall~\mu_z > 0 .$$ (43)

For each direction and each frequency group, the short-characteristics scheme described above is applied separately. For the angular integration we use the quadrature formulae of type A of Carlson (1963). In this scheme, the directions in one octant are arranged in a triangular pattern and the quadrature is invariant under rotations over multiples of $\pi /2$ around any coordinate axis. A summary of the construction procedure is given by Bruls et al. (1999). For the simulations in Sect. 4, the A4 quadrature set with three directions per octant was used.

Once the angular integration has been performed, the radiative heating rate for frequency group $\nu$ can be derived from the two equivalent expressions

$$Q^J_\nu = 4\pi\kappa_\nu\varrho~(J_\nu-B_\nu)$$ (44)

and

$$Q^F_\nu = -\nabla\cdot {\vec F}_\nu .$$ (45)

While the radiation field is defined on cell corners, the MHD solver requires cell-centered (or cell-averaged) values of the radiative heating rate. Cell-centered values for $Q^J_\nu$ are obtained by averaging over the values at the eight surrounding cell corners. In the case of $Q^F_\nu$, first average values for the components of  ${\vec F_\nu}$ are calculated on the centers of cell interfaces, then the derivatives are obtained by first-order finite differences between the averaged flux components on opposite cell interfaces.

As pointed out by Bruls et al. (1999), $Q^J_\nu$ suffers from severe accuracy problems in the optically thick regime, since then $J_\nu$ approaches $B_\nu$, and the difference of these two almost equal quantities is amplified by a factor  $\kappa_\nu~\varrho$which grows exponentially with depth. On the other hand, $Q^F_\nu$ is less accurate than $Q^J_\nu$ in the optically thin layers of the upper photosphere since small inaccuracies in the orientation of the flux vector can lead to significant errors when the divergence of the nearly constant radiative flux is determined, while the difference  $J_\nu-B_\nu$ is only slightly affected. Following a suggestion of Bruls et al. (1999), we make a smooth transition from $Q^J_\nu$ to $Q^F_\nu$ for each frequency group separately, depending on the local optical depth scale for the group considered, using $Q^F_\nu$ in regions with  $\tau_\nu \geq 0.1$ and $Q^J_\nu$ otherwise. The total (frequency-integrated) radiative heating rate is then given by

$$Q_{\rm tot} = \sum_\nu {\rm e}^{-\tau_\nu/\tau_0}~ Q^J_\nu +(1-{\rm e}^{-\tau_\nu/\tau_0})~Q^F_\nu$$ (46)

with $\tau_0=0.1$.

   
4 Simulation of a solar plage region

4.1 Simulation setup

	Figure 3: Geometrical setup of the simulation run. The vector  ${\vec B}_0$ indicates the vertical homogeneous magnetic field introduced into the hydrodynamic convection at the beginning of the magnetic phase.
Open with DEXTER

In this section we show some results of the simulation of a typical solar plage region. The simulations are part of a parameter study to investigate the properties of photospheric magneto-convection and their dependence on the average magnetic field strength. The dimensions of the computational domain are 1400 km in the vertical direction and 6000 km in each horizontal direction, with a resolution of  $100 \times 288 \times 288$ grid points. Initially, the simulation was set up as purely hydrodynamical convection, starting from a plane-parallel model of the solar atmosphere (Spruit 1974) extending from 800 km below to 600 km above the level of continuum optical depth unity at 500 nm, as initial condition (see Fig. 3). The height corresponding to  $\tau_{500} =1$ of the initial stratification was chosen as the zero level of the height coordinate z used in the following sections. For this simulation we used the artificial diffusivities described in Sect. 3.2 with  $c_{\rm shk} =1$ and $c_{\rm hyp}=0.03$. In order to stabilize the solution in the uppermost parts of the computational domain, $c_{\rm hyp}$ was increased in a layer of 200 km thickness below the top of the box and reaches a maximum value of 0.2 at the upper boundary. For the magnetic diffusivity we use the constant value $\eta = 1.1\times 10^{11}~{\rm cm^2~s^{-1}}$. Our choice of  $c_{\rm shk}$, $c_{\rm hyp}$, and $\eta$is motivated by test runs which showed that this set of parameters keeps the solution stable and well resolved without overresolving it. After convection had fully developed, a homogeneous vertical initial magnetic field of 200 G, corresponding to the average field strength in a strong solar plage region, was introduced.

4.2 Morphology and statistical properties

Figure 4: Map of (frequency-integrated) brightness ( lower right) and horizontal cuts at the the average geometrical height corresponding to optical depth unity of (clockwise from bottom left) temperature, vertical magnetic field and vertical velocity. The "mesoscale'' network of magnetic field structures is embedded in the network of granular downflows. Larger field concentrations appear dark while the brightness of small magnetic structures occasionally exceeds the brightness of granules. Most of the domain exhibits "abnormal'' granulation with small granules (compared to the "normal'' granules in the upper right corner). See also the movies provided as supplementary online material.

Open with DEXTER

The simulation results show that within a few minutes (approximately one turnover time of the convection) after the magnetic field has been introduced, the convective motions transport most of the magnetic flux into the intergranular downflow regions. During this initial phase, the magnetic field forms a network structure with maximum field strengths around 2000 G at height z=0. (It should be noted that the average  $\tau_{500} =1$ level of the fully developed convection in its statistically steady state is shifted upwards by approximately 80-100 km as a result of the turbulent pressure of the convective flow (see e.g. Stein & Nordlund 1998). The magnetic network is organized on a larger scale, which typically comprises 2-6 granules. As the simulation develops in time, the network on this scale turns out to be long-lived with a typical timescale of hours. Similar large-scale organization also appears in other convection simulations (e.g., Cattaneo et al. 2001). However, since the horizontal box size is comparable to the length scale of the pattern, it cannot be ruled out that it is an artifact.

Figure 5: Statistical properties of a layer of 100 km thickness around optical depth unity. Upper left: probability density function (PDF) of the field strength, signed with the vertical orientation of the field vector. Upper right: joint PDF of field strength and the inclination angle of ${\vec B}$ with respect to the horizontal, theta(${\vec B}$). Lower left: joint PDF of flow velocity, multiplied with the sign of its vertical component, and field strength. Lower right: joint PDF of the inclination angles of the flow, theta(${\vec v}$), and of the magnetic field, theta( ${\vec B}$). The grey-scaling indicates the probability density, with intervals of 0.5 on the $\log_{10}$ scale between greyscale levels.

Open with DEXTER

Figure 4 shows a map of the frequency-integrated emergent intensity (brightness) together with horizontal slices of temperature, vertical magnetic field, and vertical velocity at z=0 for a snapshot taken about two hours of simulated solar time after the start of the magnetic phase. The magnetic field forms elongated, sheet-like structures that extend along intergranular lanes as well as larger structures with a size of up to 1000 km ("micropores''), which are located at vertices where several downflow lanes merge. The micropores appear dark in the intensity picture, while smaller structures are usually brighter than the non-magnetic downflow lanes, their brightness occasionally exceeding the brightness of granules. A large part of the simulated area shows "abnormal'' granulation (Dunn & Zirker 1973; Title et al. 1992) with reduced granule sizes compared to the "normal'' granules in the upper right corner of the intensity map. The micropores are far from homogeneous; they show considerable small scale intensity fluctuations, which are related to localized upflows in regions of reduced field strength. While the overall shape of the magnetic network is stable, the magnetic features show a strong time dependence on small scales as magnetic flux is incessantly redistributed within the network. Consequently, the typical lifetime of the micropores is smaller than the timescale associated with the magnetic network and is comparable to the granule lifetime. The magnetic network itself is embedded in the network of granular downflows. While convective motions are effectively supressed inside strong field features, downflows occur at their edges. Basically, this result is consistent with earlier MHD simulations (e.g. Deinzer et al. 1984a,b; Knölker et al. 1991) which showed that the influx of radiation into a magnetic element drives a baroclinic flow in form of a strong downflow at its edge. It is also consistent with the observational finding that observed Stokes-V profiles in plage regions show a distinct area asymmetry (Stenflo et al. 1984; Grossmann-Doerth et al. 1989; Solanki 1989; Sigwarth et al. 1999).

Figure 5 shows some statistical properties of a layer of 100 km thickness around $\langle \tau_{500} \rangle =1$. The probability density function (PDF) for the magnetic field, multiplied with the sign of its vertical component is shown in the upper left panel. There appears to be a superposition of two components. Most of the volume considered is occupied by weak field, the probability density dropping off approximately exponentially with increasing field strength. The distribution reveals a pronounced local minimum at B=0, indicating that magnetic fields, albeit mostly weak, permeate the whole volume, and strictly field free regions are avoided. Superimposed on this exponential distribution is a Gaussian "bulge'' (the high field strength wing shows the characteristic parabolic shape on a logarithmic scale) with a maximum around 1500 G, which reflects the sheet and micropore structures in the network of concentrated magnetic field. The joint probability density function (joint PDF) of magnetic field strength and inclination angle of the field vector with respect to the horizontal plane given in the upper right panel of Fig. 5 shows that most of the strong field (above the kilogauss level) is vertical and upward directed (which is the orientation of the initial field), presumably as the result of buoyancy acting on the partially evacuated magnetic structures. The inclination angle of weak fields is much more evenly distributed. With decreasing field strength, a slight preference for upward directed fields is observed. The joint PDF of the vertical magnetic field and the flow velocity multiplied with the sign of the vertical velocity component in the lower left panel of Fig. 5 (positive velocities correspond to upflows) shows the effect of strong fields on the fluid motions: while flow velocities up to  $8~{\rm ~km\;s^{-1}}$ can be found in regions of weak field, the velocities are reduced in magnetic structures with field strengths above 1000 G. Fluid motions are not completely supressed, however, since the predominantly vertical fields permit vertical flows. Downflows are preferred inside strong field features. The lower right panel of Fig. 5 shows the joint PDF of the inclination angles of magnetic field vector and flow vector with respect to the horizontal plane. The pronounced diagonal feature indicates that flow field and magnetic field are roughly aligned in most of the volume considered. In addition to this component, there is a clear correlation of (strong) vertical magnetic field with downflows. The picture obtained here exhibits some similarities with simulations of idealized Boussinesq magneto-convection (Emonet & Cattaneo 2001), which also show a preferentially vertical orientation of strong fields embedded in downdrafts and a decrease of the velocity amplitude with increasing field strength.

	Figure 6: Left: face-on view of a thin magnetic sheet. The spreading of the field lines in the upper photosphere is clearly visible. At a depth of approximately 300 km below the surface, the sheet is disrupted. The grey, horizontal plane in dicates the height z=0. Right: magnetic map of the sheet at z=0. The arrow shows the viewing direction of the field line plot on the left.
Open with DEXTER

	Figure 7: The same sheet as in Fig. 6, from a different angle. The isosurface  $\tau _{\rm Ross}=1$ is shown as shaded surface.
Open with DEXTER

The three-dimensional field line plot given in Fig. 6 (left panel) shows a face-on view of a thin magnetic sheet (cf. magnetic map at z=0 on the right panel). The structure appears coherent around the visible surface. At a depth of 200-300 km below z=0 the sheet is disrupted by vigorous convective flows, the coherence is lost and a large part of the flux appears in form of turbulent, more or less randomly oriented field lines. We find this to be a typical property of thin magnetic flux concentrations in our simulation. Figure 7 shows the same flux sheet from a different viewing angle. The depression of the visible surface ( $\tau _{\rm Ross}=1$) associated with the magnetic sheet is clearly visible. The field lines appear to be systematically twisted, possibly as a result of a horizontal shear near the surface. In contrast, the micropore shown in Fig. 8 appears much more coherent over the whole simulated height range. While the field undergoes fragmentation and forms more concentrated strands of flux underneath the surface, a significant part of the flux remains more or less vertical in the sub-surface layers.

	Figure 8: Three-dimensional field line plot ( left) and corresponding magnetic map ( right) of a micropore. The translucent plane shows the field strength at z=0; strong fields appear white.
Open with DEXTER

The upper panel of Fig. 9 shows that the total (gas + magnetic) pressure inside strong field concentrations (here defined as fields exceeding a height-dependent critical value,  $B_{\rm c}(z)$; at a given height z, $B_{\rm c}$ is defined such that the fields with  $B > B_{\rm c}$comprise 70% of the total flux at that height) is over a large height range in balance with the external gas pressure. So, despite of the dynamic time evolution of the strong fields in the downflow network and although the typical diameters of flux concentrations are not necessarily small compared against other relevant length scales like e.g. the pressure scale height, the pressure conditions are largely consistent with the simple picture of a thin flux tube or flux sheet in horizontal magnetohydrostatic pressure balance with its surroundings. The deviations become significant in the uppermost 300 km of the computational domain where the total pressure inside strong fields clearly exceeds the external pressure. This is plausible since the field concentrations spread out with height and become nearly vertical close to the upper boundary, thus giving rise to an inwards directed magnetic curvature force which is balanced by an outwards directed magnetic pressure force. Correspondingly, the actual average field strength of the field concentrations shown in the lower panel of Fig. 9, exceeds the value predicted by the thin flux-tube approximation in the height range above $z\approx 200{-}300$ km.

Owing to the strong decrease in gas pressure with height, the plasma-$\beta$ in strong flux concentrations shown in Fig. 10 drops by more than three orders of magnitude across the vertical extent of the simulation box. While it is of the order of unity near the visible surface  $\langle \tau_{500} \rangle =1$ ( $z\approx 100$ km), it drops below 0.1 near the top of the box, which shows that the gas pressure is largely irrelevant for the internal force balance of magnetic flux concentrations in the uppermost parts of the domain.

Figure 9: Upper panel: horizontally averaged gas pressure inside (solid) and outside (dashed) strong field regions, and total (gas + magnetic) pressure in strong field regions (dotted). (See the text for the definition of strong fields used here.) Lower panel: average (rms) field strength inside (solid) and outside (dotted) strong fields as a function of height. The dashed line shows the field strength which would balance the difference between the inside and outside gas pressures.

Open with DEXTER

	Figure 10: Average plasma-$\beta$ in strong-field regions as a function of height.
Open with DEXTER

	Figure 11: Horizontally averaged (rms) horizontal velocity in regions of strong and weak magnetic field as function of height.
Open with DEXTER

	Figure 12: Detailed view of the flow and magnetic structure of a flux concentration. Upper panel: vectors of horizontal velocity superimposed on a greyscale map of the vertical magnetic field at a height of 550 km above the visible surface. Lower panel: the same at a height of 100 km. The longest arrows correspond to a velocity of  $7\;{\rm~km~s^{-1}}$.
Open with DEXTER

In the upper part of the computational domain, the structure of horizontal flows undergoes a characteristic change with height. This is illustrated by Fig. 11, which shows horizontally averaged (root-mean-square) velocities of horizontal flows as a function of height, in strong-field and weak-field regions, respectively. Horizontal flows are stronger outside magnetic field concentrations than inside below z=400 km and reach a maximum around z=200 km, where the granular upflows turn over and converge horizontally towards the downflow lanes. Above z=400 km, a different picture is obtained: the rms horizontal velocity inside magnetic structures increases with height, while the velocities outside decrease significantly. As a result, horizontal flows have larger amplitudes inside magnetic fields. As the example shown in Fig. 12 illustrates, there is a correlation between flows inside magnetic structures at large heights and converging granular flows immediately outside magnetic structures at deeper levels near the visible surface around  $\tau_{5000}=1$. The anticlockwise whirl flow at z=550 km (upper panel) corresponds to a net circulation of granular flows around the magnetic element with the same orientation at z=100 km (lower panel). It is conceivable that the shear due to such surrounding flows excites torsional Alfvén waves propagating upward along field concentrations. While such waves in the uppermost parts of the computational domain cannot be described realistically with the present closed upper boundary conditions, future versions of our code with a transmitting upper boundary will allow us to thoroughly analyze this potentially important mechanism of mechanical energy transport and assess its relevance for the energy balance in the upper layers of the solar atmosphere.

4.3 Relation between brightness and magnetic field strength

It is a well established observational fact that the brightness of magnetic structures in the photosphere strongly depends on their size (Topka et al. 1997). While small field concentrations appear as bright points in continuum-intensity pictures, larger structures like pores and sunspots are dark as a result of the reduced convective energy transport in their interior (see e.g. Zwaan 1987, for an overview of the hierarchy of magnetic elements). A particular case is the brightening of magnetic concentrations in molecular lines like the CH lines in the "G band''. Schüssler et al. (2003) and Shelyag et al. (2004) have shown that the G-band brightening is caused by a depletion of CH molecules in the partially evacuated magnetic structures.

Figure 13: Upper panel: joint probability distribution of magnetic field strength and frequency-integrated brightness. The grey-shading indicates the probability density, with level-intervals of 0.5 on the $\log_{10}$-scale. Lower panel: probability distribution functions of the magnetic field strength inside the magnetic network ( $\vert B_z\vert > 500 \;{\rm G}$), for regions brighter and darker than average, respectively. The field strength distribution on the $\tau =1$ or  $\tau = 0.1$ levels would, however, lead to a somewhat different picture since in dark structures the strong Wilson depression results in higher observed field strengths.

Open with DEXTER

A comparison of the map of magnetic field strength with the intensity picture in Fig. 4 shows that the brightest magnetic structures are typically found in regions of high field strength of the order of 2000 G at z=0. However, not all strong-field features stand out as particularly bright; larger structures with extended regions of magnetic fields above 1000 G tend to be darker than average. This impression is confirmed quantitatively by Fig. 13. The left panel shows the joint PDF of magnetic field strength at z=0 (approximately 80-100 km below the average level of optical depth unity) and emergent intensity, based on simulation data of approximately one hour simulated time. While a large part of the magnetic regions has an intensity below the average value, the largest field strengths show a trend towards increased brightness, with maximum intensities reaching or even slightly exceeding those of the brightest granules. The right panel shows the probability distribution of magnetic field strength in bright (solid line) and dark (dashed line) magnetic regions, both curves being separately normalized to unity. While the fields with below-average brightness have a rather flat distribution up to fields of 1500 G, the distribution for the bright features shows a steep increase with increasing field strength and reaches a pronounced maximum at 1700 G.

	Figure 14: Vertical field strength at z=0 ( upper panel) and intensity map ( lower panel) for a sheet-like magnetic structure. The vertical lines indicate the position of the cuts shown in Fig. 15.
Open with DEXTER

	Figure 15: Vertical cut through the sheet-like magnetic structure shown in Fig. 14. Left panel: density (grey-scaled) and magnetic field vectors projected on the veertical plane. The longest vector correponds to a field strength of 2000 G. Right panel: temperature structure and radiative flux vectors. The solid lines indicate the level  $\tau _{\rm Ross}=1$ (for the Rosseland opacity average).
Open with DEXTER

Figures 14 and 15 provide a detailed view of a typical sheet-like structure in an intergranular lane which illustrates the connection between brightness and field strength. Figure 14 shows magnetic field and intensity maps. The vertical line marks the position of the vertical cuts through this structure shown in Fig. 15. As Fig. 15 illustrates, the partial evacuation and lower temperature of the magnetic structure leads to a depression of the visible surface inside the sheet (the bold solid lines, indicating the level  $\tau _{\rm Ross}=1$). This results in radiation flowing in from the hot neighbouring granules through the sidewalls, heating up the interior of the sheet. This horizontal heating is balanced by cooling due to radiative losses in the vertical direction. The resulting net radiative heating rate in the interiour of the sheet is very small which indicates that it is nearly in radiative equilibrium and convective energy transport is not relevant. The bright appearance of the sheet is, in turn, a direct consequence of the increased temperature around the (depressed) visible surface corresponding to this radiative equilibrium.

This example confirms the basic mechanism for brightness enhancement due to channelling of radiation which was already found in earlier, more simplified models of photospheric magnetic elements (e.g. Spruit 1976; Deinzer et al. 1984a,b; Knölker et al. 1991) and demonstrates that results of idealized, two-dimensional models are indeed relevant for the explanation of some aspects of three-dimensional magneto-convection. The nature of facular brightenings close to the solar limb was studied by Keller et al. (2004) whose analysis of photospheric MHD simulations showed that faculae originate from a thin layer within granules immediately below largely transparent magnetic flux concentrations.

It should be noted that radiative heating is not the only mechanism which could affect the temperature structure of magnetic structures. In principle, Joule dissipation near the edges of the magnetic sheet, where the gradients in field strength are large, can also heat the interior of the sheet. In the solar photosphere, this effect is expected to be small compared with the radiative heat influx because of the high electric conductivity of the plasma. Since the value of the conductivity in simulations of this kind is unrealistically low - an inevitable consequence of the limited spatial resolution - Joule heating tends to be overestimated (Schüssler 1986; Hirayama 1992). However, as Fig. 16 shows, the contribution of Joule dissipation to the total heating around $\tau =1$ is negligible, indicating that the effects of limited grid resolution do not seriously affect the temperature structure obtained from the simulation.

	Figure 16: Radiative and joule heating across a flux sheet. The solid line shows the heating rate due to horizontal influx of radiation energy into the sheet along the cut indicated in Fig. 14, at a depth of 80  km below z=0. The dashed line shows the joule heating occuring at the sheet boundaries at the same depth.
Open with DEXTER

4.4 Formation of a micropore

Figure 17 illustrates the formation of a micropore. At the beginning of the time sequence shown, a small granule is surrounded by magnetic field, which is embedded in a complex structure of downflows and forms a ring around the granule. Subsequently, the initially bright granule shrinks and becomes darker as the upflow velocity in the granule decreases; the surrounding magnetic flux converges to form a dark micropore with a diameter of approximately 1000 km. In its interior, vertical motions are strongly reduced, which implies that convective energy transport is almost completely supressed. The micropore is, however, punctuated by small, non-magnetic upwellings, which penetrate it from below and carry hot plasma to the surface. While the shape of the micropore changes continuously, it more or less retains its position. Roughly 15 min after the end of the sequence shown, most of the magnetic flux of the micropore has been transported away into the magnetic network and the pore has dissolved. Apparently, the size (amount of magnetic flux) of the micropores which form in this simulation is too small to structure the convective flow around them in a way that would stabilize them and extend their lifetime beyond about 1-2 convective turnover times. The micropore formation process shown in Fig. 17 is reminiscent of the result reported by Hirzberger (2003), who observed several bright points immediately before the appearance of a micropore at the same position. It is also in qualitative agreement with the magneto-convection simulations of Bercik (2002) (see also Stein et al. 2002).

	Figure 17: Time series of vertical magnetic field ( left; light shades indicate strong fields), intensity ( middle), and vertical velocity ( right; light and dark shades indicate up- and downflows, respectively) during the formation of a micropore. An initially bright, small granule shrinks, while the magnetic field surrounding it forms a small pore which appears dark in the intensity picture.
Open with DEXTER

5 Conclusions

We have developed a 3D radiative MHD code for realistic simulations in the solar photosphere and convection zone. It includes a detailed treatment of partial ionization effects as well as non-grey radiative transfer, which makes the numerical model sufficiently realistic to allow a direct comparison with observations. Owing to the high degree of parallelization, including the radiative transfer part, the code is able to use the capabilities of large, massively parallel computers with distributed memory architecture to simulate spatially extended domains at high resolution.

The first results obtained with this code are promising. The simulation of a solar plage region with an average field strength of 200 G is in agreement with the basic picture of solar magneto-convection derived from previous simulations. But the results also suggest that a wealth of new information can be obtained by combining full three-dimensionality with high spatial resolution: while the sizes, lifetimes, and field strengths of micropores found in our simulations are consistent both with observations (Topka et al. 1997; Hirzberger 2003) and simulations (Stein et al. 2002), they appear as internally highly structured features which show significant fluctuations of brightness and field strength around the mean values and are penetrated by hot, weakly magnetized upflows. The near-surface structure of the quasi two-dimensional flux sheets in the intergranular lanes is consistent with earlier two-dimensional models (e.g. Deinzer et al. 1984a,b; Knölker et al. 1991). Our simulations, however, show that slender flux sheets which appear as coherent features around optical depth unity, get disrupted by turbulent downflows beneath the surface and are, in fact, shallow phenomena. The occurrence of strong vortical flows inside magnetic structures in the upper photosphere, which are triggered by the advection of vorticity towards the footpoints of the field concentrations in the surface layers, suggests that important information about the transport of mechanical energy in the upper layers of the solar atmosphere can be obtained from simulations of this kind.

Acknowledgements

This work has been supported by the Deutsche Forschungsgemeinschaft under grant SCHU 500/7 in the framework of the Priority Research Program "Analysis and Numerics of Conservation Laws''. A.V. and S.S. were supported by the Max-Planck-Institute for Aeronomy through doctoral stipends. A.V. would like to thank the colleagues from the University of Chicago for their hospitality and financial support during a stay at the Department of Astronomy and Astrophysics at UofC.

References

• Bercik, D. J. 2002, Ph.D. Thesis (In the text)
• Bruls, J. H. M. J., Vollmöller, P., & Schüssler, M. 1999, A&A, 348, 233 (In the text)
• Carlson, B. G. 1963, in Methods in Computational Physics, ed. B. Alder, & S. Fernbach, 1 (In the text)
• Cattaneo, F. 1999, ApJ, 515, L39 [CrossRef] (In the text)
• Cattaneo, F., Lenz, D., & Weiss, N. 2001, ApJ, 563, L91 [CrossRef] (In the text)
• Cattaneo, F., Emonet, T., & Weiss, N. 2003, ApJ, 588, 1183 [CrossRef] (In the text)
• Caunt, S. E., & Korpi, M. J. 2001, A&A, 369, 706 [EDP Sciences] [CrossRef] (In the text)
• Deinzer, W., Hensler, G., Schüssler, M., & Weishaar, E. 1984a, A&A, 139, 426 (In the text)
• Deinzer, W., Hensler, G., Schüssler, M., & Weishaar, E. 1984b, A&A, 139, 435
• Dunn, R. B., & Zirker, J. B. 1973, Sol. Phys., 33, 281 (In the text)
• Emonet, T., & Cattaneo, F. 2001, ApJ, 560, L197 [CrossRef] (In the text)
• Grossmann-Doerth, U., Schüssler, M., & Solanki, S. K. 1989, A&A, 221, 338 (In the text)
• Grossmann-Doerth, U., Schüssler, M., & Steiner, O. 1998, A&A, 337, 928 (In the text)
• Hirayama, T. 1992, Sol. Phys., 137, 33 [CrossRef] (In the text)
• Hirsch, C. 1988, Numerical Computation of Internal and External Flows (John Wiley & Son Ltd) (In the text)
• Hirzberger, J. 2003, A&A, 405, 331 [EDP Sciences] [CrossRef] (In the text)
• Irwin, A. W. 1981, ApJS, 45, 621 [CrossRef] (In the text)
• Keller, C. U., Schüssler, M., Vögler, A., & Zakharov, V. 2004, ApJ, 607, L59 [CrossRef] (In the text)
• Knölker, M., Grossmann-Doerth, U., Schüssler, M., & Weishaar, E. 1991, Adv. Space Res., 11, 285 [CrossRef] (In the text)
• Kunasz, P. B., & Auer, L. 1988, J. Quant. Spectrosc. Radiat. Transfer, 39, 67 [CrossRef] (In the text)
• Kunasz, P. B., & Olson, G. L. 1988, J. Quant. Spectrosc. Radiat. Transfer, 39, 1 [CrossRef] (In the text)
• Kurucz, R. 1993, ATLAS9 Stellar Atmosphere Programs and 2 km s^-1 grid., Kurucz CD-ROM No. 13, Cambridge, Mass.: Smithsonian Astrophysical Observatory (In the text)
• Ludwig, H. G. 1992, Ph.D. Thesis, Universität Kiel (In the text)
• Mihalas, D., Auer, L. H., & Mihalas, B. R. 1978, ApJ, 220, 1001 [CrossRef] (In the text)
• Nordlund, Å. 1982, A&A, 107, 1 (In the text)
• Nordlund, Å. 1983, Solar and Stellar Magnetic Fields: Origins and Coronal Effects, IAU Symp., 102, 79 (In the text)
• Nordlund, Å., & Stein, R. F. 1990, in Solar Photosphere: Structure, Convection, Magnetic Fields, ed. J. O. Stenflo (Kluwer: Dordrecht), IAU Symp., 138, 191 (In the text)
• Olson, G. L., & Kunasz, P. B. 1987, J. Quant. Spectrosc. Radiat. Transfer, 38, 325 [CrossRef] (In the text)
• Rast, M. P., & Toomre, J. 1993, ApJ, 419, 224 [CrossRef] (In the text)
• Schüssler, M. 1986, in Small Scale Magnetic Flux Concentrations in the Solar Photosphere, ed. W. Deinzer, M. Knölker, & H. H. Voigt, 103
• Schüssler, M. 2003, in Third International Workshop on Solar Polarization, ed. J. Trujillo Bueno, & J. Sánchez Almeida, ASP Conf. Ser., in press
• Schüssler, M., Shelyag, S., Berdyugina, S., Vögler, A., & Solanki, S. K. 2003, ApJ, 597, L173 [CrossRef]
• Shelyag, S., Schüssler, M., Solanki, S. K., Berdyugina, S., & Vögler, A. 2004, A&A, submitted (In the text)
• Sigwarth, M., Balasubramaniam, K. S., Knölker, M., & Schmidt, W. 1999, A&A, 349, 941 (In the text)
• Solanki, S. K. 1989, A&A, 224, 225 (In the text)
• Spruit, H. C. 1974, Sol. Phys., 34, 277 [CrossRef] (In the text)
• Spruit, H. C. 1976, Sol. Phys., 50, 269 [CrossRef] (In the text)
• Stein, R. F., & Nordlund, Å. 1998, ApJ, 499, 914 [CrossRef] (In the text)
• Stein, R. F., Bercik, D., & Nordlund, A. 2002, Nuovo Cimento C Geophysics Space Physics C, 25, 513 (In the text)
• Stein, R. F., & Nordlund, Å. 2003a, in Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 519 (In the text)
• Stein, R. F., & Nordlund, Å 2003b, in Modeling Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray, ASP IAU Symp., 210, 169 (In the text)
• Steiner, O., Grossmann-Doerth, U., Knölker, M., & Schüssler, M. 1998, ApJ, 49 5, 468 (In the text)
• Stenflo, J. O., Solanki, S., Harvey, J. W., & Brault, J. W. 1984, A&A, 131, 333 (In the text)
• Title, A. M., Topka, K. P., Tarbell, T. D., et al. 1992, ApJ, 393, 782 [CrossRef] (In the text)
• Topka, K. P., Tarbell, T. D., & Title, A. M. 1997, ApJ, 484, 479 [CrossRef] (In the text)
• Vögler, A. 2003, Ph.D. Thesis, Göttingen University, http://webdoc.sub.gwdg.de/diss/2004/voegler/voegler.pdf (In the text)
• Vögler, A., & Schüssler, M. 2003, Astron. Nachr., 324, 399 [CrossRef] (In the text)
• Vögler, A., Shelyag, S., Schüssler, M., et al. 2003, in Modelling of Stellar Atmospheres, ed. N. E. Piskunov, W. W. Weiss., & D. F. Gray, IAU Symp., 210, ASP Conf. Ser., 157 (In the text)
• Vögler, A., Bruls, J. H. M. J., & Schüssler, M. 2004, A&A, 421, 741 [EDP Sciences] [CrossRef] (In the text)
• Zwaan, C. 1987, ARA&A, 25, 83 (In the text)

 

  
Online Material

   
Appendix A: The equation of state

In the following, a brief outline is given how temperature and pressure are derived from  $e_{\rm int}$ and $\varrho$if the first ionization of a number of elements is taken into account. For the simulations presented here, the eleven most abundant elements in the solar photosphere were included. The internal energy per mass unit  $\varepsilon_{\rm int} = e_{\rm int}/{\varrho}$ is given by

$$\varepsilon_{\rm int}=\frac{3}{2\varrho}\left(n_{\rm e}+n_{\rm a}\right)kT+\frac{1}{\varrho}\sum n_{\rm i}^{*}\chi_{\rm i},$$ (A.1)

where the sum runs over the particle species (labeled with index i), $n_{\rm i}^{*}$ is the number density of ions of type i, and  $\chi_{\rm i}$ is the corresponding ionization energy. $n_{\rm a}=\sum n_{\rm i}$ is the number density of atoms, $n_{\rm e}$ the number density of electrons. Defining the ionization degree,  $x_{\rm i}={n_{\rm i}^{*}}/{n_{\rm i}}$, and the relative abundance, $\nu_{\rm i}={n_{\rm i}}/{n_{\rm a}}$, Eq. (A.1) can be rewritten as

$$\varepsilon_{\rm int}=\frac{3kT}{2\mu_{\rm a} m_0}\left(1+\s... ...ac{1}{\mu_{\rm a} m_0}\sum x_{\rm i} \nu_{\rm i} \chi_{\rm i},$$ (A.2)

where  $\mu_{\rm a}~$ is the mean molecular weight of the neutral gas ( $\mu_{\rm a} = 1.29$ for solar composition) and m₀ is the atomic mass unit. The ionization degrees, $x_{\rm i}$, are determined by the set of Saha equations

 

$$\frac{x_{\rm i}}{1-x_{\rm i}}\sum x_{\rm i} \nu_{\rm i}=\frac{u_{... ...\left(2\pi m_{\rm e} kT\right)^{3/2}}{h^3} \exp\left({-\chi_{\rm i}/kT}\right).$$ (A.3)

The temperature dependence of the partition functions  $u_{\rm i1},u_{\rm i0}$ are obtained from the literature (e.g., Irwin 1981). For temperatures exceeding about 16 000 K, the first ionization for all elements considered is almost complete, so that the temperature dependence of the partition functions can be neglected. In order to obtain the temperature from $\varrho$ and  $e_{\rm int}$, the nonlinear system of Eqs. (A.2) and (A.3) is solved iteratively. Once the temperature is known, the gas pressure follows from

$$p=\left(n_{\rm e}+n_{\rm a}\right)kT=\frac{\varrho}{\mu_{\rm a} m_0}\left(1+\sum x_{\rm i} \nu_{i}\right)kT.$$ (A.4)

Temperature and pressure are stored on a  $(\rho,e_{\rm int})$-grid from which the required values are interpolated during a simulation run.

Table A.1: Relative abundances and ionization energies of first ionization for the eleven most abundant elements in the solar photosphere.

   
Appendix B: Fourth-order centered difference operators

Choosing i as the index denoting the grid position along a particular spatial direction, the first and second spatial derivatives of quantity u are given by

$$\left(\frac{\partial u}{\partial x} \right)_{i}~ =\frac{1}{12\Delta x}~(-u_{i+2} ~+8u_{i+1}~-8u_{i-1}~+u_{i-2})$$ (B.1)

and

$$\left(\frac{\partial^2 u}{\partial x^2} \right)_{i}\! =\frac... ...x^2}~(-u_{i+2} ~+16u_{i+1} ~-30u_{i}~+ 16u_{i-1} ~- u_{i-2}),~$$ (B.2)

respectively, where $\Delta x$ is the grid spacing in the coordinate direction considered.

For partial differential equations corresponding to a conservation law, the centered difference scheme is equivalent to a finite volume scheme with the numerical flux

$$f_{i+1/2} = \frac{7}{12}~\left[f(u_{i+1})+f(u_{i})\right] - \frac{1}{12}~\left[f(u_{i+2})+f(u_{i-1})\right]$$ (B.3)

defined on the cell interfaces, where f(u) is the flux function of the conservation law for quantity u. Therefore, any quantitity obeying a conservation law is exactly conserved numerically as long as the total integrated flux at the domain boundaries vanishes. It is straightforward to show that, for an induction equation with constant scalar $\eta$, this discretization formally conserves the discretized form of  $\nabla\cdot{\vec B}$as a result of the symmetry of the scheme.

In a multidimensional domain, mixed derivatives are calculated by straightforward successive application of Eq. (B.1) in different coordinate directions.

   
Appendix C: Time-stepping procedure

The numerical solution of the system is advanced in time with an explicit fourth-order Runge-Kutta scheme. Defining  ${\vec U}_0$ as the vector which describes the state of the system for time t₀, i.e.

$${\vec U_0} = (\varrho, \varrho{\vec u} ,e,{\vec B})(x,y,z, t_0),$$ (C.1)

the system of partial equations can be written as

$$\frac{\partial{\vec U}} {\partial t} = {\vec R}({\vec U}),$$ (C.2)

where the vector  ${\vec R}({\vec U})$ contains the spatial derivatives and source terms in the system of equations. Then the new state  ${\vec U}_1$ for time  $t_1 = t_0 + \Delta t$ is calculated in four steps:

$${\vec U}_{\frac{1}{4}}={\vec U}_0 + \frac{\Delta t}{4}~{\vec R}({\vec U}_0),$$

$${\vec U}_{\frac{1}{3}} ={\vec U}_0 + \frac{\Delta t}{3}~{\vec R}({\vec U}_{\frac{1}{4}}),$$

$${\vec U}_{\frac{1}{2}} = {\vec U}_0 + \frac{\Delta t}{2}~{\vec R}({\vec U}_{\frac{1}{3}}),$$

and finally

$${\vec U}_1 = {\vec U}_0 + \Delta t ~ {\vec R}({\vec U}_{\frac{1}{2}}).$$ (C.3)

   
Appendix D: Numerical implementation of the diffusive terms

Figure D.1: The numerical stencil for the artificial diffusivities. Empty circles indicate the staggered grid on which the hyperdiffusive coefficients are defined. Left and middle panel: grid and stencil for the terms  $\partial _x(\nu \partial _x u)$ and  $\partial _y(\nu \partial _y u)$, respectively. Right panel: stencil for the term   $\partial _y(\nu \partial _x u)$ with mixed derivatives.

With the artificial diffusion $\nu_l$ defined in Sect. 3.2, the diffusive terms in the continuity and energy equations are given by

$$\left( \frac{\partial \varrho}{\partial t} \right)_{\rm diff}... ... \nu_l(\varrho)\frac{\partial\varrho }{\partial x_l} \right)$$ (D.1)

and

$$\left( \frac{\partial e}{\partial t} \right)_{\rm diff} = \su... ...varrho\nu_l(T)\;\frac{\partial c_p T }{\partial x_l} \right),$$ (D.2)

respectively. The diffusive terms for the vector quantities u and B have a more complicated structure. In the momentum equations we have

$$\left( \frac{\partial \varrho {\vec u}}{\partial t} \right)_{\rm diff} = \nabla\cdot {\underline{\tau}},$$ (D.3)

where  ${\underline{\tau}}$ is the symmetrized stress tensor

$$\tau_{kl} = \frac{1}{2}\varrho\left(\nu_k(u_l)\frac{\partial ... ... +\nu_l(u_k)\frac{\partial u_k}{\partial x_l} \right)\cdot$$ (D.4)

In the equation for the total energy, the dissipated energy is taken into account by a viscous heating term

$$\left( \frac{\partial e}{\partial t} \right)_{\rm visc} = \nabla\cdot \left({\vec u}\cdot{\underline{\tau}}\right).$$ (D.5)

The artifical diffusive term in the induction equation is given by

$$\left( \frac{\partial {\vec B}}{\partial t} \right)_{\rm diff} = -\nabla\times\vec{\mathcal{E}}.$$ (D.6)

The vector quantity   $\vec{\mathcal{E}}$ is defined as

$$\vec{\mathcal{E}} = \left( \begin{array}{c} \nu_y(B_z)\part... ...tial_x B_y - \nu_y(B_x)\partial_y B_x\\ \end{array} \right).$$ (D.7)

The dissipated energy is accounted for in the energy balance and appears in the energy equation as an Ohmic heating term:

$$\left( \frac{\partial e}{\partial t} \right)_{\rm ohm} = \fra... ...pi} \nabla\cdot\left({\vec B}\times\vec{\mathcal{E}} \right).$$ (D.8)

All artificial diffusion terms are discretized with second-order centered differences. The staggered grid used for the artificial diffusion is shown in Fig. D.1. At grid point i, a diffusion term of the form  $\partial _x(\nu \partial _x u)$ is calculated as

$$\left( \frac{\partial u}{\partial t} \right)_{{\rm diff},i} =... ...i - \frac{1}{2}}\cdot\frac{u_{i}-u_{i-1}}{\Delta x} \right),$$ (D.9)

where  $\nu_{i + \frac{1}{2}}$ and  $\nu_{i - \frac{1}{2}}$ are the artificial diffusion coefficient defined on the cell interfaces normal to the coordinate direction of diffusion. In a term with mixed derivatives like  $\partial _y(\nu \partial _x u)$, the hyperdiffusion coefficient is interpolated from cell interfaces onto cell centers. The inner part of such a term is then discretized as

$$\left(\nu \frac{\partial u}{\partial x}\right)_{\rm i} = \l... ...ght)\cdot \left(\frac{u_{i+1}-u_{i-1}}{2\Delta x}\right)\cdot$$ (D.10)

Intensity.mpg

MagneticField.mpg

Temperature.mpg

Velocity.mpg