Radiative transfer with scattering for domain-decomposed 3D MHD simulations of cool stellar atmospheres

W. Hayek; M. Asplund; M. Carlsson; R. Trampedach; R. Collet; B. V. Gudiksen; V. H. Hansteen; J. Leenaarts

doi:10.1051/0004-6361/201014210

Free Access

Issue		A&A Volume 517, July 2010


Article Number		A49
Number of page(s)		18
Section		Stellar atmospheres
DOI		https://doi.org/10.1051/0004-6361/201014210
Published online		03 August 2010

A&A 517, A49 (2010)

Numerical methods and application to the quiet, non-magnetic, surface of a solar-type star

W. Hayek^1,2,3 - M. Asplund¹ - M. Carlsson³ - R. Trampedach^4,2 - R. Collet¹ - B. V. Gudiksen³ - V. H. Hansteen³ - J. Leenaarts⁵

1 - Max Planck Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany
2 - Research School of Astronomy & Astrophysics, Cotter Road, Weston Creek 2611, Australia
3 - Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029 Blindern, 0315 Oslo, Norway
4 - JILA, University of Colorado, 440 UCB, Boulder, CO 80309-0440, USA
5 - Astronomical Institute, Utrecht University, Postbus 80 000, 3508 TA Utrecht, The Netherlands

Received 5 February 2010 / Accepted 28 April 2010

Abstract
Aims. We present the implementation of a radiative transfer solver with coherent scattering in the new BIFROST code for radiative magneto-hydrodynamical (MHD) simulations of stellar surface convection. The code is fully parallelized using MPI domain decomposition, which allows for large grid sizes and improved resolution of hydrodynamical structures. We apply the code to simulate the surface granulation in a solar-type star, ignoring magnetic fields, and investigate the importance of coherent scattering for the atmospheric structure.
Methods. A scattering term is added to the radiative transfer equation, requiring an iterative computation of the radiation field. We use a short-characteristics-based Gauss-Seidel acceleration scheme to compute radiative flux divergences for the energy equation. The effects of coherent scattering are tested by comparing the temperature stratification of three 3D time-dependent hydrodynamical atmosphere models of a solar-type star: without scattering, with continuum scattering only, and with both continuum and line scattering.
Results. We show that continuum scattering does not have a significant impact on the photospheric temperature structure for a star like the Sun. Including scattering in line-blanketing, however, leads to a decrease of temperatures by about 350 K below $\log_{10}\tau_{5000} \la -4$ . The effect is opposite to that of 1D hydrostatic models in radiative equilibrium, where scattering reduces the cooling effect of strong LTE lines in the higher layers of the photosphere. Coherent line scattering also changes the temperature distribution in the high atmosphere, where we observe stronger fluctuations compared to a treatment of lines as true absorbers.

Key words: radiative transfer - stars: atmospheres - Sun: atmosphere

1 Introduction

The atmospheres of late-type stars form the transition from the opaque convective envelope to the interstellar medium. Hot rising plasma transports heat to the surface, becomes transparent and looses its entropy through radiative cooling. Gravity accelerates the cooled gas back into the star, carrying kinetic energy inward and forcing the convective flow. By taking over heat transport and removing entropy, the radiation field therefore indirectly drives convection (Stein & Nordlund 1998), making radiative and hydrodynamical processes equally important at the surface. Magnetic fields have strong impact on the higher atmosphere and cause local phenomena in the surface granulation, such as spots and pores.

The classical numerical models of cool stellar atmospheres in 1D focused on a detailed description of radiative transfer, with two prominent examples being the MARCS code (Gustafsson et al. 1975) and the ATLAS code (Kurucz 1979). Assuming a plane-parallel or spherical-symmetric stratification, they include only a rudimentary treatment of convective energy transport in cool stellar atmospheres. Subsequent updates of these models (e.g., Gustafsson et al. 2008; Kurucz 1996) benefit from the largely increased computational power, refining the treatment of the strongly wavelength-dependent line opacities. Newer codes, such as PHOENIX (Hauschildt et al. 1999) can also include departures from local thermodynamic equilibrium (LTE) in the radiative transfer computation and the absorber populations. The 1D models have not only provided growing insight into the physical environment at the surface of cool stars, but have also become a standard tool for chemical abundance analyses. The wide variety of applications includes studies of galactic chemical evolution and of the origin of the elements.

The advent of fully dynamic 3D surface convection simulations has enabled a much more realistic treatment of the hydrodynamical plasma flow, deepening our understanding of convection and eliminating the need for microturbulent and macroturbulent broadening in line formation computations (see, e.g., Nordlund et al. 2009). The 3D models are capable of accurately reproducing the surface structure of the observed solar granulation with their strongly inhomogeneous surface intensities (Stein & Nordlund 1998). The velocity fields predicted by the 3D simulations lead to a close match with both the observed spectral line bisectors and the broadening of their profiles in the atmospheres of different stars (e.g., Ramírez et al. 2009; Asplund et al. 2000; Dravins & Nordlund 1990; Allende Prieto et al. 2002). Recently, impressive agreement between a new synthetic 3D model and solar observations has been found in a detailed comparison of spectral line shifts, equivalent widths and center-to-limb variations for normalized line profiles (Pereira et al. 2009b,a). In essentially all cases, this 3D model reproduced the observations with an accuracy that is comparable to the semi-empirical model of Holweger & Mueller (1974), which is traditionally used in spectroscopy of the solar photosphere.

The accuracy of the treatment of radiation in 3D, however, is still strongly limited by the available computational power. Radiative transfer easily becomes the most computationally expensive part of a simulation, since the equations must be solved for a considerably larger set of transport directions compared to hydrodynamics, and non-grey opacities must be accounted for in realistic simulations. Most of the currently existing 3D radiative (M)HD codes therefore assume LTE and capture the atmospheric height dependence of continuum and line opacities using the opacity binning method (e.g., Ludwig 1992; Nordlund 1982): the problem of computing the monochromatic radiation field for a larger number of wavelengths is reduced to the numerical solution of the radiative transfer equation for typically 5 opacity bins. Skartlien (2000) extended the opacity binning method to include coherent scattering, and showed its importance in the solar chromosphere using a 3D radiative transfer solver for parallel shared-memory architectures.

Modern large-scale computer clusters use distributed memory architectures to handle the growing complexity of scientific simulations, allowing, e.g., self-consistent MHD models of the solar chromosphere, transition region and corona (Hansteen et al. 2007; Hansteen 2004) or detailed hydrodynamical models of giant stars (Collet et al. 2007). We present a new fully MPI-parallelized radiative transfer solver with coherent scattering for the new BIFROST code for time-dependent 3D MHD simulations of cool stellar atmospheres (Gudiksen et al., in preparation).

In Sects. 2 and 3, we discuss the physics of the radiative transfer model and its implementation in the MHD code. Section 4 describes the most important continuous and line opacity sources that we include in our simulations. Section 5 describes the application of the BIFROST code to model the atmosphere of a solar-type star using radiative transfer calculations with scattering, and discusses the effects on the temperature structure.

2 Radiative transfer with scattering and the radiative flux divergence

2.1 The radiative transfer equation

Hydrodynamical simulations of cool stellar atmospheres need to cover several pressure scale heights above and below the optical surface to minimize the effect of the boundaries on the granulation flow. The exponential density stratification causes the optical depth of the plasma to span about 15 orders of magnitude from the highest to the lowest layers of the simulation. The radiative transfer problem must therefore be solved in very different physical environments: in the extremely optically thick diffusion region at the bottom of the simulation box, all photons are thermalized. At the top, the atmosphere is mostly optically thin and mainly photons in the strongest lines interact with the gas. For the bulk of the photons, the transition between these two domains is rapid; it is confined to a thin layer which appears corrugated due to the different geometrical depth variation of opacities in upflows and downflows (Stein & Nordlund 1998).

Radiative transfer is, in general, a time-dependent process, which needs to be treated simultaneously with the hydrodynamics. However, the timescale of photon propagation over a mean free path length, $t_{\lambda}=(c\chi_{\lambda})^{-1}$ , where $\chi_{\lambda}$ is the monochromatic opacity and c is the speed of light, is orders of magnitude shorter than any hydrodynamical timescale. Radiative transfer therefore decouples from the hydrodynamics and is well approximated by a time-independent problem, described by a radiative transfer equation for the monochromatic specific intensity $I_\lambda(\vec{x},\vec{\hat{n}})$ in direction $\vec{\hat{n}}$ :

