Analytical solution of beam spread function for ocean light radiative transfer

We develop a new method to analytically obtain the beam spread function (BSF) for light radiative transfer in oceanic environments. The BSF, which is defined as the lateral distribution of the (scalar) irradiance with increasing depth in response to a uni-directional beam emanating from a point source in an infinite ocean, must in general be obtained by solving the three-dimensional (3D) radiative transfer equation (RTE). By taking advantage of the highly forward-peaked scattering property of the ocean particles, we assume, for a narrow beam source, the dependence of radiance on polar angle and azimuthal angle is deliberately separated; only single scattering takes place in the azimuthal direction while multiple scattering still occurs in the polar direction. This assumption enables us to reduce the five-variable 3D RTE to a three-variable two-dimensional (2D) RTE. With this simplification, we apply Fourier spectral method to both spatial and angular variables so that we are able to analytically solve the 2D RTE and obtain the 2D BSF accordingly. Using the relations between 2D and 3D solutions acquired during the process of simplification, we are able to obtain the 3D BSF in explicit form. The 2D and 3D analytical solutions are validated by comparing with Monte Carlo radiative transfer simulations. The 2D analytical BSF agrees excellently with the Monte Carlo result. Despite assumptions of axial symmetry and spike-like azimuthal profile of the radiance in deriving the 3D BSF, the comparisons to numerical simulations are very satisfactory especially for limited optical depths (< O(5)) for single scattering albedo values typical in the ocean. The explicit form of the analytical BSF and the significant gain in computational efficiency (several orders higher) relative to RTE simulations make many forward and inverse problems in ocean optics practical for routine applications. © 2015 Optical Society of America OCIS codes: (010.4450) Oceanic optics; (010.5620) Radiative transfer; (010.3310) Laser beam transmission. References and links 1. K. M. Case and P. F. Zweifel, Linear Transport Theory (Addison-Wesley, 1967). 2. S. Q. Duntley, “Underwater lighting by submerged lasers and incandescent sources,” Tech. Rep., Scripps Institution of Oceanography, San Diego, Calif. (1971). 3. K. J. Voss and A. L. Chapin, “Measurement of the point spread function in the ocean,” Appl. Opt. 29, 3638–3642 (1990). 4. K. J. Voss, “Simple empirical model of oceanic point spread function,” Appl. Opt. 30, 2647–2651 (1991). #239808 Received 27 Apr 2015; revised 25 Jun 2015; accepted 26 Jun 2015; published 1 Jul 2015 © 2015 OSA 13 Jul 2015 | Vol. 23, No. 14 | DOI:10.1364/OE.23.017966 | OPTICS EXPRESS 17966 5. W. Hou, D. J. Gray, A. D. Weidemann, and R. A. Arnone, “Comparison and validation of point spread models for imaging in natural waters,” Opt. Express. 16(13), 9958–9965 (2008). 6. C. D. Mobley, Light and Water: Radiative Transfer In Natural Waters (Academic, 1994). 7. G. N. Plass and G. W. Kattawar, “Monte Carlo calculations of radiative transfer in the earth atmosphere-ocean system: I. flux in the atmosphere and ocean,” J. Phys. Oceanog 2(13), 139–145 (1972). 8. R. Sanchez and N. J. McCormick, “Analytic beam spread function for ocean optics applications,” Appl. Opt. 41(30), 6276–6288 (2002). 9. E. P. Zege, A. P. Ivanov, and I. L. Katsev, Image Transfer Through a Scattering Medium (Springer-Verlag, 1991). 10. G. C. Pomraning, A. K. Prinja, and J. W. VanDenburg, “An asymptotic model for the spreding of a collimated beam,” Nucl. Sci. Eng. 112, 347–360 (1992). 11. C. Börgers and E. W. Larsen, “Asymptotic derivation of the Fermi pencil-beam approximation,” Nucl. Sci. Eng. 123, 343–357 (1996). 12. W. H. Wells, “Theory of small angle scattering,” in Optics of the Sea (Interface and In-Water Transmission and Imaging), P. Halley, ed., Lecture Series no. 61 (North Atlantic Treaty Organization, Advisory Group for Aerospace Research and Development, Neuilly-sur-Seine, France, 1973). 13. R. E. Walker, Marine Light Field Statistics (Wiley, 1994). 14. J. W. McLean, J. D. Freeman, and R. E. Walker, “Beam spread function with time dispersion,” Appl. Opt. 37(21), 4701–4711 (1998). 15. W. Wei, X. Zhang, Y. Chao, and X. Zhou, “An analytical model of the power spatial distribution for underwater optical wireless communication,” Optica Applicata XLII(1), 157–166 (2002). 16. N. L. Swanson, V. M. Gehman, B. D. Billard, and T. L. Gennaro, “Limits of the small-angle approximation to the radiative transport equation,” J. Opt. Soc. Am. A 18, 385–391 (2001). 17. N. L. Swanson, B. D. Billard, V. M. Gehman, and T. L. Gennaro, “Application of the small-angle approximation to ocean water types,” Appl. Opt. 40, 3608–3613 (2001). 18. V. A. Markel, “Modified spherical harmonics method for solving the radiative transport equation,” Waves Random Media 14(1), L13–L19 (2004). 19. A. Liemert and A. Kienle, “Green’s functions for the two-dimensional radiative transfer equation in bounded media,” J. Phys. A: Math. Theor. 45, 175201 (2012). 20. Z. Xu and D. K. P. Yue, “Monte Carlo radiative transfer simulation for the near-ocean-surface high-resolution downwelling irradiance statistics,” Opt. Eng. 53, 051408 (2014). 21. T. J. Petzold, “Volume scattering functions for selected ocean water,” Tech. Rep., Scripps Institution of Oceanography, San Diego, Calif. (1972). 22. L. C. Henyey, and J. L. Greenstein,“Diffuse radiation in the galaxy,” Astrophys. J. 93, 70–83 (1941). 23. L. Florescu, J. C. Schotland, and V. A. Markel, “Single-scattering optical tomography,” Phy. Rev. E 79, 036607 (2009). 24. J. P. Boyd, Chebyshev and Fourier Spectral Methods (DOVER Publications, 2000). 25. Z. Xu, D. K. P. Yue, L. Shen, and K. J. Voss,“Patterns and statistics of in-water polarization under conditions of linear and nonlinear ocean surface waves,” J. Geophys. Res. 116, C00H12 (2011). 26. Z. Xu, X. Guo, L. Shen, and D. K. P. Yue,“Radiative transfer in ocean turbulence and its effect on underwater light field,” J. Geophys. Res. 117, C00H18 (2012).


