Radiative transfer in very optically thick circumstellar disks
M. Min^{1}  C. P. Dullemond^{2}  C. Dominik^{1}  A. de Koter^{1,3}  J. W. Hovenier^{1}
1  Astronomical institute Anton Pannekoek, University of Amsterdam,
Kruislaan 403, 1098 SJ Amsterdam, The Netherlands
2 
MaxPlanckInstitut für Astronomie Heidelberg, Königstuhl 17, 69117 Heidelberg,
Germany
3 
Astronomical institute Utrecht, University of Utrecht, PO Box 80000, 3508 TA Utrecht, The Netherlands
Received 4 December 2008 / Accepted 3 February 2009
Abstract
Aims. In this paper we present two efficient implementations of the diffusion approximation to be employed in Monte Carlo computations of radiative transfer in dusty media of massive circumstellar disks. The aim is to improve the accuracy of the computed temperature structure and to decrease the computation time. The accuracy, efficiency, and applicability of the methods in various corners of parameter space are investigated. The effects of using these methods on the vertical structure of the circumstellar disk as obtained from hydrostatic equilibrium computations are also addressed.
Methods. Two methods are presented. First, an energy diffusion approximation is used to improve the accuracy of the temperature structure in highly obscured regions of the disk, where photon counts are low. Second, a modified random walk approximation is employed to decrease the computation time. This modified random walk ensures that the photons that end up in the highdensity regions can quickly escape to the lower density regions, while the energy deposited by these photons in the disk is still computed accurately. A new radiative transfer code, MCMax, is presented in which both these diffusion approximations are implemented. These can be used simultaneously to increase both computational speed and decrease statistical noise.
Results. We conclude that the diffusion approximations allow for fast and accurate computations of the temperature structure, vertical disk structure and observables of very optically thick circumstellar disks.
Key words: radiative transfer  diffusion  stars: circumstellar matter
1 Introduction
Protoplanetary disks are the sites of planet formation and thus are the ideal environments for studying how planetary systems are formed. With the current development of new observational techniques and the increasing sensitivity of modern telescopes, the call for accurate, fast and flexible methods of interpreting these observations is increasing. Analysis tools are needed that are fast and, at the same time, that do justice to the rich detail of current observations.
An important step in translating predictions from disk properties into comparison with observations lies in solving the radiative transfer in protoplanetary disks subject to conservation of energy and thermal equilibrium. The opacity and energy balance in protoplanetary disks are dominated by their dust grains, even though they only contribute about one percent of the total disk mass. Unfortunately, determining the radiation transport is far from a trivial task, mainly because of the wide range of optical depths that have to be traced in the disk. On the one hand, the optically thin upper layers of the disk are important since they are the regions where solidstate resonances are formed that give information on the size distribution and composition of the dust grains in these layers. On the other hand, the cold midplane region of the disk is important since it contains most of the mass of the disk and thereby determines the vertical structure. Radial optical depths through this midplane can easily reach values of 10^{5}10^{6} in the visible part of the spectrum.
The required flexibility for radiative transfer can be obtained by using a Monte Carlo method (see e.g. Pinte et al. 2006; Bjorkman & Wood 2001; Niccolini et al. 2003). Its advantage is that it is accurate at all optical depths and puts virtually no restrictions on the spatial distribution of matter or on the optical properties of the dust grains. The disadvantage is that the method becomes increasingly slow for high optical depths since every interaction of a photon package has to be considered separately. Also, in the deep layers of the disk that are wellshielded from the stellar radiation, photon counts may be low, causing large errors on the derived temperature structures in these regions. Accurate values for the temperatures in the midplane regions of protoplanetary disks are crucial for computing selfconsistent vertical density distributions (Dullemond & Dominik 2004; Dullemond et al. 2007).
In this paper we discuss solutions to these problems in terms of two different formulations of the diffusion approximation for high optical depth regions. The two different formulations serve different purposes and can be used together. We present a partial diffusion approximation (PDA) which can be used to decrease statistical noise on the temperature structure and a modified random walk (MRW) procedure which is useful for decreasing computation time significantly. Both used together provide a fast and accurate radiative transfer scheme that can be used for a variety of applications. We focus on passive protoplanetary disks, i.e. where the radiation from the central star dominates, and consider only the radiative transfer through the dust component. In Sect. 2 we will outline the general radiative transfer method employed along with the implementations of the PDA and MRW. Then, in Sect. 3 we will check the accuracy of the approximations and discuss the resulting selfconsistent vertical density structures. The results are put together into recommendations for the use of the diffusion approximation in various cases in Sect. 4. Finally, in Sect. 5 we summarize our findings and draw our conclusions.
2 Computational approach
2.1 Monte Carlo radiative transfer
A flexible way of computing the transfer of radiation in an inhomogeneous, complex medium is by means of a Monte Carlo method. In this method photon packages emerging from the central star are traced through the disk, allowing them to undergo scattering, absorption and reemission events caused by the dust they encounter on their path. The Monte Carlo method for radiative transfer in dusty media has become particularly efficient since the ideas of continuous absorption and reemission have been proposed and worked out by Bjorkman & Wood (2001). The method is fully described in their paper, and we will not repeat it here.
Two major difficulties of Monte Carlo radiative transfer are:
 1.
 the statistical noise in regions shielded from radiation;
 2.
 the fact that computation time increases with roughly the square of the optical depth.
