---

Benchmark problems for continuum radiative transfer

High optical depths, anisotropic scattering, and polarisation

C. Pinte¹ - T. J. Harries¹ - M. Min² - A. M. Watson³ - C. P. Dullemond⁴ - P. Woitke^5,6 - F. Ménard⁷ - M. C. Durán-Rojas³

1 - School of Physics, University of Exeter, Stocker Road, Exeter EX4 4QL, UK
2 - Sterrenkundig Instituut Anton Pannekoek, Universiteit Amsterdam, Kruislaan 403, 1098 Amsterdam, The Netherlands
3 - Centro de Radioastronomía y Astrofsíca, Universidad Nacional Autóma de México, Morelia, Mich., Mexico
4 - Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
5 - UK Astronomy Technology Centre, Royal Observatory, Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK
6 - School of Physics & Astronomy, University of St. Andrews, North Haugh, St. Andrews KY16 9SS, UK
7 - Laboratoire d'Astrophysique de Grenoble, CNRS/UJF UMR 5571, 414 rue de la Piscine, BP 53, 38041 Grenoble Cedex 9, France

Received 19 December 2008 / Accepted 22 February 2009

Abstract
Aims. Solving the continuum radiative transfer equation in high opacity media requires sophisticated numerical tools. In order to test the reliability of such tools, we present a benchmark of radiative transfer codes in a 2D disc configuration.
Methods. We test the accuracy of seven independently developed radiative transfer codes by comparing the temperature structures, spectral energy distributions, scattered light images, and linear polarisation maps that each model predicts for a variety of disc opacities and viewing angles. The test cases have been chosen to be numerically challenging, with midplane optical depths up 10⁶, a sharp density transition at the inner edge and complex scattering matrices. We also review recent progress in the implementation of the Monte Carlo method that allow an efficient solution to these kinds of problems and discuss the advantages and limitations of Monte Carlo codes compared to those of discrete ordinate codes.
Results. For each of the test cases, the predicted results from the radiative transfer codes are within good agreement. The results indicate that these codes can be confidently used to interpret present and future observations of protoplanetary discs.

Key words: radiative transfer - circumstellar matter - accretion, accretion disks - planetary systems: protoplanetary disks - methods: numerical

1 Introduction

Dust represents an essential element in the energy balance of a variety of astrophysical objects, from the interstellar medium to the atmospheres and close circumstellar environments of numerous classes of object; from the lower mass planets and brown dwarfs, to massive stars. With the advent of high-angular resolution and high-contrast instruments, the basic structural properties (e.g., size, inclination, and surface brightness) of the circumstellar environments of the nearest and/or largest objects - discs and envelopes around young stars in nearby star-forming regions and around more distant evolved stars - are now under close scrutiny. With this unprecedented wealth of high-resolution data, from optical to radio, detailed studies of the dust content become possible and sophisticated radiative transfer (RT) codes are needed to fully exploit the data.

At short wavelengths, dust grains efficiently absorb, scatter, and polarise the starlight while at longer wavelengths dust re-emits the absorbed radiation. How much radiation is scattered and absorbed is a function of both the geometry of the circumstellar environment and the properties of the dust. In turn, the amount of absorbed radiation sets the temperature of the dust (and gas) and defines the amount of radiation that is re-emitted at longer, thermal wavelengths.

