Linking the Monte Carlo radiative transfer algorithm to the radiative transfer equation

2.1 Neumann Series solution to the RTE

The RTE can be expressed in terms of dimensionless variables and formally integrated along the ray trajectory, expressed as [22]

where 𝜏 is a dimensionless parametric distance: the optical thickness along the ray trajectory. Further, S ⁢ ( x → , s ^ ) is referred to as the source function, defined as

where 𝜔 is the particle albedo, defined as σ / β . Given the recursivity between equations (2.1) and (2.2), an asymptotic expansion concerning the dimensionless parameter relevant to scattering, the particle albedo (𝜔), leads to

where 𝑛 stands for the order of the perturbation expansion and is also related to the number of scattering events that a photon experiences [26]. The zero-th order moment is simply Beer’s law. Each subsequent moment ( I n + 1 ⁢ ( x → , s ^ ) ) in the perturbation expansion can be computed from its lower-order moment ( I n ⁢ ( x → , s ^ ) ) by

where x → ′ is the location from the incoming radiation, defined as x → ′ = x → − τ ′ ⁢ s ^ , and τ ′ is the dimensionless photon flight length. If 𝜔 is taken outside of the integration, equation (2.4) can be written as I n + 1 ⁢ ( x → , s ^ ) = ω ⁢ L ⁢ { I n ⁢ ( x → , s ^ ) } , where 𝐿 is a linear operator, accounting for the integrations on τ ′ and Ω ′ , leaving equation (2.3) in the following form:

where 𝑛 in L n represents the number of subsequent applications of the linear integral operator 𝐿 in equation (2.4).

2.2 Analytical solution of the lower-order terms

When considering a unidirectional point source (i.e., a delta function in phase space), the first term on the series ( I 0 ⁢ ( x → , s ^ ) ) becomes Beer’s law as it is obtained by solving the unperturbed RTE ( ω = 0 ). In an unbounded domain under a cylindrical coordinate system ( r → , z ), the unidirectional point source ( δ s ⁢ ( s ^ − k ^ ) ⁢ δ c ⁢ ( r → ) ⁢ δ ⁢ ( z ) ) located at the origin and directed along 𝑘 (positive 𝑧 direction), with a radiation intensity of I s , yields the following expression for I 0 ⁢ ( x → , s ^ ) :

where U ⁢ ( z ) is the Heaviside unit step function and r → is a vector in a plane perpendicular to k ^ . It should be noted that δ c ⁢ ( r → ) in equation (2.5) is a delta function of r → in polar coordinates. Additionally, δ s ⁢ ( s ^ − k ^ ) is a direction-wise delta function (on a unit sphere) defined as δ ⁢ ( μ − 1 ) ⁢ δ ⁢ ( ϕ ) , where 𝜇 is the cosine of the polar angle measured from k ^ , and 𝜙 the azimuthal angle.

The recursive formula in equation (2.4) can be used to compute the first-order scattering moment ( I 1 ⁢ ( x → , s ^ ) ) directly from equation (2.5), yielding the following integral:

As reported previously [26], the integration on equation (2.6) yields

Equations (2.5) and (2.7) will be crucial to derive the MCRT algorithm from the RTE. On the other hand, no other analytical or closed-form solutions are available for higher-order terms; however, as will be shown in upcoming sections, this is not required. Instead, a stochastic version of equation (2.4) will prove highly useful.

3.1 Sampling the open-form Green’s function

The MCRT algorithm aims at sampling f ⁢ ( x → , s ^ ) , describing the probability of photons or radiative energy propagating and absorbing at a given position ( x → ) and along a given direction ( s ^ ) after interacting, possibly multiple times, with suspended particles. As f ⁢ ( x → , s ^ ) is a summation of other PDFs ( f n ⁢ ( x → , s ^ ) ), the probability mixing method allows sampling f ⁢ ( x → , s ^ ) by sampling the terms or moments in equation (3.2). The weight preceding each term in the series defines the selection probability. In this sense, expanding the terms in equation (3.2) and nesting the higher-order terms directly within the operator (considering 𝐿 is a linear operator) yields

Nesting the terms of this open-form solution allows noticing that it can be split into two contributions (that add up to unity): one term corresponds to immediate photon absorption (preceded by ( 1 − ω ) ) and another due to absorption after a subsequent scattering event (preceded by 𝜔). It should be noted that the argument within the operator 𝐿 in the second term contains, again, two terms preceded by the same weights. This recurrent structure already points to two crucial facts about MCRT:

• 𝜔 is the criteria that allows deciding whether to sample a higher moment in the Neumann series solution, or not (i.e., if R < ω );

• the selection of higher-order moments can be done by repeatedly applying this same criteria.

3.1.1 Sampling lower-order terms

Besides defining a criteria to select which moment of f ⁢ ( x → , s ^ ) to sample, the algorithm must also consider how to sample each term. Sampling (refer to equation (3.2)) implies each term must be integrated in phase space as

P = ∫ − 1 μ ∫ 0 ϕ ∫ 0 r ∫ − ∞ z ∫ 0 2 ⁢ π f ⁢ ( x → , s ^ ) ⁢ r ⁢ d ϕ p ⁢ d z ⁢ d r ⁢ d ϕ ⁢ d μ .

The first term corresponds to the radiative energy transported by photons absorbed without experiencing a single scattering event. For a unidirectional point source (i.e., a photon bundle), the corresponding PDF is

