On sampling of scattering phase functions

TL;DR

This paper reviews and improves sampling algorithms for scattering phase functions in Monte Carlo radiative transfer, enabling accurate simulation of complex phase functions and introducing a new phase function for realistic dust modeling.

Contribution

It introduces enhanced sampling methods, including piecewise linear approximations and Gibbs sampling, for complex phase functions, and proposes a new phase function for dust scattering.

Findings

Improved sampling performance for small scattering angles.

Accurate sampling of complex analytic phase functions like Fournier-Forand.

New phase function fits realistic dust scattering data.

Abstract

Monte Carlo radiative transfer, which has been demonstrated as a successful algorithm for modeling radiation transport through the astrophysical medium, relies on sampling of scattering phase functions. We review several classic sampling algorithms such as the tabulated method and the accept-reject method for sampling the scattering phase function. The tabulated method uses a piecewise constant approximation for the true scattering phase function; we improve its sampling performance on a small scattering angle by using piecewise linear and piecewise log-linear approximations. It has previously been believed that certain complicated analytic phase functions such as the Fournier-Forand phase function cannot be simulated without approximations. We show that the tabulated method combined with the accept-reject method can be applied to sample such complicated scattering phase functions accurately. Furthermore, we introduce the Gibbs sampling method for sampling complicated approximate analytic phase functions. In addition, we propose a new modified Henyey-Greenstein phase function with exponential decay terms for modeling realistic dust scattering. Based on Monte Carlo simulations of radiative transfer through a plane-parallel medium, we also demonstrate that the result simulated with the new phase function can provide a good fit to the result simulated with the realistic dust phase function.

Figures

Equations

$$p(\Omega ,{\Omega }^{\prime })=\frac{1}{{\sigma }_{sc}(\Omega )}\frac{d{\sigma }_{sc}(\Omega )}{d{\Omega }^{\prime }}$$

$$p_{HG}(\theta |g)=\frac{1}{2}\frac{1-g^2}{(1+g^2-2g\cos \theta {)}^{3/2}}$$

$$p_D(\theta |g,\alpha )$$

$$n(r)=A{r}^{-\alpha },r\ge r_{\min }$$

$$p(\theta )=\frac{{\int }_{r_{\min }}^{\infty }n(r){\sigma }_{sc}(r)p(\theta ,r)dr}{{\int }_{r_{\min }}^{\infty }n(r){\sigma }_{sc}(r)dr}$$

$$p_{FF}(\theta )$$

$$p_{FF}^{\prime }(\theta )$$

$$\mu =F^{-1}(\xi ),\text{ and }F(\mu )={\int }_{-1}^{\mu }p({\mu }^{\prime })d{\mu }^{\prime }$$

$$\mu$$

$$\varphi$$

$$\mu =\frac{{\mu }_{k+1}-{\mu }_k}{F_{k+1}-F_k}(\xi -F_k)+{\mu }_k$$

$$p(\mu )\approx \sum _{k=1}^K{c}_k\mathbb{I}(\mu \in A_k)$$

