Issue |
A&A
Volume 511, February 2010
|
|
---|---|---|
Article Number | A81 | |
Number of page(s) | 16 | |
Section | Numerical methods and codes | |
DOI | https://doi.org/10.1051/0004-6361/200912355 | |
Published online | 16 March 2010 |
Fast and accurate frequency-dependent radiation transport for hydrodynamics simulations in massive star formation
R. Kuiper1,3 - H. Klahr1 - C. Dullemond1 - W. Kley2 - T. Henning1
1 - Max-Planck-Institut für Astronomie, Abt. Planeten- und
Sternentstehung,
Königstuhl 17, 69117 Heidelberg, Germany
2 - Universität Tübingen, Institut für Astronomie und Astrophysik,
Abt. Computational Physics, Auf der Morgenstelle 10, 72076
Tübingen, Germany
3 - Fellow of the International Max-Planck Research School for
Astronomy and Cosmic Physics at the University of Heidelberg
(IMPRS-HD), Germany
Received 20 April 2009 / Accepted 15 December 2009
Abstract
Context. Radiative feedback plays a crucial role in
the formation of massive stars. The implementation of a fast and
accurate description of the proceeding thermodynamics in pre-stellar
cores and evolving accretion disks is therefore a main effort in
current hydrodynamics simulations.
Aims. We introduce our newly implemented
three-dimensional frequency dependent radiation transport algorithm for
hydrodynamics simulations of spatial configurations with a dominant
central source.
Methods. The module combines the advantage of the
speed of an approximate flux limited diffusion (FLD) solver in the
one-temperature approach, which is valid in the static diffusion limit,
with the high accuracy of a frequency dependent first order ray-tracing
routine. The ray-tracing routine especially compensates the introduced
inaccuracies by standard approximate FLD solvers in transition regions
from optically thin to thick and yields the correct optical depths for
the frequency dependent stellar irradiation. Both components of our
module make use of realistic tabulated dust opacities. The module is
parallelized for distributed memory machines based on the message
passing interface standard. We implemented the module in the
three-dimensional high-order magneto-hydrodynamics code Pluto.
Results. We prove the viability of the scheme in a
standard radiation benchmark test compared to a full frequency
dependent Monte-Carlo based radiative transfer code. The setup includes
a central star, a circumstellar flared disk, as well as an envelope.
The test is performed for different optical depths. Considering the
frequency dependence of the stellar irradiation, the temperature
distributions can be described precisely in the optically thin, thick,
and irradiated transition regions. Resulting radiative forces onto dust
grains are reproduced with high accuracy. The achievable parallel
speedup of the method imposes no restriction on further radiative
(magneto-) hydrodynamics simulations.
Conclusions. The proposed approximate radiation
transport method enables frequency dependent radiation hydrodynamics
studies of the evolution of pre-stellar cores and circumstellar
accretion disks around an evolving massive star in a highly efficient
and accurate manner.
Key words: radiative transfer - hydrodynamics - stars: formation - circumstellar matter - accretion, accretion disks - methods: numerical
1 Introduction
Massive stars are of importance for a wide range of astrophysical problems. In spite of the fact that they are both rare and short-living, they represent the major source of radiation energy in their stellar clusters. Therefore, they act as valuable tracers of star formation rates in distant galaxies. Additionally, the radiation emitted during a massive star's lifetime influences the surroundings through various interactions such as heating and ionization of gas, evaporation of dust, and radiative forces leading to powerful stellar winds and outflows. Finally, when a massive star dies, it enriches its neighborhood with heavy elements. In this sense, massive stars are the main drivers of the morphological, dynamical, and chemical evolution of their complex environments (Zinnecker & Yorke 2007; McKee & Ostriker 2007; Henning 2008).
However, knowledge about their formation is rather poor compared to the case of low-mass star formation. Observationally, this is mostly due to their larger average distances and the fact that they are deeply embedded in dense, opaque cores especially during their early evolutionary phases. Also, their short lifetime, rareness, and complex environment pose difficulties for detailed observations. Nevertheless, current observational results, e.g. from the Spitzer Space Telescope Survey and millimeter interferometry, support the assumption that the basic concepts of star formation, including the collapse of an unstable gas and dust core, the formation of both bipolar outflows and jets, as well as accretion disks, are possibly valid throughout the whole range of stellar masses (e.g. Bally 2008; Vehoff et al. 2008; Steinacker et al. 2008; Bik et al. 2008; Kumar & Grave 2008). A good review of observational results can be found for instance in Beuther et al. (2007). Future generations of space telescopes and interferometric systems such as the Herschel Space Observatory, the James Webb Space Telescope (JWST), and the Atacama Large Millimeter Array (ALMA) will provide a deeper insight into the mechanisms of in-falling envelopes, bipolar outflows, and massive accretion disks. This will definitely put tighter constraints on current theoretical models.
From the theoretical point of view,
assuming that the formation of high-mass stars is basically a scaled-up
version of low-mass star formation, implies its own challenges:
Their rapid evolution, especially the shorter Kelvin-Helmholtz
contraction timescale, leads to the interaction of the accreting flow
of gas and
dust with the emitted radiation from the newborn
star (Shu et al. 1987).
In the extreme case of a purely spherically symmetric collapse of a
pre-stellar core (without angular momentum), the in-fall is potentially
stopped by the growing radiation pressure entirely.
In the static limit, the equilibrium between radiative and
gravitational forces at the
dust condensation radius is described by the Eddington limit ,
where L*, M*,
and
denote the stellar luminosity,
mass, and the dust opacity regarding the stellar irradiation spectrum,
respectively.
The underlying one-dimensional collapse problem is extensively studied
in analytical studies
(Shu 1977;
Hunter
1977; Penston 1969b)
as well as in numerical simulations (Foster & Chevalier 1993;
Whitworth
& Summers 1985; Hunter 1977; Penston
1969a), either in hydrodynamics or as evolution of
self-similar solutions. Also stellar radiative feedback has been
considered in one-dimensional studies (Larson & Starrfield 1971;
Kahn
1974; Wolfire
& Cassinelli 1987; Yorke & Krügel 1977).
However star formation is rarely a perfectly spherically symmetric
process.
As stated by Nakano (1989), massive star formation is at least a two-dimensional non-spherically symmetric problem, including initial angular momentum to force the formation of a circumstellar disk and polar cavities. Both will potentially help to overcome the Eddington barrier. On small scales, Turner et al. (2007) show that the so-called photon bubble instability can significantly enhance accretion and avoid the radiation pressure problem. On large scales, the formation of bipolar cavities possibly reduces the radiative force onto the accretion flow by the so-called flashlight effect (Krumholz et al. 2005b; Yorke & Sonnhalter 2002). Although Krumholz et al. (2009,2007a,2005a) claim that the Eddington barrier may be broken via a Rayleigh-Taylor instability, it is unclear if this instability is the most important one and how this instability will change in case of realistic frequency dependent radiative feedback instead of the used gray (i.e. frequency averaged) Flux Limited Diffusion (hereafter called FLD) approximation.
The important role of frequency dependent radiation transport in the formation of massive stars was already shown in Yorke & Sonnhalter (2002), but due to the huge computational overhead of the frequency dependent radiation transfer it was neither possible to study a huge amount of different initial conditions (to scan the parameter space), nor to perform high-resolution simulations of the accretion process.
Regarding these issues, we demonstrate here the usability of our newly developed fast three-dimensional frequency dependent approximate radiation transport module for numerical hydrodynamics.
The most accurate description of the physics proceeding during the collapse of a pre-stellar core would include a frequency dependent radiation transport step after each hydrodynamic timestep using a modern Monte-Carlo based or ray-tracing radiative transfer method. The CPU time needed to solve one single hydrodynamic timestep is generally orders of magnitude lower than the time for a frequency dependent radiation transport step, especially in complex geometries. Thus, this approach is not feasible with current computing technology for a large number of grid cells in more than one dimension due to the large amount of computational time needed for each radiation transport step. A more desired radiation hydrodynamics scheme should roughly spend the same CPU time on the radiation physics as on the hydrodynamics part. Sensible approximations are thus necessary to speed up the radiative transfer in such hydrodynamics studies.
In hydrostatic disk atmosphere computations, Murray et al. (1994) introduced a splitting of a two-dimensional radiation field into an irradiated and a diffuse part. Already Wolfire & Cassinelli (1986,1987) used such a splitting in one dimension to study the radiation feedback of massive stars, which has been shown to be a valid approach in one-dimensional hydrodynamics simulations by Edgar & Clarke (2003). In the present paper, we expand the method to higher dimensions and show its validity in an axially symmetric setup, consisting of a central star, a flared circumstellar disk, and an envelope. In fact, we found that it is necessary to consider the frequency dependence of the stellar irradiation feedback to reconstruct a reasonable approximation to the spatial temperature distribution.
A first approach for fast two-dimensional axially symmetric radiative transfer is, for example, the gray diffusion approximation studied by Tscharnuter & Boss (1993), which is only applicable in the optically thick limit. Ray-tracing based methods (like in Efstathiou & Rowan-Robinson 1990) show high accuracy, but as already mentioned also require much CPU time, which yields low efficiency in combination with hydrodynamics solvers. Another common approach for the description of radiative processes in circumstellar disks is the gray FLD approximation (Kley 1989; Klahr & Kley 2006; Bodenheimer et al. 1990; Klahr et al. 1999). It provides a fast method to determine the temperature evolution in the optically thick (diffusion) as well as in the optically thin (free-floating) limit, but shows stronger deviations in transition regions, demonstrated e.g. in Boley et al. (2007). In the case of gray FLD, which is still the default technique in modern radiation hydrodynamics codes, this method of course suffers strongly from the lack of frequency dependence, when compared to the accuracy of modern Monte-Carlo based codes.
In the following, we present our results in combining the advantage of the gray FLD approximation (speed) with the accuracy of frequency dependent ray-tracing. The approximation described here results in a large reduction in computing time compared to Monte-Carlo based radiative transfer. This allows us to implement this particular code in the framework of three-dimensional (magneto-) hydrodynamics simulations of circumstellar disks and in-falling envelopes on a parallel decomposed (spherical) grid. In particular, we implemented our radiation transport module into the three-dimensional high-order Godunov magneto-hydrodynamics code Pluto (Mignone et al. 2007). This tool allows us to perform high-resolution studies of accretion flows around a forming massive star.
The basic equations of the approximate radiation transport method and the numerical details are described in the following Sect. 2. The test of the accuracy of the proposed radiation transport scheme in the setup of a currently standard radiation benchmark test (Pascucci et al. 2004) is presented in Sect. 3. The parallel scaling behavior of our implementation is discussed in Sect. 4. In Sect. 5, we show the results of two standard radiation hydrodynamics shock tests performed with the FLD part of our implemented radiation transport method and the open source magneto-hydrodynamics code Pluto (Mignone et al. 2007). Conclusions as well as further studies planned are described in Sect. 6.
2 Theory and numerics of the approximate radiation transport scheme
In the following, we recapitulate the general ideas, the basic equations, and methods of the frequency dependent approximate radiation transport scheme. This section should allow the reader to follow our motivation for the newly implemented three-dimensional radiative transfer module. Every implemented formula and numerical detail is given.
The general idea of the method is to split the radiation transport into two components (Edgar & Clarke 2003; Wolfire & Cassinelli 1986; Murray et al. 1994). Stellar radiative forces will mostly act on the surrounding gas and dust mixture in the first transition region from optically thin (e.g. where the dust is evaporated) to optically thick (e.g. a massive accretion disk). This is exactly the region, where the FLD approximation is incorrect (e.g. Boley et al. 2007). To avoid this disadvantage of the FLD approximation, we first calculate the stellar radiative flux through the environment including its absorption in a corresponding first order ray-tracing routine. First order means that the spatial distribution of the radiative flux from stellar irradiation is calculated according to its frequency dependent optical depth, but re-emission of photons is shifted to a gray FLD solver. Sources for the thermal dust emission are the absorbed energy from the prior stellar irradiation step and potentially additional heating from the hydrodynamics, e.g. due to compression of the gas, accretion luminosity of sink cells, or viscous heating. In other words, this means that instead of solving the whole radiation transfer problem either in the FLD approximation or with a ray-tracing technique, we just extract the first (most important) absorption event of the stellar irradiation from the FLD solver and calculate the appropriate flux in an accurate ray-tracing manner. This splitting method allows us to consider the frequency dependence of the stellar spectrum in a very efficient manner.
In the first Sect. 2.1, we recapitulate the FLD equation for thermal dust emission. In the following two Sects. 2.2 and 2.3, we explain how this FLD solver can simply be combined with a first order ray-tracing routine (either gray or frequency dependent respectively) to include irradiation feedback from a single central object. Due to the fact that these rays are aligned with the radial coordinate axis in spherical coordinates, this kind of coordinate system is highly favored, but not required, to solve the ray-tracing step. In the last Sect. 2.4, we comment on the so-called generalized minimal residual method (GMRES), our default implicit solver algorithm for the FLD equation.
2.1 Flux limited diffusion
The thermal dust emission is solved in the FLD approximation based on a
diffusion equation for the thermal radiation energy .
Within a given spatial density
and temperature
distribution, we start from the time evolution of the internal energy
density
and thermal radiation energy density
:
with the corresponding thermal pressure P, dynamical velocity





