Fast and accurate frequencydependent radiation transport for hydrodynamics simulations in massive star formation
R. Kuiper^{1,3}  H. Klahr^{1}  C. Dullemond^{1}  W. Kley^{2}  T. Henning^{1}
1  MaxPlanckInstitut 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 MaxPlanck Research School for
Astronomy and Cosmic Physics at the University of Heidelberg
(IMPRSHD), 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 prestellar
cores and evolving accretion disks is therefore a main effort in
current hydrodynamics simulations.
Aims. We introduce our newly implemented
threedimensional 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
onetemperature approach, which is valid in the static diffusion limit,
with the high accuracy of a frequency dependent first order raytracing
routine. The raytracing 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
threedimensional highorder magnetohydrodynamics code Pluto.
Results. We prove the viability of the scheme in a
standard radiation benchmark test compared to a full frequency
dependent MonteCarlo 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 prestellar 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 shortliving, 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 lowmass 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 infalling 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 highmass stars is basically a scaledup version of lowmass star formation, implies its own challenges: Their rapid evolution, especially the shorter KelvinHelmholtz 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 prestellar core (without angular momentum), the infall 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 onedimensional 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 selfsimilar solutions. Also stellar radiative feedback has been considered in onedimensional 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 twodimensional nonspherically 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 socalled 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 socalled 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 RayleighTaylor 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 highresolution simulations of the accretion process.
Regarding these issues, we demonstrate here the usability of our newly developed fast threedimensional frequency dependent approximate radiation transport module for numerical hydrodynamics.
The most accurate description of the physics proceeding during the collapse of a prestellar core would include a frequency dependent radiation transport step after each hydrodynamic timestep using a modern MonteCarlo based or raytracing 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 twodimensional 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 onedimensional 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 twodimensional 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. Raytracing based methods (like in Efstathiou & RowanRobinson 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 (freefloating) 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 MonteCarlo 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 raytracing. The approximation described here results in a large reduction in computing time compared to MonteCarlo based radiative transfer. This allows us to implement this particular code in the framework of threedimensional (magneto) hydrodynamics simulations of circumstellar disks and infalling envelopes on a parallel decomposed (spherical) grid. In particular, we implemented our radiation transport module into the threedimensional highorder Godunov magnetohydrodynamics code Pluto (Mignone et al. 2007). This tool allows us to perform highresolution 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 magnetohydrodynamics 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 threedimensional 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 raytracing 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 reemission 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 raytracing 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 raytracing 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 raytracing 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 raytracing step. In the last Sect. 2.4, we comment on the socalled 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 , and flux of radiation energy density . The radiative heating and cooling of the gas is covered in , where c is the speed of light, the Rosseland mean opacity, and a the radiation constant. Gas and dust temperatures are assumed to be the same .
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 depends on the physics included and can additionally contain additional source terms such as accretion luminosity from sink cells or viscous heating.
In the following, we use the assumption that the gas and radiation temperature are in equilibrium (also called onetemperature 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 gasradiation equilibrium time only. The flow time is , where l and u are the characteristic length and velocity scales of the problem. The gasenergy 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 onetemperature approach yields the correct temperature of the shocked gas, but results in a steeper temperature gradient than a twotemperature 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 the specific heat capacity yield
With this relation, the Eq. (3) reduces to a standard diffusion equation
with , depending only on the ratio of internal to radiation energy. The flux of radiation energy density is determined in the FLD approximation via
with the diffusion constant . The flux limiter is chosen according to Levermore & Pomraning (1981). Scattering is neglected. In the most extreme limits the flux becomes either for highly optically thin regions (freestreaming limit) or for highly optically thick regions (diffusion limit).
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 represents the Planck mean opacity to a given temperature T and T_{*} is the star's effective temperature.
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 NewtonRaphson 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 can be calculated from the StefanBoltzmann law
with the StefanBoltzmann constant . The optical depth is, in the case of gray irradiation (see following Sect. 2.3 for specific changes due to frequency dependent irradiation), calculated along radial rays through the spherical grid via
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 t^{n} = t^{n}  t^{n1}. 

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 midfrequency , 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 , Planck mean opacities , and stellar black body spectrum as function of the frequency . 

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 overrelaxation (SOR) method, and at least as good as the wellknown improved stabilized BiConjugate 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 K_{i} 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 socalled ``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 K_{1} 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 MonteCarlo 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 solartype star, an irradiated circumstellar flared disk and an envelope. We have to choose a lowmass central star, because no benchmark for highmass 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 MonteCarlo 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 reemission) with both the MonteCarlo 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 solarlike: 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 in the midplane) to avoid divisions by zero (e.g. in the calculation of the diffusion coefficients). The normalization of the density setup is chosen to define different optical depths through the midplane of the corresponding circumstellar disk (at a visual wavelength of 550 nm):
of gas  
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, twodimensional 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 BlockJacobian and serial/local ILU preconditioning 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 10^{4} 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 MonteCarlo based comparison code RADMC
For comparison we use the MonteCarlo 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 MonteCarlo 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 raytracing 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.12.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 MonteCarlo 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 reemission of the photons (to achieve a pure absorption scenario for comparison with our first order raytracing routine), hereafter called ``onephotonlimit''. An overview of these three comparison runs is given in Table 1.
The resulting temperature distributions of these MonteCarlo 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 MonteCarlo 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 . Upper panel: radial temperature slope of the gray irradiation routine (solid line) and the MonteCarlo based comparison code (dots). Lower panel: differences between the two codes in percent. 

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 without diffusion. Upper panel: radial temperature slope of gray irradiation (dashed line), frequency dependent irradiation (solid line), and the MonteCarlo routine in the ``onephotonlimit'' (dots). Lower panel: deviations of the gray (dashed line) and frequency dependent method (solid line) from the MonteCarlo code in percent. 

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 MonteCarlo simulation in the most optically thick case .
Figure 5: Radial cut through the temperature profile in the midplane in the most optically thick case including irradiation and Flux Limited Diffusion. Upper panel: radial temperature slope of gray irradiation plus Flux Limited Diffusion (dashed line), frequency dependent irradiation plus Flux Limited Diffusion (solid line), and the corresponding MonteCarlo routine (dots). Lower panel: deviations of the gray (dashed line) and frequency dependent run (solid line) from the MonteCarlo code in percent. 

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 . The profile reproduces the turnover point at a polar angle of above the midplane from the optically thin envelope to the optically thick disk region. Solid line: Frequency dependent irradiation plus Flux Limited Diffusion. Dots: Data from the corresponding MonteCarlo comparison run. The vertical axis covers only a small temperature range from 200 to 260 K to better visualize the small deviations. 

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 MonteCarlo 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 is given by the sum of the irradiated stellar flux from the raytracing routine and the diffuse component from the FLD solver. The flux calculated in the FLD approximation is (see Sect. 2.1).
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 onedimensional cartesian grid with
an uniform grid spacing of )
leads to:
The flux at position x inside the grid cell is given by the flux entering the cell at the interface i and the absorption of this flux along the length x:
with the optical depth . The remaining integral yields:
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 does of course not contribute to the radiative force.
Figure 7: Radial cut through the radiative force profile at a polar angle of from the midplane in the most optically thick case . The profile visualizes the onset of the radiative force at the transition from the optically thin envelope to the optically thick disk region. Upper panel: radial radiative force profile of frequency dependent irradiation plus Flux Limited Diffusion (solid line with squares) as well as the result from the MonteCarlo routine (solid line). Lower panel: deviations of purely stellar (dotted line with triangles) and stellar plus thermal radiative force (solid line with squares) from the MonteCarlo run in percent. 

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 MonteCarlo 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 34% 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 MonteCarlo 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 MonteCarlo 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 raytracing part of the optically thick case (Fig. 4) stay beneath the discrepancy of the different MonteCarlo 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 socalled photon noise in the MonteCarlo 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 predescribed 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 through the midplane in the most optically thick case integrated from the outer edge of the disk inwards. 

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 (freestreaming limit) as well as in the most optically thick (diffusion limit) regions only. The apparent yet surprisingly good agreement between the MonteCarlo 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 midfrequency 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 BlockJacobian preconditioner, 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 raytracing 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 t_{2} to t_{n} (n is the number of processors) represent the wall clock time per main iteration per grid cell without the nonrecurring 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=t_{2} / t_{n}. The efficiency E is determined via E = t_{2} / (t_{n} 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 64Bit 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=t_{2} / t_{n} ( upper panel) and corresponding efficiencies ( lower panel) for a fixed (squares) and a growing problem size (stars). 

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 threedimensional radiative hydrodynamics studies of collapsing massive prestellar 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 threedimensional hydrodynamics code for spherical grids (well known in the astrophysical community, e.g. Pluto2, Flash2.5, and ZeusMP1.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), ZEUS2D (Turner & Stone 2001), as well as ZEUSMP (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 nonmoving 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 twotemperature approach read off Figs. 1012 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 gas velocity , the thermal pressure P, and the total energy E is computed using the open source magnetohydrodynamics code Pluto3 (Mignone et al. 2007). Pluto is a highorder Godunov solver code, i.e. it uses a shock capturing Riemann solver within a conservative finite volume scheme. Our default configuration consists further of a HartenLaxVanLeer (hll) Riemann solver and a socalled ``minmod'' flux limiter using piecewise linear interpolation (plm) and a RungeKutta 2 (RK2) time integration (predictorcorrectormethod, compare e.g. van Leer 1979). Therefore the total difference scheme is 2nd order accurate in time and space.
The test setup describes a piston moving with supersonic velocity through an initially uniform, cold gas. The onedimensional domain of the grid associated covers a distance of length . The isodensity 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 preheating, a preacceleration, as well as a precompression of the gas directly in front of the shock. If the temperature in this preheated region stays below the temperature of the shocked gas, it is called a subcritical radiative shock. If the temperature in the preheated 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 u_{0}, 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 nonmoving 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 ZEUSMP 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 nonmoving 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 twotemperature 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. 1012 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 preheating, preacceleration, and precompression 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 nonequilibrium 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 nonequilibrium 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 nonequilibrium 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 nonequilibrium 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 longrange tail of the preheated region, if using a nonequilibrium radiation transport method. Summing up, the deviations of the equilibrium temperature method from a nonequilibrium approach stay below 10% in the pre and postshock regions. In this dynamic diffusion shock test, the outermost tail of the preheated 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 nonequilibrium 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 magnetohydrodynamics 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 raytracing 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 raytracing 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 (freestreaming) 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 raytracing 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 raytracing 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 raytracing 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 raytracing routine. RowanRobinson (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 reemission of once absorbed stellar photons, i.e. the expansion of the raytracing 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 (MonteCarlo 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 MonteCarlo 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 raytracing 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 onetemperature 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 twotemperature 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 prestellar 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, threedimensional ideal and nonideal MHD simulations of the magnetorotational instability (MRI) in protoplanetary 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 lowmass stars. The algorithm is fully implemented in our version of the magnetohydrodynamics code Pluto (including highorder 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 MaxPlanck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRSHD). 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 magnetohydrodynamics 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 [NASA ADS] [CrossRef] [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., & RowanRobinson, 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 [NASA ADS] [CrossRef] [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 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [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]
 RowanRobinson, M. 1980, ApJS, 44, 403 [NASA ADS] [CrossRef] [Google Scholar]
 Saad, Y., & Schultz, M. H. 1986, SIAM J. Scientific Statistical Comp., 7, 856 [CrossRef] [MathSciNet] [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] [CrossRef] [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 [NASA ADS] [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 MonteCarlo 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 t^{n} = t^{n}  t^{n1}. 

Open with DEXTER  
In the text 
Figure 2: Regarding the frequency dependence of stellar irradiation feedback. Frequency dependent opacities , Planck mean opacities , and stellar black body spectrum as function of the frequency . 

Open with DEXTER  
In the text 
Figure 3: Radial cut through the temperature profile in the midplane in the most optically thin case . Upper panel: radial temperature slope of the gray irradiation routine (solid line) and the MonteCarlo based comparison code (dots). Lower panel: differences between the two codes in percent. 

Open with DEXTER  
In the text 
Figure 4: Radial cut through the temperature profile in the midplane in the most optically thick case without diffusion. Upper panel: radial temperature slope of gray irradiation (dashed line), frequency dependent irradiation (solid line), and the MonteCarlo routine in the ``onephotonlimit'' (dots). Lower panel: deviations of the gray (dashed line) and frequency dependent method (solid line) from the MonteCarlo code in percent. 

Open with DEXTER  
In the text 
Figure 5: Radial cut through the temperature profile in the midplane in the most optically thick case including irradiation and Flux Limited Diffusion. Upper panel: radial temperature slope of gray irradiation plus Flux Limited Diffusion (dashed line), frequency dependent irradiation plus Flux Limited Diffusion (solid line), and the corresponding MonteCarlo routine (dots). Lower panel: deviations of the gray (dashed line) and frequency dependent run (solid line) from the MonteCarlo code in percent. 

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 . The profile reproduces the turnover point at a polar angle of above the midplane from the optically thin envelope to the optically thick disk region. Solid line: Frequency dependent irradiation plus Flux Limited Diffusion. Dots: Data from the corresponding MonteCarlo comparison run. The vertical axis covers only a small temperature range from 200 to 260 K to better visualize the small deviations. 

Open with DEXTER  
In the text 
Figure 7: Radial cut through the radiative force profile at a polar angle of from the midplane in the most optically thick case . The profile visualizes the onset of the radiative force at the transition from the optically thin envelope to the optically thick disk region. Upper panel: radial radiative force profile of frequency dependent irradiation plus Flux Limited Diffusion (solid line with squares) as well as the result from the MonteCarlo routine (solid line). Lower panel: deviations of purely stellar (dotted line with triangles) and stellar plus thermal radiative force (solid line with squares) from the MonteCarlo run in percent. 

Open with DEXTER  
In the text 
Figure 8: Radial profile of the optical depth through the midplane in the most optically thick case integrated from the outer edge of the disk inwards. 

Open with DEXTER  
In the text 
Figure 9: Speedup factors S=t_{2} / t_{n} ( upper panel) and corresponding efficiencies ( lower panel) for a fixed (squares) and a growing problem size (stars). 

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 nonmoving 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 twotemperature approach read off Figs. 1012 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 nonmoving 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 twotemperature approach read off Fig. 8 in the original study by Ensman (1994). 

Open with DEXTER  
In the text 
Copyright ESO 2010