Improvements in the computational efficiency and convergence of the Invariant Imbedding T-matrix method for spheroids and hexagonal prisms

Siyao Zhai; R. Lee Panetta; R. Lee Panetta; Ping Yang; Ping Yang

doi:10.1364/OE.27.0A1441

1. Introduction

For atmospheric and oceanic radiative transfer simulations and remote sensing implementations, accurate computation of the single-scattering properties, namely the extinction cross section, the single-scattering albedo, and the phase matrix, of non-spherical particles is essential [1,2]. Because of the linearity of the equations governing the light-scattering process, the transformation from the incident field to the scattered field through their expansions in terms of vector spherical wave functions can be described by the Transition Matrix (T-Matrix) [3,4]. The single-scattering properties of the scattering particle can be derived from the T-matrix. The T-Matrix is an inherent property of the particle; in particular, it depends only on the particle structure, size and refractive index. The T-matrix of a homogeneous sphere reduces to a diagonal form with elements on the diagonal corresponding to the Lorenz-Mie series expansion coefficients. For nonspherical particles, the T-matrix method is computationally efficient because it is feasible to analytically average the optical properties over random orientations, as illustrated in [5,6]. Several efficient computational methods have been developed to compute the T-matrix of nonspherical and inhomogeneous particles [7,8]. The invariant-imbedding T-matrix (II-TM) method is a robust method developed to compute the single-scattering properties of non-spherical and inhomogeneous dielectric particles [9–11].

The II-TM method computes the T-Matrix of a scattering particle in a radial recurrence form starting with a certain initial inscribed sphere inside the particle and ending with the circumscribed sphere of the particle. At each recurrence level, Gauss-Legendre (GL) quadrature and matrix inversions are performed. As a result, three parameters determine the accuracy of the aforementioned numerical procedures: the radial discretization $Δ r$ ; the number of GL quadrature nodes $N_{q}$ , and the T-matrix expansion order N. In the numerical implementation of II-TM, a quadrature scheme is applied to define the geometry of the scattering particle. Any numerical errors associated with the quadrature scheme may be magnified in the matrix inversions and radial recurrence.

The focus of this study is on improving the computational efficiency and convergence of the GL quadrature. For particles with sharp edges, discontinuities occur in the quadrature method; we modify a commonly used scheme to avoid performing quadrature over such discontinuities. The example of a particle with sharp edges we consider here is inherently of considerable interest in atmospheric single-scattering studies, but the modification we make can be more generally applied. In addition, we replace the GL node generating scheme with the fastest algorithm to date [12]. These improvements substantially enhance the computational efficiency of the II-TM method for spheroids and hexagonal columns.

2. The T-matrix method

The incident and scattered electric fields can be expanded with the vector spherical wave functions (VSWFs) ${\overset{⇀}{M}}_{m n}^{(1, 3)}$ and ${\overset{⇀}{N}}_{m n}^{(1, 3)}$ in the form [13]

{\overset{⇀}{E}}_{i n c} (\overset{⇀}{r}) = E_{0} \sum_{n = 0}^{\infty} \sum_{m = - n}^{n} a_{m n} {\overset{⇀}{M}}_{m n}^{(1)} (k \overset{⇀}{r}) + b_{m n} {\overset{⇀}{N}}_{m n}^{(1)} (k \overset{⇀}{r}),

{\overset{⇀}{E}}_{s c a} (\overset{⇀}{r}) = E_{0} \sum_{n = 0}^{\infty} \sum_{m = - n}^{n} p_{m n} {\overset{⇀}{M}}_{m n}^{(3)} (k \overset{⇀}{r}) + q_{m n} {\overset{⇀}{N}}_{m n}^{(3)} (k \overset{⇀}{r}),

where k is the wavenumber,

\overset{⇀}{r}

is the position vector. Superscript 1 indicates regular wave functions in which spherical Bessel functions of the first kind are used to meet the boundary condition at coordinate center (particle center), and superscript 3 indicates radiating wave functions in which spherical Hankel functions are used to represent the field radiating from the source (particle center) to infinity. The T-matrix relates the expansion coefficients:

(\begin{matrix} p_{m n} \\ q_{m n} \end{matrix}) = E_{0} \sum_{n^{'} = 0}^{\infty} \sum_{m^{'} = - n^{'}}^{n^{'}} (\begin{matrix} T_{m n m^{'} n^{'}}^{11} & T_{m n m^{'} n^{'}}^{12} \\ T_{m n m^{'} n^{'}}^{21} & T_{m n m^{'} n^{'}}^{22} \end{matrix}) (\begin{matrix} a_{m^{'} n^{'}} \\ b_{m^{'} n^{'}} \end{matrix}) .

If index n is truncated at

n_{\max}

, the T-matrix is a

n_{\max} (n_{\max} + 2) \times n_{\max} (n_{\max} + 2)

square matrix. The two column vectors on the left and right sides of Eq. (3) are

n_{\max} (n_{\max} + 2)

column vectors. Denote the T-matrix as T and the column vectors as p and a. The above equation can be written in a concise super-matrix form:

p = T a .

