© ESO, 2013

1. Introduction

The formation, morphology, and evolution of circumstellar disks depend crucially on their temperature distribution. After their formation, angular momentum transport can be enabled by various instabilities. Purely hydrodynamic instabilities (e.g., Lin & Bodenheimer 1982; Godon & Livio 2000; Klahr & Bodenheimer 2003; Umurhan & Regev 2004; Johnson & Gammie 2005), magneto-rotational instability (Balbus & Hawley 1991; Hawley & Balbus 1991), and gravitational instability (Boss 1984) denote important paths to sustain accretion. The strength of the angular momentum transport and accretion depends on the temperature distribution of the disk. For purely hydrodynamic instabilities, this was shown, for example, in the publication series Kley & Lin (1992), Kley et al. (1993b), and Kley et al. (1993a); for magneto-rotational instability in local shearing box simulations it was shown by Hirose & Turner (2011) and Flaig et al. (2012), and for gravitational instability in, for example, Boley et al. (2007).

Also, the stability or fragmentation of circumstellar disks, which leads to the formation of planets, brown dwarfs, or binary stars, depends mainly on two temperature-dependent criteria, the Toomre criterion (Toomre 1964) and the cooling timescale (Gammie 2001). Based on these criteria, semi-analytical studies of massive accretion disks (Kratter & Matzner 2006; Kratter et al. 2008; Vaidya et al. 2009) derive the disk’s stability and its fragmentation properties. The study by Vaidya et al. (2009) includes gas and dust opacities to obtain dynamical quantities in the midplane of disks around typical massive stars and concludes that viscous heating is the dominant phenomenon in destroying dust up to roughly 10 AU from the central star. This location of the dust evaporation front is in good agreement with the dynamical disk-formation studies by Kuiper & Yorke (2013), which also take the effect of optically thick gas (disks) around massive stars into account.

Planetary companions may form in the accretion disk directly by disk fragmentation or by core-accretion (Pollack et al. 1996; Alibert et al. 2004, 2005). Based on the core-accretion scenario, statistical properties of extrasolar planets can be determined using an extrasolar planet population synthesis model (Mordasini et al. 2009a,b, 2012; Alibert et al. 2011). Effects of the disk’s temperature distribution (due to irradiation) on the outcome of such planet population synthesis models are studied in Fouchet et al. (2012).

After a companion has formed in an accretion disk, its radial drift or migration plays a major role in the further evolution of the solid body-disk interacting system. In this late stage of disk evolution, the effect of radiation on migration is an ongoing problem in modern astrophysical research (see, e.g., Bitsch et al. 2013).

Summing up, a reliable determination of the disk’s temperature distribution is a crucial ingredient of current and future (magneto-) hydrodynamical studies of all evolutionary epochs of circumstellar disks. To account for the proper temperature distribution, (magneto-) hydrodynamics studies of irradiated circumstellar disks require an accurate treatment of radiation transport.

At the same time this accurate radiation transport should desirably be as fast as the (magneto-) hydrodynamics to allow for an efficient usage of computing resources. Such fast radiation transport methods imply the acceptance of far-reaching approximations, e.g. neglecting frequency dependence or making use of the diffusion approximation. Although Harries (2011) recently combined Monte-Carlo radiation transport with hydrodynamics, they conclude that the Monte-Carlo calculations are so expensive in cpu time that the time spent for the hydrodynamics part of the simulations becomes negligible. By consuming roughly the same computing resources for the radiation transport than for the gas motion itself, approximate radiation transport methods make it possible, for example, to study the system under investigation for much longer evolution in time, run a whole set of simulations to scan the parameter space, or achieve higher resolution.