$\begin{displaymath}% \vec{\hat{n}}\cdot\vec{\nabla}I_{\lambda}(\vec{x},\vec{\hat... ...da}(\vec{x},\vec{\hat{n}})+j_{\lambda}(\vec{x},\vec{\hat{n}}), \end{displaymath}$

(1)

where $j_{\lambda}$ denotes the local emission at wavelength $\lambda$ (see, e.g., Rutten 2003; Mihalas 1978). The left-hand side of Eq. (1) is defined in the rest frame of the model atmosphere. The source and sink terms of the right-hand side are naturally described in the co-moving frame of the flowing gas. The consequent Doppler shifts are difficult to treat in 3D time-dependent simulations due to restrictions in computational power, requiring us to compute wavelength-integrated quantities in the opacity binning approximation (see below). We therefore assume a static medium, neglecting possible influences of the velocity field.

The extinction of photons is described, as customary, through the absorption coefficient $\kappa_{\lambda}$ and the scattering coefficient $\sigma_{\lambda}$ , which combine to the gas opacity,

$\begin{displaymath}% \chi_{\lambda}(\vec{x})=\kappa_{\lambda}(\vec{x})+\sigma_{\lambda}(\vec{x}), \end{displaymath}$

(2)

and give rise to the definition of the photon destruction probability

$\begin{displaymath}% \epsilon_{\lambda}(\vec{x})=\frac{\kappa_{\lambda}(\vec{x})}{\kappa_{\lambda}(\vec{x})+\sigma_{\lambda}(\vec{x})}\cdot \end{displaymath}$

(3)

Recasting the optical path ${\rm d}s=\vec{\hat{n}}\cdot {\rm d}\vec{x}$ along a ray in direction $\vec{\hat{n}}$ into the optical depth ${\rm d}\tau_{\lambda}=\chi_{\lambda}{\rm d}s$ along that direction, gives Eq. (1) the form

$\begin{displaymath}% \frac{{\rm d}I_{\lambda}}{{\rm d}\tau_{\lambda}}(\tau_{\lambda})=-I_{\lambda}(\tau_{\lambda})+S_{\lambda}(\tau_{\lambda}), \end{displaymath}$

(4)

with the source function $S_{\lambda}\equiv j_{\lambda}/\chi_{\lambda}$ . For the numerical computation, we employ the formal solution

$\begin{displaymath}% I_{\lambda}(\tau_{\lambda})=I_{\lambda}(\tau_{{\rm u},\lamb... ...au_{\lambda}}S_{\lambda}(t){\rm e}^{t-\tau_{\lambda}}{\rm d}t, \end{displaymath}$

(5)

where $I_{\lambda}(\tau_{{\rm u},\lambda})$ is the incident intensity at the upwind end of the ray at optical depth $\tau_{{\rm u},\lambda}<\tau_{\lambda}$ .

The source function $S_{\lambda}$ at optical depth $\tau_{\lambda}$ in direction $\vec{\hat{n}}$ includes local thermal radiation from the gas and coherent scattering of photons:

$\begin{displaymath}% S_{\lambda}\!=\!\frac{\sigma_{\lambda}}{4\pi\chi_{\lambda}}... ...lon_{\lambda}\right)J_{\lambda}+\epsilon_{\lambda}B_{\lambda}, \end{displaymath}$

(6)

where scattered radiation from direction $\vec{\hat{n}'}$ contributes with weight $\phi$ in the integral over the unit sphere S², $B_{\lambda}$ denotes the Planck function, and $J_{\lambda}$ is the mean intensity. The second equality holds for isotropic angular redistribution of radiation ( $\phi=1$ ). For $\epsilon_{\lambda}<1$ , the source function depends on $J_{\lambda}$ and thus, through the non-locality of radiative transfer, on radiation processes in the entire simulation domain. This turns Eq. (4) from an ordinary differential equation into an integro-differential equation.

Current limitations of available computing resources require the assumption of isotropic coherent scattering. Continuum processes in cool stellar atmospheres and very strong lines fulfill this restriction in very good or reasonable approximation, respectively, due to their weak wavelength dependence. Intermediate and weak lines are more accurately treated in complete spectral redistribution.

2.2 The radiative flux divergence and the wavelength integral

Absorption and thermal emission of radiation couples the stellar plasma with the radiation field through the transfer of heat. Photon energies in cool stars are too small to exert a significant force on the fluid compared to the gas pressure and gravity; the coupling is therefore sufficiently described by adding a radiative heating term $Q_{{\rm rad}}$ to the energy equation.

Evaluating the first moment of Eq. (1) and using the above definitions yields

$\begin{displaymath}% -\vec{\nabla}\cdot \vec{F}_{\lambda}=4\pi\chi_{\lambda}\lef... ...n_{\lambda}\chi_{\lambda}\left(J_{\lambda}-B_{\lambda}\right), \end{displaymath}$

(7)

where $\vec{\nabla}\cdot \vec{F}_{\lambda}$ is the local monochromatic radiative flux divergence. The second equality holds in the case of the coherent scattering source function (Eq. (6)). The scattering term does not contribute to heat exchange by definition, reducing radiative heating and cooling by a factor of $\epsilon _{\lambda }$ compared to the case where $S_{\lambda}=B_{\lambda}$ .

Integrating the monochromatic flux divergence in Eq. (7) over the whole wavelength spectrum of the star yields the local heating rate $Q_{{\rm rad}}$ :

$\begin{displaymath}% Q_{{\rm rad}}=-\int_{0}^{\infty}\vec{\nabla}\cdot \vec{F}_{\lambda}{\rm d}\lambda. \end{displaymath}$

(8)

In the optically thick regime, where radiative transfer is diffusive, this integral may be simplified with good accuracy by assuming the Rosseland mean opacity in the so-called gray approximation. However, gray opacities are not sufficient for a realistic treatment of the height-dependent line-blanketing above the surface, where the atmospheric structure is very sensitive to the radiation field. Atomic and molecular lines are important opacity sources in this region, changing the radiative heating and cooling compared to the simplified case of a gray atmosphere (see, e.g., Vögler et al. 2004, for a detailed discussion). The current version of the MARCS 1D atmosphere code uses the opacity sampling technique (Peytremann 1974), which approximates the spectrum through statistical sampling at ${\sim}100~000$ wavelength points. This resolution is sufficient to capture continuum absorption and line-blanketing without bias, at least in the lower parts of the atmosphere where the spectral distribution of absorbers is sufficiently widespread. Stellar atmosphere models in 3D do not allow for such a detailed treatment yet, since a single formal solution is many orders of magnitude more expensive to compute: the radiative transfer equation in our radiation-hydrodynamical model (Sect. 5.1) is solved for 240 $\times$ 240 columns and 24 transport directions, which is equivalent to ${\sim}10^{5}$ 1D calculations for each time step.

Nordlund (1982), Ludwig (1992) and Skartlien (2000) have described opacity binning techniques, where wavelength integration is performed over subsets of the spectral range before the solution of the radiative transfer equation, and the radiation field is computed for only a few mean opacities instead of the full spectrum. We will give a brief description of the technique in the following; see Skartlien (2000) for a more detailed discussion.

Integrating the radiative transfer equation (Eq. (1)) over wavelength leads to the definition of a mean opacity, mean scattering coefficient and mean absorption coefficient:

$\displaystyle \chi^{{I}} = \frac{\int\chi_{\lambda}I_{\lambda}{\rm d}\lambda}{\int I_{\lambda}{\rm d}\lambda}$	(9)
$\displaystyle \sigma^{{J}} = \frac{\int\sigma_{\lambda}J_{\lambda}{\rm d}\lambda}{\int J_{\lambda}{\rm d}\lambda}$	(10)
$\displaystyle \kappa^{{B}} = \frac{\int\kappa_{\lambda}B_{\lambda}{\rm d}\lambda}{\int B_{\lambda}{\rm d}\lambda}\cdot$	(11)

The intensity-weighted mean opacity $\chi^{{I}}$ and the mean-intensity-weighted mean scattering coefficient $\sigma^{{J}}$ depend on the unknown radiation field $I_{\lambda}$ and its angular average $J_{\lambda}$ , which must be estimated: we use 1D radiative transfer calculations on the mean stratification of the atmosphere (see Sect. 5.2), which yield approximations for $\chi^{{I}}\approx\chi^{J,{\rm 1D}}$ and $\sigma^{{J}}\approx\sigma^{J,{\rm 1D}}$ .

These three mean coefficients represent absorption, scattering and thermal emission of photons with good accuracy where the stellar atmosphere is optically thin across the spectrum. However, $\chi^{J,{\rm 1D}}$ does not ensure a correct total radiative energy flux at optical depths $\tau\gg1$ where radiative transfer is diffusive. It needs to be replaced by the Rosseland mean opacity, defined as the weighted harmonic mean

$\begin{displaymath}% \chi^{{R}}=\frac{\int\left({\rm d}B_{\lambda}/{\rm d}s\righ... ...t)\left({\rm d}B_{\lambda}/{\rm d}s\right){\rm d}\lambda}\cdot \end{displaymath}$

(12)

We consequently use a $\tau$ -weighted sum of the two quantities $\chi^{J,{\rm 1D}}$ and $\chi^{{R}}$ . The geometrical depth of the transition between the two regimes near $\tau\approx1$ varies quickly with wavelength where spectral lines are present, and it is not sufficient to consider only a single pair of mean opacities $\chi^{J,{\rm 1D}}$ and $\chi^{{R}}$ . The opacity binning method therefore defines several opacity groups, where each member reaches unit optical depth ( $\tau_{\lambda}=1$ ) at a similar geometrical depth. The integrals in Eqs. (9)-(12) are then evaluated only for a set of member wavelengths $\{\lambda_{i}\}$ in each bin i, which does not have to be continuous.

Depending on the height range of the stellar atmosphere model and the wavelength selection method, it turns out that about 5 such opacity bins are enough to capture the essence of the line-blanketing and continuum opacity and to obtain a realistic temperature structure (Vögler et al. 2004). More recent atmosphere models have been extended to 12 bins (Caffau et al. 2008). For the simulations presented in this work, we compute radiative transfer with 12 bins, where wavelengths are sorted not only by the geometrical height of the monochromatic optical surface, but also by wavelength, separating opacities in the UV, visual and infrared bands (Trampedach et al., in preparation).

It is difficult to assess the quality of the opacity binning method in realistic 3D simulations: deviations of the resulting radiative heating rates $Q_{{\rm rad}}$ from an accurate monochromatic solution have a height-dependent impact on the temperature structure (see Sect. 5), making the long-term behavior of the simulation hard to predict. The agreement of 3D model atmospheres with various observational tests indicates that opacity binning still yields a reasonable estimate for the line-blanketing.

3 The numerical implementation

The large variety of radiative transfer models for astrophysical problems inspired the development of very different analytical and numerical methods to obtain the radiation field (see, e.g., Wehrse & Kalkofen 2006, for an overview). For our given problem of computing radiative heating rates as the flux divergence $-\vec{\nabla}\cdot \vec{F}$ of a time-independent radiation field in 3D, the direct solution of Eq. (4) yields accurate results with efficient numerical schemes.

Characteristics methods, which solve the transfer problem along a discrete set of light rays to capture the anisotropy of the radiation field in the optically thin atmosphere, are a popular choice in stellar atmosphere models. Nordlund (1982) and Skartlien (2000) use Feautrier-type differential radiative transfer solvers (Feautrier 1964) for solving Eq. (4) on long characteristics. They span across the entire simulation domain, which is an obstacle for a domain-decomposed parallelization of the MHD code (see Sect. 3.2 below). Bruls et al. (1999), Vögler et al. (2005) and Muthsam et al. (2010) employ the short characteristics method (Mihalas et al. 1978; Olson & Kunasz 1987; Kunasz & Auer 1988), where the radiative transfer equation is solved on characteristics which only extend to the adjacent upwind and downwind grid layers. This method is required by our choice of iteration technique for an efficient solution of the scattering problem.

3.1 Short characteristics

The short characteristics method employs the formal solution (Eq. (5)) of the monochromatic radiative transfer equation (Eq. (4)) to compute the radiation field at the center of a three-point ray for a known source function $S_{\lambda}$ . The discretization is performed by interpolating the source function for a given wavelength $\lambda$ (or bin number) along the ray using a second-order Bézier curve (see, e.g., the discussion in Auer 2003)

$\begin{displaymath}% S(t)=(1-t)^2S_{{\rm u}}+t^2S_{0}+2t(1-t)S_{{\rm c}}, \end{displaymath}$

(13)

where $S_{{\rm u}}$ and S₀ are the upwind and local source functions and $t=(\tau-\tau_{{\rm u}})/(\tau_{0}-\tau_{{\rm u}})$ is the curve parameter. $S_{{\rm c}}$ is a control point, which shapes the interpolating curve by restricting it to the convex hull laid out by $S_{{\rm u}}$ , $S_{{\rm c}}$ and S₀. This characteristic of Bézier curves may be exploited to detect and suppress overshoots, which destabilize the numerical solution at places in the atmosphere where strong opacity and temperature gradients occur. Inserting Eq. (13) into the formal solution (Eq. (5)), evaluating the integral and reordering the terms yields the discretized expression

$\begin{displaymath}% I(\tau)=I(\tau_{{\rm u}}){\rm e}^{-(\tau-\tau_{{\rm u}})}+\Psi_{{\rm u}}S_{{\rm u}}+\Psi_{0}S_{0}+\Psi_{{\rm d}}S_{{\rm d}}. \end{displaymath}$

(14)

The shape of the three interpolation coefficients $\Psi_{{\rm u}}$ , $\Psi_{0}$ and $\Psi_{{\rm d}}$ for the upwind, center and downwind source functions depends on the control point $S_{{\rm c}}$ . Choosing

$\begin{displaymath}% S_{{\rm c}}=S_{0}-\frac{\tau_{0}-\tau_{{\rm u}}}{2}S_{0}', \end{displaymath}$

(15)

where S₀' is the centered derivative on the three-point stencil $(S_{{\rm u}},S_{0},S_{{\rm d}})$ , yields second-order interpolation. It is used where no overshoots happen and correctly reproduces the diffusion approximation at optical depths $\tau \ga 30$ (see Appendix A for the detailed shape of the $\Psi$ coefficients). In the optically thin atmosphere where $\tau \la 10^{-3}$ , a second-order expansion of the $(1-{\rm e}^{-\tau})$ terms in the $\Psi$ constants stabilizes the solver, which may therefore be implemented with single precision floating point numerics throughout the simulation domain. Optical depths $\Delta\tau$ along the characteristics are similarly computed using the Bézier interpolation technique.

The mean intensities J and the components of the flux vector $\vec{F}$ are computed by approximating the zeroth and first moment integrals by a quadrature sum over selected ray angles (``method of discrete ordinates''),

$\begin{displaymath}% J\approx \frac{1}{4\pi} \sum_{i}w_{i}I(\vec{\hat{n}}_{i}); ... ...I(\vec{\hat{n}}_{i})(\vec{\hat{n}}_{i}\cdot\vec{\hat{n}}_{j}), \end{displaymath}$

(16)

where w_i is the weight of direction $\vec{\hat{n}}_{i}$ . The best choice of quadrature depends on the expected anisotropy of the radiation field and on the quantity that needs to be computed. In our case, the components of the flux vector $\vec{F}$ need to be calculated explicitly, requiring the quadrature be invariant to rotation by $\pi/2$ to avoid directional bias. Carlson's A4 quadrature (Carlson 1963) with 3 rays per octant is an appropriate choice and represents the anisotropy with good accuracy.

Short characteristics require knowledge of the upwind intensities $I(\tau_{{\rm u}})$ for each ray direction $\vec{\hat{n}}$ , on which the sweep direction for a formal solution therefore depends. Interpolation yields all such quantities (Sect. 3.3). Shallow rays, that fail to hit the upwind layer within the grid cells, need to be extended and may cross several cells, possibly across subdomain boundaries. For the first formal solution of a simulation run, a Feautrier-type long characteristics solver delivers boundary intensity estimates; intensities from the previous iteration in the neighbor subdomains are used for all subsequent computations. Once $I(\tau)$ is known along two edges of the current layer, the remaining unknown intensities may be computed away from the boundary through vertical interpolation between the upwind layer and the current layer. It is worth noting that some long characteristics codes turn transport directions around the vertical axis with every time step to avoid numerical artefacts stemming from a fixed set of discrete ordinates. Such an effect is not observed in our short characteristics implementation. Moreover, the anisotropy of the radiation field slows down convergence of an iterative solution in optically thin parts when transport directions are turned between time steps, since the stored boundary intensities come from the previous solution (see Sect. 3.2 for further details). All ray directions are therefore kept fixed.

The discretized formal solution (Eq. (14)) in the simulation domain and averaging of the radiation field over solid angle will be abbreviated in the following using the $\Lambda$ operator, which is commonly defined through

$\begin{displaymath}% J=\Lambda S. \end{displaymath}$

(17)

$\Lambda$ is a linear matrix operator on the source functions S which represents the numerical algorithm used to compute the radiation field in the code.

3.2 The Gauss-Seidel scheme and MPI parallelization

As noted in Sect. 2.1, the coherent scattering term turns the transfer equation into an integro-differential equation for the specific intensity I. Using the $\Lambda$ operator defined in Eq. (17), the problem may be rewritten into the matrix equation

$\begin{displaymath}% \left[\mathbbm{1}-(1-\epsilon)\Lambda\right]S=\epsilon B, \end{displaymath}$

(18)

with the identity matrix $\mathbbm{1}$ . The expression represents a very large system of linear equations. Its direct solution in 3D through inversion of the operator on the left-hand side is far too complex and numerically unstable in some cases, to be of practical use. Most solvers therefore apply an iteration scheme, the choice of which depends on the structure of the $\Lambda$ operator matrix. Approximate Lambda Iteration (ALI, sometimes also called Accelerated Lambda Iteration, Cannon 1973) is a popular method to obtain a good approximation of the radiation field with fast convergence. J is computed through a formal solution and used to correct the source function. Rather than just inserting J in S, which leads to very slow convergence (or no convergence at all), a largely simplified approximate operator $\Lambda^{\ast}$ is used to compute correction values $\Delta S$ at low cost, speeding up convergence tremendously.

We employ the Gauss-Seidel scheme (Trujillo Bueno & Fabiani Bendicho 1995), an ALI method that combines the formal solution and correction steps. It mimics a tridiagonal $\Lambda^{\ast}$ operator, but the scheme does not require the expensive construction of the matrix. Source function corrections at the grid point i are obtained during a solver sweep from the expression

$\begin{displaymath}% \Delta S_{i}=\frac{(1-\epsilon_{i})J_{i}^{{\rm old/new}}+\e... ...}B_{i}-S_{i}^{{\rm old}}}{1-(1-\epsilon_{i})\Lambda_{ii}}\cdot \end{displaymath}$

(19)

$J_{i}^{{\rm old/new}}$ is the radiation field that includes the corrections in the upwind part of the simulation domain, which have already been computed during the current sweep. The dependence of $\Delta S_{i}$ on J_i in each layer for immediate correction of S_i during the sweep requires employing the short characteristics method. The denominator contains the diagonal element $\Lambda_{ii}$ of the $\Lambda$ operator, which may be computed using Eq. (14) and therefore reduces to a sum of $\Psi$ constants. Source function corrections may be applied during both upsweeps and downsweeps for faster convergence.

We tested our radiative transfer code by comparing the numerical results with an analytical solution for the case of an isothermal 1D atmosphere with constant photon destruction probability $\epsilon$ (see the discussion in Trujillo Bueno & Fabiani Bendicho 1995) and found very good agreement.

The radiation solver is parallelized using spatial domain decomposition and communication with the MPI library, adopting the virtual topology given by the MHD solver of the BIFROST code. The grid is decomposed into cuboid subdomains, allowing an arbitrary number of divisions on all three spatial axes. While this parallelization lends itself to a mixed initial and boundary value problem found in computational hydrodynamics, it is harder to apply in an efficient way to the pure boundary value problem of time-independent radiative transfer. Concurrent computation of spectral subdomains (or opacity bins) would provide a higher degree of parallelism considering the non-local dependencies in a monochromatic formal solution of our given coherent scattering problem, but such an approach would cause severe load balancing issues and suffer from node memory limitations when applying the code to very large simulations. Spatial domain decomposition may still be combined with spectral domain decomposition if radiative transfer needs to be solved for a large number of wavelengths.

Heinemann et al. (2006) have presented a domain-decomposed method based on a variant of the formal solution (Eq. (5)) on long characteristics. The solver bypasses the problem of missing incident intensities at subdomain boundaries by splitting the local and boundary contributions. While their approach efficiently solves the radiative transfer equation without scattering, the long characteristics solver would have to be combined with a different ALI scheme than Gauss-Seidel. An approximate $\Lambda^{\ast}$ operator needs a certain bandwidth around its matrix diagonal to achieve good convergence (see, e.g., the discussion in Hauschildt & Baron 2006). It is therefore more expensive to construct and invert than the diagonal operator used for the Gauss-Seidel scheme.

Our code iterates the solution, starting with the source function and subdomain boundary intensities from the previous hydrodynamical time step, until the maximum relative source function correction in the domain after the nth iteration is smaller than a preset threshold C:

$\begin{displaymath}% \max\left(\frac{\left\vert S_{i}^{n}-S_{i}^{n-1}\right\vert}{S_{i}^{n-1}}\right)\leqslant C. \end{displaymath}$

(20)

When scattering is not included, the maximum relative change of mean intensities at the boundary is used instead to test the convergence of the radiation field:

$\begin{displaymath}% \max\left(\frac{\left\vert J_{i}^{n}-J_{i}^{n-1}\right\vert}{J_{i}^{n-1}}\right)\leqslant C. \end{displaymath}$

(21)

If too few iterations are performed, the subdomain boundaries produce artifacts in the upper parts of the atmosphere, where photon mean free paths are comparable to or larger than the subdomain size. In practice, it turns out that a threshold of $C\sim10^{-3}$ yields good results in either case.

$\begin{figure} \par\includegraphics[width=7.4cm,clip]{14210f1a.eps}\hspace*{7mm} \includegraphics[width=7.55cm,clip]{14210f1b.eps} \end{figure}$

Figure 1:

Horizontal mean photon destruction probability $\epsilon$ for three bins representing continuum and weak lines (1), intermediate lines (2) and strong lines (3) in the UV, plotted as a function of geometrical depth (left) and optical depth in the respective bin (right, the average is taken over surfaces with the same vertical optical depth). The dotted line marks the stellar surface (left) and unit optical depth in each bin (right).

Open with DEXTER

$\begin{figure} \par\includegraphics[width=11.3cm,clip]{14210f2.eps} \end{figure}$

Figure 2:

Convergence of the source function for bins 1-3 (see Fig. 1) during a simulation run without domain decomposition (left column), with $2\times 2\times 2$ decomposition (center column) and $3\times 3\times 3$ decomposition (right column), and with a time step of $\Delta t=0.03$ s (upper row), $\Delta t=0.04$ s (center row) and $\Delta t=0.08$ s ( lower row). Line styles represent the same bins as in Fig. 1. Thin lines: relative source function correction $\Delta S$ after n iterations with respect to S^n-1 from the previous iteration n-1. Thick lines: relative source function correction $\Delta S$ with respect to the ``true'' solution $S^{\infty }$ . Dotted lines mark the threshold C beneath which convergence is assumed (see text).

Open with DEXTER

The convergence speed of an iterative method depends on the spectral radius $\rho$ of the operator with which corrections are computed, as the error of the solution after n iterations decreases with $\rho^{n}$ . The spectral radius approaches $\rho\approx 1-\epsilon$ for optically thick scattering media (see, e.g., the discussion in Trujillo Bueno & Fabiani Bendicho 1995). Strong scattering at high optical depths therefore leads to very poor convergence rates of the Gauss-Seidel solver, requiring hundreds of iterations in extreme situations. However, this difficulty is mostly alleviated by using the source function solution from the previous time step and the slow evolution of the plasma flow between consecutive time steps, so that the code ideally needs to fully converge the solution only once at the beginning. Domain decomposition additionally slows down convergence if the photon mean free paths cross subdomain boundaries, which is the case at continuum wavelengths in the thin atmosphere, since the subdomain boundary intensities are not initially known. Storing intensities from the previous time step again largely circumvents this problem, and the actual number of iterations per time step that is required during a simulation run depends on how fast the atmosphere evolves.

We therefore test the convergence of the solution for arbitrary time steps of our solar-type simulation using 12 opacity bins with continuum and line scattering (see Sect. 5), following a similar discussion in Skartlien (2000). The tests were run at half resolution on all axes to facilitate computation on a single core, which yields slightly faster convergence. Since the true solution S of our radiative transfer problem is unknown, we compare the approximate solution after n iterations, Sⁿ, with an approximate solution $S^{\infty }$ which we obtained after additional iterations with a lower convergence threshold of $C\sim10^{-4}$ , assuming $S^{\infty}\approx S$ with good accuracy.

We use three representative opacity bins, which cover weak, intermediate and strong opacities in the UV, with different depth-dependence of the scattering strengths. The remaining nine bins at longer wavelengths behave in a similar way. Figure 1 shows horizontal averages of the photon destruction probabilities $\epsilon$ for each bin in an arbitrary snapshot of our photospheric simulation: averages over layers with the same geometrical depth are plotted in the left panel, averages over surfaces with the same vertical optical depth are plotted in the right panel.

Figure 2 compares the convergence speed for the radiative transfer solution of the sample bins with and without domain decomposition, and with different time step lengths. Thick lines represent convergence relative to the true solution $S^{\infty }$ for each bin, thin lines show the convergence relative to the solution obtained in the previous iteration, which we use as the convergence criterion. In normal operation, the solver would stop as soon as the thin line of the currently computed opacity bin has crossed the dotted horizontal line. We caution that the number of iterations needed for a solution also depends mildly on the time stepping algorithm, since the choice of method affects the deviation of stored boundary intensities and source functions between substeps of the time integration. We therefore only analyze the behavior for the first extrapolation step of a 3rd order Runge-Kutta time stepper.

The poorer convergence speed caused by scattering at high optical depths in bin 3 is evident in all plots (thick dot-dashed line), compared to the situation in bin 1, where the photon destruction probability is larger. The small optical path lengths of bin 3 reduce the impact of domain decomposition, since the radiation field is essentially local in most parts of the simulation box. Contrary to that, bin 1 suffers most strongly from slower convergence with increasing number of subdomain divisions, as well as from some flip-flopping of $\Delta S$ . The latter is caused by high-order interpolation (see Sect. 3.3) and disappears when the solver is set to linear interpolation. High order interpolation of upwind intensities widens the domain of dependence of the short characteristics, and the effect is amplified where large path lengths in the optically thin regime cross subdomain boundaries.

Domain decomposition mildly slows down convergence, and the accuracy of the solution in bin 3 slightly deteriorates for a larger number of subdomains. Longer time steps have the same effect on that bin, causing slower convergence towards $S^{\infty }$ than indicated by the relative corrections with respect to S^n-1 (thick and thin dot-dashed lines in Fig. 2). The method devised by Skartlien (2000) exhibits similar behavior for bins with strong scattering lines.

The effect of such inaccuracies in the numerical solution of S and J on the energy transfer between the radiation field and the gas are nevertheless small or even vanish in some regions: radiative heating is reduced in the atmosphere where coherent scattering is important (see Eq. (7) and Fig. 1). Coherent scattering also effectively damps the impact of any remaining discontinuities in the radiation field across subdomain boundaries on the flux divergence in the optically thin atmosphere, so that no visible artifacts from the domain decomposition remain in the gas temperatures.

Compared to the solver proposed by Heinemann et al. (2006), it is clear that our method is not optimal for the case without scattering, since several computationally expensive formal solution and communication steps are required to obtain a radiation field that is consistent in the whole domain. It offers good performance when scattering is included, which is not considered in their method.

3.3 Interpolation and grid refinement

$\begin{figure} \par\includegraphics[width=6.6cm,clip]{14210f3a.eps}\hspace{7mm}... ...7mm} \includegraphics[width=6.6cm,clip]{14210f3d.eps} \vspace{5mm} \end{figure}$	Figure 3: Numerical diffusion of a searchlight beam with rectangular cross-section using linear (upper right), local cubic (lower left) and local cubic monotonic (lower right) interpolation, compared to the exact solution (upper left).
Open with DEXTER