Introduction
Knowledge of the response to highly collimated light beam in an anisotropic scattering medium such as turbid ocean, atmosphere and human tissue is widely used in underwater optical imaging and communication, ocean light detection and ranging (LIDAR), graphic rendering, and optical diffuse tomography. In these contexts, the spatial distribution of irradiance (usually a lateral spreading with increasing depth) in response to a narrow uni-directional light beam is defined as the beam spread function (BSF) (also called the point spread function or PSF). In a macroscopical perspective, the light propagation in a scattering medium is fully described by the (3D) radiative transfer equation (RTE). Thus, the BSF generally requires the solution of the RTE, which is a significant task. Our interest is in applications in ocean optics where, due to the highly peaked anisotropic distribution of the scattering phase function for oceanic particles, the BSF for ocean light radiative transfer is particularly difficult to solve both analytically and numerically [1].
The BSF in ocean optics has traditionally been obtained through experiments [2,3] and empirical models [4,5]. Direct measurements of spatially distributed irradiance in the field is challenging and expensive. Empirical models, while practically efficient in some applications, generally lack rigorous derivation and are often unable to reveal the underlying physics and apparent optical properties. Direct numerical or Monte Carlo solutions of the RTE are also available [6,7], but the computational effort is generally large. Computational solution of the BSF directly has also been attempted. Sanchez and McCormick [8] utilize the discrete ordinate method for directional variables and finite difference in spatial variables to obtain the numerical BSF; however, the computational cost is still comparable to solving the RTE.
The analytical solution of the RTE Green's function is fundamental and has the advantage of direct physical insights and the potential of highly computational efficiency. For the BSF, most analytical work is based on the diffusion approximation (DA) [1,9] of the RTE. However, a diffusion approximation is mainly accurate for strong (single scattering albedo ω 0 0.99) and/or isotropic scattering. For highly forward peaked scattering, asymptotic models [10,11] have also been developed. These models assume that the scattering polar angle θ of radiance is highly peaked at small values so that μ = cos θ ≈ 1. Expanding the RTE in an asymptotic series for small θ obtains a sequence of advection-diffusion equations for different orders. The asymptotic model is valid when the leading term dominates which means very small optical depth and lateral spreading. Another important approach is the small-angle approximation (SAA) [12,13]. The basic idea in SAA is again the assumption of θ 1 (μ ≈ 1), but in addition, the angle between any two ray directions is approximated by the difference of their projections in the lateral plane. With SAA, the integral operator in RTE is reduced to a convolution operator and the RTE can be solved directly in frequency domain. For a uni-directional ocean source, McLean et al [14] obtained an expression for the BSF involving two double Fourier integrals, whose quadrature still requires significant numerical effort. Further simplifications of SAA are possible (e.g., [15]), but generally leads to appreciable loss in accuracy. In general, Swanson et al [16,17] show that SAA is valid only for small optical depths. More recently, Markel [18] and Liemert and Kienle [19] developed a rigorous method for the 3D RTE Green function using the rotated-frame method. The resulting solution form is quite involved (containing multiple summations of Fourier integrals) and is computationally less practical.
We present a highly-robust and efficient analytical solution for the 3D BSF using a new azimuthal single scattering approximation to reduce the original 3D RTE to a standard 2D RTE, and then combine the convolution operator idea of SAA and the rotated-frame method to effectively solve for the (exact) 2D RTE. The azimuthal single scattering approximation is based on the assumption that the azimuthal dependence of the radiance is highly-peaked in the forward direction, which typically obtains in ocean optics. The 2D BSF solution [see Eq. (37)] contains only two summations, while the 3D BSF is related to the 2D BSF via a simple relationship [see Eq. (44)]. Comparing our analytical BSF results to Monte Carlo simulation of the original RTE [20], we confirm that the 2D BSF is exact with excellent comparisons between the two for all parameter ranges. For the 3D BSF, the comparisons are very satisfactory for single scattering albedo values typical in the ocean. As expected, deviations occur only for large optical depths (> O(5)) at lateral positions close to the beam, and for very large single scattering albedo (up to 0.9).
Due to the relative simple form(s) of the 2D and 3D analytical solution(s), the present BSF is very amenable to numerical evaluation, and more efficient than direct Monte Carlo simulation especially for single scattering albedo that is not small. This approach provides great potentials for practical applications of both forward and inverse ocean optics problems, particularly in areas of optical oceanography includes predicting the LIDAR returns and the light scattering in the ocean turbulence. The method is espeecially suitable to deal with light scattering by turbulence of a turbid ocean, which is usually obtained by solving 3D RTE with temporally/spatially varying inherent optical properties. Considering the fact that the BSF is equivalent to the Green's function of RTE, one can immediately obtain the scattered light field by the turbid ocean turbulence using the first-order approximation and the present analytic approach.

