1. Introduction

Knowledge of light scattering by particles is inevitably necessary in numerous disciplines including atmospheric and oceanic radiative transfer, remote sensing, climate science, astrophysics, planetary sciences, and bio-optics. In the case of spherical particles, the rigorous Lorenz-Mie theory (commonly known as the Mie theory) was well established over 100 years ago [1,2], and is regarded as a major milestone in physics and optics. In this theory, the method of separation of variables in spherical coordinates is used to expand the incident and scattered electromagnetic waves in vector spherical harmonics.

Many particles in nature have nonspherical shapes. The spheroidal shape model can largely capture the nonspherical effects, even of many very complex natural shapes, on the optical properties. In the case of a sphere, only one degree of freedom, namely, particle size, is needed to determine the particle morphology. In the case of a spheroid, two degrees of freedom, namely, particle maximum dimension and aspect ratio (the ratio of the major-axis length to the minor-axis length), are needed to determine the particle morphology. The overall shapes of many natural particles such as bacteria and microweeds are approximately spheroidal [3]. For practical applications, dust aerosol particles can be approximated as spheroids in simulating their optical properties. For example, Dubovik et al. [4] shows that simulation based on the spherical model substantially deviates from the measured angular distribution of scattered light, but the spheroidal model provides a quite realistic result compared to measurements. The particle nonspherical effect on the optical properties of dust aerosol has a significant impact on remote sensing of dust aerosol properties and assessment of aerosol radiative forcing involved in climate study. Yang et al. [5] simulates the spectrum at the top of the atmosphere under dusty conditions by assuming the dust aerosol particles are spherical or spheroidal. As dust aerosol particles are exclusively nonspherical, the simulation results suggest that assuming these particles to be “equivalent” spheres would cause significant errors or even misleading results in remote sensing and radiative transfer simulations involved in climate study.

The extended boundary condition T-matrix method (EBCM) [6,7], invariant imbedding T-matrix method (IITM) [8,9] and several geometric optics methods [10–14] are commonly used for light-scattering computations involving spheroidal particles. EBCM is applicable only to spheroidal particles with relatively small size parameters (defined as major semiaxis multiplied by modified wavenumber) and is not capable of computing the optical properties of coarse-mode aerosol particles in ultraviolet to near infrared wavelengths. IITM can be applied to particles with moderate size parameters but is prohibitively computer-intensive when the size parameter is large. In practical applications, it is necessary to consider spheroids with various aspect ratios to model the nonsphericity effect of realistic particles. The applicability of EBCM and IITM degrades if the aspect ratio of spheroids deviates substantially from unity. For instance, the EBCM can be applied to a case with a size parameter more than 100 and aspect ratio 0.5, but is only applicable up to size parameter 20 if the aspect ratio is 0.1 [7]. For large particles, the geometric optics method can be used, which, however, may have unknown errors at some scattering angles, and the errors cannot be quantified because a benchmark reference is not available. In principle, the accuracy of the geometric optics method increases with increasing particle size. However, in the case of a sphere, the geometric optics method fails to simulate the backscattering peak corresponding to an optical phenomenon referred to as the “glory” even in the case of a large water drop, as compared to the Lorenz-Mie theory. Thus, we conjecture that the geometric optics method is not accurate in some cases for spheroids no matter how large the scattering particles are. Note that the physical-geometric optics method has been applied to accurately solve light scattering by large particles with planar facets ([15]; and the relevant references cited therein).

In the history of light scattering research, a Lorenz-Mie type solution was sought for the scattering of light by a spheroid without a limit on particle size and aspect ratio. In this line, the most pronounced achievement is the milestone work by Asano and Yamamoto [16], who published a theoretical solution for light scattering by a spheroid in terms of vector spheroidal harmonics. However, their computational program based on knowledge of spheroidal functions (and the relevant numerical schemes), and computer coding techniques at that time, was limited to small-size-parameter particles (up to about 30). Although other researchers (e.g., [17]) tried to improve the Asano and Yamamoto work, there has been very little progress in numerical capabilities for practical applications.