At every time step, the hydrodynamical solver updates mass densities and internal energy densities. These quantities are used to look up tabulated opacities, bin-integrated Planck functions and photon destruction probabilities at every grid point. In general, the characteristics grid needed to represent the anisotropy of the radiation field does not coincide with the hydrodynamical mesh, requiring the interpolation of $\chi$ , S and the upwind intensities $I_{{\rm u}}$ during the formal solution.

The accuracy of this interpolation strongly influences the overall accuracy of the solver, and there is a large choice of possible methods (see, e.g., the discussion in Auer 2003). Linear interpolation is fast and avoids instabilities produced by interpolation overshoots, but yields poor estimates where the radiation field is not well-resolved, e.g. between granules and intergranular lanes at the optical surface. It also amplifies the numerical diffusion effect of short characteristics, where lateral diffusion artificially transports radiation away from the beam.

To illustrate this behavior, we repeat the searchlight test of Kunasz & Auer (1988), where a rectangular light beam is cast through an empty 3D box with a 100³ mesh and zero opacity. Any diffusion of radiation away from the beam results in a broadening of the beam profile at the surface and can only stem from the interpolation of unattenuated upwind intensities. The light source illuminates the bottom of the 3D box, where it initially covers an area of 30² mesh points; it is slanted with an angle of $\theta=28.1^{\circ}$ off the vertical and an azimuth of $\phi=45.0^{\circ}$ . The upper left panel in Fig. 3 shows the beam profile at the top of the 3D box expected from an exact solution of the unattenuated transfer problem through vacuum; note that the finite resolution of the surface in the plot leads to a slightly widened profile. The upper right panel shows the broadening of the beam profile caused by 100 consecutive linear interpolations applied for the numerical transfer through the box. Although the area-integrated intensity is conserved with good accuracy, limited by the machine precision, the beam is visibly widened through numerical diffusion. The lower left panel in Fig. 3 shows the result when using local cubic interpolation for the transport problem. The broadening is reduced, but the overshooting cubic polynomials produce ringing and negative intensities. We therefore use the local cubic monotonic interpolation scheme of Fritsch & Butland (1984), which effectively suppresses overshoots by using weighted harmonic mean derivatives, in consecutive 1D-1D interpolation on horizontal planes, and local quadratic interpolation on vertical cell walls (see Appendix B for further details). The lower right panel in Fig. 3 shows the result from the searchlight test, where the beam profile is conserved to a satisfactory degree. Numerical diffusion is reduced and reaches a level which renders the computed flux divergences comparable to those obtained with long characteristics codes: although upwind intensities do not need interpolation along the beam, diffusion affects the local flux divergences when transfered from the slanted long characteristics grid back to the hydrodynamical grid.

