Fast and accurate frequency-dependent radiation transport for hydrodynamics simulations in massive star formation

R. Kuiper; H. Klahr; C. Dullemond; W. Kley; T. Henning

doi:10.1051/0004-6361/200912355

Home

All issues

Volume 511 (February 2010)

A&A, 511 (2010) A81

Full HTML

Free Access

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

A&A 511, A81 (2010)

Fast and accurate frequency-dependent radiation transport for hydrodynamics simulations in massive star formation

R. Kuiper^1,3 - H. Klahr¹ - C. Dullemond¹ - W. Kley² - T. Henning¹

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 $L_*/M_* = (4 \pi G c) / \kappa_*$ , where L_*, M_*, and $\kappa_*$ 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 $E_{\rm R}$ . Within a given spatial density $\rho(\vec{x})$ and temperature $T(\vec{x})$ distribution, we start from the time evolution of the internal energy density $E_{\rm int}$ and thermal radiation energy density $E_{\rm R}$ :

$\displaystyle \partial_t E_{\rm int} + \vec{\nabla} \cdot \left(E_{\rm int} ~ \vec{u} \right)$	=	$\displaystyle - P ~ \vec{\nabla} \cdot \vec{u} + \Lambda$	(1)
$\displaystyle \partial_t E_{\rm R} + \vec{\nabla} \cdot \left(E_{\rm R} ~ \vec{u} \right)$	=	$\displaystyle - \vec{\nabla} \cdot \vec{F} - \Lambda$	(2)

with the corresponding thermal pressure P, dynamical velocity $\vec{u}$ , and flux of radiation energy density $\vec{F}$ . The radiative heating and cooling of the gas is covered in $\Lambda = \rho c \kappa_{\rm R} \left(a T^4 - E_{\rm R} \right)$ , where c is the speed of light, $\kappa_{\rm R}$ the Rosseland mean opacity, and a the radiation constant. Gas and dust temperatures are assumed to be the same $\left(T_{\rm gas} = T_{\rm dust} = T\right)$ .

We add up the Eqs. (1) and (2) and solve the transport term $\vec{\nabla} \left(E_{\rm R} \vec{u} \right)$ separately, if necessary, in the dynamical problem. If we consider the transport of the internal energy $\vec{\nabla} \cdot \left(E_{\rm int} ~ \vec{u} \right)$ already during the corresponding hydrodynamics step (operator splitting), the remaining terms yield

$\begin{displaymath} \partial_t \left(E_{\rm int} + E_{\rm R}\right) = - \vec{\nabla} \cdot \vec{F} + Q^+, \end{displaymath}$

(3)

where the source term $Q^+ = - P \vec{\nabla} \cdot \vec{u} + \ldots$ 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 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 $t_{\rm flow} \approx l / u$ , where l and u are the characteristic length and velocity scales of the problem. The gas-energy equilibrium time is of the order of $t_{\rm eq} \approx (\kappa \rho c)^{-1}$ , since the energy exchange rate is $\dot{E}_{\rm R} \approx \kappa \rho c E_{\rm R}$ . Thus, the requirement that $t_{\rm flow} \gg t_{\rm eq}$ reduces to the requirement that $\tau \gg u / c$ , where $\tau \approx \kappa \rho l$ 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

$\begin{displaymath} E_{\rm R} = a T^4. \end{displaymath}$

(4)

This expression plus the internal energy density $E_{\rm int} = c_{\rm V} \rho T$ with the specific heat capacity $c_{\rm V}$ yield

$\begin{displaymath} \partial_t E_{\rm int} = c_{\rm V} \rho \partial_t T = \frac{c_{\rm V} \rho}{4 a T^3} \partial_t E_{\rm R}\cdot \end{displaymath}$

(5)

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

$\begin{displaymath} \partial_t E_{\rm R} = f_{\rm c} \left( - \vec{\nabla} \cdot \vec{F} + Q^+ \right) \end{displaymath}$

(6)

with $f_{\rm c} = \left(c_{\rm V} ~ \rho / \left(4 a T^3 \right) + 1\right)^{-1}$ , depending only on the ratio of internal to radiation energy. The flux $\vec{F}$ of radiation energy density is determined in the FLD approximation via

$\begin{displaymath} \vec{F} = - D \vec{\nabla} E_{\rm R} \end{displaymath}$

(7)

