## Abstract

The propagation of light through turbid media is of fundamental interest in a number of areas of optical science including atmospheric and oceanographic science, astrophysics and medicine amongst many others. The angular distribution of photons after a single scattering event is determined by the scattering phase function of the material the light is passing through. However, in many instances photons experience multiple scattering events and there is currently no equivalent function to describe the resulting angular distribution of photons. Here we present simple analytic formulas that describe the angular distribution of photons after multiple scattering events, based only on knowledge of the single scattering albedo and the single scattering phase function.

©2011 Optical Society of America

## 1. Introduction

Many natural media, including the oceans, clouds, nebulae, skin and other living tissues, exhibit strong effects of multiple scattering. Optical text books often discuss single scattering theory in depth, but give only a cursory introduction to multiple scattering [1]. There are specialist books on multiple scattering [2–5] that treat the subject in much greater detail, but there remains a need for simple analytical relationships to describe the phenomenon. We present a number of new formulations relating multiple scattering signals to well-understood single scattering parameters that go some way to satisfying this demand.

The scattering phase function $\tilde{\beta}\left(\theta \right)$ describes the angular distribution of scattered photons in the case of single scattering. However, in many cases photons undergo more than one scattering event and there is a need to develop effective phase functions ${\tilde{\beta}}_{ms}\left(\theta \right)$ that describe the angular distribution of photons after multiple scattering events [6]. The asymmetry parameter for the single scattering phase function is given by

*n*order scattering has

^{th}*g*. They also showed that the resulting

_{n}= g_{1}^{n}*n*order phase function takes the form of a HG phase function.

^{th}Multiple scattering is a 3D process that cannot, in the general case, be reduced to 2D by assuming azimuthal symmetry. After two scattering events at *θ _{1}* and

*θ*, a photon will propagate at an angle

_{2}*θ*to the original incident direction, where

_{n}*ψ*is the azimuthal angle of scattering for the second scattering event relative to the first. After two scattering events the photon will have a direction of propagation between

*θ*and

_{1}- θ_{2}*θ*(Fig. 1 ). Here we assume that the medium is infinite, isotropic and homogeneous.

_{1}+ θ_{2}The purpose of this paper is to determine if the results of Pfeiffer and Chapman [7] are valid for the general case of phase functions that are not defined on the asymmetry parameter and to develop a robust method for determining an effective scattering phase function for the case of multiple scattering. In the first instance Monte Carlo simulations are used to provide fully parameterised light fields that can be used to test the above hypotheses and to establish general relationships that are applicable across the range of sciences where multiple scattering is a feature. Here we define effective scattering phase functions in the same sense as Pfeiffer and Chapman [7] as describing the angular distribution of scattered photons for n^{th} order scattering, as determined from Monte Carlo simulations. Furthermore, we define an effective phase function for multiple scattering which describes the angular distribution of photons after all orders of scattering have been simulated (limited only by the number of photons used in the Monte Carlo simulations). We believe this name represents a logical extension of the classical, single scattering phase function.

## 2. Methods

The forward directed Monte Carlo code has been extensively tested previously and used in both studies of marine environment [8,9] and for correction of measurement errors of ocean optics instrumentation [10,11]. Absorption and scattering parameters as well as the scattering phase functions were treated as statistical properties, which makes it possible to decide the fate of each virtual photon using random numbers. A pseudo-random number generator of tested uniformity [12] and a period of 2^{127} (KISS by George Marsaglia) was used. In the setup used for this study, the modelled light source was a parallel light beam placed in a virtually unbounded space (no virtual photons reached the 3D model walls which were several thousands of optical depths from the light source). The beam is assumed to be vertical for simplicity. For each run, one billion (10^{9}) virtual photons were used. The zenith angle of the virtual photon and the scattering order were recorded for each scattering event. This made possible the estimation of phase functions for each scattering order as well as the effective phase function (the latter using all scattering events).