After the T-matrix is obtained, it is straightforward to compute the single-scattering properties (see [14] for mathematical details).

3. Improving the convergence of Gaussian quadrature in II-TM

In the II-TM method, the T-matrix of a particle is obtained recursively, starting with the T-matrix $T (r_{0})$ of an inscribed sphere of the non-spherical particle, successively calculating T-matrices for a series of concentric spheres, and ending at the circumscribed sphere of the particle (r = R). The T-matrix of the inscribed sphere $T (r_{0})$ can be obtained with the Lorenz-Mie theory [9]. The T-matrix recurrence formula that relates the T-matrix $T (r_{0})$ for sphere of radius $r_{p}$ to the T-matrix $T (r_{p - 1})$ for the sphere of smaller radius $r_{p - 1}$ is given in super-matrix notation [9] as:

\begin{array}{l} T (r_{p}) = Q_{11} (r_{p}) + [I + Q_{12} (r_{p})] {[I - T (r_{p - 1}) Q_{22} (r_{p})]}^{- 1} \cdot \\ T (r_{p}) [I + Q_{21} (r_{p})] . \end{array}

_{$Q_{i j} (r_{p}), i, j = 1, 2$}are all

n_{\max} (n_{\max} + 2) \times n_{\max} (n_{\max} + 2)

square super-matrices given by

Q_{11} (r_{p}) = i k J^{T} (r_{p}) Q (r_{p}) J (r_{p}),

Q_{12} (r_{p}) = i k J^{T} (r_{p}) Q (r_{p}) H (r_{p}),

Q_{21} (r_{p}) = i k H^{T} (r_{p}) Q (r_{p}) J (r_{p}),

Q_{22} (r_{p}) = i k H^{T} (r_{p}) Q (r_{p}) H (r_{p}),

where

J (r_{p})

and

H (r_{p})

are

n_{\max} (n_{\max} + 2) \times n_{\max} (n_{\max} + 2)

diagonal square super-matrices whose diagonal elements are defined by

{\bar{\bar{J}}}_{n} (r_{p})

and

{\bar{\bar{H}}}_{n} (r_{p})

[9]:

{\bar{\bar{J}}}_{n} (r_{p}) = (\begin{matrix} j_{n} (k r_{p}) \\ 0 \\ 0 \end{matrix} \begin{matrix} 0 \\ \frac{1}{k r} \frac{\partial}{\partial r} r j_{n} (k r_{p}) \\ \sqrt{n (n + 1)} \frac{j_{n} (k r_{p})}{k r_{p}} \end{matrix}),

{\bar{\bar{H}}}_{n} (r_{p}) = (\begin{matrix} h_{n}^{(1)} (k r_{p}) \\ 0 \\ 0 \end{matrix} \begin{matrix} 0 \\ \frac{1}{k r} \frac{\partial}{\partial r} r h_{n}^{(1)} (k r_{p}) \\ \sqrt{n (n + 1)} \frac{h_{n}^{(1)} (k r_{p})}{k r_{p}} \end{matrix}),

where

j_{n}

is the spherical Bessel function of the first kind,

h_{n}^{(1)}

is the spherical Hankel function of the first kind. The super-matrix

Q (r_{p})

in Eqs. (6)-(9) is given by,

Q (r_{p}) = Δ r {[I - Δ r U (r_{p}) g (r_{p}, r_{p})]}^{- 1} U (r_{p}) .

In Eq. (12), the elements of

g (r_{p}, r_{p})

contain spherical Bessel and Hankel functions of the first kind, the elements of

U (r_{p})

contain zenith and azimuth integrals on the spherical layers that are used to define the particle geometry (see Eqs. (23)-(25) of [9]). Two inversions must be performed to obtain the T-matrix

T (r_{p})

in each iteration (Eqs. (5) and (12)). Equation (5) was first derived through discretizing integral equations of T-matrix [13]. Figure 1 shows the geometries. The origin of the coordinate system is at the particle center.

Fig. 1 Schematic of the radial recurrence. Particle in this schematic configuration is a particle with arbitrary geometry marked in green. Coordinate origin is on the geometric center of the particle. Two spherical shells are indicated by dashed circles. Arrows indicate the two radii $r_{p}$ and $r_{p + 1}$ .

Download Full Size | PDF

The matrix elements of $U (r_{p})$ for each iteration are surface integrals in zenith and azimuth directions on the spherical shell. The nonzero elements in the U-matrix can be written as [10]