with the diffusion constant $D = \lambda c / \left(\kappa_{\rm R} ~ \rho \right)$ . The flux limiter $\lambda$ is chosen according to Levermore & Pomraning (1981). Scattering is neglected. In the most extreme limits the flux becomes either $\vec{F} = - c E_{\rm R} \vec{\nabla} E_{\rm R} / \vert\vec{\nabla} E_{\rm R}\vert$ for highly optically thin regions (free-streaming limit) or $\vec{F} = - c \vec{\nabla} E_{\rm R} / \left(3 ~ \kappa_{\rm R} ~ \rho \right)$ for highly optically thick regions (diffusion limit).

Summing up, we can describe the thermal/diffuse radiative processes via:

$\begin{displaymath} \partial_t E_{\rm R} = - f_{\rm c} \left( \vec{\nabla} \cdot \vec{F} - Q^+ \right) \end{displaymath}$

(8)

with

$\displaystyle f_{\rm c}$	=	$\displaystyle \left(c_{\rm V} ~ \rho / \left(4 ~ a ~ T^3 \right) + 1 \right)^{-1}$	(9)
$\displaystyle \vec{F}$	=	$\displaystyle - D ~ \vec{\nabla} E_{\rm R}$	(10)
D	=	$\displaystyle \lambda ~ c / \left(\kappa_{\rm R} ~ \rho \right)$	(11)
$\displaystyle \lambda$	=	$\displaystyle \left(2+R \right) / \left(6+3R+R^2 \right)$	(12)
R	=	$\displaystyle \vert\vec{\nabla} E_{\rm R}\vert / \left(\kappa_{\rm R} ~ \rho ~ E_{\rm R} \right)$	(13)
Q⁺	=	$\displaystyle - P ~ \vec{\nabla} \cdot \vec{u} + \mbox{ additional source terms.}$	(14)

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 $\vec{F}_*$ , released along rays in the radially outward direction. $E_{\rm R}$ is the energy density of the radiation emitted by the circumstellar material, and does not include the stellar irradiation. The additional radiation power $- \vec{\nabla} \cdot \vec{F}_*$ from this irradiation is therefore added to the diffusion Eq. (8) as a source term

$\begin{displaymath} \partial_t E_{\rm R} = - f_{\rm c} \left( \vec{\nabla} \cdot \vec{F} + \vec{\nabla} \cdot \vec{F}_* - Q^+ \right) \end{displaymath}$

(15)

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:

$\begin{displaymath} a ~ T^4 = E_{\rm R} + \frac{\kappa_{\rm P}(T_*)}{\kappa_{\rm P}(T)} \frac{\vert\vec{F}_*\vert}{c}\cdot \end{displaymath}$

(16)

Here $\kappa_{\rm P}(T)$ 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 $E_{\rm R}\left(\vec{x}\right)$ as the solution vector), the corresponding dust temperature is calculated using Eq. (16). Since the opacity $\kappa_{\rm P}\left(T\right)$ 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 $\kappa_{\rm P}\left(T_*\right) / \kappa_{\rm P}\left(T\right)$ 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

$\begin{displaymath} \vec{F}_*\left(r\right) = \vec{F}_*\left(R_*\right) ~ \left(\frac{R_*}{r}\right)^2 \exp\left(-\tau\left(r\right)\right). \end{displaymath}$

(17)

The (boundary) flux at the stellar surface $\vec{F}_*\left(R_*\right)$ can be calculated from the Stefan-Boltzmann law

$\begin{displaymath} \vert\vec{F}_*\left(R_*\right)\vert = \sigma_{\rm SB} T_*^4 \end{displaymath}$

(18)

with the Stefan-Boltzmann constant $\sigma_{\rm SB} = 0.25 a c$ . The optical depth $\tau$ 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

$\begin{displaymath} \tau\left(r\right) = \int_{R_*}^r \kappa_{\rm P}\left(T_*\right) ~ \rho\left(r\right) {\rm d}r. \end{displaymath}$

(19)

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.

$\begin{figure} \includegraphics[width=9cm]{12355fig1.eps} \end{figure}$	Figure 1: Schematic flow chart of the radiation module. Exponents declare the iteration (timestep) number n. The actual timestep corresponds to d tⁿ = tⁿ - t^n-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 $\nu_i$ , 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.

