## Abstract

An analytic solution to the problem of determining photon direction after successive scatterings in an infinite, homogeneous, isotropic medium, where each scattering event is in accordance with a two-term Henyey-Greenstein phase function, is presented and compared against Monte Carlo simulation results. The photon direction is described by a probability density function of the dot product of the initial direction and the direction after multiple scattering events, and it is found that such a probability density function can be represented as a weighted series of one-term Henyey-Greenstein phase functions.

© 2008 Optical Society of America

1. Introduction

Many fields in biomedical optics, remote sensing, and astrophysics rely upon the knowledge of the propagation of photons within a scattering medium. The exact distribution of photons within such a medium is difficult to determine, and simplifications are used to render the problem tractable. One such simplification is to use the Henyey-Greenstein phase function [1] to describe the single-particle scattering angle of photons. A linear combination of Henyey-Greenstein phase functions is often used to provide a better fit to real-world data [2].

In many applications, Monte Carlo simulation is performed to determine the exact photon trajectories and distributions for arbitrary media [3]. However, such simulation is computationally expensive and may not be suitable for large orders of scattering.

In imaging applications, the information carrying ability of a photon can be related to the photon’s direction **Ω**, a unit vector of the instantaneous trajectory [4]. A photon that has not deviated far from its initial direction will carry more information than a photon that has wandered far afield. Thus, knowledge of a photon’s direction **Ω** with respect to its initial direction **Ω**
_{0} after successive scattering events is important. For an isotropic medium, the photon’s direction with respect to its initial direction is conveniently defined by the cosine of the angle between them, hereafter referred to as the scattering cosine, cos *θ*=**Ω**
_{0}·**Ω**.

While performing Monte Carlo simulations of successive order, multiple scattering using Photon Transport Simulator software previously developed [4], it was noticed that the mean scattering cosine <cos *θ*> followed a definite pattern as the scattering order *n* was increased. Specifically, it was found that with a one-term Henyey-Greenstein (OTHG) phase function, the mean scattering cosine decreased by a multiplier of the anisotropy factor *g* at each successive scattering order. Such behavior is consistent with reported features of the OTHG phase function [5][6].

If a single scattering phase function is expanded using the Legendre polynomials [7], then the phase function for the *n*
^{th} scattering can also be expanded using the Legendre polynomials. This implies that the probability density function of the scattering cosine after multiple scatterings, when the single particle scattering is defined by a linear combination of Henyey-Greenstein phase functions, is also a linear combination of Henyey-Greenstein phase functions.

We derive an expression for the probability density function of the scattering cosine after multiple scatterings, when the single particle scattering is defined by a two-term Henyey-Greenstein (TTHG) phase function, and we compare these results against Monte Carlo simulations.

## 2. Single-particle probability density distributions

#### 2.1 One-term Henyey-Greenstein phase function

The one-term Henyey-Greenstein phase function p_{hg}(cos *θ*) describes the probability density that a photon with direction **Ω**
_{0} will scatter to direction **Ω** with scattering cosine cos *θ* [8].

When expanded in terms of Legendre polynomials [5], it is

Parameter *g*, also referred to as the anisotropy factor, has the property that it is also the mean cosine of the scattering angle <cos *θ*>. By definition, g can vary between -1 and 1, and it can be readily seen that when *g*=-1, there is complete backscattering (**Ω**=-**Ω**
_{0}); when *g*=0, there is isotropic scattering (**Ω** is unrelated to **Ω**
_{0}); and when *g*=1, there is no scattering (**Ω**=**Ω**
_{0}).

The OTHG phase function does not have a physical basis, but is often used to provide an empirical fit. It has been shown to have good agreement in applications where photons are collected in the far-field and have undergone successive order, multiple scattering.

## 2.2 Two-term Henyey-Greenstein phase function

For many applications, the OTHG phase function does not provide sufficient accuracy. In these cases a two-term Henyey-Greenstein phase function is often used. Such a phase function has been shown to accurately represent photon transport through tissue [9] and the ocean [10].