We add up the Eqs. (1)
and (2) and solve
the transport term
separately, if necessary, in the dynamical problem. If we consider the
transport of the internal energy
already during the corresponding hydrodynamics step (operator
splitting), the remaining terms yield
where the source term

In the following, we use the assumption that the gas and
radiation temperature are in equilibrium (also called one-temperature
radiation transport), which is e.g. a widely used approach in radiation
hydrodynamics simulations of circumstellar disks. In a purely FLD
radiation transport method, this assumption is justified for optically
thick regions (e.g. deeply inside an accretion disk). The
usage of our split radiation
scheme guarantees also the correct gas temperature in regions of
dominating stellar irradiation
(like an optically thin envelope or disk atmosphere). In residual
regions, shielded from the irradiation behind the optically thick
circumstellar disk, the gas and radiation temperature will be the same
in systems where the flow time is long compared to the
gas-radiation equilibrium time only. The flow time is ,
where l and u are the
characteristic length and velocity scales of the problem. The
gas-energy equilibrium time is of the order of
,
since the energy exchange rate is
.
Thus, the requirement that
reduces to the requirement that
,
where
is the optical depth of the system. This is the regime that Mihalas & Mihalas (1984)
identifies as static diffusion. Thus,
the scheme is viable in the static diffusion limit. This is consistent
with our test results for the radiative shock problem presented in
Sect. 5.
In this special case of a strong shock in a shielded, optically thin
region,
the one-temperature approach yields the correct temperature of the
shocked gas, but results in a steeper temperature gradient than a
two-temperature radiation transport method in the
upstream direction of the radiation flux.
Consideration of such small
effects at outer regions of the domain goes beyond the scope of this
research project, which focuses on the details of stellar radiative
feedback on the accretion flow onto a massive star. Additionally, Krumholz et al. (2007b)
show that massive star formation is generally in the static diffusion
limit.
In the end, we will benefit from the speedup of the radiation transport
solver by cutting the numbers of unknown variables in half. Moreover,
the stiffness of the set of equations is much less if the above local
equilibrium assumption is made.
Expressing both energies on the left hand side of
Eq. (3)
in terms of temperature allows us to
derive a relation between the time derivatives of the energies. The
radiation energy density in absence of irradiation (see following
Sect. 2.2
for the case including irradiation) is given by
This expression plus the internal energy density


With this relation, the Eq. (3) reduces to a standard diffusion equation
with


with the diffusion constant




Summing up, we can describe the thermal/diffuse radiative
processes via:
with
2.2 Irradiation
When including stellar irradiation, some of the equations have to be
modified.
The irradiation from a central object is treated as an additional flux ,
released along rays in the radially outward direction.
is the energy density of the radiation emitted by the circumstellar
material, and does not include the stellar irradiation. The additional
radiation power
from
this irradiation is therefore added to the diffusion Eq. (8) as a source term
and has to be added to the left hand side of Eq. (4), which is used to calculate either the thermal radiation density or the corresponding dust temperature T under the assumption that the dust is in equilibrium with the combined stellar and diffuse radiation field:
Here