$\begin{figure} \par\includegraphics[width=12cm,clip]{12355_02.eps} \end{figure}$	Figure 2: Regarding the frequency dependence of stellar irradiation feedback. Frequency dependent opacities $\kappa \left (\nu \right )$ , Planck mean opacities $\kappa _{\rm P}\left (T_\right )$ , and stellar black body spectrum $B_\nu \left (\nu ,T_\right )$ as function of the frequency $\nu$ .
Open with DEXTER

Each dot in Fig. 2 represents the mid-frequency of the corresponding frequency bin. Including the frequency dependence, we have to replace the Planck mean opacity $\kappa _{\rm P}\left (T_*\right )$ in Eqs. (16) and (19) with the frequency dependent opacities $\kappa \left (\nu \right )$ and the stellar flux $\vec{F}_*$ in Eqs. (15) and (16) with the corresponding sum over the frequency dependent fluxes $\vec{F}_*\left(\nu\right)$ , each calculated via Eq. (17).

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 $B_\nu \left (\nu ,T_*\right )$ in the frequency dependent interval:

$\begin{displaymath} \vert\vec{F}_*\left(R_*, \nu_i\right)\vert = \frac{c}{4} \in... ...+ \nu_{i+1}\right)/2} B_\nu \left( \nu, T_*\right) {\rm d}\nu. \end{displaymath}$

(20)

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 $A {\vec x} = {\vec b}$ 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 $K_i = span\{\vec{b}, A \vec{b}, A^2 \vec{b}, ..., A^{i-1} \vec{b}\}$ . In each solver iteration i the GMRES method increments the subspace K_i used with an additional basis vector $A^{i-1} {\vec b}$ and approximates the solution of the system of linear equations by the vector ${\vec x}_i$ which minimizes the norm of the residual $r = \vert A \vec{x}_i - {\vec b}\vert$ . 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 $O\left(i^2\right)$ . 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 K₁ but with the last approximate solution $\vec{x_n}$ . Due to the growing computational effort with $O\left(i^2\right)$ 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 $\left(T_* = 5800~\mbox{K}\right)$ and the stellar radius is fixed to 1 solar radius $\left(R_* = 1~R_\odot\right)$ . The disk ranges from $r_{\rm min} = 1$ AU up to $r_{\rm max} = 1000$ 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:

$\begin{displaymath}\rho\left(r, z\right) = \rho_0 ~ f_1\left(r\right) ~ f_2\left(r, z\right) \end{displaymath}$

(21)

with the radially and vertically dependent functions

$\begin{displaymath}f_1\left(r\right) = \frac{r_{\rm d}}{r} \end{displaymath}$

(22)

and

$\begin{displaymath}f_2\left(r, z\right) = \exp \left(- \frac{\pi}{4} \left( \frac{z}{h\left(r\right)}\right)^2 \right) \end{displaymath}$

(23)

making use of the following abbreviations

$\displaystyle h\left(r\right)$	=	$\displaystyle z_{\rm d} \left( \frac{r}{r_{\rm d}} \right)^{1.125},$	(24)
$\displaystyle r_{\rm d}$	=	$\displaystyle \frac{r_{\rm max}}{2} = 500~\mbox{ AU} \mbox{, and}$	(25)
$\displaystyle z_{\rm d}$	=	$\displaystyle \frac{r_{\rm max}}{8} = 125~\mbox{ AU.}$	(26)

The lowest density is limited to a relative factor of 10^-100 compared to the highest density (at $r_{\rm min}$ in the midplane) to avoid divisions by zero (e.g. in the calculation of the diffusion coefficients). The normalization $\rho_0$ of the density setup is chosen to define different optical depths $\tau_{550~{\rm nm}}$ through the midplane of the corresponding circumstellar disk (at a visual wavelength of 550 nm):

$\tau_{550~{\rm nm}}$	$\rho_0 [\mbox{g}~\mbox{cm}^{-3}]$	$M_{\rm tot}$ of gas $[M_\odot]$
0.1	$8.321 \times 10^{-21}$	$1.1 \times 10^{-5}$
100	$8.321 \times 10^{-18}$	$1.1 \times 10^{-2}$

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 $180\hbox{$^\circ$ }$ .

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 $\left( \vec{\nabla} E_{\rm R} = 0 \right)$ , i.e. radiative flux over these boundaries is prohibited. The radial outer boundary is defined as a constant Dirichlet boundary corresponding to $T_0 = 14.7~\mbox{ K}$ $\left( E_{\rm R} = a ~ T_0^4 \right)$ .