Variables reduction of three-dimensional radiative transfer equation
In general, the beam spread function (BSF) is the spatial distributions of the light intensity in response to a narrow incident beam. In this work, we define the BSF as the scalar irradiance which is the function of the lateral spatial variables (x, y) and depth variable z for an uni-directional point source inside a homogenous three-dimensional horizontally-unbounded ocean column of large depth (Fig. 1). Here all position variables are normalized by the beam attenuation coefficient to represent the optical distances. We use E (3) (x, y, z) to represent the scalar irradiance, where for later reference (and in contrast to the two-dimensional case), the superscript denotes the irradiance is in three dimensions (3D). We assume the uni-directional light source is placed at the origin illuminates in positive z-direction. The scalar irradiance is axisymmetric with respect to the z-axis, and it is convenient to consider the problem in cylindrical coordinates The scalar irradiance is determined by the radiance L(r, z, θ , φ − φ r ); L is the function of both the spatial variables (r, φ r , z), and directional variables (θ , φ ); θ is the polar angle and φ is the azimuthal angle. Introducing a new azimuthal variable ψ = φ − φ r , the scalar irradiance can be defined as: The radiance L is a function of five independent variables and is the solution of the 3D radiative transfer equation (RTE): where δ (·) is the Dirac delta function; Φ(θ , ψ ; θ , ψ) is the phase function with arbitrary shape. So far, no assumption on the phase function and the solution has been made yet. The closedform solution of Eq. (2) is generally very difficult to obtain [18,19].
Now assumption can be made in the context of the typical ocean environment (e.g., case 1 water [6]) where light scattering tends to be highly forward peaked. The typical phase functions, such as the Petzold [21], and the Henyey-Greenstein (H-G) [22] with the asymmetric parameter g > 0.9, are highly peaked in the forward direction. Under the geometry of an uni-directional beam that emanates from a point source into the ocean, the scattering of light tends to spread more energy in the polar direction than in the azimuthal direction; the lateral spreading of a (narrow) incident beam is small relative to the forward direction. Therefore, a simplifying assumption is possible based on the highly peaked scattering property and is called 'azimuthal single scattering approximation'. In this approximation, we assume that in the azimuthal direction the scattering is relatively less important than in the polar direction, and therefore only single scattering occurs in the azimuthal direction. This assumption can be expressed mathematically by deliberately separating out azimuthal angle ψ in the radiance L. The dependence of L on azimuthal angle ψ is thus nontrivial only for |ψ| 1, especially for locations not close to the incident beam. According to the characteristics of single scattering [23], the radiance L now can be written as: Figure 1 demonstrates the geometry of the azimuthal single scattering approximation. It is seen that under such approximation, scattering of light is assumed to take place within the plane OABO r which is identical to the plane ψ = 0 and π. For forward-peaked scattering, the radiance is small except near θ = 0, and we write sin(θ ) ≈ sin(θ − θ ) and cos(θ ) ≈ cos(θ − θ ). We replace L in Eq. (2) with Eq. (3) and perform the integration over ψ to obtain: Note that the 3D phase function Φ in Eq. (2) is now replaced by a 2D phase function p(θ − θ ). The function p can be expressed explicitly by expanding Φ in spherical harmonics, and upon integration over ψ, to yield: where P n (x) is the Legendre polynomial and the coefficient g n is given by: By defining the azimuthal-averaged radiance I(r, z, θ ): the Eq. (4) can be written in two-dimensional form: We change the variable r to x and extending the original range of r from (0, ∞) to (−∞, ∞); at the same time, we extend range of the polar angle θ from (0, π) and [−π, π]. We are able to obtain the standard two-dimensional RTE with uni-directional point source: (9) We define scalar irradiance for the two-dimensional RTE: The relation between the original 3D scalar irradiance E (3) (r, z) and the 2D E (2) (x, z) is simply: Note that Eq. (11) is new and is an immediate result of the single scattering approximation in azimuthal direction we introduce. Thus, for conditions such as typical ocean where the scattering is highly focussed in the forward direction, the 3D BSF can be obtained directly and explicitly from the 2D BSF.