After solving the diffusion Eq. (15) implicitly (with
the radiation energy
as the solution vector), the corresponding dust temperature is
calculated using Eq. (16).
Since the opacity
depends on the temperature the right hand side of Eq. (16)
also depends on the solution of T. This
makes an iterative procedure based on the Newton-Raphson method
necessary to find the solution.
The ratio of opacities
corresponds to the relation of emission from and absorption by dust.
The stellar radiative flux as a function of distance r
from the central object is calculated by
The (boundary) flux at the stellar surface

with the Stefan-Boltzmann constant


In the case of the radiation benchmark test by Pascucci et al. (2004), the optical depth between the stellar surface and the inner disk (the inner boundary of the spherical computational domain) is assumed to be negligible.
A flow chart of the radiation module described so far is shown for a static problem such as the benchmark test of Sect. 3 in Fig. 1.
![]() |
Figure 1: Schematic flow chart of the radiation module. Exponents declare the iteration (timestep) number n. The actual timestep corresponds to d tn = tn - tn-1. |
Open with DEXTER |
2.3 Frequency dependent irradiation
Taking into account the frequency dependence of the stellar flux, we
consider a fixed number of frequency bins,
characterized by their mid-frequency ,
instead of the gray irradiation. For the radiation benchmark test we
use the opacity tables of Draine
& Lee (1984), including 61 frequency bins,
see Fig. 2.
![]() |
Figure 2:
Regarding the frequency dependence of stellar irradiation feedback.
Frequency dependent opacities |
Open with DEXTER |




