This article has an erratum: [https://doi.org/10.1051/0004-6361/201015025e]
Volume 520, September-October 2010
|Number of page(s)||3|
|Section||Numerical methods and codes|
|Published online||05 October 2010|
On the modified random walk algorithm for Monte-Carlo radiation transfer
T. P. Robitaille
Spitzer Postdoctoral Fellow, Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA, 02138, USA
Received 20 May 2010 / Accepted 12 September 2010
Min et al. (2009) presented two complementary techniques that use the diffusion approximation to allow efficient Monte-Carlo 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
The problem of Monte-Carlo 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 Monte-Carlo 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. Monte-Carlo 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 re-emitted. 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.
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
where is the mass extinction coefficient ( ). Therefore, Eq. (2) can be rewritten as
Re-arranging this equation, one obtains
which can be simplified, since:
The source function for this emissivity is
which is expected for thermal emission from dust in LTE in the optically thick regime.
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 frequency-independent constant for any given dust type. Equation (9) then becomes:
M09 suggest solving this through an iterative method, but in fact, this can be solved exactly, by re-arranging for as for Eq. (4) in Sect. 2.1, giving:
Substituting this back into Eq. (10) gives
The in the integrals cancel out (since they are frequency independent), so that
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.
The diffusion coefficient is given by M09 as
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
The terms can be simplified, and the term can factored out as it is not dependent on frequency:
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.
The mean opacity to absorption in the diffusion region is
This is a mean frequency-dependent 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
As before, this can be simplified by canceling the terms:
which is the standard Planck mean mass absorption coefficient.
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
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 over-estimated 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/re-emission, and does not increase the computing time noticeably.
If the starting criterion is met, a sphere of radius R0 and centered on the current photon position can be set up, with R0
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:
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 pre-compute the sum very
accurately for a range of y values, and to then interpolate for y given
in the Monte-Carlo code. Once the value of y is determined, one can compute the distance traveled to exit the diffusion sphere, using
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
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 .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 pre-computed in Monte-Carlo 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. Acknowledgements
I 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.
Copyright ESO 2010
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.