U_{m n m^{'} n^{'}}^{i j} = A_{m n m^{'} n^{'}} \int_{0}^{π} d θ \sin θ F_{m m^{'}} (r, θ) K_{m n m^{'} n^{'}}^{i j} (θ), i, j = 1, 2,

U_{m n m^{'} n^{'}}^{33} = A^{'}_{m n m^{'} n^{'}} \int_{0}^{π} d θ \sin θ F^{'}_{m m^{'}} (r, θ) d_{0 m}^{n} (θ) d_{0 m^{'}}^{n^{'}} (θ),

where the integrations in the azimuth direction

F_{m m'} (r, θ)

and

F'_{m m'} (r, θ)

are

F_{m m^{'}} (r, θ) = \int_{0}^{2 π} d φ e^{- i (m - m^{'}) φ} [ε (r, θ, φ) - 1],

F^{'}_{m m^{'}} (r, θ) = \int_{0}^{2 π} d φ e^{- i (m - m^{'}) φ} \frac{[ε (r, θ, φ) - 1]}{ε (r, θ, φ)},

and the constants

A_{m n m' n'}

and

A'_{m n m' n'}

are

A_{m n m^{'} n^{'}} = {(k r)}^{2} {(- 1)}^{m + m^{'}} {[\frac{2 n + 1}{4 π n (n + 1)}]}^{1 / 2} {[\frac{2 n^{'} + 1}{4 π n^{'} (n^{'} + 1)}]}^{1 / 2},

A^{'}_{m n m^{'} n^{'}} = \sqrt{n n^{'} (n + 1) (n^{'} + 1)} A_{m n m^{'} n^{'}} .

The terms

K_{m n m' n'}^{i j} (θ)

appearing in the integrand of Eq. (13) are defined by

\begin{array}{l} K_{m n m^{'} n^{'}}^{i j} (θ) = \\ (\begin{matrix} π_{m n} (θ) π_{m^{'} n^{'}} (θ) + τ_{m n} (θ) τ_{m^{'} n^{'}} (θ) & - i [π_{m n} (θ) τ_{m^{'} n^{'}} (θ) + τ_{m n} (θ) π_{m^{'} n^{'}} (θ)] \\ i [π_{m n} (θ) τ_{m^{'} n^{'}} (θ) + τ_{m n} (θ) π_{m^{'} n^{'}} (θ)] & π_{m n} (θ) π_{m^{'} n^{'}} (θ) + τ_{m n} (θ) τ_{m^{'} n^{'}} (θ) \end{matrix}), \end{array}

where

π_{m n}

and

τ_{m n}

are obtained from the Wigner-d function

d_{0 m}^{n} (θ)

π_{m n} (θ) = \frac{m}{\sin θ} d_{0 m}^{n} (θ),

τ_{m n} (θ) = \frac{d}{d θ} d_{0 m}^{n} (θ) .

The U-matrix elements are different for different particles. Below we outline the formulations for the spheroid and hexagonal column. These shapes are used as basic particle structures in modelling atmospheric ice crystals and other aerosols. A spheroid has two important kinds of symmetries that allow computational simplifications: rotational symmetry with respect to the z-axis

ε (r, θ, φ) = ε (r, θ)

and mirror symmetry with respect to the xy-plane

ε (r, θ) = ε (r, π - θ)

. Also,

ε (r, θ) - 1 = {\begin{matrix} 0, θ_{2} < θ, \\ ε - 1, 0 < θ < θ_{2}, \end{matrix}

where

θ_{2}

is the angle at which the spherical shell intersects the spheroid (Fig. 2).

Fig. 2 Procedures in the sequence of expanding spherical shells in the recurrence. In the case of the spheroid, $θ_{1} = 0$ . $θ_{2}$ is the angle where the spherical shell intersects the spheroid.

Download Full Size | PDF

The following properties of the Wigner-d functions allow us to simplify the formulations for spheroid [15]

\begin{array}{l} [π_{m n} (θ) π_{m n^{'}} (θ) + τ_{m n} (θ) τ_{m n^{'}} (θ)] \sin θ \\ = \frac{d}{d θ} (τ_{m n'} (θ) d_{0 m}^{n} (θ) \sin θ) + n' (n' + 1) \sin θ d_{0 m}^{n} (θ) d_{0 m}^{n'} (θ) \\ = \frac{d}{d θ} (τ_{m n} (θ) d_{0 m}^{n'} (θ) \sin θ) + n (n + 1) \sin θ d_{0 m}^{n} (θ) d_{0 m}^{n'} (θ), \end{array}

[π_{m n} (θ) τ_{m n^{'}} (θ) + τ_{m n} (θ) π_{m n^{'}} (θ)] \sin θ = m \frac{d}{d θ} (d_{0 m}^{n} (θ) d_{0 m}^{n'} (θ)) .

Substituting Eqs. (23) and (24) into Eqs. (13) and (14), followed by integration by parts, we obtain the simplified formulations:

\begin{array}{l} U_{m n m' n'}^{11} = U_{m n m' n'}^{22} \\ = A_{m n m' n'} 2 π δ_{m m'} [ε - 1] \cdot \\ [c_{n n'} τ_{n n'} (θ_{2}) d_{0 m}^{n} (θ_{2}) \sin θ_{2} + n' (n' + 1) \int_{0}^{θ_{2}} d θ \sin θ d_{0 m}^{n} (θ) d_{0 m}^{n'} (θ)], \\ c_{n n'} = 1 + (^{- 1) n + n'}, \end{array}

U_{m n m' n'}^{33} = A'_{m n m' n'} 2 π δ_{m m'} \frac{ε - 1}{ε} c_{n n'} \int_{0}^{θ_{2}} d θ \sin θ d_{0 m}^{n} (θ) d_{0 m}^{n'} (θ),

\begin{array}{l} U_{m n m' n'}^{12} = - U_{m n m' n'}^{21} \\ = A_{m n m' n'} 2 π δ_{m m'} [ε - 1] c'_{n n'} (- i m) [d_{0 m}^{n} (θ_{2}) d_{0 m}^{n'} (θ_{2}) - d_{0 m}^{n} (0) d_{0 m}^{n'} (0)], \\ c'_{n n'} = 1 + (^{- 1) n + n' + 1}, \end{array}

The original integrals in the zenith direction in Eqs. (13) and (14) are reduced to

\int_{0}^{θ_{2}} d θ \sin θ d_{0 m}^{n} (θ) d_{0 m}^{n'} (θ)

in Eqs. (25) and (26). The hexagonal column also has a rotational symmetry about the z-axis, but it is no longer continuous and is instead discrete. As a consequence, the azimuth integrals

F_{m m'} (r, θ)

and

F'_{m m'} (r, θ)

cannot be decoupled from the zenith integrals. The values of the azimuth integrals are dependent on

θ

and r because the spherical shell intersects a hexagonal column in different ways as

θ

and r change. Figure 3 illustrates the three different situations as the recurrence progresses in the case of a long column.

Fig. 3 Procedures in the sequence of expanding spherical shells in the recurrence for a hexagonal column. $θ_{2}$ is the angle where the spherical shell intersects the column.

Download Full Size | PDF

Using the mirror symmetry of $ε (r, θ) = ε (r, π - θ)$ , we can simplify Eqs. (13) and (14):

\begin{array}{l} U_{m n m' n'}^{11} = U_{m n m' n'}^{22} = A_{m n m' n'} c_{n n'} \int_{θ_{1}}^{θ_{2}} d θ \sin θ F_{m m'} (r, θ) K_{m n m' n'}^{i i} (θ), \\ c_{n n'} = 1 + (^{- 1) n + n'}, \end{array}

U_{m n m' n'}^{33} = A'_{m n m' n'} c_{n n'} \int_{θ_{1}}^{θ_{2}} d θ \sin θ F'_{m m'} (r, θ) d_{0 m}^{n} (θ) d_{0 m'}^{n'} (θ),

\begin{array}{l} U_{m n m' n'}^{12} = - U_{m n m' n'}^{21} \\ = A_{m n m' n'} c'_{n n'} \int_{θ_{1}}^{θ_{2}} d θ \sin θ F_{m m'} (r, θ) K_{m n m' n'}^{i j} (θ), \\ c'_{n n'} = 1 + (^{- 1) n + n' + 1} . \end{array}

_{$F_{m m'} (r, θ)$}and

F'_{m m'} (r, θ)

in this case are piecewise smooth functions in

θ

. Figure 3 shows

F_{m m'} (r, θ)

of some index

m m'

at a certain

r_{p}

. A discontinuity (in the first derivative) occurs when the spherical shell intersects the hexagonal surface. As a result, the U-matrix integrand

F_{m m'} (r, θ) K_{m n m' n'}^{i i} (θ)

is no longer smooth in the interval

[θ_{1}, θ_{2}]

. The location of the discontinuity

θ = θ_{j u m p}

is given in Eq. (16) of [10], where

r_{p} \sin θ_{j u m p} = \frac{\sqrt{3} a}{2}

.

Fig. 4 An example $F_{m m'} (r, θ)$ as function of $θ$ for hexagonal column.

Download Full Size | PDF

The surface integrations in the zenith direction (Eqs. (25)~(30)) are naturally computed using Gauss-Legendre (GL) quadrature in the II-TM numerical implementations:

\int_{0}^{π} (d θ \sin θ) (F_{m m'} (r, θ) K_{m n m' n'}^{i j} (θ)) \approx \sum_{l = 1}^{N_{q}} w_{l} F {}_{m m'}{(r, θ_{l})} K_{m n m' n'}^{i j} (θ_{l}) .

In addition to the T-matrix truncation order N and the radial discretization

Δ r

, the number of quadrature nodes

N_{q}

is an important parameter controlling the accuracy of the II-TM method. For a standard Gaussian Legendre quadrature such as Eq. (31), the optimal nodes for quadrature of order

N_{q}

are given by the roots of the Legendre polynomials

P_{N_{q}} (x) = P_{N_{q}} (c o s θ)

. The GL quadrature is optimal if the integrand

(F_{m m'} (r, θ) K_{m n m' n'}^{i j} (θ))

can be well approximated by polynomials of order

2 N_{q} - 1

or less. Generally, this means that the various orders of derivative of the integrand are smooth. Discontinuities or singularities in the integrand or its derivatives will reduce the smoothness of the integrand. Higher order polynomials are required to approximate the fine scale oscillations of the integrand thus

N_{q}

can go to very large value. We use Gaussian quadrature on the following two integrals to illustrate our point that convergence rate of the Gaussian quadrature depends on behaviour of the integrand.

\int_{0}^{1} f_{1} (x) d x \approx \sum_{n = 1}^{N_{q}} w_{n} f_{1} (x_{n}), f_{1} (x) = {\begin{matrix} 1, 0 \leq x \leq 0.5, \\ 2 (1 - x), 0.5 < x \leq 1. \end{matrix}

\int_{0}^{1} f_{2} (x) d x \approx \sum_{n = 1}^{N_{q}} w_{n} f_{2} (x_{n}), f_{2} (x) = {\begin{matrix} 1, 0 \leq x \leq 0.5, \\ 0.5 [1 + \cos 2 π (x - 0.5)], 0.5 < x \leq 1. \end{matrix}

Fig. 5 $f_{1} (x)$ and $f_{2} (x)$ .

Download Full Size | PDF

Function $f_{2} (x)$ is smooth inside the [0,1] interval while $f_{1} (x)$ possesses a discontinuity in its first-order derivative at $x = 0.5$ . The behaviour of $f_{1} (x)$ is similar to that of $F_{m m'} (θ)$ in Fig. 4; they each have a discontinuity in the first-order derivative within their integration range. The two integrals of $f_{1} (x)$ and $f_{2} (x)$ have the same exact value (0.75), and we can demonstrate the convergence rate of Gaussian quadrature in terms of the rate of decrease of relative error in the numerical approximations. Table 1 shows the magnitude of the relative errors with respect to increasing the number of nodes $N_{q}$ for the two integrands. Apparently $f_{2} (x)$ shows a faster rate of convergence.

Table 1. Convergence rate of Gaussian quadrature for 2 different integrands.

View Table | View all tables in this article

Back to the II-TM method, the spheroid has the simplest integrand (in Eqs. (25) and (26)) among all other shapes:

\int_{0}^{θ_{2}} d θ \sin θ d_{0 m}^{n} (θ) d_{0 m}^{n'} (θ),

where the integrand is the product of Wigner-d functions

d_{0 m}^{n} (θ) d_{0 m}^{n'} (θ)

. We expect rapid convergence of the quadrature with respect to

N_{q}

, because the integrand

d_{0 m}^{n} (θ) d_{0 m}^{n'} (θ)

can be relatively well approximated by polynomials. Figure 6 shows the convergence rate of the outputs (

Q_{e x t}

, asymmetry factor (g factor),

P_{11} (180^{\circ})

) with respect to

N_{q}

for these different size parameters.

Q_{e x t}

is the extinction efficiency,

P_{11} (180^{\circ})

is the phase function in the exact backscattering direction and g factor is the asymmetry factor.

N_{q} = 1000

is chosen as the reference value and outputs of

N_{q} = 200, 400, 600, 800

are compared against it to demonstrate the convergence rate. T-matrix truncation order N and radial discretization

Δ r

are fixed. The T-matrix truncation order N for a particle with size parameter x is determined according to the formula

N = x + 4.05 x^{1 / 3} + 5,

in effect adding three more terms to the familiar choice used for spherical particles [16]

N_{M I E} = x + 4.05 x^{1 / 3} + 2,

No explicit formula has been published for calculations of the T-matrix for non-spherical particles, and we simply have made a choice that appears to increase the reliability of the results without greatly increasing computational cost. The increment in radius

Δ r

is determined according to the formula,

\frac{Δ r}{R - r_{0}} = 0.1,

_{$r_{0}$}is the radius of the homogeneous sphere inscribed inside the non-spherical particle, and is the starting radius for the radial recurrence. R is the radius of the circumscribed sphere. The value 0.1 has been found to be sufficient for particles without small scale features such as surface roughness; in the present cases where the particle shapes are regular spheroids and hexagonal columns, 0.1 should be sufficient.

The spheroid shape is indicated in the figures. x indicates the particle size parameter defined with respect to the radius of the circumscribed sphere of the particle.

Fig. 6 Convergence rate of the outputs ( $Q_{e x t}$ , asymmetry factor (g factor), $P_{11} (180^{\circ})$ ) with respect to number of quadrature nodes $N_{q}$ . Six panels correspond to 6 different spheroidal shapes indicated with size parameter x and aspect ratio a/b.

Download Full Size | PDF

The convergence rates of different sizes show mixed behaviors. Generally, we can conclude that the convergence rate for spheroid is fast. In most cases, the outputs agree to more than 12 significant digits. The obstacle for spheroid is the intense oscillation of high-order Wigner-d functions for larger spheroids. Once $N_{q}$ gets large enough (>200), the sampling rate is sufficiently large and we have rapid convergence. The integrals in the hexagonal column formulation (in Eqs. (28)~(30)) take the following form:

\int_{θ_{1}}^{θ_{2}} d θ \sin θ F_{m m'} (r, θ) K_{m n m' n'}^{i j} (θ) .

It was pointed out in Fig. 4 that the function

F_{m m'} (r, θ)

in the above integral has a discontinuity in its first derivative when the spherical shell crosses the hexagonal column. We expect this discontinuity to undermine the convergence of the quadrature if we do the quadrature in the interval

[θ_{1}, θ_{2}]

which contains the discontinuity. Figure 6 shows the convergence rate of the outputs (the extinction efficiency

Q_{e x t}

, the asymmetry factor (the g factor),

P_{11} (180^{o})

with respect to

N_{q}

.

N_{q} = 1000

is again used as the reference value. The hexagonal column shape is indicated in the figure and x is the size parameter.

Fig. 7 Convergence rate of the outputs ( $Q_{e x t}$ , asymmetry factor (g factor), $P_{11} (180^{\circ})$ ) with respect to number of quadrature nodes $N_{q}$ . Six panels correspond to 6 different hexagonal column shapes indicated with size parameter x and aspect ratio a/b.

Download Full Size | PDF

Compared to spheroids (Fig. 6), the convergence rate is much slower. The discontinuity in $F_{m m'} (r, θ)$ reduces the convergence rate of the quadrature. Note that there is only one discontinuity $θ_{j u m p}$ in $[θ_{1}, θ_{2}]$ . A similar situation was encountered in the EBCM method for hexagonal columns [17]. To address the problem of discontinuities, a natural solution is to split the interval and perform quadrature on the two intervals where the integrand is smooth (Eq. (41).

\begin{array}{l} \int_{θ_{1}}^{θ_{2}} d θ \sin θ F_{m m'} (r, θ) K_{m n m' n'}^{i j} (θ) = \\ \int_{θ_{1}}^{θ_{j u m p}} d θ \sin θ F_{m m'} (r, θ) K_{m n m' n'}^{i j} (θ) + \int_{θ_{j u m p}}^{θ_{2}} d θ \sin θ F_{m m'} (r, θ) K_{m n m' n'}^{i j} (θ) . \end{array}

After making this modification, the results are shown in Fig. 8. Compared to Fig. 6, the convergence rate is accelerated. At this rate, we can expect to reach convergence to 16 digits at about

N_{q} = 1000

. This justifies the use

N_{q} \approx 1000

in practical applications.

Fig. 8 Convergence rate of the outputs ( $Q_{e x t}$ , asymmetry factor (g factor), $P_{11} (180^{\circ})$ ) with respect to number of quadrature nodes $N_{q}$ . Six panels correspond to 6 different hexagonal column shapes indicated with size parameter x and aspect ratio a/b. Unlike Fig. 7, the integration range is split into 2 smooth intervals so the convergence rate is accelerated (see the text).

Download Full Size | PDF

From Figs. 6 and 8, we can conclude that, with respect to the number of quadrature points $N_{q}$ , to reach convergence for common particle sizes and aspect ratios, $N_{q} \approx 1000$ is sufficient for hexagonal columns and $N_{q} \approx 600$ is sufficient for spheroids.

4. Accelerating the node and weight generating scheme for Gaussian quadrature in II-TM

As stated above, the optimal nodes for Gaussian-Legendre quadrature of order L are the roots of Legendre Polynomials $P_{L} (x)$ . There are no closed formulas for location of the roots so they have to be generated numerically along with the corresponding weights. Different algorithms exist for generating the nodes and weights of Gauss-Legendre quadrature [18,19].

The current algorithm in the II-TM is based on the Newton’s method (Newton-Raphson) to find the roots of $P_{L} (x)$ . Denote the L node-weight pairs ${x_{L, k}, w_{L, k}}, k \in [0, L]$ and they satisfy $P_{L} (x_{L, k}) = 0$ . The weights are given by $w_{L, k} = \frac{2 (1 - x_{L, k}^{2})}{{[(n + 1) P_{L + 1} (x_{L, K})]}^{2}}$ . For large L, some initial value $x_{L, k, 0}$ for $x_{L, k}$ is given and then Newton-Raphson iterations are performed to obtain a true value for $x_{L, k}$ up to some precision. A key in the design of fast algorithms is proper choice of asymptotic expansions (large L) of the Legendre polynomials. Considerable research effort has gone into the search for suitable asymptotic expansions [12,19–21].

The state-of-the-art for Gaussian-Legendre quadrature is the Bogaert’s algorithm [12,22]. In their algorithm, the Newton-Raphson iteration is avoided by providing asymptotic expansions directly for $x_{L, k}$ and $w_{L, k}$ . This algorithm is much faster than previous methods. In the following discussions, a commonly used algorithm is referred to as “Newton” and the new algorithm is referred to as “Bogaert”. We first illustrate the difference in the performance of the algorithms in the simple case of the model function $f_{1} (x)$ of the previous section (Fig. 5 and Table 1). Table 2 shows the CPU time required to compute the numerical quadrature using the two algorithms. The Bogaert algorithm is clearly much faster. Obviously, the Bogaert algorithm is much faster.

Table 2. CPU time required to run the same code section with Newton and Bogaert algorithms.

View Table | View all tables in this article

We now turn to implementation of the Bogaert method in conjunction with the II-TM. For L>100, the Bogaert asymptotic expansions are used; for L<100, tabulated values for nodes and weights are used. Table 3 shows the acceleration for hexagonal columns and spheroids. II-TM is implemented with the Massage Passing Interface (MPI) parallelization and 100 CPUs were used for the experiments. Significant acceleration is achieved for hexagonal columns of all sizes and aspect ratios. Moderate acceleration is achieved for spheroid in terms of CPU time. However, in two cases, the Bogaert algorithm slightly slows down the method in terms of wall clock time. As mentioned above, the spheroid and hexagon have different symmetries, and correspondingly the code structures of the II-TM implementations are different. As a result, the respective MPI parallelizations are different. The result is that the corresponding acceleration rates associated with the Bogaert algorithm are different. The optimization of the Bogaert algorithm in the II-TM MPI implementations for a wide range of particle geometries is subject of future work.

Table 3. Acceleration for hexagonal column and spheroids of various sizes and aspect ratios.

View Table | View all tables in this article

With the accelerated code, we compute the single scattering-properties of large hexagonal columns with size parameter (defined with the circumscribed sphere radius) x = 225 and aspect ratio = 1. Two cases are considered. One has the index of refraction $m = 1.308 + i 1.43 \times 10^{- 9}$ which corresponds to ice refractive index at 650 nm wavelength. Another one has the index of refraction $m = 1.276 + i 0.413$ which corresponds to ice refractive index at 12 $μ m$ wavelength. Random orientation is assumed. Figures 9 and 10 show comparison between the II-TM method and a physical geometric optics method (PGOM) [23]. Very good agreement is achieved between the two methods.

Fig. 9 Phase matrix elements for hexagonal column with size parameter x = 225 (defined with the circumscribed sphere radius). Comparison is between PGOM (blue) and II-TM (red dashed). Particle refractive index is $m = 1.308 + i 1.43 \times 10^{- 9}$ .

Download Full Size | PDF

Fig. 10 Phase matrix elements for hexagonal column with size parameter x = 225 (defined with the circumscribed sphere radius). Comparison is between PGOM (blue) and II-TM (red dashed). Particle refractive index is $m = 1.276 + i 0.413$ .

Download Full Size | PDF

5. Summary

We studied the performance of the angular quadrature in the II-TM method for spheroid and hexagonal column. The precision of the angular quadrature is of fundamental importance to the II-TM method and we want to achieve convergence by choosing an appropriate number of quadrature nodes $N_{q}$ . In the case of a hexagonal column, the previous formulation [9] is modified to avoid a discontinuity in the integration domain, and the convergence rate is largely accelerated. Based on convergence test results from spheroids and hexagonal columns with various sizes, we conclude that $N_{q} \approx 600$ for spheroid and $N_{q} \approx 1000$ for hexagonal column is sufficiently large for the improved numerical implementation. A new node and weight generating algorithm (Bogaert algorithm) is implemented in the II-TM method, and considerable acceleration in computation time is achieved for hexagonal columns and spheroids. With the improved II-TM method, we compute the single-scattering properties of two large hexagonal columns with size parameter x = 225 and compare the results against the physical geometric optics method [23]. Good agreement is found between the two methods.

Funding

National Science Foundation (AGS-1826936); Texas A&M David Bullock Harris Chair in Geosciences endowment funds (02-512231-10000).

Acknowledgment

We gratefully acknowledge the computational support of the Texas A&M High Performance Computing Center.

References

1. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley & Sons, 1983).

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

3. P. C. Waterman, “Matrix formulation of electromagnetic scattering,” Proc. IEEE 53(8), 805–812 (1965). [CrossRef]

4. P. C. Waterman, “The T-matrix revisited,” J. Opt. Soc. Am. A 24(8), 2257–2267 (2007). [CrossRef] [PubMed]

5. N. G. Khlebtsov, “Orientational averaging of light-scattering observables in the J-matrix approach,” Appl. Opt. 31(25), 5359–5365 (1992). [CrossRef] [PubMed]

6. M. I. Mishchenko, “Light scattering by randomly oriented axially symmetric particles,” J. Opt. Soc. Am. A 8(6), 871–882 (1991). [CrossRef]

7. T. Wriedt, “Review of the null-field method with discrete sources,” J. Quant. Spectrosc. Radiat. Transf. 106(1-3), 535–545 (2007). [CrossRef]

8. M. I. Mishchenko and L. D. Travis, “Capabilities and limitations of a current FORTRAN implementation of the T-matrix method for randomly oriented, rotationally symmetric scatterers,” J. Quant. Spectrosc. Radiat. Transf. 60(3), 309–324 (1998). [CrossRef]

9. L. Bi, P. Yang, G. W. Kattawar, and M. I. Mishchenko, “Efficient implementation of the invariant imbedding T-matrix method and the separation of variables method applied to large nonspherical inhomogeneous particles,” J. Quant. Spectrosc. Radiat. Transf. 116, 169–183 (2013). [CrossRef]

10. L. Bi and P. Yang, “Accurate simulation of the optical properties of atmospheric ice crystals with the invariant imbedding T-matrix method,” J. Quant. Spectrosc. Radiat. Transf. 138, 17–35 (2014). [CrossRef]

11. T. Wriedt and Y. Eremin, The Generalized Multipole Technique for Light Scattering. (Springer, 2018).

12. I. Bogaert, “Iteration-Free Computation of Gauss–Legendre Quadrature Nodes and Weights,” SIAM J. Sci. Comput. 36(3), A1008–A1026 (2014). [CrossRef]

13. B. R. Johnson, “Invariant imbedding T matrix approach to electromagnetic scattering,” Appl. Opt. 27(23), 4861–4873 (1988). [CrossRef] [PubMed]

14. M. Mischenko, L. Travis, and A. Lacis, Scattering, absorption, and emission of light by small particles, (Cambridge University Press,1970).

15. W. R. C. Somerville, B. Auguié, and E. C. Le Ru, “Simplified expressions of the T-matrix integrals for electromagnetic scattering,” Opt. Lett. 36(17), 3482–3484 (2011). [CrossRef] [PubMed]

16. W. J. Wiscombe, “Improved Mie scattering algorithms,” Appl. Opt. 19(9), 1505–1509 (1980). [CrossRef] [PubMed]

17. F. M. Kahnert, J. J. Stamnes, and K. Stamnes, “Application of the extended boundary condition method to particles with sharp edges: a comparison of two surface integration approaches,” Appl. Opt. 40(18), 3101–3109 (2001). [CrossRef] [PubMed]

18. A. Glaser, X. Liu, and V. Rokhlin, “A Fast Algorithm for the Calculation of the Roots of Special Functions,” SIAM J. Sci. Comput. 29(4), 1420–1438 (2007). [CrossRef]

19. G. H. Golub and J. H. Welsch, “Calculation of Gauss Quadrature Rules,” Math. Comput. 23(106), 221–230 (1969). [CrossRef]

20. N. Hale and A. Townsend, “Fast and Accurate Computation of Gauss–Legendre and Gauss–Jacobi Quadrature Nodes and Weights,” SIAM J. Sci. Comput. 35(2), A652–A674 (2013). [CrossRef]

21. A. Townsend, T. Trogdon, and S. Olver, “Fast computation of Gauss quadrature nodes and weights on the whole real line,” IMA J. Numer. Anal. 36(1), 337–358 (2014).

22. A. Townsend, “The race to compute high-order Gauss-Legendre quadature,” SIAM J. Sci. Comput. (2015).

23. B. Sun, P. Yang, G. W. Kattawar, and X. Zhang, “Physical-geometric optics method for large size faceted particles,” Opt. Express 25(20), 24044–24060 (2017). [CrossRef] [PubMed]

	$f_{1} (x)$	$f_{2} (x)$
	Relative error = (Estimate – Exact)/Exact
$N_{q}$	order of magnitude of relative error	order of magnitude of relative error
10	$10^{- 3}$	$10^{- 9}$
100	$10^{- 5}$	$10^{- 9}$
1000	$10^{- 7}$	$10^{- 9}$
10000	$10^{- 9}$	$10^{- 9}$

$N_{q}$	$f_{1} (x)$		CPU time acceleration (Newton/Bogaert)
$N_{q}$	CPU time(sec) Newton	CPU time(sec) Bogaert	CPU time acceleration (Newton/Bogaert)
10	<0.001	<0.001	–
100	<0.001	<0.001	–
1000	0.026	<0.001	–
10000	1.643	0.001	1643

Shape	Size parameter	Aspect ratio	CPU time acceleration (Newton/Bogaert)	Wall clock time acceleration (Newton/Bogaert)
Hexagonal columns	x = 30	0.5	7.9	283.7
	x = 30	0.75	9.5	432
	x = 65	0.5	2.2	34.1
	x = 65	0.75	2.3	34.6
	x = 100	0.5	1.7	8.6
	x = 100	0.75	1.6	9.1

spheroids	x = 30	0.5	1.1	0.9
	x = 30	0.75	1.03	1.1
	x = 65	0.5	1.03	1.0
	x = 65	0.75	1.2	1.1
	x = 100	0.5	1.1	0.8
	x = 100	0.75	1.03	1.0

	$f_{1} (x)$	$f_{2} (x)$
	Relative error = (Estimate – Exact)/Exact
$N_{q}$	order of magnitude of relative error	order of magnitude of relative error
10	$10^{- 3}$	$10^{- 9}$
100	$10^{- 5}$	$10^{- 9}$
1000	$10^{- 7}$	$10^{- 9}$
10000	$10^{- 9}$	$10^{- 9}$

$N_{q}$	$f_{1} (x)$		CPU time acceleration (Newton/Bogaert)
$N_{q}$	CPU time(sec) Newton	CPU time(sec) Bogaert	CPU time acceleration (Newton/Bogaert)
10	<0.001	<0.001	–
100	<0.001	<0.001	–
1000	0.026	<0.001	–
10000	1.643	0.001	1643

Improvements in the computational efficiency and convergence of the Invariant Imbedding T-matrix method for spheroids and hexagonal prisms

Abstract

1. Introduction

2. The T-matrix method

3. Improving the convergence of Gaussian quadrature in II-TM

4. Accelerating the node and weight generating scheme for Gaussian quadrature in II-TM

5. Summary

Funding

Acknowledgment

References

Cited By

Figures (10)

Tables (3)

Equations (39)

Optics Express