f 0 ⁢ ( x → ) = e − z ⁢ U ⁢ ( z ) ⁢ δ ⁢ ( r → ) r ⁢ δ ⁢ ( ϕ ) ⁢ δ ⁢ ( μ )

Sampling this term implies determining the propagation direction ( s → ) and displacement. More specifically, the spherical (direction-wise) and radial (position-wise) delta functions, in equation (2.5), point to the fact that s → and 𝑟 (the propagation direction and the radial position, respectively), remain unchanged. Thus, computing the photon propagation along 𝑧 follows from

(3.4) P 0 ⁢ ( x → ) = ∫ − ∞ z U ⁢ ( z ) ⁢ e − z ⁢ d z .

On the other hand, the dimensionless axial position (𝑧) can be interpreted as the photon flight length times the extinction coefficient ( l 0 ⁢ β ), sampled between 0 and ∞ (due to the presence of the step function at 0, U ⁢ ( z ) ). Thus, solving equation (3.4) and sampling the probability space (𝑃), then leaves

(3.5) l 0 ⁢ ( x → ) = − 1 β ⁢ ln ⁡ ( R ) ,

where 𝑅 is a uniformly distributed random number, between 0 and 1.

The second term in equation (3.2), f 1 , corresponds to the radiative energy transported by photons absorbed after experiencing a single scattering event. The PDF for this second term reads

f 1 ⁢ ( z , r , ϕ p , Ω ) = 1 r ⁢ p ⁢ ( μ ) 1 − μ 2 ⁢ e − z ⁢ e − 1 − μ 1 − μ 2 ⁢ r ⁢ δ ⁢ ( ϕ − ϕ p ) ⁢ U ⁢ ( z r 2 + z 2 − μ ) .

The first thing to notice is that the PDF is symmetric both in space ( ϕ p ) and direction (𝜙) such that one of these variables should be sampled from a uniform distribution (e.g., ϕ = 2 ⁢ π ⁢ R ); however, the term δ ⁢ ( ϕ − ϕ p ) also implies that ϕ = ϕ p . Thus, one random number fixes both variables. This leaves a dependence between 𝜇, 𝑟, 𝑧 as follows:

(3.6) P 1 = ∫ − 1 μ ∫ − ∞ z ∫ 0 r 1 r ⁢ p ⁢ ( μ ) 1 − μ 2 ⁢ e μ ⁢ r 1 − μ 2 − z ⁢ e − r 1 − μ 2 ⁢ r ⁢ d r ⁢ d z ⁢ d μ .

Before proceeding with the selection of 𝜇, 𝑟 and 𝑧, a change of variables is required. Figure 1 highlights the change of spatial variables required to recast equation (3.6) in terms of the photon flight lengths ( l i ), more commonly used by the traditional MCRT algorithm.

Figure 1

Sample trajectory of the first scattering event, highlighting the variables involved in equation (3.6).

Thus, the arguments in the exponential can be recast, by noticing (Figure 1) that

l 0 = z − μ ⁢ r 1 − μ 2 (initial photon flight length) , l 1 = r 1 − μ 2 (subsequent photon flight length) .

Then the integrand in equation (3.6) becomes

(3.7) P 1 = ∫ − 1 μ ∫ 0 r ∫ − ∞ z p ⁢ ( μ ) ⁢ e − l 0 ⁢ e − l 1 ⁢ d l 0 ⁢ d l 1 ⁢ d μ ,

as the denominator 1 r cancels out and 1 / 1 − μ 2 is absorbed during the change of variables (as the determinant of the jacobian relating l 0 and l 1 to 𝑟 and 𝑧). Moreover, noticing that, under the change of variables, the step function U ⁢ ( z r 2 + z 2 − μ ) becomes U ⁢ ( l 0 > 0 ) , the integration limits are analogous to equation (3.4) so that equation (3.7) can be expressed as

(3.8) P 1 = ∫ 0 l 0 e − l 0 ⁢ d l 0 ⋅ ∫ − 1 μ p ⁢ ( μ ) ⁢ d μ ⋅ ∫ 0 l 1 e − l 1 ⁢ d l 1 ,

yielding separable factors that point to the fact that this multi-dimensional PDF can be evaluated as three uncorrelated stochastic processes, i.e., as a Markov chain of events. Thus, sampling l 0 and l 1 is done as in equation (3.5), while 𝜇 is sampled by integrating and inverting the applicable phase function ( p ⁢ ( μ ) ) as μ = P ⁢ ( R ) − 1 .

3.2 MCRT algorithm

Thus, as presented in Figure 2, sampling 𝑓 implies evolving the random walker’s location ( P → ) and direction D → by sorting the photon flight length (equation (3.5)), then sorting the possibility of extending the random walk (if R < ω ) by a subsequent scattering and extinction event (through equation (3.9)). In Figure 2, P s − 1 represents the phase function after being integrated and inverted (from which μ s is sampled), while f → represents the coordinate change that corrects the direction of propagation D → . This should be noted as the sampled variables refer to a coordinate system centered at the random walker’s current position and direction ( μ = 0 points along the propagation direction, before scattering). Thus, a change of coordinates must be made to account for the random walker’s evolution on the global coordinate system.

Figure 2

MCRT algorithm for the evolution of a random walker in phase space (position P → and direction D → ) in an unbounded domain. The 𝑅’s represent uniformly distributed random numbers between 0 and 1, 𝐿 is the photon flight length, μ s is the cosine of the scattering angle.