To get a reliable understanding of the structure and evolution of these ``dusty'' objects, be it the evolution of dust grain sizes, the temperature dependent chemistry, or simply the density profiles, it is highly desirable to model not only the integrated fluxes (i.e., the spectral energy distributions, hereafter SED), but also the resolved brightness and/or polarisation profiles when available. This can only be done by solving the radiative transfer (hereafter RT) problem in media that can have large optical depths and/or complex geometries and compositions. Recent studies of circumstellar discs are based on detailed comparisons of high-quality data sets, combining various kinds of observation (SED, multiple wavelength scattered light images, polarisation map, infrared or millimetre visibilities) to the predictions of RT codes (e.g. Pinte et al. 2008a; Doucet et al. 2007; Steinacker et al. 2006; Tannirkulam et al. 2008; Watson & Stapelfeldt 2004; Pinte et al. 2008b; Wolf et al. 2003; Pontoppidan et al. 2007; D'Alessio et al. 2006; Pinte et al. 2007; Glauser et al. 2008; Wood et al. 2002; Fitzgerald et al. 2007). Such studies will become more and more common with the advent of new instruments like VLT/SPHERE, Gemini/GPI, JWST, Herschel and ALMA, and validating the accuracy of RT codes is of particular importance.

Analytical solutions do not exist for wavelength-dependent radiative transfer and sophisticated numerical methods must be used. Testing the reliability of RT computations requires in that case to compare the solutions to well-defined problems by independent codes. Such a work has been done by Ivezic et al. (1997) for a 1D spherical geometry and by Pascucci et al. (2004, hereafter P04) in a 2D disc configuration. The later work compared in detail the calculations of five radiative codes and has been used as a reference to validate newly developed RT codes (e.g. Ercolano et al. 2005; Pinte et al. 2006; Harries et al. 2004). The test cases were however limited to relatively modest optical depths (midplane opacity in the V band), orders of magnitude smaller than the actual optical depths of protoplanetary discs for which radiative transfer codes are generally used in the literature. High optical depths represent challenging calculations for radiative transfer codes, with potential convergence issues, and additional tests are required to confidently trust the results of radiative codes in this regime. Furthermore, calculations in P04 were done assuming isotropic scattering and restricted to SEDs. With the advent of high-resolution observations of young stellar objects, validating the calculations of resolved surface brightness of these objects is now also crucially needed.

In this paper, we extend the work of P04 to realistic, very optically thick discs with anisotropic scattering. We perform a comparison of seven RT codes in a well defined 2-dimensional disc configuration, with simple dust properties. The prediction of four Monte Carlo codes (temperature structure in the disc, the emergent SED as well as monochromatic scattered light images, and polarisation maps for different disc opacities and viewing angles) are compared in the case of anisotropic scattering (Sects. 2 to 4). Additional comparisons with discrete ordinate codes, in the case of isotropic scattering and without scattering, are presented in Sect. 5.

2 Radiative transfer modelling

2.1 The radiative transfer problem

Solving the RT problem in dusty environments aims at determining the (polarised) specific intensity at each point and direction of the volume and at each wavelength . This intensity is obtained by solving the stationary transfer equation.

In the case of randomly oriented dust particles, the radiative transfer equation can be written adopting the Stokes formalism:

where is the Stokes vector, with I representing the total intensity, Q and U the linearly polarised intensities, and V the circularly polarised intensity. , and are the absorption, scattering and extinction opacities, respectively. s is the length along the direction of propagation. is the scattering (or Mueller) matrix describing the changes in the Stokes vector when the light is scattered from the direction to the direction . is the Planck function and is the unitary Stokes vector representing unpolarised emission .

Computation of the thermal emission requires to determine the dust temperature structure . This temperature is determined by writing that the dust is in radiative equilibrium. If we assume that the dust is at the local thermodynamic equilibrium and that there is no more sources of energy than the radiation field, the temperature is obtained by solving the implicit equation:

where is the mean specific intensity (i.e. the specific intensity averaged over all solid angles).

The system of Eqs. (1) and (2) completely defines the RT problem when the dust optical properties ( , , ) and sources of radiation (initial conditions for Eq. (1)) are given. It is important to note that opacities can depend on the temperature, in which case solving simultaneous Eqs. (1) and (2) requires an iterative scheme. Most of the time however, dust opacities do not vary much with temperature and can be assumed to be constant. We will make this assumption in the following analysis.

Additionally, radiative transfer plays an integral role in the physics of the disc. It alters the density structure via hydrostatic equilibrium (e.g. Walker et al. 2004) and impacts the dust content, and hence the opacity. For instance differential dust sublimation at the inner edge can destroy some of the dust grains, the temperature of a dust grain depending on its size (e.g., Tannirkulam et al. 2007) and composition (e.g., Woitke 2006). Similarly the formation of ice mantles around the grains in the cold, outer regions of the disc also affects the grain opacities. Including any of these effects requires an iterative approach. In this paper, we restrict ourselves to the benchmark of radiative transfer solvers and keep the density structure and dust properties fixed.

2.2 Numerical methods

2.2.1 The Monte Carlo method

Anisotropic scattering by dust grains precludes the use of direct methods to solve the continuum RT equation and Monte Carlo methods are commonly used instead. They solve the RT equation by stochastically propagating ``photon packets'' through the dusty environment. The transport of packets is governed by scattering, absorption and re-emission events that are controlled by the optical properties of the medium (opacity, albedo, scattering phase function, etc) and by the temperature distribution. Upon leaving the model boundaries, ``photon packets'' are used to build an SED and/or synthetic images.

The Monte Carlo scheme estimates physical quantities by statistical means, which potentially leads to noisy results when the number of packets sampling some regions and/or directions in the model becomes low. Several variance reduction techniques have been developed to improve the sampling of the Monte Carlo method: forced first scattering (Cashwell & Everett 1959), peel-off techniques (Yusef-Zadeh et al. 1984), estimation of the mean specific intensity (Lucy 1999), immediate re-emission (Bjorkman & Wood 2001), and importance weighting schemes (Juvela 2005). These various techniques allowed the Monte Carlo method to progressively become competitive against grid-based methods, and it is now more and more commonly used to solve the continuum RT problem.

Despite these techniques, Monte Carlo methods can become computationally expensive when the optical depth becomes very large. In our disc configuration, this leads to two major difficulties:

• because gradients of opacity are oriented toward the disc surface, packets entering the disc tend to escape after a few interactions. Very few packets penetrate the central regions of the disc, leading to a noisy temperature structure close to the disc midplane;
• for edge-on configurations, the flux in the near- and mid-infrared is dominated by packets originating from the dense central regions of the disc. These packets will be scattered towards the observer in the surface layers of the disc, where the probability of interaction is very low (due the very low density). Packets escaping the dense regions (after a large number of interactions) almost always go through the surface layers without interacting with the dust grains and large number of packets are required to converge SEDs or images in this regime.

In the next paragraphs, we describe two schemes that, when coupled with a Monte Carlo approach, significantly improve the efficiency of the codes. They have helped to overcome the previously mentioned difficulties and to efficiently solve the test cases presented in this paper.

2.2.2 Diffusion approximation

In the deep regions of the disc, solving the complete RT equation is not necessary. The radiation field becomes isotropic and the source function becomes equal to the Planck function. The behaviour of the density of energy is in that case properly described by the diffusion theory:

where the diffusion coefficient is defined as with the local Rosseland opacity.

In this case, Monte Carlo methods can be efficiently coupled to diffusion approximation methods. The model can be divided into two regions. The first one, that corresponds the surface of the disc represents all parts of the model volume where the optical depth in any direction is smaller than a given threshold. In this region, the temperature structure is computed with a Monte Carlo method, eventually including acceleration schemes such as the immediate re-emission concept and/or estimation of the mean specific intensity.

In the rest of the model volume, i.e. in the central regions of the disc, the temperature structure can be solved using the diffusion approximation. The temperature structure at the edge of the diffusion approximation region (initial condition for Eq. (3)) is given by the solution of the Monte Carlo calculations.

The optical depth threshold which defines these two regions must be set high enough to ensure that the radiation field inside the diffusion approximation region is isotropic and dominated by the local emission.

For an optimal efficiency, this method must avoid calculating the complete propagation of photon packets inside the diffusion approximation region. This can be done by using various methods. The propagation of the packets can be calculated in a faster way by using a modified random walk procedure (Min et al. 2009; Fleck & Canfield 1984). This method combines multiple interaction steps in one computation, while still keeping track of the energy deposited in an accurate sense. In this way, the interaction between the optically thick regions and the upper layers of the disc is properly computed. In addition, the number of steps taken can be easily adjusted to the local radiation field and density gradient, making it a highly flexible method. A ``mirror'' condition can also be used: when a packet enters the diffusion approximation region, it is sent back with the same energy and wavelength but with an opposite direction vector. Although not rigorously exact, this method (used by TORUS and MCFOST) was found to be very accurate when compared to the full Monte Carlo solution (see Sect. 4.1).

Details and tests on the accuracy of diffusion approximation methods are presented in Min et al. (2009).

2.2.3 Ray-tracing

When the specific intensity is known in each point of the model, direct ray-tracing methods using the formal solution of the RT equation can be used to produce observables.

Ray-tracing has for instance been used in combination with Monte Carlo methods to produce SEDs and emission maps in the infrared and millimetre regimes, where scattering can be considered isotropic in some cases (Dullemond & Dominik 2004; Wolf 2003). When the scattering is isotropic, only the mean specific intensity and temperature structure are required to calculate the source function. These two quantities can be easily estimated with a Monte Carlo method (Lucy 1999) and do not require large amount of memory to be stored. After an initial Monte Carlo run computing the total mean intensity and temperature structure in the disc, SEDs and/or maps can then be produced by integrating the source function on rays originating from the observer. Which such a method, the Monte Carlo method is used only to estimate the specific intensity and not the images and/or SEDs. The resulting noise in the observables is much lower as it only reflects the noise in the mean specific intensity and no longer the noise associated to the production of the observables themselves.

This method combining Monte Carlo and ray-tracing can be extended to any wavelength if the angular dependence of the scattering component of the source function is preserved in the calculations. The Monte Carlo method produces all the information needed to perform such calculations, as it can give an estimate of the specific intensity, with its complete angular dependence, and not only of the mean specific intensity. However, storing the full spatial, angular and wavelength dependence of the radiation field requires large amounts of memory which is currently beyond computational capacities.

This difficulty can be overcome by invoking successive monochromatic Monte Carlo runs, which removes the need to store the wavelength dependence of the specific intensity. An initial multi-wavelength Monte Carlo run calculates the temperature structure in the disc. The SED is then constructed wavelength by wavelength with successive monochromatic Monte Carlo runs that estimate the specific intensity at each point of the model.

In MCFOST, the specific intensity is then saved for a set of angular directions (method 1). At the end of each Monte Carlo run, the scattering emissivity in any direction is calculated from the specific intensity and resolved maps and/or integrated fluxes for any inclinations are finally obtained by ray-tracing. This step (monochromatic Monte Carlo run + ray-tracing run) is repeated over all wavelengths, without storing the specific intensity at the previous wavelength.

A slightly different method is adopted in MCMax, where the scattering emissivity in a given set of directions is stored instead of the specific intensity itself (method 2). The scattering emissivity is the product of the specific intensity by the scattering phase function, i.e. the last term in Eq. (1). Each time a packet crosses a cell, its contributions to the scattering emissivity in the chosen directions are calculated by multiplying the packet energy by the local phase function. At the end of each monochromatic Monte Carlo run, maps and/or fluxes at the chosen inclinations are produced by integrating the source function via a ray-tracing method.

Method 1 avoids the expensive calculations of the scattering emissivity each time a packet crosses a cell but requires a larger amount of memory to store the angular dependence of the specific intensity. If the radiation field is stored for a few specific wavelengths, it also allows to produce scattered light images and emission maps at any inclinations, by only running additional ray-tracing calculations. However, with this method, the angular sampling of the radiation field must be performed with care, especially when scattering is very anisotropic. This issue is not encountered with method 2, which has the same, almost perfect, angular sampling of the radiation field as classical Monte Carlo methods.

2.3 Codes description

2.3.1 MCFOST

MCFOST is a 3D continuum and line radiative transfer code based on the Monte Carlo method (Pinte et al. 2006). Temperature structures are calculated using the immediate re-emission concept of Bjorkman & Wood (2001) but with a continuous deposition of energy to estimate the mean intensity (Lucy 1999). The code uses a spherical or cylindrical grid, with an adaptive mesh refinement at the inner edge (based on the opacity gradient) so as to properly sample the inner radius of the disc.

Several improvements have been implemented on top of the original algorithm presented in Pinte et al. (2006). In very optically thick parts of the model, the temperature structure is calculated with a diffusion approximation method, using the Monte Carlo calculations as limit conditions. Equation (3) is solved as an asymptotic limit of the time dependent diffusion equation, via an implicit scheme to ensure stability and accuracy. The transition between the Monte Carlo and diffusion approximation regions is set to an optical depth of 1 000 at the wavelength where the stellar emission peaks. To avoid calculating the propagation of photon packets inside the diffusion approximation domain, packets are mirrored at the boundaries of the Monte Carlo domain.

The temperature structure and radiation field estimated by the Monte Carlo runs are used to produce images, polarisation maps, and SEDs with a ray-tracing method, where the emerging flux is obtained by calculating the formal solution of the radiative equation along rays.

2.3.2 MCMax

The Monte Carlo radiative transfer code MCMax is based on the scheme of immediate re-emission as proposed by Bjorkman & Wood (2001). For the temperature structure the method of continuous absorption by Lucy (1999) is implemented. The photons are traced in 3D on a spherical coordinate grid while the geometry of the system is set to be cylindrically symmetric. For the optically thick regions a modified random walk procedure is applied in order to make multiple interaction steps in a single computation (see Min et al. 2009; Fleck & Canfield 1984). This has the advantage that the computational speed is increased significantly while the temperature structure is still computed with high accuracy. After the Monte Carlo procedure a partial diffusion approximation is used for these regions in the disc that received too few photons to determine a reliable temperature structure (see Min et al. 2009). All the observables are constructed by integrating the formal solution using ray-tracing. In this way noise on the observables is reduced significantly.

The spatial grid at the inner edge of the disc is set in such a way that the optical depth for both the local radiation field and the stellar radiation are sampled logarithmically.

2.3.3 Pinball

Pinball is a Monte-Carlo code that calculates scattered-light images; it does not calculate the equilibrium temperature or include thermal re-emission. An earlier version and a pair of simple test cases were described by Watson & Henney (2001). The current version includes polarisation.

2.3.4 TORUS

TORUS is a 3D continuum and line radiative transfer code based on the Monte Carlo method (Harries 2000; Kurosawa et al. 2004; Harries et al. 2004). Radiative equilibrium is computed using the continuous absorption algorithm from Lucy (1999). Calculations are performed on a 2D, cylindrical adaptive-mesh grid. Storing the opacity information on an adaptive mesh has particular advantages for the problem considered here, since it allows an adequate sampling of the inner edge of the disc, where the opacity gradient is very steep. The temperature structure in the central regions of the disc is computed with a diffusion approximation method.

For scattered light images, the enforced scattering concept (Cashwell & Everett 1959) as well as the peel-off technique (Yusef-Zadeh et al. 1984) are implemented to reduce the variance.

The descriptions of the codes ProDiMo, RADMC and RADICAL, that have performed the test cases with isotropic scattering and without scattering, are presented in Sect. 5.1.

3 Benchmark problem

All the codes in this paper that calculate the thermal equilibrium have successfully reproduced the P04 benchmark, from optically thin configurations to optical depths of 100 in the optical. The test cases presented here are complementary to those in P04 and are restricted to optical depths higher than 1 000. The full description of the benchmark problem, including tabulated values for the disc density and the dust properties, as well as the results for all codes are presented on the webpage http://www.astro.ex.ac.uk/people/cpinte/benchmark/. This should allow additional codes to compare their results with the ones presented in this paper.

3.1 System geometry

The model consists of a dusty disc surrounding a central star and located at a distance of 140 pc.

We consider an axisymmetric flared density structure with a Gaussian vertical profile . We use power-law distributions for the surface density and the scale height h(r) = h₀  (r/r₀)^1.125 where r is the radial coordinate in the equatorial plane and h₀ = 10 AU is the scale height at the radius r₀ =100 AU. The disc extends from an inner cylindrical radius  AU to an outer limit  AU. The edges of the disc are assumed to be sharp, i.e. vertical: there is nothing inside and outside and the density is defined by the previously mentioned power-laws between them. The dust disc mass is the only parameter varied and takes 4 different values: , , and .

This configuration was chosen because it represents a more difficult problem to solve than the test case presented by P04: the disc extends much closer to the star, 0.1 AU instead of 1 AU and the radial gradient of density is much steeper with a slope of surface density of -1.5 instead of 0.125, leading to much higher densities close to the inner edge of the disc, and hence much higher disc optical depths.

The star is defined as a uniformly radiating blackbody sphere at a temperature of 4 000 K and with a radius of 2 solar radii.

3.2 Dust properties

Dust grains are defined as homogeneous and spherical particles with a single size of 1 m and are composed of astronomical silicates (Weingartner & Draine 2001). The grain mass density is fixed to 3.5 g/cm³.

The dust optical properties: extinction and scattering opacities (Fig. 1), scattering phase functions, and Mueller matrices are calculated using the Mie theory. The resulting midplane optical depth in I band (0.81 m), from the star to the observer, is ranging from to when the disc mass varies from to . For simplicity, in the following, we will label the different models , , and .

In the simplifying case of Mie scattering, the matrix becomes block-diagonal with only 4 non-zero elements:

where the individual elements only depend on the scattering angle and not on the azimuthal angle. The first element S₁₁, also known as the ``phase function'' is plotted in Fig. 2 for the wavelength of 1 m used to calculate scattered light images and polarisation maps.

A spherical grain of 1 m at a wavelength of 1 m is in the middle of a resonance region, where constructive and destructive interference within the dust grain results in phase (and polarisation) functions with strong oscillations. These effects are not observed in the case of a grain size distribution (where the oscillations corresponding to different grain sizes are averaged) or in the case of more naturally shaped particles, like aggregates. However, we choose these dust properties as they represent a better test case for the radiative transfer codes. The oscillations in the elements of the Mueller matrix must be seen in the synthetic maps allowing a more detailed comparison between codes. For comparison an Henyey-Greenstein phase function with the same asymmetry parameter is over-plotted in Fig. 2.

All the calculations presented in Sect. 4 are done assuming anisotropic scattering and use the previously presented Mueller matrix. To compare the results of the Monte Carlo codes with those obtained from discrete ordinate codes, the temperature structures and SEDs are also calculated in the isotropic case (i.e. S₁₁ is constant), using the opacities calculated with the Mie theory (Sect. 5).

3.3 Maps and SEDs

SEDs, images and polarisation maps are calculated at 10 inclinations equally spaced in cosine, i.e. for , 0.15, ..., 0.85 and 0.95. This corresponds to inclination angles ranging from 18.2 to 87.1 from pole-on. Scattered light images and polarisation maps are calculated at 1 m. The pixel scale is 25.61 mas.pixel^-1 (i.e. 251 pixels for a physical size of 900 AU at a distance of 140 pc). This is roughly a factor 2 smaller than the pixel scale of the WFPC and ACS cameras on-board the Hubble Space Telescope.

4 Results

4.1 Temperature structures

Figures 3 and 4 show the temperature distributions calculated by the different codes, as well as the difference between codes. Overall, the agreement is very good with difference almost always smaller than 10%.

Figure 3 presents the temperature along the disc midplane for the and cases. Very close to the inner edge, where the disc is directly heated by the star, the agreement between codes is excellent with maximum differences of the order of 1%. At large radii (>100 AU), the disc becomes optically thin at optical wavelengths in the vertical direction, and the midplane is heated by the stellar light that is scattered in the surface layers of the discs. In these regions, the peak-to-peak differences between codes remain below 1.5 and 3% for and cases, respectively. This shows that all codes deal smilarly with the redistribution of energy by anisotropic scattering. This is confirmed by the vertical cut at a radius of 200 AU (Fig. 4, right panel), where the differences are of the order of 1%, except at the turnover point from optically thin to optically thick (the place where the temperature suddenly drops) where differences reach 5%.

Differences in the radial temperature profile are larger between the inner edge and 1 AU, i.e. in regions were the stellar radiation does not penetrate, even via scattering. In these regions, the heating mechanism is the dust re-emission by the upper layers of the disc, which represent the most difficult case for RT codes. Nevertheless, the peak-to-peak differences between codes remain limited, smaller than 4% for the  case. For the cases, differences are most of the time smaller than 10%, except in a very small regions between  AU and  AU where the maximum difference is 20%. Differences remain also very small in the vertical direction as shown in the left panel of Fig. 4. They are below 5% from the midplane up to the disc surface, where they become smaller than 1%.

4.2 SEDs

The emerging spectral energy distributions for the and models are presented in Fig. 5 for different inclinations ranging from an almost face-on ( ) to an almost edge-on disc ( ).

The shape of the SED is strongly dependent on the inclination, moving from a stellar photosphere plus disc excess for low inclinations to a double-bumped SED, typical of very close to edge-on systems, when the stellar photosphere is obscured by the disc. The disc being optically thick at short wavelengths, the visible and near-IR stellar light is blocked and the emission in this wavelength range is dominated by the stellar scattered light coming from the disc. At longer wavelengths (>10-12 m) the dust emission dominates, resulting in a steep positive slope and a double-bumped SED. Not surprisingly, the transition between these two characteristic SEDs also depends on the disc opacity. For instance, for an inclination of 75.5, the star is still seen directly in the case, whereas it is already strongly obscured in the case.

It is interesting to note that for the most edge-on case, the flux in the optical is almost independent of the disc opacity. Similarly, in the most optically thick case, there is only small differences in the SEDs at 81.4 and 87.1. Indeed, when the star is completely obscured, the flux is dominated by stellar light that has scattered on the upper layers of the outer disc. This scattered light is mainly a function of the dust properties and of the scattering geometry but not of the optical depth in the line of sight. In these cases, optical and near-infrared photometry cannot be used to get estimates of the extinction by dust of the central object.

At low inclinations, the characteristic 9.8 m amorphous silicate feature is seen in emission. At higher inclinations, for example for the case, the feature is now observed in absorption over the continuum emission of the disc. It should be noted that, at very high inclinations, although the dip is roughly centred on the silicate feature, it is not associated to it as demonstrated by the exceptional breadth of the feature. The disc is now optically thick in the silicate feature but also in the adjacent continuum, and the absorption feature from the silicates vanishes.

At long wavelengths (>500 m), the disc become optically thin in most of its parts and the emerging flux no longer depends on the system inclination.

Overall, the differences between codes are small in all cases, which is illustrated by the thinness of the grey envelopes around the lines in Fig. 5, representing the range of the results obtained by the different codes. Figure 6 gives a more detailed view of the differences between the codes for the various inclinations and optical depths.

The left column show the results for the results for the case. When the star is seen directly (i=18.5 and ), the agreement between the three codes is excellent, with peak to peak difference smaller than 5% over the whole wavelength range. Closer to edge-on, differences remain smaller than 10%, except at very short wavelengths (<1 m) where the contribution from scattered light is dominating the SED. As the wavelength becomes shorter, the scattering becomes more forward throwing and calculations are very sensitive to the angular sampling of the scattering phase function of the codes. As a result, the agreement between codes becomes worse at short wavelengths. Differences remain smaller than 15% down to m and significant differences are only observed around m.

For the case, the disc is seen at a grazing incidence and the optical depth from the star to the observer is varying strongly across the stellar disc, from 2 at the top of the stellar surface to 200 at the bottom. In this case, using a point source for the star does not provide the correct result, and special care must be taken to resolve the stellar photosphere. In the case of a uniformly radiating sphere, as presented in this benchmark, the origin of the photon packets is uniformly distributed on the stellar surface, and a uniform distribution in the cosine of the angle between the photon direction and the normal to the surface at the point of origin is used to set the initial propagation direction of packets.

For the case (right panel of Fig. 6), the results are similar to the ones for the lower optical depth case, but with larger differences in the near- and mid-infrared. This part of the SED is one of the most challenging wavelength range for RT codes, where different contributions (thermal emission from the inner disc seen directly or through the outer disc, direct or scattered stellar light, scattered thermal emission from the disc) can dominate the emerging flux depending on the system geometry and dust opacities. Furthermore, most of the flux in this wavelength regime is coming from the inner edge of the disc and the output spectrum is very sensitive to the grid resolution adopted by the different codes. MCMax and MCFOST present very good agreement over all the wavelength range, including in the near- and mid-infrared regime, and for all inclinations. TORUS shows slightly larger differences, probably due to a lower spatial resolution at the inner edge. We note that the spatial resolution of the TORUS code is limited by the maximum cell depth in the AMR grid, which is currently set to 30, correponding to a dynamical range of 2³⁰ (a limit dictated by numerical precision in the photon path integrator). These differences are maximum around m and vary from 15% in the low inclination case up to 30% when the inclination is increasing.

4.3 Scattered light images and polarisation maps

Figure 7 presents the scattered light images of the most optically thick disc for and . The synthetic maps clearly display the effects of the anisotropy of the scattering. The oscillations in the phase function (Fig. 2) are directly observed in the maps.

The three panels on the right present the flux obtained by the various codes, and the corresponding differences, along horizontal and vertical cuts in the images. The codes agree to within 10% where the flux is significant. TORUS shows slightly larger deviations that are due to a larger Monte carlo noise in the images. All codes predict the same oscillations as a function of the position, indicating that the implementation of the scattering phase function is correct in all codes. At greater radii, larger differences are observed due to the various grid geometries used by the different codes. For instance, some codes use a spherical grid whereas other codes use a cylindrical grid with a vertical cut-off in density. The sampling of the density structure at the outer edge of the model domain is then slightly different for each code. This results in systematic differences between codes, but only in the regions of the synthetic maps where the flux is extremely low.

The vertical cut (# 3) samples the disc ``dark lane'', corresponding to the optically thick midplane. In this region, the flux is dominated by photons that have scattered several times in the disc before reaching the observer. The very good agreement between codes in the dark lane indicates that all of them deal properly with multiple scattering.

Figure 8 presents the polarisation maps for the same configurations as in the Fig. 7. The maps show complex patterns due to the strong variations of the elements of the Mueller matrix with the scattering angle. All codes produce very similar maps (see horizontal and vertical cuts), at both inclination angles, with differences smaller than 5 points of polarisation degree in the central regions of the maps, i.e. in regions where the flux is large enough to allow resolved polarisation measurement in actual observations. At greater radii, differences become slightly larger, but they remain limited to 10 points of polarisation degree. TORUS shows larger deviations, due to a larger Monte Carlo noise in the simulations, which biases the polarisation degree towards larger values. In regions where the polarised flux is large, the agreement between TORUS and the other codes is also very good. For the highest inclination, the results from TORUS are not shown due to a low signal-to-noise.

5 Comparison with discrete ordinate codes

All results presented thus far were obtained using Monte Carlo codes. To further compare the predictions of radiative transfer codes, we present, in this section, results obtained with discrete ordinate codes and discuss their limitations compared to those of Monte Carlo codes. Discrete ordinate codes solve the radiative equation along predetermined sets of directions. This integration can be performed with ``long'' or ``short characteristics'' and various schemes can be used to iterate between the temperature structure and specific intensity (or its moments), like the Accelerated Lambda Iteration (ALI) or Variable Eddington Tensor (VET) methods (see Mihalas & Mihalas 1984).

The benchmark problem presented in this section is the same as in the previous sections but with the following modifications:

• the star is now considered as a point source, but with the same spectrum and luminosity as previously defined;
• the scattering is assumed either to be isotropic (but with the same opacities as before), either to be negligible (the scattering opacity is set to 0).

These tests have been defined to allow a larger number of codes to reproduce the calculations. They correspond to an extension of the P04 benchmark to optical depths higher than 100, but without any further complexity. In this paper, they have been calculated by at least one of the previously tested Monte Carlo code and by two additional discrete ordinate codes and one additional Monte Carlo code.

5.1 Code description

5.1.1 ProDiMo

ProDiMo is an acronym for Protoplanetary Disk Model (Woitke et al. 2009) which consistently solves the chemistry, the heating/cooling balance of the gas, the dust radiative transfer and the vertical stratification of protoplanetary discs, mainly for the purpose of interpreting far IR to mm gas emission lines.

For realistic gas models, it is essential to calculate the dust temperature structure in the disc as well as the transport of UV photons including scattering, which drive the photo-chemistry. Furthermore, radiative pumping by continuum radiation changes the non-LTE population of atoms and molecules and have an important impact on the cooling rates. These strong physical couplings necessitate to solve a full 2D dust radiative transfer as one module in a global iterative procedure. It is this radiative transfer module inside ProDiMo that participates in this benchmark test. Its basic task is to provide and for the gas modelling - it is not meant for the interface to dust observations (SEDs, scattered light images etc.).

ProDiMo solves the frequency-dependent 2D dust continuum radiative transfer of irradiated discs by means of a simple, ray-based, long-characteristic method. From each grid point in the disc, a limited number of rays (here 172) are traced backwards, while solving the radiative transfer equation with isotropic scattering. The setup of the ray directions is critical for the optically thin parts of the disc, in particular at near IR wavelengths where the illumination originates from small and far hot regions. This is done in a manual fashion in ProDiMo to ensure there are more rays pointing toward the hot inner regions than toward the cooler interstellar side. One central ray pointing toward the star is reserved and covers the solid angle occupied by the star.

Instead of a treatment with a large number of wavelength grid points, ProDiMo uses a coarse wavelength grid (here ) from 100 nm to 1000 m, and treat the opacities, intensities and source functions with band means, e.g. where . One radiation transfer iteration takes about 4 s for a low resolution grid, and about 70 s for a  grid on a single-processor 2.66 GHz Linux machine, which is comparable to the computational efforts taken to solve the disc chemistry.

In order to solve the condition of radiative equilibrium and the scattering problem, a simple -type iteration is applied. The source functions are pre-calculated on the grid points and fixed during one iteration. During the ray-tracing, the opacities and source functions are interpolated from the grid point values. After having solved all rays from all points in all frequency bands, the mean intensities are updated and the dust temperatures are re-calculated. Without further accelerations, this  iteration converges only for problems up to a midplane optical depth of about 100. In order to accelerate the convergence, we apply the procedure of Auer (1984) known as ``Ng''-iteration. This enables us to solve radiative transfer problems up to , depending on the geometry of the model.

For higher optical depths, we apply an approximate procedure following an idea of C.P. Dullemond which consists in reducing the dust density in the central midplane regions in the following way. For every vertical column (considering the downward direction) we do not increase the dust density any further once a certain critical optical depth at 1 m is reached ( ), provided that the radial optical depth is also > . With this ``trick'', we can manage test problems up to with this code.

5.1.2 RADMC

RADMC is a Monte-Carlo based continuum radiative transfer code for 2-D axisymmetric configurations such as circumstellar discs and envelopes. The basic algorithm is that of Bjorkman & Wood (2001), but with a continuous deposition of energy instead of the discrete deposition as described in the original paper. In this way the temperature profile is also smooth in the very optically thin regions of the model. The temperature corrections are not computed every time a photon package enters a cell and leaves some energy, but only if the energy deposited since the last temperature update is larger than some threshold value. In this way the not so cheap temperature update is done only when needed. Of course, the higher this threshold, the faster the code is but the less reliable the result. Typically a threshold of a few percent is taken, meaning that if the energy deposited in the cell has increased by more than a few percent of the energy as it was during the last temperature update, then a new temperature update is done. Once the main Monte Carlo process is over, the spectra and images from RADMC are made using a post-processing step: the ray tracing program RAYTRACE uses the dust temperature and isotropic scattering source function computed by RADMC to calculate the formal solution of the transfer equation along rays through the model. This yields images at every discrete wavelength bin. By integrating over the images one obtains a flux at each wavelength, i.e. a spectrum. Care is taken to arrange the pixels of the images such that all the flux is captured, both from the very outer regions and from the very inner regions. This is done using ``circular images'', described in detail in Dullemond & Turolla (2000). Because the spectra and images are computed as a post-processing step, as opposed to the more classical photon collection during the Monte Carlo process itself, it is hard to include full non-isotropic scattering. To keep the flexibility to view the object from every angle and at every wavelength without having to repeat the RADMC run, one has to store the entire scattering source function , where  and are the local directional coordinates. Such a 5-dimensional array is extremely large and requires of the order of a gigabyte of disc space, which is not very practical. Instead one could prescribe before calling RADMC at which angle you wish to view the object, removing the need to store the source function also as a function of , , meaning we have a 3-D array which is much easier to store. This is done in several other codes in this paper. This is not implemented in RADMC.

5.1.3 RADICAL-VET

The code RADICAL is a classical discrete ordinate method for radiative transfer, i.e. it is not based on a Monte Carlo approach and is therefore completely deterministic. It is based on methods that are routinely used in models of stellar atmospheres, with some adaptions. The algorithm is that of Variable Eddington Tensors (VET), which is a multi-dimensional version of the method of Variable Eddington Factors (VEF) described in the book by Mihalas & Mihalas (1984). A 1-D version of this method, with special application to the kind of continuum radiative transfer problems encountered in protoplanetary discs, was described in Dullemond et al. (2002). For such 1-D geometries the method is extremely accurate and efficient. It works well and converges quickly for optical depths ranging from small (1) to extremely large (10⁶). The RADICAL-VET code is a 2-D version of this algorithm. For 2-D or 3-D geometries the method has some numerical difficulties related to the computation of the flux-mean opacity in regions of extremely low flux (e.g. the midplane of a passive irradiated disc). In practice, however, these difficulties are not fatal, although they could lower the reliability of the method in such flow-flux regions. For further details on the VET method the reader is referred to Dullemond et al. (2002), even though this paper describes only the 1-D version of this method.

5.2 Results and discussion

The midplane temperature profiles are presented in Figs. 9 and 10 for the cases with isotropic scattering and no scattering respectively. Convergence was not properly reached for some codes for the highest mass case with isotropic scattering and we present the results for the case instead.

The agreement between Monte Carlo codes is again very good in the isotropic case (Fig. 9), with peak-to-peak difference smaller than a few percents in the case and smaller than 15% in the case. Differences with discrete ordinate codes are significantly largers.

This paper states an important verification test for the method of reducing the dust density implemented in ProDiMo and described in Sect. 5.1.1. It shows that the error of this procedure in comparison to the results of the latest Monte Carlo codes combined with diffusion solvers is smaller than about 20% in the midplane regions. We note however that the relative precision in general is much better than in the midplane. The average difference between ProDiMo and MCFOST temperature structures, is smaller than 5% (1 sigma-deviations) in all cases.

The benchmark test has revealed three major problems quite typical for ray-based codes. First, the limited number of fixed rays causes some small artifacts of the temperature structure in the optically thin distant midplane - these problems can be solved by using a larger number of rays which is, however, computationally expensive. Second, the precise temperature determination in the optically thick core of the disc is hard with simple discrete ordinate codes like ProDiMo. The numerical solution of the radiative transfer equation has always some discretisation errors superimposed, and in an optically thick situation it is the small difference that determines the next temperature iteration. This eventually limits the quality of the ``forecast'' by the Ng-iteration, and disables the convergence for very optically thick problems. Third, the results close to the midplane suffer from the very large gradients present close to the inner boundary, and the results depend on the numerical details, e.g. how to interpolate the source function.

Comparisons with the RADICAL code illustrates some of the difficulties of the VET method at very high optical depths. RADICAL overestimates the midplane temperature in the central regions of the disc by about 20% and 40% for the low and high mass cases respectively (Fig. 10). The 2-D geometry introduces a number of difficulties which make the 2-D VET algorithm less stable and less reliable than its 1-D version. The main problem is the choice of discrete angular coordinates at each grid point. For the formal transfer RADICAL-VET uses the method of Short Characteristics. The choice of the angular distribution of these short characteristics is essential for the reliability of the result. In Dullemond & Turolla (2000) a good choice was described, but in the end no choice is perfect and the reliability of the results may depend on this. Another difficulty is that the VET method uses flux-mean opacities which are the generalisations of Rosseland mean opacities. By using flux-mean opacities instead of Rosseland mean opacities we may expect to get the true result instead of just an approximate result. But in regions where the flux is extremely small, such as in the midplane of an extremely optically thick disc, taking the average of the opacity based on a quantity that is nearly zero is dangerous and may lead to large errors. In spite of these caveats the VET algorithm gives reasonable results for most problems.

6 Summary

We have presented solutions for the continuum radiative transfer in high-opacity circumstellar discs. The problems have optical depths up to 10⁶ and include anisotropic scattering. They represent realistic configurations for discs around low mass stars and validate the use of the codes to model current and future observations of discs.

We have compared the results of four independent Monte Carlo codes for the temperature structure, SEDs, scattered light images and/or polarisation maps, in the case of anisotropic scattering. Overall, the agreement between codes is very good. In the most optically thick case, SEDs agree within 20% over almost all of the wavelength range. Differences become larger only at wavelengths shorter than 0.2 m for edge-on configurations, i.e. when the flux is extremely low and not observed in practice. Pixel-to-pixel differences in high-resolution scattered light images remain limited to 10% and the polarisation maps do not differ by more than 5 points of polarisation degree in regions where the polarisation can be effectively measured by observations. Each observation (SED, image or polarisation map) was reproduced by at least three of the codes, providing robust solutions to test other RT codes.

The benchmark problems represent challenging test cases for RT codes. The convergence of Monte Carlo methods alone become extremely slow for the most optically thick cases and specific numerical schemes are required to efficiently compute the temperature structure, emerging SEDs and images, for instance combining a Monte Carlo approach with ray-tracing and/or diffusion approximation methods.

Comparisons between Monte Carlo codes and discrete ordinate codes, in the cases with isotropic scattering or without scattering, show relatively large differences, that increase with optical depth. Ray-tracing and VET methods were successfully compared to Monte Carlo methods up to moderate optical depths ( , P04) but the convergence of such codes seem to become delicate in the test cases presented here. They provide good approximate solutions but must be used with care at high optical depths, when the goal is to perform detailed comparisons with observations.

Acknowledgements

C. Pinte and F. Ménard would like to thank G. Duchêne for fruitful discussions on MCFOST. Some of the computations presented in this paper were performed at the Service Commun de Calcul Intensif de l'Observatoire de Grenoble (SCCI) and on the University of Exeter's SGI Altix ICE 8200 supercomputer. C. Pinte acknowledges the funding from the European Commission's Seventh Framework Program as a Marie Curie Intra-European Fellow (PIEF-GA-2008-220891).

References

• Auer, L. H. 1984, Difference equations and linearization methods for radiative transfer. (Methods in Radiative Transfer), 237 (In the text)
• Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615 [NASA ADS] [CrossRef] (In the text)
• Cashwell, E., & Everett, C. 1959, A practical manual on the Monte Carlo Method for random walk problem (New York: Pergamon Press) (In the text)
• D'Alessio, P., Calvet, N., Hartmann, L., Franco-Hernández, R., & Servín, H. 2006, ApJ, 638, 314 [NASA ADS] [CrossRef]
• Doucet, C., Habart, E., Pantin, E., et al. 2007, A&A, 470, 625 [NASA ADS] [CrossRef] [EDP Sciences]
• Dullemond, C. P., & Dominik, C. 2004, A&A, 417, 159 [NASA ADS] [CrossRef] [EDP Sciences]
• Dullemond, C. P., & Turolla, R. 2000, A&A, 360, 1187 [NASA ADS] (In the text)
• Dullemond, C. P., van Zadelhoff, G. J., & Natta, A. 2002, A&A, 389, 464 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Ercolano, B., Barlow, M. J., & Storey, P. J. 2005, MNRAS, 362, 1038 [NASA ADS] [CrossRef]
• Fitzgerald, M. P., Kalas, P. G., Duchêne, G., Pinte, C., & Graham, J. R. 2007, ApJ, 670, 536 [NASA ADS] [CrossRef]
• Fleck, Jr., J. A., & Canfield, E. H. 1984, Journal of Computational Physics, 54, 508 [NASA ADS] [CrossRef]
• Glauser, A. M., Ménard, F., Pinte, C., et al. 2008, A&A, 485, 531 [NASA ADS] [CrossRef] [EDP Sciences]
• Harries, T. J. 2000, MNRAS, 315, 722 [NASA ADS] [CrossRef]
• Harries, T. J., Monnier, J. D., Symington, N. H., & Kurosawa, R. 2004, MNRAS, 350, 565 [NASA ADS] [CrossRef]
• Ivezic, Z., Groenewegen, M. A. T., Men'shchikov, A., & Szczerba, R. 1997, MNRAS, 291, 121 [NASA ADS] (In the text)
• Juvela, M. 2005, A&A, 440, 531 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Kurosawa, R., Harries, T. J., Bate, M. R., & Symington, N. H. 2004, MNRAS, 351, 1134 [NASA ADS] [CrossRef]
• Lucy, L. B. 1999, A&A, 345, 211 [NASA ADS] (In the text)
• Mihalas, D., & Mihalas, B. W. 1984, ApJ, 283, 469 [NASA ADS] [CrossRef] (In the text)
• Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009, ArXiv e-prints
• Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797 [NASA ADS] [CrossRef] [EDP Sciences]
• Pinte, C., Fouchet, L., Ménard, F., Gonzalez, J.-F., & Duchêne, G. 2007, A&A, 469, 963 [NASA ADS] [CrossRef] [EDP Sciences]
• Pinte, C., Ménard, F., Berger, J. P., Benisty, M., & Malbet, F. 2008a, ApJ, 673, L63 [NASA ADS] [CrossRef]
• Pinte, C., Padgett, D. L., Ménard, F., et al. 2008b, A&A, 489, 633 [NASA ADS] [CrossRef] [EDP Sciences]
• Pontoppidan, K. M., Dullemond, C. P., Blake, G. A., et al. 2007, ApJ, 656, 991 [NASA ADS] [CrossRef]
• Steinacker, J., Chini, R., Nielbock, M., et al. 2006, A&A, 456, 1013 [NASA ADS] [CrossRef] [EDP Sciences]
• Tannirkulam, A., Harries, T. J., & Monnier, J. D. 2007, ApJ, 661, 374 [NASA ADS] [CrossRef] (In the text)
• Tannirkulam, A., Monnier, J. D., Harries, T. J., et al. 2008, ApJ, 689, 513 [NASA ADS] [CrossRef]
• Walker, C., Wood, K., Lada, C. J., et al. 2004, MNRAS, 351, 607 [NASA ADS] [CrossRef] (In the text)
• Watson, A. M., & Henney, W. J. 2001, Rev. Mex. Astron. Astrofis, 37, 221 [NASA ADS] (In the text)
• Watson, A. M., & Stapelfeldt, K. R. 2004, ApJ, 602, 860 [NASA ADS] [CrossRef]
• Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296 [NASA ADS] [CrossRef] (In the text)
• Woitke, P. 2006, A&A, 460, L9 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, submitted (In the text)
• Wolf, S. 2003, Comp. Phys. Commun., 150, 99 [NASA ADS] [CrossRef]
• Wolf, S., Padgett, D. L., & Stapelfeldt, K. R. 2003, ApJ, 588, 373 [NASA ADS] [CrossRef]
• Wood, K., Wolff, M. J., Bjorkman, J. E., & Whitney, B. 2002, ApJ, 564, 887 [NASA ADS] [CrossRef]
• Yusef-Zadeh, F., Morris, M., & White, R. L. 1984, ApJ, 278, 186 [NASA ADS] [CrossRef] (In the text)

Footnotes

...

Because the grains are randomly oriented, the nett dust thermal emission is not polarised.

... perfect

Only limited by the numerical precision.

... difference

Over 500 representative points in the grid.

All Figures

---

Copyright ESO 2009