We significantly expand Asano-Yamamoto’s work. In particular, we develop an innovative separation of variables method (SVM) in spheroidal coordinates in conjunction with using the Wigner-d functions instead of associated Legendre polynomials to expand the angular spheroidal functions to solve the light scattering problem by spheroids. In numerical computation, we adopt an advanced numerical computational capability by van Buren et al. [18–20] to compute the radial spheroidal functions. The most pronounced achievement of this study is that the applicable size parameter of the present SVM can be more than 600 and the new method is not restricted by extreme aspect ratios. The description of the current SVM is detailed in Section 2. We show some computational examples using SVM in Section 3, and Section 4 summarizes this study.

2. Methodology

Spheroidal shapes include prolate and oblate spheroids. In Fig. 1(a), the spheroids are in the particle coordinate (o-x′y′z′) system where the z′ axis is the symmetry rotational axis. We define the aspect ratio as $a/c$, where a is the semiaxis length on the x′oy′ plane, and c is the semiaxis length along the z′ axis. $a/c \,<\, 1$ and $a/c \,>\, 1$ correspond to prolate and oblate spheroids, respectively.

 

Fig. 1. Illustration of spheroidal shapes and coordinates. (a) Prolate and oblate spheroids are in the particle coordinate system (o-x′y′z′). a is the semi-axis length on the x′oy′ plane, and c is the semi-axis length along the z′ axis. The aspect ratio is defined as $a/c$. (b) Oblate spheroidal coordinates ($\eta ,\; \xi $) on the Cartesian coordinate plane yoz. F₁ and F₂ are the two foci. d is the semifocal distance. (c) Laboratory (o-xyz) and particle (o-x′y′z′) coordinates. $\alpha $ and $\beta $ are Euler angles. ${\hat{e}^i}$ is the direction of incidence along the z axis. ${\hat{e}^s}$ is the scattering direction. $\theta $ and $\varphi $ are scattering zenith and azimuth angles with respect to o-xyz. $\theta ^{\prime}$ and $\varphi ^{\prime}$ are scattering zenith and azimuth angles with respect to o-x′y′z′.

Download Full Size | PDF

Light scattering by a sphere is analytically formulated using spherical coordinates in the Lorenz-Mie theory. Analogously, light scattering by a spheroid can be analytically solved in spheroidal coordinates, as illustrated in Fig. 1(b). If we assume a Cartesian coordinate system $({x,\; y,\; z} )$ shares its origin with a spheroidal coordinate system $({\xi ,\; \eta ,\; \varphi } )$, the relation between the two coordinate systems is expressed as [21]

Orientation of an arbitrary particle relative to the laboratory coordinate system can be characterized by three Euler angles, $\alpha $, $\beta $ and $\gamma $, which describe the rotation of the particle coordinate relative to a fixed laboratory coordinate [22]. When $\alpha = \beta = \gamma = 0$, the particle and laboratory coordinates coincide. Because spheroids have rotational symmetry, only angles $\alpha $ and $\beta $ are needed, as defined in Fig. 1(c). We define the angle of incidence $\zeta $ as the angle between the spheroid’s symmetry rotational axis and the direction of incidence. According to Fig. 1(c), $\zeta $ is the angle between the laboratory coordinate z axis and the particle coordinate (z′) axis. Therefore, we have $\zeta = \beta $.

2.1 Separation of variable method in spheroidal coordinates

Asano and Yamamoto [16], and Voshchinnikov and Farafonov [17] propose different SVM formulations for light scattering in the spheroidal coordinates. In particular, in [16], the electromagnetic fields are expanded in terms of vector spheroidal functions (VSFs). The formulation in [17] expresses the electric and magnetic fields as the derivatives of scalar potentials, and the scalar potentials can be expanded in terms of spheroidal functions. In both formulations, the expansion coefficients of the scattered field are derived by solving a linear equation system. The previous computational capabilities based on these two formulations have similar applicable size parameter and aspect ratio ranges. This study is based on expanding the formulation by [16], which is most similar to the SVM formulation of the Lorenz-Mie theory for light scattering by spheres in the spherical coordinate system.

In a light scattering event, the incident, scattered and internal electromagnetic fields can be expanded in terms of VSFs, following [16], in the form

The VSFs $\mathbf{M}$ and $\mathbf{N}$ can be derived from the scalar spheroidal function $\psi $ [16],