The boundary condition, given by Eq. (18), is now
determined for each frequency bin by the integral over the
corresponding part of the stellar black body Planck
function
in the frequency dependent interval:
2.4 The generalized minimal residual method
The FLD Eqs. (8)
or (15)
are solved in our newly implemented approximate radiation transport
module for hydrodynamics in an implicit fashion via a generalized
minimal residual (GMRES) solver. The GMRES solver is a Krylov
subspace (KSP) method for solving a system of linear equations
via
an approximate inversion of the large but sparse matrix A
and is an advancement of the minimal residual (MinRes) method. The
GMRES method was developed in 1986 and is described in Saad & Schultz (1986).
This method is much better than the conjugate gradient (CG) or
successive over-relaxation (SOR) method, and at least as good as the
well-known improved stabilized Bi-Conjugate Gradient (BiCGstab) method
used, e.g., in Yorke &
Sonnhalter (2002).
The general idea of minimal residual solvers based on the
KSP method is the following: The ith KSP
is defined as .
In each solver iteration i the GMRES method
increments the subspace Ki
used with an additional basis vector
and approximates the solution of the system of linear equations by the
vector
which minimizes the norm of the residual
.
This method converges monotonically and theoretically reaches the exact
solution after performing as many iterations as the column number of
the matrix A (which equals the number of
grid cells). Of course, in practice the iteration is already stopped
after reaching a specified relative or absolute tolerance of the
residual, which normally takes only a small number of iterations. The
computation of each iteration grows like
.
In the current implementation, we use the so-called ``GMRES restarted''
by default.
GMRES restarted never performs all the iterations to reach the exact
solution.
After a priorily fixed number of internal iterations n,
the solver starts a second time in the first subspace K1
but with the last approximate solution
.
Due to the growing computational effort with
this approach generally results in a speedup of the
computation.
Table 1: Overview of performed Monte-Carlo runs for comparison.
Table 2: Overview of runs using the proposed approximate radiation transport module.
The radiation transport module is parallelized for distributed memory machines, using the message passing interface (MPI) language. The results of a detailed parallel performance test of the whole radiation transport module, including this GMRES solver is presented in Sect. 4.
3 Frequency dependent test of the approximate radiation transport
The approximate radiation transport introduced in the previous section can now be tested for realistic dust opacities in a standard benchmark test for irradiated circumstellar disk models. The setup of the following comparison is adopted from Pascucci et al. (2004) and includes a central solar-type star, an irradiated circumstellar flared disk and an envelope. We have to choose a low-mass central star, because no benchmark for high-mass stars was performed so far. However, the tests should not depend on the actual size and luminosity of the central star. For comparison, we choose a standard full frequency dependent Monte-Carlo based radiation transport code. The comparison is done for two different (low and high) optical depths taken from the original radiation benchmark test. To test each of the components (gray and frequency dependent irradiation as well as Flux Limited Diffusion) of the proposed approximate radiative transfer separately, we perform several test runs with and without the different physical processes (absorption and re-emission) with both the Monte-Carlo based (see Table 1) and the approximate (see Table 2) radiative transfer code.
3.1 Setup
3.1.1 Physical setup of the star, the disk, and the envelope
The stellar parameters are solar-like: The effective temperature of the
star is 5800 K and
the stellar radius is fixed to 1 solar radius
.
The disk ranges from
AU up to
AU.
Although the numerical setup of the gas density is done in
spherical coordinates, the analytic setup of the gas density, as
described in the original benchmark test, is given in cylindrical
coordinates:
![]() |
(21) |
with the radially and vertically dependent functions
![]() |
(22) |
and
![]() |
(23) |
making use of the following abbreviations
![]() |
= | ![]() |
(24) |
![]() |
= | ![]() |
(25) |
![]() |
= | ![]() |
(26) |
The lowest density is limited to a relative factor of 10-100 compared to the highest density (at



![]() |
![]() |
![]() ![]() |
0.1 |
![]() |
![]() |
100 |
![]() |
![]() |
The opacity tables used are the same as in the original benchmark test (Pascucci et al. 2004) taken from Draine & Lee (1984). They are displayed in Fig. 2.
3.1.2 Numerical setup of the approximate radiation transfer code
The runs of the approximate radiative transfer code are performed on a
radially stretched, polar uniform, spherical, two-dimensional grid. The
grid consists of 60 cells in the radial direction times 61
cells in the polar direction (plus additional cells for the storage of
boundary conditions). The polar range covers the full spatial setup of .
Stretching of the radial grid dimension by an additional 10% from one cell to the next is applied. The implicit diffusion Eq. (15) is solved via the GMRES method (see Sect. 2.4) after parallel/global Block-Jacobian and serial/local ILU pre-conditioning in the framework of the version 2.3.3 of the open source parallel solver library PETSc (Portable, Extensible Toolkit for Scientific Computation). More detailed information about this solver library can be found in Balay et al. (2001,2004).
The gradient of the radiation energy is zero at the inner
radial and both polar boundaries ,
i.e. radiative flux over these boundaries is prohibited. The radial
outer boundary is defined as a constant Dirichlet boundary
corresponding to
.
The timestep used for the FLD solver is 104 s. The main iteration (circle of irradiation and FLD steps) is stopped when the relative change of the temperature in each cell is smaller than 0.01%, leading to 3 main iterations in the ``purely absorption'' runs and more than 600 main iterations in the ``irradiation plus FLD'' runs.
3.1.3 Numerical setup of the Monte-Carlo based comparison code RADMC
For comparison we use the Monte-Carlo based radiation transfer code
RADMC described in Dullemond
& Turolla (2000) and Dullemond & Dominik
(2004).
The general solver method of the code is based on
Bjorkman & Wood
(2001). The Monte-Carlo runs are performed on a
grid assuming symmetry to the disk midplane. The grid is stretched in
both directions. One million photons are used.
Scattering can be handled by this code, but is simply switched off as it is neglected in the radiation transport module described here. Scattering would increase the temperature in the irradiated parts (up to an optical depth of about unity) by about 2% in the optically thin envelope up to a maximum of 19% in the optically thick inner rim of the disks midplane due to higher extinction. The more effectively shielded outer regions of the disk would be about 10% cooler. For more massive and luminous stars the effect of scattering would decrease due to stronger forward scattering, which is included in our ray-tracing routine per definition.
3.2 Configurations of runs performed
The following comparison of the results is divided into three parts
(each of the proposed components of the theoretical Sects. 2.1-2.3
are tested independently): First, we study a pure absorption
scenario without diffusion in the optically thin and thick case.
Afterwards we include diffusion effects.
Therefore, we perform three different runs of the Monte-Carlo based
code:
a full run for the optically thin case ,
a full run for the optically thick case
,
and an additional run for the optically thick case
with excluded re-emission of the photons (to achieve a pure absorption
scenario for comparison with our first order ray-tracing routine),
hereafter called ``one-photon-limit''.
An overview of these three comparison runs is given in Table 1.
The resulting temperature distributions of these Monte-Carlo runs are compared with the results of the approximate radiative transfer runs including different components of our module: We discuss five different configurations for the optically thin and thick setup including gray and frequency dependent irradiation plus potential diffusion. An overview of the physics applied and the resulting deviations of these runs from the comparison data is given in Table 2. These results are discussed and illustrated in detail in the following subsections.
3.3 Results
3.3.1 Gray absorption in an optically thin disk
In the most optically thin case ,
diffusion effects should be negligible. Therefore, we can test the
validity of the routines described in Sect. 2.2
without running the diffusion routine. Compared to the full Monte-Carlo
simulation from RADMC also the deviation in the most ``difficult''
region of the midplane (due to having the highest absorption) stays
below 2% (see Fig. 3).
In this optically thin limit, the diffusion effects are indeed negligible: An additionally performed full run of the approximate radiation transfer module with frequency dependent irradiation plus Flux Limited Diffusion shows a variation in the radial temperature slope from the pure gray irradiation run below 1%.
![]() |
Figure 3:
Radial cut through the temperature profile in the midplane in the most
optically thin case |
Open with DEXTER |
3.3.2 Gray and frequency dependent absorption in an optically thick disk
In the most optically thick case ,
we test a pure absorption case to distinguish the deviations introduced
by the FLD approximation in the full run from
the deviations introduced by the irradiation component of our module.
![]() |
Figure 4:
Radial cut through the temperature profile in the midplane in the most
optically thick case |
Open with DEXTER |
3.3.3 Gray and frequency dependent absorption plus diffusion in an optically thick disk (the complete problem)
In this and the following section, we show the comparison results
between the complete approximate radiation transport method with the
combination of (gray or frequency dependent) irradiation and Flux
Limited Diffusion and the corresponding Monte-Carlo simulation in the
most optically thick case .
![]() |
Figure 5:
Radial cut through the temperature profile in the midplane in the most
optically thick case |
Open with DEXTER |
![]() |
Figure 6:
Polar cut through the temperature profile at r=2 AU
of the frequency dependent irradiation plus Flux Limited Diffusion run
for the most optically thick case |
Open with DEXTER |
Also in this case, we found that the frequency dependent irradiation is
necessary to achieve
the accuracy needed for realistic radiative feedback in hydrodynamics
studies (e.g. in massive star formation or in irradiated accretion
disks). Gray irradiation combined with the FLD approximation leads to
deviations up to 38.4% in the resulting temperature slope (see
Fig. 5).
Including frequency dependent irradiation and FLD the radial
temperature profiles agree within 0.7% at the inner boundary up to a
maximum deviation of 11.1% at roughly .
The comparison between the complete radiation scheme and the
corresponding Monte-Carlo run is shown in radial and polar temperature
profiles in Figs. 5
and 6
respectively. The turnover point (at a polar angle from the midplane of
)
from the optically thin envelope to the optically thick disk region is
reproduced very well.
Further remarks and discussion of these results are given in the last Sect. 3.4 at the end of the comparison section.
3.3.4 Radiative force
Stellar radiative feedback on the dynamics of the environment plays a crucial role in the formation of massive stars. The heating will probably prevent further fragmentation of the cloud by enhancing the Jeans mass (e.g. Krumholz et al. 2007a). Furthermore, the dusty environment feels the radiative force when absorbing the radiation due to momentum conservation, which potentially stops the accretion process for highly luminous massive stars. Therefore, it is necessary to compute the correct radiative flux (and its derivative, the radiative force) in addition to the temperature distribution in such simulations.
The radiative force density is, accordingly to
Mihalas & Mihalas
(1984), given by
In our split radiation transport the total radiative flux




In discretized space the opacity and density are constant over
a single grid cell and the irradiated stellar radiative flux
is calculated at the cell interfaces. To calculate the mean radiative
force inside this cell it is necessary to integrate the radiative force
over the cell volume (e.g. a simple ansatz of averaging only
the stored fluxes at the interfaces towards the cell center would lead
to unphysically high radiative forces for
,
i.e. for the case that most of the flux is absorbed on a length scale
much smaller than the grid size). Integrating the above formula
Eq. (27)
over space (for simplicity here: an one-dimensional cartesian grid with
an uniform grid spacing of
)
leads to:
The flux at position x inside the grid cell is given by the flux

with the optical depth

Finally, the mean radiative force is therefore given by the difference of the left and the right flux into and/or out of the cell respectively
which equals in continuous space the derivative of the radiative flux
Combining Eqs. (27) and (33) reflects the fact that without emission the radiative flux is given by the differential equation
In other words, the relationship in Eq. (34) shows that the original Eq. (27) and the derived Eq. (33) are indeed identical (in continuous physical space), whereas the left hand side of Eq. (34) gives the expression easily available in discretized space. In our code we use a grid in spherical coordinates, thus the derivative of the stellar irradiation is given by
So the purely geometrical dilution of the flux in the radially outward direction

![]() |
Figure 7:
Radial cut through the radiative force profile at a polar angle of |
Open with DEXTER |
We compute the radiative force in the setup of the radiation benchmark test with gray and frequency dependent irradiation as well as the thermal radiative force for the most optically thick case and compare our results with the corresponding Monte-Carlo based run. The result is visualized in Fig. 7. The peak position is reproduced very well. The absolute value of the peak is underestimated by 3-4% only. Behind the absorption peak the radiative force smoothly drops down and the relative deviations grow radially outwards, but stay below 10%. The fraction of the radiative force, resulting from the thermal flux, is relatively small and is most important at and directly after the absorption peak (where most of the thermal radiation is emitted). This fraction will presumably be higher in denser environments. At larger radii (r > 300 AU) where the disk becomes highly optically thin for its own thermal radiation, the radiative force resulting from the thermal flux is negligible (cp. Figs. 8 and 7).
The gray approximation and the corresponding frequency
dependent run show only small deviations (<
at the outer part of the disk for radii roughly larger than
200 AU (after the absorption of the stellar irradiation). With
the setup of Pascucci
et al. (2004) the gray approximation leads to a
higher radiative force than the frequency dependent ones, but in
general the difference of both methods depends strongly on the stellar
luminosity. The black body spectrum of more luminous stars will shift
to higher frequencies (cp. Fig. 2).
Further remarks, explanations, and detailed discussion of the resulting radiative force and temperature profiles are given in the following Sect. 3.4.
3.4 Remarks and analysis
3.4.1 The Monte-Carlo comparison code
For comparison and interpretation of the results, we should mention
that in the original radiation benchmark test
(Pascucci et al. 2004)
the different Monte-Carlo codes themselves differ in the radial
temperature profile through the midplane of the circumstellar disk in
the most optically thick case
by 5% in most of the region between 1.2 to 200 AU and up to
15% towards the outer border of the computational domain in the radial
direction. This means that the deviations in the optically thin case
(Fig. 3)
as well as the frequency dependent ray-tracing
part of the optically thick case (Fig. 4) stay beneath
the discrepancy of the different Monte-Carlo
solutions. As expected, the direct stellar irradiation is determined
highly accurately, when considering the frequency dependence.
Therefore, the errors introduced by using the FLD approximation can be
limited in the test throughout the irradiated regions.
The influence of the so-called photon noise in the Monte-Carlo
method is illustrated in the highly optically thin regions
(
from the midplane) in Fig. 6, where the
temperature should actually be independent for large polar angles (as
displayed by the solid line).
3.4.2 The Flux Limited Diffusion approximation
A special feature of the setup of the original benchmark test (Pascucci et al. 2004)
is the fact that even in the most optically thick case ,
which is defined for a wavelength of 550 nm, the pre-described disk is
locally optically thin for the radiation from thermal dust emission.
Integrating the corresponding local optical depth
from the outer edge of the disk through the midplane towards the center
yields a final optical depth of
,
see Fig. 8.
![]() |
Figure 8:
Radial profile of the optical depth |
Open with DEXTER |
The plot shows clearly the low optical depth for the thermal component
of the radiation field especially in the outer part of the disk, which
results in an underestimation of the temperature in the transition
region at roughly 200 AU
due to the overestimation of the radiative flux in the outward
direction by the FLD approximation.
The FLD approximation is known to be valid in the most
optically thin (free-streaming limit) as well as in the most optically
thick (diffusion limit) regions only. The apparent yet surprisingly
good agreement between the Monte-Carlo based runs and the herein
proposed radiation transport module also in the intermediate region of
the flared disk atmosphere (see Fig. 6) is due to
the newly implemented direct irradiation routine which yields the
correct flux and depth of penetration for the different frequency bins
of the stellar irradiation
spectrum (see Sect. 3.3.2).
The slight underestimation of the temperature at
(see Fig. 5)
in the disk midplane in the most optically thick case is most likely a
result of an intermediate region (transition from optically thick to
optically thin) in the outward radial direction, which is shielded from
the direct irradiation and is not in good agreement with the FLD
approximation.
3.4.3 The frequency dependence
Approximating the frequency dependence of the stellar spectrum by gray
(frequency averaged) Planck mean opacities results in an
overestimation of the optical depth in the infrared part
and an underestimation of the absorption in the ultraviolet part of the
stellar spectrum (see Fig. 2). This leads to
an under- or overestimation of the radiative force onto dust grains at
the first absorption peak depending on the stellar luminosity due to
less absorption of the most energetic UV photons or stronger absorption
of the photons at the peak of the stellar black body spectrum.
Secondly, the gray approximation results in an underestimation of the
resulting temperature deeply inside the disk due to absorption of the
infrared photons already at inner disk radii.
The quantitative amount of the resulting deviations depends in general
on the specific stellar and dust properties. Due to the steep decay of
the stellar black body spectrum at high frequencies, the difference of
the gray and frequency dependent radiative force turns out to be very
small in this specific setup. On the other hand,
consideration of the frequency dependence is essential to limit the
deviations in the resulting temperature profile to
less than 11.1% (compared to 38.4% for gray irradiation plus
FLD, see Fig. 5).
These issues are well illustrated in Fig. 2.
The figure shows the frequency dependent opacities of
Draine & Lee (1984),
the approximated frequency averaged value of the Planck mean opacity
regarding the stellar effective surface temperature
as well as the black body spectrum of the central star (to visualize
the amount of radiative flux which is emitted per frequency bin). Each
dot in the figure marks the mid-frequency of the correspondingly chosen
frequency bin. The effects on the resulting temperature profile and
radiative force cannot be generalized easily.
They depend on the underlying dust model and strongly on the properties
of the central star, which yield a shift of the peak position of the
black body spectrum in Fig. 2 according to
Wien's
displacement law .
In the specific setup of Pascucci
et al. (2004) the Planck mean opacity at the black
body peak position is higher than the frequency dependent ones, leading
to a slightly higher radiative force. The strong overestimation of the
opacity in the infrared regime leads to the huge discrepancy of the
gray approximation
in the radial temperature profile through the disk.
4 Parallel performance of the approximate radiation transport module
The parallelization of the radiation transport scheme and the GMRES solver (see Sect. 2.4 for details of how the solver works) are taken care of by the PETSc library (Portable, Extensible Toolkit for Scientific computation), see also Balay et al. (1997). To test the parallel speedup of the implemented radiation transport module we perform two tests with an extended version of the circumstellar disk setup introduced in Sect. 3.1.
We adopt the most optically thick setup for
and expand it to three dimensions assuming axial symmetry. All runs
include frequency dependent irradiation and Flux Limited Diffusion.
The tests run for 10 main iterations, which consume the main
computational effort for the approximate solver (later on, near
equilibrium, the internal iterations needed decrease strongly). The
number of internal iterations of the approximate
implicit solver is fixed to 100 to guarantee the same amount
of
computation in all runs during this benchmark test.
Due to the parallel Block-Jacobian pre-conditioner, the number of
needed internal iterations (for a specified accuracy) normally
increases with increasing number of processors. The precise value for
the increase is strongly problem dependent (Smith, member of the
development team of the PETSc library, private communication).
The parallel domain (the linear system of equations) is split in the azimuthal and polar direction only, which insures good speedup and efficiency. Decomposing the domain in the radial direction would decrease the parallel performance due to the fact that in the ray-tracing routine it is necessary to compute and therefore communicate the flux from the central sink cell to the outer boundary from the inside outward. Since the knowledge of the flux at the inner cell interface is needed to compute the flux at the outer interface, this method is hardly parallelizable as a domain decomposition.
The measured times t2 to tn (n is the number of processors) represent the wall clock time per main iteration per grid cell without the non-recurring initialization and finalization of the code. We perform runs with 2 up to 64 processors due to the fact that single job submission is not available on the cluster we use. Cases, in which the local cache size would exceed the parallel decomposed problem size, are not taken into account, i.e. no misleading super linear speedup for high number of used processors is shown here. The speedup S is calculated as the ratio of the ``serial'' run time compared to the wall clock time used by the parallel run: S=t2 / tn. The efficiency E is determined via E = t2 / (tn n) = S/n. Each run for a specific grid and a specific number of processors is performed three times and averaged afterwards, but the differences of the resulting run times are negligible.
All tests are performed on a 64-Bit Opteron cluster consisting of 80 nodes with two CPUs each.
4.1 The constant grid test
The runs during this test are performed on a grid consisting of grid
cells. Each processor covers therefore a
subdomain, depending on the number of processors used. This test shows
the speedup one can gain when running a fixed problem on more and more
processors. Therefore, the parallel efficiency declines stronger than
in the following growing grid test due to the fact that with the usage
of an increasing number of processors one lower the amount of
computation and
increase the amount of communication per single CPU. This means that
the granularity (ratio of computation to communication) of the parallel
problem drops
strongly with an increasing number of processors.
The resulting speedup factors and efficiencies of this test are shown
in Fig. 9.
![]() |
Figure 9:
Speedup factors S=t2
/ tn (
upper panel) and corresponding efficiencies |
Open with DEXTER |
4.2 The growing grid test
The runs during this test are performed on a grid consisting of 32,768*n
grid cells (e.g. a
grid), where n is again the number of processors
used. Each processor covers therefore a subdomain containing
32 768 grid cells respectively during all runs,
independent of the number of processors used.
This test shows the efficiency achievable when using more processors to
run a bigger problem.
This is a realistic setup for our three-dimensional radiative
hydrodynamics studies of collapsing massive
pre-stellar cores, which are at the current limit of our clusters. The
resulting speedup factors and efficiencies are shown in Fig. 9.
4.3 Parallel performance results
The constant grid test shows a clear speedup for the fixed problem, so it is simply possible to compute a fixed problem faster by using more processors. The growing grid case shows a high efficiency of more than 95% during all runs. This speedup seems to be higher than the speedup of any actual freely available three-dimensional hydrodynamics code for spherical grids (well known in the astrophysical community, e.g. Pluto2, Flash2.5, and Zeus-MP1.5), which we tested during the summer 2007 in an accretion disk setup in hydrostatic equilibrium. In this sense, the resulting parallel speedup is high enough for the planned integration of this radiation transfer module into a (magneto-) hydrodynamics framework.
Finally, we can recommend the usage of such an implicit modern KSP solver method leading to fast convergence of the radiative diffusion problem (at least for the setup discussed) while simultaneously offering high parallel efficiency. Admittedly, attention should be paid to the general fact that such an approximate solver method strongly depends on the physical problem at hand as well as on the specified accuracy or abort criterion.
5 Radiative hydrodynamics shock tests
To test the approximate radiation transport scheme also in dynamical interaction with a streaming fluid, we perform two standard radiative shock tests. We adopt the setup of the supercritical and subcritical radiative shock tests for the VISPHOT code in Ensman (1994). These radiative shock tests were already repeated in tests of the TITAN code (Sincell et al. 1999a,b), ZEUS-2D (Turner & Stone 2001), as well as ZEUS-MP (Hayes et al. 2006). Analytic approximations for this kind of problem are given by Zel'Dovich & Raizer (1967).
![]() |
Figure 10: Radiative supercritical shock. Resulting density, pressure, velocity, and temperature distributions at four different snapshots in time. Dashed lines represent the adiabatic runs, solid lines the radiative ones. The time snapshots are taken ( from left to right) at 860 s, 4000 s, 7500 s, and 13 000 s after launching. Mostly horizontal lines at the lower border of the graphics refer to the initial setup. The snapshots at 4000 s are additionally marked by open circles for every 10th grid cell to illustrate the resolution used. The spatial axes are retranslated into the non-moving frame used in the visualization by Ensman (1994) for the sake of comparison. The filled dots in the four profiles at 13 000 s represent characteristic sample points in the two-temperature approach read off Figs. 10-12 in the original study by Ensman (1994). |
Open with DEXTER |
To follow the motion of the gas, we solve (additionally to the
radiation transport equations of Sect. 2) the
equations of compressible hydrodynamics in conservative form
with the acceleration source term from the radiation transport
![]() |
(39) |
derived in Sect. 3.3.4. The evolution of the gas density


The test setup describes a piston moving with supersonic
velocity through an initially uniform, cold gas. The one-dimensional
domain of the grid associated covers a distance of length .
The iso-density in the domain is fixed to
.
For testing purposes, the gas is set to be completely ionized,
thus the mean molecular weight is
.
An ideal equation of state is used with
.
The opacity is fixed to a constant value of
.
The initial temperature drops linearly from 85 K at the
starting position
of the piston to 10 K at the outer boundary.
The velocity
of the piston is used to determine the strength of the shock.
While the piston moves through the domain the radiative energy
from the shocked gas will stream upwards leading to a pre-heating, a
pre-acceleration, as well as a pre-compression of the gas directly in
front of the shock.
If the temperature in this pre-heated region stays below the
temperature of the shocked gas, it is called a subcritical radiative
shock.
If the temperature in the pre-heated region equals the temperature of
the shocked gas, it is called a supercritical radiative shock. The
smallest piston velocity leading to a supercritical shock defines the
critical velocity
.
Ensman (1994) used a Lagrangian grid moving with the piston velocity. This setup is translated into an Eulerian grid by setting the initial velocity in the whole domain as well as the permanent velocity at the outer boundary to the negative of the piston velocity u0, compressing the gas at the inner reflective boundary, which represents the moving piston. In the visualization of our results (see Figs. 10 and 11) the spatial axes are retranslated into the non-moving frame used in the visualization by Ensman (1994) for the sake of comparison. The spherical coordinate system used at large radii (to achieve a planar geometry) by Ensman (1994) is translated into cartesian coordinates. We use 512 uniform grid cells to cover the spatial extent of the grid. These grid adjustments were also used in the test of the ZEUS-MP code (Hayes et al. 2006).
In this radiative shock test setup we are able to check the dynamical behaviour of the radiation module and the hydrodynamics. On the other hand, it implies a much easier treatment of radiation transport (due to the fact that the optical depth in front of the shock is practically constant) than the prior static but frequency dependent benchmark test of Pascucci et al. (2004). Therefore, only the FLD routine (as described in Sect. 2.1) is needed to perform the problem runs.
We study this shock scenario in purely adiabatic as well as radiative hydrodynamics simulations.
5.1 Radiative supercritical shock
![]() |
Figure 11: Radiative subcritical shock. Resulting density, pressure, velocity, and temperature distributions at five different snapshots in time. Dashed lines represent the adiabatic runs, solid lines the radiative ones. The time snapshots are taken ( from left to right) at 350 s, 5400 s, 17 000 s, 28 000 s, and 38 000 s after launching. The snapshot at 350 s mostly refers to the initial setup. The snapshots at 5400 s are additionally marked by open circles for every 10th grid cell to illustrate the resolution used. The spatial axes are retranslated into the non-moving frame used in the visualization by Ensman (1994) for the sake of comparison. The filled dots in the temperature profile at 38 000 s represent characteristic sample points in the two-temperature approach read off Fig. 8 in the original study by Ensman (1994). |
Open with DEXTER |

We choose the piston velocity to be the one used originally by
Ensman (1994), .
In some of the subsequent tests mentioned above a slightly lower piston
velocity
for
the supercritical radiative shock was used by the authors.
Figure 10 displays the resulting density, pressure, velocity, and temperature distributions for four different times (same as Figs. 10-12 in Ensman (1994)) as well as the initial setup. The analytic limit of a maximum jump in density by four in adiabatic shocks is reproduced. The effects of pre-heating, pre-acceleration, and pre-compression are clearly visible. The position of the peak temperature and the sharp declines in the density, the pressure, and the velocity at the location of the shock front fit the one by Ensman (1994). That means, the propagation speed of the shock is reproduced correctly. Although the non-equilibrium gas temperature spike, the maximum temperature, cannot be reproduced with an equilibrium radiation transport method (the width of this gas temperature spike would always be less than a mean free path across), the value of the peak temperature is in good agreement with the radiation temperature of the corresponding non-equilibrium run by Ensman (1994) aside from minor geometrical effects of the modified coordinate system as already discussed in Hayes et al. (2006). Furthermore, the equilibrium temperature distribution resembles the temperature distribution of the non-equilibrium run by Ensman (1994) in regions, where the gas temperature equals the radiation temperature (which is the most part of the domain in such a supercritical shock). As expected, the temperature profile calculated by our equilibrium diffusion method sharply declines at the boundary to the non-equilibrium region predicted by Ensman (1994). The outermost data points by Ensman (1994) in the pressure, the velocity, and the temperature distributions in Fig. 10 display the long-range tail of the pre-heated region, if using a non-equilibrium radiation transport method. Summing up, the deviations of the equilibrium temperature method from a non-equilibrium approach stay below 10% in the pre- and post-shock regions. In this dynamic diffusion shock test, the outermost tail of the pre-heated region, where the gas and the radiation temperature is not in equilibrium, can naturally not be observed in an equilibrium approach. In the context of massive star formation, Krumholz et al. (2007b) already show that massive star formation is generally in the static diffusion limit.
5.2 Radiative subcritical shock
In a radiative subcritical shock (
)
the gas in front of the shock region is heated to a temperature lower
than the temperature of the shocked gas region only.
In a non-equilibrium radiation transport method the gas temperature
sharply
drops down (after the narrow temperature spike similar to the
supercritical shock), whereas the radiation energy penetrates the cold
unshocked gas region and declines more smoothly (see also Fig. 8 in Ensman 1994).
An equilibrium radiation transport is expected to first ignore the narrow gas temperature spike. Secondly, it drops down sharply when arriving the region, where the gas and radiation temperature are out of equilibrium. In other words, the equilibrium radiation transport resembles the decline of the gas temperature distribution.
We choose the piston velocity to be the one used originally by
Ensman (1994), .
Figure 11
displays the resulting density,
pressure, velocity, and temperature distributions for five different
times (same as Fig. 8
in Ensman 1994).
Given the constraints of the equilibrium radiation transport discussed
above
the resulting distributions fully satisfy the predictions.
6 Conclusions and outlook
We introduced a fast and accurate radiation transport module, feasible for prospective usage in hydrodynamics and magneto-hydrodynamics simulations of circumstellar environments. The basic physical principle of the proposed radiation transport scheme is the splitting of the full radiative transfer problem into an accurate first order ray-tracing component and an approximate treatment of second order effects (Wolfire & Cassinelli 1986; Edgar & Clarke 2003; Murray et al. 1994). The method combines thereby the advantages of the Flux Limited Diffusion approximation (speed) and a first order frequency dependent ray-tracing routine (accuracy). The thermal radiation by dust grains is solved in the static diffusion limit via a FLD solver. The FLD approximation yields the correct solution in the optically thin (free-streaming) and optically thick (diffusion) limit. But in the transition regions, where also the most important feedback due to stellar radiative forces occurs, the FLD approximation is not valid anymore. This first impact of stellar irradiation is therefore accurately solved in the ray-tracing routine.
Our implementation of the FLD method uses a MPI parallelized implicit GMRES restarted solver yielding high speedup factors comparable to modern hydrodynamics solvers. The ray-tracing routine, which calculates the distribution of stellar irradiation, benefits from the usage of spherical coordinates, in which the outgoing stellar rays are aligned with the radial coordinate axis. Expanding this first order ray-tracing to frequency dependent irradiation is therefore convenient and results in higher accuracy of the temperature slope deeply inside the disk due to lower optical depth of the infrared part of the stellar spectrum. A minor flaw of the proposed radiative transfer module is the neglecting of scattering. But in the radiation benchmark test scattering has only a marginal influence on the resulting temperature distribution: Scattering would increase the temperature in the irradiated parts (up to an optical depth of about unity) due to higher extinction from about 2% in the optically thin envelope up to a maximum of 19% in the optically thick inner rim of the disks midplane. The more effectively shielded outer regions of the disk would be about 10% cooler. Secondly, for more massive and luminous stars the effect of scattering will decrease due to stronger forward scattering of the most energetic UV photons, which is automatically included in our frequency dependent first order ray-tracing routine. Rowan-Robinson (1980) stated that neglecting the effect of scattering except for the first scattering of the starlight yield a good approximation. Another flaw is the frequency averaged instead of frequency dependent re-emission of once absorbed stellar photons, i.e. the expansion of the ray-tracing routine to second and higher order would result in a more and more realistic (but also computationally more expensive) radiative transfer.
We tested the approximate radiation transport scheme for
hydrodynamics within the setup of a well established standard benchmark
test for (Monte-Carlo and finite differences based) radiative transfer
codes (Pascucci
et al. 2004), consisting of a central star, a
circumstellar flared disk, and an envelope. The deviations of the
resulting temperature slopes stay below a maximum of 11.1% compared to
a frequency dependent
Monte-Carlo based radiation transport method, which requires much more
CPU time. This deviation in temperature will be reflected in the
associated hydrodynamics in roughly 6% deviation in the
corresponding sound speed .
Neglecting the frequency dependence of the stellar irradiation however
results in significant higher
deviations up to 38.4% in the temperature distribution. Also the
computed (frequency dependent) radiative force onto
dust grains is in good agreement (below 3% at the inner rim and below
5% throughout the first absorption peak up to 10% relative deviation at
the outer boundary, where the radiative force gets negligible) with the
comparison radiation transport method likely due to the usage of direct
stellar ray-tracing feedback. The overall resulting deviations
potentially passed to the hydrodynamics solver are indeed smaller than
the lack of knowledge of the appropriate opacities. Summing up, the
accuracy of the proposed radiation transport method legitimates the
approximations needed for adequate speed in
combination with hydrodynamics.
In shock simulations, we show that the one-temperature approach, which is valid in the static diffusion limit, is sufficiently accurate for the astrophysical applications we have in mind. Krumholz et al. (2007b) already show that massive star formation is generally in the static diffusion limit. Nevertheless an extension to a two-temperature approach is straight forward, if desired for other astrophysical applications in the future.
Furthermore, we performed parallel performance tests of the radiative transfer module with different grid and parallelization setups. The parallel speedup of the solver imposes no restriction on large scale radiative (magneto-) hydrodynamics simulations.
Therefore, we conclude that frequency dependent radiative feedback in hydrodynamics studies of the formation of a massive star becomes feasible with the herein proposed approximate radiation transport method. Hydrodynamics simulations of monolithic pre-stellar core collapse models in massive star formation, including direct frequency dependent stellar feedback, is now one of our primary goals for this radiation transport module. Secondly, three-dimensional ideal and non-ideal MHD simulations of the magneto-rotational instability (MRI) in proto-planetary accretion disks with radiative cooling are intended (Kuiper et al. 2009). Thirdly, the formation of circumstellar disks and the evolution of their morphology in terms of potential planet formation can be studied also around low-mass stars. The algorithm is fully implemented in our version of the magneto-hydrodynamics code Pluto (including high-order Godunov methods for hydrodynamics, full implementation of the viscosity stress tensor, and additionally external forces, see Mignone et al. 2007).
AcknowledgementsThis research has been supported by the International Max-Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). Special thanks goes to Fabian Janssen for benchmarking different hydrodynamics codes on the PIA cluster during summer 2007. R. Kuiper would like to thank Gerhard Hoffmann for his permanent helpful suggestions and the discussions about numerical details, especially the parallelization, of the developed radiation transport module as well as Jürgen Steinacker for helpful conversation about the radiation benchmark test of Pascucci et al. (2004). Author R. Kuiper further thanks A. Mignone for the main development of the open source magneto-hydrodynamics code Pluto as well as his great support in code specific questions. Authors H. Klahr, C. P. Dullemond and W. Kley have been supported in part by the Deutsche Forschungsgemeinschaft DFG through grant DFG Forschergruppe 759 ``The Formation of Planets: The Critical First Growth Phase''. Parallel performance tests have been performed on the rio cluster of the Max Planck computing center in Garching.
References
- Balay, S., Gropp, W. D., McInnes, L. C., et al. 1997, Birkhäuser Press, 163 [Google Scholar]
- Balay, S., Buschelman, K., Gropp, W. D., et al. 2001, PETSc Web page, http://www.mcs.anl.gov/petsc [Google Scholar]
- Balay, S., Buschelman, K., Eijkhout, V., et al. 2004, PETSc Users Manual, Tech. rep., PETSc, revision 2.1.5 [Google Scholar]
- Bally, J. 2008, in Massive Star Formation: Observations confront Theory, ed. H. Beuther, H. Linz, & T. Henning, 387, 158 [Google Scholar]
- Beuther, H., Churchwell, E. B., McKee, C. F., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 165 [Google Scholar]
- Bik, A., Lenorzer, A., Thi, W. F., et al. 2008, in Massive Star Formation: Observations confront Theory, ed. H. Beuther, H. Linz, & T. Henning, ASP Conf. Ser., 387, 78 [Google Scholar]
- Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615 [NASA ADS] [CrossRef] [Google Scholar]
- Bodenheimer, P., Yorke, H. W., Rozyczka, M., et al. 1990, ApJ, 355, 651 [NASA ADS] [CrossRef] [Google Scholar]
- Boley, A. C., Durisen, R. H., Nordlund, Å., et al. 2007, ApJ, 665, 1254 [NASA ADS] [CrossRef] [Google Scholar]
- Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [Google Scholar]
- Dullemond, C. P., & Turolla, R. 2000, A&A, 360, 1187 [NASA ADS] [Google Scholar]
- Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Edgar, R., & Clarke, C. 2003, MNRAS, 338, 962 [NASA ADS] [CrossRef] [Google Scholar]
- Efstathiou, A., & Rowan-Robinson, M. 1990, MNRAS, 245, 275 [NASA ADS] [Google Scholar]
- Ensman, L. 1994, ApJ, 424, 275 [NASA ADS] [CrossRef] [Google Scholar]
- Foster, P. N., & Chevalier, R. A. 1993, ApJ, 416, 303 [NASA ADS] [CrossRef] [Google Scholar]
- Hayes, J. C., Norman, M. L., Fiedler, R. A., et al. 2006, ApJS, 165, 188 [NASA ADS] [CrossRef] [Google Scholar]
- Henning, T. 2008, in Massive Star Formation: Observations confront Theory, ed. H. Beuther, H. Linz, & T. Henning, ASP Conf. Ser., 387, 452 [Google Scholar]
- Hunter, C. 1977, ApJ, 218, 834 [NASA ADS] [CrossRef] [Google Scholar]
- Kahn, F. D. 1974, A&A, 37, 149 [NASA ADS] [Google Scholar]
- Klahr, H. H., & Kley, W. 2006, A&A, 445, 747 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Klahr, H. H., Henning, T., & Kley, W. 1999, ApJ, 514, 325 [NASA ADS] [CrossRef] [Google Scholar]
- Kley, W. 1989, A&A, 208, 98 [NASA ADS] [Google Scholar]
- Krumholz, M. R., Klein, R. I., & McKee, C. F. 2005a, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 8271 [Google Scholar]
- Krumholz, M. R., McKee, C. F., & Klein, R. I. 2005b, ApJ, 618, L33 [NASA ADS] [CrossRef] [Google Scholar]
- Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007a, ApJ, 656, 959 [NASA ADS] [CrossRef] [Google Scholar]
- Krumholz, M. R., Klein, R. I., McKee, C. F., et al. 2007b, ApJ, 667, 626 [NASA ADS] [CrossRef] [Google Scholar]
- Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009, Science, 323, 754 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
- Kuiper, R., Flock, M., & Klahr, H. H. 2009, in Cosmic Magnetic Fields: From Planets, to Stars and Galaxies, ed. K. G. Srassmeier, A. G. Kosovichev, & J. E. Beckman, Proc. IAU Symp., 259, 103 [Google Scholar]
- Kumar, M. S. N., & Grave, J. M. C. 2008, in Massive Star Formation: Observations confront Theory, ed. H. Beuther, H. Linz, & T. Henning, 387, 323 [Google Scholar]
- Larson, R. B., & Starrfield, S. 1971, A&A, 13, 190 [NASA ADS] [Google Scholar]
- Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
- McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 [NASA ADS] [CrossRef] [Google Scholar]
- Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
- Mihalas, D., & Mihalas, B. W. 1984 (Oxford University Press: Dover Publications edition 1999) [Google Scholar]
- Murray, S. D., Castor, J. I., Klein, R. I., et al. 1994, ApJ, 435, 631 [NASA ADS] [CrossRef] [Google Scholar]
- Nakano, T. 1989, ApJ, 345, 464 [NASA ADS] [CrossRef] [Google Scholar]
- Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793 [Google Scholar]
- Penston, M. V. 1969a, MNRAS, 145, 457 [NASA ADS] [Google Scholar]
- Penston, M. V. 1969b, MNRAS, 144, 425 [NASA ADS] [CrossRef] [Google Scholar]
- Rowan-Robinson, M. 1980, ApJS, 44, 403 [NASA ADS] [CrossRef] [Google Scholar]
- Saad, Y., & Schultz, M. H. 1986, SIAM J. Scientific Statistical Comp., 7, 856 [Google Scholar]
- Shu, F. H. 1977, ApJ, 214, 488 [NASA ADS] [CrossRef] [Google Scholar]
- Shu, F. H., Lizano, S., & Adams, F. C. 1987, in Star forming regions, ed. M. Peimbert, & J. Jugaku, 115, 417 [Google Scholar]
- Sincell, M. W., Gehmeyr, M., & Mihalas, D. 1999a, Shock Waves, 9, 391 [NASA ADS] [CrossRef] [Google Scholar]
- Sincell, M. W., Gehmeyr, M., & Mihalas, D. 1999b, Shock Waves, 9, 403 [NASA ADS] [CrossRef] [Google Scholar]
- Steinacker, J., Chini, R., Nielbock, M., et al. 2008, in Massive Star Formation: Observations confront Theory, ed. H. Beuther, H. Linz, & T. Henning, ASP Conf. Ser., 387, 270 [Google Scholar]
- Tscharnuter, W. M., & Boss, A. P. 1993, in Protostars and Planets III, ed. E. H. Levy, & J. I. Lunine, 921 [Google Scholar]
- Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 [NASA ADS] [CrossRef] [Google Scholar]
- Turner, N. J., Quataert, E., & Yorke, H. W. 2007, ApJ, 662, 1052 [NASA ADS] [CrossRef] [Google Scholar]
- van Leer, B. 1979, JCP, 32, 101 [NASA ADS] [Google Scholar]
- Vehoff, S., Nürnberger, D. E. A., Hummel, C. A., et al. 2008, in Massive Star Formation: Observations confront Theory, ed. H. Beuther, H. Linz, & T. Henning, ASP Conf. Ser., 387, 444 [Google Scholar]
- Whitworth, A. P., & Summers, D. 1985, MNRAS, 214, 1 [NASA ADS] [CrossRef] [Google Scholar]
- Wolfire, M. G., & Cassinelli, J. P. 1986, ApJ, 310, 207 [CrossRef] [Google Scholar]
- Wolfire, M. G., & Cassinelli, J. P. 1987, ApJ, 319, 850 [NASA ADS] [CrossRef] [Google Scholar]
- Yorke, H. W., & Krügel, E. 1977, A&A, 54, 183 [NASA ADS] [Google Scholar]
- Yorke, H. W., & Sonnhalter, C. 2002, ApJ, 569, 846 [NASA ADS] [CrossRef] [Google Scholar]
- Zel'Dovich, Y. B., & Raizer, Y. P. 1967 (New York: Academic Press. Dover Publications, edition 2002) [Google Scholar]
- Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Table 1: Overview of performed Monte-Carlo runs for comparison.
Table 2: Overview of runs using the proposed approximate radiation transport module.
All Figures
![]() |
Figure 1: Schematic flow chart of the radiation module. Exponents declare the iteration (timestep) number n. The actual timestep corresponds to d tn = tn - tn-1. |
Open with DEXTER | |
In the text |
![]() |
Figure 2:
Regarding the frequency dependence of stellar irradiation feedback.
Frequency dependent opacities |
Open with DEXTER | |
In the text |
![]() |
Figure 3:
Radial cut through the temperature profile in the midplane in the most
optically thin case |
Open with DEXTER | |
In the text |
![]() |
Figure 4:
Radial cut through the temperature profile in the midplane in the most
optically thick case |
Open with DEXTER | |
In the text |
![]() |
Figure 5:
Radial cut through the temperature profile in the midplane in the most
optically thick case |
Open with DEXTER | |
In the text |
![]() |
Figure 6:
Polar cut through the temperature profile at r=2 AU
of the frequency dependent irradiation plus Flux Limited Diffusion run
for the most optically thick case |
Open with DEXTER | |
In the text |
![]() |
Figure 7:
Radial cut through the radiative force profile at a polar angle of |
Open with DEXTER | |
In the text |
![]() |
Figure 8:
Radial profile of the optical depth |
Open with DEXTER | |
In the text |
![]() |
Figure 9:
Speedup factors S=t2
/ tn (
upper panel) and corresponding efficiencies |
Open with DEXTER | |
In the text |
![]() |
Figure 10: Radiative supercritical shock. Resulting density, pressure, velocity, and temperature distributions at four different snapshots in time. Dashed lines represent the adiabatic runs, solid lines the radiative ones. The time snapshots are taken ( from left to right) at 860 s, 4000 s, 7500 s, and 13 000 s after launching. Mostly horizontal lines at the lower border of the graphics refer to the initial setup. The snapshots at 4000 s are additionally marked by open circles for every 10th grid cell to illustrate the resolution used. The spatial axes are retranslated into the non-moving frame used in the visualization by Ensman (1994) for the sake of comparison. The filled dots in the four profiles at 13 000 s represent characteristic sample points in the two-temperature approach read off Figs. 10-12 in the original study by Ensman (1994). |
Open with DEXTER | |
In the text |
![]() |
Figure 11: Radiative subcritical shock. Resulting density, pressure, velocity, and temperature distributions at five different snapshots in time. Dashed lines represent the adiabatic runs, solid lines the radiative ones. The time snapshots are taken ( from left to right) at 350 s, 5400 s, 17 000 s, 28 000 s, and 38 000 s after launching. The snapshot at 350 s mostly refers to the initial setup. The snapshots at 5400 s are additionally marked by open circles for every 10th grid cell to illustrate the resolution used. The spatial axes are retranslated into the non-moving frame used in the visualization by Ensman (1994) for the sake of comparison. The filled dots in the temperature profile at 38 000 s represent characteristic sample points in the two-temperature approach read off Fig. 8 in the original study by Ensman (1994). |
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.