The timestep used for the FLD solver is 10⁴ 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 $60 \times 31$ 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 $\tau_{550~{\rm nm}} = 0.1$ , a full run for the optically thick case $\tau_{550~{\rm nm}} = 100$ , and an additional run for the optically thick case $\tau_{550~{\rm nm}} = 100$ 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 $\tau_{550~{\rm nm}} = 0.1$ , 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%.

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_03.eps} \end{figure}$	Figure 3: Radial cut through the temperature profile in the midplane in the most optically thin case $\tau_{550~{\rm nm}} = 0.1$ . Upper panel: radial temperature slope of the gray irradiation routine (solid line) and the Monte-Carlo 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 $\tau_{550~{\rm nm}} = 100$ , 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.

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_04.eps} \end{figure}$

Figure 4:

Radial cut through the temperature profile in the midplane in the most optically thick case $\tau_{550~{\rm nm}} = 100$ without diffusion. Upper panel: radial temperature slope of gray irradiation (dashed line), frequency dependent irradiation (solid line), and the Monte-Carlo routine in the ``one-photon-limit'' (dots). Lower panel: deviations of the gray (dashed line) and frequency dependent method (solid line) from the Monte-Carlo code in percent.

Open with DEXTER

Therefore, we run the Monte-Carlo code only until every initial photon is absorbed, neglecting re-emission or scattering events. In this scenario, we are able to probe the absorption routines of Sects. 2.2 and 2.3 in detail and determine the improvement by considering the frequency dependence of the stellar irradiation. The resulting temperature profiles through the midplane for the case of gray and frequency dependent irradiation as well as the corresponding deviations from the Monte-Carlo comparison run are shown in Fig. 4. Indeed, we found that in order to limit the deviation in the absorption part of the radiation module to less than 5% it is essential to account for the frequency dependence of the stellar irradiation. That is due to the fact that the infrared part of the stellar spectrum has a lower optical depth than the UV part. Neglecting the frequency dependence results roughly in a 57% cooler disk at large radii and a 10% hotter inner rim compared to the Monte-Carlo based radiative transfer code.

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 $\tau_{550~{\rm nm}} = 100$ .

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_05.eps} \vspace*{4.5mm} \end{figure}$

Figure 5:

Radial cut through the temperature profile in the midplane in the most optically thick case $\tau_{550~{\rm nm}} = 100$ 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 Monte-Carlo routine (dots). Lower panel: deviations of the gray (dashed line) and frequency dependent run (solid line) from the Monte-Carlo code in percent.

Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_06.eps} \end{figure}$

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 $\tau_{550~{\rm nm}} = 100$ . The profile reproduces the turnover point at a polar angle of $\theta \approx 19\hbox {$^\circ $ }$ 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 Monte-Carlo 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 $r \approx 200~\mbox{ AU}$ . 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 $\theta \approx 19\hbox {$^\circ $ }$ ) 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

$\begin{displaymath} \vec{f}_{\rm rad} = \rho \kappa \frac{\vec{F}_{\rm tot}}{c}\cdot \end{displaymath}$

(27)

In our split radiation transport the total radiative flux $\vec{F}_{\rm tot}$ is given by the sum of the irradiated stellar flux $\vec{F}_*$ from the ray-tracing routine and the diffuse component $\vec{F}$ from the FLD solver. The flux calculated in the FLD approximation is $\vec{F} = -D \vec{\nabla} E_{\rm R}$ (see Sect. 2.1).

In discretized space the opacity and density are constant over a single grid cell and the irradiated stellar radiative flux $\vec{F}_*$ 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 $\tau \gg 1$ , 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 $\Delta x$ ) leads to:

$\displaystyle \vec{f}_*$	=	$\displaystyle \frac{1}{\Delta x} \int_0^{\Delta x}\rho \kappa \frac{\vec{F}_*}{c} {\rm d}x$	(28)
	=	$\displaystyle \frac{1}{c \Delta x} ~ \rho ~ \kappa \int_0^{\Delta x} \vec{F}_* {\rm d}x.$	(29)

The flux at position x inside the grid cell is given by the flux $\vec{F}_*^i$ entering the cell at the interface i and the absorption of this flux along the length x:

$\displaystyle \vec{F}_*\left(x\right) = \vec{F}_*^i {\rm e}^{-\tau} \nonumber$

with the optical depth $\tau = \kappa \rho x$ . The remaining integral yields:

$\displaystyle \vec{f}_*$	=	$\displaystyle \frac{1}{c \Delta x} \rho \kappa \vec{F}_*^i \int_0^{\Delta x} \exp\left(-\kappa \rho x \right) {\rm d}x$	(30)
	=	$\displaystyle - \frac{1}{c \Delta x} \vec{F}_*^i \left(\exp\left(-\kappa \rho \Delta x \right) - 1 \right)\cdot$	(31)

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

$\begin{displaymath} \vec{f}_* = - \frac{1}{c} \frac{\vec{F}_*^{i+1} - \vec{F}_*^i}{\Delta x}, \end{displaymath}$

(32)

which equals in continuous space the derivative of the radiative flux

$\begin{displaymath} \vec{f}_* = - \frac{1}{c} \partial_x \vec{F}_* \cdot \end{displaymath}$

(33)

Combining Eqs. (27) and (33) reflects the fact that without emission the radiative flux is given by the differential equation

$\begin{displaymath} \partial_x \vec{F}_* = - \rho \kappa \vec{F}_*. \end{displaymath}$

(34)

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

$\begin{displaymath} \vec{f}_* = - \frac{1}{c r^2} \partial_r \left(r^2 \vec{F}_*\right) \cdot \end{displaymath}$

(35)

So the purely geometrical dilution of the flux in the radially outward direction $\left(\vec{F}_*(r) \propto r^{-2}\right)$ does of course not contribute to the radiative force.

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_07.eps} \end{figure}$

Figure 7:

Radial cut through the radiative force profile at a polar angle of $\theta \approx 27\hbox {$^\circ $ }$ from the midplane in the most optically thick case $\tau _{550{\rm nm}} = 100$ . 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 Monte-Carlo 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 Monte-Carlo 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 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 (< $5\%)$ 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 $\tau_{550~{\rm nm}} = 100$ 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 ( $\vert\theta\vert > 30\hbox{$^\circ$ }$ 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 $\tau _{550{\rm nm}} = 100$ , 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 $\tau_{\rm R}\left(r\right) = \kappa_{\rm R}\left(T\right) ~ \rho\left(r\right) \Delta r$ from the outer edge of the disk through the midplane towards the center yields a final optical depth of $\tau_{\rm R} \approx 0.5$ , see Fig. 8.

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_08.eps} \end{figure}$	Figure 8: Radial profile of the optical depth $\tau _{\rm R}$ through the midplane in the most optically thick case $\tau_{550~{\rm nm}} = 100$ 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 $r \approx$ 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 $r \approx 200~\mbox{ AU}$ (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 $\nu_{\rm peak} \propto T$ . 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 $\tau_{550~\mbox{nm}} = 100$ 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 t₂ to t_n (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=t₂ / t_n. The efficiency E is determined via E = t₂ / (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 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 $64 \times 64 \times 256$ grid cells. Each processor covers therefore a $\left(64 \times 64 \times 256\right) / n$ 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.

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_09.eps} \end{figure}$	Figure 9: Speedup factors S=t₂ / t_n ( upper panel) and corresponding efficiencies $E=t_2 / \left (t_n n\right ) = S/n$ ( 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 $64 \times 64 \times 8n$ 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).

$\begin{figure} \par\includegraphics[width=17cm]{12355_10.eps} \end{figure}$

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

$\displaystyle \partial_t \rho + \vec{\nabla} \cdot \left(\rho ~ \vec{u}\right)$	=	0	(36)
$\displaystyle \partial_t \left(\rho ~ \vec{u}\right) + \vec{\nabla} \left(\rho ~ \vec{u} ~ \vec{u}\right)$	=	$\displaystyle - \vec{\nabla} ~ P + \rho ~ \vec{a}$	(37)
$\displaystyle \partial_t E + \vec{\nabla} \cdot \left(\left(E + P\right) ~ \vec{u}\right)$	=	$\displaystyle \rho ~ \vec{u} \cdot \vec{a}$	(38)

with the acceleration source term from the radiation transport

$\begin{displaymath}\vec{a} = - \frac{\vec{\nabla} \cdot \vec{F_*}}{\rho ~ c} \vec{e}_r - \kappa_{\rm R} \frac{D ~ \vec{\nabla} E_{\rm R}}{c}, \end{displaymath}$

(39)

derived in Sect. 3.3.4. The evolution of the gas density $\rho$ , the gas velocity $\vec{u}$ , the thermal pressure P, and the total energy E is computed using the open source magneto-hydrodynamics code Pluto3 (Mignone et al. 2007). Pluto is a high-order 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 Harten-Lax-VanLeer (hll) Riemann solver and a so-called ``minmod'' flux limiter using piecewise linear interpolation (plm) and a Runge-Kutta 2 (RK2) time integration (predictor-corrector-method, 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 one-dimensional domain of the grid associated covers a distance of length $l_0=7 \times 10^{10}~\mbox{ cm}$ . The iso-density in the domain is fixed to $\rho_0 = 7.78 \times 10^{-10}~\mbox{g~cm}^{-3}$ . For testing purposes, the gas is set to be completely ionized, thus the mean molecular weight is $\mu = 0.5$ . An ideal equation of state is used with $\gamma = 5/3$ . The opacity is fixed to a constant value of $\kappa = 0.4~\mbox{cm}^2~\mbox{g}^{-1}$ . The initial temperature drops linearly from 85 K at the starting position of the piston to 10 K at the outer boundary. The velocity $u_0 > c_{\rm s}$ 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 $u_{\rm c}$ .

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₀, 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

$\begin{figure} \par\includegraphics[width=17cm]{12355_11.eps} \end{figure}$

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

Higher velocity of the piston results in a stronger pre-heating of the gas directly in front of the shocked gas. When a piston velocity higher than the critical velocity $u_{\rm c}$ is used, the peak temperature of the pre-heated region (which is equal to the temperature of the shocked gas) will not increase any more, but the enlargement of the pre-heated region will extend further out. The gas temperature of the pre-heated region will presumably be in equilibrium with the radiation temperature for the most part of the domain. At one point of the sloping tail of the temperature distribution in the upstream direction both temperatures will differ. In an equilibrium diffusion (also called one-temperature diffusion) code like ours the temperature will decline at this point sharply (see also Fig. 15 in Ensman 1994, which compares an adiabatic, an equilibrium, a non-equilibrium and a full radiation transport method for this setup). The reason for this sharp decline is (according to Ensman 1994) that the radiation can only penetrate through the cold environment, if it heats this environment up to the equilibrium temperature. Once the radiation energy is absorbed the radiative flux is zero. In a non-equilibrium radiation transport method the radiation energy can penetrate further into the environment simply by heating it to less than the equilibrium temperature.

We choose the piston velocity to be the one used originally by Ensman (1994), $u_0 = 20~\mbox{km}~\mbox{s}^{-1}$ . In some of the subsequent tests mentioned above a slightly lower piston velocity $\left(u_0 = 16~\mbox{km}~\mbox{s}^{-1}\right)$ 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 ( $c_{\rm s} < u_0 < u_{\rm c}$ ) 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), $u_0 = 6~\mbox{km}~\mbox{s}^{-1}$ . 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 $c_{\rm s} \propto \sqrt{T}$ . 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).

Acknowledgements

This 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 [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., & 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 [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]
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 [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 Monte-Carlo runs for comparison.

Table 2: Overview of runs using the proposed approximate radiation transport module.

All Figures

$\begin{figure} \includegraphics[width=9cm]{12355fig1.eps} \end{figure}$	Figure 1: Schematic flow chart of the radiation module. Exponents declare the iteration (timestep) number n. The actual timestep corresponds to d tⁿ = tⁿ - t^n-1.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=12cm,clip]{12355_02.eps} \end{figure}$	Figure 2: Regarding the frequency dependence of stellar irradiation feedback. Frequency dependent opacities $\kappa \left (\nu \right )$ , Planck mean opacities $\kappa _{\rm P}\left (T_\right )$ , and stellar black body spectrum $B_\nu \left (\nu ,T_\right )$ as function of the frequency $\nu$ .
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_03.eps} \end{figure}$	Figure 3: Radial cut through the temperature profile in the midplane in the most optically thin case $\tau_{550~{\rm nm}} = 0.1$ . Upper panel: radial temperature slope of the gray irradiation routine (solid line) and the Monte-Carlo based comparison code (dots). Lower panel: differences between the two codes in percent.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_04.eps} \end{figure}$	Figure 4: Radial cut through the temperature profile in the midplane in the most optically thick case $\tau_{550~{\rm nm}} = 100$ without diffusion. Upper panel: radial temperature slope of gray irradiation (dashed line), frequency dependent irradiation (solid line), and the Monte-Carlo routine in the ``one-photon-limit'' (dots). Lower panel: deviations of the gray (dashed line) and frequency dependent method (solid line) from the Monte-Carlo code in percent.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_05.eps} \vspace*{4.5mm} \end{figure}$	Figure 5: Radial cut through the temperature profile in the midplane in the most optically thick case $\tau_{550~{\rm nm}} = 100$ 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 Monte-Carlo routine (dots). Lower panel: deviations of the gray (dashed line) and frequency dependent run (solid line) from the Monte-Carlo code in percent.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_06.eps} \end{figure}$	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 $\tau_{550~{\rm nm}} = 100$ . The profile reproduces the turnover point at a polar angle of $\theta \approx 19\hbox {$^\circ $ }$ 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 Monte-Carlo 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

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_07.eps} \end{figure}$	Figure 7: Radial cut through the radiative force profile at a polar angle of $\theta \approx 27\hbox {$^\circ $ }$ from the midplane in the most optically thick case $\tau _{550{\rm nm}} = 100$ . 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 Monte-Carlo 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 Monte-Carlo run in percent.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_08.eps} \end{figure}$	Figure 8: Radial profile of the optical depth $\tau _{\rm R}$ through the midplane in the most optically thick case $\tau_{550~{\rm nm}} = 100$ integrated from the outer edge of the disk inwards.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.8cm]{12355_09.eps} \end{figure}$	Figure 9: Speedup factors S=t₂ / t_n ( upper panel) and corresponding efficiencies $E=t_2 / \left (t_n n\right ) = S/n$ ( lower panel) for a fixed (squares) and a growing problem size (stars).
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=17cm]{12355_10.eps} \end{figure}$	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