The basic mesh on which radiative transfer is computed is imposed by the MHD solver. This is usually not critical in the optically thin upper atmosphere and the optically thick interior, where radiative transfer is simple and may even be over-resolved. The opposite is the case in the transition region around the optical surface, where opacities drop rapidly due to their strong temperature dependence and cause a runaway cooling effect (Stein & Nordlund 1998). For a solar simulation, 1D tests performed by Nordlund & Stein (1991) indicate that a vertical spacing of ${\la}10$ km is desirable at this atmospheric height. Using a non-linear vertical grid with the finest resolution around the surface, this is easily achievable in 3D for modern MPI-based domain decomposed radiative hydrodynamics codes. However, for large coronal simulations or in the case of giant stars, where the spatial scales needed to resolve hydrodynamics and radiation transport exhibit much larger disparity than in the Sun, finding the optimal grid leads to a conflict. Besides the larger simulation size, too small length intervals $\Delta x$ drastically increase the stiffness of the hydrodynamical equations, where the stability-limited time steps of the transport and diffusion terms scale with $\Delta x$ and $\Delta x^2$ , respectively, and quickly become exceedingly small. In extreme cases, both effects may increase computation times of a model beyond tractability.

A fully adaptive mesh for computing radiative transfer would yield optimal results without affecting the stiffness of the equations, but is difficult to realize in a characteristics method. We achieve partial adaptivity by inserting horizontal layers in the hydrodynamical mesh for the radiative transfer computation, reducing optical path lengths without reducing the time steps. The refinement is based on the maximum vertical gradient of the Rosseland mean opacity in each layer and reassessed in regular intervals. While inserting additional layers slows down convergence of the Gauss-Seidel method (see Sect. 3.2), this is again overcome by storing the source function from the previous time step.

3.4 Numerical flux divergences

Having established a method for numerically computing radiative transfer with coherent scattering in a decomposed simulation domain, we now need to obtain flux divergences $\vec{\nabla}\cdot\vec{F}$ , a derivative of the radiation field.

The right-hand side of Eq. (7) involves only local quantities that are defined on the cell centers of the hydrodynamical mesh, where $Q_{{\rm rad}}$ is eventually needed, and therefore seems a natural choice. The expression $\chi(J-S)$ is numerically stable in the optically thin regime, where round-off errors of a possibly vanishing difference between J and S are attenuated by the exponential outward decrease of the opacity $\chi$ . At the same time, $\chi$ amplifies round-off errors of (J-S) beneath the optical surface, where the radiation field thermalizes ( $J\approx B$ , also in the scattering case since $\epsilon>0$ ): the flux divergence again vanishes, but the finite machine precision prevents complete cancellation of the terms.

It is possible to stabilize a short characteristics solver in the whole simulation domain by subtracting S₀ from the discretized formal solution (Eq. (14)), which yields the modified integration constant $\tilde{\Psi}_{0}=\Psi_{0}-1$ . Using this equation, one obtains a monochromatic $Q_{{\rm rad},\lambda}(\vec{x},\vec{\hat{n}})$ along each ray. We note, however, that this leads to a deviation between the radiative energy $\iint Q_{{\rm rad},\lambda}(\vec{x},\vec{\hat{n}}){\rm d}\Omega {\rm d}V$ that is emitted by the gas in the simulation volume V per time unit, and the outgoing radiative flux computed from the specific intensities at the surface: the expressions are not equivalent anymore in their discretized form, and numerical errors affect the two values in a different way.

The discretized flux divergence $\vec{\nabla}\cdot\vec{F}$ on the left-hand side of Eq. (7) using finite difference quotients is stable in the optically thick regime, but its accuracy deteriorates outward: round-off errors quickly become significant, as the internal energy per gas volume decreases exponentially (see also the discussion in Bruls et al. 1999).

Adopting the approach presented in Bruls et al. (1999) and Vögler et al. (2005), we combine both expressions through exponential bridging in each vertical column of the simulation domain as a function of bin optical depth to benefit from their respective advantages. We slightly reduce the transition range between the regimes by a squared exponent, resulting in the expression:

$\begin{displaymath}% Q_{{\rm rad}}={\rm e}^{-\left(\tau/\tau_{0}\right)^2}Q_{{\r... ...{-\left(\tau/\tau_{0}\right)^2}\right)Q_{{\rm rad}}^{\vec{F}}, \end{displaymath}$

(22)

where $\tau_{0}=0.1$ , $Q_{{\rm rad}}^J=4\pi\chi(J-S)$ and $Q_{{\rm rad}}^F=-\vec{\nabla}\cdot \vec{F}$ , representing the two sides of Eq. (7). The total radiative energy computed with this expression delivers a consistent surface flux, since $Q_{{\rm rad}}^F$ contributes most of the radiative heating.

Following Vögler et al. (2005), we compute radiative transfer on cell corners to improve the accuracy of $Q_{{\rm rad}}^{\vec{F}}$ . Radiative fluxes $\vec{F}$ are averaged over cell corners surrounding each face before computing difference quotients, while $Q_{{\rm rad}}^J$ is averaged over all eight cell corners surrounding each grid point. Both expressions use exactly the same stencil and exhibit very good agreement around the threshold optical depth $\tau_{0}$ in our solar-type simulation.

Flux divergences are computed only on the hydrodynamical grid. Additional layers that are possibly inserted by the radiative transfer solver just serve to stabilize the computation and may simply be omitted when computing $Q_{{\rm rad}}$ , since conservation of the radiative energy flux through the hydrodynamical cell surfaces must hold.

4 Absorption and scattering opacity sources in the Sun

$\begin{figure} \par\includegraphics[width=12cm,clip]{14210f4.eps} \end{figure}$

Figure 4:

Wavelength and depth dependence of the continuum photon destruction probabilities $\epsilon _{\lambda }^{{\rm c}}$ (upper panel), the van Regemorter (1962) line photon destruction probabilities $\epsilon _{\lambda }^{{\rm l}}$ (center panel) and the total photon destruction probabilities $\epsilon _{\lambda }$ (lower panel) for the mean solar-type stratification. The zero point on the depth axis marks the stellar surface.

Open with DEXTER

A complete description of radiative transfer in stellar atmospheres requires a detailed wavelength-resolved treatment of numerous radiative absorption and emission processes, collisions with neutral atoms, electrons and ions in the plasma, as well as an evaluation of the feedback of the radiation field on the level populations of the interacting particles. The complexity of the resulting problem vastly exceeds current computational resources. We therefore restrict all of the underlying thermodynamical plasma states to LTE, neglecting the effects of radiation on the excitation and ionization of atoms and photo-dissociation of molecules. The cross-sections and level populations needed for the absorption and scattering coefficients then depend only on the gas density $\rho$ and the temperature T. Microscopic plasma thermodynamics is treated with the Mihalas-Hummer-Däppen equation of state (EOS) for stellar envelopes (Hummer & Mihalas 1988; Mihalas et al. 1988; Däppen et al. 1988; Mihalas et al. 1990) and used in tabulated form. The solar chemical composition for the 15 elements included in the EOS and for the opacities is taken from the abundances of Asplund et al. (2005).

4.1 Continuum opacity

The most important continuous opacity sources are various transitions of hydrogen atoms, their ions and molecules. The H^- ionization opacity dominates the solar continuum around the optical surface in the visual band; the large temperature sensitivity of the weakly bound second electron in the hydrogen atom causes runaway radiation cooling and the strong temperature gradient found at the top of the granules in the Sun (Stein & Nordlund 1998). Most solar continuum photons originate from this very thin layer. Among many other processes, photoionization of metals contributes significantly to the continuous opacity at shorter wavelengths. Table D.1 gives an overview of all sources considered in this work; our data is mostly identical to those used in the latest MARCS models (see Table 1 in Gustafsson et al. 2008), but includes additional bound-free data from the Opacity Project and the Iron Project (see Trampedach et al., in prep., for further details). We also include opacities of the second ionization stage for many metals, allowing 3D models to extend deeper into the convection zone than their 1D counterparts, which is a requirement for correctly simulating surface granulation.

The upper panel in Fig. 4 shows the wavelength and depth dependence of the continuum photon destruction probabilities $\epsilon _{\lambda }^{{\rm c}}$ for the mean stratification of our 3D model, including all continuous absorption and scattering opacity sources considered here. Continuum scattering has a significant contribution mostly above the surface, photons thermalize beneath at almost all wavelengths. Note that the narrow features at the short-wavelength end are the scattering resonances of the Lyman series; Lyman lines are nevertheless treated as true absorbers if line scattering is not included in the simulations. The Rayleigh scattering tail of H I contributes mostly to the UV continuum opacity in the upper solar photosphere due to its comparatively small cross-section and strong wavelength dependence ( $\sigma_{\lambda}\sim\lambda^{-4}$ ). The importance of elastic scattering on neutral hydrogen is outweighed by thermalizing processes closer to the surface and at short wavelengths. Electron scattering is wavelength-independent in the spectral range considered here, and becomes significant in the upper photosphere, where metals are the most important electron donors. It is mostly notable red-ward of the 1.644 $\mu$ m edge of H^- bound-free, before H^- free-free absorption takes over. Rayleigh scattering on He I atoms only gives minor contributions to the UV continuum opacity in the upper photosphere. The scattering opacity of H₂ molecules is negligible. Rayleigh and electron scattering are treated as isotropic, neglecting their weak $(1+\cos^{2}\theta)$ anisotropy, where $\theta$ is the scattering angle away from the incident direction (see, e.g., Mihalas 1978).

Between 5000 Å and 1.644 $\mu$ m, the strong H^- bound-free absorption opacity thermalizes the photons. Its dominance slightly decreases in the cool outermost layers owing to the lack of free electrons to form the ion.

4.2 Line opacity

Spectral line absorption and scattering are important processes which dictate the near-radiative equilibrium found in the solar photosphere. The heating/cooling effect of this line-blanketing forces the flatness of the observed temperature gradient, balancing the adiabatic dynamical gradient; see the discussion in Sect. 5.4. Spectral lines are particularly significant opacity sources at short wavelengths where many radiative bound-bound transitions of metals lie.

We obtain line opacities from extensive opacity sampling tables provided by B. Plez (2008, priv. comm.) as part of the MARCS collaboration. The data are based on VALD with some modifications; see Gustafsson et al. (2008) for further details. The original line data combine scattering and absorption contributions in a total opacity, which is sampled with ${\sim}100~000$ wavelengths and tabulated for a range of temperatures and pressures. The tables assume Saha ionization equilibrium and Boltzmann level populations to obtain the absorber density fractions. Departures from LTE, e.g. through radiative ionization, are neglected.

Following Skartlien (2000), we estimate the importance of scattering in line transitions by computing a photon destruction probability $\epsilon _{\lambda }^{{\rm l}}$ for every line opacity sample, using the van Regemorter (1962) formula (see Appendix C). We assume all scattering atoms to be neutral, accounting for the large contribution of Fe I to the line-blanketing (Anderson 1989), and all transitions to be permitted, in which case the assumptions of the van Regemorter (1962) formula yield reasonable estimates. Only electrons are taken into account for collisional de-excitation. The estimated photon destruction probability $\epsilon _{\lambda }^{{\rm l}}$ is then a function of wavelength, temperature and electron pressure, and independent of the actual transition. It may therefore also be applied in cases where the line opacity sample includes several transitions (see the discussion in Appendix C). Line transitions are treated as independent two-level processes without taking the coupling of the respective level populations into account, which is a reasonable assumption for resonance lines.

The center and lower panels in Fig. 4 show the wavelength and depth dependence of the estimated $\epsilon _{\lambda }^{{\rm l}}$ of spectral lines and the total photon destruction probabilities $\epsilon _{\lambda }$ , including all considered continuous and line processes. It is clear that collisional de-excitation dominates beneath the surface and at the longest wavelengths. Resonant line scattering becomes important towards optical and shorter wavelengths at increasing depth.

With the exception of very strong lines, line scattering is generally not coherent due to the Doppler shifts in the moving gas, which are not accounted for in our calculations. The two-level approximation probably gives a reasonably realistic picture of strong permitted lines, but departures from the LTE populations of the atomic levels are still neglected. The important Fe I opacity deviates from the LTE estimate in higher layers (see Fig. 7 in Short & Hauschildt 2005), thereby affecting the overall magnitude of the line-blanketing in these regions. Moreover, the accuracy of the opacity sampling method itself deteriorates outwards, where fewer and fewer lines contribute to the opacity. The van Regemorter approximation assumes resonant line scattering and consequently produces poorer estimates for all non-resonant lines. In summary, we should expect to obtain an order-of-magnitude estimate for the effects of scattering on the atmospheric structure. A more detailed picture requires a full treatment of the departures from LTE level populations and velocity fields, which is still out of reach for time-dependent 3D simulations.

$\begin{figure} \par\includegraphics[width=8cm,clip]{14210f5.eps} \end{figure}$

Figure 5:

Horizontal average heating rate per unit mass around the stellar surface at an arbitrary time step of the simulation. Boxes show $\langle q_{{\rm rad}} \rangle$ when computed on the hydrodynamical grid; the vertical resolution reaches 7 km around the peak. The solid line shows the result after inserting two additional horizontal layers in each hydrodynamical cell. The upper panel gives the average deviation between the two cases in the same units.

Open with DEXTER

5 The effects of scattering on the photospheric temperature structure of a solar-type star

5.1 The 3D hydrodynamical surface convection model

To investigate the effects of scattering on the atmosphere of a solar-type star, we conduct time-dependent radiative hydrodynamical simulations of the quiet surface, neglecting the effects of magnetic fields. We solve the fully compressible Navier-Stokes equations, the mass conservation equation and the energy equation, along with the time-independent radiative transfer equation (Eq. (4)); see, e.g., Stein & Nordlund (1998) and Nordlund et al. (2009) for further details. Our $240\times240\times226$ model covers a horizontal area of 6 Mm $\times$ 6 Mm at a constant resolution of 25 km, and extends approximately 700 km above and 2.8 Mm below the surface. The vertical resolution reaches 7 km around the radiative cooling peak and decreases in the optically thick and thin parts of the simulation; radiative transfer is thus resolved well enough that only ${\sim}3$ % of the rays would be affected by overshoots (see Sect. 3.1). We test the accuracy of the vertical resolution using the adaptive refinement, inserting two extra layers before each computation of radiative transfer. Local differences between the two calculations reach ${\sim}3$ $\times$ 10¹⁰ erg g^-1 s^-1, owing to the strong sensitivity of the heating rate per unit mass, $q_{{\rm rad}}\equiv Q_{{\rm rad}}/\rho$ , to the local temperature gradients in the highly inhomogeneous granulation flow. On the average, however, the change in radiative flux divergence is negligible (see the upper panel of Fig. 5), and the radiation field is well resolved on the hydrodynamical grid. Note the difference between the magnitude of the cooling peaks in Figs. 5 and 6: the 1D calculation is based on the mean structure; in the 3D case, the average over each depth layer in the 3D box is taken and thus includes lateral inhomogeneities produced by the granulation flow.

$\begin{figure} \par\includegraphics[width=7.5cm,clip]{14210f6a.eps}\hspace*{3mm} \includegraphics[width=7.5cm,clip]{14210f6b.eps} \vspace*{6mm} \end{figure}$

Figure 6:

Left: heating rates $q_{{\rm rad}}$ per unit mass as a function of monochromatic optical depth at 5000 Å, computed on the 1D mean structure with full opacity sampling for three cases: without scattering (solid line), with continuum scattering (dashed line), and with continuum and line scattering (dot-dashed line). The upper panel shows the deviations of the latter two cases from the computation without scattering (same axis units as in the lower panel). Right: same computation, but using mean opacities and scattering albedos in 12 bins for the radiative transfer computations.

Open with DEXTER

Horizontal boundaries are periodic to mimic an infinitely extended atmosphere, vertical boundaries at the top and bottom of the simulation box are open to minimize the interference with the granulation flow. Mass conservation is ensured at the bottom by keeping the gas pressure constant; the underlying convection zone is mimicked by setting the entropy of the inflowing gas. The upper atmosphere is stabilized by setting internal energies to a slowly evolving average at the top.

We approximate the wavelength integral (Eq. (8)) with 12 opacity bins to account for the depth-dependence and wavelength-dependence of the absorption and scattering coefficients. The simulation box extends far into the optically thin atmosphere with $\langle \tau_{5000}\rangle$ $\approx$ 10^-6, where irradiation $I_{{\rm top}}^{-}$ from above is negligible. Rosseland optical depths at the bottom typically reach $\langle \tau_{{\rm Ross}}\rangle$ $\approx$ 10⁷, where radiative transfer is entirely diffusive and the radiation field is completely thermalized. We therefore set the diffusion approximation $I_{{\rm bot}}^{+}=B_{{\rm bot}}+{\rm d}B/{\rm d}\tau$ for all ingoing intensities at the bottom.

The three simulations discussed in Sect. 5.3 have mean effective temperatures $T_{{\rm eff}}$ between 5804 K and 5811 K with average temporal fluctuations of about 13 K; they are thus slightly hotter than the Sun. For our purposes, there is no need to exactly reproduce the solar $T_{{\rm eff}}$ . The simulations yield time-series of snapshots spanning ${\sim}1$ h of stellar time each, covering several granule lifetimes ( $t\sim10$ min) and several periods of the dominant p-mode ( $t\sim5$ min). Our simulation box covers about 10 granules with typical sizes of the order of ${\sim}1$ Mm, allowing us to obtain a statistically meaningful sample of the surface flow in terms of the ergodic hypothesis. The model without scattering was computed with a coarser radiation time step of 0.2 s, keeping the radiation field constant during the intermediate hydrodynamical calculations. The slow evolution of the flow field and the locality of the Planck source function allow such reduction of the computation time in very good approximation.

5.2 Scattering in the 1D mean stratification

We first test the importance of scattering in the 1D mean stratification of our 3D model (the S=B case, see Sect. 5.3) by comparing the wavelength-integrated $q_{{\rm rad}}$ , using the full opacity-sampled spectrum. Radiative transfer was computed in 1D using a direct block matrix Feautrier-type solver with coherent scattering (for a detailed description see, e.g., Rutten 2003) and 4th order Radau quadrature for the integral over the polar angle. The left-hand side of Fig. 6 shows $q_{{\rm rad}}$ without scattering and S=B, with continuum scattering only, and with both continuum and line scattering (lower panel), as well as the deviations from the first case (upper panel).

Continuum scattering seems to have very little impact on $q_{{\rm rad}}$ for the given mean structure; the cooling is slightly stronger near the surface. This behavior is expected from the mostly large photon destruction probabilities $\epsilon _{\lambda }^{{\rm c}}$ shown in the upper panel of Fig. 4.

The differences are slightly larger when scattering is included in the line-blanketing: the small heating bump, where cool uprising gas is heated from beneath by hot granules (see the discussion in Stein & Nordlund 1998), and the cooling peak beneath the surface both slightly weaken, since the fraction of scattered photons in the line-blanketing does not contribute to heat exchange (cf. the right-hand side of Eq. (7)). The upper atmosphere, however, now shows slight heating of the mean structure.

We repeat the same test with the binned opacities, computing 1D radiative transfer with and without scattering for 12 mean opacities, photon destruction probabilities and bin-integrated Planck functions. The right panels of Fig. 6 compare again the three different cases. The binning has been optimized for matching sampled and binned $q_{{\rm rad}}$ in the S=B case (solid lines in the lower panels of Fig. 6). The continuum scattering calculation with opacity bins underestimates the cooling beneath the surface. The disparity increases further when line scattering is included; the relative deviations reach 7.5% in the cooling peak (dot-dashed lines in Fig. 6). However, the overall impact of scattering radiative transfer on the temperature structure of the 3D atmosphere above $\tau_{5000} \ga 10^{-3}$ is small (see Sect. 5.3 and Fig. 7), the same binning setup was therefore adopted for all three simulations. Higher up in the atmosphere, at $\tau_{5000} \la 10^{-3}$ , opacity binned radiative transfer shows slightly stronger heating of the gas.

5.3 Scattering in the mean 3D model

In order to assess the effects of continuum and line scattering, we perform three independent simulation runs: the first one treats radiation without scattering by adding all scattering opacity to the absorption opacity and assuming a Planck source function S=B. The second one includes continuum scattering in the source function and only adds line scattering opacity to the absorption opacity, and the third one includes scattering both in the continuum and in the line-blanketing. All three time series start from the same initial snapshot and span the exact same amount of simulation time. Snapshots are taken at regular intervals of $\Delta t_{{\rm sim}}=10$ s. We consider time steps at $t_{{\rm sim}}>8$ min after the initial snapshot to allow the atmosphere to adjust to any changes in the radiative heating rates. Exploiting the tight correlation between gas temperature T and vertical optical depth $\tau$ (Stein & Nordlund 1998), we interpolate the 3D temperature cube at each time step of the series onto surfaces with the same optical depth, using a reference $\tau$ -scale at 5000 Å. We then compute the average temperature of each surface in the 3D cube, which yields a 1D mean temperature profile for every snapshot. These profiles are finally averaged over time, and we obtain a very robust characteristic $T-\tau$ relation.

$\begin{figure} \par\includegraphics[width=7.8cm,clip]{14210f7.eps} \end{figure}$	Figure 7: Horizontal and temporal average of the mean temperature structure as a function of optical depth at 5000 Å without scattering (solid line), with continuum scattering (dashed), and with continuum and line scattering (dot-dashed). The upper panel shows the deviation from the first case.
Open with DEXTER

Figure 7 compares the resulting horizontal and temporal mean temperature profiles. The simulations without scattering and with continuum scattering have practically identical stratifications, as expected from the continuum photon destruction probabilities $\epsilon _{\lambda }^{{\rm c}}$ (Fig. 4) and the 1D test presented in the previous section; continuum scattering is therefore insignificant for the atmospheric stratification in solar-type stars.

The effects of scattering on line-blanketing in and below the photosphere are also rather weak (dot-dashed line in Fig. 7). The gas temperatures above $\tau_{5000} \ga 10^{-2}$ deviate up to 40 K from the stratification without scattering, resulting in a slightly steeper temperature gradient around the surface ( $\tau_{5000}=1$ ). Since our adopted binning setup overestimates the deviations for the 1D mean structure (right-hand side of Fig. 6), the impact of line scattering is probably even smaller at $\tau_{5000} \ga 10^{-2}$ . The temperature structure in the lower photosphere is thus hardly affected by scattering. The opposite is the case in the high photosphere and above ( $\tau_{5000} \la 10^{-4}$ ), where we observe temperatures that are about 350 K lower, resulting in a significantly steeper mean gradient.

5.4 Comparison of the 1D and 3D calculations and with other model atmospheres

$\begin{figure} \par\includegraphics[width=11.5cm,clip]{14210f8.eps} \end{figure}$

Figure 8:

Horizontal averages of the radiative heating rates $q_{{\rm rad}}$ for the continuum scattering case (circles) and the continuum and line scattering case (diamonds) as a function of simulation time t, at optical depths $\tau _{5000}=10^{-5}$ (upper panel), $\tau _{5000}=10^{-4}$ (center panel) and $\tau _{5000}=10^{-3}$ (lower panel). Dashed lines show the spatial and temporal averages for the continuum scattering case, where line scattering is treated as true absorption; dot-dashed lines show the spatial and temporal averages for the line scattering case. Dotted lines indicate zero heating.

Open with DEXTER

The effects of line scattering on the temperature structure of the 3D model seem to be opposite of 1D hydrostatic models in radiative equilibrium, where heating of the highest layers rather than cooling is observed. Indeed, the 1D calculations on the mean 3D atmosphere exhibit slight heating in this region when scattering is included (Fig. 6). The temperature gradient would therefore become shallower if the 1D calculations were iterated under the assumption of radiative equilibrium (see, e.g., the discussion in Rutten 2003).

The total radiative flux divergence includes several components: hot radiation from deeper layers at short wavelengths dominates the heating of the gas; the steep outward ${\rm d}B_{\lambda}/{\rm d}T$ gradient causes a positive growing (J-S) split. The effect declines in higher layers due to the rapidly decreasing opacity (cf. Eq. (7)). Strong LTE lines may heat or cool the higher atmosphere (since $J\approx B$ in deeper parts), depending on the spectral region and local temperature gradient, which determine the sign of the (J-S) split. Including coherent scattering in line-blanketing effectively reduces both radiative heating and cooling in high layers through the outwards decreasing $\epsilon^{{\rm l}}$ (see Fig. 1). As a consequence, strong resonance lines become unimportant for the temperature structure in high layers, and radiative heating at shorter wavelengths decreases.

In the 1D mean atmosphere, scattering-weakened line cooling shifts the total $q_{{\rm rad}}$ slightly towards positive values. The behavior of the 3D case can be understood by considering the dynamical nature of our 3D models. Following a derivation in Mihalas & Mihalas (1984), we insert the continuity equation

$\begin{displaymath}% \frac{D\rho}{{\rm D}t}+\rho\vec{\nabla}\cdot \vec{u}=0, \end{displaymath}$

(23)

where $\rho$ is the gas density, $D/{\rm D}t$ is the material derivative and $\vec{u}$ is the gas velocity, into the energy equation,

$\begin{displaymath}% \frac{De}{{\rm D}t}+\frac{P}{\rho}\vec{\nabla}\cdot \vec{u}=q_{{\rm rad}}, \end{displaymath}$

(24)

where e is the internal energy per unit mass, P is the gas pressure, and $q_{{\rm rad}}$ is the radiative heating rate per unit mass; we omit the viscous dissipation term for simplicity. The resulting expression,

$\begin{displaymath}% \frac{De}{{\rm D}t}-\frac{P}{\rho^{2}}\frac{D\rho}{{\rm D}t}=q_{{\rm rad}}, \end{displaymath}$

(25)

is the first law of thermodynamics. An upflowing (downflowing) gas parcel cools (heats) through expansion (compression) represented by the $D\rho/{\rm D}t$ term in Eq. (25), and is exposed to radiative heating through the $q_{{\rm rad}}$ term. Equation (25) is equivalent to the expression

$\begin{displaymath}% T\frac{Ds}{{\rm D}t}=q_{{\rm rad}}, \end{displaymath}$

(26)

where T is the gas temperature and s is the entropy per unit mass, and it is immediately clear that gas motion is adiabatic when $q_{{\rm rad}}\rightarrow 0$ . In the photosphere of the 3D simulation, temperatures are not affected by scattering. In the upper atmosphere, below $\tau_{5000}\approx10^{-4}$ , scattering strongly reduces the line-blanketing. Small or vanishing heating rates $q_{{\rm rad}}$ cause the temperature stratification to steepen towards an adiabatic gradient.

Figure 8 compares the time evolution of radiative heating rates at three optical depths in the atmosphere, averaged over surfaces of constant optical depth to approximately account for vertical gas motion. The plot shows a sequence of snapshots taken at regular simulation time intervals of 10 s; the zero point on the abscissa is arbitrary. At $\tau _{5000}=10^{-3}$ and $\tau _{5000}=10^{-4}$ (lower and center panels), the continuum scattering case (circles) and the continuum and line scattering case (diamonds) exhibit similar positive heating rates on the average (dashed and dot-dashed lines) and thus similar average temperatures (Fig. 7). Line scattering radiative transfer produces slightly stronger mean heating at $\tau _{5000}=10^{-3}$ , but fluctuates with lower amplitude. At $\tau _{5000}=10^{-5}$ , $q_{{\rm rad}}$ practically vanishes on the time average in the line scattering case, but there is still significant radiative heating with line scattering as true absorption. Note the dynamical variation of the sequences: contrary to 1D hydrostatic models, where the radiation field is time-independent by definition, the evolution of the 3D simulations produces fluctuating radiative heating.

$\begin{figure} \par\includegraphics[width=12cm,clip]{14210f9.eps} \end{figure}$	Figure 9: Temperature histograms at three different geometrical heights z above the optical surface, integrated over each simulation run. Solid lines show the radiative transfer computation with continuum scattering, dot-dashed lines the case where continuum and line scattering are included. Each temperature distribution is normalized.
Open with DEXTER

Wedemeyer et al. (2004) presented 3D radiation-hydrodynamical simulations of the solar atmosphere that include a chromosphere, using radiative transfer without scattering and solving the equation only for the Rosseland mean opacity to suppress radiative cooling by strong LTE lines. They found an increasing asymmetry of the gas temperature distribution with increasing height above the surface, and a bifurcation in the chromosphere. Wedemeyer et al. (2004) further observed that treating strong spectral lines as true absorption with the opacity binning method reduces the amplitude of temperature fluctuations, which are caused by outward propagating acoustic waves, resulting in unrealistically low maximum temperatures in high layers. Skartlien (2000) investigated scattering radiative transfer in the chromosphere, comparing radiative heating with and without scattering, and came to the conclusion that including line scattering reduces this damping effect of LTE lines.

Our simulations do not include a chromosphere; the internal energy at the top boundary is set to a slowly evolving mean instead. In the line scattering case, where radiative transfer has only weak influence on the gas, the temperature gradient is sensitive to this boundary condition and thus not well-constrained. However, this does not compromise our conclusions, since the boundary is free to adapt to any upward or downward shift in the mean energies of the gas beneath.

Figure 9 shows temperature distributions of the simulations with continuum scattering and with continuum and line scattering at three different heights above the surface. Our simulations do not reach the same geometrical heights as those of Skartlien (2000) and Wedemeyer et al. (2004), and we use a more realistic radiative transfer treatment with 12 opacity bins. We find a similarly growing asymmetry in the temperature distribution of the line scattering simulation in the outer layers (cf. Fig. 7 in Wedemeyer et al. 2004). Treating strong lines as absorbers shifts the mean temperature upward and removes the high temperature tail of the distribution, in qualitative agreement with the findings of Skartlien (2000) and Wedemeyer et al. (2004).

Figure 10 shows horizontal and temporal averages of the relative temperature fluctuations, which we define as

$\begin{displaymath}% \frac{\Delta T_{{\rm rms}}}{\langle T\rangle}=\frac{\sqrt{\langle (T-\langle T\rangle )^{2}\rangle}}{\langle T\rangle} \end{displaymath}$

(27)

in every geometrical depth layer (cf. Eq. (2) and Fig. 9 in Wedemeyer et al. 2004). The comparison between the cases with continuum scattering and with continuum and line scattering confirms the damping of temperature fluctuations through line absorption. Note the decreasing $\Delta T_{{\rm rms}}$ at the top of the simulation, which is induced by the hydrodynamical boundary conditions.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{14210f10.eps} \end{figure}$	Figure 10: Horizontal and temporal averages of the relative temperature fluctuations $\Delta T_{{\rm rms}}/\langle T \rangle$ as a function of atmospheric height, computed with continuum scattering (solid line) and with continuum and line scattering (dot-dashed line).
Open with DEXTER

We conclude that line scattering is an important ingredient for model atmospheres of solar-type stars that include a chromosphere; while gray radiative transfer reduces damping through strong LTE lines, it cannot produce a realistic temperature structure.

Anderson (1989) presented simplified, and Short & Hauschildt (2005) presented full 1D non-LTE line-blanketing calculations, respectively, for hydrostatic model atmospheres of solar-type stars. The departures of line-blanketing from LTE through iron-group elements heat up the atmosphere in the height range $10^{-6} \la \tau_{5000} \la 10^{-2}$ . Our 3D model predicts a predominant temperature decrease as we discussed above. However, it is not clear how departures from LTE in the absorber populations through the ionization balance etc. would affect the atmospheric structure in our 3D simulations, making a direct comparison with the 1D models difficult.

Doppler shifts may have a significant influence on line absorption in higher layers, which expose line cores to hot radiation from deeper in. Vögler et al. (2004) estimated the effects to be insignificant in the photosphere, but their work was based on 1D tests. The large scattering albedo of strong resonance lines, however, should reduce the impact of Doppler shifts higher up.

6 Conclusions

We presented a 3D radiative transfer method with coherent scattering for time-dependent (M)HD simulations of stellar atmospheres with the new BIFROST code (Gudiksen et al., in prep.). The simulations are parallelized through domain-decomposition to take advantage of large-scale computer clusters. The solver is based on short characteristics and the Gauss-Seidel scheme for an iterative computation of the radiation field and the radiative flux divergence in the whole simulation domain. We use monotonic interpolation to reduce the numerical diffusion effect of short characteristics and represent the source function integral with Bézier polynomials to suppress interpolation overshoots. A partial grid refinement scheme is included to improve the resolution of the radiative transfer computation where strong vertical opacity gradients occur. The wavelength integral is treated in the opacity binning approximation, using 12 bins that divide the opacity spectrum by formation height and wavelength.

The effects of coherent scattering on the temperature structure of a solar-type star are investigated with 3D time-dependent hydrodynamical simulations of magnetically quiet surface convection, including Rayleigh scattering and electron scattering in the continuum and estimated line scattering using the van Regemorter formula. While continuum scattering processes are not important for the mean temperature stratification, we find lower temperatures in the upper atmosphere when scattering is included in the line-blanketing. 3D radiative-hydrodynamical atmospheres thus show the opposite behavior of 1D hydrostatic atmospheres in radiative equilibrium, where scattering in strong lines effectively heats the outer layers.

3D LTE models of solar surface convection have been very successful at reproducing various observational tests, and our results indicate that the solar photosphere is indeed well represented when scattering is not included in radiative transfer. It therefore seems that a refined treatment of the line-blanketing through, e.g., opacity distribution functions or opacity sampling will be the next significant step to improve the realism of 3D radiative-hydrodynamical model atmospheres. Scattering radiative transfer is nevertheless an important ingredient of consistent 3D MHD models of the solar chromosphere, transition region and corona.

While it is not unexpected to see only small differences in the photospheres of solar-type stars when scattering is taken into account, this is likely to change for the much less dense atmospheres of giants, where the importance of Rayleigh scattering increases. The case of metal-poor giants is particularly interesting in that respect, owing to their significance for understanding galactic chemical evolution and the origin of the elements.

Acknowledgements

This research project has been supported by a Marie Curie Early Stage Research Training Fellowship of the European Community's Sixth Framework Programme under contract number MEST-CT-2005-020395: The USO-SP International School for Solar Physics.

Appendix A: Bézier interpolation of source functions and opacities

The discrete formal solution (Eq. (14)) of the radiative transfer equation (Eq. (4)) requires interpolating the source function $S(\tau)$ along the short characteristic. While linear interpolation never overshoots, its accuracy is not sufficient in optically thick media, since the discrete expression is not equivalent to the diffusion approximation. Second-order interpolation significantly improves the accuracy, but suffers from strong overshoots where $\Delta S/\Delta\tau$ gradients change rapidly between the upwind and downwind halves of the characteristic. In extreme cases, this can even destabilize the solver and produce spikes in the local flux divergences.

Bézier-type interpolation techniques allow for a direct detection and suppression of such overshoots by virtue of a control point in the polynomial which shapes its curve (see Auer 2003). A second-order Bézier polynomial may be written in the parameterized form

$\begin{displaymath}% S(t)=S_{{\rm u}}(1-t)^{2}+S_{0}t^{2}+2S_{{\rm c}}t(1-t), \end{displaymath}$

(A.1)

where $S_{{\rm u}}$ and S₀ are the source functions at the upwind end and the center point of the characteristic, between which interpolation is needed, $t=(\tau-\tau_{{\rm u}})/(\tau_{0}-\tau_{{\rm u}})$ is the curve parameter, and $S_{{\rm c}}$ is the control point. A Beziér curve is always bounded by the convex hull of the three points $S_{{\rm u}}$ , $S_{{\rm c}}$ and S₀. Using the abbreviations

$\begin{displaymath}% \Delta\tau_{{\rm u}}=\tau_{0}-\tau_{{\rm u}}; \hspace{0.5cm} \Delta\tau_{{\rm d}}=\tau_{{\rm d}}-\tau_{0} \end{displaymath}$

(A.2)

for the optical depths along the characteristic and choosing the control point

$\begin{eqnarray*}S_{{\rm c}} &=& S_{0}-\frac{\Delta\tau_{{\rm u}}}{2}S'_{0}\\ ... ..._{{\rm d}}}\frac{S_{{\rm d}}-S_{0}}{\Delta\tau_{{\rm d}}}\right) \end{eqnarray*}$

yields second-order interpolation of S, which now also depends on the source function $S_{{\rm d}}$ at the downwind end. Introducing the three functions

$\begin{eqnarray*}&& U_{0}(t) = 1-{\rm e}^{-t} \\ && U_{1}(t) = t-U_{0}(t) \\ && U_{2}(t) = t^{2}-2U_{1}(t), \end{eqnarray*}$

and evaluating the integral of the Bézier polynomial results in the familiar second-order integration constants for Eq. (14),

$\begin{eqnarray*}&& \Psi_{{\rm u}} = U_{0}(\Delta\tau_{{\rm u}})+\frac{U_{2}(\De... ...lta\tau_{{\rm d}}(\Delta\tau_{{\rm u}}+\Delta\tau_{{\rm d}})}\\ \end{eqnarray*}$

(cf. Eqs. (8) and (9) in Kunasz & Auer 1988). If the source functions $S_{{\rm u}}$ , S₀ and $S_{{\rm d}}$ have an extremum at S₀, choosing $S_{{\rm c}}=S_{0}$ enforces S'₀=0, yielding the constants

$\begin{eqnarray*}&& \Psi_{{\rm u}} = U_{0}(\Delta\tau_{{\rm u}})+\frac{U_{2}(\De... ...au_{{\rm u}})}{\Delta\tau_{{\rm u}}^2}\\ && \Psi_{{\rm d}} = 0. \end{eqnarray*}$

Overshoots are avoided by limiting $S_{{\rm c}}$ to the range of the data points: $\min(S_{{\rm u}},S_{0})\leq S_{{\rm c}}\leq\max(S_{{\rm u}},S_{0})$ . If $S_{{\rm c}}$ lies outside these limits, choosing $S_{{\rm c}}=S_{{\rm u}}$ results in the constants

$\begin{eqnarray*}&& \Psi_{{\rm u}} = U_{0}(\Delta\tau_{{\rm u}})-\frac{U_{2}(\De... ..._{{\rm u}})}{\Delta\tau_{{\rm u}}^{2}}\\ && \Psi_{{\rm d}} = 0. \end{eqnarray*}$

Note that, contrary to the first two cases, suppressing such overshoots leads to discontinuous left-hand and right-hand derivatives at S₀.

Optical depths are computed in an analogue fashion to avoid negative results. A second-order Bézier polynomial $\chi(s)$ interpolates opacities over the path length $\Delta s$ along the ray; integration over s yields the optical depth interval

$\begin{displaymath}% \Delta\tau=\int_{0}^{\Delta s}\chi(\sigma){\rm d}\sigma=\fr... ...elta s}{3}\left(\chi_{{\rm u}}+\chi_{0}+\chi_{{\rm c}}\right), \end{displaymath}$

(A.3)

where the control point $\chi_{{\rm c}}$ is selected according to the same criteria as discussed above for $S_{{\rm c}}$ .

Appendix B: Local cubic monotonic interpolation

The radiative transfer solver uses local cubic interpolation for interpolating data from the hydrodynamical grid onto the characteristics grid. This choice of method improves the accuracy compared to linear interpolation, while ensuring local control of the interpolating polynomial to reduce artifacts. In addition to being a one-pass algorithm, the method also exhibits good computational performance through its instruction-per-data ratio, which is well-suited for modern multi-core CPUs, where high computation speeds are contrasted with slow memory access.

2D interpolation is approximated by consecutive 1D interpolation using a cubic polynomial

f(t)=at³+bt²+ct+d,

(B.1)

with the curve parameter $t\in[0,1]$ . The coefficients a-d depend on the adjacent data points f₁ and f₂ and their derivatives f'₁ and f'₂. Inserting the data and reordering the terms, the polynomial may be rewritten in the form

$\begin{displaymath}% f(t)=\alpha(t)f_{1}+\beta(t)f_{2}+\gamma(t)f'_{1}+\delta(t)f'_{2}, \end{displaymath}$

(B.2)

where the interpolation weights $\alpha$ , $\beta$ , $\gamma$ and $\delta$ now depend on the parameter t:

$\begin{eqnarray*}&& \alpha(t) = 2t^{3}-3t^{2}+1\\ && \beta(t) = 3t^{2}-2t^{3}\\... ...ight)\Delta x\\ && \delta(t) = \left(t^{3}-t^{2}\right)\Delta x \end{eqnarray*}$