The TTHG phase function p_{hg}
_{2}(cos *θ*) is comprised of two weighted OTHG phase functions, each with different anisotropy factors *g*
* _{α}*,

*g*

*, such that*

_{β}*α*+

*β*=1.

When expanded in terms of Legendre polynomials, it is

## 3. Successive scattering probability density distributions

#### 3.1 One-term successive scattering probability density function

The *n*
^{th} order phase function describing the scattering cosine after *n* scattering events for a OTHG phase function has been shown to be [5]

By inspection, it can been seen that the probability density function p_{hg},_{n}(cos *θ*, *g*) has such a form for *n*=0, 1 and ∞. For *n*=0, the initial photon direction, *g*
^{0}=1. For *n*=1, the first scattering event, *g*
^{1}=*g*. For, n=∞, fully scattered, *g*
^{∞}=0. It is evident that the mean scattering cosine after *n* successive scatterings is *g*
* ^{n}*.

## 3.2 Two-term successive scattering probability density function

In a similar manner, the *n*
^{th} order phase function describing the scattering cosine for a TTHG phase function can be expanded as

The TTHG phase function p_{hg}
_{2}(cos *θ*) is comprised of the summation of two OTHG phase functions with weights *α* and *β*. At the first scattering event, a photon has probability *α* of scattering according to p_{hg}(cos *θ*, *g*
* _{α}*) and probability

*β*of scattering according to p

_{hg}(cos

*θ*,

*g*

*).*

_{β}At the second scattering event, a photon that scattered according to p_{hg}(cos *θ*, *g*
* _{α}*) still has probability

*α*of scattering according to p

_{hg}(cos

*θ*,

*g*

*) again. The probability that the photon will scatter twice in succession according to p*

_{α}_{hg}(cos

*θ*,

*g*

*) is*

_{α}*α*

^{2}. In like manner, the probability that a photon will scatter twice in succession according to p

_{hg}(cos

*θ*,

*g*

*) is*

_{β}*β*

^{2}. The probability that a photon will scatter first according to p

_{hg}(cos

*θ*,

*g*

*) and then according to p*

_{α}_{hg}(cos

*θ*,

*g*

*) is*

_{β}*αβ*. The probability density function p

_{hg2},

_{2}(cos

*θ*) is the weighted sum of the four individual contributions.

$$+\beta \alpha {p}_{\mathrm{hg},\beta \alpha}(\mathrm{cos}\theta ,{g}_{\beta \alpha})+{\beta}^{2}{p}_{\mathrm{hg},\beta \beta}(\mathrm{cos}\theta ,{g}_{\beta})$$

By equating *g*
* _{α}*=

*g*

*=*

_{β}*g*and comparing to the OTHG phase function, it is easily shown that

*g*

*=*

_{αβ}*g*

*=*

_{βα}*g*

_{α}*g*

*.*

_{β}Equations (8) to (11) show simplified forms for p_{hg}
_{2},* _{n}*(cos

*θ*) for

*n*=0 to 3, based on the above analysis.

$$+{\beta}^{2}{p}_{\mathrm{hg}}(\mathrm{cos}\theta ,{{g}_{\beta}}^{2})$$

$$+3\alpha {\beta}^{2}{p}_{\mathrm{hg}}(\mathrm{cos}\theta ,{g}_{\alpha}{{g}_{\beta}}^{2})+{\beta}^{3}{p}_{\mathrm{hg}}(\mathrm{cos}\theta ,{{g}_{\beta}}^{3})$$

It can be seen that after *n* successive scatterings, the probability that a photon has scattered according to a particular sequence of applications of p_{hg}(cos *θ*, *g*
* _{α}*) and p

_{hg}(cos

*θ*,

*g*

*) is*

_{β}*α*

^{a}*β*

*, where*

^{b}*a*+

*b*=

*n*, and that binomial coefficients specify the number of possible paths to get to order

*n*that contain only

*i*occurrences of

*α*branches. The final form of p

_{hg}

_{2},

*(cos*

_{n}*θ*) is thus:

It can be seen that p_{hg}
_{2},* _{n}*(cos

*θ*) is a weighted summation of

*n*+1 terms of OTHG phase functions. In a similar manner, the mean scattering cosine after

*n*successive scatterings is

## 4. Monte Carlo simulation

#### 4.1 Method

In order to test probability density function p_{hg2},* _{n}*(cos

*θ*) and confirm correct behavior with p

_{hg},

*(cos*

_{n}*θ*), the Photon Transport Simulator software was modified to record photon direction information after each scattering event. Each photon was scattered up to 100 times according to either the OTHG or TTHG phase function. After each scattering event, a photon density map (a table of scattering cosine versus scattering order) was updated. Once all photons had been launched, the photon density map was post-processed to determine the scattering cosine probability density distribution as a function of scattering order. All simulations were run with 10

^{8}photons on a 3 GHz PC operating under Windows XP.

## 4.2 Photon scattering

For the OTHG phase function, the azimuthal angle *φ* is uniformly distributed between 0 and 2π, while the scattering angle *θ* follows the p_{hg}(cos *θ*) probability distribution.

To determine cos *θ*, p_{hg}(cos *θ*) was sampled with a uniform random variable *ξ*∈[0..1].

For the TTHG phase function, the probability that the photon is scattered according to p_{hg}(cos *θ*, *g*
* _{α}*) is

*α*and according to p

_{hg}(cos

*θ*,

_{g}

*) is*

_{β}*β*=1-

*α*. Therefore, first a uniform random variable

*χ*∈[0..1] was used to determine which phase function to use [11], and then the selected phase function was sampled with a uniform random variable

*ξ*.

## 4.3 Recording of results

A two-dimensional table, comprised of 100 scattering orders *n* and 804 bins of scattering cosine cos *θ*, was created to record photon density with respect to scattering order and scattering cosine. The cos *θ* bins were more finely grained as cos *θ*→1 in order to capture photons sufficiently accurately when *g*→1.

After each scattering event, a single bin of the photon density map was incremented by one. The bin selected was based on the dot product of the photon’s new direction **Ω** with the photons initial direction **Ω**
_{0} and the number of scattering events *n* for that photon.

## 5. Results

#### 5.1 One-term successive scattering

Monte Carlo simulation of a OTHG phase function using a wide variety of anisotropy factors *g* showed good correlation between p_{hg},* _{n}*(cos

*θ*) and the recorded values. Representative results are shown in Figs. 1 and 2.

Figure 1 shows the Monte Carlo simulation results for a OTHG phase function with *g*=0.9 at scattering orders of *n*=1..4. This anisotropy factor was selected as it is often used to represent many common biological media. It can be seen that the probability density distribution p_{hg},* _{n}*(cos

*θ*) from Eq. (5) closely matches the density distribution obtained from the simulation.

Figure 2 shows results for higher scattering orders of *n*=10, 20, .., 50, and has an expanded vertical scale. As expected, as *n* increases, the scattering cosine more closely becomes uniform and p_{hg},* _{n}*(cos

*θ*) approaches a value of ½.

The correlation coefficient *r* between the Monte Carlo simulation data and p_{hg},* _{n}*(cos

*θ*) is greater than 0.9999 for

*n*≤10 and greater than 0.92 for

*n*≤50. Limited testing shows that simulating larger numbers of photons improves

*r*by reducing statistical fluctuations.

## 5.2 Two-term successive scattering

Monte Carlo simulation of a TTHG phase function with a range of weights *α*, *β*, and anisotropy factors *g*
* _{α}*,

*g*

*showed good correlation between p*

_{β}_{hg}

_{2},

*(cos*

_{n}*θ*) and the recorded values. Representative results are shown in Figs. 3 and 4.

Figure 3 shows the Monte Carlo simulation results for a TTHG phase function with *α*=0.9, *g*
* _{α}*=0.9,

*β*=0.1,

*g*

*=0 at scattering orders*

_{β}*n*=1..4 and 10. This phase function represents scattering through human dermis [9].