Here, we extend our previous comparison test of radiation transport methods (Kuiper et al. 2010b) by considering three types of approximate radiation transport methods and expanding the setup of the disks to very high optical depth. As physical test setup, we choose the standard radiation benchmark tests by Pascucci et al. (2004) and Pinte et al. (2009). We perform simulations for a wide range of optical depths of the disk’s midplane from τ_550 nm = 0.1 up to τ_810 nm = 1.22 × 10⁺⁶. For comparison, we use results from the Monte-Carlo code RADMC (Dullemond & Turolla 2000; Dullemond & Dominik 2004a).

We determine the accuracy both of the gray and frequency-dependent hybrid radiation transport method, introduced in Kuiper et al. (2010b), and of the gray flux-limited diffusion (FLD) approximation, which denotes a standard technique in current astrophysical research. In the hybrid method, stellar irradiation is computed in a ray-tracing (RT) step and the thermal (re-)emission of dust is shifted to a gray FLD solver. Previous studies using the hybrid radiation transport approach covered the formation of massive stars and their radiation pressure feedback (Kuiper et al. 2010a, 2011; Kuiper & Yorke 2013), along with the stability of radiation-pressure-dominated outflow cavities around forming high-mass stars (Kuiper et al. 2012). The latter study emphasizes the importance of the long-range effect of stellar irradiation to compute the stability of the outflow cavities. It shows that the usage of the gray FLD approximation yields an underestimate of the radiative force in the expanding cavity shell by 1−2 orders of magnitude (depending on the size of the cavity and the stellar parameters). These results were later confirmed by radiation-hydrodynamics simulations of Harries et al. (2012) and radiation transport models of Owen et al. (2012). Another example of the importance of such radiative long-range effects was already given in three-dimensional radiation-hydrodynamics studies of proto-planetary disks by Boley et al. (2007).

For a historical overview of the development of approximate radiation transport methods we refer the reader to the introduction of Kuiper et al. (2010b).

2. Methods

The equilibrium temperature for static disk configurations is computed using different kinds of radiation transport methods, including fast approximate radiation transport solvers and a highly accurate Monte-Carlo code.

As approximate radiation transport solvers, we use the gray FLD approximation, the gray irradiation + gray FLD approximation, and the frequency-dependent irradiation + gray FLD approximation. We have described the details of these techniques and their actual numerical implementation in our radiation transport module called MAKEMAKE in Kuiper et al. (2010b) and Kuiper et al. (2012). Due to their speed, these approximate radiation transport solvers are ideally suited to be used in (magneto-) hydrodynamics simulations as well.

As a comparison code, we use the Monte-Carlo based radiation transfer code RADMC, which is described in Dullemond & Turolla (2000) and Dullemond & Dominik (2004a). The general solver method of RADMC is based on Bjorkman & Wood (2001). The code is also described briefly in the original radiation benchmark tests of Pascucci et al. (2004) and Pinte et al. (2009). While scattering can be handled by this code, it is simply switched off because it is also neglected in the approximate radiation transport methods for (magneto-) hydrodynamics simulations. For the given disk and envelope configurations of moderate optical depth τ_550 nm = 100, scattering would increase the temperature in the irradiated parts by about 2% in the optically thin envelope and up to a maximum of 19% in the optically thick inner rim of the disk’s midplane (up to an optical depth of about unity) due to higher extinction. The more effectively shielded outer regions of the disk would be about 10% cooler.

3. Physical setup, numerical configuration, and runs performed

3.1. Physical setup