Analytical solution of two-dimensional scalar irradiance
We solve the standard 2D RTE (9) by using Fourier spectral method [24] in combination with the rotated frame method [18,19]. We expand the 2D radiance I(x, z, θ ) in terms of Fourier integrals: where the coefficients are:Î k x =k cos θ k , k z =k sin θ k . Equation (9) becomes: We define a new variable α = θ − θ k and write the radiance aŝ Similarly, the phase function p can be written as The source function in the spectral domain becomes Substituting I and p into Eq. (14), we obtain Truncation at some finite number of modes −N < m < N, we obtain a algebraic linear system for the unknown modal amplitudes: Accordingly, A and B are symmetric tri-diagonal and diagonal matrices, respectively, with elements: A n,m = δ n,m+1 + δ n+1,m , and B n,m = β n δ n,m , where δ n,m is the Kronecker delta function and β n =1 − ω 0 p n . To solveX analytically, we symmetrize the matrix A by defining We obtain, upon formal inversion:X where I is the identity matrix and D is a symmetric matrix with elements: D is Hermitian and can be decomposed as where the eigenmatrix V has real elements V nm , and Σ is diagonal with (real) eigenvalue λ m . Solving Eq. (27), we obtain the solutionX: Applying Eq. (10), we have The 2D scalar irradiance in physical space E (2) (x, z) can be obtained after inverse Fourier transform in cylindrical coordinates (ρ, θ ρ ): z=ρ cos θ ρ , x=ρ sin θ ρ . We obtain: To perform the θ k integral in Eq. (33), we use the identity where J n (x) is the bessel function of the first kind, to obatin To simplify this result, we invoke the properties of the eigenvalues and eigenmatrix of D: Considering the real part of the results, we obtain the exact 2D scalar irradiance expressed in the form: where K n (x) is the modified Bessel function of the second kind and (39) is the hypergeometric function and I n (x) is the modified Bessel function of the first kind. Following the same approach, we can derive the series form of the solution for the ballistic source. The only differences between the complete solution and ballistic solution are the eigenvalues and eigenvectors. The eigenvalues λ b m and eigenvectors V b nm for ballistic solution can be obtained analytically: and Therefore, the final solution of the scattered light fields for 2D BSF from an unidirectional beam source is E Change the variable of (ρ, θ ρ ) to (x, z), and we obtain the expression of the beam spread function obtained from a standard two-dimensional radiative transfer.
It is worth mentioning that H R n (ρ, λ ) and H I n (ρ, λ ) are solutions after inverse Fourier transforms for angular variables. Though they have explicit forms, they must be treated carefully. The two functions are both substractions of two terms. The two functions are non-singular but both terms in the function are singular for certain n, e.g., when n is even both 1 F 2 and csc(nπ/2) in H R n are singular; whilst n is odd 1 F 2 and sec(nπ/2) in H I n are singular. The accurate approach to evaluate these two functions is to expand the singular terms at n; the singular terms are crossed out and the non-singular terms in the expansion can be calculated easily. A less accurate but simpler approach is to take the limit of n numerically using the fact that both integrals are continuous with n.
To summarize the method to solve the two-dimensional standard RTE analytically, the radiance is expanded to Fourier series for both spatial and angular variables. The eigenvalues and eigenvectors must be calculated and stored in advance. The scalar irradiance as the function of depth and lateral position is obtained by the summation of terms containing the eigenvalues and elements of eigenvectors.
The solution of 3D beam spread function of light radiative transfer in an infinite ocean is obtained based on the solution of 2D RTE through Eq. (11) and (43); therefore, the 3D BSF for the scattered light fields can be expressed as:

Convergence evaluation
The validation of the present theory is made by comparisons with Monte Carlo radiative transfer simulation in two and three dimensions for arbitrary scattering phase function [25,26]. The codes are parallelized for high-performance computing platforms. All the results of BSF obtained by the Monte Carlo simulation for comparisons are performed by launching O(10 11 ) photons. Before validation of the (exact) 2D BSF solution [Eq. (37)], we first evaluate the convergence of Eq. (37) with increasing Fourier modes N. For simplicity, we utilize the 2D Henyey-Greenstein (H-G) phase function [22] in the evaluations: where Φ = θ − θ is scattering angle and g is the asymmetric parameter. The dependence of the solution on N is investigated with N varying from 40 to 800. The relative error (RE) of the scattered scalar irradiance is defined as RE = |E (2) s,N − E (2) s,∞ | |E (2) where E (2) s,∞ is the best estimate of the true value using Richardson extrapolation on the calculated data. Figure 2 shows RE for three different single scattering albedos ω 0 at different horizontal positions x (at a given depth). Exponential convergence RE∼ e −κN is clearly obtained. For given x, κ is independent of ω 0 (although RE is somewhat greater for smaller ω 0 ). With increasing x (at a given depth), κ is also somewhat greater.
It is worth noting that the convergence rate is related to the delta function angular source of the BSF. In many practical applications, the source angular spread is not infinitesimal, and one might approximate the delta function with a generalized delta function, say, of the form: As ε →0, the original delta function BSF is recovered. For small ε > 0, the converge of Eq.
(37) is improved. This is illustrated in Fig. 3 comparing the delta function 2D BSF with that using Eq. (37), for different ε. The deviation between the finite ε and ε=0 results are primarily for smaller x, as expected. As ε is decreased, the result approaches the ε=0 result smoothly; while the number of Fourier modes N required also increases.