with the grid spacing $\Delta x$ between the two data points. The shape of the curve is defined by the derivatives f'₁ and f'₂. A natural choice is the mean of the left-handed and right-handed difference quotients $f'_{{\rm L}}$ and $f'_{{\rm R}}$ at both end points. An unweighted arithmetic mean leads to wiggles and overshoots where strong gradients appear. We therefore adopt a recipe by Fritsch & Butland (1984), which uses a weighted harmonic mean

$\begin{displaymath}% f'=\left\{ \begin{array}{cc} \frac{f'_{{\rm L}}f'_{{\rm R}}... ...space{0.5cm} f'_{{\rm L}}f'_{{\rm R}} \le 0 \end{array}\right. \end{displaymath}$

(B.3)

with the weighting factor

$\begin{displaymath}% \alpha=\frac{1}{3}\left(1+\frac{\Delta x_{{\rm R}}}{\Delta x_{{\rm L}}+\Delta x_{{\rm R}}}\right), \end{displaymath}$

(B.4)

which depends on the grid spacing $\Delta x_{{\rm L}}$ and $\Delta x_{{\rm R}}$ on the left and right sides of the data point. The weighted harmonic mean biases f' towards the smaller of the two difference quotients $f'_{{\rm L}}$ and $f'_{{\rm R}}$ where strong gradients occur, effectively suppressing overshoots.