$$c_{k}=\left\{\begin{array}[]{rl}(F_{k+1}-F_{k})/(\mu_{k+1}-\mu_{k})&\text{if the CDF is known},\\ (p(\mu_{k+1})+p(\mu_{k}))/2s&\text{otherwise}\end{array}\right.$$

$$p(\mu )\approx \sum _{k=1}^K(a_k\mu +b_k)\mathbb{I}(\mu \in A_k)$$

$$\frac{1}{2}{a}_k{\mu }^2+b_k\mu -\frac{1}{2}{a}_k{\mu }_k^2-b_k{\mu }_k=\xi$$

$$p(\mu ,u|\alpha ,g)\propto \mathbb{I}(u<1+\alpha {\mu }^2)\frac{1-g^2}{2(1+g^2-2g\mu {)}^{3/2}}$$

$$p(u|\mu ,\alpha ,g)$$

$$p(\mu |u,\alpha ,g)$$

$$f(\cos \theta |g,a,b,k,k^{\prime })=C\frac{(1-g^2)(1+a\exp (-k\theta )+b\exp (-k^{\prime }(\pi -\theta ))}{(1+g^2-2g\cos \theta {)}^{3/2}}$$

$$p(\mu ,u|g,a,b,k,k^{\prime })\propto \mathbb{I}(u<1+a\exp (-k\theta )+b\exp (-k^{\prime }(\pi -\theta )))\frac{1-g^2}{2(1+g^2-2g\mu {)}^{3/2}}$$

Full text

K****eywords Monte Carlo radiative transfer  $\cdot$ Sampling method

1 Introduction

Monte Carlo radiative transfer (MCRT) [1, 2, 3] is a widely used numerical algorithm for studying transport of radiation in interstellar mediums [4], and the planetary atmosphere [5]. The MCRT primarily relies on producing pseudo random numbers to evaluate the integral form of the radiative transfer equation with statistical simulation. Considerable efforts have been devoted to produce sequences of random numbers with high qualities in order to implement a Monte Carlo simulation. In MCRT, a key step is to sample scattering directions from the phase function, which plays a fundamental role in understanding the physical properties of the medium [6, 7, 8, 9, 10, 11, 12, 13, 14]. ]. Hence, sampling of scattering phase functions is related to the accuracy of the simulated results in MCRT, and consequently, efficient sampling algorithms are particularly essential and in high demand.

For an approximate analytic phase function, it is effective to conduct simulations using the direct method [15, 2] by performing an inverse transform of the cumulative distribution function (CDF) corresponding to the scattering phase function. For realistic phase functions such as the Mie’s phase function, an analytic form of their inverse CDF may not exist or computing their values may be too expensive. Therefore, sampling has to be conducted in different ways. For example, the tabulated method [16, 17] is often used for sampling such complicated phase functions by constructing a look-up table, where random numbers are obtained by linear interpolation. As an approximate sampling method, the performance of the tabulated method is dependent on the number of interpolation points. However, to sample realistic phase functions (with peaks near the forward and backward directions) accurately, a considerably greater number of points are required to obtain a reliable result [18]. Poor sampling near forward and backward directions has a large negative impact on the simulation accuracy of radiation transport [19, 18]. The accept-reject method [15, 20] is another available method. For the Rayleigh phase function, Frisvad [21] reviewed several sampling methods and concluded that the simple accept-reject method was the most efficient. For realistic phase functions [7, 8, 9, 10, 13] with a highly forward scattering peak, the simple accept-reject method poses certain limitations owing to the high rejection rate [22].

In this paper, we propose several improvements for the widely used simple tabulated method for sampling realistic phase functions. Although some of these improvements may have been previously designed for general purpose [22] and may be known to some researchers in the radiative transfer community, there have been no detailed descriptions on them in the literature, particularly on how to correctly sample scattering phase functions with these methods. Moreover, several complicated approximate phase functions such as the Fournier-Forand (FF) phase function [23, 24] have often been considered to be impossible to sample accurately. In this work, we show that these phase functions could be accurately simulated with the proposed hybrid sampling method. In addition, Gibbs sampling [15], a Markov chain sampling method, is also suggested for generating random numbers for sampling complicated approximate analytic phase functions. As examples, detailed Gibbs sampling procedures are described for the Draine phase function [25] and a new modified Henyey-Greenstein phase function with two exponential decay terms. This new phase function is also validated based on Monte Carlo simulations. This work is expected to contribute to the fields of astronomy and earth sciences by providing detailed descriptions on how to sample difficult scattering phase functions in applications of Monte Carlo radiative transfer.

2 Phase function

Light scattering by small particles [26] is characterised by the scattering cross section and the scattering phase function. The scattering phase function is defined as the normalised differential scattering cross section [26]:

[TABLE]

where $\Omega$ and $\Omega^{\prime}$ are the incident and scattered directions of light. Under this definition, the scattering phase function $p(\Omega,\Omega^{\prime})$ describes the angular distribution of scattered light for a specific wavelength. If the medium is assumed to be isotropic, the scattering phase function is a function of the angle $\theta$ between $\Omega$ and $\Omega^{\prime}$ only, i.e., $\cos\theta$.

The scattering phase function for small particles can be obtained numerically by solving Maxwell’s equations. Various numerical methods [26] have been proposed to address the scattering problem. However, many approximate analytic phase functions are still widely employed owing to their simplicity and usefulness. One such phase function proposed by Henyey-Greenstein n[27] (HG) to model scattering of dust can be expressed as

[TABLE]

where $-1\leq g\leq 1$. This single parameter phase function has been used to simulate the scattering properties of various particles. The parameter $g$ can be determined by the first moment $g=<\cos\theta>$. Although the HG phase function is quite successful, it fails to reproduce observations. Various numerical and experimental studies have shown that the HG phase function underestimates the possibility of small and large angle scattering. To address this problem, other scattering phase functions [28, 29, 30] were established such as the one proposed by Draine:

[TABLE]

with the normalization constant $C=\frac{1}{1+\alpha(1+2g^{2})/3}$. The Draine phase function tends to the Rayleigh phase function at long wavelengths and reduces to the phase function proposed by Cornette & Shanks[31] when $\alpha=1$.

A more sophisticated approach for modelling realistic phase functions can be established by assuming a particle size distribution and computing the ensemble average. For example, the power-law distribution for atmospheric and marine particles [24] can be expressed as:

[TABLE]

where $\alpha$ is a parameter, $r$ is the radius of the particle volume-equivalent sphere, $r_{\min}(>0)$ is the lower bound of $r$, and $A=(\alpha-1)r_{\min}^{\alpha-1}$ is the normalisation constant. It should be noted that $r_{\min}$ is necessary to ensure that the distribution is reasonable, which is often ignored in the literature. Then, the ensemble average scattering phase function can be obtained by calculating the ratio of ensemble average differential scattering cross section to the ensemble average scattering cross section:

[TABLE]

Based on the power-law particle size distribution, Fournier and Forand presented a two-parameter approximate analytic phase function [24] where each particle scatters light with an anomalous diffraction approximation. The latest form of the FF phase function can be written as

[TABLE]

where $\nu=\frac{3-\alpha}{2}$ and $\delta=\frac{4}{3(m-1)^{2}}\sin^{2}(\theta/2)$, $\alpha$ is the exponent of the power-law distribution, $m$ is the real part of the relative refractive index of the particle. The FF phase function can successfully model a marine environment. However, the value of $p_{FF}(\theta)$ at $\theta=0^{\circ}$ can not be defined, which is rather peculiar. Hence, we suggest a regularised FF phase function:

[TABLE]

where $\mathbb{P}(\theta\leq\theta_{0})$ represents the probability that the scattering angle is smaller than or equal to $\theta_{0}$ and $\mathbb{I}(x\in S)$($S$ is a set.) is the indicator function, which is 1 if $x\in S$ and 0 otherwise .

3 Sampling Methods

3.1 Inverse CDF method

Monte Carlo simulations are based on the generalisation of random samples that are distributed according to a certain distribution. Samples from nonuniform distributions can be generated by transforming sequences of uniform random numbers on $(0,1)$. If an inverse CDF exists, a continuous random variable $\mu$ can be generated by taking the transform of a uniform random variable $\xi$ on $(0,1)$,

[TABLE]

where $\mu$ is a random variable with the desired distribution, $p(\mu)$ denotes the probability density (in this setting, the scattering phase function) and $F(\mu)$ denotes the CDF. For example, the HG phase function can be sampled easily with the inverse CDF technique:

[TABLE]

where $\mu=\cos\theta$, and $\xi_{1}$ and $\xi_{2}$ are random numbers uniformly distributed on $(0,1)$.

3.2 Accept-reject method

It is sometimes difficult or even impossible to obtain an analytical form of the inverse CDF for certain phase functions. In such cases, the accept-reject method can be used, in which simulate the desired density $f_{X}$ (scattering phase function) is simulated using the similar but simpler density function $f_{Y}$.

In the accept-reject method, we first sample a random number $x$ from density $f_{Y}$, then draw a uniform random number $\xi\sim\mathcal{U}(0,1)$, and accept $x$ if $\xi<f_{X}(x)/(Mf_{Y}(x))$ ( $M$ can be chosen as the maximum of $f_{X}/f_{Y}$.). Otherwise, we return to the first step. Iterations of this procedure will produce a sequence of random numbers from the density $f_{X}$. The average number of iterations for generating a random variate is proportional to $M$ ($M$ is larger than 1.). In general, if $f_{X}(x)$ can be written as $f_{X}(x)=Cg(x)f_{Y}(x)$, where $f_{Y}(x)$ is another density function and $C$ is a positive constant that makes $g(x)\in[0,1]$, random numbers from $f_{X}(x)$ can be generated by the following procedure:

3.3 Tabulated method

It is often difficult to simulate common phase functions for realistic particles using the simple accept–reject method owing to the highly forward scattering peak. The tabulated method is usually applied for such cases where a look-up table is constructed and then random numbers are generated using linear interpolation. In the following section, we present several tabulated methods for accurately sampling random scattering directions in Monte Carlo radiative transfer.

In the tabulated method, the range $[0,1]$ of the CDF for the scattering phase function is divided uniformly into subintervals. Hence, the sampling space is also divided into subsets. First, we draw a uniform random number $\xi\sim\mathcal{U}(0,1)$ first. Then, the subinterval where $\xi$ is located is determined, $\xi\in[F_{k},F_{k+1})$ ($\{F_{k}\}$ represents values of the CDF.). A random number $\mu$ from the desired scattering phase function can be found by linear interpolation:

[TABLE]

where $\mu_{k}$($-1=\mu_{1}<\mu_{2}<...<\mu_{K+1}=1$) represents the value corresponding to $F_{k}$. The index of the interval can be found using the bisection search algorithm[32]. It should be noted that in this widely used tabulated method, the scattering phase function is actually approximated by a piecewise constant function:

[TABLE]

Here, $\mathbb{I}(\cdot)$ is the indicator function, $A_{k}$ represents the interval $(\mu_{k},\mu_{k+1})$ and

[TABLE]

with $s=\sum_{k}c_{k}(\mu_{k+1}-\mu_{k})$. As an example, Fig. 1 illustrates a normalized histogram using samples generated with the tabulated method and piecewise constant approximation (PCA). However, the quality of random numbers generated from the tabulated method above depends on the number of subintervals. Moreover, PCA sometimes is not accurate enough to expand a regular phase function near the forward and backward scattering regions. The phase function can be better approximated by a piecewise linear function(PLA) that agrees well with the actual value at each of the points $-1=\mu_{1}<\mu_{2}<...<\mu_{K+1}=1$. Based on this observation, improvements can be made by extending PCA to PLA or the piecewise log-linear approximation (PLLA). For PLA, the corresponding approximate phase function can be represented as

[TABLE]

where $a_{k}=(p(\mu_{k+1})-p(\mu_{k}))/(\mu_{k+1}-\mu_{k})$ and $b_{k}$.The sampling procedure is similar to the one described above, except the following quadratic equation is solved:

[TABLE]

The procedure for sampling is shown in Algorithm 1. Fig. 2 illustrates a normalised histogram of $\mu$ sampled from PLA of a HG phase function. And PLA outperforms PCA, as shown in Fig. 3 and Fig. 4. Further acceleration, as shown in Algorithm 2, can be made if we select the values of the CDF to be uniform on $[0,1],$[17, 18]. By using such a setting, the index of subintervals can be found more easily. The trade-off would be intervals must be fine enough to approximate the scattering phase function near the backward scattering region.

Another improvement was established in this study by dividing the entire interval of $\mu\in[-1,1]$ into equal probable subintervals similar to the tabulated method; however, for each subinterval, the accept-reject method is applied instead of interpolation. In particular, if the analytic form of the phase function is known, the phase function can be sampled accurately using this approach. The procedure for the hybrid tabulated method combined with the accept-reject method is shown in Algorithm 3. Fig. 5 compares the sample quality of the hybrid and the simple tabulated method with normalised histograms of random $\cos\theta$ generated from the FF phase function. As shown in Fig. 5, the FF phase function, which has long been considered impossible to exactly sample, can be sampled accurately using the hybrid tabulated method.

3.4 Gibbs sampling

Tabulated methods are simple and easy to implement; however, the accuracy of the simulation depends on the number of interpolation points. Moreover, tabulated methods are computationally expensive, and they consume large amounts of memory to be implemented on GPUs. In the following section, Gibbs sampling, a sampling method using Markov chains is introduced for generating random numbers from certain complicated analytic scattering phase functions.

The Gibbs sampling method relies on building a joint density with conditional distributions that are easy to simulate. Unlike other methods, auxiliary random variables are required in for Gibbs sampling. Random numbers are iteratively generated from a sequence of conditional distributions. The limit distribution of this sequence converges to the desired density

As an example, the Draine phase function is sampled using the Gibbs sampling method. The joint density required to be sample is

[TABLE]

where $u$ is an auxiliary random variable. The two conditional densities are

[TABLE]

The sampling procedure is summarised in Algorithm 4.

The Fig.  6 compares the Draine phase function with different parameters and Fig. 7 shows normalised histograms for the Draine phase function samples generated using the Gibbs sampling method. The autocorrelation functions of random scattering angles are plotted in Fig. 8; from this figure, we can see that the autocorrelation of Lag-1 is considerably below 0.1. As expected, the Gibbs sampling method is an extremely powerful method that allows us to sample complicated analytical phase functions.

4 Simulation and Results

Analytic phase functions such as the Draine phase function can provide good approximations to realistic dust scattering at long wavelengths. For a higher frequency case, the true phase function of particles exhibits strong peaks (both forward and backward) and modelling using these analytic phase functions can result in noticeable errors. Considering two peaks of the true phase functions at the forward and backward directions, a new phase function is proposed with exponential decay terms:

[TABLE]

where $g$, $a$, $b$, $k$ and $k^{\prime}$ are certain parameters, and $C$ is the normalisation constant that can be determined numerically.

Fig. 9 compares the ensemble averaged phase function (modelled with a rough spheroid [12] ) and the new phase function. The sampling process for the new phase function is similar, and the joint density is

[TABLE]

The sampling procedure is summarised in Algorithm 5. Fig. 10a illustrates the normalised histogram of the sample from the new phase function using the Gibbs sampling method and the autocorrelation function plot of samples from this new phase function is shown in Fig. 10b. The figure shows no obvious short-term autocorrelation in the scattering angle sample, which indicates the independence of these samples.

As a more realistic test, a backward Monte Carlo [5] simulation is conducted to compute the reflected and transmitted radiance through a plane-parallel medium. Each photon is traced from the detector to a light source in a backward manner. The radiance is calculated by summing the contribution of each photon. The albedo of the bottom boundary is set to 0, and only volume scattering events are considered here. The new phase function is used to model dust scattering and simulated using the Gibbs sampling method. The measured scattering phase function obtained from the Amsterdam Light Scattering Database is used as the true dust scattering phase function, sampled using Algorithm 1.

The left panel of Fig. 11 illustrates the relative error of the simulated transmitted radiance with the HG phase function and the new phase function. The error with the new phase function is below 2% for nearly all angles and approximately $1/2$ of the error with the HG phase function. The right panel of Fig. 11 illustrates the relative error for the simulated reflected radiance with these two phase functions. The result simulated with the new phase function also outperforms that with the HG phase function, except for only a small range of angles. For observation angles from $-30^{\circ}$ to $30^{\circ}$, the relative error corresponding to the new phase function is roughly constant with a value of approximately 1%. Conversely, the relative error corresponding to the HG phase function reaches a maximum at the backward direction and is approximately five times larger than the relative error corresponding to the new phase function. Consequently, results simulated with the new phase function fit the realistic dust scattering results nearly uniformly better than those simulated with the widely used HG phase function.

5 Discussion and Conclusion

In this paper, several new sampling methods were proposed for simulating the scattering phase function in the MCRT algorithm. We show that the commonly used tabulated method actually samples a PCA of the exact scattering phase function, which may not be a good choice for modelling realistic dust scattering. Improvements are made by extending the PCA to a PLA or a PLLA. This modification requires minor revisions in the code of the MCRT algorithm and is not computationally expensive. Furthermore, a hybrid method combining the tabulated method and the accept-reject method is proposed by substituting the interpolation procedure with an accept-reject procedure. As an example, the regularized FF function is sampled exactly with this hybrid method. In addition, the Gibbs sampling method is also recommended for generating random numbers from certain complicated analytic phase functions; this has previously not been used in MCRT. Based on the HG phase function, a new phase function with two exponential decay terms is also proposed to better model realistic dust scattering. We also demonstrate that this new phase function can be sampled efficiently using the Gibbs sampling method. The performance of this new phase function is validated for a plane-parallel medium simulated using the backward MCRT algorithm. The results simulated with the new phase function fit the realistic dust scattering model considerably better than those simulated with the HG phase function, particularly near the backscattering region. It is recommended that this new phase can be used for modelling radiation transport in dusty planetary atmospheres and interstellar mediums.

Acknowledgments

We would like to thank anonymous reviewers and the review editor for their helpful comments. Part of this work was supported by the National Natural Science Foundation of China (NSFC) (Grant No. 41705010) and the Fundamental Research Funds for the Central Universities (DUT19RC(4)039).