Invariant-imbedding T-matrix method

Invariant-imbedding T-matrix method 4 Light scattering by spherical particles can be analytically solved using the LorenzMie theory (van de Hulst, 1...

Download PDF

669KB Sizes 0 Downloads 90 Views

Report

PDF Reader
Full Text

Invariant-imbedding T-matrix method

4

Light scattering by spherical particles can be analytically solved using the LorenzMie theory (van de Hulst, 1957; Bohren and Huffman, 1983). Moreover, the scattering by a spherical particle or a spherical shell has been extensively studied (e.g., Aden and Kerker, 1951; Kattawar and Hood, 1976; Wiscombe, 1980; Toon and Ackerman, 1981; Hovenac and Lock, 1992; Nussenzveig, 1992; Gouesbet and Grehan, 2011). The Lorenz-Mie theory can solve the scattering of light by a spherical particle for the entire practical size range. Importantly, the cpu time usage for the calculation in the case of a spherical particle is so minimal that it does not pose a concern given the present computational resources. For nonspherical particles, however, only scattering by small particles can be obtained using numerically accurate methods. Here, a numerically accurate method refers to a method that gives converged results with a sufficiently fine resolution. In Section 2.3, several frequently used accurate methods, such as FDTD and DDA, are introduced for the computation of the scattering of light by nonspherical particles. With an increasing particle size, not only do the cpu time usage and the computer memory demand significantly increase, but the convergence rates for these methods also substantially degrade. Moreover, the single-scattering properties of a particle calculated by any numerically accurate method depend on the incident direction with respect to the orientation of the scattering particle. Accordingly, the scattering properties of a particle under the random orientation condition can only be obtained by numerically averaging the results associated with different orientations. The number of orientations necessary for convergence drastically increases with increased particle size. Consequently, these numerically accurate methods are only efficient for light scattering by small particles. The T-matrix method, however, is a semianalytical method expressed in terms of vector spherical wave functions. The T-matrix of a particle is also independent of the incident direction. Accordingly, the random orientation condition can be analytically accounted for once the T-matrix is given. The T-matrix concept is given in Chapter 3, including introduction of the vector spherical wave function and analytical computation for random orientations using the T-matrix. A commonly used method, the extended boundary condition method (EBCM), which uses a surface integral equation and matrix inversion, is introduced in Section 3.3.2. The EBCM has been extensively used to calculate light scattering by particles with axial symmetry, such as spheroids, cylinders, and Chebyshev particles (Mishchenko et al., 2002). Over a period of time, the EBCM was almost a synonym of the “T-matrix method.” In principle, the EBCM is a general theory to calculate the T-matrix and can be applied to any symmetric or

Invariant Imbedding T-matrix Method for Light Scattering by Nonspherical and Inhomogeneous Particles https://doi.org/10.1016/B978-0-12-818090-7.00004-8 © 2020 Elsevier Inc. All rights reserved.

146

Invariant Imbedding T-matrix Method

asymmetric particle, but this method is usually applied to the computation for particles with axial symmetry (e.g., Mugnai and Wiscombe, 1986; Mishchenko and Travis, 1998; Havemann and Baran, 2001). One reason for this is that effort of computing the T-matrix of a particle with axial symmetry significantly reduces due to the T-matrix symmetry (Schulz et al., 1999; Kahnert et al., 2001). Another reason is related to numerical stability. The computation of the T-matrix includes a matrix inversion process, that is, T ¼ (RgQ)Q1. Matrices RgQ and Q are both calculated from the surface integrals introduced in Section 3.3.2. For nonspherical particles with large size parameters or largely deformed geometries, such as those having extreme aspect ratios, the eigenvalues of matrix Q differ by many orders of magnitude, and the matrix Q becomes ill-conditioned due to the fact that the surface integrals introduce the spherical Hankel function of the first kind in different arguments into the matrix, and instability of the matrix inversion process becomes worse if the dimension of the matrix Q becomes larger. Several efforts have been made to improve the numerical stability of the T-matrix computation. A perturbation approach was developed for the application of EBCM, which replaces the matrix inversion problem with an implicit T-matrix equation, and the implicit equation can be solved iteratively (Kahnert and Rother, 2011). However, the perturbation method only works for the particle geometry slightly deviating from a reference geometry. A straightforward improvement for the numerical stability results from using the extended precision variables to replace double precision variables in the computation and inversion of matrix Q (Mishchenko and Travis, 1994). Computations using extended precision may be roughly five to eight times slower than using double precision, but the applicable particle sizes can be two or three times larger if extended precision is used. Another method uses the unitary property of the T-matrix introduced in Section 3.4.1 and iteratively obtains stable results for strongly deformed geometries (e.g., Lakhtakia et al., 1984), where a large difference between the smallest and largest dimensions is needed to characterize the particle. However, this unitary equation is only applicable to particles without absorption. A method using discrete sources can also be implemented in the EBCM to improve numerical stability for prolate and oblate particles. One method is called the iterative EBCM (IEBCM) (Iskander et al., 1983, 1989), and another is the null field with discrete sources (e.g., Doicu et al., 2006). However, although the numerical stability is attained by such methods, the computational complexity and cpu time increase compared with the original EBCM. An alternative method is required not only to retain numerical stability but also to preserve numerical efficiency, so the method can be applied to particles with regular or extreme geometries and large size parameters. The invariant-imbedding T-matrix method (IITM) is a volume integral method that uses the invariant-imbedding technique to iteratively obtain the T-matrix of an arbitrary particle. The volume integral equation of the electric field has the form of the Fredholm equation of the second kind, which is usually well conditioned (Press et al., 2007). The IITM shows numerical stability in matrix inversion and

Invariant-imbedding T-matrix method

147

can be used to obtain the T-matrix of particles with large size parameters or extreme geometries. The invariant-imbedding technique was introduced by Ambarzumian (1943). One of the most famous applications of the invariant-imbedding technique is reported in Radiative Transfer by Chandrasekhar (1960). The general concept of invariant imbedding is discussed in detail in An Introduction to Invariant Imbedding by Bellman and Wing (1975). The concept of “invariant imbedding” was first applied in light scattering to obtain the T-matrix of an arbitrary particle by Johnson (1988). Unfortunately, this paper received little attention from the light scattering community and was cited only several times in the literature before 2013. Johnson’s paper was fortunately revisited, and the IITM was implemented as a powerful light scattering computational program by Bi et al. (2013). The IITM shows stability not only for particles with axial symmetry but also for particles without axial symmetry (Bi and Yang, 2014; Yang et al., 2019). Moreover, it is straightforward to extend IITM from homogeneous to inhomogeneous particles and from particles with symmetry to particles without symmetry. The current computational range of the IITM has been applied to particle geometries with axial symmetry, finite-fold rotational symmetry (such as a hexagonal column with sixfold rotational symmetry), or even with no symmetry (Bi and Yang, 2014; Liu et al., 2015; Sun et al., 2017; Yang et al., 2019). This section presents a comprehensive introduction to the IITM. Section 4.1 briefly derives the volume integral equation from Maxwell’s equations and reorganizes the expansion of the electromagnetic fields using the vector spherical wave functions. In addition, the free space dyadic Green function is revisited using eigenfunction expansions and can easily be generalized to other sets of eigenfunction expansions, such as expansion using both the outgoing and incoming vector spherical wave functions used in the Debye series study. The invariant-imbedding technique is introduced in Section 4.2. Section 4.3 gives the concept of the IITM, and discussions about several aspects of the IITM are given in Section 4.4.

4.1

Electromagnetic volume integral equation

4.1.1 Volume integral equation Here, only dielectric particles are considered, and the particles are assumed to be isotropic and have a linear response to an applied field. Maxwell’s equations in such a medium for time-harmonic fields can be given as follows: ! ! r E r ¼ iωμH r ,

(4.1.1a)

! ! r H r ¼ iωεE r ,

(4.1.1b)

148

Invariant Imbedding T-matrix Method

Fig. 4.1 A dielectric scattering particle. The region occupied by the particle is denoted as V1, while the exterior region is denoted as V0.

where ε and μ are the permittivity and the permeability in the medium, respectively. These are the general equations for the electromagnetic field in a medium. For light scattering of a particle in a surrounding medium as in Fig. 4.1, the permittivity of the !

particle is ε1 r , which may vary inside the particle, and the permittivity of the

surrounding medium is ε. Also, we assume that both the particle and the surrounding medium are nonferromagnetic and thus μ1 ¼ μ ¼ μ0. According to Maxwell’s equations in Eq. (4.1.1), we have the following equations: ! ! ! r r E r k2 E r ¼ 0, r 2 V0 ,

(4.1.2a)

! ! ! r r E r k12 E r ¼ 0, r 2 V1 ,

(4.1.2b)

! ! ! E r ¼ Einc r + Esca r , pﬃﬃﬃﬃﬃﬃﬃ where k ¼ ω εμ0 and k1 ¼ ω

(4.1.3)

rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ! ε1 r μ0 ; the space is divided into the region occu-

pied by the particle V1 and the region exterior to the particle V0. It is evident that the incident field always satisfies the equation ! ! r r Einc r k2 Einc r ¼ 0:

(4.1.4)

Eq. (4.1.2) can be rewritten in a compact form (Mishchenko et al., 2002): ! ! ! r r E r k2 E r ¼ j r ,

(4.1.5)

Invariant-imbedding T-matrix method

149

where ! j r ¼

(

!

0, h r 2 V0 , i ! ! ! k2 ε1 r =ε 1 E r , r 2 V1 :

(4.1.6)

The free space dyadic Green function satisfies the following equation: $ ! ! $ ! ! $ ! ! r r G r , r 0 k2 G r , r 0 ¼ I δ r r 0 ,

(4.1.7)

!

where the operator r acts only on the vector r . We use the vector Green theorem to derive the volume integral equation. It is rewritten here as (Morse and Feshbach, 1953; Tsang et al., 2000; Mishchenko et al., 2002) ð

þ dV½A ðr r BÞ B ðr r AÞ ¼ dS^ ns ðB r A A r BÞ,

V

S

(4.1.8) where volume V is enclosed by surface S and A and B are arbitrary vectors. Using Eq. (3.3.19), the identity of the incident field at infinity is þ nh i $ h $! ! io ! ! ! ! E r ¼ dS^ ns r Einc r G r , r 0 + Einc r r G r , r 0 : inc !0

S∞

(4.1.9)

Next, we can insert the total electric field as the vector A as ! A¼E r ,

(4.1.10a)

and B in this case is $ ! ! ! B ¼ G r , r0 a ,

(4.1.10b)

Substituting Eq. (4.1.10) into Eq. (4.1.8) and integrating it over all space, the vector Green theorem applied to Eqs. (4.1.5)–(4.1.7) and Eq. (4.1.9) can give the following expression: ð $ !0 ! ! ! ! inc !0 E r ¼E r + j r G r , r 0 d3 r V1

150

Invariant Imbedding T-matrix Method

þ

dS^ ns

nh i $ h $ ! ! io ! ! ! ! r Esca r G r , r 0 + Esca r r G r , r 0 :

S∞

(4.1.11) Using the Sommerfeld radiation condition with respect to the dyadic Green function given in Eq. (3.3.22), the surface integration at infinity in Eq. (4.1.11) goes to zero. Using the symmetric property of the dyadic Green function given in Eq. (3.3.5), Eq. (4.1.11) gives the volume integral equation for the electric field: ð $ ! ! ! ! ! ! ! E r ¼ Einc r + u r 0 G r , r 0 E r 0 d3 r 0 ,

(4.1.12a)

V1

ð ! ! $ ! ! ! ! Esca r ¼ u r 0 G r , r 0 E r 0 d 3 r 0 ,

(4.1.12b)

V1

h i ! ! where u r 0 ¼ k2 ε1 r 0 =ε 1 .

4.1.2 Equivalence between volume integral equation and surface integral equation In Section 3.3.2, the surface integral equation of the electric field is given in Eq. (3.3.23) using the vector Green theorem and is reorganized as follows: þ n h $ i h $ ! !io ! ! ! ! E r ¼ dS0 n^s iωμ0 H r 0 G r 0 , r + E r 0 r0 G r 0 , r , sca !

S

(4.1.13) !

!

where the field point r is the region outside the scattering particle, that is, r 2 V0 ; S encloses the scattering particle and n^s is the outward unit normal to the surface; and ! r0 acts on the variable r 0 . Both the volume integral equation shown in Eq. (4.1.12) and surface integral equation shown in Eq. (4.1.13) are derived from Maxwell’s equations and the vector Green theorem, which should be equivalent. Using Gauss’s theorem, the surface integral term in Eq. (4.1.13) becomes the volume integral as follows: ð n h $ i h $ ! !io ! ! ! ! ! ! Esca r ¼ d3 r 0 r0 iωμ0 H r 0 G r 0 , r + E r 0 r0 G r 0 , r : V1

(4.1.14)

Invariant-imbedding T-matrix method

151

This divergence operator identity will be useful in what follows: r ðA BÞ ¼ ðr AÞ B A ðr BÞ:

(4.1.15)

Using Maxwell’s equations in Eqs. (4.1.1), (4.1.2) and Eq. (4.1.7), Eq. (4.1.14) can be further transformed: ð n $ $ o ! ! ! ! ! ! ! ! ! ! ! Esca r ¼ d 3 r 0 k12 E r 0 G r 0 , r k2 E r 0 G r 0 , r E r 0 δ r r 0 : V1

(4.1.16)

! ! ! Since the field point is outside region V1, the integral over E r 0 δ r r 0 becomes zero. Also, the free space dyadic Green function has the symmetry property: $ ! ! $ ! ! G r 0, r ¼ GT r , r 0 :

(4.1.17)

Eq. (4.1.16) becomes the volume integral equation: ð $ ! ! ! ! r ¼ k2 ðε1 =ε 1ÞG r , r 0 E r 0 d3 r 0 : E sca !

(4.1.18)

V

The equivalence between the volume integral equation and the surface integral equation is proved using the Green theorem. In the far field or r ! ∞, the free space dyadic Green function is given in Eq. (3.3.8). Also, for an arbitrary vector A, one has $ I r^r^ A ¼ ^ r r^ A:

(4.1.19)

The volume integral and the surface integral equations can be given in the far field as ! Esca r

r!∞

¼

exp ðikr Þ ik3 ikr 4π

ð

1 m2

$ ! ! ! I r^r^ E r 0 exp ik^ r r 0 d3 r 0 ,

V1

(4.1.20) E r sca !

r!∞

þ nh i ωμ h io exp ðikr Þ k2 ! ! ! ¼ , r^ d2 r 0 n^s E r 0 0 r^ n^s H r 0 ikr 4π k S

(4.1.21) where m2 ¼ ε1/ε is the square of the relative refractive index of the particle.

152

Invariant Imbedding T-matrix Method

4.1.3 Free space dyadic green function revisited The free space dyadic Green function is given using the so-called Ohm-Rayleigh method in Section 3.3.1. The method first expands a corresponding δ function using vector spherical wave functions and then finds the expansion of the free space dyadic Green function. This method is used to get the expansion of the free space dyadic Green function in various coordinate systems, such as Cartesian, cylindrical, and spherical coordinates (Tai, 1994). Contour integration is used to get the corresponding expansions. When the expansion eigenfunctions change, the contour integration may be complicated. Here, we want to use a general technique to expand the Green function so the eigenfunctions can be any combination as long as their Wronskian is not equal to zero, which is introduced by Tai (1994). Since the expansion of the scalar Green function is straightforward using this technique, we shall start with the expansion of the scalar Green function and then find the expansion of the dyadic Green function using the relation between the scalar and the dyadic Green functions. Note that only the expansion in spherical coordinates is discussed in this book and the expansion in other coordinate systems can be found in Morse and Feshbach (1953) and Tai (1994). The direct expansion of the dyadic Green function using this general technique can also be found in Morse and Feshbach (1953). For convenience, the scalar Green function and its expression given in Eq. (3.3.1) are restated here. The scalar Green function for the scalar Helmholtz equation satisfies the equation ! ! ! ! ! ! r2 g r , r 0 + k2 g r , r 0 ¼ δ r r 0 ,

(4.1.22)

The explicit expression for the scalar Green function is ( Jackson, 1998)

!

!

g r , r0

i h ! ! exp ik r r 0 ¼ ! ! : 4π r r 0

(4.1.23)

We will first expand the scalar Green function using scalar spherical wave functions and then discuss how to expand the free space dyadic Green function in terms of the scalar expansions. A general technique for computing the scalar Green function is given in Chapter 7 of Morse and Feshbach (1953). Using this technique, not only can the Green function in different coordinate systems be expanded, but the Green function can also have different basis functions. Here, the scalar Green function is given in spherical coordinates using spherical Bessel functions and spherical Hankel functions of the first kind as the basis set. Since two angular coordinates (θ, φ) have a finite domain, the spherical harmonics function associated with angular coordinates forms a complete set of two-dimensional eigenfunctions, and the eigenvalues are discrete. The Green function can be expanded in terms of the spherical harmonics, and the inhomogeneous equation related to the radial coordinate can be solved using the two

Invariant-imbedding T-matrix method

153

solutions of the corresponding homogeneous equations. The spherical harmonics satisfy the orthogonality relation:

γ 0mn

2

ð

∗n0 m0 ðθ, φÞnm ðθ, φÞsin θdθdφ ¼ δnn0 δmm0 ,

(4.1.24a)

4π

γ 0mn

2

¼

2n + 1 ðn mÞ! : 4π ðn + mÞ!

(4.1.24b)

The spherical harmonics are defined in Eqs. (3.2.11), (3.2.12). The scalar Green function can be expanded using the spherical harmonics and the undetermined functions as ∞ X n X ! ! g r , r0 ¼ Rn ðr, r 0 Þmn ðθ0 , φ0 Þnm ðθ, φÞ,

(4.1.25)

n ¼ 0 m ¼ n

where the functions Rn(r, r0 ) and nm ðθ0 , φ0 Þ are undetermined. The radial function Rn(r, r0 ) can be shown later to be only related to index n. Substituting Eq. (4.1.25) into Eq. (4.1.22) and using Eqs. (3.2.5), (3.2.8a), the equation for the scalar Green function is

1 ∂ 2∂ nð n + 1Þ 0 2 0 r R nm ðθ0 , φ0 Þnm ðθ, φÞ ð r, r Þ + k ð r, r Þ R n n r 2 ∂r ∂r r2 n ¼ 0 m ¼ n ∞ X n X

¼

δðr r 0 Þδðθ θ0 Þδðφ φ0 Þ : r 2 sin θ

(4.1.26)

Using the orthogonality of spherical harmonics in Eq. (4.1.24), we can decouple Eq. (4.1.26) into two equations: 2 nm ðθ0 , φ0 Þ ¼ γ 0mn ∗nm ðθ0 , φ0 Þ, "

# " # 2 ∂ nðn + 1 Þ δð r r 0 Þ 0 Rn ðr, r Þ + 1 Rn ðr, r 0 Þ ¼ + : 2 2 ðkr Þ ∂ðkr Þ ∂ðkr Þ ðkr Þ ðkr Þ2 ∂2

(4.1.27) (4.1.28)

Eq. (4.1.28) is an inhomogeneous equation for the radial function, and the homogeneous equation is actually the spherical Bessel equation based on Eq. (3.2.8b). Before giving the solution of Eq. (4.1.28), let us first discuss the general form of a second-order ordinary differential equation (ODE). Homogeneous and corresponding inhomogeneous second-order ODEs are stated as d 2 y ðx Þ dyðxÞ + qðxÞ ¼ 0, + pðxÞ 2 dx dx

(4.1.29a)

154

Invariant Imbedding T-matrix Method

d 2 y ðx Þ dyðxÞ + qðxÞ ¼ r ðxÞ, + pð x Þ 2 dx dx

(4.1.29b)

and y1(x) and y2(x) are two linearly independent solutions, that is, the Wronskian is not zero, or Δðy1 , y2 Þ ¼ y1 y02 y2 y01 6¼ 0:

(4.1.30)

Accordingly, the solution of the inhomogeneous second-order ODE can be explicitly expressed using the two independent solutions from the homogeneous second-order ODE (Morse and Feshbach, 1953): ðx yðxÞ ¼ y2 ðxÞ dt

ð r ðtÞy1 ðtÞ r ðtÞy2 ðtÞ + y1 ðxÞ dt , Δðy1 ðtÞ, y2 ðtÞÞ Δ y ð 1 ðtÞ, y2 ðtÞÞ x

(4.1.31)

where the lower and upper limit choices depend on the two independent solutions y1(x) and y2(x) and the boundary condition of y(x). Note that the solution in Eq. (4.1.31) exists only when the coefficient of the second-order derivative in Eq. (4.1.29) is unity. Comparing Eq. (4.1.28) and Eq. (4.1.29), the inhomogeneous term here is the right-hand side of Eq. (4.1.28). Accordingly, the solution of Eq. (4.1.28) can be given as Rn ðr, r 0 Þ ¼ z2n ðr Þ

ðr du a

δðu r 0 Þz1n ðuÞ ðkr 0 Þ2 Δðz

1n ðuÞ, z2n ðuÞÞ

+ z1n ðr Þ

ðb du r

δðu r 0 Þz2n ðuÞ

, ðkr 0 Þ2 Δðz1n ðuÞ, z2n ðuÞÞ (4.1.32)

where the choices that the first limit is less than r and the second limit is larger than r are allowed because boundary conditions are not determined. For r > r0 , the second integral in Eq. (4.1.32) becomes zero, while for r < r0 , the first integral becomes zero. Accordingly, the solution in radial coordinates can be written as 0

Rn ðr, r Þ ¼

k ðkr 0 Þ2 Δðz1n ðr 0 Þ, z2n ðr 0 ÞÞ

(

z1n ðr Þz2n ðr 0 Þ; r r 0 , z2n ðr Þz1n ðr 0 Þ; r r 0 ,

(4.1.33)

where the property kδ(kx) ¼ δ(x) since k > 0 is used. For our light scattering formulation, the boundary condition requires the radial function to be finite at the origin (r ¼ 0) and also to the radiation condition at infinity (r ! ∞) as given in Eq. (3.2.35). As a result, the following two independent homogeneous solutions exist: z1n ðr Þ ¼ jn ðkr Þ, z2n ðr Þ ¼ hðn1Þ ðkr Þ, and the corresponding Wronskian is given in Eq. (3.2.32b) as

(4.1.34)

Invariant-imbedding T-matrix method

Δ jn ðkr Þ, hðn1Þ ðkr Þ ¼

155

i ðkr Þ2

:

(4.1.35)

Consequently, the radial solution can be explicitly written as ( 0

Rn ðr, r Þ ¼ ik

jn ðkr Þhðn1Þ ðkr 0 Þ; r r 0 ,

(4.1.36)

hðn1Þ ðkr Þjn ðkr 0 Þ; r r 0 :

Then,

!

!

g r , r 0 ¼ ik

∞ X n X n ¼ 0m ¼ n

¼ ik

2

γ 0mn ∗nm ðθ0 , φ0 Þnm ðθ, φÞ

∞ n X 2n + 1 X n¼ 0

4π

( ð1Þ

m ¼ n

m

8 < jn ðkr Þhðn1Þ ðkr 0 Þ; r r 0 : hð1Þ ðkr Þj ðkr 0 Þ; r r 0 n n

Rgψ mn ðkr, θ, φÞψ mn ðkr 0 , θ0 , φ0 Þ; r r 0 , ψ mn ðkr, θ, φÞRgψ mn ðkr 0 , θ0 , φ0 Þ; r r 0 : (4.1.37)

The definition of the dyadic Green function in Eq. (4.1.7) is from the wave equation of the electric field in Eq. (4.1.5). Only in this section do we denote it as the free space $ ! ! electric dyadic Green function G e r , r 0 for distinction. Correspondingly, it is straightforward to define the dyadic Green function associated with the wave equation of the magnetic field as (Tai, 1994) h$ i $ ! ! $ ! ! ! ! r r G m r , r 0 k2 G m r , r 0 ¼ r I δ r r 0 ,

(4.1.38)

$ ! ! where G m r , r 0 is the free space magnetic dyadic Green function. The

explicit expressions for the two dyadic Green functions using the scalar Green function are

$ ! ! $ 1 ! ! G e r , r 0 ¼ I + 2 rr g r , r 0 , k h$ i $ ! ! ! ! Gm r , r 0 ¼ r I g r , r 0 :

(4.1.39a)

(4.1.39b)

After some algebraic steps, the free space magnetic dyadic Green function has this symmetry property: $ T ! ! $ ! ! Gm r 0, r ¼ Gm r , r 0 :

(4.1.40)

156

Invariant Imbedding T-matrix Method

It is evident that the two dyadic Green functions are related by $ ! ! $ ! ! r Ge r , r 0 ¼ Gm r , r 0 ,

(4.1.41a)

$ ! ! $ ! $ ! ! ! r G m r , r 0 I δ r r 0 ¼ k2 G e r , r 0 :

(4.1.41b)

We have two facts to point out in relation to introducing the free space magnetic dyadic Green function in this section. The scalar Green function is expanded using $ ! ! the scalar spherical wave functions, and the explicit expression for G e r , r 0 contains the gradient operator in Eq. (4.1.39a). The longitudinal component of the vector spherical wave function (Rg)Lmn would inevitably appear in the process of expanding $ ! ! G e r , r 0 , even though the final expansion may not explicitly contain the longitudinal component. Moreover, the longitudinal component is associated with the field $ ! ! point in the source region. The curl operator in Eq. (4.1.39b) for G m r , r 0 automatically excludes the longitudinal component. Consequently, the derivation of $ ! ! $ ! ! G m r , r 0 is much easier than the one of G e r , r 0 . A comparison of the degree of complexity using the two methods is shown in Chapter 5 of Tai (1994). The expansion of the free space dyadic Green function using the Ohm-Rayleigh method in Section 3.3.1 contains an isolated singular term in Eq. (3.3.14), which can be given $ ! ! in terms of G m r , r 0 and Eq. (4.1.41b). For convenience, the cases r < r0 and r > r0 are denoted by superscript minus and plus signs, respectively. Therefore, we can rewrite the scalar Green function as follows: ∞ X n X 2n + 1 ! ! ð1Þm + 1 Rgψ mn ðkr, θ, φÞψ mn ðkr 0 , θ0 , φ0 Þ, g r , r 0 ¼ ik 4π n ¼ 0 m ¼ n

(4.1.42a) ∞ X n X 2n + 1 ! ! ð1Þm + 1 ψ mn ðkr, θ, φÞRgψ mn ðkr 0 , θ0 , φ0 Þ, g + r , r 0 ¼ ik 4π n ¼ 0 m ¼ n

(4.1.42b) ! ! ! ! ! ! g r , r 0 ¼ g + r , r 0 U ðr r 0 Þ + g r , r 0 Uðr 0 r Þ,

(4.1.43)

where the step function U is defined as ( 0

U ðr r Þ ¼

1, r > r 0 0, r < r 0

,

(4.1.44a)

Invariant-imbedding T-matrix method

U ðr 0 r Þ ¼

0, r > r 0 1, r < r 0

,

rU ðr r 0 Þ ¼ rU ðr 0 r Þ ¼ r^δðr r 0 Þ:

157

(4.1.44b) (4.1.45)

Using the explicit expression for the free space magnetic dyadic Green function in Eq. (4.1.39b) and Eq. (4.1.43), we have the expression $ ! ! $ ! ! $ ! ! 0 0 G m r , r 0 ¼ G m+ r , r 0 U ðr r 0 Þ + G m r , r U ðr r Þ

$ ! ! ! ! + r^ I δðr r 0 Þ g + r , r 0 g r , r 0 : h$ i $ ! ! ! ! G m r , r 0 ¼ r I g r , r 0 :

(4.1.46) (4.1.47)

The scalar Green function at r ¼ r0 is continuous as a result of the reciprocity shown in Eq. (3.3.2). Consequently, the free space magnetic dyadic Green function in Eq. (4.1.46) is reduced to $ ! ! $ ! ! $ ! ! 0 0 G m r , r 0 ¼ G m+ r , r 0 Uðr r 0 Þ + G m r , r U ðr r Þ:

(4.1.48)

Then, one can have h h $ ! ! $ ! ! i $ ! ! i 0 U ðr 0 r Þ r G m r , r 0 ¼ r G m+ r , r 0 Uðr r 0 Þ + r G m r , r n h$ $ io ! ! ! !0 + δðr r 0 Þ r^ G m+ r , r 0 G : m r , r

(4.1.49)

The tangential discontinuity of the Green function at r ¼ r0 can be stated as (Tai, 1994, Section 4.1 and p. 128) h$ $ i $ ! ! ! !0 r^ G m+ r , r 0 G ¼ I r^r^ δðΩ Ω0 Þ, m r , r

(4.1.50)

where δ(Ω Ω0 ) is a two-dimensional delta function at the spherical surface r ¼ r0 . Using Eq. (4.1.41b) and Eqs. (4.1.49), (4.1.50), the free space electric dyadic Green function can be given as $ ! ! $ ! ! i $ ! ! i 1h 1h 0 G e r , r 0 ¼ 2 r G m+ r , r 0 U ðr r 0 Þ + 2 r G U ðr 0 r Þ m r , r k k ! ! r^r^δ r r 0 , (4.1.51a) k2

158

Invariant Imbedding T-matrix Method

! ! δ r r 0 ¼ δðr r 0 ÞδðΩ Ω0 Þ,

(4.1.51b)

! ! where δ r r 0 is a three-dimensional delta function. Accordingly, the electric dyadic Green function can be obtained from the magnetic dyadic Green function and the isolated singular term. The explicit expression for $ ! ! G m r , r 0 in Eq. (4.1.47) can be written as ! ! $ ! ! ! ðxÞ ! !0 ðyÞ ! !0 ðzÞ ! !0 0 r , r x^ + G r , r y^ + G r , r ^z, G ¼ G m r , r m m m

(4.1.52a)

h i ! ! ! ð cÞ ¼ r g r , r 0 c^ , c^ ¼ x^, y^, ^ z, G m

(4.1.52b)

where x^, y^, and ^ z are the unit vectors in Cartesian coordinates and x^ ¼ sin θ cos φ^ r + cos θ cos φ^ θ sin φ^ φ,

(4.1.53a)

y^ ¼ sin θ sin φ^ r + cos θ sin φ^ θ + cos φ^ φ,

(4.1.53b)

z ¼ cos θ^ ^ r sin θ^ θ:

(4.1.53c)

In addition to the standard vector spherical wave functions, we can define the following ones for convenience in expanding the free space dyadic Green function: r ðΨ mn c^Þ 1 cÞ ¼ r Nðmn , k k

(4.1.54a)

r r ðΨ mn c^Þ 1 cÞ ¼ r Mðmn , k k2

(4.1.54b)

cÞ Mðmn ¼ γ mn

cÞ Nðmn ¼ γ mn

1 cÞ Lðmn ¼ γ 0mn ðrΨ mn Þ c^ ¼ Lmn c^, k 0

(4.1.54c)

where the definitions of γ mn and γ mn are given in Eq. (3.2.39) and the type of the spher(c) (c) ical Bessel function used in M(c) mn, Nmn, or Lmn depends on Ψ mn and can be spherical Bessel functions jn, spherical Neumann functions yn, or spherical Hankel functions of (2) (1) the first kind h(1) n or the second kind hn . For instance, Ψ mn uses hn in Eq. (4.1.54), (c) (c) (c) (c) (c) (c) that is, Ψ mn ¼ ψ mn, then Mmn ¼ Mmn, Nmn ¼ Nmn, and Lmn ¼ Lmn. Using the recurrence relations of the spherical Bessel functions listed in Eq. (3.2.29) and the associated Legendre polynomials listed in Eq. (3.2.24), after some straightforward though tedious algebra, the defined VSWF in Eq. (4.1.54) can be presented using the standard VSWF (Tai, 1994):

Invariant-imbedding T-matrix method

159

xÞ Aðmn i Bm + 1n Bm1n ¼ + ðn + mÞðn m + 1Þ 2nðn + 1Þ γ m + 1n γ mn γ m1n 1 Am + 1n + 1 Am1n + 1 ðn m + 1Þðn m + 2Þ 2ðn + 1Þð2n + 1Þ γ m + 1n + 1 γ m1n + 1 1 Am + 1n1 Am1n1 ðn + m 1Þðn + mÞ + , (4.1.55a) 2nð2n + 1Þ γ m + 1n1 γ m1n1 yÞ Aðmn 1 Bm + 1n Bm1n ¼ ðn + mÞðn m + 1Þ 2nðn + 1Þ γ m + 1n γ mn γ m1n i Am + 1n + 1 Am1n + 1 + ðn m + 1Þðn m + 2Þ + 2ðn + 1Þð2n + 1Þ γ m + 1n + 1 γ m1n + 1 (4.1.55b) i Am + 1n1 Am1n1 + ðn + m 1Þðn + mÞ , 2nð2n + 1Þ γ m + 1n1 γ m1n1 zÞ Aðmn im Bmn 1 n m + 1 Amn + 1 n + m Amn1 , ¼ + + γ mn nðn + 1Þ γ mn 2n + 1 n + 1 γ mn + 1 n γ mn1

(4.1.55c)

where

if A ¼ M, B ¼ N, if A ¼ N, B ¼ M,

(4.1.56)

and the Mmn and Nmn are defined in Eq. (3.2.100). Also (Tai, 1994), xÞ Lðmn 1 ¼ ½Ψ m + 1n + 1 ðn m + 1Þðn m + 2ÞΨ m1n + 1 2ð2n + 1Þ γ 0mn

+

1 ½Ψ m + 1n1 ðn + m 1Þðn + mÞΨ m1n1 , 2ð2n + 1Þ

yÞ Lðmn i ¼ ½Ψ m + 1n + 1 + ðn m + 1Þðn m + 2ÞΨ m1n + 1 0 2ð2n + 1Þ γ mn i ½Ψ m + 1n1 + ðn + m 1Þðn + mÞΨ m1n1 , 2ð2n + 1Þ zÞ Lðmn nm+1 n+m Ψ mn + 1 + Ψ mn1 : ¼ 0 2n + 1 2n + 1 γ mn

(4.1.57a)

(4.1.57b)

(4.1.57c)

In Eqs. (4.1.55), (4.1.57), m can be any integer: positive, negative, or zero. Similar formulas for M(c) and N(c) are given by Tai (1994), but different definitions are used.

160

Invariant Imbedding T-matrix Method

Using r < r0 as an example, the corresponding free space magnetic dyadic Green function can be expressed as (Tai, 1994) ∞ X n cÞ X ! RgMðmn ðkr, θ, φÞ 2n + 1 ðcÞ ! !0 2 G r , r ð1Þm + 1 ψ mn ðkr 0 , θ0 , φ0 Þ ¼ ik m γ 4π mn n ¼ 0 m ¼ n ∞ X n X ð1Þm RgMmn ðkr, θ, φÞ ðcÞ RgNmn ðkr, θ, φÞ ðcÞ 2 αmn + βmn , ¼ ik γ mn γ mn 4π n ¼ 0 m ¼ n

c ¼ x,y, z,

(4.1.58)

where xÞ αðmn ¼

1 0 ψ ðn m 1Þðn mÞψ 0m1n1 2n m + 1n1 1 0 ψ m + 1n + 1 ðn + m + 1Þðn + m + 2Þψ 0m1n + 1 , 2ð n + 1Þ i 0 ψ + ðn m 1Þðn mÞψ 0m1n1 2n m + 1n1 0 i ψ m + 1n + 1 + ðn + m + 1Þðn + m + 2Þψ 0m1n + 1 , + 2ð n + 1Þ

(4.1.59a)

yÞ αðmn ¼

zÞ αðmn ¼

nm 0 n+m+1 0 ψ mn1 + ψ mn + 1 , n n+1

(4.1.59c)

2n + 1 0 ψ m + 1n + ðn + m + 1Þðn mÞψ 0m1n , 2nðn + 1Þ

(4.1.60a)

2n + 1 0 ψ ðn + m + 1Þðn mÞψ 0m1n , 2nðn + 1Þ m + 1n

(4.1.60b)

xÞ βðmn ¼i

yÞ βðmn ¼

(4.1.59b)

zÞ βðmn ¼ i

2n + 1 mψ 0 : nðn + 1Þ mn

(4.1.60c)

The primed arguments of ψ in Eqs. (4.1.59), (4.1.60) are suppressed and are indicated by ψ 0 instead. Using Eq. (4.1.52), the free space magnetic dyadic Green function can be written as $

! !0 G m r , r

∞ X n i X ð1Þm RgMmn ðkr, θ, φÞ h ðxÞ yÞ zÞ αmn x^ + αðmn ^z y^ + αðmn ¼ ik γ mn 4π n ¼ 1 m ¼ n 2

+

i RgNmn ðkr, θ, φÞ h ðxÞ yÞ zÞ βmn x^ + βðmn ^z : y^ + βðmn γ mn

(4.1.61)

Based on the definitions of Eq. (4.1.54) and the standard vector spherical wave functions, we have the following equations:

Invariant-imbedding T-matrix method

161

Mmn MðcÞ ! c^ ¼ mn k r , γ mn γ mn

(4.1.62a)

Nmn NðcÞ ! L ð cÞ c^ ¼ mn k r + 2 0mn , γ mn γ mn γ mn

(4.1.62b)

Mmn ! k r ¼ 0, γ mn

(4.1.63a)

Nmn ! k r ¼ nðn + 1Þψ mn : γ mn

(4.1.63b)

Using the aforementioned equations, we can have xÞ yÞ zÞ βðmn z¼ ^ x^ + βðmn y^ + βðmn

2n + 1 Mmn ðkr 0 , θ0 , φ0 Þ , nð n + 1Þ γ mn

(4.1.64a)

xÞ yÞ zÞ αðmn z¼ ^ x^ + αðmn y^ + αðmn

2n + 1 Nmn ðkr 0 , θ0 , φ0 Þ : nð n + 1Þ γ mn

(4.1.64b)

The free space magnetic dyadic Green function at r < r0 can be rewritten as ∞ X n X $ ! ! 0 G ð1Þm ½RgMmn ðkr, θ, φÞNmn ðkr 0 , θ0 , φ0 Þ ¼ ik2 m r , r n ¼ 1 m ¼ n

+ RgNmn ðkr, θ, φÞMmn ðkr 0 , θ0 , φ0 Þ:

(4.1.65)

The expansion at r > r0 following the same steps can be given as ∞ X n X $ ! ! ð1Þm ½Mmn ðkr, θ, φÞRgNmn ðkr 0 , θ0 , φ0 Þ G m+ r , r 0 ¼ ik2 n ¼ 1 m ¼ n

+ Nmn ðkr, θ, φÞRgMmn ðkr 0 , θ0 , φ0 Þ:

(4.1.66)

Using the relation between the free space electric and magnetic dyadic Green function given in Eq. (4.1.51), the free space electric dyadic Green function can be written as $

!

!

Ge r , r 0

! ! ∞ X n δ r r0 X + ik ð1Þm ¼ ^ r r^ k2 n ¼ 1 m ¼ n

Mmn ðkr, θ, φÞRgMmn ðkr 0 , θ0 , φ0 Þ + Nmn ðkr, θ, φÞRgNmn ðkr 0 , θ0 , φ0 Þ, r > r 0 , RgMmn ðkr, θ, φÞMmn ðkr 0 , θ0 , φ0 Þ + RgNmn ðkr, θ, φÞNmn ðkr 0 , θ0 , φ0 Þ, r < r 0 : (4.1.67)

162

Invariant Imbedding T-matrix Method

This section presents a general technique to obtain the expansion of the free space dyadic Green function. The crucial step is to choose two independent solutions of the spherical Bessel equation as shown in Eq. (3.1.34). It can be any combination, such (2) (1) as (jn(kr), yn(kr)), (h(1) n (kr), hn (kr)), or (jn(kr), hn (kr)). Using the corresponding Wronskian of the chosen two functions as shown in Eq. (4.1.35), the scalar Green function can be obtained as shown in Eq. (4.1.37). Then, the free space magnetic dyadic Green function can be obtained using the scalar Green function as shown in Eqs. (4.1.38)–(4.1.66). Finally, the free space electric Green function can be obtained after some algebra series as shown in Eq. (4.1.67). In this section, we use (jn(kr), h(1) n (kr)) as the basis to give the expansion. Using other independent combinations can be easily realized using the same steps because all the definition and recurrence relations are the same for different spherical Bessel functions.

4.1.4 Vector spherical wave function expansion Several quantities related to the angular functions and the radial functions will be introduced, so the discussion can be performed in a compact way. By analyzing the angular functions in the form of vector spherical wave functions, the following angular quantities are defined as Xmn ðθ, φÞ ¼ ðAmn ðθ, φÞ Bmn ðθ, 0 φÞ Cmn ðθ, φÞÞ 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ n 0 0 nðn + 1Þd0m ðθÞ B C ¼ ð1Þm βn exp ðimφÞ@ iπ mn ðθÞ τmn ðθÞ A, 0 0 τmn ðθÞ iπ mn ðθÞ (4.1.68) where sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2n + 1 , βn ¼ 4πnðn + 1Þ

(4.1.69) 0

1 0 Amn ðθ, φÞ ¼ ð1Þm βn exp ðimφÞ@ iπ mn ðθÞ A, τmn ðθÞ 0

0

(4.1.70a)

1

B C Bmn ðθ, φÞ ¼ ð1Þm βn exp ðimφÞ@ τmn ðθÞ A, iπ mn ðθÞ 0 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 n ðθ Þ nðn + 1Þd0m B C Cmn ðθ, φÞ ¼ ð1Þm βn exp ðimφÞ@ A, 0 0

(4.1.70b)

(4.1.70c)

Invariant-imbedding T-matrix method

163

where the real functions π mn(θ) and τmn(θ) and their symmetric relations are defined in Eq. (3.2.59), the real function dn0m(θ) is defined in Eq. (3.2.79a), and the corresponding symmetric relations are defined in Eq. (3.2.78). Using the symmetric relations, one can obtain A∗mn ¼ ð1Þm Amn ,B∗mn ¼ ð1Þm Bmn ,C∗mn ¼ ð1Þm Cmn ,

(4.1.71a)

X∗mn ¼ ð1Þm Xmn :

(4.1.71b)

Using Eq. (3.2.59) and Eq. (3.3.38) for the functions π mn(θ) and τmn(θ) and using Eq. (3.2.76a) for the function dn0m(θ), their orthogonality relations are 1 2 1 2 1 2

ðπ

sin θdθ½π mn ðθÞπ mn0 ðθÞ + τmn ðθÞτmn0 ðθÞ ¼

0

ðπ

nð n + 1Þ δnn0 , 2n + 1

sin θdθ½π mn ðθÞτmn0 ðθÞ + τmn ðθÞπ mn0 ðθÞ ¼ 0,

(4.1.72a)

(4.1.72b)

0

ðπ 0

0

n n ðθÞd0m ðθ Þ ¼ sin θdθd0m

1 δnn0 : 2n + 1

(4.1.72c)

Moreover, the radial quantities can be defined as 0

1 0 ha ðn; kr Þ 0 hb ðn; kr Þ A, Hn ðkr Þ ¼ @ 0 hc ðn; kr Þ 1 0 ja ðn; kr Þ 0 jb ðn; kr Þ A, Jn ðkr Þ ¼ @ 0 jc ðn; kr Þ

(4.1.73a)

0

(4.1.73b)

where 1 hð1Þ ðkr Þ n ð1Þ B 1 d krh ðkr Þ C ha ðn; kr Þ n B C C, @ hb ðn; kr Þ A ¼ B kr dðkr Þ B C @ ð 1 Þ hc ðn; kr Þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ hn ðkr Þ A nð n + 1Þ kr 0 1 jn ðkr Þ 0 1 ja ðn; kr Þ B 1 d½krjn ðkr Þ C C @ jb ðn; kr Þ A ¼ B B kr dðkr Þ C: @ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ A jc ðn; kr Þ jn ðkr Þ nð n + 1Þ kr 0

1

0

(4.1.74a)

(4.1.74b)

164

Invariant Imbedding T-matrix Method

The regular and outgoing vector spherical wave functions can be written in a compact way as follows: ð RgMmn ðkr, θ, φÞ RgNmn ðkr, θ, φÞ Þ ¼ ð Amn ðθ, φÞja ðn; kr Þ Bmn ðθ, φÞjb ðn; kr Þ + Cmn ðθ, φÞjc ðn; kr Þ Þ

(4.1.75a)

¼ Xmn ðθ, φÞJn ðkr Þ, ð Mmn ðkr, θ, φÞ Nmn ðkr, θ, φÞ Þ ¼ ð Amn ðθ, φÞha ðn; kr Þ Bmn ðθ, φÞhb ðn; kr Þ + Cmn ðθ, φÞhc ðn; kr Þ Þ

(4.1.75b)

¼ Xmn ðθ, φÞHn ðkr Þ: The free space Dyadic green function in Eq. (4.1.67) can also be written in a compact way: ! !0 δ r r $ ! ! $ ! ! 0 G r , r 0 ¼ ^ + G r r^ 0 r , r , k2

(4.1.76a)

8 1 0 > RgMmn kr 0 , θ0 , φ0 > > > 0 @ A > Mmn ðkr, θ, φÞ Nmn ðkr, θ, φÞ > , r>r , > > < RgN kr0 , θ0 , φ0

∞ n X X $ ! !0 G 0 r , r ¼ ik ð1Þm > > > n ¼ 1 m ¼ n >

> > > > :

mn

1 0 Mmn kr 0 , θ0 , φ0 0 A RgMmn ðkr, θ, φÞ RgNmn ðkr, θ, φÞ @ , r
ð4:1:76bÞ $ ! ! $ T ! ! $ ! ! The expansion for G 0 r , r 0 has the symmetric property, G 0 r 0 , r ¼ G 0 r , r 0 ,

as stated in Eq. (3.3.5) and a result of reciprocity (Tai, 1994; Tsang et al., 2000). Using Eqs. (4.1.68), (4.1.75), the free space Dyadic green function can be written in a simple form ∞ X n X $ ! ! Xmn ðθ, φÞ Gn ðr, r 0 Þ X{mn ðθ0 , φ0 Þ, G0 r , r 0 ¼

(4.1.77)

n ¼ 1 m ¼ n

where the superscript “†” represents the matrix transpose and complex conjugate, and ( 0

Gn ðr, r Þ ¼ ik

Hn ðkr ÞJTn ðkr 0 Þ, r > r 0 , Jn ðkr ÞHTn ðkr 0 Þ, r < r:

Gn(r, r0 ) is discontinuous at r ¼ r0 and GTn (r0 , r) ¼ Gn(r, r0 ).

(4.1.78)

Invariant-imbedding T-matrix method

4.2

165

Concept of the invariant-imbedding T-matrix method

Spherical coordinates are used in this book to obtain the T-matrix method. Accordingly, the imbedding process is performed in the radial coordinate in the form of infinitesimal spherical shells. Consequently, we need to integrate the angular coordinates and find the invariant-imbedding formula associated with the radial coordinate.

4.2.1 Radial integral Using the volume integral equation associated with the electric field in (4.1.12a) and the free space dyadic Green function (4.1.76a), the volume integral equation can be transformed as 3 ! ð $ u r r^r^ ! ! ! ! ! inc ! 4I + 5E ! r r + u r 0 G 0 r , r 0 E r 0 d3 r 0 : ¼ E 2 k 2

$

(4.2.1)

V1

The singular term in the Green function is considered in Eq. (4.2.1). Accordingly, we can define the electric field excluding the contribution from the source region as !

E0 r

and

3 2 ! h i1 u r r^r^ $ ! ! ! 4 5E ! r ¼ Z r E r : E0 r ¼ I + k2

(4.2.2)

The coefficient Z is written in matrix form in Eq. (4.2.2) using spherical coordinates and can be given as 0 1 ! ε=ε1 r 0 0 ! B C Z r ¼@ 0 1 0 A: 0 0 1

(4.2.3)

! ! It is evident that E0 r E r , if r > rc, where rc is the radius of the smallest circumscribedsphere of the scattering particleas shown in Fig. 4.2. Replacing the elec!

tric field E r

!

with the electric field E0 r

in Eq. (4.2.1), the volume integral

equation can be rearranged as ð $ ! ! ! ! ! ! ! ! E0 r ¼ Einc r + u r 0 G 0 r , r 0 Z r 0 E0 r 0 d3 r 0 : V1

(4.2.4)

166

Invariant Imbedding T-matrix Method

Fig. 4.2 Schematic figure for the invariant-imbedding method. In panel (A), the shaded region is the dielectric scattering particle, and ri and rc are the radii of the largest inscribed and the smallest circumscribed spheres of the particle. In panel (B), the region where the radius is between ri and rc can be divided into spherical shells, and light scattering can be computed using the invariant-imbedding technique. rn1 and rn are the radii of layers (n 1) and n, respectively.

In Eq. (3.2.51a), the incident field can be expanded using the regular vector spherical wave functions. Substituting Eq. (4.1.77) into Eq. (4.2.4) and using a single incident component ð RgMmn RgNmn Þ, we can rewrite the aforementioned equation as

! E0 r

m 0 n0

¼ Xm0 n0 ðΩÞJn0 ðkr Þ +

∞ X n X

Xmn ðΩÞ

n ¼ 1m ¼ n

2

ð

ð rc (

fdr 0 r 02 Gn ðr, r 0 Þ

0

! ! ! 4 dΩ u r 0 X{mn ðΩ0 Þ Z r 0 E0 r 0 0

m0 n0

3) 5 ,

(4.2.5)

4π

where the angular arguments (θ, φ) are replaced by the argument Ω for convenience. We can organize all the angular-related functions into a new function: ð h ! ! ! Fmnm0 n0 ðr Þ ¼ r 2 dΩ u r X{mn ðΩÞ Z r E0 r

i

m 0 n0

:

(4.2.6)

4π

Accordingly, Eq. (4.2.5) can be rewritten as

! E0 r

m0 n0

¼ Xm0 n0 ðΩÞJn0 ðkr Þ +

∞ X n X n ¼ 1 m ¼ n

Xmn ðΩÞ

ð rc

dr 0 Gn ðr, r 0 Þ Fmnm0 n0 ðr 0 Þ:

0

(4.2.7)

Invariant-imbedding T-matrix method

167

When r > rc, this satisfies the relation r > r0 . Substituting Eq. (4.1.78) into Eq. (4.2.7), we can obtain ! E0 r

m 0 n0

¼X

m0 n0

ðΩÞJ ðkr Þ +

∞ X n X

n0

n ¼ 1 m ¼ n

Xmn ðΩÞ Hn ðkr Þik

ð rc 0

dr 0 JTn ðkr 0 Þ Fmnm0 n0 ðr 0 Þ: (4.2.8)

We can only consider a single incident component ð RgMmn RgNmn Þ in Eq. (4.2.5). Using Eq. (4.1.75b), Xmn(Ω)Hn(kr) in the second term of Eq. (4.2.8) means a single component associated with the outgoing vector spherical wave functions. The physical meaning of Eq. (4.2.8) is evident. The electric field in Eq. (4.2.8) is composed of two terms, specifically, the incident field and the scattered field. According to the definition of the T-matrix in Eq. (3.2.53), the T-matrix can be represented as the following integral: Tmnm0 n0 ¼ ik

ð rc 0

dr 0 JTn ðkr 0 Þ Fmnm0 n0 ðr 0 Þ:

(4.2.9)

Eq. (4.2.8) is physically simplified to the following form:

! E0 r

m0 n0

¼ Xm0 n0 ðΩÞJn0 ðkr Þ +

∞ X n X

Xmn ðΩÞ Hn ðkr ÞTmnm0 n0 :

(4.2.10)

n ¼ 1 m ¼ n

The electric field inside the smallest circumscribed sphere of the scattering particle is still necessary according to Eqs. (4.2.6), (4.2.9). Using Eqs. (4.2.5), (4.2.6), the angular integral function Fmnm0 n0 has the following relation: Fmnm0 n0 ðrÞ ¼ Umnm0 n0 ðrÞJn0 ðkr Þ +

∞ X

n0 0 X

n00 ¼ 1 m0 0 ¼ n0 0

Umnm0 0 n0 0 ðrÞ

ð rc

dr 0 Gn0 0 ðr, r0 Þ Fm0 0 n0 0 m0 n0 ðr0 Þ,

0

(4.2.11)

where ð ! Umnm0 n0 ðr Þ ¼ r 2 dΩu r X{mn ðΩÞ Zðr, ΩÞ Xm0 n0 ðΩÞ:

(4.2.12)

Eq. (4.2.9) and Eqs. (4.2.11), (4.2.12) are the key equations to obtain the T-matrix of a dielectric scattering particle using the invariant-imbedding T-matrix method. Eqs. (4.2.5)–(4.2.12) earlier are given in component form for clarity. For succinctness, we next use matrix notation to replace all the equations and quantities. As stated in Section 3.2.3, the indices (m,n) are replaced by one index l. Once the truncation term N associated with the vector spherical wave functions is determined, the matrix or vector dimension L associated with the index l is given in Eq. (3.2.54). Accordingly, matrices J and H have dimension 3L 2L and are block diagonal; matrix F has

168

Invariant Imbedding T-matrix Method

dimension 3L 2L; matrices U and G have dimension 3L 3L, and G is block diagonal; and the dimension of the T-matrix T is 2L 2L. Eqs. (4.2.9), (4.2.11) can be written in compact matrix forms as T ¼ ik

ð rc

dr 0 JT ðkr 0 ÞFðr 0 Þ,

(4.2.13)

0

Fðr Þ ¼ Uðr ÞJðkr Þ + Uðr Þ

ð rc

dr 0 Gðr, r 0 ÞFðr 0 Þ:

(4.2.14)

0

4.2.2 Invariant-imbedding method: Differential form and difference form The invariant-imbedding method is comprehensively introduced in general formulas and applications by Bellman and Wing (1975). Here, we only present the invariantimbedding technique according to the T-matrix computations given in Eqs. (4.2.13), (4.2.14) in our convenience. The upper limit rc in Eqs. (4.2.13), (4.2.14) is a constant. If we consider light scattering by a particle whose smallest circumscribed sphere radius is rn, the upper limit now is a variable, 0 rn rc. The integral forms in Eqs. (4.2.13), (4.2.14) remain unchanged with the increase of the radial length, which is called the principle of invariance. We can imbed the radial length from rn to rn + Δr and find the relation between T(rn + Δr) and T(rn) or the expression for dTdrðrnn Þ. The imbedding method is similar to a perturbation method so that we can obtain the perturbation equation dTdrðrnn Þ. The final T-matrix can be obtained by iteratively imbedding the radial length rn from 0 to rc, as shown in Fig. 4.3. The invariant-imbedding method transforms a boundary condition problem to an initial condition problem. In the

Fig. 4.3 Invariant-imbedding steps. The invariant-imbedding algorithm starts from the inscribed sphere of a scattering particle and keeps imbedding spherical shells until the whole particle is included in the imbedded region.

Invariant-imbedding T-matrix method

169

following part, we will give two invariant-imbedding formulas for the T-matrix: the derivative form dTdrðrnn Þ, which is derived from the integral equation, and the difference form for T(rn + Δr) and T(rn), which is derived from the quadrature expressions associated with the integral equations. Moreover, we will prove that the difference form is reduced to the derivative form when Δr ! 0. For clarity, we rewrite Eqs. (4.2.13), (4.2.14) and change the upper limit from the constant rc to an intermediate variable rn as Tðrn Þ ¼ ik

ð rn

dr 0 JT ðkr 0 ÞFðrn j r 0 Þ,

(4.2.15)

0

Fðrn j r Þ ¼ Uðr ÞJðkr Þ + Uðr Þ

ð rn

dr 0 Gðr, r 0 ÞFðrn j r 0 Þ,

(4.2.16)

0

where the rn dependence in the matrix F is explicitly shown. It is apparent that Eq. (4.2.16) is a Fredholm integral equation of the second kind, where the forcing term U(r)J(r) and the integral kernel G(r, r0 ) are known. A general invariant-imbedding algorithm for the Fredholm integral equation of the second kind is given in Section 12.9 by Bellman and Wing (1975). Since our integral kernel in Eq. (4.1.78) and the forcing term are explicitly known, the invariant-imbedding form for the T-matrix is presented in the following part. If we differentiate Eq. (4.2.16) with respect to the variable rn, we can obtain ∂Fðrn j r Þ ¼ Uðr ÞJðkr ÞikHT ðkr n ÞFðrn j rn Þ + Uðr Þ ∂rn

ð rn

dr 0 Gðr, r 0 Þ

0

∂Fðrn j r 0 Þ , ∂rn (4.2.17)

where the condition r < rn for G(r, rn) is used in the first term of the right-hand side. Comparing the forcing term and the integral kernel in Eqs. (4.2.16), (4.2.17), it is evident that ∂Fðrn j r Þ ¼ Fðrn j r ÞikHT ðkr n ÞFðrn j rn Þ: ∂rn

(4.2.18)

Using Eq. (4.2.16), we can obtain Fðrn j rn Þ ¼ Uðrn Þ½Jðkr n Þ + Hðkr n ÞTðrn Þ,

(4.2.19)

where the condition rn > r0 for G(rn, r0 ) is used in the second term of the right-hand side and Eq. (4.2.15) is also used. Differentiating Eq. (4.2.15) with respect to the variable rn, we can obtain dTðrn Þ ¼ ikJT ðkr n ÞFðrn j rn Þ + ik dr n

ð rn 0

dr 0 JT ðkr 0 Þ

∂Fðrn j r 0 Þ : ∂rn

(4.2.20)

170

Invariant Imbedding T-matrix Method

Using Eq. (4.2.15), Eqs. (4.2.18), (4.2.19), Eq. (4.2.20) can be presented as dTðrn Þ ¼ ik JT ðkr n Þ + Tðrn ÞHT ðkr n Þ Uðrn Þ½Jðkr n Þ + Hðkr n ÞTðrn Þ: dr n

(4.2.21)

Eq. (4.2.21) is the invariant-imbedding equation for the T-matrix in a derivative form ( Johnson, 1988). A similar treatment with respect to the derivative form is also given by Doicu and Wriedt (2018). Before obtaining the invariant-imbedding equation for the T-matrix in a difference form, we first need to discretize the integrations into quadrature expressions for Eqs. (4.2.13), (4.2.14) as shown in (e.g., Press et al., 2007, p. 989) T ¼ ik

n0 X

ωj JT kr j F rj ,

(4.2.22)

j¼ 1

Fðri Þ ¼ Uðri ÞJðkr i Þ + Uðri Þ

n0 X

ωj G ri , rj F rj , i ¼ 1, …,n0 ,

(4.2.23)

j¼ 1

where n0 is the total order of discretization, or total layers, {rj} are the abscissas, and {ωj} are the corresponding weights of a quadrature rule. Eq. (4.2.23) is actually a set of linear algebraic equations associated with function F(ri). The solution can be algebraically obtained by solving the linear equations because the algebraic matrix for the Fredholm equation of the second kind is usually well conditioned (Press et al., 2007, p. 990). However, the dimension of the matrix or vector is too large to directly compute since Eq. (4.2.23) itself is a matrix equation. Alternatively, the imbedding form is employed to iteratively obtain the solution. Because the imbedding form is derived from the Fredholm equation of the second kind, the matrix expression in the imbedding form still remains well conditioned. When the layer number is n, where 1 n n0, the corresponding T-matrix and radial functions are the same as the ones in Eqs. (4.2.22), (4.2.23) except that the upper limits of the summations have to be n instead of n0, due to the principle of invariance. We can rewrite Eqs. (4.2.22), (4.2.23) for a changing layer number n and explicitly present the layer number dependence as follows: TðnÞ ¼ ik

n X

ωj JT kr j F nj rj ,

(4.2.24)

j¼ 1

Fðnj ri Þ ¼ Uðri ÞJðkr i Þ + Uðri Þ

n X

ωj G ri , rj F nj rj , i ¼ 1,…,n,

(4.2.25)

j¼ 1

where T(n) ¼ T(rn). The difference form of the invariant-imbedding method for the T-matrix was proposed by Johnson (1988). Using similar steps to the derivation of the derivative form, the difference form is given in the following part.

Invariant-imbedding T-matrix method

171

Using Eqs. (4.1.78), (4.2.25), we can get the expression for F(n j rn) when i ¼ n as follows: " Fðnj rn Þ ¼

ω1 n Qðrn Þ

Jðkr n Þ + Hðkr n Þik

n1 X

ωj J ðkr n ÞF nj rj T

# ,

(4.2.26)

j¼ 1

where Qðrn Þ ¼ ωn ½I ωn Uðkr n ÞGðrn , rn Þ1 Uðrn Þ:

(4.2.27)

Then, we can get the expression for F(nj ri) when 1 i < n as follows: Fðnj ri Þ ¼ Uðri ÞJðkr i Þ + Uðri Þ

n1 X

ωj G ri , rj F nj rj + Uðri ÞJðkr i Þωn ikHT ðkr n ÞFðnj rn Þ,

j¼ 1

(4.2.28a)

or, n1 X ωj G ri , rj F nj rj : Fðnj ri Þ ¼ Uðri ÞJðkr i Þ I + ωn ikHT ðkr n ÞFðnj rn Þ + Uðri Þ j¼ 1

(4.2.28b) Comparing the forcing term in Eq. (4.2.28b) and the forcing term in Eq. (4.2.25) and also using the expression for F(n 1j ri) in Eq. (4.2.25), it is evident that we can define the following relation: Fðnj ri Þ ¼ Fðn 1j ri ÞðI + pn Þ,

(4.2.29)

where pn ¼ ωn ikHT ðkr n ÞFðnj rn Þ:

(4.2.30)

Using Eqs. (4.2.24), (4.2.29), Eq. (4.2.26) can be rewritten as Fðnj rn Þ ¼ ω1 n Qðrn Þ½Jðkr n Þ + Hðkr n ÞTðn 1ÞðI + pn Þ,

(4.2.31)

Using Eqs. (4.2.30), (4.2.31), we can explicitly write down pn ¼ Qhj ðrn Þ + Qhh ðrn ÞTðn 1ÞðI + pn Þ,

(4.2.32a)

I + pn ¼ ½I Qhh ðrn ÞTðn 1Þ1 I + Qhj ðrn Þ ,

(4.2.32b)

or,

172

Invariant Imbedding T-matrix Method

where Qhj ¼ ikHT QJ,

(4.2.33)

Qhh ¼ ikHT QH:

(4.2.34)

Accordingly, the T-matrix for n layers in Eq. (4.2.24) can be written as TðnÞ ¼ ikωn JT ðkr n ÞFðnj rn Þ + Tðn 1Þ½I + pn :

(4.2.35)

Using Eqs. (4.2.31), (4.2.32), Eq. (4.2.35) can be written in an iterative form as TðnÞ ¼ Qjj ðrn Þ + I + Qjh ðrn Þ Tðn 1Þ½I Qhh ðrn ÞTðn 1Þ1 I + Qhj ðrn Þ , (4.2.36) where Qjj ¼ ikJT QJ,

(4.2.37)

Qjh ¼ ikJT QH:

(4.2.38)

Eq. (4.2.36) is the invariant-imbedding T-matrix method in difference form. The difference form directly gives a recurrence relation, which is computationally efficient, so the difference form is used in the following T-matrix computation. Since all Q-related functions are only related to the quadrature point rn, they can be treated as the contributions of the reflection and transmission associated with the corresponding single imbedding layer. Qjj(rn) represents the exterior to exterior reflection of the single layer and Qhh(rn) the interior to interior reflection of the single layer; I + Qhj(rn) represents the transmission from the exterior to interior direction and I + Qjh(rn) the transmission of the single layer from the interior to exterior direction. Moreover, T(n 1) represents the reflection from all the (n-1) layers. The matrix inversion in Eq. (4.2.36) can be decomposed into summations using a Taylor expansion as ½I Qhh ðrn ÞTðn 1Þ1 ¼ I +

∞ X

½Qhh ðrn ÞTðn 1Þp :

(4.2.39)

p¼ 1

The physical interpretation of the difference form can be given in terms of the scattering orders between the single layer and the interior layers as follows: First-order scattering: Qjj(rn) Second-order scattering: [I + Qjh(rn)]T(n 1)[I + Qhj(rn)] Third-order scattering: [I + Qjh(rn)]T(n 1)Qhh(rn)T(n 1)[I + Qhj(rn)] … pth-order scattering: [I + Qjh(rn)]T(n 1)[Qhh(rn)T(n 1)]p[I + Qhj(rn)] …

Invariant-imbedding T-matrix method

173

(A)

(B) Fig. 4.4 Physical interpretation of the invariant-imbedding expression of the T-matrix. The left panel of (A) shows reflection and transmission for a single layer for exterior and interior incidence, and the right panel shows the reflection for the layers inside the single layer. (B) A schematic of the physical interpretation. The dash-dot arrows represent the incident light.

which are also shown in Fig. 4.4. Now, we prove that the difference form is the same as the derivative form in the limit Δr ! 0. We use the layer depth between layer n-1 and layer n as the weight, that is, ωn ¼ Δr ¼ rn rn1 :

(4.2.40)

174

Invariant Imbedding T-matrix Method

Only retaining terms in which are first order in Δr, we can obtain the following expressions: Qðrn Þ ΔrUðrn Þ,

(4.2.41)

Qhj ikΔrHT UJ, Qhh ikΔrHT UH,

(4.2.42a)

Qjj ikΔrJT UJ, Qjh ikΔrJT UH,

(4.2.42b)

½I Qhh ðrn ÞTðn 1Þ1 I + ikΔrHT ðrn ÞUðrn ÞHðrn ÞTðn 1Þ:

(4.2.43)

Substituting Eqs. (4.2.41)–(4.2.43) into Eq. (4.2.36) and retaining first-order terms in Δr, the invariant-imbedding equation in the derivative form can be obtained in the limit Δr ! 0 as follows ( Johnson, 1988): dTðnÞ TðnÞ Tðn 1Þ ¼ lim Δr!0 dr Δr ¼ ik JT ðrn Þ + Tðrn ÞHT ðrn Þ Uðrn Þ½Jðrn Þ + Hðrn ÞTðrn Þ:

(4.2.44)

4.2.3 T-matrix computation As stated in Eq. (4.1.78), function G(r, r0 ) is discontinuous when r approaches r0 from above or below. The integration in Eq. (4.2.14) is replaced by a quadrature expression in Eq. (4.2.23). The quadrature method should be carefully chosen because high-order quadrature methods may give worse results than low-order quadrature methods when applied to an integrand with a discontinuity (Farrington, 1961), here caused by G(r, r0 ). Therefore, the trapezoidal rule is used in computation to give the quadrature expression in Eq. (4.2.23). Accordingly, the discontinuous function G(r, r0 ) at r ¼ r0 can be stated as the average of two limiting values approaching r0 from above and below: Gn ðr, r Þ ¼

ik Hn ðkr ÞJTn ðkr Þ + Jn ðkr ÞHTn ðkr Þ : 2

(4.2.45)

In the imbedding process, a recurrence relation with respect to the T-matrix is given in Eq. (4.2.36). As stated in Section 3.3.3, the T-matrix of a homogeneous spherical particle can be analytically obtained using the Lorenz-Mie theory, and the coefficients are given in Eq. (3.3.41). To optimally save on computational time, the initial value for Eq. (4.2.36) is taken to be the T-matrix at r ¼ ri as shown in Fig. 4.3, and the imbedding or iterative region is from ri to rc. Eq. (4.2.36) is the equation for iteratively computing the T-matrix of an arbitrary particle. As shown for the physical interpretation in Fig. 4.4, the T-matrix T(n) associated with n layers is determined by the T-matrix T(n 1) associated with (n 1) layers and Qjj, Qjh, Qhj, and Qhh, which only depend on the characteristics of the nth layer. Radial functions H(krn), J(krn), and G(rn, rn) only depend on the radius

Invariant-imbedding T-matrix method

175

of the nth layer and not on any particle properties, such as shape and refractive index. Accordingly, the only particle-related quantity is the matrix U(rn), and using its definition in Eq. (4.2.12), it can be explicitly expressed as sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 0 ð2n + 1Þð2n0 + 1Þ Umnm0 n0 ðr Þ ¼ k2 r 2 ð1Þm + m nðn + 1Þn0 ðn0 + 1Þ

where

1 4π

ðπ 0

sin θdθ

ð 2π 0

ε1 ðr, ΩÞ 1 umnm0 n0 ðr, ΩÞ , dφ exp ½iðm0 mÞφ ε (4.2.46)

0

1 i½π mn ðθÞτm0 n0 ðθÞ π mn ðθÞπ m0 n0 ðθÞ 0 B + τmn ðθÞτm0 n0 ðθÞ C + τmn ðθÞπ m0 n0 ðθÞ B C B C B i½π ðθÞτ 0 0 ðθÞ C π mn ðθÞπ m0 n0 ðθÞ B mn C mn 0 umnm0 n0 ðr, ΩÞ ¼ B +τ ðθÞπ 0 0 ðθÞ C, 0 0 + τ ð θ Þτ ð θ Þ mn mn B C mn mn B C pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ n 0 B C n 0 0 @ nðn + 1Þn ðn + 1Þd0m ðθÞd0m0 ðθÞ A 0 0 ½ε1 ðr, ΩÞ=ε

ð4:2:47Þ where ε1(r, Ω) and ε are the permittivity of the scattering particle and the surrounding medium, respectively, and the refractive index m of the scattering particle relative to the medium has the relation m2(r, Ω) ¼ ε1(r, Ω)/ε. All symmetries associated with the T-matrix caused by the geometric symmetries of the scattering particle introduced in Section 3.4.2 depend on the symmetries of matrix U because the only shape-related quantities in the recurrence relation of Eq. (4.2.26) are all involved in matrix U. Moreover, the permittivity ε1(r, Ω) in matrix U is the only shape-related quantity. The application of matrix symmetries can reduce the computational burdens. Especially, axially or discrete rotational symmetries can decouple the T-matrix, so the computation can be significantly reduced. We will discuss the expressions of matrix U under some shape symmetries. The argument θ for functions π mn, τmn, and dn0m and the arguments (r, Ω) or (r, θ, φ) for the permittivity ε1 in the following expressions are suppressed for clarity. If the scattering particle has horizontal mirror symmetry, the particle permittivity is described as ε1 ðr, π θ, φÞ ¼ ε1 ðr, θ, φÞ:

(4.2.48)

Using the symmetric relations in Eq. (4.1.72), matrix U can be written as sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð2n + 1Þð2n0 + 1Þ m + m0 2 2 Umnm0 n0 ðrÞ ¼ k r ð1Þ nðn + 1Þn0 ðn0 + 1Þ ð 2π n ðπ o ε 1 2 1 sinθdθ dφ exp ½iðm0 mÞφ 1 e umnm0 n0 ðr, ΩÞ , ε 4π 0 0 (4.2.49)

176

Invariant Imbedding T-matrix Method

0

e umnm0 n0

cmnm0 n0 ie cmnm0 n0 0 B ½π mn π m0 n0 + τmn τm0 n0 ½ π mn τ m0 n0 + τmn π m0 n0 B B B ie cmnm0 n0 B cmnm0 n0 0 B ½π mn π m0 n0 + τmn τm0 n0 ¼ B ½π mn τm0 n0 + τmn π m0 n0 B B cmnm0 n0 B pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ n n0 B @ nðn + 1Þn0 ðn0 + 1Þd0m d0m0 0 0 ½ε1 =ε

1 C C C C C C C, C C C C A

(4.2.50)

where cmnm0 n0 ¼ 1 + ð1Þn + m + n

0

cemnm0 n0 ¼ 1 ð1Þn + m + n

+ m0 0

,

+ m0

(4.2.51a)

:

(4.2.51b)

If the scattering particle has vertical mirror symmetry, we generally assume the line of symmetry to be along the direction specified by azimuthal angle φ0 (or φ0 + π), and the particle permittivity can be described as ε1 ðr, θ, 2π + 2φ0 φÞ ¼ ε1 ðr, θ, φÞ:

(4.2.52)

Matrix U can be expressed as sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð2n + 1Þð2n0 + 1Þ exp ½iðm0 mÞφ0 Umnm0 n0 ðr Þ ¼ k2 r2 ð1Þ nðn + 1Þn0 ðn0 + 1Þ

ð φ0 + π

ð 1 π ε1 ðr, ΩÞ sin θdθ dφ cos ½iðm0 mÞðφ φ0 Þ 1 umnm0 n0 ðr, ΩÞ , 2π 0 ε φ0 m + m0

(4.2.53)

where umnm0 n0 is given in Eq. (4.2.47). If the scattering particle has N-fold (N is a finite integer) rotational symmetry, the particle permittivity satisfies the relation ε1 ðr, θ, φ + j2π=N Þ ¼ ε1 ðr, θ, φÞ, j ¼ 1, 2, …,N 1:

(4.2.54)

Matrix U can be partially decoupled as sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð2n + 1Þð2n0 + 1Þ Umnm0 n0 ðr Þ ¼ δðm + lNÞm0 k r ð1Þ nðn + 1Þn0 ðn0 + 1Þ 2 2

N 4π

ðπ 0

lN

ð 2π n ε o 1 1 umnm0 n0 , sin θdθ N dφ exp ½ilNφ ε 0

(4.2.55)

Invariant-imbedding T-matrix method

177

where l is an integer, that is, if mod(m 0 m, N) ¼ 0 (mod is the modular arithmetic function), δ(m+lN)m0 ¼ 1; lN ¼ l ∗ N. Accordingly, matrix U can be decomposed into N components, and the components can be formally described as

ðUmnm0 n0 Þj , j ¼ 0,1, 2, …,N 1

mod ðm0 m, N Þ ¼ 0, mod ðm, N Þ ¼ j:

(4.2.56)

If the scattering particle has axial symmetry, the particle permittivity is azimuthally independent. The matrix U can be decoupled with component m as sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð

ð2n + 1Þð2n0 + 1Þ 1 π ε1 ðr, θÞ 2 2 1 umnmn0 ðr, θÞ , Umnm0 n0 ðr Þ ¼ δmm0 k r sin θdθ nðn + 1Þn0 ðn0 + 1Þ 2 0 ε (4.2.57) where 1 π mn ðθÞπ mn0 ðθÞ i½π mn ðθÞτmn0 ðθÞ 0 C B + τmn ðθÞτmn0 ðθÞ + τmn ðθÞπ mn0 ðθÞ C B C B C B i½π mn ðθÞτmn0 ðθÞ π mn ðθÞπ mn0 ðθÞ C: B 0 umnmn0 ðr, θÞ ¼ B + τ ðθÞπ 0 ðθÞ C + τmn ðθÞτmn0 ðθÞ mn mn C B pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ n C B 0 n 0 0 @ nðn + 1Þn ðn + 1Þd0m ðθÞd0m ðθÞ A 0 0 ½ε1 ðr, θÞ=ε 0

ð4:2:58Þ Accordingly, the matrix U can be decomposed into an infinite number of components as ðUmnmn0 Þm , m ¼ 0, 1, 2, …

(4.2.59)

For a particle with spherical symmetry, the particle permittivity is independent of the angular variables. Using the symmetry given in Eq. (4.1.72), matrix U is completely decoupled as Umnm0 n0 ðr Þ ¼ δmm0 δnn0 k2 r 2

ε ε1 ðr Þ 1 diag 1 1 : ε1 ð r Þ ε

(4.2.60)

There are similar discussions for matrix U in Bi et al. (2013). The layer-related matrices H(krn), J(krn), and G(rn, rn) are radial functions and block diagonal. Accordingly, if matrix U is decoupled such as when it is associated with a particle with N-fold rotational symmetry or axial symmetry, the T-matrix also has the same decoupled property. The symmetric property of the T-matrix associated with the particle shape is comprehensively discussed in Section 3.4.2.

178

Invariant Imbedding T-matrix Method

We have discussed the decoupled T-matrix here according to the symmetry of matrix U. If the particle has N-fold rotational symmetry, the recurrence relation for the T-matrix can be decomposed into N components as h i ðTðnÞÞj ¼ Qjj ðrn Þ j + I + Qjh ðrn Þ j ðTðn 1ÞÞj h i1 h i I ðQhh ðrn ÞÞj ðTðn 1ÞÞj I + Qhj ðrn Þ j ,

(4.2.61)

where the jth component is defined as mod(m m0 , N) ¼ 0 and j ¼ mod(m, N) in Eq. (4.2.56). As stated in Section 3.4.2, only components for j ¼ 0, 1, 2, …, int[N/2] are independent, and the remaining components can be obtained using the T-matrix symmetry. If the particle has axial symmetry, the recurrence relation for the T-matrix can be decomposed into infinite components according to index m as h i ðTðnÞÞm ¼ Qjj ðrn Þ m + I + Qjh ðrn Þ m ðTðn 1ÞÞm i 1 h I ðQhh ðrn ÞÞm ðTðn 1ÞÞm I + Qhj ðrn Þ m ,

(4.2.62)

where component m is defined as the one in Eq. (4.2.57). Only components for m ¼ 0, 1, 2, … are independent, and the remaining components can be obtained using the T-matrix symmetry. If the particle has spherical symmetry, the recurrence relation for the T-matrix is completely decoupled as h i ðTðnÞÞn ¼ Qjj ðrn Þ n + I + Qjh ðrn Þ n ðTðn 1ÞÞn 1 h i I ðQhh ðrn ÞÞn ðTðn 1ÞÞn I + Qhj ðrn Þ n ,

(4.2.63)

where all quantities are m independent and n ¼ 1, 2, …; all the Q and the matrices

11 Tn ð nÞ 0 T-matrix are 2 2 diagonal matrices, for instance, ðTðnÞÞn ¼ ; 0 T22 n ð nÞ and the analytical expressions are given in Eq. (3.3.41).

4.2.4 Spherical validation of differential form using Lorenz-Mie coefficients Similarly to Eq. (4.2.63), the derivative form of the IITM with respect to a spherical particle given in Eq. (4.2.21) can be explicitly given by 2 dT11 n ðr Þ ¼ ikðkr Þ2 m2 1 ja + T11 n ðr Þha , dr

(4.2.64a)

2 2 2 i h dT22 n ðr Þ ¼ ikðkr Þ2 m2 1 jb + T22 + jc + T22 n ðr Þhb n ðr Þhc =m , dr

(4.2.64b)

Invariant-imbedding T-matrix method

179

where k is the wave number in the surrounding medium; m here is the refractive index relative to the surrounding medium, that is, m2 ¼ ε1(r)/ε; n is the index to denote the expansion of VSWF; and the functions j and h are defined in terms of the spherical Bessel function and the spherical Hankel function of the first kind as in Eq. (4.1.74), respectively. For a homogeneous spherical particle with relative refractive index m and radius r, the Lorenz-Mie coefficients are given in Eq. (3.3.40), and the corresponding RiccatiBessel functions are given in Eq. (3.3.41). We use the conventional notation x ¼ kr for spherical particles. After some algebra, the derivatives of the Lorenz-Mie coefficients can be given by dbn ðxÞ ¼ i 1 m2 ðζ n ðxÞ bn ðxÞξn ðxÞÞ2 : (4.2.65a) dx 2 2 dan ðxÞ ¼ i 1 m2 ζ 0n ðxÞ an ðxÞξ0n ðxÞ + nðn + 1Þ jn ðxÞ an ðxÞhðn1Þ ðxÞ =m2 , dx (4.2.65b) It is evident that 8 ξn ðxÞ ¼ xha ðxÞ; < ζ n ðxÞ ¼ xja ðxÞ, 0 ζp ð x Þ ¼ xj ð x Þ, ξ0n ðxÞﬃ ¼ xhb ðxÞ; pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ : n ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ b nðn + 1Þhðn1Þ ðxÞ ¼ hc ðxÞ: nðn + 1Þjn ðxÞ ¼ jc ðxÞ,

(4.2.66)

22 For a spherical particle, T11 n ¼ bn, and Tn ¼ an. Accordingly, Eq. (4.2.65) is the same as Eq. (4.2.64). This verifies the differential form of the IITM for a spherical particle using the Lorenz-Mie coefficients.

4.3

Application of the IITM to arbitrary particle morphologies

4.3.1 Stability and memory restriction The key equation for computing the T-matrix is Eq. (4.2.14) with respect to function F(r), which is a Fredholm equation of the second kind. The integral can be approximated by the choice of a quadrature rule in Eq. (4.2.23). The algebraic solution associated with n0 quadrature points is usually well conditioned because a diagonally dominant matrix inversion is encountered (Press et al., 2007, Chapter 19). The equation in Eq. (4.2.23) is itself a matrix equation. The matrix of a matrix equation, however, invokes super large matrices, which can become extremely time-consuming for large particles. Alternatively, the imbedding algorithm is employed to iteratively get the solution, as shown in Eq. (4.2.36). Moreover, the invariant-imbedding presentation of the T-matrix remains well conditioned. Accordingly, IITM can be applied to arbitrary particle morphologies.

180

Invariant Imbedding T-matrix Method

Even without the concern of an ill-conditioned matrix, the efficiency of the IITM depends on the memory and computational time requirements. The matrix computations involve matrix inversion, multiplication, and addition. A large matrix requires both large memory requirements and long computational times for the matrix operations. Apparently, we can have the following relations using the properties of the spherical Bessel function introduced in Section 3.2.1, as follows: n, n0 !∞ ðωn UGÞmnm0 n0 ! 0,

Qjj

mnm0 n0

n, n0 !∞ , Qjh mnm0 n0 , ðQhh Þmnm0 n0 , Qhj mnm0 n0 , Tmnm0 n0 ! 0:

(4.3.1) (4.3.2)

Accordingly, the matrices in the recurrence relation can be truncated at the term related to the current imbedding radius. Naturally, the truncation term increases with the imbedding process. We assume that the truncation term for the vector spherical wave functions is N, that is, ∞ X n X n ¼ 1 m ¼ n

!

N X n X

:

(4.3.3)

n ¼ 1 m ¼ n

As described in Section 3.2.2, index (m, n) can be replaced by index l on the basis of a unique correspondence. For clarity, we only consider the maximum truncation term N0 associated with the radius of the smallest circumscribed sphere of the scattering particle. The dimension of the vector spherical wave function is L ¼ N0(N0 + 2). For an arbitrary particle morphology, the dimensions of matrices U and G are 3L 3L, and the dimensions of matrices H and J are 3L 2L. The dimensions of matrices Qjj, Qjh, Qhh, Qhj, and T are 2L 2L. For a particle with N-fold rotational symmetry, however, the dimensions of the matrices are roughly decoupled as L/N. Note that the decoupled dimension depends on the remainder in Eq. (4.2.56) and is not exactly L/N. Accordingly, the dimensions for (int[N/2] + 1)-independent T-matrix components are roughly 2 NL 2 NL . For a particle with axial symmetry, the matrix is decoupled according to index m. Since only index values m > 0 are independent, the dimension for m ¼ 0 is N0, and the dimension for m ¼ 1, 2, …, N0 is (N0 m + 1). For instance, the dimensions for T-matrix component m ¼ 0 is 2N0 2N0, and the dimensions for T-matrix components with m > 0 are 2(N0 m + 1) 2(N0 m + 1), where the details are given in Table 4.1. Fig. 4.5 shows the memory requirements of the T-matrix associated with the size parameter. The truncation term N0 uses x +4.05x1/3 + 8 (where x ¼ krc) as an example to show the requirements. We will discuss the truncation term at the end of this section. The T-matrix is a complex matrix, and double precision is assumed in the computation. With an increased particle size, the memory requirements for the T-matrix significantly increase for particles with arbitrary morphologies. The memory requirements in the computation are much larger than shown in Fig. 4.5 because the transitional matrices, such as Qjj, Qjh, Qhh, Qhj, and U, have memory requirements similar to the T-matrix.

Invariant-imbedding T-matrix method

181

Table 4.1 Matrix dimension and total element number for an arbitrary morphology, finite-fold symmetric morphology, and axially symmetric morphology. L 5 N0(N0 + 2) Arbitrary morphology N-Fold symmetric morphology Axially symmetric morphology

Dimension Independent elements Dimension Independent elements Dimension for m ¼ 1, 2, …, N0 Dimension for m¼0 Independent elements

U

Qjj, Qjh, Qhh, Qhj, T

3L 3L 9L2

2L 2L 4L2

3 NL 3 NL

2 NL 2 NL

2

2

9 NL 2 ð int½N=2 + 1Þ

4 NL 2 ð int½N=2 + 1Þ

3(N0 m + 1) 3(N0 m + 1) 3N0 3N0

2(N0 m + 1) 2(N0 m + 1) 2N0 2N0

h i h i 9 N02 + N0 ðN0 + 16Þð2N0 + 1Þ 4 N02 + N0 ðN0 + 16Þð2N0 + 1Þ

Note that N indicates the morphology symmetry and N0 is the truncation term in the expansion using the vector spherical wave functions.

Memory requirement of T-matrix (G)

200

Arbitrary morphology 6-Fold rotationally symmetric morphology Axially symmetric morphology

150

100

50

0 0

50

100 Size parameter krc

150

200

Fig. 4.5 The memory requirement for the T-matrix associated with the size parameter in the unit of gigabytes. The truncation term N0 uses x + 4.05x1/3 + 8 (where x ¼ krc) as an example to show the memory requirement. The memory requirements at x ¼ 200 for arbitrary morphology, sixfold symmetric morphology, and axially symmetric morphology are roughly 175G, 19G, and 2.5G, respectively.

According to the analysis of the memory requirement for the T-matrix, two parallelization methods, Open Multi-Processing (OpenMP) and Message Passing Interface (MPI), are used to reduce the computational time and the memory restriction. OpenMP supports shared-memory programming. Even though it cannot reduce the memory restriction, the computational time can be significantly reduced. Also,

182

Invariant Imbedding T-matrix Method

OpenMP is easy to implement. On the other hand, MPI supports distributed-memory programming. Not only can it distribute the memory requirements to multiple nodes to significantly reduce memory restrictions, but also it can reduce the computational time. However, during the computation of the MPI, message passing in different processes is usually necessary, which can take time and increase the complexity of implementation. Based on the shape symmetry and the particle size, we use different parallelizations to improve the computational efficiency as shown in Table 4.2. For arbitrary morphologies, OpenMP is used for small particles (the estimated size parameter is x ¼ krc < 50). For large particles, MPI with message passing has to be used to reduce the memory requirements for each node. For finite-fold symmetric morphologies, the T-matrix is decomposed into several components shown in Eq. (4.2.56). We assign the independent components to different nodes and use OpenMP for the computation of the corresponding component. Accordingly, OpenMP and MPI without message passing are used for sizes x ¼ krc < 100. For large particles, MPI with message passing during the recurrence relation computation must be used. For axially symmetric morphologies, the T-matrix is decomposed into components according to index m. Each m-component can be assigned to a process, and no message passing is necessary during the recurrence computation. Since the required memory for axially symmetric morphologies is much smaller than for other morphologies, MPI without message passing can be usually used to calculate the T-matrix for intermediate particle sizes.

4.3.2 Extension to arbitrary morphology The stability of the invariant-imbedding method enables extensive applications not only for axially symmetric morphologies, finite-fold symmetric morphologies, and arbitrary morphologies but also for both homogeneous and inhomogeneous particles. As mentioned in Section 4.2.3, the imbedding process starts at radius ri of the largest inscribed sphere of the scattering particle and ends at radius rc of the smallest circumscribed sphere of the scattering particle as shown in Fig. 4.2A. Then, we can determine the quadrature points from radius ri to radius rc, so the intersected area between the scattering particle and the spherical shell at certain quadrature points can be determined. Consequently, the angular integrations in matrix U can be computed so that matrices Qjj, Qjh, Qhj, and Qhh can be obtained in terms of the spherical Bessel functions. The T-matrix can be updated step by step using the recurrence relation until the Table 4.2 Parallelization methods used for different morphologies at different sizes.

x < 50 50 < x < 100 100 < x

Arbitrary morphology

N-Fold symmetric morphology

Axially symmetric morphology

OpenMP MPI MPI

MPI + OpenMP MPI + OpenMP MPI

MPI MPI MPI

The numbers given in the table are rough or estimated.

Invariant-imbedding T-matrix method

183

quadrature point reaches radius rc. Technically, we need to find the starting radius ri, the ending radius rc, and the intersected area between the scattering particle and the spherical shell at a quadrature point to compute matrix U. The angular integrations in matrix U as given in Section 4.2.3 are computed using Gaussian quadrature for both the zenith (cosθ) and azimuthal (φ) integration. For simple shapes, such as a spheroid, a cylinder, or a hexagonal column, the starting and ending radii and the intersected area at certain quadrature points can be analytically given. For a complex shape, these parameters are not analytically obtained. We will discuss how to proceed with the invariant-imbedding algorithm for a complex shape with a given origin. Without loss of generality, we can assume that the complex shape is constructed by faceted meshes. Each facet can be specified by vertices and a normal direction. Moreover, the faceted meshes are assumed to not intersect. Similar discretization using faceted meshes has been used extensively to emulate a complex shape in many fields such as in technology and computer sciences.

4.3.2.1 Inscribed sphere and circumscribed sphere of a complex shape For a given origin at a complex shape, the radius rc of the smallest circumscribed sphere of the morphology can be easily obtained by comparing the distances between the origin and the vertices. In contrast, obtaining radius ri is not as apparent as obtaining radius rc. We can first find the minimum distance between the origin and one facet (a polygon) and then exhaust all the facets to find the least distance. The following gives the steps to obtain the minimum distance from the origin to a facet. The equation for a plane including a facet or a polygon can be formally expressed as ax + by + cz ¼ d:

(4.3.4)

We assume that point p0 is in the plane given in Eq. (4.3.4) and has the minimum distance from the origin. In Cartesian coordinates, the coordinates of point p0 can be expressed as p0 ¼

d ða, b, cÞ: a2 + b2 + c 2

(4.3.5)

Whether the point is on the boundary or not can be easily determined, and if the point is not on the boundary, whether the point is inside or outside the polygon can be determined using the ray-casting algorithm (Roth, 1982). If the point is on the boundary or inside the polygon, the minimum distance from the origin to the facet is straightforwardly computed by pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ jp0 j ¼ jd j= a2 + b2 + c2 :

(4.3.6)

If the point is outside the polygon, we can calculate the minimum distance from the point to the polygon. The polygon is composed of line segments so that the minimum

184

Invariant Imbedding T-matrix Method

distance from the point to the polygon can be obtained by comparing the minimum distances from the point to the line segments. The two vertices of a line segment are expressed by points p1 ¼ (x1, y1, z1) and p2 ¼ (x2, y2, z2). Any point on the line can be parameterized as p ¼ p1 + tðp2 p1 Þ,

(4.3.7)

where t is an unknown variable. The point on the line including points p1 and p2 that has the minimum distance from point p0 can be described as ð p p0 Þ ð p2 p1 Þ ¼ 0 tm ¼

ð p0 p1 Þ ð p2 p1 Þ : ð p2 p1 Þ ð p2 p1 Þ

(4.3.8a) (4.3.8b)

Accordingly, the point that has the minimum distance from point p0 is 8 tm 0; < p1 , pm ¼ p1 + tm ðp2 p1 Þ, 0 < tm < 1; : p2 , 1 tm :

(4.3.9)

Afterward, we can exhaust all the line segments to find the minimum distance from the point p0 to the polygon so that the minimum distance from the origin to the polygon can be easily obtained. Then, we can exhaust all the polygons to find the minimum distance from the origin to the polygons. The right minimum distance is the radius of the largest inscribed sphere.

4.3.2.2 Intersection between spherical shell at certain radius and scattering particle Gaussian quadrature is used for the angular integrations in matrix U. The zenith integration (μ ¼ cos θ) is replaced by Gaussian quadrature, so the spherical shell at a certain radius is decomposed into a series of circles according to the cosine of the zenith angle. For one of the circles with certain μ, the next step is to find the azimuthal intervals that can be obtained by locating the intersections between this circle and the facets. Then, the integrations for the zenith and azimuthal angles can be worked out so that matrix U can be given. The algorithm is also described by Bi et al. (2013). Without loss of generality, the circle at imbedding radius r0 and zenith angle θ0 is studied to find the azimuthal intervals. Any arbitrary point on the circle in Cartesian coordinate can be described as (r0 sin θ0 cos φ, r0 sin θ0 sin φ, r0 cos θ0). The equation of the plane containing a facet is given in Eq. (4.3.4). The intersection between the circle and the facet satisfies the following equation: C cos ðφ φ0 Þ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ , A2 + B2

(4.3.10a)

Invariant-imbedding T-matrix method

185

A ¼ ar 0 sin θ0 , B ¼ br 0 sin θ0 , C ¼ d cr 0 cos θ0 :

(4.3.10b)

B A sin φ0 ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ , cos φ0 ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ : 2 2 2 A +B A + B2

(4.3.10c)

The intersection points can be described as pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 8 2 + B2 ; no intersection, jCj > pAﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ > > ﬃ < 2 + B2 ; φ , j C j ¼ A 0

φ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ C > > : φ0 arccos pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ , jCj < A2 + B2 : 2 2 A +B

(4.3.11)

Since no intersection and only one intersection point do not contribute to the azimuthal integration, we shall deal with the one intersection case as no intersection. If there are two intersection points, the ray-casting algorithm (Roth, 1982) can be used to determine if the two intersection points are inside the facet or not. We can check all the facets to find all the intersection points. There could be no intersection or a finite number of intersection points. If there is no intersection between the circle and all facets, the circle could be inside or outside the scattering particle. An arbitrary point on the circle can be picked to determine if the circle is inside the particle or not. For convenience, the chosen point is p ¼ (r0 sin θ0, 0, r0 cos θ0) at 0-degree azimuthal angle, and the center of the circle is o ¼ (0, 0, r0 cos θ0). We can find all the intersection points between ray op and all the facets. If there is no intersection, point p and the circle must be outside the particle. If there are N intersections, the distances from point o to the intersection points can be sorted as d1 < d2 ⋯ < dN :

(4.3.12)

The outward normal to the facets corresponding to the order of Eq. (4.3.12) can be given as n^1 ! n^2 ⋯ ! n^N :

(4.3.13)

The sign function is defined as sign(x) ¼ x/jx j (the x ¼ 0 case does not occur here). Accordingly, we can have the following sign function between vector op and the normal to the facets as signðop n^1 Þ ! signðop n^2 Þ⋯ ! signðop n^N Þ:

(4.3.14)

The sign series in Eq. (4.3.14) alternates between 1 and 1. Then, we can conclude that the circle inside or outside the particle by 8

signðop n^1 Þ ¼ 1, outside, > > op < d , j j > 1 > signðop n^1 Þ ¼ 1, inside, < (4.3.15) dN < jopj, outside,

> > ^ ð Þ ¼ 1, outside, sign op n > i > : di < jopj < di + 1 , signðop n^Þ ¼ 1, inside:

186

Invariant Imbedding T-matrix Method

Whether the circle is inside or outside the particle is determined if there is no intersection between the circle and the facets. If there are M intersection points, the points can be sorted based on the azimuthal angles as follows: φ 1 ! φ2 ⋯ ! φ M ! φ1 :

(4.3.16)

All the azimuthal angles are arranged from 0 to 2π. The corresponding tangential vectors of the circle at the intersection points and outward normal to the facets stated in the order of Eq. (4.3.16) can be given as ! ! ! ν 1 ! ν 2⋯ ! ν M, ! ν i ¼ ð sin θ0 sin φi ,

sin θ0 cos φi , cos θ0 ð cos φi sin φi ÞÞ, i ¼ 1, 2,…, M, (4.3.17a)

n^1 ! n^2 ⋯ ! n^M :

(4.3.17b)

Accordingly, we can have the following sign functions corresponding to Eq. (4.3.16) as follows: ! ! ! ! sign ν 1 n^1 ! sign ν 2 n^2 ⋯ ! sign ν N n^N ! sign ν 1 n^1 :

(4.3.18)

The sign series in Eq. (4.3.18) alternates between 1 and 1. The signs equal to 1 and 1 mean the circle enters the particle and exits the particle, respectively. Accordingly, the corresponding azimuthal range [φi, φi+1] can be determined to be inside or outside the particle as 8 < Inside, if sign ! ν i n^i ¼ 1, ½ φi , φi + 1 : Outside, if sign ! ν i n^i ¼ 1:

(4.3.19)

Of course, only the range inside the particle is considered in the integration. The same steps are used to exhaust all the quadrature zenith points. Now, the radii of the largest inscribed sphere and the smallest circumscribed sphere of the particle and the intersected area between the spherical shell and the particle are solved in this section.

4.3.3 Discussion about the truncation term and imbedding step size The T-matrix method is a semianalytical method to obtain the light scattering properties of a complex shape. The term “analytical” here refers to the expansion using the vector spherical wave functions. Numerically, the computation of the T-matrix must truncate the VSWF from an infinite to a finite number of terms based on the particle size, as described in Eq. (4.3.3). The IITM is an imbedding method to obtain the

Invariant-imbedding T-matrix method

187

T-matrix of a scattering particle. The trapezoidal rule is used in the integration computation since the integrand is discontinuous. Accordingly, the step size can affect the computational accuracy. The truncation term N and the step size Δρ ¼ kΔr are discussed simply in this section. The truncation term N for spherical particles is given by Wiscombe (1980), and the expression usually used for spheres is N ¼ x + 4:05x1=3 + 2,

(4.3.20)

where x ¼ kr, r is the radius of the sphere, and Eq. (4.3.20) is obtained by fitting formula N x + cx1/3 suggested by Khare (1976). The physical interpretation is clear although Eq. (4.3.20) is an estimated formula. The first term x in Eq. (4.3.20) corresponds to a geometric term according to the localization principle described by van de Hulst (1957, p. 208), the second term x1/3 accounts for the edge effect (Khare, 1976; Wiscombe, 1980), and the remaining parameters are given by fitting a large amount of data. For nonspherical particles, the semianalytical T-matrix method is usually employed in light scattering. Fitting a formula for nonspherical particles is not realistic because not only the nonspherical particles have different morphology but also the memory and time requirements restrict the feasible amount of computation. When EBCM is used for a particle with axial symmetry, the truncation term is tested by determining the convergence of the T-matrix at component m ¼ 0 (Mishchenko, 1993). The convergence criterion is established for particles with axial symmetry. The IITM uses the imbedding method to recursively obtain the T-matrix of the particle. The same convergence criterion for the IITM can be used to determine the truncation term at every imbedding step when the particle has axial symmetry. Moreover, the convergence criterion for a particle with finite-fold rotational symmetry or complex shape still needs to be explored, since the IITM has excellent stability and can be applied to arbitrary morphologies. However, because of the enormous number of imbedding steps, the application of the convergence criterion on every imbedding step is not computationally efficient. Due to the outstanding stability of the IITM, the truncation term at an imbedding step is given as an input parameter described as h i 1=3 Ni ¼ int xi + 4:05xi + N0 ,

(4.3.21)

where variable N0 is an input parameter, xi ¼ krc, and rc is the radius of the smallest circumscribed sphere of the imbedded shape. Empirically, for small particles (krc < 100) and simple morphologies such as a spheroid, a cylinder, a hexagonal column, or a hexahedron, N0 ¼ 8 can be used to get convergent results, and N0 ¼ 20 is an adequate number for small particles. For large particles or complex morphologies, N0 should be increased to 30 or larger depending on the particle size and the complexity. The advantage for IITM is that a larger truncation term can give better convergence, but the trade-off is an increase of computational time. Another conservative formula for the truncation term can be used as h i Ni ¼ int xc + 4:05x1=3 + N0 , c

(4.3.22)

188

Invariant Imbedding T-matrix Method

where xc ¼ krc and the truncation term Ni at an imbedding step is the same as the one for the whole scattering particle. This formula gives a comfortable degree of convergence even for a particle with complex morphology (Liu et al., 2015). However, it takes a longer time than the one in Eq. (4.3.21). In Chapter 5, we will show that the formula in Eq. (4.3.22) can also be applied to large size particles and still give good convergence. The value of N0 is important because an insufficient truncation term can give an incorrect result. The recurrence relation for the T-matrix in Eq. (4.2.36) is derived from the quadrature expression of the integrations in Eqs. (4.2.22), (4.2.23). Using Eqs. (4.2.15), (4.2.16), the quadrature expressions of the integrations from radius rn1 to radius rn can be described as Tðrn Þ ¼ ik

ð rn1

dr 0 JT ðr 0 ÞFðrn j r 0 Þ + ikωn JT ðrn ÞFðrn j rn Þ,

(4.3.23)

0

Fðrn j r Þ ¼ Uðr ÞJðr Þ + Uðr Þ

ð rn1

dr 0 Gðr, r 0 ÞFðrn j r 0 Þ + Uðr ÞGðr, rn ÞFðrn j rn Þ:

0

(4.3.24) The recurrence relation between the T-matrix T(rn) at radius rn and the T-matrix T(rn1) at radius rn1 is actually obtained from the quadrature expression of the integration from radius rn1 to radius rn. Accordingly, the error originates from the quadrature expression. For the trapezoidal rule, the error can be usually expressed as follows (Atkinson, 1989): error ¼

ðkr c kr i Þ3 d2 F ðkr ξ Þ , 12n20 dðkr Þ2

(4.3.25)

where F (rξ) represents the integrand in Eqs. (4.3.23), (4.3.24) and rξ 2 [ri, rc]. Although Eq. (4.3.25) shows that more steps n0 can give more accurate results, the computational time is linearly increased with the number of steps. Moreover, the step number n0 should be increased or the step size Δρ should be decreased when the k(rc ri) is large. In Section 5.1, we will use the Lorenz-Mie coefficients as benchmarks to discuss the step size.

Invariant-imbedding T-matrix method

Invariant-imbedding T-matrix method

Recommend Documents