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 |
Radiative transfer with scattering for domain-decomposed 3D MHD simulations of cool stellar atmospheres
Numerical methods and application to the quiet, non-magnetic, surface of a solar-type star
W. Hayek1,2,3 - M. Asplund1 - M. Carlsson3 - R. Trampedach4,2 - R. Collet1 - B. V. Gudiksen3 - V. H. Hansteen3 - J. Leenaarts5
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
.
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,
,
where
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
in direction
:
where


The extinction of photons is described, as customary, through the absorption coefficient
and the scattering coefficient
,
which combine to the gas opacity,
![]() |
(2) |
and give rise to the definition of the photon destruction probability
![]() |
(3) |
Recasting the optical path



with the source function

where


The source function
at optical depth
in direction
includes local thermal radiation from the gas and coherent scattering of photons:
where scattered radiation from direction







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
to the energy equation.
Evaluating the first moment of Eq. (1) and using the above definitions yields
where



Integrating the monochromatic flux divergence in Eq. (7) over the whole wavelength spectrum of the star yields the local heating rate
:
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



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:
The intensity-weighted mean opacity






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,
does not ensure a correct total radiative energy flux at optical depths
where radiative transfer is diffusive. It needs to be replaced by
the Rosseland mean opacity, defined as the weighted harmonic mean
We consequently use a








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
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
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
.
The discretization is performed by interpolating the source function for a given wavelength
(or bin number) along the ray using a second-order Bézier curve (see, e.g., the discussion in Auer 2003)
where





The shape of the three interpolation coefficients




![]() |
(15) |
where S0' is the centered derivative on the three-point stencil







The mean intensities J and the components of the flux vector are
computed by approximating the zeroth and first moment integrals by a
quadrature sum over selected ray angles (``method of discrete
ordinates''),
![]() |
(16) |
where wi is the weight of direction