2.2 Diffusion approximation for optically thick regions
2.2.1 Decreasing photon noise: partial diffusion approximation (PDA)
In regions in the disk where the mean free path of the photons is much smaller than the length scale over which density, ,
and temperature, T, change, the energy transport is by radiative diffusion. This is described by the radiative diffusion equation (see Wehrse et al. 2000 who provide an extension of the classical one dimensional result by Rosseland 1924 to three dimensions)
where E is the local energy density, c is the speed of light, and D is the diffusion coefficient which in the Rosseland approximation is given by
where is the Rosseland mean free path, and is the Rosseland mean opacity given by
(3) 
Here is the wavelength dependent mass extinction coefficient, is the Planck function, and is the frequency of radiation. The Rosseland diffusion coefficient is an approximate measure. In case scattering cannot be neglected, the diffusion coefficient has to be adjusted. The expression for the diffusion coefficient we use is given in Sect. 2.2.3.
We can use the diffusion equation to improve the accuracy of the temperatures in regions of high optical depth. The way we implement this is that after the Monte Carlo procedure we solve the time independent version of Eq. (1), i.e. we put
in regions of low photon counts. In this case we get the diffusion approximation for the temperature in the cell
Equation (4) results in a system of linear equations which can be solved if we set appropriate boundary conditions. Since we wish to solve the diffusion equation only in a limited region of the disk (the regions around the midplane), the boundary conditions are simply set by the temperatures as determined by the Monte Carlo procedure at the edge of the region with low photon counts. We refer to this method as the Partial Diffusion Approximation (PDA). The accuracy of the results of this approximation will be discussed in detail later in this paper. Note that the PDA is only used after the full Monte Carlo run has already finished. Thus observables that are directly obtained from the escaping Monte Carlo photon packages are not influenced by it. The main use of the PDA is to obtain a reliable temperature structure for, for example, iteratively solving the vertical hydrostatic density distribution.
The PDA assumes that the probability that a photon escapes from the optically thick region without interactions is zero. In reality though, a photon emitted at a very long wavelength encounters a very low opacity, and thus the disk can become optically thin for these photons. This will cause these regions to cool more efficiently than predicted by the PDA. Thus, in general the PDA is expected to overestimate the temperature slightly. We will return to this point when we discuss the accuracy of the PDA.
2.2.2 Increasing computational speed: modified random walk (MRW)
Monte Carlo radiative transfer computations can be very slow when implemented in a straightforward manner, especially for geometrical configurations where a photon package has a large number of interactions with the dust before it leaves the system. For protoplanetary disks this occurs especially when a single photon package ``gets lost'' in the optically thick regions of the disk, i.e. the disk midplane. It can then take a very high number of interactions before the photon package finally leaves the system. Although only very few photons get to these regions, in some configurations of rather massive protoplanetary disks these photons completely dominate the computation time. We can try to overcome this problem by letting photons make multiple interaction steps in a single computation. This can be done by using the radiative diffusion equation as outlined below.
In the regions of radiative diffusion the photons travel by a random walk. This random walk is basically a solution to the radiative diffusion equation (Eq. (1)).
In the Monte Carlo procedure we wish to know the distance a photon travels before leaving a given region in the disk. For the Monte Carlo procedure, space is subdivided in regions of constant density and temperature and we can apply the method of Fleck & Canfield (1984). The solution to Eq. (1) in the case of an infinite homogeneous medium is simply
where is the fraction of the energy that has diffused to position in a time t, and .
Equation (5) is only valid for an infinite, homogeneous medium while the cells in the disk are of finite size. The way to solve this is outlined by Fleck & Canfield (1984) and consists of setting an absorbing boundary condition to the diffusion equation at r=R_{0}, with R_{0} being the radial boundary of the homogeneous region. This ensures that the time it takes for the photon to reach the boundary of the region is the time when the photon crosses the boundary for the first time. The solution to this is a Modified Random Walk (MRW) and is given by
To determine the likelihood that the photon is still inside a sphere with radius R_{0} after a time t we have to integrate from 0 to R_{0}
with
(8) 
In the Monte Carlo procedure we now draw a random number, P_{0}, between 0 and 1 and solve P(t)=P_{0} for t. This now gives us the time the photon has traveled before reaching the edge of the sphere. The energy that is absorbed by the dust grains per unit mass during this random walk to the edge of the sphere is equal to
(9) 
where is the energy of the photon package, and is the average mass absorption coefficient encountered by the photons
(10) 
where is the emission coefficient, and is the wavelength dependent mass absorption coefficient. The emission coefficient is proportional to the probability that a photon leaves an interaction event (scattering or reemission) with frequency . An iterative way to obtain is given in Sect. 2.2.3. The above equation averages the frequencydependent mass absorption coefficient weighted with the average distance traveled at that particular frequency. This procedure can then be repeated until the photon is close enough to the cell's edge where the regular Monte Carlo method can take over.
Since the MRW is only valid in optically very thick regions, an important aspect of the above method is a criterion to start the MRW procedure. A good criterion for this is that the smallest distance to the cell's edge is larger than a few times the Rosseland mean free path (Fleck & Canfield 1984), i.e.
The parameter can be adjusted for speed or accuracy. Values of close to unity ensure fast radiative transfer but with relatively large errors, while higher values of ensure more accurate results but at the cost of a slightly increased runtime. In addition to this criterion we also require the photon package to have experienced several tens of interactions in a single cell before the MRW is employed. Note that there is always a finite chance that a photon gets emitted at a wavelength where the optical depth is low and the photon can escape the cell without interactions. The criterion set in Eq. (11) has to be chosen such that the fraction of such photons is small. We study the effects of different values of in Sect. 3.2.
2.2.3 The diffusion coefficient
In this section we derive the diffusion coefficient used in the modified random walk procedure and the partial diffusion approximation. We start with the diffusion coefficient as given by Fleck & Canfield (1984) for the modified random walk
where d is the distance a photon package travels, is the mean free path, and is the mean of the pathlength squared. For monochromatic radiative transport at frequency we have that and .
Since we have absorption and reemission in addition to scattering, we have to consider a broad range of wavelengths. Therefore, we have to average the mean free path over all frequencies using a proper weighting function. The proper weighting function is the energy probability distribution leaving an interaction (absorption/reemission or scattering) event which is proportional to the emission coefficient .
This emission coefficient obeys the following integral equation
where is the wavelength dependent mass scattering coefficient. Equation (13) can be solved iteratively. When an guess of is available, an improved estimate can be obtained by substituting this guess into the right hand side of Eq. (13). By repeating this procedure several times the solution is easily obtained. A good initial guess in most cases is to ignore scattering and take .
Usually the procedure described above converges in a few iterations. When is computed
and
are given by
and
Note that when there is no scattering, i.e. , Eqs. (12)(15) reduce to the Rosseland diffusion coefficient given by Eq. (2).
The formulae above hold for the case of isotropic scattering, which we consider in this paper. However, simple corrections can be made for anisotropic scattering. As shown by Ishimaru (1978) in the diffusion domain similar results are obtained when
and the asymmetry parameter g are both replaced by
and g' given by
(16) 
The simplest form is by taking the solution for isotropic scattering, i.e. g'=0, with .
2.3 Description of the code: MCMax
The radiative transfer code MCMax is based on the Monte Carlo method outlined by Bjorkman & Wood (2001) to solve the temperature structure in the disk. The observables are obtained by analytic integration of the computed source function in order to avoid noise in the observables. The partial diffusion approximation as well as the modified random walk as discussed above are implemented as options. The code solves the radiative transfer in 3 dimensions but it is assumed that the geometry is axisymmetric. Sophisticated regridding algorithms make sure that the optical depth to the stellar and local radiation field is properly sampled so that the temperature structures of extreme density distributions can be solved accurately.
The scattering of radiation by the dust grains is treated in a complete sense. This means that the full angledependent Mueller matrix is used for each scattering event. In this way, also polarization maps of the disks can be obtained. A comparison of the results for anisotropic scattering with other codes will be presented in Pinte et al. (2009). In this study we assume isotropic scattering.
In addition to solving the radiative transfer equations subject to the constraint of radiative equilibrium, the MCMax code can also be used to solve the vertical scale height of the disk self consistently by iterating the radiative transfer computations and subsequently solving vertical hydrostatic equilibrium equations as outlined in Sect. 2.5.
2.4 Disk setup
For the density structure of the disk we follow partly the benchmark setup as defined by Pascucci et al. (2004). However, we modify the setup to more massive disks with higher optical depths. We deviate from the Pascucci et al. (2004) setup on three important points.
 1.
 The surface density of the disk, ,