In spheroidal coordinates, $\psi = \psi ({\eta ,\xi ,\varphi ,d} )$, and can be written as the product of three equations, $\psi \equiv R({\xi ,d} )S({\eta ,d} )\mathrm{\Phi }(\varphi )$. Equation (5) can thus be separated into three ordinary differential equations [21],

On the boundary of a spheroid, the components of the incident, scattered and internal electromagnetic fields must satisfy the following boundary conditions [16]:

Utilizing the orthogonality of the terms in VSFs with different degrees and orders, and Eq. (7), the expansion coefficients of the scattered and internal fields can be obtained by solving a linear equation system involving known parameters ${f_{mn}}$ and ${g_{mn}}$, and terms with angular and radial spheroidal functions.

The expansion coefficients of the scattered field contain the optical property information of the scattering particle. The far-field scattering properties including extinction (${C_{\textrm{ext}}}$) and scattering (${C_{\textrm{sca}}}$) cross sections, and amplitude scattering matrix elements, can be computed using the expansion coefficients. As an example, the amplitude scattering matrix $\mathbf{S}$ is defined as

Detailed formulations of the linear equation system and expressions for the extinction and scattering cross sections can be found in [16].

2.2 Spheroidal functions

The above description of the SVM suggests that accurate computation of the angular and radial spheroidal functions is a prerequisite to obtain the expansion coefficients of the scattered field.

Conventionally, ${S_{mn}}$ can be expanded in terms of associated Legendre polynomials $P_n^m$ [21],

Here we adopt the Hodge [24] method to calculate the approximate values of ${\lambda _{mn}}$ for $n = m,\; m + 2,\; m + 4,\; \ldots $. In [24], the calculation of ${\lambda _{mn}}$ is converted to a problem of calculating the eigenvalues of a symmetric tridiagonal matrix $\lambda _{mn}^{\prime}$. Because the matrix has a finite dimension, $\lambda _{mn}^{\prime}$ values are not identical to ${\lambda _{mn}}$. We then utilize the Bouwkamp [25] method to correct $\lambda _{mn}^{\prime}$ and obtain accurate ${\lambda _{mn}}$ values. Once an accurate ${\lambda _{mn}}$ is attained, the expansion coefficients $f_r^{mn}$ at all possible r can be calculated using the recurrence relations.

With the expansion Eq. (11), the SVM formulation has many series summation operations involving $f_r^{mn}$ . In Fig. 2, we show the absolute value of $f_r^{mn}$ based on the expansion in terms of $P_{m + r}^m$, which has substantial variations at different expansion orders r. The large difference among the expansion coefficients may cause errors in summing the series in the numerical computation due to overflow and underflow problems. Alternatively, we consider the expansion of ${S_{mn}}$ in terms of Wigner-d functions,

 

Fig. 2. The absolute values of the expansion coefficients ($f_r^{mn}$) of ${S_{mn}}$ in terms of associated Legendre polynomials and Wigner-d functions, with $m = 100$, $n = 200$ and $kd = 500$.

Download Full Size | PDF

Wigner-d functions are used in T-matrix methods to replace the associated Legendre polynomials [9,26], and are found to make the T-matrix algorithm more numerically stable [26]. Figure 2 compares the absolute values of $f_r^{mn}$ in the expansion in terms of $d_{m0}^n$ and $P_n^m$. Obviously, the variation of $f_r^{mn}$ in terms of $d_{m0}^n$ is much smaller than that in terms of $P_n^m$. With the ${S_{mn}}$ expansion Eq. (12), the numerical implementation of the SVM can avoid overflow and underflow problems to some extent.

${S_{mn}}$ is normalized as

Because we use Wigner-d functions rather than associated Legendre polynomials to expand ${S_{mn}}$, the expressions of equations in the SVM formulation involving the expansion coefficients are different from those in [16]. For example, an angular function ${A_{mn}}(\eta )$ [16] with $m \ge 1$ in the SVM formulation is defined as

Other angular functions can be derived in a similar manner based on the orthogonality of Wigner-d functions. The expressions of the angular functions and their expansion coefficients for $m \ge 1$ are given in Supplement 1. The expressions of the expansion coefficients for $m = 0$ are identical to those in [16].