Validation of 2D analytical BSF
We now validate the present 2D theory for ocean radiative transfer BSF by comparing the closed-form analytical solution with Monte Carlo simulations. We assume a homogeneous ocean with phase function as given earlier [Eq.(45)]. The only two factors that affect the lateral profiles of irradiance are the single scattering albedo ω 0 and optical depth z. Figure 4 shows the comparisons for broad ranges of ω 0 and z. For fixed optical depth z, Fig. 4(a), shows excellent agreement between analytical and simulation results for different ω 0 over all values of x. Similarly, for fixed ω 0 and varying z [ Fig. 4(b)] the comparison is excellent.

Validation of 3D analytical BSF
We now compare the 3D BSF given by the (approximate) analytical form [Eq. (44)] against (3D) Monte Carlo result. We again assume homogeneous ocean with 3D H-G phase function [22]: where Θ is the scattering angle. Figure 5 shows the comparisons for different single scattering albedos ω 0 and optical depths z. The comparisons are overall very satisfactory, showing the general validity of the azimuthal single scattering approximation we propose. For smaller optical depth, the results are especially good (for all lateral distance r) where multiple azimuthal scattering is expected to be limited. For large z, the deviations are primarily for small r, where multiple azimuthal scattering is expected to be more important. Comparing the results for different ω 0 , the deviations are somewhat increased for increasing ω 0 (for larger z at small r), again as expected, since greater ω 0 is associated with enhanced multiple azimuthal scattering.

Conclusion
The main contribution of the present work is the development of the analytical 3D BSF for homogenous ocean environment using a new single scattering approximation for the azimuthal direction (only). This reduces the 3D beam spread problem to that in 2D. We obtain the exact analytical 2D BSF by combining the Fourier spectral approach and the rotated-frame method [18,19], and obtaining for the first time an explicit expression in terms of sums of Fourier modes expressed in special functions (without inverse Fourier integrals). The present solution (involving a double Fourier sum) is a substantial improvement over existing small-angle approximations (containing four nested Fourier integrals [14]) in terms of computational simplicity and efficacy.
We compare the new 2D and 3D analytical BSFs against direct Monte Carlo simulations. The results are excellent over broad ranges of optical depths and single scattering albedos. The overall good comparisons for the 3D BSF support the general validity of the azimuthal single scattering approximation, with (small) deviations limited, as expected, to large optical depths and small lateral distances for extremely turbid water (single scattering albedos greater that ∼0.9).
The analytical BSFs are found to be computationally efficient with exponential convergence (error ∼ e −κN ) with increasing number of Fourier modes N (compared with ∼ N −1/2 convergence rate of standard Monte Carlo simulation using N realizations).
Our results obtained under the assumption of azimuthal single scattering should be valid in the ocean environment where scattering of light is highly forward-peaked. Comparisons and validation of these theoretical results against direct ocean measurements would be very desirable.