Quadratic interpolation uses only one of the two derivatives f'₁ and f'₂, depending on the interpolation parameter t. The interpolation coefficients are

	$\displaystyle t\le\frac{1}{2}: \hspace{2.5cm} t>\frac{1}{2}:$
	$\displaystyle \alpha(t) =1-t^2 \hspace{1.8cm} \alpha(t) = (1-t)^2$
	$\displaystyle \beta(t) = t^2 \hspace{2.4cm} \beta(t) = t(2-t)$
	$\displaystyle \gamma(t) = t(1-t)\Delta x \hspace{1.2cm} \gamma(t) = 0$
	$\displaystyle \delta(t) = 0 \hspace{2.4cm} \delta(t) = t(t-1)\Delta x.$	(B.5)

Appendix C: Line scattering with the van Regemorter formula

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{14210fc1.eps} \vspace*{3mm} \end{figure}$	Figure C.1: Dependence of the tabulated Gaunt factor integral $P(\Delta E_{\lambda },T)$ for collisions of electrons with neutral atoms in the van Regemorter (1962) formula on the transition energy $\Delta E_{\lambda }$ and the gas temperature T. The dashed lines mark the boundaries for typical values of $\Delta E_{\lambda }/kT$ found in solar surface convection simulations (Sect. 5.1).
Open with DEXTER

The photon destruction probabilities in line transitions may be estimated using the semi-empirical van Regemorter (1962) formula to obtain electron collision rates, following the discussion in Skartlien (2000). Neglecting other contributions from, e.g., collisions with neutral hydrogen atoms, the de-excitation rate for electron collisions according to this formula is given by

$\begin{displaymath}% C_{21}\sim\lambda^3N_{{\rm e}}T^{-1/2}A_{21}P(\Delta E_{\lambda},T), \end{displaymath}$

(C.1)

where $\lambda$ is the transition wavelength, $N_{{\rm e}}$ is the electron density, T is the gas temperature, and A₂₁ is the Einstein coefficient for the corresponding spontaneous radiative transition. The function $P(\Delta E_{\lambda },T)$ abbreviates the velocity integral over the empirically calibrated Gaunt factor of the scattered electron, and depends on the transition energy $\Delta E_{\lambda }$ and the gas temperature T. We adopt the tabulated data for neutral atoms of van Regemorter (1962), see Fig. C.1.

The photon destruction probability for a two-level atom is given by

$\begin{displaymath}% \epsilon_{\lambda}=\frac{\kappa_{\lambda}}{\kappa_{\lambda}... ...gma_{\lambda}}=\frac{C_{21}}{C_{21}+A_{21}+B_{21}B_{\lambda}}, \end{displaymath}$

(C.2)

where B₂₁ is the rate for induced de-excitation, and $B_{\lambda}$ is the Planck function. Neglecting the induced de-excitation term, Eq. (C.2) simplifies to

$\begin{displaymath}% \epsilon_{\lambda}\approx\frac{1}{1+A_{21}/C_{21}}\cdot \end{displaymath}$

(C.3)

$\epsilon _{\lambda }$ is independent of the actual transition after inserting the van Regemorter formula Eq. (C.1), and thus only a function of $\lambda$ , $N_{{\rm e}}$ and T.

Line opacities in stellar spectra often combine contributions from many transitions at a given wavelength. The total monochromatic photon destruction probability of an opacity sample at wavelength $\lambda$ is given by the sum over all transitions,

$\begin{displaymath}% \epsilon^{{\rm l}}_{\lambda}=\frac{\sum_{i}\kappa^{{\rm l}}_{\lambda,i}}{\sum_{i}\chi^{{\rm l}}_{\lambda,i}}\cdot \end{displaymath}$

(C.4)

Inserting Eq. (C.3), thus assuming the above mentioned approximations, yields

$\begin{displaymath}% \epsilon^{{\rm l}}_{\lambda}\approx\frac{\sum_{i}\epsilon_{... ...da,i}}{\sum_{i}\chi^{{\rm l}}_{\lambda,i}}=\epsilon_{\lambda}, \end{displaymath}$

(C.5)

where the equality holds since $\epsilon _{\lambda }$ is independent of the actual transition i. The absorption and scattering contributions $\kappa^{{\rm l}}_{\lambda}$ and $\sigma^{{\rm l}}_{\lambda}$ to each opacity sample $\chi^{{\rm l}}_{\lambda}$ are then isolated using $\epsilon^{{\rm l}}_{\lambda}$ and added to the coefficients of the continuum processes (see Table D.1).

Appendix D: Continuum opacity sources

Table D.1: Continuum opacity sources.

References