$\begin{figure} \par\includegraphics[width=17cm]{12355_11.eps} \end{figure}$	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

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.

[1] Balay, S., Gropp, W. D., McInnes, L. C., et al. 1997, Birkhäuser Press, 163 [Google Scholar]

[2] Balay, S., Buschelman, K., Gropp, W. D., et al. 2001, PETSc Web page, http://www.mcs.anl.gov/petsc [Google Scholar]

[3] Balay, S., Buschelman, K., Eijkhout, V., et al. 2004, PETSc Users Manual, Tech. rep., PETSc, revision 2.1.5 [Google Scholar]

[4] Bally, J. 2008, in Massive Star Formation: Observations confront Theory, ed. H. Beuther, H. Linz, & T. Henning, 387, 158 [Google Scholar]

[5] 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]

[6] 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]

[7] Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615 [NASA ADS] [CrossRef] [Google Scholar]

[8] Bodenheimer, P., Yorke, H. W., Rozyczka, M., et al. 1990, ApJ, 355, 651 [NASA ADS] [CrossRef] [Google Scholar]

[9] Boley, A. C., Durisen, R. H., Nordlund, Å., et al. 2007, ApJ, 665, 1254 [NASA ADS] [CrossRef] [Google Scholar]

[10] Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [NASA ADS] [CrossRef] [Google Scholar]

[11] Dullemond, C. P., & Turolla, R. 2000, A&A, 360, 1187 [NASA ADS] [Google Scholar]

[12] Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[13] Edgar, R., & Clarke, C. 2003, MNRAS, 338, 962 [NASA ADS] [CrossRef] [Google Scholar]

[14] Efstathiou, A., & Rowan-Robinson, M. 1990, MNRAS, 245, 275 [NASA ADS] [Google Scholar]

[15] Ensman, L. 1994, ApJ, 424, 275 [NASA ADS] [CrossRef] [Google Scholar]

[16] Foster, P. N., & Chevalier, R. A. 1993, ApJ, 416, 303 [NASA ADS] [CrossRef] [Google Scholar]

[17] Hayes, J. C., Norman, M. L., Fiedler, R. A., et al. 2006, ApJS, 165, 188 [NASA ADS] [CrossRef] [Google Scholar]

[18] Henning, T. 2008, in Massive Star Formation: Observations confront Theory, ed. H. Beuther, H. Linz, & T. Henning, ASP Conf. Ser., 387, 452 [Google Scholar]

[19] Hunter, C. 1977, ApJ, 218, 834 [NASA ADS] [CrossRef] [Google Scholar]

[20] Kahn, F. D. 1974, A&A, 37, 149 [NASA ADS] [Google Scholar]