as chosen by Pascucci et al. (2004) increases with the distance to the central star R as
.
This implies that the vertical optical depth through the disk increases with distance, putting the most optically thick regions far away from the star. We consider a radial dependence of the surface density which is a power law with
(which is closer to the theoretically predicted value, see Dullemond et al. 2007).
 2.
 The outer disk radius in the Pascucci et al. (2004) setup is chosen very large,
AU, while protoplanetary disks are thought to have smaller sizes. We therefore adopt a value of
AU.
 3.
 For the vertical scale height of the disk we choose a parameterization which is close to the case of vertical hydrostatic equilibrium. In Sect. 3.4 we will show the results for the iterated solution of the vertical density distribution.
Here is the dust density in the disk, is the total dust mass, R is the distance to the center of the star (in cylindrical coordinates), z is the height above the midplane, p is the power of the scaleheight powerlaw, h_{1} is the scaleheight of the disk at 1 AU, and and are the inner and outer radius of the disk. For the values of the parameters see Table 1. The radius of the inner edge of the disk is chosen to correspond roughly to the dust condensation temperature.
By integrating Eq. (18) through the midplane (z=0) from
to
we can compute the radial optical depth
This corresponds, for the parameters of the disks we have chosen, to an optical depth at 550 nm of for the model with .
We take a very simplistic description of the grain properties. The dust grains in the disk are homogeneous spherical particles with a radius of m composed of socalled ``astronomical silicate'' (Draine & Lee 1984). The scattering by the grains is assumed to be isotropic. We would like to note here that although the grains have in principle a nonisotropic phase function, the errors on the temperature structure and emerging spectra made by assuming isotropic scattering are generally small. So although both codes used for the calculations can in principle handle anisotropic scattering, we choose to use isotropic scattering for the computations in this paper.
Table 1: The parameters for the model setup (as used in Eq. (18)).
2.5 Vertical structure
At first we will discuss the results using the above parameterized version of the density structure. In Sect. 3.4 we solve the vertical structure using hydrostatic equilibrium computations and iterate on the vertical structure until it sufficiently converges. For these computations we kept the surface density fixed at
.
The vertical density profile is then set by the balance between pressure and gravity and is given by the solution of the following differential equation (see e.g. Dullemond 2002).
where G is the gravitational constant, and the pressure with the Boltzmann constant, the mean molecular weight (for a H_{2}, He mixture) and m_{u} is the proton mass. The procedure is now to first compute the temperature structure for an initial guess of the density distribution. These temperatures and their derivatives can then be used to solve Eq. (20), and this can be iterated several times until the convergence criterion is met.
In order to judge if the vertical structure is sufficiently converged, a convergence criterion is needed. For this we used the temperature structure of the disk. The error in the temperature structure can be computed using the statistics from the MonteCarlo procedure. This allows us to compute the difference between two iterations as compared to the statistical errors. When the average difference between the temperature structures of two consecutive iterations is less than we conclude that the model is converged.
3 Results
3.1 Results with partial diffusion approximation
In this section we will compare the results with and without the partial diffusion approximation (PDA) as outlined in Sect. 2.2. We will compare the results of the runs with 10^{5} photon packages with and without radiative diffusion with those obtained using 10^{8} photon packages without the partial diffusion approximation or modified random walk procedure. We will refer to the latter model as the reference model. We have used the PDA everywhere where in the 10^{5} photon model the number of photons visiting a cell is below 30 or the relative error in the temperature as derived from the photon statistics is higher than 5%. The relative difference in the computed temperature structure compared to the run with 10^{8} photon packages is shown in Fig. 1. In this figure negative errors correspond to the case where the temperature using the PDA is lower than that of the reference model, while positive errors represent higher temperatures. It is clear that the PDA improves the agreement with the high photon run considerably, although there are still regions with relatively large errors. These errors are partly caused by errors in the temperature at the boundary of the region where the PDA is used. Most of the diffusion occurs vertically. Since the cells at the boundary are used to set the boundary conditions of the PDA, an error in the temperature in these cells propagates all the way to the midplane, causing the vertical lines of large errors seen in Fig. 1.
The largest errors in the derived temperatures without the use of the PDA are (which lies outside the limits of errors shown in Fig. 1). This is in the regions of the disk where there are no photon packages in the 10^{5} photon run. By using the PDA we can reduce this error to .
Figure 1: The relative error in the temperature computed by the model using 10^{5} photon packages as compared to the reference model using 10^{8} photon packages as a function of position in the disk. The left panel shows the model without the radiative diffusion approximation, the right panel with the radiative diffusion approximation. The black contour line encloses the region where the number of photon packages entering a cell is below 30, which corresponds to an error in the derived temperature of 5%. The dust mass in these models is . 