The setup is adopted from the original standard radiation benchmark tests Pascucci et al. (2004) and Pinte et al. (2009). It describes a slightly flared circumstellar disk within an optically thin envelope around a central proto-star. Such a flaring is assumed to be a natural consequence of the heating (e.g., Kenyon & Hartmann 1987; Chiang & Goldreich 1997), but generally geometries of circumstellar disks can have a large variety (Dullemond & Dominik 2004b). The setup is axially symmetric. Hence, the radiation transport is done in 2D. The gas density distribution of the circumstellar disk is given in cylindrical coordinates (r,z) by $\mathit{\rho }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathit{,}$(1)with the radially and vertically dependent functions ${\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{r}}_{\mathrm{d}}\mathit{/}\mathit{r}\mathrm{and}{\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}\exp \left(\mathrm{-}\frac{\mathit{\pi }}{\mathrm{4}}{\left(\frac{\mathit{z}}{\mathit{h}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\right)}^{\mathrm{2}}\right)$(2)making use of the abbreviations h(r) = z_d(r/r_d)^1.125, r_d = 500 AU, and z_d = 125 AU. The normalization ρ₀ of the density setup is chosen to define different optical depths τ through the midplane of the corresponding circumstellar disk. We compute problems scanning seven orders of magnitude in optical depth.

Fig. 1 Visualization of the disk configuration and final temperature for the case τ550 nm = 100. Solid lines denote iso-density contours. Colors denote the final temperature distribution.

The central star is assumed to radiate as a black body. In the Pascucci et al. (2004) radiation benchmark test (τ_550 nm ≤ 100), the star has a radius of R_∗ = 1 R_⊙ and a temperature of T_∗ = 5800 K. In the Pinte et al. (2009) radiation benchmark test (τ_810 nm ≥ 1000), the star has a radius of R_∗ = 2 R_⊙ and a temperature of T_∗ = 4000 K. A visualization of the disk setup and the final temperature distribution is given in Fig. 1.

The opacity table used is the same as in the original benchmark test of Pascucci et al. (2004) taken from Draine & Lee (1984). These opacities are shown in Fig. 2. The dust-to-gas mass ratio is set constant to 1%. Scattering is neglected.

Fig. 2 Frequency-dependent opacities from Draine & Lee (1984). The solid line denotes an exponential fit to the long wavelength regime.

3.2. Numerical configuration

In the Pascucci et al. (2004) radiation benchmark test (τ_550 nm ≤ 100), the computational domain goes from R_min = 1 AU up to R_max = 1000 AU. In the Pinte et al. (2009) radiation benchmark test (τ_810 nm > 1000), the computational domain goes from 0.1 AU up to 400 AU. In both setups, the polar angle covers the full domain from 0 to 180 degrees.

The boundary conditions in the polar direction are chosen so as not to allow any radiative flux over the polar axis. The boundary condition at the outer radial direction is chosen to allow for a radiative flux out of the computational domain computed in the optically thin limit. The boundary condition for the thermal radiation field at the inner radial direction is in the RT simulations given by a low value of the thermal radiation field (the total radiation field at the inner disk rim is dominated by the irradiation). The boundary condition for the thermal radiation field at the inner radial direction is in the simulations using only FLD determined from the surface temperature of the central star assuming an optically thin limit between the stellar surface and the inner disk rim (as in the original benchmark tests).

The grids in spherical coordinates are chosen to resolve the local optical depths of the circumstellar disks in the polar direction in order to guarantee the computation of correct cooling properties. For all cases up to τ_810 nm = 1.22 × 10^+ 3, the grid consists of n_r × n_θ = 128 × 360 cells. For τ_810 nm = 1.22 × 10^+ 4, the grid has n_r × n_θ = 64 × 2880 grid cells. For τ_810 nm = 1.22 × 10⁺⁶, the grid has n_r × n_θ = 16 × 10⁵ grid cells.

We note that results of the same accuracy are obtained by also using 128 × 360 grid cells for the higher optical depths runs. However, they make use of a numerical “trick” similar to the one described in Pinte et al. (2009) for the ProDiMo code, where the density distribution was altered to limit the optical depth. Here, we promote a modified algorithm of the same idea: the resolution of a 128 × 360 grid is sufficient also for high optical depth runs, if an upper limit of the optical depth of τ_max ≤ 1 is applied to each individual grid cell during the radiation transport. Although we focus in the following discussion on the results of the high resolution grids only, a solution of a simulation run using this numerical trick is shown for comparison for the most optically thick case of τ_810 nm = 1.22 × 10⁶ in Fig. 7.

The numerical configuration of the Monte-Carlo code RADMC is given in Pascucci et al. (2004) and Pinte et al. (2009) with the difference that we do not consider scattering here.

3.3. Runs performed

We computed the equilibrium temperature distribution in the disk and envelope using different radiation transport methods. We studied the gray FLD approximation, the gray irradiation + gray FLD approximation, and the frequency-dependent irradiation + gray FLD approximation. As a reference solution, we used results obtained with the Monte-Carlo code RADMC. We performed simulations for various optical depths of the disk’s midplane from τ_550 nm = 0.1 up to τ_810 nm = 1.22 × 10⁺⁶. An overview of runs performed is given in Table 1.

4. Simulation results

We compared the resulting temperature profiles through the disk’s midplane for the different radiation transport methods depending on the optical depth of the circumstellar disk. The final accuracy of the different approximate radiation transport methods was computed as the deviation to the comparison code results by RADMC. These deviations and the standard deviation of the comparison code with respect to the other radiation transport codes participating in the original benchmark tests of Pascucci et al. (2004) and Pinte et al. (2009) are given in the last column of Table 1. Visualizations of the individual temperature profiles and the deviations to the reference code are given in Figs. 3 to 7.

Fig. 3 Temperature profiles (upper panel) in the midplane of the circumstellar disk for the case of low optical depth τ550 nm = 0.1. In the lower panel the deviations of the three different radiation transport methods to the comparison code result are displayed. Solid lines denote results for the frequency-dependent RT + gray FLD method. Dotted lines denote results for the gray RT + gray FLD method. Dashed lines denote results for the gray FLD method. Pluses “+” denote results for the comparison code RADMC. In this optical thin case, the results for gray RT + gray FLD are almost identical to the frequency-dependent RT + gray FLD. Here, analytical estimates for the optically thin approximation are given: crosses “×” mark the analytical estimate by Spitzer (1978) for irradiated regions far away from the star, circles “o” mark the analytical estimate in the gray and isotropic approximation (T ∝ r−1/2).

Fig. 4 Same as Fig. 3 for simulation runs with τ550 nm = 102.

Fig. 5 Same as Fig. 3 for simulation runs with τ810 nm = 1.22 × 103. In this fairly optical thick case, the results for gray RT + gray FLD are almost identical to the frequency-dependent RT + gray FLD, small differences are visible at the innermost rim of the disk only.

Fig. 6 Same as Fig. 3 for simulation runs with τ810nm = 1.22 × 104. In this highly optical thick case, the results for gray RT + gray FLD are identical to the frequency-dependent RT + gray FLD.

Fig. 7 Same as Fig. 3 for simulation runs with τ810 nm = 1.22 × 106. In this highly optical thick case, the results for gray RT + gray FLD are identical to the frequency-dependent RT + gray FLD. Here, the dotted lines represent the results for the frequency-dependent RT + gray FLD method in a low-resolution simulation if an upper limit of the optical depth per grid cell of τmax ≤ 1 is used (see Sect. 3.2 for a description of this numerical trick).

5. Discussion

In the following, we discuss the simulation results and their implications, depending on the different regimes of optical depth. Afterwards, we discuss the validity of the various approximations used in the different radiation transport methods described above.

5.1. Different regimes of optical depth

We begin our discussion with the two limiting cases of very optically thin and thick media. The regime of moderate optical depth is addressed afterwards.

5.1.1. The most optically thin regime

In the optically thin limit, the temperature profile can be accurately estimated analytically. As given in Spitzer (1978), the local dust temperature in irradiated regions far away from the central star (r ≫ R_∗) is presented by $\mathit{T}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\left(\frac{{\mathit{R}}_{\min }}{\mathrm{2}\mathit{r}}\right)}^{\frac{\mathrm{2}}{\mathrm{4}\mathrm{+}\mathit{\beta }}}{\mathit{T}}_{\min }\mathit{,}$(3)with the exponent β of the opacity power law in the long wavelength regime and the temperature T_min at the inner rim of the computational domain at R_min. Usually, β can be approximated as β ≈ 2. The exponent of the opacity power law in the long wavelength regime (Draine & Lee 1984) is determined to be β = 2.0508 (see Fig. 2). In the limiting case of gray and isotropic (i.e., non-irradiation) radiation fields, the local dust temperature is expected to follow $\mathit{T}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\sqrt{\frac{{\mathit{R}}_{\min }}{\mathit{r}}}{\mathit{T}}_{\min }\mathit{.}$(4)The analytical estimate (3)by Spitzer (1978) and the gray, isotropic relation (4)are shown together with the simulation results for the most optically thin case in Fig. 3.