Short characteristics require knowledge of the upwind intensities
for each ray direction
,
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
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 operator, which is commonly defined through

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 operator defined in Eq. (17), the problem may be rewritten into the matrix equation
with the identity matrix




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
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
![]() |
(19) |





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
(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
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:
![]() |
(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:
![]() |
(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

![]() |
Figure 1:
Horizontal mean photon destruction probability |
Open with DEXTER |
![]() |
Figure 2:
Convergence of the source function for bins 1-3 (see Fig. 1) during a simulation run without domain decomposition (left column), with
|
Open with DEXTER |
The convergence speed of an iterative method depends on the spectral radius
of the operator with which corrections are computed, as the error of the solution after n iterations decreases with
.
The spectral radius approaches
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, Sn, with an approximate solution
which we obtained after additional iterations with a lower convergence threshold of
,
assuming
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
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
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 .
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
than indicated by the relative corrections with respect to Sn-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
![]() |
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 , S and the upwind intensities
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 1003 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 302 mesh points; it is slanted with an angle of
off the vertical and an azimuth of
.
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 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
drastically increase the stiffness of the hydrodynamical equations,
where the stability-limited time steps of the transport and diffusion
terms scale with
and
,
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
,
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
is eventually needed, and therefore seems a natural choice. The expression
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
.
At the same time,
amplifies round-off errors of (J-S) beneath the optical surface, where the radiation field thermalizes (
,
also in the scattering case since
): 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 S0 from the discretized formal solution (Eq. (14)), which yields the modified integration constant
.
Using this equation, one obtains a monochromatic
along each ray. We note, however, that this leads to a deviation between the radiative energy
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
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:
![]() |
(22) |
where




Following Vögler et al. (2005), we compute radiative transfer on cell corners to improve the accuracy of
.
Radiative fluxes
are averaged over cell corners surrounding each face before computing difference quotients, while
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
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
,
since conservation of the radiative energy flux through the hydrodynamical cell surfaces must hold.
4 Absorption and scattering opacity sources in the Sun
![]() |
Figure 4:
Wavelength and depth dependence of the continuum photon destruction probabilities
|
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
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
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 (
).
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
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 H2 molecules is negligible. Rayleigh and electron scattering are treated as isotropic, neglecting their weak
anisotropy, where
is the scattering angle away from the incident direction (see, e.g., Mihalas 1978).
Between 5000 Å and 1.644 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
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
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
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
of spectral lines and the total photon destruction probabilities
,
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.
![]() |
Figure 5: Horizontal average
heating rate per unit mass around the stellar surface at an arbitrary
time step of the simulation. Boxes show
|
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
model covers a horizontal area of 6 Mm
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
% 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
1010 erg g-1 s-1, owing to the strong sensitivity of the heating rate per unit mass,
,
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.
![]() |
Figure 6:
Left: heating rates
|
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
10-6, where irradiation
from above is negligible. Rosseland optical depths at the bottom typically reach
107,
where radiative transfer is entirely diffusive and the radiation field
is completely thermalized. We therefore set the diffusion approximation
for all ingoing intensities at the bottom.
The three simulations discussed in Sect. 5.3 have mean effective temperatures
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
.
The simulations yield time-series of snapshots spanning
h of stellar time each, covering several granule lifetimes (
min) and several periods of the dominant p-mode (
min). Our simulation box covers about 10 granules with typical sizes of the order of
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
,
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
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
for
the given mean structure; the cooling is slightly stronger near
the surface. This behavior is expected from the mostly large photon
destruction probabilities
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
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
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
,
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
s. We consider time steps at
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
(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
-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
relation.
![]() |
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
(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
deviate
up to 40 K from the stratification without scattering, resulting
in a slightly steeper temperature gradient around the surface (
). 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
.
The temperature structure in the lower photosphere is thus hardly
affected by scattering. The opposite is the case in the high
photosphere and above (
), 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
![]() |
Figure 8:
Horizontal averages of the radiative heating rates
|
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
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
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
(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
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
![]() |
(23) |
where



![]() |
(24) |
where e is the internal energy per unit mass, P is the gas pressure, and

is the first law of thermodynamics. An upflowing (downflowing) gas parcel cools (heats) through expansion (compression) represented by the


![]() |
(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



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
and
(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
,
but fluctuates with lower amplitude. At
,
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.
![]() |
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
![]() |
(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

![]() |
Figure 10:
Horizontal and temporal averages of the relative temperature fluctuations
|
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
.
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.
AcknowledgementsThis 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
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
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
![]() |
(A.1) |
where





![]() |
(A.2) |
for the optical depths along the characteristic and choosing the control point

yields second-order interpolation of S, which now also depends on the source function


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

(cf. Eqs. (8) and (9) in Kunasz & Auer 1988). If the source functions




Overshoots are avoided by limiting





Note that, contrary to the first two cases, suppressing such overshoots leads to discontinuous left-hand and right-hand derivatives at S0.
Optical depths are computed in an analogue fashion to avoid negative results. A second-order Bézier polynomial
interpolates opacities over the path length
along the ray; integration over s yields the optical depth interval
![]() |
(A.3) |
where the control point


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)=at3+bt2+ct+d, | (B.1) |
with the curve parameter
![$t\in[0,1]$](/articles/aa/full_html/2010/09/aa14210-10/img201.png)
![]() |
(B.2) |
where the interpolation weights





with the grid spacing



![]() |
(B.3) |
with the weighting factor
![]() |
(B.4) |
which depends on the grid spacing




Quadratic interpolation uses only one of the two derivatives f'1 and f'2, depending on the interpolation parameter t. The interpolation coefficients are
![]() |
||
![]() |
||
![]() |
||
![]() |
||
![]() |
(B.5) |
Appendix C: Line scattering with the van Regemorter formula
![]() |
Figure C.1:
Dependence of the tabulated Gaunt factor integral
|
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
where




The photon destruction probability for a two-level atom is given by
where B21 is the rate for induced de-excitation, and




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
is given by the sum over all transitions,
![]() |
(C.4) |
Inserting Eq. (C.3), thus assuming the above mentioned approximations, yields
![]() |
(C.5) |
where the equality holds since





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 [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 [Google Scholar]
- Bautista, M. A. 1999, A&AS, 137, 529 [Google Scholar]
- Bell, K. L. 1980, J. Phys. B, 13, 1859 [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 [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 [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 [Google Scholar]
- Peach, G. 1970, MmRAS, 73, 1 [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 [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 [Google Scholar]
- Wishart, A. W. 1979, MNRAS, 187, 59P [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Table D.1: Continuum opacity sources.
All Figures
![]() |
Figure 1:
Horizontal mean photon destruction probability |
Open with DEXTER | |
In the text |
![]() |
Figure 2:
Convergence of the source function for bins 1-3 (see Fig. 1) during a simulation run without domain decomposition (left column), with
|
Open with DEXTER | |
In the text |
![]() |
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 |
![]() |
Figure 4:
Wavelength and depth dependence of the continuum photon destruction probabilities
|
Open with DEXTER | |
In the text |
![]() |
Figure 5: Horizontal average
heating rate per unit mass around the stellar surface at an arbitrary
time step of the simulation. Boxes show
|
Open with DEXTER | |
In the text |
![]() |
Figure 6:
Left: heating rates
|
Open with DEXTER | |
In the text |
![]() |
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 |
![]() |
Figure 8:
Horizontal averages of the radiative heating rates
|
Open with DEXTER | |
In the text |
![]() |
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 |
![]() |
Figure 10:
Horizontal and temporal averages of the relative temperature fluctuations
|
Open with DEXTER | |
In the text |
![]() |
Figure C.1:
Dependence of the tabulated Gaunt factor integral
|
Open with DEXTER | |
In the text |
Copyright ESO 2010
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.