A&A 401, 405-418 (2003)
DOI: 10.1051/0004-6361:20021853
J. Steinacker1 - Th. Henning2 - A. Bacmann3 - D. Semenov1
1 - Astrophysical Institute and University Observatory
(AIU),
University of Jena,
Schillergässchen 2-3, 07745 Jena, Germany
2 -
Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
3 -
European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching, Germany
Received 23 September 2002 / Accepted 9 December 2002
Abstract
We present the new grid-based code STEINRAY which has been developed
to solve the full 3D continuum radiative
transfer problem generally arising in the analysis of star-forming
regions, matter around evolved stars, starburst galaxies,
or tori around active galactic nuclei.
The program calculates the intensity emerging from these complicated
structures using a combination of step-size controlled ray-tracing and
adaptive multi-wavelength photon transport grids.
Along with a 2nd order finite differencing approach, the grids are
optimized to reduce the discretization error and provide
global error control.
The full wavelength-dependent problem is
solved without any flux approximation, and for arbitrary scattering
properties of the dust.
The program is designed to provide spatially resolved images, visibilities, and spectra
of complex dust distributions even without any symmetry for wavelengths ranging
from the UV to FIR and allows for
multiple internal and external sources.
In this paper, the algorithm is described and the capabilities of the
code are illustrated by the treatment of 1D and 3D test cases.
Analyzing the properties of typical cosmic dust, we
discuss the wavelength range for which
the time-consuming solution on adaptive grids can be omitted.
The temperature is calculated self-consistently using standard
accelerated -iteration.
Key words: radiative transfer - methods: numerical - accretion, accretion disks - ISM: dust extinction
Nearly all information about astrophysical sources comes from analyzing the radiation we receive from the objects. Therefore, radiative transport (thereafter RT) is one of the most fundamental processes considered in astrophysics whenever radiation is altered on its way from the source to the telescope. The alteration can range from small perturbations like the reddening of stellar light due to interstellar extinction or the change of polarization due to the Faraday effect, to almost complete shielding at short wavelengths and thermal re-emission at infrared and sub-millimeter wavelengths in the case of deeply-embedded star-forming regions, envelopes, and evolved stars, or tori around active galactic nuclei.
In view of the complex structure of the matter around these objects revealing filamentary, disk- or ring-like distributions, simple 1D approximations often fail to describe the observed images and spectral energy distributions, making multi-dimensional RT mandatory. But even for symmetric dust distributions, asymmetric source distributions or boundary conditions like, e.g., a single star near the dust cloud boundary have to be treated with 3D RT. This will be even more necessary in the near future, as the oncoming interferometers such as VLTI, Keck, and ALMA will reach resolutions revealing details at milli-arcsecond level, hence resolving complicated 3D structures like planetary gaps (Bryden et al. 1999; D'Angelo et al. 2002), density waves (Pfalzner et al. 2000), interaction streamers (Günther & Kley 2002), and disk warps (Bouvier et al. 1999).
But it is not only for explaining complex observational data that flexible three-dimensional RT programs are required. 3D Smooth-Particle Hydrodynamical (SPH) and Magneto-Hydrodynamical (MHD) simulations are now available and produce time-dependent density and temperature distributions of dust and gas in star-forming regions, proto-planetary disks, and in dust tori around AGNs. Without a 3D RT code, it is impossible to predict at which wavelengths or by which telescope these 3D features can be detected. Moreover, only with a 3D RT program can the validity of the approximate RT used in 3D SPH/MHD simulations be verified.
In line RT, the existing codes often use coarse approximations for the spatial resolution and the scattering to be able to handle the numerically intensive calculation of the level population (e.g. Folini & Walder 1999; Uitenbroek 2001; and references therein).
For 3D continuum RT, different methods have been applied to solve the complex transport equation, but only a few codes are available. The advantage of Monte-Carlo methods is the easy treatment of complex density distributions, complicated scattering functions, and polarization. Monte-Carlo methods require high numbers of photons to cover re-emission in all directions and to treat optically thick regions, which can be reduced using elaborated concepts like immediate reradiation. The major drawback, however, is that there is no global error control when using Monte-Carlo schemes. They have been applied for 3D configurations in the papers of Egan & Shipman (1995), Wolf et al. (1999), Wolf & Henning (1999) incorporating polarization, and Gordon et al. (2001), Misselt et al. (2001) considering also transient heating.
Integrating the RT equation along the considered line of sight by using a ray-tracer is a very flexible method for treating complex density configurations. It provides full error control, but requires the implicit re-calculation of the step size and becomes very time-consuming when the optical depth is high, so that approximations have to be used. The first solution was published by Yorke (1986) solving the wavelength-dependent problem in flux-limited approximation and without self-consistent temperature iteration. Well-established in hydrodynamical simulations, solutions on fixed grids can also be performed in RT using a finite difference or short characteristics discretization. However, resolution of all relevant structures for all wavelengths is very hard to achieve for most astrophysical applications. Given that the intensity depends on wavelength, three spatial and two directional coordinates, a decent resolution of 100 points in each variable already leads to solution vectors with 1012 components, which is not possible to store simultaneously even in supercomputers nowadays. Together with the complicated integro-differential type of the RT equation which makes application of common sparse matrix solvers impossible, a solution of the 3D RT can only be achieved by using sophisticated adaptive grids and fast solution algorithms. Two 3D continuum RT grid codes are available. Stenholm et al. (1991) solved the transport equation on a regular grid using 1st order finite differencing with low spatial resolution. Folini et al. (2003) have developed a code to solve the optically thick, stationary, wavelength decoupled, and unpolarized NLTE radiative transfer problem for moving media of a given density.
This publication is the second in a series presenting the new grid-based 3D RT code STEINRAY which incorporates the use of multi-wavelength adaptive photon transport grids to provide global error control of the solution. The underlying grid generation mechanism was described in Steinacker et al. (2002b). In this paper, the full discretization and solution strategy is given including the self-consistent temperature iteration, as well as simple test cases to show the capabilities of the code. In another paper, we compare the results for a standard accretion disk in a 2D benchmark project with three other RT codes (Pascucci et al. 2003). Further publications will deal with applications to circumbinary disks, warped disks, gaps in disks, pre-stellar cores, and dust tori around AGNs.
In Sect. 2, the continuum RT equation is described and the solution strategy is outlined. The optical properties of typical cosmic dust are briefly reviewed in Sect. 3 and limiting wavelengths are derived to decide when time-consuming calculations on the adaptive multi-wavelengths photon transport grids have to be performed. The solution of the transport equation in the case of substantial scattering is described in Sect. 4. Section 5 describes the discretization of the temperature distribution. The code is tested in Sect. 6 for simple 1D cases and a complex 3D application. We summarize our findings in Sect. 7.
We describe the radiation field by
the total specific intensity
,
where
gives the location in space,
is the direction of the radiation, and
its wavelength.
Starting with the boundary values, we can calculate
the transport of the radiation through the considered
domain by solving the stationary 3D RT equation
![]() |
(2) |
The physical quantities describing the efficiency to absorb and scatter the
incident radiation by an ensemble of dust grains can be written as
![]() |
(3) |
Two source terms of radiation are explicitely given in
Eq. (1), namely scattering into the line of
sight and re-radiation by the dust particles.
The scattering into the direction
is described by the probability
that radiation is scattered from the direction
into
,
with the solid angle
.
The second source is the re-radiation of the dust particles
.
Intensity and dust temperature are not independent. The radiation field
determines the temperature and in turn, the dust re-emission contributes to
the radiation field.
This couples the partial integro-differential RT equation to the
local energy balance equation describing how a dust particle is heated
by the source radiation and the radiation of all other particles.
The balance
equation for the energy density in local thermal equilibrium
at point
is
A simultaneous treatment of the 3D RT and hydrodynamical (HD) equations in one code is currently beyond the capabilities of nowadays computers. HD simulations commonly use an approximate RT to calculate radiative heating, while in turn, RT codes can use the derived densities and the heating sources to calculate the radiation field at a given time or in a stationary picture.
In view of the substantial computational effort to solve the
3D transport equation, it is mandatory to use any approximation
that is allowed by the physical conditions.
With vanishing scattering integral, (1) becomes a 1st order differential
equation which can be solved without problems using ray-tracing.
To make use of this approximation, and
since the operators in the integro-differential Eq. (1) are linear in ,
we can (following e.g. Efstathiou & Rowan-Robinson 1990)
split the total specific intensity into an unprocessed passing source
component I* that includes also the thermal contribution from
the dust, and a processed component I of radiation that has encountered scattering
As intended by the split (5), the first transfer Eq. (7) can be transformed to a path integral and thus easily be calculated using the formal exponential solution. For an empirical optical data set, the numerical solution is conveniently derived for all wavelengths e.g. using fifth-order Runge-Kutta with adaptive stepsize control (Press et al. 1992) as ray-tracing routine. This provides error control for the solution, and can be used for all optical depths from thin to thick regions. Moreover, it allows to treat multiple external and internal sources of radiation.
The second Eq. (8) still has an integro-differential
form and requires a separate, more sophisticated treatment.
The third Eq. (9) allows to update the temperature
from an intensity that has been calculated assuming a fixed temperature
using a standard accelerated -iteration (Ng 1974).
In order to decide where the scattering integral becomes negligible
and when the solution of the RT equation I* can be calculated
by solving (7) easily,
the optical properties of the dust grains have to be investigated.
For the sake of comparison, we concentrate on homogeneous spherical particles
of a uniform size (radius m) consisting of amorphous carbon
(Preibisch et al. 1993; Rouleau & Martin 1991)
or astro-silicate (Draine & Lee 1984; Draine 1985).
Their optical properties are calculated
by the Mie theory for spheres (Mie 1908) using the modified code of
Barber & Hill (1990).
More complicated dust models like mixtures of (in)homogeneous spherical particles of various sizes and compositions can be treated by the code without difficulty, but at expense of longer computational time.
To illustrate the typical behavior of the optical properties, we show in
Fig. 1
the ratio of the scattering and absorption efficiency factors
of dust particles consisting of amorphous carbon for different
sizes a as a function of wavelength.
![]() |
Figure 1: Ratio of scattering and absorption efficency factor as function of the wavelength. The different curves correspond to different particle sizes a. Two solid lines at the values 2 and -2 limit the regions where scattering and absorption is considered to dominate, respectively. |
Open with DEXTER |
![]() |
Figure 2: Limiting wavelength for which scattering can be neglected, plotted as a function of particle size a for amorphous carbon (dotted) and astro-silicate (solid). In the parts above the lines, absorption is more than a factor of 100 stronger than scattering. |
Open with DEXTER |
Depending on a given particle size a and wavelength ,
the program
switches from direct solution of the scattering
Eqs. (7)-(9) for
to the easier ray-tracing solution of (7) for
with the limiting wavelength
.
If scattering is important, beside the scattering cross section,
the phase function
describing the probability of
the incoming radiation to be scattered from the direction
into
also has to be specified. This is important, as for isotropic scattering,
when
,
(1) can be simplified.
Integrating over the solid angle, the scattering integral
disappears and the equation becomes an ordinary differential equation
![]() |
(13) |
![]() |
Figure 3:
Scattering phase function in polar coordinates for homogenous
spherical silicate particle of size ![]() ![]() ![]() |
Open with DEXTER |
![]() |
Figure 4:
Scattering phase function in polar coordinates for homogenous
spherical silicate particle of size
![]() ![]() ![]() |
Open with DEXTER |
Anisotropic scattering can become very important.
If the dust grains are, e.g., of size m,
in the optically thin case, there will be just a few scattering events
before the radiation escapes the object, and the scattered radiation
is strongly peaked in the forward direction.
Assuming isotropic scattering, the scattered radiation will be more
homogeneous and substantially different from the anisotropic case.
We just mention as one prominent example the atmosphere of an accretion
disk where the matter has low density
and just a few scattering
events occur before the radiation escapes the object.
In particular, interpreting images obtained for edge-on accretion disks
at short wavelengths is impossible without using the correct anisotropic
phase function.
In the program, we always use the correct anisotropic phase function.
If more than one dust size is considered, the phase function may differ substantially between the different sizes. For efficient mixing of the dust grains (assuming the same density distribution for all sizes), a mean scattering function may be defined to handle the transport of the scattered part of the radiation approximatively.
In the following, we describe how the integro-differential Eq. (8) is discretized and solved in the program.
For the direction coordinates ,
we choose spherical coordinates
,
while Cartesian coordinates x,y,z are best suited
for spatial coordinates in 3D RT problems with complex density distributions.
Hence, (8) turns into
For the spatial discretization, we use the adaptive optimized
multi-wavelength photon transport grids presented in Steinacker et al. (2002b).
These grids are obtained with a refinement grid generator
to minimize the 1st order discretization error and
guarantee global error control for the solution of the radiative
transfer problem on the grid. For each wavelength grid point, a spatial grid is
calculated, so that the number
of grid points per grid
depends on the wavelength. The three Cartesian coordinates are
xfs,yfs,zfs,
respectively, used at wavelength grid point f where
.
The numbering of the grid points by s is caused by the refinement
procedure which generates the adaptive grids. They will be renumbered with
respect to the boundary
conditions and the considered direction.
The intensity is therefore discretized by a solution vector with
components,
reaching values around 106 for typical astrophysical applications.
Each
vector component of the intensity can be written as
![]() |
(15) |
The notation for a use of a 2nd order
finite differencing scheme on an adaptively
refined grid is shown in Fig. 5a.
![]() |
Figure 5: Panel a) Sketch of a refined grid and the location of the seven point 2nd order finite differencing stencil centered at point A. In the upper part, the points are labeled according to their numbering with respect to point A along the x, y, and z-direction, respectively. Point B marks the center of a stencil moved downward with respect to A, with the point C indicating a stencil point which is not located on the grid. Panel b) Refined grid and two cell deep shadow grid for a direction that is indicated in the small coordinate system at the bottom. |
Open with DEXTER |
In 2nd order approximation, the derivatives with respect to the point
(xfs,yfs,zfs) and suppressing the s index
read as
Inserting (16)-(23) into the radiative
transfer Eq. (14), we obtain
![]() ![]() |
||
![]() ![]() |
(24) |
![]() |
![]() |
![]() |
(25) |
![]() |
![]() |
![]() |
(26) |
![]() |
![]() |
![]() |
(27) |
Vf,g0,0,0 | ![]() |
![]() |
(29) |
Uf,g0,0,0 | ![]() |
![]() |
|
![]() |
|||
![]() |
(30) |
The boundary values are defined on a shadow grid that is not part of the adaptive photon grid. The boundary intensities are determined by the interstellar radiation field and nearby sources. As we use a 2nd order differencing scheme, the shadow grid will be 2 cells deep. Figure 5b shows such a shadow grid for a ray direction with positive components indicated by the arrow in the bottom right corner. The cell size of the shadow grid will not influence the solution and is chosen to be 1/10 of the computational domain.
Since intensities of radiation propagating in all directions are
calculated, the stencil will also be moved through the grid
in other directions than
the one following the orientation of the axis.
Equation (28) is also valid also for the down-winding case
(one of the components of direction vector
is negative)
if you re-order the grid points with respect to the direction vector
,
as can be shown by considering
rf0<rf-1<rf-2 with
.
A major numerical concern in 3D calculations is to reduce
the amount of calculations within the innermost loop of the
RT solver as much as possible.
In this loop, (7) and (8) are solved
by -iteration.
We therefore pre-calculate the positions
of all stencils for all directions. There are eight different direction
types to be considered ranging from all three Cartesian components
of
being positive to all being negative.
Beside that, the standard refinement used in generating the grids
guarantees that we can switch from a floating point storage of
the grid points to integer values as multiples of the smallest
occurring grid cell size. This reduces storage space of the grids and
increases the speed of the solution substantially.
Note that the solution scheme (28) has been derived for
a Cartesian grid, but equivalent schemes could be derived for
other coordinate systems.
![]() |
Figure 6: Examples of temperature grids used by STEINRAY. Depending on the density distribution, equidistant Cartesian a), spherical grids with logarithmic radial points b), strongly overlapping nested Cartesian grids c), cylindrical grids with logarithmic height griding d), and adaptive grids with standard refinement e) are used. The number of grid points has been reduced to keep the figure readable. |
Open with DEXTER |
The overall solution strategy can be summarized as follows: we specify a density distribution of the dust, its size distribution, absorption and scattering properties, as well as the scattering phase function. We define the location and the spectrum of the radiation sources, and the boundary conditions. We specify the limiting wavelength for which we can neglect scattering.
To start the computation, we analyze the density
distribution and determine minima, maxima and strong gradients.
Using the concept of penetration depths (see Steinacker et al. 2002b),
we construct adaptive photon transport grids for each wavelength.
For the temperature, we create another grid depending on the symmetry
of the density and source distribution. On this grid, we determine
a start temperature by illuminating each cell by the sources and
then each cell by each cell.
With this temperature T0, the unprocessed intensity I* is calculated
by solving (7). If scattering is important, also I is determined
by solving (8). In the latter case, we start with the scattered intensity
from the formerly treated wavelength and calculate the source term (6).
Then, the intensity is updated solving (8) with this source term
and this -iteration is continued until the intensity change
drops below a given value
(typically we choose
).
We use the standard technique proposed in Ng (1974)
to accelerate convergence by using former solutions.
When
has been calculated, we update the temperature using (9),
and continue this outer
-iteration until the temperature has converged
to the self-consistent value.
During the iterations, we check the global energy to be conserved (and equal
to the start source energy).
It is possible to derive a priori guesses as to where the radiation
transport needs fine spatial grids to resolve the physics correctly
and hence to deduce the emerging radiation. Knowing the density
distribution and the optical properties of the dust particles,
the optical depth throughout a grid cell can be used
to derive a wavelength-dependent first order finite differencing criterion for the
grid generation (see Steinacker et al. 2002b). In most cases, this leads
to a refinement in the region around .
For the temperature, this is hardly possible, as the temperature couples
to all wavelengths via (9) and hence, a clear definition of the
layer is not evident.
In regions where the radiation field is dominated by the source
contribution, the temperatures are expected to be high with high gradients.
Considering e.g. an accretion disk, strong temperature gradients
at the inner boundary and in the disk atmosphere are related to strong
density gradients.
For shielded regions, radiation is dominated by re-emitted long-wavelength
dust emission and the temperature gradients should be small.
Therefore, adaptively refining grids similar to the density grids obtained in Steinacker et al. (2002b) could be used providing a global error control also for the temperature. However, the use of two different adaptive grids is very time-consuming. It has to be kept in mind that both the ray-tracing solver and the scattering-integral solver on the photon grids need the temperature at arbitrary points within the computational domain, and this very often. An 8-point interpolation is an acceptable effort, but with an adaptive grid, a tree-code search algorithm is needed to find the corresponding grid cell, which is too much of an effort especially for the scattering integral solver. In view of the expected smooth behavior (compared to the radiation field) of the temperature, simpler grids are used in STEINRAY with the drawback of having no a-priori error control for the temperature. Some of them are plotted in Fig. 6, namely an equidistant Cartesian grid (a), a spherical grid with radial points being logarithmically equidistant (b), strongly overlapping nested Cartesian grids (c), cylindrical grids with logarithmic height griding (d), and adaptive grids with standard refinement (e) for cases where the same grid for radiation and temperature can be used as the gradients within the object are smooth.
The outer temperature iteration requires to start with a given initial
temperature distribution
and then to improve
this distribution using (9)
to obtain the correct temperatures
.
As this outer iteration requires to run the entire solution process
several times, it is essential to have quick
convergence and a good initial guess.
In the Appendix, we derive a scheme to calculate the initial
temperature distribution T0 for a given density distribution and optical
dust properties.
T0 is obtained neglecting scattering but considering both
stellar-cell and cell-cell illumination.
For deeply embedded
sources, the transformation of short-wavelength photons to infrared
radiation takes place in a thin dust layer, and as it matters little
for the overall temperature if the radiation is scattered there a
few times, T-T0 will be small and convergence fast.
If one deals with a configuration that has an extended atmosphere, namely
a region that can be seen from outside with an optical depth of the order
of 1 for wavelengths where scattering is important (like e.g. an accretion
disk), T0 and T will differ substantially in visible regions and thus
the images calculated from them will deviate too.
As we have pointed out in Sect. 3, forward scattering dominates the
scattering processes at wavelengths of strong scattering.
An improved T1 can be calculated using this effect:
in the source illumination of each cell, only the absorption
instead of the complete extinction is considered, arguing that the
radiation is not nearly lost but scattered forward in first approximation.
Testing RT codes is difficult as even in simple 1D cases, for realistic density distributions, the RT equation already gets too complicated to be solved analytically. Benchmark studies based on well-defined RT problems can aid here to compare the results of different codes. It may be pointed out here explicitly, though, that there may be situations where the majority of codes fails, e.g. due to resolution problems of standard algorithms. Only by incorporating error control is it possible to infer reliability of the technique. Hence, error control in STEINRAY will be one of the issues we will discuss in the following.
For 1D continuum RT, a benchmark was published by Ivezic et al. (1997). They compared results of three codes with similar algorithms for an application in spherical geometry. For 2D continuum RT, Pascucci et al. (2003) have defined a benchmark based on simulations of a standard accretion disk around a solar-type star. The project involved four separately developed codes, two based on the Monte-Carlo technique and two based on grids, each with very different implementation. STEINRAY has been tested in this collaboration in good overall agreement with the other three codes.
There is no benchmark project for 3D continuum radiative transfer yet. Hence, we present first results of the code two-folded: we show simple test cases where the solution is either exactly or approximately known, and a more complex 3D scenario to show the capabilities. For applications to circumbinary disks, warped disks, disks with planets, and AGN dust tori, as well for a 3D benchmark test including images, we refer to the later papers of this series.
In order to compare with an analytical solution, we start with the
most simple case possible, namely radiation passing through a
homogeneous dust medium with a constant density and no dust emission
radiation.
We follow the radiation up to a distance of 1000 AU from the source that is emitting
radiation just along the x-axis. For this so-called slab geometry, we assume
a constant dust number density of
1 m-3, spherical silicate dust particles of size
m, and a wavelength of
m so that scattering
can be neglected.
The optical depth
![]() |
(31) |
![]() |
Figure 7: Relative error of the intensity in a homogenous illuminated cloud slab without continuum re-emission of the dust. The straight line shows the error when using a constant stepsize tracer, the noisy line gives the error of a 5th-order Runge Kutta-solver. |
Open with DEXTER |
In spherical geometry, the intensity will spatially depend on the
radius only, and we can use the scaling argument by Ivezic et al. (1995)
to reduce the amount of free parameters in the 1D RT problem by
normalizing the occurring parameters.
We use the benchmark provided by
Ivezic et al. (1997) to test the code in spherical geometry. From their
models, we choose the scaled example of a central source with T*=2500 K,
a temperature at the inner boundary of T1=800 K, an outer
radius normalized to the inner boundary of Y=1000,
density power law index p=0 (constant density),
and optical depth at the
reference wavelength
m to be
.
For the sake of comparison, we also use isotropic scattering.
We show the adaptive photon transport grid for the scattered radiation
at a wavelength of
m in Fig. 8.
![]() |
Figure 8:
Adaptive photon transport grid for the scattered part of
the radiation passing through a dust envelope with spherically symmetric
density distribution, for the wavelength
![]() |
Open with DEXTER |
The upper panel of
Fig. 9 shows the temperature we obtained (solid line)
along the z-axis as a function
of the normalized radius and the benchmark results (crosses).
![]() |
Figure 9:
Comparison with the 1D benchmark results of Ivezic et al. (1997)
for T1=800 K, Y=1000, p=0, ![]() ![]() ![]() ![]() ![]() |
Open with DEXTER |
Tests with a 2D RT benchmark have been given in Pascucci et al. (2003) comparing the results obtained with STEINRAY for a 2D standard accretion disk with three other 2D and 3D RT codes. They show very good overall agreement in both the spectrum and the temperature distribution. We therefore proceed to full 3D RT. But here, no benchmark has been defined sofar. And moreover, there is no benchmark yet that gives images for multi-dimensional continuum RT. So, to show the capabilities of the code, we just give here images of a complex test case and refer to further applications and benchmark runs in forthcoming papers.
The example we choose is a spherically symmetric cloud core that contains
two young stars having formed at slightly different times, close to the
core center. We assume the stars to be massive (
),
with a luminosity of
,
and a surface temperature of T=20 000 K.
The emerging radiation is assumed to have already formed cavities free of dust, and all
remnants like disks that may have formed the stars are destroyed. We choose a slight time
shift in the evolution to have cavities with different sizes 1000 AU and
2000 AU, respectively, and at distances of 1500 and 2000 AU from the core center.
The stars are moved away from the cavity centers
towards the core center to distances of 1500AU and 1000 AU, respectively. This is to take into account the density gradient
within the core and its influence on the cavity formation.
The density within the core is assumed to
be a Gaussian of the form
with
m-3and r0=2000 AU.
Figure 10 illustrates the density distribution in the z=0-plane.
![]() |
Figure 10: Cut through the density distribution of the 3D test case at z=0. The two massive stars are located in cavities that have formed within a dense core with a density maximum in-between the two stars, their position is indicated by peaks. |
Open with DEXTER |
For
m, the configuration is optically thin and the
dust emission at the surface of both
cavity shells accumulates so that they are clearly visible at all
viewing angles.
The re-emission function
![]() |
(32) |
![]() |
Figure 11: Temperature distribution of the cavity configuration in the x-y-plane. a) Single star located on the x-axis at x=-1500 AU. b) Two stars located on the x-axis at x=1000 AU and x=-1500 AU. |
Open with DEXTER |
![]() |
Figure 12:
Images of the cavity configuration at
![]() ![]() ![]() ![]() ![]() |
Open with DEXTER |
At wavelength of
m, the inner part of the core is
optically thick.
Part of the emission from the dust around the larger cavity is obscured
by the dense core.
The same is valid for the image at
m. As the absorption
is stronger here, the smaller cavity is deeply embedded and emission is
only seen from parts where the density of the dense core drop due to the
Gaussian profile. In both pictures, emission from the star is too weak to
be seen.
In the left upper panel of Fig. 12, the stars can be seen at
m (view along the y-axis, stars at x=1000 AU and -1500 AU,
total extend 10 000 AU). Still it is the dust close to the surface of the
cavity that is dominating the emission, but only the emission in the outer
parts can be seen as the extinction drops with the density.
The emission from the cavity surfaces is not entirely smooth, showing a slight numerically caused ring pattern. This is an effect of the finite grid size of 100 AU of the temperature grid. Doubling the number of grid points in each direction (100 to 200) would increase the size of the temperature array by a factor of 8 and the computational time by a factor of 64. The use of an adaptive temperature grid reduces the finite grid cell errors but increases the runtime as the interpolation between the adaptive photon transport grid for the scattered radiation and the temperature grid becomes more time-consuming.
The numerically by far most difficult part is to produce images in the
optical and UV. We have chosen
m, as the scattering
is maximal in this wavelength range. With a deeply embedded pair of
stars like in our configuration, hardly any scattered light can be
expected to leave the core before being absorbed. As the images showed
just the reddened stars, we decided to lower the density by a factor
of 50 to make escape of scattered light possible.
We show two images at this wavelength in Fig. 13
with an inclination of the view by
with respect to the x-axis,
![]() |
Figure 13:
Images of the cavity configuration at
![]() ![]() |
Open with DEXTER |
The right image is calculated for the same viewing angle and wavelength,
but with a realistic scattering function as described in Sect. 3.
At
m, forward scattering dominates.
Therefore, only the front part of the surface of the larger cavity that is close to the star
and that has substantial density can scatter the radiation to the observer.
Additionally there is some diffuse scattered radiation coming from the front border
of the dense center of the core caused by the other star.
This comparison shows that a wrong approximation of the scattering will alter the resulting images substantially. Animations varying the viewing angle to support visualization of the 3D structure can be found under http://www.astro. uni-jena.de/Users/stein/Ani/anin.htm.
Both ray-tracing and finite differencing on the transport grids
allow error control of the solution. Nevertheless, we point out that in contrast
to the intensity calculations, for the
temperature iteration, it is impossible to give an a priori error limit.
Here, we rely on usual a posteriori methods, namely to check the global
energy conservation when integrating the out-coming flux over a closed surface
around the source of energy. The main sources of numerical errors are the
interpolations of the temperatures for the ray-tracing and the transport grids.
As they have to be carried out often, a quadratic or higher polynomial interpolation is
too time-consuming. The same is valid for temperature grids with cell numbers
exceeding 106. For the simulations, global energy conservation within 5%
was used as convergence criterion for the outer -iteration.
In this paper, we presented the solution algorithm implemented
in the grid-based code STEINRAY, designed to
solve the 3D continuum radiative transfer problem for the intensity emerging
from objects in star-forming regions, evolved stars with envelopes,
starburst galaxies, and AGNs.
The method is a combination of step-size controlled ray-tracing and
solution on adaptive multi-wavelength photon transport grids, on which the
finite-differencing discretization error is minimized.
We briefly analyzed the optical properties of typical cosmic dust grains,
and discussed the wavelength range for which
the time-consuming solution on adaptive grids has to be used.
For the temperature, we have presented and discussed possible grids on which the
temperature distribution can be calculated self-consistently.
Aside from 2D benchmark comparisons presented elsewhere, we have tested the
code with simple 1D cases and illustrated the capabilities by treating a complex 3D test
case. The temperature is calculated self-consistently using standard accelerated
-iteration.
Acknowledgements
J. S. thanks F. Evans for valuable comments in the course of improving the adaptive photon transport grids. We thank D. Folini for an excellent referee report helping to improve the paper at various points.
Computer time can be saved substantially when starting the solution of the RT equation with a temperature that is already close to the correct value. The temperature iteration in the RT solution scheme will converge quickly when an initial temperature is used that is derived for the case that the high-energetic source radiation has been absorbed and re-distributed by the dust particles already. Here, we derive a scheme to calculate the initial temperature for a given dust configuration. The main approximation is that the influence of scattering is neglected for the heating process. In view of the strong forward scattering in the wavelength regions where scattering plays a role, this assumption is reasonable to derive an initial temperature distribution. We concentrate here on stellar illumination, but a corresponding scheme can be derived for other radiation sources like evolved stars or the central engine of an active galactic nucleus.
The power emitted by a star with the luminosity
that is received by the dust particle k
of size a with the absorption efficiency
at a distance rk through vacuum is (Evans 1994)
![]() |
(33) |
![]() |
(34) |
![]() |
(35) |
![]() |
(36) |
In the presence of extincting matter between star and dust,
the intensity along a ray
from the star
to cell K is damped according to
![]() |
(37) |
![]() |
= | ![]() |
|
![]() |
![]() |
(38) |
![]() |
(39) |
![]() |
(40) |
![]() |
(41) |
![]() |
= | ![]() |
|
![]() |
|||
![]() |
(42) |
![]() |
(43) |
![]() | |||
![]() |
(44) |
Numerically,
and
are obtained from ray-tracing.
For large numbers of cells, ray-tracing should be truncated when the
intensity has dropped below a given limit characterizing a small
fraction of the thermal energy in the cell. In the optically thick
case, only cells located within a few mean free paths of the photons
will contribute to the cell-cell illumination.
We calculate
![]() |
(45) |
![]() |
(46) |
![]() ![]() | (47) | ||
![]() |