Open with DEXTER 
The errors when using only 10^{5} photon packages and the PDA are still considerable. This is because by using such low number of photons, the region where the diffusion approximation is employed is large. When we increase the number of photon packages to 10^{6} we can decrease the size of this region significantly (see Fig. 2). The maximum error in this case is , but only in a very small region of the disk. This model still runs in a few minutes on a normal desktop computer and because of the high gain in accuracy seems to be the best tradeoff between speed and accuracy for most cases.
The diffusion approximation is only valid in regions with very high optical depths. The above errors are caused by the Monte Carlo noise and the intrinsic errors caused by the assumption entering the PDA. To study the size of the intrinsic errors we have used the partial diffusion approximation to compute the temperatures in the 10^{8} photon run in regions that have less than 30 000 or 3000 photons. These regions are comparable to the regions where the PDA was used for the 10^{5} and the 10^{6} photon runs respectively. However, now the temperatures at the boundary of the diffusion region are determined to a very high accuracy. Therefore, we can now see the intrinsic error made by the PDA. For the region corresponding to the 10^{5} photon run (limit of 30 000 photons) the largest relative error in the temperature compared to the reference model is still (see Fig. 3). For the region corresponding to the 10^{6} photon run (limit 3000 photons) the largest relative error is .
As mentioned in Sect. 2.2.1 the PDA always tends to overestimate the temperature. In principle this causes violation of energy conservation. When the region where the PDA is used is sufficiently shielded from the observer, this is no problem. However, one can think of exotic density distributions where large, moderately optically thick regions are shielded from stellar radiation but not from the observer. In these cases, the errors in the temperature distribution can be high and will result in errors in the computed spectra and images when these are obtained from the obtained temperature distribution. This will be visible because in these cases the integrated total luminosity will be higher than the stellar luminosity. Since the PDA is only used after the Monte Carlo run, observables directly obtained from escaping photons will not be influenced. For raytracing the observables, like spectra, images and visibilities, the noise in the temperature structure usually provides no problems. Thus, when raytracing observables it is better to use the temperature as obtained before using the PDA. For determining the vertical scale height it is required that the temperature distribution is smooth and that the temperature is well defined throughout the entire disk. Therefore, temperatures obtained by the PDA are used when computing the vertical scaleheight of the disk (or other dynamical or geometrical properties of the disk that depend on temperature and require a smooth temperature profile).
3.2 Results with modified random walk approximation
In this section we study the accuracy of the Modified Random Walk (MRW) procedure. We have run the radiative transfer code using 10^{8} photon packages with and without the MRW. We used and . The difference in execution time is shown in Table 2. In this table we present the average runtime per photon package for the various models with and without the modified random walk (MRW) method relative to the model without the MRW (which was 0.22 ms in our case on a normal desktop computer). For this comparison we have also added disks with extreme mass ( ). Although such disks are too massive to be stable in reality, it represents the results one would get for more extreme density distributions (for example a disk with the same mass but with a steeper power law index of the surface density). It is clear that the MRW speeds up the model for massive disks.
We compare the two cases with different values for again with the reference model using 10^{8} photon packages. The results are shown in Fig. 4. In the left panel we show the case for . The relative error in the temperature in this case varies between and . The errors are concentrated in the innermost regions. The second case, with , is shown in the right panel of Fig. 4 and has comparable errors. In this case the errors vary between and . More importantly, in this case the region where the errors occur is much smaller.
We conclude from the above that the temperatures are, on average, systematically underestimated by the MRW. Part of the reason for this is that for a small circumscribing sphere, the diffusion approximation used to derive the random walk procedure tends to underestimate the time it takes to leave this sphere. When the radius of the sphere is increased the approximation gets better. We can do this by increasing which ensures that the MRW is only activated when the circumscribing sphere is large enough. For we find that the errors are below everywhere in the disk, and that the region with the largest errors is very small and located near the midplane of the disk. A second reason for the errors on the temperature for low values of is that the direction of propagation of the photon package after leaving the random walk sphere is not well defined. When the radius of the sphere is increased, the travel direction at the edge becomes less important and thus the temperature estimate becomes more accurate.
Figure 2: Same as the right panel of Fig. 1 but now comparing a model using 10^{6} photon packages with the reference model using 10^{8} photon packages. 

Open with DEXTER 
Figure 3: The relative error in the temperature computed by the model using 10^{8} photon packages with PDA for the region with ``low'' photon counts as compared to the reference model using 10^{8} photon packages as a function of position in the disk. Everywhere in the disk where less than 30 000 photon packages ( left panel) or 3000 photon packages ( right panel) entered we used the PDA. This region is enclosed by the black contour line. The dust mass in these models is . 

Open with DEXTER 
Figure 4: The relative error in the temperature using 10^{8} photon packages and the MRW compared to the reference model using 10^{8} photon packages as a function of position in the disk. The left panel shows the case where , the right panel shows the case for . The dust mass in these models is . 