Allende Prieto, C., Asplund, M., García López, R. J., & Lambert, D. L. 2002, ApJ, 567, 544 [NASA ADS] [CrossRef] [Google Scholar]
Allison, A. C., & Dalgarno, A. 1969, AD, 1, 91 [NASA ADS] [CrossRef] [Google Scholar]
Anderson, L. S. 1989, ApJ, 339, 558 [NASA ADS] [CrossRef] [Google Scholar]
Asplund, M., Nordlund, Å., Trampedach, R., Allende Prieto, C., & Stein, R. F. 2000, A&A, 359, 729 [NASA ADS] [Google Scholar]
Asplund, M., Grevesse, N., & Sauval, A. J. 2005, in Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis, ed. T. G. Barnes, & F. N. Bash, ASP Conf. Ser., 336, 25 [Google Scholar]
Auer, L. 2003, in Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 3 [Google Scholar]
Bautista, M. A. 1997, A&AS, 122, 167 [NASA ADS] [CrossRef] [EDP Sciences] [PubMed] [Google Scholar]
Bautista, M. A. 1999, A&AS, 137, 529 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
Bell, K. L. 1980, J. Phys. B, 13, 1859 [NASA ADS] [CrossRef] [Google Scholar]
Bell, K. L., & Berrington, K. A. 1987, J. Phys. B, 20, 801 [NASA ADS] [CrossRef] [Google Scholar]
Bell, K. L., Hibbert, A., & Berrington, K. A. 1988, J. Phys. B, 21, 2319 [NASA ADS] [CrossRef] [Google Scholar]
Borysow, A., Jørgensen, U. G., & Fu, Y. 2001, JQSRT, 68 [Google Scholar]
Broad, J. T., & Reinhardt, W. P. 1976, Phys. Rev. A, 14, 2159 [NASA ADS] [CrossRef] [Google Scholar]
Bruls, J. H. M. J., Vollmöller, P., & Schüssler, M. 1999, A&A, 348, 233 [NASA ADS] [Google Scholar]
Caffau, E., Ludwig, H., Steffen, M., et al. 2008, A&A, 488, 1031 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cannon, C. J. 1973, ApJ, 185, 621 [NASA ADS] [CrossRef] [Google Scholar]
Carlson, B. G. 1963, Methods in Computational Physics, 1, Statistical Physics, ed. B. Alder, S. Fernbach, & M. Rotenberg (New York and London: Academic Press), 1 [Google Scholar]
Collet, R., Asplund, M., & Trampedach, R. 2007, A&A, 469, 687 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Däppen, W., Mihalas, D., Hummer, D. G., & Mihalas, B. W. 1988, ApJ, 332, 261 [NASA ADS] [CrossRef] [Google Scholar]
Doyle, R. O. 1968, ApJ, 153, 987 [NASA ADS] [CrossRef] [Google Scholar]
Dravins, D., & Nordlund, A. 1990, A&A, 228, 203 [NASA ADS] [Google Scholar]
Feautrier, P. 1964, Comptes Rendus Académie des Sciences, 258, 3189 [NASA ADS] [Google Scholar]
Fritsch, F. N., & Butland, J. 1984, SIAM J. Sci. Stat. Comput., 5, 300 [CrossRef] [MathSciNet] [Google Scholar]
Gavrila, M. 1967, Phys. Rev., 163, 147 [NASA ADS] [CrossRef] [Google Scholar]
Gustafsson, M., & Frommhold, L. 2001, ApJ, 546, 1168 [NASA ADS] [CrossRef] [Google Scholar]
Gustafsson, M., & Frommhold, L. 2003, A&A, 400, 1161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gustafsson, B., Bell, R. A., Eriksson, K., & Nordlund, A. 1975, A&A, 42, 407 [NASA ADS] [Google Scholar]
Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hansteen, V. H. 2004, in Multi Wavelength Investigations of Solar Activity, ed. A. Stepanov, E. Benevolenskaya, & E. Kosovichev, IAU Symp., 223, 385 [Google Scholar]
Hansteen, V. H., Carlsson, M., & Gudiksen, B. 2007, in The Physics of Chromospheric Plasmas, ed. P. Heinzel, I. Dorotovic, & R. J. Rutten, ASP Conf. Ser., 368, 107 [Google Scholar]
Hauschildt, P. H., & Baron, E. 2006, A&A, 451, 273 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hauschildt, P. H., Allard, F., & Baron, E. 1999, ApJ, 512, 377 [NASA ADS] [CrossRef] [Google Scholar]
Heinemann, T., Dobler, W., Nordlund, Å., & Brandenburg, A. 2006, A&A, 448, 731 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Holweger, H., & Mueller, E. A. 1974, Sol. Phys., 39, 19 [NASA ADS] [CrossRef] [Google Scholar]
Hummer, D. G., & Mihalas, D. 1988, ApJ, 331, 794 [NASA ADS] [CrossRef] [Google Scholar]
John, T. L. 1975, MNRAS, 172, 305 [NASA ADS] [Google Scholar]
John, T. L. 1995, J. Phys. B, 28, 23 [NASA ADS] [CrossRef] [Google Scholar]
Jørgensen, U. G., Hammer, D., Borysow, A., & Falkesgaard, J. 2000, A&A, 361, 283 [NASA ADS] [Google Scholar]
Karzas, W. J., & Latter, R. 1961, ApJS, 6, 167 [NASA ADS] [CrossRef] [Google Scholar]
Kunasz, P., & Auer, L. H. 1988, JQSRT, 39, 67 [NASA ADS] [CrossRef] [Google Scholar]
Kurucz, R. L. 1979, ApJS, 40, 1 [NASA ADS] [CrossRef] [Google Scholar]
Kurucz, R. L. 1996, in Stellar Surface Structure, ed. A. Wehlau, D. F. Gray, & J. Rice, IAU Symp., 176, 523 [Google Scholar]
Kurucz, R. L., van Dishoeck, E. F., & Tarafdar, S. P. 1987, ApJ, 322, 992 [NASA ADS] [CrossRef] [Google Scholar]
Langhoff, P. W., Sims, J., & Corcoran, C. T. 1974, Phys. Rev. A, 10, 829 [NASA ADS] [CrossRef] [Google Scholar]
Ludwig, H.-G. 1992, Ph.D. Thesis, Universität Kiel [Google Scholar]
Mihalas, D. 1978, Stellar atmospheres, 2nd edn, ed. J. Hevelius (San Francisco: W. H. Freeman and Co.) [Google Scholar]
Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics, ed. D. Mihalas, & B. W. Mihalas (New York: Oxford University Press) [Google Scholar]
Mihalas, D., Auer, L. H., & Mihalas, B. R. 1978, ApJ, 220, 1001 [NASA ADS] [CrossRef] [Google Scholar]
Mihalas, D., Däppen, W., & Hummer, D. G. 1988, ApJ, 331, 815 [NASA ADS] [CrossRef] [Google Scholar]
Mihalas, D., Hummer, D. G., Mihalas, B. W., & Däppen, W. 1990, ApJ, 350, 300 [NASA ADS] [CrossRef] [Google Scholar]
Muthsam, H. J., Kupka, F., Löw-Baselli, B., et al. 2010, New Astron., 15, 460 [NASA ADS] [CrossRef] [Google Scholar]
Nahar, S. N. 1995, ApJS, 101, 423 [NASA ADS] [CrossRef] [Google Scholar]
Nahar, S. N. 1996, Phys. Rev. A, 53, 1545 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Nahar, S. N. 1998, Phys. Rev. A, 58, 3766 [NASA ADS] [CrossRef] [Google Scholar]
Nahar, S. N. 2000, ApJS, 126, 537 [NASA ADS] [CrossRef] [Google Scholar]
Nahar, S. N. 2002, Phys. Rev. A, 65, 052702 [NASA ADS] [CrossRef] [Google Scholar]
Nahar, S. N., & Pradhan, A. K. 1991, Phys. Rev. A, 44, 2935 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Nahar, S. N., & Pradhan, A. K. 1993, J. Phys. B Atom. Mol. Phys., 26, 1109 [NASA ADS] [CrossRef] [Google Scholar]
Nahar, S. N., & Pradhan, A. K. 1994a, J. Phys. B Atom. Mol. Phys., 27, 429 [NASA ADS] [CrossRef] [Google Scholar]
Nahar, S. N., & Pradhan, A. K. 1994b, Phys. Rev. A, 49, 1816 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Nahar, S. N., & Pradhan, A. K. 1997, ApJS, 111, 339 [NASA ADS] [CrossRef] [Google Scholar]
Nordlund, A. 1982, A&A, 107, 1 [NASA ADS] [Google Scholar]
Nordlund, A., & Stein, R. F. 1991, in Stellar Atmospheres - Beyond Classical Models, NATO ASIC Proc., 341, 263 [Google Scholar]
Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Liv. Rev. Sol. Phys., 6, 2 [Google Scholar]
Olson, G. L., & Kunasz, P. B. 1987, JQSRT, 38, 325 [NASA ADS] [CrossRef] [Google Scholar]
Peach, G. 1970, MmRAS, 73, 1 [NASA ADS] [Google Scholar]
Pereira, T. M. D., Asplund, M., & Kiselman, D. 2009a, A&A, 508, 1403 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pereira, T. M. D., Kiselman, D., & Asplund, M. 2009b, A&A, 507, 417 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Peytremann, E. 1974, A&A, 33, 203 [NASA ADS] [Google Scholar]
Ramírez, I., Allende Prieto, C., Koesterke, L., Lambert, D. L., & Asplund, M. 2009, A&A, 501, 1087 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ramsbottom, C. A., Bell, K. L., & Berrington, K. A. 1992, J. Phys. B Atom. Mol. Phys., 25, 1443 [NASA ADS] [CrossRef] [Google Scholar]
Rutten, R. J. 2003, Radiative Transfer in Stellar Atmospheres, 8th edn. (Utrecht University lecture notes) [Google Scholar]
Seaton, M. J., Yan, Y., Mihalas, D., & Pradhan, A. K. 1994, MNRAS, 266, 805 [NASA ADS] [CrossRef] [Google Scholar]
Short, C. I., & Hauschildt, P. H. 2005, ApJ, 618, 926 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Skartlien, R. 2000, ApJ, 536, 465 [NASA ADS] [CrossRef] [Google Scholar]
Stancil, P. C. 1994, ApJ, 430, 360 [NASA ADS] [CrossRef] [Google Scholar]
Stein, R. F., & Nordlund, A. 1998, ApJ, 499, 914 [NASA ADS] [CrossRef] [Google Scholar]
Trujillo Bueno, J., & Fabiani Bendicho, P. 1995, ApJ, 455, 646 [NASA ADS] [CrossRef] [Google Scholar]
van Regemorter, H. 1962, ApJ, 136, 906 [NASA ADS] [CrossRef] [Google Scholar]
Victor, G. A., & Dalgarno, A. 1969, J. Chem. Phys., 50, 2535 [NASA ADS] [CrossRef] [Google Scholar]
Vögler, A., Bruls, J. H. M. J., & Schüssler, M. 2004, A&A, 421, 741 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Wedemeyer, S., Freytag, B., Steffen, M., Ludwig, H., & Holweger, H. 2004, A&A, 414, 1121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Wehrse, R., & Kalkofen, W. 2006, A&A Rev., 13, 3 [NASA ADS] [Google Scholar]
Wishart, A. W. 1979, MNRAS, 187, 59P [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table D.1: Continuum opacity sources.

All Figures

$\begin{figure} \par\includegraphics[width=7.4cm,clip]{14210f1a.eps}\hspace*{7mm} \includegraphics[width=7.55cm,clip]{14210f1b.eps} \end{figure}$	Figure 1: Horizontal mean photon destruction probability $\epsilon$ for three bins representing continuum and weak lines (1), intermediate lines (2) and strong lines (3) in the UV, plotted as a function of geometrical depth (left) and optical depth in the respective bin (right, the average is taken over surfaces with the same vertical optical depth). The dotted line marks the stellar surface (left) and unit optical depth in each bin (right).
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=11.3cm,clip]{14210f2.eps} \end{figure}$	Figure 2: Convergence of the source function for bins 1-3 (see Fig. 1) during a simulation run without domain decomposition (left column), with $2\times 2\times 2$ decomposition (center column) and $3\times 3\times 3$ decomposition (right column), and with a time step of $\Delta t=0.03$ s (upper row), $\Delta t=0.04$ s (center row) and $\Delta t=0.08$ s ( lower row). Line styles represent the same bins as in Fig. 1. Thin lines: relative source function correction $\Delta S$ after n iterations with respect to S^n-1 from the previous iteration n-1. Thick lines: relative source function correction $\Delta S$ with respect to the ``true'' solution $S^{\infty }$ . Dotted lines mark the threshold C beneath which convergence is assumed (see text).
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=6.6cm,clip]{14210f3a.eps}\hspace{7mm}... ...7mm} \includegraphics[width=6.6cm,clip]{14210f3d.eps} \vspace{5mm} \end{figure}$	Figure 3: Numerical diffusion of a searchlight beam with rectangular cross-section using linear (upper right), local cubic (lower left) and local cubic monotonic (lower right) interpolation, compared to the exact solution (upper left).
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=12cm,clip]{14210f4.eps} \end{figure}$	Figure 4: Wavelength and depth dependence of the continuum photon destruction probabilities $\epsilon _{\lambda }^{{\rm c}}$ (upper panel), the van Regemorter (1962) line photon destruction probabilities $\epsilon _{\lambda }^{{\rm l}}$ (center panel) and the total photon destruction probabilities $\epsilon _{\lambda }$ (lower panel) for the mean solar-type stratification. The zero point on the depth axis marks the stellar surface.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8cm,clip]{14210f5.eps} \end{figure}$	Figure 5: Horizontal average heating rate per unit mass around the stellar surface at an arbitrary time step of the simulation. Boxes show $\langle q_{{\rm rad}} \rangle$ when computed on the hydrodynamical grid; the vertical resolution reaches 7 km around the peak. The solid line shows the result after inserting two additional horizontal layers in each hydrodynamical cell. The upper panel gives the average deviation between the two cases in the same units.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=7.5cm,clip]{14210f6a.eps}\hspace{3mm} \includegraphics[width=7.5cm,clip]{14210f6b.eps} \vspace{6mm} \end{figure}$	Figure 6: Left: heating rates $q_{{\rm rad}}$ per unit mass as a function of monochromatic optical depth at 5000 Å, computed on the 1D mean structure with full opacity sampling for three cases: without scattering (solid line), with continuum scattering (dashed line), and with continuum and line scattering (dot-dashed line). The upper panel shows the deviations of the latter two cases from the computation without scattering (same axis units as in the lower panel). Right: same computation, but using mean opacities and scattering albedos in 12 bins for the radiative transfer computations.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=7.8cm,clip]{14210f7.eps} \end{figure}$	Figure 7: Horizontal and temporal average of the mean temperature structure as a function of optical depth at 5000 Å without scattering (solid line), with continuum scattering (dashed), and with continuum and line scattering (dot-dashed). The upper panel shows the deviation from the first case.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=11.5cm,clip]{14210f8.eps} \end{figure}$	Figure 8: Horizontal averages of the radiative heating rates $q_{{\rm rad}}$ for the continuum scattering case (circles) and the continuum and line scattering case (diamonds) as a function of simulation time t, at optical depths $\tau _{5000}=10^{-5}$ (upper panel), $\tau _{5000}=10^{-4}$ (center panel) and $\tau _{5000}=10^{-3}$ (lower panel). Dashed lines show the spatial and temporal averages for the continuum scattering case, where line scattering is treated as true absorption; dot-dashed lines show the spatial and temporal averages for the line scattering case. Dotted lines indicate zero heating.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=12cm,clip]{14210f9.eps} \end{figure}$	Figure 9: Temperature histograms at three different geometrical heights z above the optical surface, integrated over each simulation run. Solid lines show the radiative transfer computation with continuum scattering, dot-dashed lines the case where continuum and line scattering are included. Each temperature distribution is normalized.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{14210f10.eps} \end{figure}$	Figure 10: Horizontal and temporal averages of the relative temperature fluctuations $\Delta T_{{\rm rms}}/\langle T \rangle$ as a function of atmospheric height, computed with continuum scattering (solid line) and with continuum and line scattering (dot-dashed line).
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{14210fc1.eps} \vspace*{3mm} \end{figure}$	Figure C.1: Dependence of the tabulated Gaunt factor integral $P(\Delta E_{\lambda },T)$ for collisions of electrons with neutral atoms in the van Regemorter (1962) formula on the transition energy $\Delta E_{\lambda }$ and the gas temperature T. The dashed lines mark the boundaries for typical values of $\Delta E_{\lambda }/kT$ found in solar surface convection simulations (Sect. 5.1).
Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.