As expected, the gray FLD approximation resembles the slope of the gray isotropic assumption, while the frequency-dependent and gray hybrid methods resemble the Spitzer estimate for irradiated regions. In this limiting case of very low optical depth, the frequency-dependent and gray irradiation routines give identical results.

5.1.2. The most optically thick regime

In general, configurations including regions of higher optical depth in otherwise optically thin environments denote the numerically most challenging problems. Typically very accurate radiation transport methods such as RT and Monte-Carlo have to consume the majority of cpu time to reach a converged temperature in the optically thick regions. While photon propagation in the optically thick regions alone is most easily described in the diffusion limit, the (flux-limited) diffusion approximation fails to account for shadowing effects in multi-dimensional problems due to the loss of angular information in its derivation.

The simulation results for the most optically thick case are shown in Fig. 7. The RT + gray FLD runs are in overall good agreement with the comparison results from RADMC. The maximum deviation of 43% corresponds to the maximum difference between the various radiation transport codes of the original benchmark test (Pinte et al. 2009). As expected in this limiting case of very high optical depth, the frequency-dependent and gray irradiation routines give identical results.

The gray FLD approximation fails to reproduce the shadow behind the optically thick inner disk rim. In this method, the stellar radiative flux diffuses around this obstacle and heats the disk through the optically thin atmosphere at larger radii. The resulting midplane temperature is up to a factor of 2.8 higher than in the comparison run.