Open with DEXTER 
There is no effect of the errors on the temperature structure on the observables. This is because the regions where the MRW is activated are sufficiently shielded from the observer. However, note that in contrast to the PDA, the MRW is used during the Monte Carlo procedure and thus, in principle, also influences the observables directly obtained from escaping Monte Carlo photon packages. These errors will be discussed below.
The MRW makes the dependence of the computation time on the mass of the disk much weaker. In addition to the cases presented in this paper, we have performed computations for extreme density distributions (e.g. with ) where the decrease in computation time due to the random walk procedure is up to a few orders of magnitude.
3.3 Observables: the effect on the SEDs
To test the effects of the MRW on the predicted observables we have considered the spectral energy distribution (SED), visibility curves and images. We will discuss here only the results on the SEDs. For the case where we solve the hydrostatic equilibrium (Sect. 3.4) we will also show the resulting visibility curves. Note that the temperature computed using the PDA is not used in the raytracing procedure for computing the observables. Instead we use the temperature directly derived from the Monte Carlo run including the MRW procedure.
For the SEDs we looked at the entire spectrum from 0.1 m up to 1 cm. The resulting errors for the high mass disk seen under an inclination of 45 are plotted in Fig. 5. For the low mass disk the effects of using the random walk procedure is negligible, as it is only used in a very small region of the disk. For the high mass disk the effects of using the MRW with can be as large as 1% for a disk seen under moderate inclinations. For inclinations close to edge on (i.e. inclination higher than 80), the differences can be a few (up to 4) percent (not shown). For the error made using the MRW is much smaller. For moderate inclinations it is always lower than 0.2%. For nearly edge on disks the error can be up to 0.5% (not shown). The increasing error at the long wavelength side of the SED is caused by the fact that at these wavelengths one starts to see through the disk, so the temperatures in the midplane are more important. The reason that the errors are much smaller in the case of is that, although the errors on the temperature are of the same order, the region where this error is made is much smaller and deeper hidden in the disk. We conclude that the effect of the MRW procedure on the SED is very small.
Table 2: The relative average runtime per photon package for thevarious models.
Figure 5: The relative differences between the SEDs computed with the random walk procedure using two different values of as compared to the reference model without the random walk procedure. The results shown are for the high mass disk under an inclination of 45. 

Open with DEXTER 
3.4 Solving the vertical density distribution
The vertical structure of the disk is mostly set by the midplane temperature, and therefore it is important to study the effects of the PDA and MRW approximations on the iterated vertical density distribution. We expect the effects of the errors of the temperature to translate into somewhat smaller errors in the vertical structure since the scale height of the disk is roughly proportional to . In this section we first analyse how the obtained structure of the disk is influenced, and after this evaluate what the consequences are for simulated observations. As a convergence criterion we use the difference between two iterations in terms of standard deviations, , integrated over the entire disk where we use the statistical error in the temperatures to compute a normalized standard deviation.
We consider here three cases:
 Case 1:
 using 10^{5} photon packages, applying the modified random walk (), and partial diffusion approximation (using a minimum number of photons of 30).
 Case 2:
 using 10^{6} photon packages, applying the modified random walk (), and partial diffusion approximation (using a minimum number of photons of 30).
 Case 3:
 using 2 10^{7} photon packages, without applying the modified random walk. We do apply the partial diffusion approximation but with a minimum number of photons of only 1 (the partial diffusion approximation is only used to assure that no cells exist with a zero temperature, which would cause the vertical scale height to collapse).
For the low mass disk ( ) the results of the vertical height of the disk is shown in Fig. 6. In this figure we show the height of the surface of the disk where an optical depth of unity is reached in the direction perpendicular to the midplane at a wavelength of m. In the first case (using 10^{5} photon packages) the diffusion region used for the PDA is quite large, causing an overestimate of the midplane temperature and therefore an overestimate of the height of the disk above the diffusion region. This region lies just behind the inner rim. By using 10^{6} photon packages (Case 2) the diffusion region is sufficiently small to overcome this problem, and no significant differences are observed between Case 2 and 3. In all cases convergence was reached within 9 iterations.
Figure 6: The vertical surface height of the low mass, i.e. , disk. The ordinate is the height of the surface of the disk where an optical depth of unity is reached when looking from the top at a wavelength of m divided by the radius, R. All three cases described in the text are displayed. The solid black line represents Case 1, the dashed line Case 2, and the solid grey line Case 3. 

Open with DEXTER 
Figure 7: The differences in terms of (see text) between the temperature structures computed in the 26th iteration and all previous iterations. Clearly, the structure as in the 26th iteration is close to that already computed in the 19th and 12th iteration. This oscillatory behavior is caused by ``waves'' propagating along the surface of the disk. 

Open with DEXTER 
Figure 8: Same as Fig. 6 but for a dust mass of . The solid black line represents Case 1, the dashed line Case 2, and the solid grey line Case 3. The curves for Case 2 and 3 represent oscillatory solutions (see main text). 

Open with DEXTER 
Figure 9: The resulting SEDs as computed for the two disk masses after iteration on the vertical structure (solid lines). The corresponding SEDs for the parameterized disk structures are also shown (dotted lines). The distance to the object is 150 pc. 

Open with DEXTER 
For the higher mass disk it is much harder to reach convergence. For this disk mass, for both Case 2 and 3 periodic fluctuations are observed with a period of 7 to 8 iterations (see Fig. 7). These fluctuations are caused by instabilities in the disk structure as described by D'Alessio et al. (1999); Dullemond (2000); Watanabe & Lin (2008) and cause the convergence criterion to fail. As can be seen from Fig. 7 the solution is periodic, implying that the procedure is converged but to an oscillatory solution. The question currently remains if in reality these fluctuations are suppressed by viscous heating, turbulence or hydrodynamical effects causing the disk to respond slower than the heating and cooling timescale. For Case 1 (10^{5} photon packages) these fluctuations are not observed and convergence is reached in 7 iterations (see Fig. 8). However, the resulting vertical structure is too high in the region just behind the inner rim. This is, again, caused by the overestimate of the temperature at the midplane due to the PDA. We conclude that for both cases using 10^{6} photon packages is sufficient to obtain a reliable vertical density profile. In this case the diffusion region, where the PDA is employed, is deep enough in the disk to ensure the diffusion approximation to be valid and not cause large errors.
Figure 10: The relative differences between the SEDs computed for the three different cases defined in the text. The upper panel shows the results for the dust mass of , the lower for a dust mass of . 

