Issue 
A&A
Volume 520, SeptemberOctober 2010



Article Number  A70  
Number of page(s)  3  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201015025  
Published online  05 October 2010 
On the modified random walk algorithm for MonteCarlo radiation transfer
(Research Note)
T. P. Robitaille
Spitzer Postdoctoral Fellow, HarvardSmithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA, 02138, USA
Received 20 May 2010 / Accepted 12 September 2010
Abstract
Min et al. (2009) presented two complementary techniques that
use the diffusion approximation to allow efficient MonteCarlo
radiation transfer in very optically thick regions: a modified random
walk and a partial diffusion approximation. In this note, I show that
the calculations required for the modified random walk method can be
significantly simplified. In particular, the diffusion coefficient and
the mass absorption coefficients required for the modified random walk
are in fact the same as the standard diffusion coefficient and the
Planck mean mass absorption coefficient.
Key words: radiative transfer  diffusion  circumstellar matter  methods: numerical
1 Introduction
The problem of MonteCarlo radiation transfer in very optically thick regions  such as in the midplane of circumstellar disks  is challenging. Without any approximations, photon packets can get trapped for millions of interactions, increasing the required computational time by several orders of magnitude. Min et al. (2009, hereafter M09) presented two complementary methods to greatly improve the efficiency of MonteCarlo radiation transfer codes in very optically thick regions: a modified random walk (MRW) and a partial diffusion approximation (PDA). The MRW prevents photons from getting stuck in very optically thick regions, and the PDA allows temperatures to be calculated in regions that see few or no photons.
The essence of the MRW method is that instead of computing thousands to millions of individual absorption or scattering events for a single photon in these optically thick regions, one can make use of the solution to the diffusion approximation inside small regions to propagate the photon efficiently. MonteCarlo radiation transfer codes propagate photons in a grid made up of cells of constant density and temperature. Therefore if the mean optical depth to the edge of a cell is much larger than unity, one can set up a sphere whose radius is smaller than the distance to the closest wall, inside which the density will be constant, and travel to the edge of a sphere in a single step using the diffusion approximation.
A probability distribution function is used to sample the true distance traveled to exit this sphere (since the photon would follow a random walk inside the sphere, rather than moving in a straight line). This true distance, which depends on the radius of the sphere and the local diffusion coefficient D, can then be used along with the mass absorption coefficient to compute the total amount of energy deposited in the dust during the diffusion. This is required in order to compute the temperature in the cell accurately. Finally, the photon is emitted from a random position on the surface of the sphere with a frequency sampled from the Planck function.
M09 provide equations for D, , and the dust emission coefficient , taking into account that photons can be both scattered and absorbed and reemitted. In M09, the suggested algorithm is to first calculate iteratively, and to then use to compute D and . In this note, I show that does not need to be solved iteratively, but can be solved directly, and I use this solution to show that D and can in fact very easily be computed, resulting in both a simpler implementation of the MRW, and in some cases performance gains.
2 Derivation
2.1 Emission coefficient
In the presence of isotropic scattering, the emissivity of dust in local thermodynamic equilibrium (LTE) is given by
where is the frequency of the radiation, is the mass absorption coefficient, is the Planck function at the temperature T of the dust, is the mass scattering coefficient, and is the mean intensity of the radiation field. Assuming that the radiation is isotropic, , where is the intensity of the radiation. In the optically thick regime, , where is the source function. Therefore,
The source function is defined as the ratio of the total emissivity to the total extinction, which in this case is
(3) 
where is the mass extinction coefficient ( ). Therefore, Eq. (2) can be rewritten as
Rearranging this equation, one obtains
(5) 
which can be simplified, since:
(6) 
Therefore,
The source function for this emissivity is
(8) 
which is expected for thermal emission from dust in LTE in the optically thick regime.
2.2 Comparison with the M09 emission coefficient
In Eq. (13) of their paper, M09 wrote the emission coefficient for dust in an optically thick region as:
Their original equation included a term instead of , because they make use of the Bjorkman & Wood (2001) temperature correction method for determining the dust temperature. I use here for consistency, but throughout the derivation, can be replaced by with no impact on the final result.
Equation (9) is similar to Eq. (4),
but includes an extra term which is the ratio of two integrals. It is
not clear why this term was included by M09, because other than
possibly changing the dust temperature  which would affect
 the addition of scattering should not affect the thermal emissivity if