5.1.3. The regime of moderate optical depth

The frequency-dependent and gray irradiation routines differ only in the case of moderate optical depth, shown in Fig. 4. In this intermediate regime, the disk’s midplane is optically thick for the higher frequency part of the stellar irradiation spectrum, but optically thin in the long wavelength regime. Therefore, only the frequency-dependent irradiation routine is able to resemble the heating of the irradiated midplane at larger radii correctly. In the case of gray irradiation, the long wavelength flux is already absorbed at the inner disk rim and the temperature at larger radii follows the slope of the gray FLD method.

5.2. Validity of approximations

5.2.1. The gray flux-limited diffusion approximation for thermal dust emission

In all our approximate radiation transport simulations, the thermal (re-)emission of dust grains is computed in the gray FLD approximation. This approximation reduces cpu time significantly (up to several orders of magnitude for a high-resolution 3D scheme) to allow the usage of the radiation transport method within (magneto-) hydrodynamics simulations. If the stellar irradiation is accounted for within a frequency-dependent RT step, the results of the approximate radiation transport simulations are in accordance with the deviations given by the differences from the various high-level radiation transport methods from the original radiation benchmark tests (Pascucci et al. 2004; Pinte et al. 2009). Hence, we conclude that the gray FLD approximation is an efficient and accurate method for solving the thermal dust (re-)emission of circumstellar disks. In general, the FLD method suffers from the issues explained in Sect. 5.2.3.

5.2.2. The gray approximation for the stellar irradiation

The frequency-dependent treatment of the stellar irradiation spectrum improves the results for configurations of moderate optical depth (for the particular test case of a Sun-like star and τ_550 nm ~ 100), in which the disk’s midplane is optically thick for the high-frequency part of the stellar spectrum, while it is optically thin for a substantial part of the long wavelength regime. In the limiting cases of very low and very high optical depth, the gray irradiation is in agreement with the frequency-dependent runs.