Open with DEXTER 
Figure 11: The visibility curves for the low mass, i.e. , disk and a baseline of 10 m. The distance to the object is 150 pc. All three cases described in the text are shown. The solid black line represents Case 1, the dashed line Case 2, and the solid grey line Case 3. 

Open with DEXTER 
We have computed the spectra and visibility curves of the resulting density and temperature structures. The spectra of the iterated disks (case 1) together with those of the corresponding parametrized disks are shown for disks viewed at an inclination of 45 in Fig. 9. In order to judge the effects of the overestimate of the disk height just behind the inner rim we show the differences in the SEDs between the different cases we considered in Fig. 10 for a disk viewed at an inclination of 15 and 85. For the low mass case it is clear that the largest error is made in the 110 m region. The emission in this spectral region predominantly comes from those radii in the disk where the height of the disk is increased because of the use of the PDA (see Fig. 6), causing these regions to emit more efficiently when the PDA is employed. The errors at long wavelengths, as seen in the high mass case, are caused by the use of the MRW, as discussed above for the noniterated case. In all cases we studied the differences are smaller than 5%, independent of disk inclination. Therefore, we conclude that for modeling SED observations, it suffices to use the very fast, i.e. 10^{5} photon package, method with the PDA. Also, the effects of the oscillations in the vertical structure of the massive disk have only a very small impact (less than 5%) on the SEDs.
For the visibilities things can be different. In Fig. 11 we show the visibility curves in the Nband, as would be observed by e.g. the MidInfrared Interferometric Instrument (MIDI) at the Very Large Telescope Interferometer (VLTI), for the low mass disk at a baseline of 100 m (the disk is set at a distance of 150 pc). It is already clear that the overestimate in the vertical height for Case 1 results in a slightly different visibility profile. For the high mass disk we have no convergence of the vertical structure. For Case 1 we get very different results as compared to either Case 2 or 3 which give visibility curves that are very close together (not shown). In Fig. 12 we show the visibility curves for Case 2 after different numbers of iterations. It is clear that the oscillatory behavior seen in the surface height of the disk can be observed in the visibility curves.
The oscillations found in the vertical structure of the disk do not appear to be caused by numerical problems. The waves originate in the region that just emerges from the shadow of the inner rim, and is thus again directly illuminated by the central star. In our case this is around 89 AU. From there onwards they travel in until they fade out inside the shadow of the inner rim. The fact that it doesn't seem to be a numerical artifact means that if there are no processes to stabilize the vertical structure (like those discussed by D'Alessio et al. 1999; Dullemond 2000), the oscillations in the visibilities mentioned above could in principle be observed. The period with which these oscillations occur can be estimated from the time it takes for a sound wave to propagate from the midplane to the surface of the disk. Around 8 AU, where the waves originate, this timescale is 2 years. The period of the oscillations is around 7 iterations and thus in time is estimated to be roughly 7 times the vertical timescale, i.e. 14 years. We should note that the timescale for a wave to travel from the midplane to the surface of the disk is a function of the distance to the central star. For example, at 3 AU the timescale is only 5 months. It is hard to estimate from simple considerations if this will have an increasing or a decreasing effect on the strength of the waves. Computations using radiation hydrodynamics have to be performed to study these effects in more detail.
4 Recommendations for use and implementation
In this section we summarize the results discussed above and make recommendations for use and implementation of the PDA and the MRW procedures for various cases.
Figure 12: The visibility curves for the high mass, i.e. , disk and a baseline of 100 m. The distance to the object is 150 pc. The visibilities computed after three different numbers of iterations for Case 2 are shown. The solid black line is after 19 iterations, the dashed line after 23 iterations, and the solid grey line after 26 iterations. 