The Henyey-Greenstein [13] scattering phase function is calculated from

*g*is the mean cosine of the scattering phase function. Mobley [14] showed that

*g*= 0.924 gave a best-fit to the ‘particle phase function’ derived from Petzold’s [15] measurements of the scattering phase function of seawater. Fournier and Forand [16,17] gave an alternative analytic expression for the scattering phase function based upon a population of Mie scatterers with a hyperbolic (Junge) particle size distribution (slope

*μ*) and a real refractive index,

*n*,

_{r}*n*and

*μ*, and integrating Eq. (4), Mobley et al [18] were able to obtain a formulation of the Fournier-Forand (FF) scattering phase function using only the particle backscattering ratio,

*B*=

_{p}*b*, as input. This formulation has been used widely in ocean optics and is adopted here. For the purpose of this paper the key difference between HG and FF scattering phase functions is that the former are defined on the asymmetry parameter, g, while the latter are not. This difference facilitates testing the general applicability of the Pfeiffer and Chapman [7] result for any arbitrary scattering phase function.

_{bp}/ b_{p}## 3. Results

Monte Carlo simulations were performed for HG scattering phase functions with a range of asymmetry values (0.5 < g_{1} < 0.99). In every case the angular distribution of *n ^{th}* order scattering was well described by a HG function with

*g*, closely reproducing the results of Pfeiffer and Chapman [7] (Fig. 2 ). To test the general validity of this result, we repeated the simulations with FF scattering phase functions that are used to describe scattering in marine environments and are not specifically parameterised on the asymmetry parameter [16,17]. Figure 3 shows that

_{n}= g_{1}^{n}*g*is true for all orders of scattering up to and including 50th order scattering (the result is only limited by the number of photons used in the simulation). We note that this result is consistent with taking the mean of both sides of Eq. (2) over all values of

_{n}= g_{1}^{n}*ψ*.

The single scattering albedo, *ω*, is defined as the ratio of scattering, *b*, to attenuation, *c*, where *c = a + b* and *a* is the absorption coefficient. From Monte Carlo simulations we were able to resolve the final angular distribution of photons that had experienced multiple scattering and hence could determine an effective scattering phase function ${\tilde{\beta}}_{ms}\left(\theta \right)$. Figure 4
shows that this effective scattering phase function has an asymmetry parameter *g _{ms}*, that is related to the single scattering asymmetry parameter

*g*through

_{1}*g*< 1 and 0 ≤ ω ≤ 1.

_{1}The probability of single scattering into an angular bin centred on *θ _{1}* is given by $\tilde{\beta}\left({\theta}_{1}\right)\mathrm{sin}{\theta}_{1}\mathrm{\Delta}{\theta}_{1}$. The phase function for

*n*order scattering can be calculated from the single scattering phase function by iterating through each order of scattering

^{th}*n,*using

*θ*determined for each combination of

_{n}*θ*,

_{1}*θ*,

_{2}*ψ*using Eq. (2). The validity of this parameterisation is demonstrated in Fig. 5a where

*n*order HG phase functions are accurately reproduced up to order 50 using Eq. (6) and using the result of Pfeiffer and Chapman [7] for HG functions. Figure 5b shows that Eq. (6) accurately reproduces

^{th}*n*order FF phase functions obtained from Monte Carlo simulations up to at least 50th order.

^{th}The effective scattering phase function for multiple scattering can be calculated from *n ^{th}* order phase functions using Eq. (7)

*n*order phase functions obtained from Monte Carlo simulations.

^{th}## 4. Conclusions