5.2.3. The flux-limited diffusion approximation for the stellar source term

In general, treating the stellar source term within the FLD approximation involves two major problems. First, the FLD approximation assumes a locally isotropic radiation field. This assumption yields the wrong radial slope of the temperature even in highly optically thin regions, as depicted in Fig. 3. Dust grains in an optically thin region around a star absorb photons mostly from the direction towards the star. Hence, this scenario is highly anisotropic.

Second, in the derivation of the FLD approximation, the momentum expansion of the radiation transport equation is closed by relating the radiative flux to the gradient of the radiation energy. This assumption prevents the FLD approximation from accounting for shadows behind radiative obstacles, such as the optically thick inner rim of circumstellar disks. The radiative flux in FLD simulations tends to diffuse from the stellar surface around the optically thick inner disk rim and heats the disk midplane at larger radii to unrealistically high values.

Hence, we conclude that the FLD approximation results in large errors if used for (stellar) irradiation sources. In radiation-(magneto-) hydrodynamic studies of accretion disks, the strong increase of the midplane temperature would result, for example, in an overestimate of the disk stability based on the Toomre criterion (Toomre 1964).

5.2.4. The hybrid approximate radiation transport method

In general, the hybrid approximate radiation transport method using a frequency-dependent RT step for stellar irradiation and a gray FLD solver for thermal dust emission results in high accuracy for the whole range of optical depths. By using only a small fraction of the cpu time needed for the more sophisticated comparison code, this method fully agrees with cases of optical depths τ_550 nm ≤ 100 as well as with the high optical depth cases. Only in the intermediate regime of a total optical depth of about τ_810 nm ~ 1000 does the approximate method overestimate the temperature in the disk midplane behind the most optically thick inner rim by about 48%, compared to only 20% deviation of the various original benchmark codes in Pinte et al. (2009).

6. Summary

We checked the reliability of approximate radiation transport methods using flared disk configurations of standard radiation benchmark tests. We compared the resulting temperature distributions for a wide range of optical depths from a fairly optically thin case (τ_550 nm = 0.1) up to a highly optically thick case (τ_810 nm = 1.22 × 10⁺⁶). The outcome for the various methods and different regimes of optical depth are summarized in Fig. 8.

Fig. 8 Maximum temperature deviation in the disk’s midplane as function of the optical depth of the circumstellar disk. The solid line denotes the results for the frequency-dependent RT + gray FLD method. The dotted line denotes the results for the gray RT + gray FLD method. The dashed line denotes the results for the gray FLD method. The pluses “+” denote the deviations of the results of the various high-level radiation transport codes used in the original benchmark tests of Pascucci et al. (2004) and Pinte et al. (2009). These values are read off Figs. 4 and 5 in Pascucci et al. (2004) and Fig. 10 in Pinte et al. (2009).

Approximating the stellar irradiation spectrum by a gray Planck opacity is in accordance with a frequency-dependent treatment in the optically thin and thick limit. For intermediate optical depths, a gray approximation of the stellar irradiation spectrum yields a slightly hotter inner rim and a slightly cooler midplane of the disk at larger radii.

The gray FLD approximation fails to compute an appropriate temperature profile in all regimes of optical depth; the maximum deviations to the comparison runs are 50% in the optically thin and up to 280% in the optically thick limit. For low optical depth, the isotropic assumption within the FLD method yields a too steep decrease of the radial temperature slope. For higher optical depths, the FLD approximation does not render the shadow behind the optically thick inner rim of the circumstellar disk, yielding artificial heating at larger disk radii. This issue of neglected shadows is expected to be even more problematic for increasing stellar luminosity.