Open with DEXTER 
4.1 Passively heated disks
For radiative transfer through a fixed density structure the recommended setup for computations depends on the observables one wants to obtain. For computing SEDs or optical/infrared images it suffices to use only the MRW procedure, the PDA is then not required. Since the errors made by the MRW procedure are localized strongly in the high density regions, one can use a low value of (e.g. 15). Errors are then smaller than 1.5%. When also millimeter images or visibilities are required a higher value of is recommended (e.g. 10).
When accurate temperatures are required, e.g. for computations of chemical reactions or of the vertical structure, the PDA can be switched on. The number of photon packages used and the limit when the PDA is used can be tuned to the desired accuracy and computational speed. One needs to keep in mind that the diffusion approximation itself is not without error, so the use of a sufficient number of photon packages is advised. In order to keep the statistical error in the temperature structure below 5% a lower limit of the number of photon packages per cell has to be set to approximately 30.
There are roughly two types of situations which require different considerations for the limits of the diffusion approximation and the number of photons needed which we will refer to as local accuracy and global accuracy.
4.1.1 Local accuracy
The first, and most difficult case is when the temperature is needed at each location in the disk with high accuracy, e.g. 5%. In this case, the PDA has to be restricted to a very small region deep inside the disk to decrease the intrinsic error of the diffusion approximation. In order to do so, a sufficient number of photon packages (e.g. 10^{7}) have to be used, and the limiting number of packages has to be set to 30 to avoid statistical errors.
4.1.2 Global accuracy
The second situation is when it is most important that the temperature structure is not too noisy, i.e. the fluctuations from cell to cell need to be small. This can be the case when one wants to look at more global disk characteristics in which the temperature plays a role like, for example, when solving for the vertical structure of the disk. In this case, a larger error in the temperature can be allowed in favor of a more stable and smooth temperature structure. Consequently, fewer photon packages can be used (e.g. 10^{6}) and the limit for using the PDA can be set to 30. Although the errors on the temperature structure can in this case be about 10%, the noise level is much lower due to the smoothing caused by the PDA. Also, in this case the for the MRW procedure can be set to a relatively low value (e.g. ) since the temperatures in the region where the MRW is actually employed will be largely overwritten by the PDA.
4.2 Viscously heated disks
Although in this paper we have not considered the case of viscous heating, the MRW procedure does allow for complete transport of the energy released by viscous heating through the disk. When the MRW is not employed it is not feasible to transport a significant number of photon packages from the midplane, where the energy is released, to the surface layers. The computations we performed on the high mass disk showed that without the MRW procedure the computation time for a photon package to reach the surface layers was roughly a minute on a normal desktop computer. By using the MRW procedure with this was reduced, in our case, to less than a second. This still does not allow for a very large number of photon packages, but at least several thousands can be used. The energy released by viscous heating ensures that the regions shielded from stellar radiation also receive sufficient photon packages to determine an accurate temperature. Thus in this case the PDA can be turned off. More details on this subject are beyond the scope of this paper and will be left for a future study.
5 Conclusions
We have presented two efficient implementations of the diffusion approximation in the high density regions of massive protoplanetary disks: the partial diffusion approximation and the modified random walk procedure.
First, the partial diffusion approximation (PDA) can be used to increase the accuracy of the temperature structure in the highly obscured regions. The PDA is mainly useful when the vertical scaleheight of the disk has to be computed selfconsistently. In order to converge the vertical scaleheight computations it is important to have a temperature structure that is stable throughout the disk. The partial diffusion approximation is therefore ideal to reduce the photon noise, inherent to MonteCarlo radiative transfer computations, in the midplane regions.
Second, the modified random walk (MRW) procedure increases the computational speed significantly, while errors are small. When the modified random walk procedure is not employed the computation time required is a strong function of the mass of the disk. When it is employed, this dependence is much weaker allowing for computations of disks with much more extreme optical depths, like for example in disks with steeper density gradients.
The combination of the above two methods allows one to perform computations for massive protoplanetary disks with high accuracy and within reasonable time. The modified random walk procedure ensures that enough photons can be emitted in order to get reasonable photon counts in most of the disk. The partial diffusion approximation can then be employed to correct the temperatures in those regions where errors are still large due to poor statistics.
We have performed computations to iteratively compute the vertical scaleheight of the disk. For this we have studied carefully the effects of both approximations mentioned above. The effects of using the MRW (with ) on the vertical structure are very small, and have no impact on the observables. When very few photon packages are used, the region where one has to apply the PDA is large. This causes overestimates of the temperatures in a relatively large region of the disk, resulting in a local overestimate of the vertical height. Although this has no effect on the SED, the visibilities obtained with interferometric measurements are affected. We conclude that for computations of spatially unresolved observations, only a small number of photon packages is enough. When spatially resolved observations are to be simulated, one needs to resort to higher numbers of photon packages in order to shrink the region where the PDA has to be applied.
We also find that the vertical structure of massive disks converges to an oscillatory solution. ``Waves'' traveling the surface of the disk cause the density structure to fluctuate. These waves originate at the location where the disk emerges from the shadow of the inner rim and from there travel inwards. When the damping mechanisms for these waves are not efficient, this can be important when considering spatially resolved computations or observations of flaring protoplanetary disks.
Acknowledgements
M. Min acknowledges financial support from the Netherlands Organisation for Scientific Research (NWO) through a Veni grant. We would like to thank an anonymous referee for positive and constructive feedback on the paper.
Appendix A: Benchmark test of MCMax
The new code MCMax has been tested for accuracy. We did this by comparing the computed spectra and temperature structures to those obtained by the independently developed radiative transfer code RADMC (first used in Dullemond & Dominik 2004). Both codes have been checked to reproduce the benchmark results of Pascucci et al. (2004). In addition to this both codes are involved in benchmark computations for higher mass disks by Pinte et al. (2009) and no discrepancies with other codes have been found so far.
We tested the codes MCMax and RADMC first for the benchmark setup of Pascucci et al. (2004). All results where reproduced without any problems. The differences with the benchmark SEDs are comparable with those presented by Pascucci et al. (2004) between the different codes used in the benchmark computations. For the temperature structures we compared the results with those computed using the RADICAL code and found that the differences are in all cases less than .
Figure A.1: The relative error in the temperature computed by MCMax and RADMC. The left panel shows the model with a dust mass of , while the right panel shows the model with of dust. 

Open with DEXTER 
Figure A.2: The resulting SEDs as computed for the two disk masses. The SEDs are shown for two disk inclination angles. The solid line gives the resulting spectrum from MCMax while the dotted line shows the result as computed by RADMC. The distance to the object is 150 pc. 

Open with DEXTER 
As an additional test case we performed computations for the 10^{3} and 10^{6} solar mass dust disks given the density setup as described in Sect. 2.4 using 10^{8} photon packages for both codes. By using this number of photon packages we ensure that the photon count is very high everywhere in the disk. For this run we did not use the random walk or the diffusion approximation to ensure that the results are not influenced by the assumptions made in these approximations.
The comparison of the computed temperature structures is shown in Fig. A.1. As the error we plot
(A.1) 
For the model with of dust the differences are below 6% everywhere in the disk. All differences seem to be caused by interpolating between the different spatial grids used in both codes. For the model with of dust, the differences are slightly larger. In this case we can distinguish three types of differences:
 1.
 Those caused by the interpolation of the spatial grids.
 2.
 Those caused by the fact that the inner boundary in the RADMC code is smoothed, while in the MCMax code it is sharp.
 3.
 Errors caused by photon noise and numerical implementation of the methods.
