EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 520 (September-October 2010) A&A, 520 (2010) A70 Full HTML
Free Access

This article has an erratum: [https://doi.org/10.1051/0004-6361/201015025e]

Issue
A&A
Volume 520, September-October 2010
Article Number A70
Number of page(s) 3
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201015025
Published online 05 October 2010
A&A 520, A70 (2010)

On the modified random walk algorithm for Monte-Carlo radiation transfer
(Research Note)

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

Abstract
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

1 Introduction

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 $\bar{\kappa}$ 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, $\bar{\kappa}$, and the dust emission coefficient $\eta_\nu$, taking into account that photons can be both scattered and absorbed and re-emitted. In M09, the suggested algorithm is to first calculate $\eta_\nu$ iteratively, and to then use $\eta_\nu$ to compute D and $\bar{\kappa}$. In this note, I show that $\eta_\nu$ does not need to be solved iteratively, but can be solved directly, and I use this solution to show that D and $\bar{\kappa}$ 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

\begin{displaymath} \eta_\nu = \kappa_\nu~B_\nu(T)+ \sigma_\nu~J_\nu \end{displaymath} (1)

where $\nu$ is the frequency of the radiation, $\kappa_\nu$ is the mass absorption coefficient, $B_\nu(T)$ is the Planck function at the temperature T of the dust, $\sigma_\nu$ is the mass scattering coefficient, and $J_\nu$ is the mean intensity of the radiation field. Assuming that the radiation is isotropic, $J_\nu = I_\nu$, where $I_\nu$ is the intensity of the radiation. In the optically thick regime, $I_\nu=S_\nu$, where $S_\nu$ is the source function. Therefore,

\begin{displaymath} \eta_\nu = \kappa_\nu~B_\nu(T) + \sigma_\nu~S_\nu. \end{displaymath} (2)

The source function is defined as the ratio of the total emissivity to the total extinction, which in this case is

\begin{displaymath}S_\nu = \frac{\eta_\nu}{\chi_\nu}, \end{displaymath} (3)

where $\chi_\nu$ is the mass extinction coefficient ( $\chi_\nu = \kappa_\nu + \sigma_\nu$). Therefore, Eq. (2) can be rewritten as

\begin{displaymath} \eta_\nu = \kappa_\nu~B_\nu(T) + \sigma_\nu~\frac{\eta_\nu}{\chi_\nu}\cdot \end{displaymath} (4)

Re-arranging this equation, one obtains

\begin{displaymath}\eta_\nu = \frac{\kappa_\nu~B_\nu(T)}{1 - \sigma_\nu/\chi_\nu}, \end{displaymath} (5)

which can be simplified, since:

\begin{displaymath}\frac{\kappa_\nu}{1-\sigma_\nu/\chi_\nu} = \frac{\kappa_\nu}{... ...\chi_\nu} = \frac{\kappa_\nu}{\kappa_\nu/\chi_\nu} = \chi_\nu. \end{displaymath} (6)

Therefore,

\begin{displaymath} \eta_\nu = \chi_\nu~B_\nu(T). \end{displaymath} (7)

The source function for this emissivity is

\begin{displaymath}S_\nu = \frac{\eta_\nu}{\chi_\nu} = B_\nu(T), \end{displaymath} (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:

\begin{displaymath} \eta_\nu = \kappa_\nu~B_\nu(T) ~\frac{\displaystyle{\int_0^\... ...eta_\nu~{\rm d}\nu}}+ \frac{\sigma_\nu}{\chi_\nu}\eta_\nu\cdot \end{displaymath} (9)

Their original equation included a ${\rm d}B_\nu(T)/{\rm d}T$ term instead of $B_\nu(T)$, because they make use of the Bjorkman & Wood (2001) temperature correction method for determining the dust temperature. I use $B_\nu(T)$ here for consistency, but throughout the derivation, $B_\nu(T)$ can be replaced by ${\rm d}B_\nu(T)/{\rm d}T$ 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 $B_\nu(T)$ - the addition of scattering should not affect the thermal emissivity if $\kappa_\nu$ is held constant. However, even with this extra term, Eq. (9) can in fact be simplified. One can set

\begin{displaymath} \epsilon = \frac{\displaystyle{\int_0^\infty\eta_\nu~\frac{\... ...~{\rm d}\nu}}{\displaystyle{\int_0^\infty\eta_\nu~{\rm d}\nu}} \end{displaymath} (10)

which is a frequency-independent constant for any given dust type. Equation (9) then becomes:

\begin{displaymath}\eta_\nu = \kappa_\nu~B_\nu(T)~\epsilon+\frac{\sigma_\nu}{\chi_\nu}~\eta_\nu \end{displaymath} (11)

M09 suggest solving this through an iterative method, but in fact, this can be solved exactly, by re-arranging for $\eta_\nu$ as for Eq. (4) in Sect. 2.1, giving:

\begin{displaymath}\eta_\nu = \chi_\nu~B_\nu(T)~\epsilon. \end{displaymath} (12)

Substituting this back into Eq. (10) gives

\begin{displaymath}\epsilon = \frac{\displaystyle{\int_0^\infty\chi_\nu~B_\nu(T)... ...yle{\int_0^\infty\chi_\nu~B_\nu(T)~\epsilon~~{\rm d}\nu}}\cdot \end{displaymath} (13)

The $\epsilon$ in the integrals cancel out (since they are frequency independent), so that

\begin{displaymath}\epsilon = \frac{\displaystyle{\int_0^\infty\kappa_\nu~B_\nu(... ...\rm d}\nu}} = \frac{\bar{\kappa}_{\rm P}}{\bar{\chi}_{\rm P}}, \end{displaymath} (14)

where $\bar{\kappa}_{\rm P}$ and $\bar{\chi}_{\rm P}$ 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

\begin{displaymath} \eta_\nu = \chi_\nu~B_\nu(T)~\frac{\bar{\kappa}_{\rm P}}{\bar{\chi}_{\rm P}}\cdot \end{displaymath} (15)

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

\begin{displaymath}D = \frac{\langle d^2\rangle}{6~\langle {d}\rangle} = \frac{1... ...ystyle{\int_0^\infty\frac{\eta_\nu}{\rho\chi_\nu}~{\rm d}\nu}} \end{displaymath} (16)

where $\langle d\rangle$ is the mean free path of the photons, $\langle d^2\rangle$ is the mean of the path lengths squared, and $\rho$ is the density of the dust. Substituting either Eq. (7) or (15) into the above gives

\begin{displaymath}D = \frac{1}{3}~ \frac{\displaystyle{\int_0^\infty\frac{\chi_... ...\infty\frac{\chi_\nu B_\nu(T)}{\rho\chi_\nu}~{\rm d}\nu}}\cdot \end{displaymath} (17)

The $\chi_\nu$ terms can be simplified, and the $\rho$ term can factored out as it is not dependent on frequency:

\begin{displaymath}D = \frac{1}{3\rho}~ \frac{\displaystyle{\int_0^\infty\frac{B... ...d}\nu}}{\displaystyle{\int_0^\infty B_\nu(T)~{\rm d}\nu}}\cdot \end{displaymath} (18)

The second term can easily be recognized as $1/\bar{\chi}_{\rm R}$, the inverse of the Rosseland mean mass extinction coefficient. Thus,

\begin{displaymath} D = \frac{1}{3\rho\bar{\chi}_{\rm R}}\cdot \end{displaymath} (19)

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

\begin{displaymath} \bar{\kappa} = \frac{\displaystyle{\int_0^\infty\frac{\eta_\... ...yle{\int_0^\infty \frac{\eta_\nu }{\chi_\nu}~{\rm d}\nu}}\cdot \end{displaymath} (20)

This is a mean frequency-dependent mass absorption coefficient weighted by the probability of emission at a given frequency $\eta_\nu$, and by the mean free path at each frequency, $1/\rho\chi_\nu$ (where the $\rho$ term cancels out). Substituting Eqs. (7) or (15) into Eq. (20) gives

\begin{displaymath}\bar{\kappa} = \frac{\displaystyle{\int_0^\infty\frac{\chi_\n... ..._0^\infty \frac{\chi_\nu~B_\nu(T)}{\chi_\nu}~{\rm d}\nu}}\cdot \end{displaymath} (21)

As before, this can be simplified by canceling the $\chi_\nu$ terms:

\begin{displaymath}\bar{\kappa} = \frac{\displaystyle{\int_0^\infty B_\nu(T)~\ka... ...le{\int_0^\infty B_\nu(T)~{\rm d}\nu}} = \bar{\kappa}_{\rm P}, \end{displaymath} (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:

\begin{displaymath}d_{\rm min} > \frac{\gamma}{\rho\bar{\chi}_{\rm R}}\cdot \end{displaymath} (23)

The $\gamma$ parameter controls the balance of speed versus accuracy. If $\gamma$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 $\gamma$ 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 $\zeta\in[0,1]$ and uses it to solve for y in the following equation:

\begin{displaymath}\zeta = 2~\sum_{n=1}^{\infty}~(-1)^{n+1}~y^{n^2}. \end{displaymath} (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 pre-compute the sum very accurately for a range of y values, and to then interpolate for y given $\zeta$ in the Monte-Carlo code. Once the value of y is determined, one can compute the distance traveled to exit the diffusion sphere, using

\begin{displaymath}ct = -\ln{y}~\left(\frac{R_0}{\pi}\right)^2~\frac{1}{D} \end{displaymath} (25)

where D is the diffusion coefficient given by $D = 1/ 3\rho\bar{\chi}_{\rm R}$. 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

\begin{displaymath}E = E_\gamma~ct~\rho~\bar{\kappa_{\rm P}} \end{displaymath} (26)

where $E_ \gamma $ is the energy of the photon, $\rho$ is the density, and $\bar{\kappa}_{\rm P}$ is the Planck mean mass absorption coefficient.

Finally, the emergent spectrum from an optically thick region is given by $I_\nu=S_\nu=B_\nu(T)$, 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 $dB_\nu(T)/{\rm d}T$ rather than simply $B_\nu(T)$.

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 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.


References

  1. Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615 [NASA ADS] [CrossRef] [Google Scholar]
  2. Lucy, L. B. 1999, A&A, 344, 282 [NASA ADS] [Google Scholar]
  3. 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 (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.