Figure 4 shows the Monte Carlo simulation results for a TTHG phase function with *α*=0.96, *g*
* _{α}*=0.9,

*β*=0.04,

*g*

*=-0.26 at scattering orders*

_{β}*n*=1..4 and 10. This phase function has significant backscattering and matches the first three moments of the seawater-Petzold particle phase function [12].

It can be seen that the p_{hg2},* _{n}*(cos

*θ*) from Eq. (12) closely matches the density distribution obtained from the simulations. The correlation coefficient r is greater than 0.995 for all results shown.

## 6. Discussion and conclusions

The results from Monte Carlo simulations agree with the analytic probability distribution functions p_{hg},* _{n}*(cos

*θ*) and p

_{hg2},

*(cos*

_{n}*θ*) and demonstrate that the probability density function for the n

^{th}scattering cosine can be calculated directly, without resorting to numerical simulation.

The fact that linear combinations of single-particle OTHG phase functions can be expanded as linear combinations of Legendre polynomial series provides the theoretical basis for the fact that the n^{th} order TTHG phase function can be exactly represented as a weighted sum of OTHG phase functions.

We have shown that for a TTHG phase function, the number of weighted terms in the probability density function p_{hg2},* _{n}*(cos

*θ*) is equal to the scattering order

*n*plus one - as the complexity of the solution is linear with scattering order, it is possible to calculate p

_{hg2},

*(cos*

_{n}*θ*) rapidly, even as

*n*becomes large. This is an improvement over the 2

^{n}terms suggested by Eq. (6).

Future work will investigate applying similar methods to larger linear combinations of Henyey-Greenstein phase functions.

## References and links

**1. **L.G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. **93**, 70–83 (1941). [CrossRef]

**2. **C. D. Mobley, L. K. Sundman, and E. Boss, “Phase function effects on oceanic light fields,” Appl. Opt. **41**, 1035–1050 (2002). [CrossRef] [PubMed]

**3. **L. Wang and H. Wu, *Biomedical Optics: Principles and Imaging* (John Wiley and Sons, 2007), *Ch. 3*.

**4. **N. Pfeiffer and G.H. Chapman, “Monte Carlo Simulations of the Growth and Decay of Quasi-Ballistic Photon Fractions with Depth in an Isotropic Medium,” Proc. SPIE **5695**, 136–147 (2005). [CrossRef]

**5. **I. Turcu and R. Bratfalean, “Narrowly peaked forward light scattering on particulate media I. Assessment of the multiple scattering contributions to the effective phase function,” J. Opt. A: Pure Appl. Opt.10, (2008). [CrossRef]

**6. **W. E. Vargas and G. A. Niklasson, “Forward-scattering ratios and average pathlength parameter in radiative transfer models,” J. Phys.: Condens. Matter **9**, 9083–9096 (1997). [CrossRef]

**7. **H. C. van de Hulst, *Multiple Light Scattering*, Vol 1 (Academic, 1980).

**8. **T. Binzoni1, T. S. Leung, A. H. Gandjbakhche, D. Rüfenacht, and D. T. Delpy, “The use of the Henyey-Greenstein phase function in Monte Carlo simulations in biomedical optics,” Phys. Med. Biol. **51**, N313–N322 (2006). [CrossRef]

**9. **S. L. Jacques, C. A. Alter, and S. A. Prahl, “Angular Dependence of HeNe laser Light Scattering by Human Dermis,” Laser Life Sci. **1**, 309–333 (1987).

**10. **V. I. Haltrin
, “Two-term Henyey-Greenstein light scattering phase function for seawater,” in *IGARSS ’99: Proceeding of the International Geoscience and Remote Sensing Symposium*,
T. I. Stein
, ed. (IEEE, 1999), pp. 1423–1425.

**11. **S. A. Prahl, *Light Transport in Tissue*, App. A1, PhD Thesis, (University of Texas at Austin, 1988).

**12. **R. A. Leathers and N. J. McCormick, “Ocean inherent optical property estimation from irradiances,” Appl. Opt. **36**, 8685–8698 (1997). [CrossRef]