is held constant. However, even with this extra term, Eq. (9) can in fact be simplified.
One can set
which is a frequencyindependent constant for any given dust type. Equation (9) then becomes:
(11) 
M09 suggest solving this through an iterative method, but in fact, this can be solved exactly, by rearranging for as for Eq. (4) in Sect. 2.1, giving:
(12) 
Substituting this back into Eq. (10) gives
(13) 
The in the integrals cancel out (since they are frequency independent), so that
(14) 
where and are the Planck mean mass absorption and extinction coefficients respectively. Therefore, Eq. (9) can solved exactly rather than through an interative procedure, with the solution given by
This is very similar to Eq. (7), but includes an extra constant multiplicative term. However, the derivations in the following sections are valid regardless of whether the emissivity is calculated using Eqs. (7) or (15), as in all cases multiplicative constants to the emissivity cancel out.
2.3 Diffusion coefficient
The diffusion coefficient is given by M09 as
(16) 
where is the mean free path of the photons, is the mean of the path lengths squared, and is the density of the dust. Substituting either Eq. (7) or (15) into the above gives
(17) 
The terms can be simplified, and the term can factored out as it is not dependent on frequency:
(18) 
The second term can easily be recognized as , the inverse of the Rosseland mean mass extinction coefficient. Thus,
This means that the expression for the diffusion coefficient is in fact the same whether or not scattering is included.
2.4 Mean opacity to absorption
The mean opacity to absorption in the diffusion region is
This is a mean frequencydependent mass absorption coefficient weighted by the probability of emission at a given frequency , and by the mean free path at each frequency, (where the term cancels out). Substituting Eqs. (7) or (15) into Eq. (20) gives
(21) 
As before, this can be simplified by canceling the terms:
(22) 
which is the standard Planck mean mass absorption coefficient.
3 Implementation
In this section, I summarize the M09 MRW algorithm, in the light of the new equations derived in Sect. 2.
A good criterion for deciding whether to start the MRW procedure was
suggested by M09, and consists in determining whether the distance to
the closest cell wall is greater than a few times the Rosseland mean
free path:
(23) 
The parameter controls the balance of speed versus accuracy. If is set too low, then the sphere used for the MRW may be optically thin at long wavelengths, and therefore the photons might have escaped directly rather than diffuse if the MRW had not been used. This leads to overestimated temperatures, as described in more detail in M09. If is set too high, then the performance gains from using the MRW are decreased, since the photons still need to undergo a significant number of interactions. The starting criterion can be checked after each scattering or absorption/reemission, and does not increase the computing time noticeably.
If the starting criterion is met, a sphere of radius R_{0} and centered on the current photon position can be set up, with R_{0}
being at most the distance to the closest grid cell wall. The diffusion
approximation can then be solved exactly inside the sphere (see M09 for
details). The diffusion solution leads to the following algorithm:
first, one samples a random number
and uses it to solve for y in the following equation:
(24) 
As y tends to 1, the sum needs to be computed to higher and higher values of n
in order to preserve a constant numerical accuracy. The most efficient
way to carry out this sampling is to precompute the sum very
accurately for a range of y values, and to then interpolate for y given
in the MonteCarlo code. Once the value of y is determined, one can compute the distance traveled to exit the diffusion sphere, using
(25) 
where D is the diffusion coefficient given by . The energy absorbed by the dust grains (which is needed to calculate the temperature in both the Lucy 1999 and Bjorkman & Wood 2001 methods) during this MRW is then
(26) 
where is the energy of the photon, is the density, and is the Planck mean mass absorption coefficient.
Finally, the emergent spectrum from an optically thick region is given by , so the frequency of the photons exiting the diffusion sphere should be sampled from the Planck function at the local dust temperature. If the Bjorkman & Wood (2001) temperature correction method is used, the frequency of the photons should be sampled from rather than simply .
4 Summary
In this note, I have shown that the MRW equations presented by M09 for the emission coefficient of dust, the diffusion coefficient, and the mean opacity to absorption can be simplified considerably. The emission coefficient defined by M09 does not need to be solved iteratively, but instead is directly given by Eq. (15). The expression for the diffusion coefficient including scattering is in fact identical to that without scattering, and is given by Eq. (19). Finally, the mean opacity to absorption is simply the Planck mean mass absorption coefficient. All of these values are directly related to the Planck and Rosseland mean mass extinction and absorption coefficients, which are usually already precomputed in MonteCarlo radiation transfer codes. Thus, the number of calculations involved with the MRW can be greatly reduced, and the MRW technique can be implemented into existing codes with very little effort. An overview of the algorithm, including caveats, is given in Sect. 3. AcknowledgementsI wish to thank the referee for a careful review, and for insightful comments that improved this research note. I also wish to thank Barbara Whitney, Kenny Wood, and Katharine Johnston for useful discussions. Support for this work was provided by NASA through the Spitzer Space Telescope Fellowship Program, through a contract issued by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA.
References
 Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 1999, A&A, 344, 282 [NASA ADS] [Google Scholar]
 Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009, A&A, 497, 155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Copyright ESO 2010
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.