[21] Klahr, H. H., & Kley, W. 2006, A&A, 445, 747 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[22] Klahr, H. H., Henning, T., & Kley, W. 1999, ApJ, 514, 325 [NASA ADS] [CrossRef] [Google Scholar]

[23] Kley, W. 1989, A&A, 208, 98 [NASA ADS] [Google Scholar]

[24] 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]

[25] Krumholz, M. R., McKee, C. F., & Klein, R. I. 2005b, ApJ, 618, L33 [NASA ADS] [CrossRef] [Google Scholar]

[26] Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007a, ApJ, 656, 959 [NASA ADS] [CrossRef] [Google Scholar]

[27] Krumholz, M. R., Klein, R. I., McKee, C. F., et al. 2007b, ApJ, 667, 626 [NASA ADS] [CrossRef] [Google Scholar]

[28] 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]

[29] 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]

[30] 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]

[31] Larson, R. B., & Starrfield, S. 1971, A&A, 13, 190 [NASA ADS] [Google Scholar]

[32] Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]

[33] McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 [NASA ADS] [CrossRef] [Google Scholar]

[34] Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] [Google Scholar]

[35] Mihalas, D., & Mihalas, B. W. 1984 (Oxford University Press: Dover Publications edition 1999) [Google Scholar]

[36] Murray, S. D., Castor, J. I., Klein, R. I., et al. 1994, ApJ, 435, 631 [NASA ADS] [CrossRef] [Google Scholar]

[37] Nakano, T. 1989, ApJ, 345, 464 [NASA ADS] [CrossRef] [Google Scholar]

[38] Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]

[39] Penston, M. V. 1969a, MNRAS, 145, 457 [NASA ADS] [Google Scholar]

[40] Penston, M. V. 1969b, MNRAS, 144, 425 [NASA ADS] [CrossRef] [Google Scholar]

[41] Rowan-Robinson, M. 1980, ApJS, 44, 403 [NASA ADS] [CrossRef] [Google Scholar]

[42] Saad, Y., & Schultz, M. H. 1986, SIAM J. Scientific Statistical Comp., 7, 856 [CrossRef] [MathSciNet] [Google Scholar]

[43] Shu, F. H. 1977, ApJ, 214, 488 [NASA ADS] [CrossRef] [Google Scholar]

[44] Shu, F. H., Lizano, S., & Adams, F. C. 1987, in Star forming regions, ed. M. Peimbert, & J. Jugaku, 115, 417 [Google Scholar]

[45] Sincell, M. W., Gehmeyr, M., & Mihalas, D. 1999a, Shock Waves, 9, 391 [NASA ADS] [CrossRef] [Google Scholar]

[46] Sincell, M. W., Gehmeyr, M., & Mihalas, D. 1999b, Shock Waves, 9, 403 [NASA ADS] [CrossRef] [Google Scholar]

[47] 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]

[48] Tscharnuter, W. M., & Boss, A. P. 1993, in Protostars and Planets III, ed. E. H. Levy, & J. I. Lunine, 921 [Google Scholar]

[49] Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 [NASA ADS] [CrossRef] [Google Scholar]

[50] Turner, N. J., Quataert, E., & Yorke, H. W. 2007, ApJ, 662, 1052 [NASA ADS] [CrossRef] [Google Scholar]

[51] van Leer, B. 1979, JCP, 32, 101 [NASA ADS] [CrossRef] [Google Scholar]

[52] 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]

[53] Whitworth, A. P., & Summers, D. 1985, MNRAS, 214, 1 [NASA ADS] [CrossRef] [Google Scholar]

[54] Wolfire, M. G., & Cassinelli, J. P. 1986, ApJ, 310, 207 [NASA ADS] [CrossRef] [Google Scholar]

[55] Wolfire, M. G., & Cassinelli, J. P. 1987, ApJ, 319, 850 [NASA ADS] [CrossRef] [Google Scholar]

[56] Yorke, H. W., & Krügel, E. 1977, A&A, 54, 183 [NASA ADS] [Google Scholar]

[57] Yorke, H. W., & Sonnhalter, C. 2002, ApJ, 569, 846 [NASA ADS] [CrossRef] [Google Scholar]

[58] Zel'Dovich, Y. B., & Raizer, Y. P. 1967 (New York: Academic Press. Dover Publications, edition 2002) [Google Scholar]

[59] Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481 [NASA ADS] [CrossRef] [Google Scholar]