The frequency-dependent RT + gray FLD approximation denotes a highly accurate treatment of the radiation field in the optically thin and thick regime. Only in the intermediate regime of moderate optical depth τ_810 nm ~ 1000 are the deviations of the radiation transport codes from the original benchmark tests by Pinte et al. (2009) lower than for the hybrid approximate radiation transport method. The high accuracy of the hybrid approximate radiation transport is attended by very low computational costs compared with the very expensive Monte-Carlo calculations.

Hence, we close this comparison study by pointing out that the frequency-dependent RT + gray FLD method is ideally suited for follow-up (magneto-) hydrodynamical studies of circumstellar accretion disks.

Acknowledgments

We thank Cornelis Dullemond for providing the RADMC reference simulation data and for fruitful discussions of approximate radiation transport solvers in the context of passively irradiated circumstellar disks. Author R.K. acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG) via the collaborative research project SFB 881 “The Milky Way System” in sub-project B2. Both authors thank the DFG for support via the priority program SPP 1573 “Physics of the ISM” as well as the Baden-Württemberg Foundation for support via contract research (grant P- LS-SPII/18) via their program Internationale Spitzenforschung II.

References

All Tables

All Figures

	Fig. 1 Visualization of the disk configuration and final temperature for the case τ550 nm = 100. Solid lines denote iso-density contours. Colors denote the final temperature distribution.
In the text

	Fig. 2 Frequency-dependent opacities from Draine & Lee (1984). The solid line denotes an exponential fit to the long wavelength regime.
In the text

Fig. 3 Temperature profiles (upper panel) in the midplane of the circumstellar disk for the case of low optical depth τ550 nm = 0.1. In the lower panel the deviations of the three different radiation transport methods to the comparison code result are displayed. Solid lines denote results for the frequency-dependent RT + gray FLD method. Dotted lines denote results for the gray RT + gray FLD method. Dashed lines denote results for the gray FLD method. Pluses “+” denote results for the comparison code RADMC. In this optical thin case, the results for gray RT + gray FLD are almost identical to the frequency-dependent RT + gray FLD. Here, analytical estimates for the optically thin approximation are given: crosses “×” mark the analytical estimate by Spitzer (1978) for irradiated regions far away from the star, circles “o” mark the analytical estimate in the gray and isotropic approximation (T ∝ r−1/2).

In the text

	Fig. 4 Same as Fig. 3 for simulation runs with τ550 nm = 102.
In the text

	Fig. 5 Same as Fig. 3 for simulation runs with τ810 nm = 1.22 × 103. In this fairly optical thick case, the results for gray RT + gray FLD are almost identical to the frequency-dependent RT + gray FLD, small differences are visible at the innermost rim of the disk only.
In the text

	Fig. 6 Same as Fig. 3 for simulation runs with τ810nm = 1.22 × 104. In this highly optical thick case, the results for gray RT + gray FLD are identical to the frequency-dependent RT + gray FLD.
In the text

Fig. 7 Same as Fig. 3 for simulation runs with τ810 nm = 1.22 × 106. In this highly optical thick case, the results for gray RT + gray FLD are identical to the frequency-dependent RT + gray FLD. Here, the dotted lines represent the results for the frequency-dependent RT + gray FLD method in a low-resolution simulation if an upper limit of the optical depth per grid cell of τmax ≤ 1 is used (see Sect. 3.2 for a description of this numerical trick).

In the text

Fig. 8 Maximum temperature deviation in the disk’s midplane as function of the optical depth of the circumstellar disk. The solid line denotes the results for the frequency-dependent RT + gray FLD method. The dotted line denotes the results for the gray RT + gray FLD method. The dashed line denotes the results for the gray FLD method. The pluses “+” denote the deviations of the results of the various high-level radiation transport codes used in the original benchmark tests of Pascucci et al. (2004) and Pinte et al. (2009). These values are read off Figs. 4 and 5 in Pascucci et al. (2004) and Fig. 10 in Pinte et al. (2009).

In the text