Multiple scattering is a phenomenon that affects the propagation of light through turbid media. Examples are easily found in marine and atmospheric optics, biomedical optics, astrophysics and plasma physics amongst others. Derivation of approximate analytical solutions for multiple scattering is an active area of research, with one potential benefit (under limited circumstances) being the ability to vastly reduce computation time in comparison to standard Monte Carlo approaches [18]. Equations (5), (6), and (7) are, to the best of our knowledge, new developments that provide excellent fits to simulated data and are consistent with rational analysis of the multiple scattering problem. The results presented in this paper are of general relevance and provide a simple framework in which the relationship between single scattering events and multiple scattering effects can be easily parameterised and understood. We note that the asymmetry parameter (also known as the mean cosine) is a particularly useful descriptor of the scattering phase function and further effort to parameterise popular phase functions (such as Fournier-Forand) in terms of *g* would be well placed. It is anticipated that these results will facilitate new insights into the effect of multiple scattering on a variety of applications such as generation of remote sensing signals, impact on optical measurements and interpretation of scattering signals from turbid media.

## Acknowledgments

The authors would like to thank the Royal Society of Edinburgh International Exchange Programme for funding a number of visits between laboratories. This work was also supported by the Natural Environment Research Council through the award of an Advanced Fellowship to D. M. (NE/E013678/1) and by statutory task I.3.3 of Institute of Oceanology PAS.

## References and links

**1. **H. C. van de Hulst, *Light Scattering by Small Particles* (Dover, 1981).

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

**3. **H. C. van de Hulst, *Multiple Light Scattering, Vol. 2* (Academic Press, 1980).

**4. **M. I. Mishchenko, L. D. Travis, and A. A. Lacis, *Multiple Scattering of Light by Particles: Radiative Transfer and Coherent Backscattering* (Cambridge Univ. Press, 2006).

**5. **A. A. Kokhanovsky, *Cloud Optics* (Springer, 2006).

**6. **M. Sydor, “Statistical treatment of remote sensing reflectance from coastal ocean water: proportionality of reflectance from multiple scattering to source function b/a,” J. Coast. Res. **235**, 1183–1192 (2007). [CrossRef]

**7. **N. Pfeiffer and G. H. Chapman, “Successive order, multiple scattering of two-term Henyey-Greenstein phase functions,” Opt. Express **16**(18), 13637–13642 (2008). [CrossRef] [PubMed]

**8. **P. Flatau, J. Piskozub, and J. R. V. Zaneveld, “Asymptotic light field in the presence of a bubble-layer,” Opt. Express **5**(5), 120–124 (1999). [CrossRef] [PubMed]

**9. **Z. Otremba and J. Piskozub, “Modelling the bidirectional reflectance distribution function (BRDF) of seawater polluted by an oil film,” Opt. Express **12**(8), 1671–1676 (2004). [CrossRef] [PubMed]

**10. **J. Piskozub, A. R. Weeks, J. N. Schwarz, and I. S. Robinson, “Self-shading of upwelling irradiance for an instrument with sensors on a sidearm,” Appl. Opt. **39**(12), 1872–1878 (2000). [CrossRef]

**11. **D. McKee, J. Piskozub, and I. Brown, “Scattering error corrections for in situ absorption and attenuation measurements,” Opt. Express **16**(24), 19480–19492 (2008). [CrossRef] [PubMed]

**12. **G. Marsaglia, *The Diehard Battery of Tests of Randomness* (1995), http://www.stat.fsu.edu/pub/diehard/.

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

**14. **C. D. Mobley, *Light and Water: Radiative Transfer in Natural Waters*. (Academic, 1994).

**15. **T. J. Petzold, T. J. SIO Ref. 72–78, Scripps Institute of Oceanography (U. California, 1972).

**16. **G. R. Fournier and J. L. Forand, “Analytic phase function for ocean water,” Proc. SPIE **2258**, 194–201 (1994). [CrossRef]

**17. **J. L. Forand and G. R. Fournier, “Particle distributions and index of refraction estimation for Canadian waters,” Proc. SPIE **3761**, 34–44 (1999). [CrossRef]

**18. **I. Turcu and M. Kirillin, “Quasi-ballistic light scattering – analytical models versus Monte Carlo simulations,” J. Phys.: Conf. Ser. **182**, 012035 (2009). [CrossRef]