Flammer [21] describes several formulas to compute radial spheroidal functions. The computation of first- and second-kind radial spheroidal functions $R_{mn}^{(1 )}$ and $R_{mn}^{(2 )}$ involves series summation, in which each term contains the product of $f_r^{mn}$ and a Bessel function for $R_{mn}^{(1 )}$ or a Neumann function for $R_{mn}^{(2 )}$. In the SVM formulation, we need $R_{mn}^{(1 )}$ and the third-kind radial spheroidal functions, $R_{mn}^{(3 )} = R_{mn}^{(1 )} + iR_{mn}^{(2 )}$, and their first derivatives with respect to $\xi $. The series summation may diverge when $kd$ and $n - m$ values are large, especially for $R_{mn}^{(2 )}$ and $dR_{mn}^{(2 )}/d\xi $. In our numerical implementation of the SVM, we utilize the approach and Fortran subroutines developed by van Buren et al. [18–20] (subroutine programs are available on https://github.com/MathieuandSpheroidalWaveFunctions) to compute $R_{mn}^{(1 )}$ and $R_{mn}^{(2 )}$, and their first derivatives. The van Buren et al. method is accurate and numerically stable for $kd$ up to several thousand. The Wronskian of the radial spheroidal function has a simple expression [21],

2.3 Random orientation-averaging

The SVM optical property computational result is for a spheroid with fixed orientation. In practical applications such as remote sensing of a plume of aerosols, we may need random orientation-averaged single-scattering properties, including the extinction efficiency (${\textrm{Q}_e}$), single-scattering albedo (SSA), asymmetry factor (g), and scattering phase matrix (P).

Random orientation-averaged ${\textrm{Q}_e}$ and SSA of a spheroid are computed as

Note that the SVM computation is conducted in the particle coordinate system. For an orientation with Euler angle $\beta $, the incident electric field in the particle coordinate system is expressed as

The integrals in Eqs. (19)–(21) are evaluated numerically. In particular, Gaussian quadrature is used to compute the integrals with respect to $\beta $. The numbers of $\beta $ and $\alpha $ are large enough to guarantee converged results.

The SVM formulation described above includes infinite series, such as the expansions of the angular spheroidal function and electromagnetic fields. Every infinite series must be truncated at a finite order, so the SVM numerical implementation truncates the expansion of angular spheroidal functions, the expansion of the electromagnetic field in Eqs. (2) and (3), and the expansion coefficients of the scattered field. While the exact results at infinite order are unknown, and the truncation criteria below are preliminary, results in Section 3 show that single-scattering properties computed by SVM are very similar to the results of other well-tested methods.

In our computations, summation of $f_r^{mn}({kd} )$ stops when the absolute value reaches a very small number related to the computer precision. For the expansion of the electromagnetic field, we define the minimum expansion degree n to obtain the converged result as N. Our numerical experiment indicates that N has positive correlations with the size parameter, aspect ratio extremality and real part of refractive index. The convergence is determined by the criterion [26]

3. Example results

In this section, we show random orientation-averaged single-scattering properties computed by the SVM, and comparisons with EBCM and IITM, and a separate comparison with the improved geometric optics method (IGOM; [11,14]).

In the computation of the following results, we express the particle size parameter as k multiplied by the major semiaxis. The size parameters of prolate and oblate spheroids are defined as kc and ka, respectively. Two refractive indices $m = 1.31$ and $m = 1.53 + i0.001$ are considered in the computation. Refractive index $1.31$ is for water ice at visible wavelengths, and $1.53 + i0.001$ is for mineral dust at visible to near infrared wavelengths. The refractive index of the surrounding medium is assumed to be 1.

Figure 3 shows ${\textrm{Q}_e}$ and g of nonabsorptive spheroids computed by EBCM, IITM and SVM. The aspect ratios of the prolate and oblate spheroids are 0.5 and 2.0. EBCM can be applied to size parameters up to 117 and 111 for prolate and oblate spheroids. IITM can be applied to a size parameter up to 500 for both shapes. The applicable size parameter range of the SVM is more than 600 for both shapes. At the size parameters where all three methods are applicable, ${\textrm{Q}_e}$ and g computed by the three methods are identical up to the fourth decimal digit.

 

Fig. 3. Extinction efficiency (${\textrm{Q}_e}$) and asymmetry factor (g) computed by EBCM, IITM and SVM. (a) and (c) are for a prolate spheroid with aspect ratio 0.5. (b) and (d) are for an oblate spheroid with aspect ratio 2. The refractive index is $m = 1.31$. The vertical dashed lines in each plot mark the maximum applicable size parameters of EBCM ($kc = 117$ for prolate and $ka = 111$ for oblate) and IITM ($kc = 500$ and $ka = 500$ for both prolate and oblate).

Download Full Size | PDF

Figure 4 shows ${\textrm{Q}_e}$, SSA and g of absorptive spheroids computed by EBCM, IITM and SVM. The prolate and oblate spheroids have extreme aspect ratios 0.1 and 10. The maximum applicable size parameters of EBCM and IITM become much smaller for cases with extreme aspect ratios. For the cases in Fig. 4, EBCM is applicable to size parameters up to 20 and 21 for prolate and oblate shapes. IITM is applicable up to size parameter 83. The SVM is not restricted by extreme aspect ratios, and is still applicable for size parameters more than 500.

 

Fig. 4. Extinction efficiency (${\textrm{Q}_e}$), single-scattering albedo (SSA) and asymmetry factor (g) computed by EBCM, IITM and SVM. (a), (c) and (e) are for a prolate spheroid with aspect ratio 0.1. (b), (d) and (f) are for an oblate spheroid with aspect ratio 10. The refractive index is $m = 1.53 + i0.001$. The vertical dashed lines in each plot mark the maximum applicable size parameters of EBCM ($kc = 21$ for prolate and $ka = 20$ for oblate) and IITM ($kc = 83$ and $ka = 83$ for both prolate and oblate).

Download Full Size | PDF

Figure 5 compares the scattering phase matrix elements computed by EBCM, IITM and SVM for an oblate spheroid with size parameter 111, which is the largest size parameter where all the three methods are applicable with aspect ratio 2 and refractive index 1.31. The results by the three methods are consistent with noticeable slight differences in the ${P_{12}}$ and ${P_{43}}$ elements. At scattering angles from 150 to 180°, the scattering phase matrix elements have strong oscillations, and the three methods are still consistent.

 

Fig. 5. Scattering phase matrix elements computed by EBCM, IITM and SVM for an oblate spheroid with size parameter (ka) 111, aspect ratio 2 and refractive index 1.31. (a) ${P_{11}}$; (b) ${P_{12}}/{P_{11}}$; (c) ${P_{22}}/{P_{11}}$; (d) ${P_{33}}/{P_{11}}$; (e) ${P_{43}}/{P_{11}}$; (f) ${P_{44}}/{P_{11}}$. The inset plots in (a), (c)-(f) show the values at scattering angles from 150 to 180°.

Download Full Size | PDF

Figure 6 compares the scattering phase matrix elements computed by EBCM, IITM and SVM for a prolate spheroid with size parameter 21, which is the largest size parameter where all three methods are applicable with aspect ratio 0.1 and refractive index 1.53 + i0.001. The results by the three methods are identical at most of the scattering angles. ${P_{12}}$ and ${P_{43}}$ elements by IITM are slightly different from those by EBCM and SVM.

 

Fig. 6. The same as Fig. 5, but for a prolate spheroid with size parameter (kc) 21, aspect ratio 0.1 and refractive index 1.53 + i0.001.

Download Full Size | PDF

In Figs. 5 and 6, the computations by all three methods converge with respect to their numerical computational parameters and procedures, such as the dimension of the T-matrix in EBCM and IITM, and the expansion degree truncation and integration for random orientation averaging in the SVM. However, some results do not converge to exactly the same values, containing slight differences in the scattering phase matrices. Although EBCM and IITM are considered to be numerically exact methods, their detailed algorithms differ. The numerical implementation of the SVM can also incur numerical errors such as in solving the linear equation system and computing the series Eq. (9). Therefore, the results by the three methods may have different numerical errors, which are manifested in their scattering phase matrix results.

Figures 7 and 8 compare scattering phase matrix elements computed by SVM and IGOM with large size parameters. In Fig. 7, the spheroid has a size parameter of 600 and is nonabsorptive. Although IGOM uses the geometric optics approximation, its result is quite consistent with SVM for such a large size parameter. At scattering angles close to 180°, the IGOM result does not show a backscattering peak, which is obvious in the SVM result. In Fig. 8, the spheroid has a size parameter of 500, aspect ratio 10, and some absorption. The IGOM result is still largely consistent with the SVM result.

 

Fig. 7. Scattering phase matrix elements computed by SVM and IGOM for a prolate spheroid with size parameter (kc) 600, aspect ratio 0.5 and refractive index 1.31.

Download Full Size | PDF

 

Fig. 8. The same as Fig. 7, but for an oblate spheroid with size parameter (ka) 500, aspect ratio 10 and refractive index 1.53 + i0.001.

Download Full Size | PDF

In SVM numerical computations, for aspect ratios that are not much greater or smaller than unity, such as 2 and 0.5, we only use double-precision floating-point real variables (64 bits) and can achieve accurate computation for size parameters up to 600. For extreme aspect ratios, such as 10 and 0.1, it is necessary to use extended double-precision floating-point real variables (128 bits) to compute radial spheroidal functions for large size parameters. The SVM developed in this study is numerically stable. CPU time is mainly spent in the computation of the spheroidal functions, solving the linear equations and calculating the amplitude scattering matrix. Storage of the constants in the linear equations uses most of the memory space. We use the Open Multi-Processing (OpenMP) technique to parallelize the computer program written in Fortran.

With one 48-core node on the Grace cluster in Texas A&M University, computing the SVM result in Fig. 7 for a spheroid with size parameter 600 and aspect ratio 0.5 takes only a few minutes for one particle orientation at a single scattering angle, and in total about 4 hours for the random orientation-averaged result at a single scattering angle. The computation of the IGOM result in Fig. 7 only needs about 8 seconds with one core. For large size parameters, an economical approach may be to use IGOM to compute scattering phase matrices excluding backscattering, and to use SVM to compute scattering phase matrices at backscattering directions, and extinction efficiency and single-scattering albedo.

In Fig. 5 for an oblate spheroid with size parameter 111 and aspect ratio 2.0, the EBCM computation takes 717 seconds with a single core. The IITM computation takes 404 seconds with 96 cores. With 48 cores, the SVM computation takes 18 seconds for one particle orientation at a single scattering angle, and 115 seconds for the random orientation-averaged result at a single scattering angle. In Fig. 6 for a prolate spheroid with size parameter 21 and aspect ratio 0.1, the EBCM computation takes 6 seconds with a single core. The IITM computation takes 19 seconds with 96 cores. With a single core, the SVM computation takes 42 seconds for the random orientation-averaged result at a single scattering angle.

To compute random orientation-averaged single-scattering properties, EBCM and IITM only need to compute the T-matrix for a single particle orientation. Then, the random orientation-averaged single-scattering properties can be computed analytically with this T-matrix [26]. The SVM must solve linear equation systems for each Euler angle $\beta $. In addition, a majority of total computational time is spent in computing the scattering phase matrix using the expansion coefficients of the scattered field. In the applicable size parameter and aspect ratio ranges of all three methods, the EBCM is the most efficient for randomly oriented particles, whereas the SVM is the most efficient for a single orientation. The SVM can also be used to compute the T-matrix of a particle [29], and computation of the random orientation-averaged single-scattering properties would be much more efficient using SVM.

4. Summary

We develop a new Lorenz-Mie theory-type SVM in spheroidal coordinates in conjunction with using the Wigner-d function to expand the angular spheroidal functions. In numerical implementation of the present SVM, we utilize an advanced radial spheroidal function computational algorithm. The present SVM has much larger applicable size parameter and aspect ratio ranges than previous spheroidal coordinate-based SVMs. The computational results by SVM are consistent with those by two T-matrix methods, namely EBCM and IITM. The SVM is applicable to greater size parameters than the T-matrix methods, and is not restricted when the aspect ratio deviates much from unity. The geometric optics method is found to be consistent with the SVM for spheroids with large size parameters.