The resulting spectral energy distributions for a disk at a distance of 150 pc at two different inclinations are shown in Fig. A.2. The differences between the RADMC and MCMax models are small, although for the high inclination cases differences can be observed even on a loglog scale. Since the shape of the spectra is identical, the most likely cause for this is the spatial grid chosen by both codes.
We conclude that the accuracy of the codes is very good, and the results compare very well.
References
 Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615 [NASA ADS] [CrossRef]
 D'Alessio, P., Cantó, J., Hartmann, L., Calvet, N., & Lizano, S. 1999, ApJ, 511, 896 [NASA ADS] [CrossRef]
 Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [NASA ADS] [CrossRef] (In the text)
 Dullemond, C. P. 2000, A&A, 361, L17 [NASA ADS]
 Dullemond, C. P. 2002, A&A, 395, 853 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Dullemond, C. P., & Dominik, C. 2004, A&A, 417, 159 [NASA ADS] [CrossRef] [EDP Sciences]
 Dullemond, C. P., Hollenbach, D., Kamp, I., & D'Alessio, P. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 555
 Fleck, Jr., J. A., & Canfield, E. H. 1984, J. Comput. Phys., 54, 508 [NASA ADS] [CrossRef] (In the text)
 Ishimaru, A. 1978, Wave propagation and scattering in random media, Vol. I, Single scattering and transport theory, Research supported by the US Air Force, NSF, and NIH (New York: Academic Press, Inc.), 267 (In the text)
 Niccolini, G., Woitke, P., & Lopez, B. 2003, A&A, 399, 703 [NASA ADS] [CrossRef] [EDP Sciences]
 Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, in press (In the text)
 Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797 [NASA ADS] [CrossRef] [EDP Sciences]
 Rosseland, S. 1924, MNRAS, 84, 525 [NASA ADS] (In the text)
 Watanabe, S.I., & Lin, D. N. C. 2008, ApJ, 672, 1183 [NASA ADS] [CrossRef]
 Wehrse, R., Baschek, B., & von Waldenfels, W. 2000, A&A, 359, 780 [NASA ADS] (In the text)
All Tables
Table 1: The parameters for the model setup (as used in Eq. (18)).
Table 2: The relative average runtime per photon package for thevarious models.
All Figures
Figure 1: The relative error in the temperature computed by the model using 10^{5} photon packages as compared to the reference model using 10^{8} photon packages as a function of position in the disk. The left panel shows the model without the radiative diffusion approximation, the right panel with the radiative diffusion approximation. The black contour line encloses the region where the number of photon packages entering a cell is below 30, which corresponds to an error in the derived temperature of 5%. The dust mass in these models is . 

Open with DEXTER  
In the text 
Figure 2: Same as the right panel of Fig. 1 but now comparing a model using 10^{6} photon packages with the reference model using 10^{8} photon packages. 

Open with DEXTER  
In the text 
Figure 3: The relative error in the temperature computed by the model using 10^{8} photon packages with PDA for the region with ``low'' photon counts as compared to the reference model using 10^{8} photon packages as a function of position in the disk. Everywhere in the disk where less than 30 000 photon packages ( left panel) or 3000 photon packages ( right panel) entered we used the PDA. This region is enclosed by the black contour line. The dust mass in these models is . 

Open with DEXTER  
In the text 
Figure 4: The relative error in the temperature using 10^{8} photon packages and the MRW compared to the reference model using 10^{8} photon packages as a function of position in the disk. The left panel shows the case where , the right panel shows the case for . The dust mass in these models is . 

Open with DEXTER  
In the text 
Figure 5: The relative differences between the SEDs computed with the random walk procedure using two different values of as compared to the reference model without the random walk procedure. The results shown are for the high mass disk under an inclination of 45. 

Open with DEXTER  
In the text 
Figure 6: The vertical surface height of the low mass, i.e. , disk. The ordinate is the height of the surface of the disk where an optical depth of unity is reached when looking from the top at a wavelength of m divided by the radius, R. All three cases described in the text are displayed. The solid black line represents Case 1, the dashed line Case 2, and the solid grey line Case 3. 

Open with DEXTER  
In the text 
Figure 7: The differences in terms of (see text) between the temperature structures computed in the 26th iteration and all previous iterations. Clearly, the structure as in the 26th iteration is close to that already computed in the 19th and 12th iteration. This oscillatory behavior is caused by ``waves'' propagating along the surface of the disk. 

Open with DEXTER  
In the text 
Figure 8: Same as Fig. 6 but for a dust mass of . The solid black line represents Case 1, the dashed line Case 2, and the solid grey line Case 3. The curves for Case 2 and 3 represent oscillatory solutions (see main text). 

Open with DEXTER  
In the text 
Figure 9: The resulting SEDs as computed for the two disk masses after iteration on the vertical structure (solid lines). The corresponding SEDs for the parameterized disk structures are also shown (dotted lines). The distance to the object is 150 pc. 

Open with DEXTER  
In the text 
Figure 10: The relative differences between the SEDs computed for the three different cases defined in the text. The upper panel shows the results for the dust mass of , the lower for a dust mass of . 

Open with DEXTER  
In the text 
Figure 11: The visibility curves for the low mass, i.e. , disk and a baseline of 10 m. The distance to the object is 150 pc. All three cases described in the text are shown. The solid black line represents Case 1, the dashed line Case 2, and the solid grey line Case 3. 

Open with DEXTER  
In the text 
Figure 12: The visibility curves for the high mass, i.e. , disk and a baseline of 100 m. The distance to the object is 150 pc. The visibilities computed after three different numbers of iterations for Case 2 are shown. The solid black line is after 19 iterations, the dashed line after 23 iterations, and the solid grey line after 26 iterations. 

Open with DEXTER  
In the text 
Figure A.1: The relative error in the temperature computed by MCMax and RADMC. The left panel shows the model with a dust mass of , while the right panel shows the model with of dust. 

Open with DEXTER  
In the text 
Figure A.2: The resulting SEDs as computed for the two disk masses. The SEDs are shown for two disk inclination angles. The solid line gives the resulting spectrum from MCMax while the dotted line shows the result as computed by RADMC. The distance to the object is 150 pc. 

Open with DEXTER  
In the text 
Copyright ESO 2009