Monte Carlo Radiative Transfer - Basic Principles¶

Radiative transfer describes how electromagnetic radiation (light) propagates through a medium. Since there are a large number of light-matter interactions that have the possibility of affecting the light propagation, solving the radiative transfer problem is in most astrophysical cases only possible by means of numerical calculations. While there are many different numerical techniques available to tackle this, Monte Carlo techniques have become a successful and elegant tool particularly for radiative transfer problems in supernovae.

Monte Carlo Radiative Transfer methods track a sufficiently large number of photons (light particles) as they propagate through the supernova ejecta. The initial properties of these photons are randomly (in a probabilistic sense) assigned in accordance with the macroscopic properties of the radiation field (see Energy Packet Initialization) and in a similar manner the decisions about when, where and how the photons interact with the surrounding material are made (see Propagation). Given a large enough sample, these photons behave as a microcosm of all of the light traveling through the ejecta – that is, based on the behavior of these photons, we can draw conclusions about the propagation of light through the ejecta as a whole (see Volume-Based Estimators). This is eventually used to determine the actual steady-state plasma properties (see Updating Plasma and Convergence) and the emitted spectrum (see Spectrum Generation).

Random Sampling Basics¶

Both the initialization of these photons and their propagation through the ejecta are modeled through probabilistic processes. This involves assigning probabilities to the occurrence of certain events or properties. For example, during isotropic scattering, finding a photon scattering into any given direction is equally likely. During the Monte Carlo simulation, assignments according to these probabilities have to be frequently performed. For this purpose, so-called Random Number Generators are available. These produce (pseudo-) random numbers \(z\) uniformly distributed on the interval \([0,1]\). The challenge now lies in using these numbers to sample any physical process involved in the Radiative transfer calculation. From a probability theory point of view, this just implies finding a mapping between the probability distribution governing the physical process and the one underlying the Random Number Generator. This process is typically referred to as random sampling.

Inverse transformation method¶

Note

This is a very superficial and sloppy description of Random sampling. More detailed and rigorous accounts are found in the standard literature, for example in [KalosWhitlock08].

The simplest and most-used technique in Monte Carlo radiative transfer applications is referred to as the inverse transformation method and involves the cumulative distribution function. In general, a random process captured by the random variable \(X\) is governed by a probability density \(\rho_X(x)\) (the continuous counterpart to discrete probabilities), with \(\rho_X(x) \mathrm{d}x\) describing the probability of the variable taking values in the interval \([x, x+\mathrm{d}x]\). The cumulative distribution function in turn describes, as the name suggests, the probability of the variable taking any value between \(-\infty\) and \(x\):

\[f_X(x) = \int_{-\infty}^x \mathrm{d}x \rho_X(x)\]

Since the probability density is by definition always positive, the cumulative distribution function is monotonically increasing. This constitutes the basis for the inverse transformation function. Consider two random variables, \(X\) and \(Y\). A mapping between those may be established by equating their cumulative distribution functions. Numbers \(y\) distributed according to one of the probability densities (in this case \(\rho_Y(y)\)) may then be used to sample the other process by

\[x = f_X^{-1}\left(f_Y(y)\right).\]

For the Random Number Generators described above, the cumulative distribution function is trivial, namely \(f_Z(z) = z\). However, the inverse distribution sampling method relies on finding the analytic inverse function of the cumulative distribution function governing the physical processes to be sampled. If this is not possible, other sampling methods, such as von-Neumann rejection sampling techniques, have to be used.

Examples¶

A few examples are provided to illustrate the random sampling process using the inverse transformation technique.

Isotropic Scattering¶

Consider the case of an isotropic scattering. Here, all scattering angles are equally likely. Thus, the corresponding (normalized) probability density and the cumulative distribution function are given by

\[\begin{split}\rho_{\mu}(\mu) &= \frac{1}{2}\\ f_{\mu}(\mu) &= \frac{1}{2} (\mu - 1).\end{split}\]

This leads to the sampling rule

\[\mu = 2 z - 1.\]

Next Interaction Location¶

The probability of a photon interacting after covering an optical depth \(\tau\) is given by (see Packet Propagation)

\[\begin{split}\rho_{\tau}(\tau) &= \exp(-\tau)\\ f_{\tau}(\tau) &= 1 - \exp(-\tau).\end{split}\]

With the inverse transformation method, the optical depth to the next interaction location may then be sampled by

\[\tau = - \mathrm{ln}(1 - z)\]

which is equivalent to

\[\tau = - \mathrm{ln}z.\]