Reprint of: Numerical methods in electromagnetic scattering theory

Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

www.elsevier.com/locate/jqsrt

Reprint of

Numerical methods in electromagnetic scattering theory F. Michael Kahnert∗ Norwegian Institute for Air Research, P.O. Box 100, N-2027 Kjeller, Norway Received 14 May 2002; received in revised form 19 August 2002; accepted 19 August 2002

Abstract An overview is given over some of the most widely used numerical techniques for solving the electromagnetic scattering problem that start from rigorous electromagnetic theory. In particular, the theoretical foundations of the separation of variables method, the 5nite-di6erence time-domain method, the 5nite-element method, the method of lines, the point matching method, the method of moments, the discrete dipole approximation, and the null-5eld method (or extended boundary condition method) are reviewed, and the advantages and disadvantages of the di6erent methods are discussed. Aspects concerning the T matrix formulation and the surface Green’s function formulation of the electromagnetic scattering problem are addressed. ? 2003 Elsevier Science Ltd. All rights reserved.

Contents 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2. The electromagnetic scattering problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. Di6erential equation methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. The separation of variables method (SVM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1. Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2. Applications and practical issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3. Strengths and drawbacks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. The 5nite di6erence time domain method (FDTD) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1. Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2. Practical issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3. Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.4. Strengths and drawbacks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3. The 5nite element method (FEM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1. Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ∗

Tel.: +47-63-89-8035; fax: +47-63-89-8050. E-mail address: [email protected] (F.M. Kahnert).

0022-4073/03/$ - see front matter ? 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0022-4073(02)00321-7 This article is a reprint of a previously published article. For citation purposes, please use the original publication details; J Quant Spectrosc Radiat Transfer 2003;79-80:775-824.**DOI of original item: doi:10.1016/S0022-4073(02)00321-7 1791

777 777 780 780 780 783 783 784 784 785 786 787 787 787

776

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

Nomenclature CPU DDA DMF FDTD FEM GMT GOA GPMM MoL MoM NFM PMM SGF SVM

central processing unit discrete dipole approximation discretised Mie formalism 5nite-di6erence time-domain method 5nite-element method generalised multipole technique geometrical optics approximation generalised point-matching method method of lines method of moments null-5eld method point-matching method surface Green’s function separation of variables method

3.3.2. Practical issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3. Strengths and drawbacks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4. The method of lines (MoL) and the discretised Mie formalism (DMF) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5. The point-matching method (PMM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1. Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2. Underlying assumptions of the PMM and the GPMM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.3. Practical issues and applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.4. Strengths and drawbacks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. Volume-integral equation methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1. The volume-integral equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1. The singular kernel problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2. The method of moments (MoM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3. The discrete dipole approximation (DDA) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4. Practical issues and applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5. Strengths and drawbacks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. Surface-integral equation methods: The null-5eld method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1. The surface-integral equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2. The null-5eld method (NFM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3. NFM: Convergence issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4. NFM: Practical issues and applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5. NFM: Strengths and drawbacks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6. Other methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7. Intercomparisons of methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8. The T matrix formulation of the electromagnetic scattering problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.1. Ensembles of randomly oriented particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.2. Particles with point-group symmetries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9. Surface Green’s function (SGF) formulation of the electromagnetic scattering problem . . . . . . . . . . . . . . . . . . . . . . . 10. Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

788 788 789 793 793 794 796 796 796 796 797 798 799 801 802 802 802 803 805 806 807 808 808 810 811 812 814 816

Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 817 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 817

1792

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

777

1. Introduction Electromagnetic scattering by particles is a 5eld of active research with high relevance for such diverse 5elds as atmospheric science, oceanography, astronomy, and engineering sciences, with speci5c applications in remote sensing, meteorology, ocean optics, climate research, scattering by interplanetary dust grains, bio-optical imaging, antenna theory, particle sizing technology, and coating technology. The historical works by Mie [1] and Debye [2] considered electromagnetic scattering by spherical particles. Comprehensive descriptions of Mie theory can be found in the classical books by Stratton [3] and by Born and Wolf [4]. The history of works in the 5rst half of the 20th century on light scattering by particles is reviewed in the book by van de Hulst [5]. The purpose of this article is to provide a concise overview over some of the most widely used modern techniques for solving the electromagnetic scattering problem for nonspherical particles. This overview will be restricted to those techniques that are based on rigorous electromagnetic theory. Speci5c emphasis will be placed on reviewing the most important aspects concerning the theoretical foundations of di6erent methods, on pointing out possible interrelations, and on discussing the strengths and drawbacks of di6erent methods in view of practical applications. A comprehensive review of rigorous and approximate methods for performing light scattering computations for nonspherical particles and their applications is given in a recent book edited by Mishchenko, Hovenier, and Travis [6]. A recent paper by Wriedt [7] contains, besides a concise review of elastic light scattering by nonspherical particles, a discussion of extensions of Mie theory to ferromagnetic, coated, inhomogeneous, excentrically multilayered, and chiral spheres. This article is organised as follows. In Section 2, the theoretical basics of the electromagnetic scattering problem are reviewed and some notation is introduced. In Section 3 an overview is given over the separation of variables method (SVM), the 5nite-di6erence time-domain method (FDTD), the 5nite-element method (FEM), the method of lines (MoL), and the point-matching method (PMM) as typical exponents of di6erential equation methods. Section 4 reviews two volume-integral equation methods, namely, the method of moments (MoM) and the discrete dipole approximation (DDA). In Section 5 the null-5eld method (NFM) is discussed as an example for surface-integral equation methods. A few references to other light scattering methods are given in Section 6. Inter-comparisons of light scattering codes are discussed in Section 7. In Sections 8 and 9 some aspects pertaining to the T matrix formulation and to the surface Green’s function (SGF) formulation of the electromagnetic scattering problem, respectively, are discussed. Concluding remarks are given in Section 10. 2. The electromagnetic scattering problem Fig. 1 shows a schematic view of the electromagnetic scattering problem. An incident electromagnetic 5eld Einc , Hinc interacts with a particle that occupies a bounded region − in three dimensional space R3 . The total 5elds Etot and Htot in the surrounding medium + = R3 \− are equal to the sum of the incident and the scattered 5elds, i.e. Etot = Einc + Esca ; Htot = Hinc + Hsca . The total 5elds inside the particle are called the internal 5elds Eint ; Hint . The problem is to determine the unknown scattered and internal 5elds in terms of the known incident 5elds. The origin of the coordinate system is usually chosen to be an inner point of − .

1793

778

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

Fig. 1. Schematic view of the electromagnetic scattering problem.

To determine the unknown 5elds, one has to solve Maxwell’s equations 4 1 @D(r; t) ∇ × H(r; t) = ; J(r; t) + c0 c0 @t 1 @B(r; t) ; ∇ × E(r; t) = − c0 @t

(1) (2)

∇ · D(r; t) = 4(r; t);

(3)

∇ · B(r; t) = 0;

(4)

where cgs-units have been chosen, and where c0 denotes the speed of light in vacuo. In the electromagnetic scattering problem, we assume that there are no sources in space, i.e.  = 0 and J = 0 ∀r; t. Maxwell’s equations provide eight equations for the twelve unknown 5eld components. In addition to Maxwell’s equations one needs material equations to determine the 5elds. For low enough 5eld strengths the material equations can usually be assumed to be linear, i.e. D(r; t) = (r)E(r; t);

(5)

B(r; t) = (r)H(r; t);

(6)

where the electric permittivity and the magnetic permeability are for isotropic media scalar functions, otherwise tensors. With these assumptions Maxwell’s equations simplify to (r) @E(r; t) ∇ × H(r; t) = ; (7) c0 @t

(r) @H(r; t) ; (8) ∇ × E(r; t) = − c0 @t ∇ · E(r; t) = 0;

(9)

∇ · H(r; t) = 0:

(10)

1794

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

779

The magnetic permeability and the permittivity (r) = 0 r (r) (where r is called the relative permittivity, and where 0 denotes the permittivity in vacuo) are in general functions of the position vector r. (However, in many applications it is assumed that the surrounding medium is homogeneous and isotropic). Eqs. (7)–(10) constitute an over-determined system of eight equations for six unknown 5eld components. Thus one has to 5nd solenoidal (i.e. divergence-free) solutions to the Maxwell curl equations (7) and (8). Maxwell’s curl equations are hyperbolic di6erential equations that need to be solved subject to an initial condition. By taking the curl of Eq. (8) and the time-derivative of Eq. (7), one can substitute the latter equation into the former, which leads to the wave equation for the electric 5eld given by ∇ × ∇ × E(r; t) +

(r) (r) @2 E(r; t) = 0: @t 2 c02

(11)

Assuming a time-harmonic electric 5eld E(r; t) = exp(−i!t)E(r) one obtains the vector Helmholtz equation given by ∇ × ∇ × E(r) −

!2 (r) (r)E(r) = 0: c02

(12)

Using the vector identity ∇ × ∇ × A = ∇(∇ · A) − ∇2 A and Eq. (9), the vector Helmholtz equation becomes [∇2 + k 2 (r)]E(r) = 0;

(13)

k 2 (r) = !2 (r) (r)=c02 :

(14)

Analogously to Eq. (13) one derives the vector Helmholtz equation for the position-dependent part H(r) of the magnetic 5eld, [∇2 + k 2 (r)]H(r) = 0:

(15)

The vector Helmholtz equation is an elliptic di6erential equation that needs to be solved subject to boundary conditions. The tangential components of the total 5elds have to be continuous across the particle surface S (where S = @ is the boundary between the spatial domains − and + ). Further, the tangential components of the electric and magnetic 5elds have to approach zero as 1=r as the distance r from the origin approaches in5nity. This boundary condition at in5nity is called the radiation condition. Those numerical techniques that are classi5ed as di6erential equation techniques are based on numerically solving either Eqs. (7) and (8) subject to an initial condition, or Eqs. (13) and (15) subject to boundary conditions on the particle surface and at in5nity. In the following the separation of variables method (SVM), the 5nite-di6erence time-domain method (FDTD), the 5nite-element method (FEM), the point-matching method (PMM), as well as the method of lines (MoL) will be reviewed as typical exponents of di6erential equation methods.

1795

780

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

3. Di!erential equation methods 3.1. The separation of variables method (SVM) 3.1.1. Principle The basic idea of the SVM is to make a separation ansatz for the solution to the scalar Helmholtz equation [∇2 + k 2 ] (r) = 0;

(16)

and to obtain a set of di6erential equations for each component function from the scalar Helmholtz equation. From the set of solution functions to these di6erential equations one can construct solenoidal vector wave functions that solve the vector Helmholtz equation. The incident, scattered, and internal 5elds are expanded in suitable vector wave functions, and the expansion coeOcients are determined by enforcing the boundary conditions on the particle surface. In spherical coordinates, e.g., Eq. (16) becomes
 1 @2 @2 @ @ 1 1 sin  + r + + k 2 = 0: (17) r 2 sin  @ @ r @r 2 r 2 sin2  @2 With the separation ansatz (r; ; ) = R(r)()();

(18)

one obtains three ordinary di6erential equations 
d 2 (rR) n(n + 1) 2 0= rR; + k − dr 2 r2 

d m2 1 d ; sin  + n(n + 1) − 0= sin  d d sin2  0=

d2  + m2 : d2

(19) (20) (21)

where m2 and n(n+1) are separation constants. The solutions R(r)=zn( j) (r) to the radial equation are spherical Bessel functions (zn(1) = jn ) spherical Neumann functions (zn(2) = nn ), and spherical Hankel (4) (2) functions of the 5rst kind (zn(3) = h(1) n = jn + inn ) and of the second kind (zn = hn = jn − inn ), where n ∈ N. The solutions () = Pnm () to the polar equation are associated Legendre functions, where m = −n; −n + 1; : : : ; n. Finally, the solutions () = exp(im) to the azimuthal equation are harmonic functions. From the set of solution functions ( j) n; m (r; ; )

= zn( j) (r)Pnm () exp(im)

(22)

to the scalar Helmholtz equation in spherical coordinates one can construct a set of vector wave functions a

Ln;( j)m = ∇ · (a

( j) n; m );

(23)

1796

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824 a

Mn;( j)m = ∇ × (a

( j) n; m );

a

Nn;( j)m = 1=k∇ × a Mn;( j)m :

781

(24) (25)

If we choose the vector a used in the construction of the spherical vector wave functions to be either the position vector a = r or any constant vector (e.g. a Cartesian unit vector a = x; ˆ a = y, ˆ or a = z), ˆ then one can show that the spherical vector wave functions satisfy the vector Helmholtz equation in spherical coordinates. However, only the functions a Mn;( j)m and a Nn;( j)m are solenoidal. Due to Eqs. (9) and (10), the functions a Ln;( j)m are not suitable for representing the solution to the electromagnetic scattering problem, because they are not divergence-free. The incident, scattered and internal electric 5elds can now be approximated by an expansion in spherical vector wave functions according to Einc (k0 r) =

Ncut  n  n=0 m=−n

sca

E (k0 r) =

Ncut  n  n=0 m=−n

int

E (kr) =

Ncut  n  n=0 m=−n

(Ncut ) (1) (1) cut ) [an;(Nm; M Mn; m (k0 r) + an; m; N Nn; m (k0 r)];

(26)

(Ncut ) (3) (3) cut ) [pn;(Nm; M Mn; m (k0 r) + pn; m; N Nn; m (k0 r)];

(27)

[cn;(Nm;cutM) Mn;(1)m (kr) + cn;(Nm;cutN) Nn;(1)m (kr)]:

(28)

Here the wavenumber inside the particle has been denoted by k and the wavenumber in the surrounding medium has been denoted by k0 . For the expansion of the incident and internal 5elds one employs spherical vector wave functions of the 5rst kind (which are regular at the origin), whereas the scattered 5eld is expanded in spherical vector wave functions of the third kind (which satisfy the radiation condition). In an analogous way, one can expand the incident, scattered and internal magnetic 5elds. These expansions are then substituted into the boundary conditions on the particle surface nˆ+ × [Einc (k0 r) + Esca (k0 r) − Eint (kr)]|r∈@ = 0;

(29)

nˆ+ × [Hinc (k0 r) + Hsca (k0 r) − Hint (kr)]|r∈@ = 0;

(30)

which express the continuity of the tangential 5eld components across the boundary surface @. Here nˆ+ denotes the outward pointing normal vector to the boundary surface. The boundary conditions cut ) yield a system of linear equations for determining the unknown expansion coeOcients pn;(Nm; $ and (Ncut ) cn; m; $ ($ = M; N ) of the scattered and internal 5elds, respectively, in terms of the known expansion cut ) coeOcients an;(Nm; $ of the incident 5eld. For spherical particles this system of linear equations can be inverted analytically. This inversion leads to the well-known solution [1,2] known as Mie theory. In this special case the expansion coeOcients are independent of the truncation parameter Ncut . Another way of deriving the Mie solution is the classical approach based on using Debye potentials [4]. The procedure for deriving a set of vector wave functions can be applied to any other coordinate system in which the scalar Helmholtz equation becomes separable. In spheroidal coordinates, e.g., the spherical vector wave functions will be replaced by spheroidal vector wave functions [8,9]. Subsequently one can substitute an expansion of the 5elds in the set of spheroidal vector wave functions

1797

782

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

into the boundary condition on the particle surface for determining the unknown coeOcients. This was done for spheroidal particles by Oguchi [10], Asano and Yamamoto [11] and by Sinha and McPhie [12]. However, for a spheroidal particle the resulting system of linear equations can only be inverted numerically. The origin for this is that the spheroidal vector wave functions are not orthogonal, even though they are complete and even though the scalar spheroidal wave functions are orthogonal. The original approach taken, e.g., by Asano and Yamamoto [11] was based on employing spheroidal vector wave functions r Mn;( j)m and r Nn;( j)m that are constructed by replacing in Eqs. (24) and (25) the scalar spherical wave functions by corresponding scalar spheroidal wave functions, and by choosing a = r. Farafonov [13] and Voshchinnikov and Farafonov [14] persued a di6erent approach by splitting up the 5elds into an axisymmetric component, which can be treated by use of Abraham’s potentials, and a non-axisymmetric component, which is expanded in a mixed set of spheroidal vector wave functions r Mn;( j)m ; r Nn;( j)m ; zˆMn;( j)m , and zˆNn;( j)m . Asano and Yamamoto’s approach [11] is equivalent to describing the electromagnetic 5eld by means of Debye potentials, which is similar to the Mie solution for spheres. The approach by Farafonov [13] and Voshchinnikov and Farafonov [14] is equivalent to representing the electromagnetic 5eld by a combination of Debye potentials and Hertz potentials, which are employed for solving, respectively, the scattering problem for spheres and in5nite circular cylinders. This ansatz results in a larger Rexibility for treating both mildly aspherical spheroids as well as highly elongated prolate or strongly Rattened oblate spheroids. Asano and Yamamoto [11] implemented expansions for the radial spheroidal wave functions of the second kind that are only valid for mildly aspherical spheroids [8]. This severely limited the range of aspect ratios accessible to their implementation. Farafonov [15] was able to obtain numerical results for the radiation scattered by prolate spheroids with an aspect ratio a=b as small as 1/100 and by an oblate spheroid with an aspect ratio as large as 100. Here b and a denote the size of the spheroid along its symmetry axis and perpendicular to its symmetry axis, respectively. Voshchinnikov and Farafonov [16] propose a more Rexible numerical treatment of the radial spheroidal functions based on using di6erent series expansions. This approach has the potential of treating large, elongated particles. Voshchinnikov and Farafonov [14] estimate that their SVM implementation is faster than the SVM formulation by Asano and Yamamoto [11] by 1-2 orders of magnitude, depending on the spheroid’s aspect ratio. In applications one is usually interested in the scattered electric 5eld far away from the scattering particle. Far-5eld scattering is described by the (2 × 2) amplitude scattering matrix, which relates the two transverse components of the incident electric 5eld to the two transverse components of the scattered electric 5eld in the far-5eld region. Once the expansion coeOcients of the scattered 5eld have been determined, it is straightforward to compute the amplitude scattering matrix and the optical cross sections [11]. From the amplitude scattering matrix one can compute the phase matrix as described in Ref. [6]. In the incoherent-scattering case, the phase matrix of an ensemble of independently scattering particles can be obtained by adding up the phase matrices of each individual particle in the ensemble. The ensemble-averaged phase matrix can then be used as input to radiative transfer computations [17,18]. If one is interested in computing the scattered 5eld for di6erent orientations of the particle with respect to the incident 5eld, then a remaining drawback of the SVM in its conventional formulation is that one needs to repeat the entire numerical procedure for each orientational angle. In particular, computing the phase matrix and other optical properties of ensembles of randomly oriented spheroids

1798

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

783

requires us to compute the scattered 5eld for a large number of discrete orientational angles followed by numerical integration over these angles [19]. Such computations can place high demands on CPU-time resources. Historically this general problem in electromagnetic scattering computations was 5rst overcome in the T matrix formulation (see Section 5). Mishchenko [20] developed a method for analytically calculating the orientationally averaged optical properties for ensembles of randomly oriented particles from the electromagnetic scattering solution in the T matrix formulation. In Ref. [21] these results were exploited in the SVM by reformulating the method in such a way as to compute a T matrix in spheroidal coordinates. The T matrix was obtained by taking as incident 5elds the spheroidal vector wave functions instead of the commonly used plane waves. 3.1.2. Applications and practical issues Extensions of the SVM exist for coated [22–26] and chiral spheroids [27] as well as to systems of two and more spheroids [28–31]. Spheroidal particles are frequently used as proxies for nonspherical particles in sensitivity studies with regard to questions relevant to atmospheric optics [32–35] and astrophysical applications [14,36,37], to name but a few examples. A review of the application of the SVM to spheroidal particles has been presented by Ciric and Cooray [38]. As mentioned earlier, an important practical issue in the SVM concerns the numerical computation of the spheroidal wave functions. Usually one expands the radial and angular spheroidal functions in a suitable set of other functions that can easily be computed numerically [8]. Our ability to accurately compute the spheroidal wave function critically depends on the development of suitable series expansions [16] as well as on reliable routines for computing the expansion coeOcients [39]. Currently, accurate numerical computations of spheroidal wave functions have been reported [39] for an inter-focal distance f of the spheroids up to f = 40', where ' denotes the wavelength of light. 3.1.3. Strengths and drawbacks An important advantage of the SVM is its high numerical accuracy, making it a method suitable for benchmark computations. Another advantage is that the SVM can be used to compute a T matrix [21]. This allows us to calculate the optical properties of an ensemble of randomly oriented particles analytically and thus very eOciently from the T matrix. The SVM solution automatically satis5es the radiation condition. The spheroidal vector wave functions suit the geometry of spheroidal particles. As a result, Voshchinnikov [26] could demonstrate that the expansion of the 5eld scattered by a spheroid converges at a speed that only depends on the size parameter kas of the spheroid (where as denotes the spheroid’s semi-major axis), but not on the particle’s aspect ratio. A disadvantage of the SVM is that for large size parameters of the spheroid and/or large refractive indices |m| the system of linear equations to be solved becomes large, and ill-conditioning problems may occur. The computational complexity (i.e. the number of numerical operations in the algorithm) of the SVM scales between N ∼ O(x3 ) and N ∼ O(x4 ), where x = k0 r denotes the particle’s size parameter, i.e. a typical size r multiplied with the wavenumber k0 . As a 5nal remark, it is worth noting that traditionally the term SVM has been restricted to solving the electromagnetic scattering problem for particles the surface of which coincides with

1799

784

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

a coordinate-hypersurface in one of the eleven coordinate systems in which the scalar Helmholtz equation is separable [40]. However, the essential assumption in the SVM is actually that the incident, scattered, and internal 5elds can be approximated by a 5nite expansion on the boundary @, and that the coeOcients in the expansion of the surface 5elds are identical with the corresponding coeOcients of the expansion of the incident and scattered 5elds in + and the internal 5elds in − , respectively. This assumption is not only satis5ed for separable boundary surfaces [41,42]. Thus the SVM can be generalised to non-separable boundary surfaces. We will refer to this as the generalised SVM. 3.2. The
(31)

(32)

(33)

and similar equations are obtained for the y and z components of the induction law, and for the x; y, and z components of Faraday’s law [Eq. (7)]. c0 denotes the propagation velocity of the electromagnetic 5eld. Eq. (33) illustrates that the solution algorithm in the FDTD is iterative. The x component of the magnetic 5eld is computed at a time step n + 1=2 in terms of the magnetic and electric 5elds at earlier time steps. Thus the essence of the FDTD algorithm is to numerically solve an initial-value

1800

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

785

problem by marching a plane wave or pulse source, which is “switched on” at some initial time, through discrete time steps over a 5nite discretised spatial domain that includes the particle. The scattering particle is represented in the FDTD by its permittivity , which in general can vary in space. To each grid cell an average value of , which occurs in Faraday’s law [Eq. (7)], is assigned. Yang and Liou [45] recommend to employ the Maxwell–Garnett [46] rule  y+Sy  z+Sz  x+Sx (i; V j; k) − 1 (x ; y ; z  ) − 1 1 = (34) d x dy d z (i; V j; k) + 2 SxSySz x (x ; y ; z  ) + 2 y z for evaluating the mean permittivity in each cell. The incident wave given by the initial condition is usually chosen to be either a plane wave or a pulse. If a sinusoidal plane wave is chosen, then the iteration process is continued until the transient solution for the 5elds has converged to a steady-state solution. If a pulse source is chosen, then the time-iteration is continued until the 5eld values have decayed below a given threshold value. The near 5eld solution one obtains for the 5elds inside the 5nite computational domain of the spatial grid mesh is Fourier-transformed into the frequency domain and subsequently propagated to the far 5eld. The far-5eld propagation is performed by employing either a surface-integration technique [47–49] or, for non-conducting particles, a volume-integration technique [50]. From the frequency-domain far 5eld solution one can compute the optical properties of the particle, such as the amplitude scattering matrix, the extinction cross section, the scattering cross section, and the absorbption cross section [45]. From the amplitude scattering matrix, one can obtain the phase matrix [6]. 3.2.2. Practical issues It should be emphasised again that the spatial grid mesh, in which the time-evolution of the electromagnetic 5eld is iteratively computed from an initial state, comprises a
(36)

with the pseudo-di6erential one-wave operators  1=2 2 2 2 2 2 (@ =@y + @ =@z ) 1 @ c @ ± ± 1− : Lˆx = @x c @t @2 =@t 2

(37)

Let’s suppose that part of the boundary @D of the computational domain is located at x = 0, and that the scattered 5eld is propagating in the positive x-direction (i.e. the particle is located somewhere

1801

786

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

in the spatial region x ¡ 0). The condition for the scattered 5eld component Ezsca that no reRections should occur in the negative x-direction (i.e. back into the computational domain) at x=0 is expressed by the one-way wave equation sca Lˆ− x Ez |x=0 = 0:

(38)

The pseudo-di6erential one-way wave operator Lˆ− x can be approximated, e.g. by performing a Taylor expansion up to second order 1 @ 1 c(@2 =@y2 + @2 =@z 2 ) @ − + : Lˆ− x ≈ @x c @t 2 @=@t

(39)

With this approximation, the one-way wave equation can be discretised and used in the FDTD algorithm as an approximate condition for minimising reRections at the arti5cial boundary @D back into the computational domain. A more recent method is based on constructing a multi-layered medium at the arti5cial boundary @D consisting of uniaxial layers that minimises reRections back into the computational domain D [53,54]. This so-called uniaxial perfectly matched layer has been demonstrated to cause reRections that are smaller by three orders of magnitude as compared to the approximate analytical methods [53,55,56]. The good performance of the uniaxial perfectly matched layer results in a higher computational stability and allows one to use a smaller free space between the particle and the boundary, which places lower demands on computer memory and central processing unit (CPU) time. Another practical issue in the FDTD that has received considerable attention is the method of spatial discretisation. In a rectangular grid mesh a curved particle surface is represented by a staircase of rectangular cells. To avoid possible computational artifacts associated with this, various non-rectangular grids have been devised [57–62]. However, rectangular grids are signi5cantly easier to implement and not limited to certain particle geometries. In many cases, particularly for non-ferromagnetic and non-conducting particles, numerical inaccuracies associated with a rectangular grid can be avoided [49,50] by using the Maxwell–Garnett rule [46] given in Eq. (34) to compute the mean values of the dielectric constant at the grid points. Finally, the size of the discrete volume cells and of the discrete time steps are critical for achieving a suOciently high numerical accuracy. The spatial increments Sx; Sy, and Sz should not be larger than '=20, where ' denotes the wavelength of the incident wave. The time increment should satisfy the Courant–Friedrichs–Levy condition [63] given by c0 St 6 [1=Sx2 + 1=Sy2 + 1=Sz 2 ]−1=2 :

(40)

3.2.3. Applications The FDTD is a well-tested method that has been applied to such diverse applications as antenna scattering [64], analysis of microstrip circuits [65], absorption in tissue [66], or scattering by atmospheric ice crystals [50], irregular aggregates [67], multi-faceted concave and convex particles [68], spheres with nonspherical inclusions [45], hexagonal columns with air bubble and soot inclusions [45], or clusters of particles with or without inclusions [45], to name but a few examples. For 2D scattering by in5nite circular cylinders in a 5xed orientation, Yang and Liou [49] have reported results for size parameters up to 60. 3D scattering results for spheres of size parameters up to 40 have been reported by Sun et al. [69]. Sun and Fu [70] have shown that for spheres

1802

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

787

of high refractive indices up to values of 7:1499 + 2:914i and size parameters up to 6 their modern FDTD implementation yields results for scalar optical properties (extinction and absorption eOciency) with an error less than 4%, and for the scattering phase function with an error less than 5%. 3.2.4. Strengths and drawbacks The FDTD algorithm is conceptually simple and straightforward to implement. One of its main advantages is the high Rexibility in terms of particle geometries. The FDTD is applicable to inhomogeneous, anisotropic, and arbitrarily shaped particles. Due to the time-iterative solution process, one avoids the problem of solving large systems of linear equations, as encountered in many other numerical techniques. A drawback is that one needs to perform computations over a spatial domain that is larger than the particle (in contrast to, e.g., volume integral equation methods). This can result in time consuming computations, especially for larger particles. Another disadvantage is that computations need to be repeated for each orientation of the particle with respect to the incident wave, which makes this method rather ineOcient for applications involving ensembles of randomly oriented particles. The number of numerical operations N in the FDTD algorithm increases approximately with the fourth power of the particle size parameter x; N ∼ O(x4 ). 3.3. The
1803

788

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

3.3.2. Practical issues Both in the FEM and in the FDTD the scattering problem is solved in a 5nite spatial region only. In the FDTD this leads to the challenge of 5nding suitable methods for ensuring that the propagating transient waves are not reRected at the arti5cial boundary @D back into the computational domain D. In the frequency-domain FEM, by contrast, one is not dealing with a transient propagating wave, but with the steady-state spatial distribution E(r) and H(r) of the electric and magnetic 5elds. Whereas Maxwell’s equations in the time-domain are hyperbolic di6erential equations that are solved as an initial value problem, Helmholtz’ equation is an elliptic di6erential equation that is solved as a boundary value problem. The solution has to satisfy not only the boundary conditions on @, but also the boundary condition at in5nity. Thus one challenge in the FEM is to 5nd suitable ways of ensuring that the FEM-solution computed in the 5nite computational domain satis5es the radiation condition. One way to emulate the radiation condition is to derive approximate boundary conditions at the edge of the computational domain [75]. Another way to enforce the radiation condition is achieved by combining the FEM with a surface-integral method [76,77] which can, however, destroy the banded structure of the coeOcient matrix. Probably the most elegant way to enforce the radiation condition in the FEM is the so-called unimoment method [71]. In this approach the computational domain D is taken to be a sphere of radius R containing the particle. In the medium outside D (r ¿ R), the scattered 5eld is expanded in outgoing spherical vector wave functions, just like in Eq. (27). The spherical vector wave functions of the third kind automatically satisfy the radiation condition. The FEM-solution inside D and the expansion coeOcients of the 5elds outside D are matched on the surface @D of the computational domain. Another important practical issue is the choice of the grid mesh. Whereas rectangular meshes are quite popular in the FDTD because of their Rexibility and generality, it is more common in the FEM to use a grid mesh that conforms with the particle surface, such as the semiannual conformal mesh [71]. This allows to naturally enforce the boundary conditions on the particle surface. Moreover, the spherical shape of the semi-annual conformal mesh allows easy use of the unimoment method for enforcing the radiation condition. The size of the cells should not be larger than '=20. 3.3.3. Strengths and drawbacks The FEM is conceptually simple and straightforward to implement. It is in principle applicable to arbitrarily-shaped and/or inhomogeneous particles, although particle-conforming grid meshes need to be adapted to the individual particle shapes. The fact that the coeOcient matrix is banded is an advantage over volume-integral equation methods. The singular-kernel problem of volume-integral equation methods and the related discussions about suitable approximations of the polarisability as a function of permittivity is avoided in the FEM (see Section 4). However, FEM computations need to be carried out over a computational domain that is larger than the particle. This is a disadvantage compared to integral equation methods. Another major drawback is the necessity to repeat computations for each new orientational angle of the particle with respect to the incident 5eld, which results in high CPU-time requirements when applying the FEM to ensembles of randomly oriented particles. The accuracy of computations strongly depends on the spatial discretisation and is therefore diOcult to control.

1804

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

789

The computational complexity of the FEM is N ∼O(x7 ) in connection with Gaussian elimination, and N ∼O(x4 ) in connection with the conjugate gradient method. 3.4. The method of lines (MoL) and the discretised Mie formalism (DMF) The MoL [78,79] is a method frequently applied for solving partial di6erential equations in guided-wave theory. Its application to scattering problems is referred to as the DMF [80,81]. The DMF is not quite as simple as the formalisms of most other di6erential equation techniques, but it possesses a few methodical advantages, it reveals some interesting aspects on the relation between methods based on approximating the solution to the scattering problem and methods based on approximating the di6erential operator in the vector Helmholtz equation, and it can be used to derive an interesting new formulation of the solution to the electromagnetic scattering problem. The DMF is based on discretising not all three spatial dimensions, as in the FEM, but (in spherical coordinates) only the polar and (for non-axisymmetric particles) the azimuthal coordinate. This leads to a di6erential equation for the radial coordinate that can be brought into Bessel’s form. Among all solution functions one picks out those that represent outgoing waves. Thus, in contrast to the FEM, the DMF analytically satis5es the radiation condition. After expanding the general solution in this basis and after determining the unknown expansion coeOcients by enforcing the boundary conditions, one obtains a surface Green’s function solution for the electromagnetic scattering problem that is valid not only for separable, but for all star-shaped geometries. To illustrate some general aspects of this formalism it is enough to restrict oneself to the problem of scattering by axisymmetric particles, and to consider the scalar problem of determining the solution 0 in + of the scalar Helmholtz equation subject to a Dirichlet boundary condition on the boundary surface @, and subject to the radiation condition. For the axisymmetric problem, the potential 0 can be expanded in a Fourier series  0(r; ; ) = 0(l) (r; ) exp(il): (41) l

In + each Fourier component 0(l) satis5es the scalar Helmholtz equation [∇2 + k02 ]0(l) = 0

(42)

where ∇2 = r

@ l2 @2 1 @ sin  − r+ : 2 @r sin  @ @ sin2 

(43)

The main conceptual di6erence to other 5nite-di6erence techniques, such as the FDTD or the FEM, is that only some but not all spatial coordinates are discretised in the di6erential equation. For an axisymmetric problem formulated in spherical coordinates, only the polar variable is discretised, as shown in Fig. 2, but not the radial variable. Thus @0 1 ≈ [0(r; i+1 ) − 0(r; i )] (44) @ i h

1805

790

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

Fig. 2. Discretisation of the polar coordinate in the MoL.

and similarly for second derivatives. With this scheme Eq. (42) is transformed into  

 2 d 2 2 2 (l) |0(l) = |0 ; h r 2 r + k0 r 1 − P dr

(45)

where the vector |0(l) contains the components 0(l) (r; 0 ); 0(l) (r; 1 ); : : : ; 0(l) (r; Nd +1 ), where Nd denotes the number of division points used in the discretisation of the polar interval [0; ]. The matrix P(l) contains the 5nite-di6erence quotients from the discretisation of the polar derivatives in Helmholtz’ equation [82,83]. Since only derivatives up to second order occur, this matrix is band triagonal. It can be symmetrised with a suitable transformation matrix Z according to (l) P(l) → Psym = Z−1 · P(l) · Z

(46)

and diagonalised (l) − ((l) ) · X(l) = 0; (Psym

(47)

(l) , and where the column vectors of where the diagonal matrix ((l) contains the eigenvalues of Psym (l) (l) the matrix X are the eigenvectors of Psym . Eq. (45) now becomes  

 2 1 (l) d 2 2 |U (l) = |0 ; (48) r 2 r + k0 r 1 − 2 ( dr h

1806

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

|U (l) = (X(l) )−1 · Z−1 · |0(l) ;

791

(49)

which can be brought into Bessel’s form. The radiating solutions are Ui(l)

H (1) (l) (k0 r) vi ; = √ k0 r

vi(l)2 =

'i(l) 1 + : 4 h2

(50)

Thus the 5rst remarkable result is that the radiation condition can automatically be satis5ed [82,83]. This is a distinct advantage of the DMF over other 5nite-di6erence techniques, which is a direct consequence of the fact that the radial coordinate has not been discretised in this method. At this point one has obtained a set of general radial solution functions that satisfy the boundary condition at in5nity. To obtain a concrete solution |0(l) (r) in + belonging to a particular Dirichlet potential |0S(l) given on the boundary surface S = @, one expands |0(l) in a suitable set of expansion functions and determines the unknown expansion coeOcients by enforcing the boundary condition. As a set of expansion functions, one can choose [41] H5(l) (r) · |x5(l)

(51)

   H (1) (k r) 0  (l) v5 √ ; H5(l) (r) = diag  k0 r 

(52)

where |x5(l) denote the eigenvectors in the matrix X(l) in Eq. (47). Thus one expands the desired solution |0(l) (r) in the surrounding medium according to (l)

|0 (r) =

Ncut 

a5(l; Ncut ) H5(l) (r) · |x5(l) ;

(53)

5=0

where Ncut 6 Nd , i.e. the truncation index of the sum approximation (and thus the number of unknown expansion coeOcients) can be chosen smaller than the number of division points (which is equal to the dimension of the space spanned by the functions H5(l) (r) · |x5(l) ). The expansion given in Eq. (53) is substituted into the Dirichlet boundary condition |0(l) (R ) = |0S(l) (R ) ;

(54)

where R denotes the parameterisation of the boundary surface in spherical coordinates. This results in a system of linear equations for determining the unknown expansion coeOcients a5(l; Ncut ) . The formal inversion of this equation system yields a5(l; Ncut ) =

Ncut 

(l) −1 (l) {A(l) @ }5; 6 x6 |0S (R ) ;

5 = 0; : : : ; Ncut ;

(55)

6=0

where the geometry factors (l) (l) (l) {A(l) @ }5; 6 = x6 |H5 (R )|x5

(56)

contain the information about the geometry of the boundary surface @.

1807

792

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

We now recall that the superscript l denotes expansion orders in a Fourier expansion with respect to the azimuthal coordinate . Thus  exp(il)|0(l) (r) (57) |0(r; ) = |0S(l) (R )

l =

2

0

d exp(−il )|0S (R ;  ) :

(58)

Eqs. (53) and (55) in conjunction with Eqs. (57) and (58) yield the solution |0(r; ) in + in the form  2 (Ncut ) |0(r; ) = d G@ (r; ; R ;  )|0S (R ;  ) (59) 0

with the surface Green’s function de5ned by (Ncut ) (r; ; R ;  ) G@

=

Ncut  l

(l) −1 (l) (l) exp [il( −  )]{A(l) @ }5; 6 H5 (r)|x5 x6 |:

(60)

5;6=0

Thus the second remarkable result is that the DMF applied to Helmholtz’ equation results in a surface Green’s function solution of the electromagnetic scattering problem [41]. As shown in Ref. [42] one can also use the DMF approach to derive a surface Green’s function in such a way that the sum approximation of the solution converges on @ in the least squares sense, and uniformly, in + . The surface Green’s function formulation will be discussed in more detail in Section 9. A third interesting point is the following. The derivation of the formalism started with an approximation of the di6erential operator of the Helmholtz equation [Eq. (45)]. Therefore it has earlier been assumed that the relative convergence phenomenon occurring in methods based on approximating the solution function with a sum expansion does not occur in the MoL [78]. However, as 5rst noted by Rother [41], this is not true. As can be seen in Eq. (53), the original approximation of the di6erential operator can be reformulated as a sum approximation of the solution. The unknown expansion coeOcients are determined by enforcing the boundary conditions on the boundary surface. When Ncut = Nd , the conventional MoL is actually equivalent to the point matching method (see Section 3.5). A fourth interesting point to note is that numerical experiments have revealed [41,83] that in the limit of an increasing number of lines Nd the radial and angular functions converge to the spherical Hankel functions of the 5rst kind and to the associated Legendre functions, respectively N →∞

|x5(l) d→ |P5(l) ;

(61)

H (1) (l) (k0 r) Nd →∞ (1) v5 √ → h5 (k0 r): (62) k0 r With this result, Eq. (60) can be rewritten as an expansion of the surface Green’s function in the basis of radiating spherical wave functions. In spherical coordinates this numerically observed phenomenon could not yet be shown mathematically. However, in Cartesian coordinates one can show analytically that the expansion functions in the DMF converge to the corresponding eigenfunctions of Helmholtz’ equation in Cartesian coordinates [83]. This aspect, which in earlier discussions about the nature of the MoL [84,85] has

1808

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

793

been overlooked, was 5rst noted by Rother [83,41]. Thus in the limit of an increasing number of discretisation lines the MoL leads to an expansion of the solution in spherical functions for arbitrary star-shaped surfaces, where the expansion coeOcients of the solution on the boundary surface and in the surrounding medium are identical. However, as noted in Section 3.1, this is just the essence of the generalisation of the SVM. Consequently, when applied in a coordinate systems in which the scalar Helmholtz equation is separable the MoL is just a complicated discretised form of the generalised SVM. Thus, in practical applications the MoL can only o6er advantages over the SVM if the boundary value problem is formulated in more complex coordinate systems in which the Helmholtz equation is not separable, or if applied to other di6erential equations (e.g. Sturm-Liouville). The DMF has been applied to a few selected nonspherical geometries, such as spheroids [80,81] and “potato”-shaped hydrometeors [83]. The main methodical advantages of the DMF are that the radiation condition is automatically satis5ed, and that the formalism yields a surface Green’s function formulation of the scattering problem. As will be discussed in Section 9, this can be done in such a way that the resulting approximation of the solution has mathematically very well-de5ned convergence properties. However, nothing can be gained from applying the DMF rather than the generalised SVM in all cases in which the generalised SVM is applicable. 3.5. The point-matching method (PMM) 3.5.1. Principle In the PMM [10] the electric and magnetic 5elds are represented as truncated expansions in spherical vector wave functions as given in Eqs. (26)–(28), and analogous expansions for the magnetic 5elds. Subsequently a 5nite set of points on the boundary surface @ is chosen (where, in contrast to the conventional applications of the SVM, the boundary surface is not necessarily assumed to be separable). The boundary conditions given in Eqs. (29) and (30) in conjunction with the corresponding 5eld expansions are enforced in the 5nite number of matching points. This yields a system of algebraic equations K 

Ai; j dj = bi ;

i = 1; : : : ; L;

(63)

j=1 cut ) where K = 4(Ncut + 1)2 is the number of unknown coeOcients cn;(Nm;cut$) ; pn;(Nm; $ (n = 0; : : : ; Ncut ; m=−n; : : : ; n; $ =M; N ) in the truncated expansion of the internal and scattered 5elds, and where L=6 is the number of matching points with L ¿ K. (The factor 1/6 results from the fact that the boundary conditions provide 6 equations). The known inhomogeneity bi depends on the orientation, polarisacut ) tion and frequency of the incident wave (and thus on the known expansion coeOcients an;(Nm; $ ). The coeOcient matrix Ai; j contains spherical vector wave functions evaluated at the matching points that lie on the boundary @. Thus the matrix elements Ai; j contain the information about the particle’s geometry and the wave numbers k and k0 in − and + , respectively. The unknown expansion cut ) coeOcients cn;(Nm;cut$) and pn;(Nm; $ are contained in the vector with elements dj . In the earliest application of the PMM [10], the number of matching points was chosen such that the number of equations equals the number of unknowns (K = L), and the system of linear equations given in Eq. (63) was solved with standard techniques. However, the results obtained with

1809

794

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

this approach are often inaccurate and numerically unstable. Also, the accuracy strongly depends on the choice of the matching points. In the so-called generalised PMM (GPMM) [86], one chooses more than twice as many equations as unknowns (L ¿ 2K). In this case the system of linear equations given in Eq. (63) is over-determined, and the unknown expansion coeOcients are determined by minimising the norm 2 K  Ai; j dj − bi (64) j=1 in the least-squares sense. The results obtained by this approach are numerically stable and more accurate than in the ordinary PMM [87]. The large discrepancy in numerical performance between the ordinary PMM and the GPMM is quite puzzling, considering that both methods are essentially based on the same idea. At 5rst glance, the di6erence between both approaches appears to lie in technical aspects concerning the two di6erent numerical methods employed to solve the system of linear equations given in Eq. (63). However, this notion is not correct. Interestingly enough, it turns out that the ordinary PMM and the GPMM are based on fundamentally di6erent mathematical assumptions. It is worth discussing these assumptions and their validity in more detail, since the implications are also relevant for some other electromagnetic scattering methods that employ an expansion of the electromagnetic 5elds. 3.5.2. Underlying assumptions of the PMM and the GPMM By expanding the 5elds on the particle surface according to Eqs. (26)–(28) one assumes that the expansion functions are linearly independent and are complete (but not necessarily orthogonal) on the boundary surface @. Also, one usually assumes that the expansion functions are linearly independent. For the case of scalar spherical wave functions detailed discussions can be found in the literature about the conditions under which these assumptions are satis5ed [88–90]. The scalar spherical wave functions can be shown to be linearly independent and complete on the boundary surface if @ satis5es certain smoothness assumptions (Lyapunov surface). In addition, the linear independence and completeness of the scalar spherical functions of the 5rst kind requires that k 2 is not an eigenvalue of the interior Dirichlet problem. The proof relies on the existence and uniqueness of the solution to the exterior and interior Dirichlet problem of the Helmholtz equation. It is possible to relax the smoothness assumptions and to generalise the proof to spherical vector wave functions. By choosing an equal number of equations as unknowns, the PMM essentially invokes the assumption that the expansions of the 5elds given in Eqs. (26)–(28) converge in any given matching point. In other words, the PMM assumes uniform convergence of the expansions of the 5elds on the non-spherical boundary surface @. This is one of the main assumptions of the so-called Rayleigh hypothesis. Following the de5nition by Millar [88], the Rayleigh hypothesis states that the expansion coeOcients of the 5elds can be determined such that (1) the expansions of the 5elds given in Eqs. (26)–(28) converge uniformly on the boundary surface @; and

1810

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

795

(2) the expansion of the scattered 5eld given in Eq. (27) converges uniformly everywhere in the surrounding medium + ; and cut ) (3) the expansion coeOcients pn;(Nm; $ are independent of the truncation parameter Ncut . The 5rst assumption is made in the ordinary PMM. Several examples are known for which this assumption is clearly not satis5ed [88]. It should be noted that the third assumption in the Rayleigh hypothesis, which has to do with Sommerfeld’s requirement of 5nality [91], is in general not satis5ed for non-spherical boundaries. However, this assumption is not made in the PMM. It should further be noted that the appropriateness of the term “Rayleigh hypothesis” as de5ned by Millar [88] may be debatable in view of the assumptions made in Rayleigh’s original work. However, due to the frequent use of this term in the literature it is probably best to adhere to this term. Incidentally, the early objections to Rayleigh’s work provoked considerable controversy [92–95]. It has been questioned if it was appropriate to expand the scattered 5eld in + in the basis of outgoing wave functions alone. In the immediate vicinity of concave sections on the boundary surface @, thus the argument, the scattered 5eld should in part consist of inwardly propagating waves. However, Millar [88] noted that this is not the origin of the problems related to Rayleigh’s hypothesis. Since Helmholtz’ equation treats the scattering problem in the frequency-domain, one is not dealing with propagating waves at all. The use of “outgoing” wave functions for the expansion of the scattered 5eld is exclusively motivated by the fact that these wave functions satisfy the radiation condition. In the GPMM, an over-determined system of linear equations is solved by least-squares error minimisation. Thus this method is based on the assumption that the expansions of the 5elds given in Eqs. (26)–(28) converge on the boundary surface @ in the least-squares sense. This is a considerably weaker assumption compared to the assumption of uniform convergence contained in the Rayleigh hypothesis. The expansion coeOcients of the 5elds can always be determined such that the expansions of the 5elds converge on the nonspherical boundary surface in the least-squares sense [96]. In addition, the following theorem applies [88,96]. cut ) If the expansion coeOcients pn;(Nm; $ of the scattered 5eld are determined such that the expansion sca of E converges on @ in the least-squares sense, then the same expansion converges uniformly in the surrounding medium + . Thus one can conclude that the GPMM is based on mathematically sound assumptions, and that the solution obtained with this method has mathematically clearly de5ned convergence properties. By contrast, the ordinary PMM is based on assumptions that are in general not justi5ed. However, it is important to keep in mind that a numerical technique is not guaranteed to be numerically eOcient just because it is mathematically sound. Nor is a numerical method necessarily unpractical because a general mathematical proof for its convergence is lacking (A good example for the latter point is the null-5eld method, cf. Section 5). Nevertheless, the case of the PMM and the GPMM provides a good example in which the di6erent experiences gained in numerical applications clearly correspond to what we expect from the investigation of the mathematical justi5cations of either method. The above discussion has shown that the frequently made ansatz of approximating the 5elds by 5nite expansions, which are substituted into the boundary conditions, does not yet contain any assumptions about the convergence of these expansions. Only when we specify a method for determining the unknown expansion coeOcients (e.g. inverting a square matrix in the PMM or performing least-squares error minimisation in the GPMM) then we invoke speci5c assumptions about the kind

1811

796

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

of convergence of the expansions. Di6erent methods will yield di6erent sets of expansion coeOcients for the same expansions, as the comparison of the PMM and the GPMM shows. 3.5.3. Practical issues and applications The early implementations of the PMM/GPMM have been applied to perfectly conducting and lossy dielectric spheres and spheroids [86,97]. An important practical problem with the GPMM is that the use of spherical vector wave functions in the expansions of the 5elds can lead to numerical problems when dealing with particles signi5cantly departing from spherical shape. This problem can be overcome by applying the generalised multipole technique (GMT) originally proposed by Hafner [98,99] and applied in various electromagnetic scattering problems [100–102]. In order to handle more elongated shapes, this method employs multiple spherical expansions about several expansion origins that are located at appropriate positions in the interior region − . The GMT has been applied to the GPMM [87,103,104]. Thus it has been possible to apply the GPMM-GMT to elongated axisymmetric particles, such as hemispherically capped cylinders, 5nite circular cylinders, non-convex “peanut”-shaped particles, and spheroids with axial ratios up to 20 and particle sizes as large as 75' [104]. The GMT has also been applied in methods other than the GPMM (see Section 5). It is also known under di6erent names, such as multiple multipole technique [105], discrete sources method [106], 5ctitious sources method [107], or Yasuura method [108]. A comprehensive review of the GMT literature is given by Wriedt [109]. Extensions of the GMT to anisotropic scatterers [110] and to coated axisymmetric objects [111] exist. 3.5.4. Strengths and drawbacks The well-de5ned convergence properties of the GPMM have already been emphasised, as well as the fact that the radiation condition is automatically satis5ed. In principle, it should be possible to use the GPMM for computing a T matrix, and to use the T matrix for analytically calculating the averaged optical properties of ensembles of randomly oriented particles. This could be achieved in complete analogy to the corresponding approach taken in the SVM [21]. However, to the best of the author’s knowledge analytical orientational averaging has not yet been performed in connection with the GPMM. The GPMM is less Rexible in terms of applications to di6erent particle shapes as compared to methods such as the FDTD, the FEM, or volume-integral equation methods (see Section 4). The origin lies in the expansions of the 5elds in spherical vector wave functions, which become increasingly unsuitable the more the particle’s geometry departs from that of a sphere. As we have seen, elongated particles require the use of specially adapted GPMM implementations with higher computer-code complexity. 4. Volume-integral equation methods 4.1. The volume-integral equation The vector Helmholtz equation (12) is a homogeneous di6erential equation with, in general, non-constant coeOcients (r) (hereafter, it is assumed that the particle is non-ferromagnetic, i.e.

1812

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

= 1). By formally introducing a volume-current density  −i! 0 [ r (r) − 1]E(r) : r ∈ − ; J(r) = 0 : r ∈ + :

797

(65)

Eq. (12) can be rewritten as (∇ × ∇ × 1ˆ − k02 )E(r) = i! 0 J(r):

(66)

Thus the homogeneous di6erential equation with non-constant coeOcients has been transformed into an inhomogeneous di6erential equation with constant coeOcients, which is much easier to solve. A standard solution technique is to 5nd the Green’s function, i.e. the solution to the corresponding di6erential equation with a Dirac-delta-inhomogeneity (∇ × ∇ × 1ˆ − k02 )G0 (r; r ) = 1:3 (r − r )

(67)

subject to the radiation condition. For the inhomogeneous vector Helmholtz equation with a Diracdelta source given in (67), the solution subject to the radiation condition is known analytically. It is given by   exp(ik0 |r − r |) 1  G0 (r; r ) = 1ˆ + ∇∇ : (68) k0 4|r − r | This Green’s function is known as the free-space dyadic Green’s function. Much e6ort has been expended into obtaining series representations of this Green’s function in di6erent coordinate systems [112]. With the free-space dyadic Green’s function, the solution to Eq. (66) [or equivalently, Eq. (12)] can be written as  inc E(r) = E + i! 0 d 3 r  G0 (r; r ) · J(r ); (69) −

inc

where the incident 5eld E is a special solution to the homogeneous Helmholtz equation obtained by setting J = 0 in Eq. (66). Substituting into Eq. (69) the de5nition of J given in Eq. (65) one obtains  inc 2 E(r) = E + k0 d 3 r  [ r (r ) − 1]G0 (r; r ) · E(r ): (70) −

In a review paper, Lakhtakia and Mulholland [113], underlined the common conceptual roots of two volume-integral equation techniques, the method of moments (MoM) and the discrete dipole approximation (DDA). They developed the MoM-formalism from Eq. (70), and the DDA-formalism from the equivalent Eq. (69). 4.1.1. The singular kernel problem Before reviewing the theoretical basics of the MoM and the DDA, an important practical issue needs to be discussed. Inspection of Eq. (68) shows that r→r


G0 (r; r ) → O(|r − r |−3 ):

(71)

The singularity of the free-space dyadic Green’s function requires us to perform the volume integration with utmost care. Suppose we want to evaluate V d 3 r  G0 (r0 ; r ) · J(r ) where V ⊂ − :

1813

798

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

r0 is assumed to be an inner point of V , i.e. there exists a surrounding V0 ⊂ V such that r0 ∈ V0 and r0 ∈ @V0 (see Fig. 3). Then the volume-integral over V can be split up as follows [114]   3    d r G0 (r0 ; r ) · J(r ) = d 3 r  G0 (r0 ; r ) · J(r ) V \ V0

V



+

V0

1 + 2 k0

d 3 r  [G0 (r0 ; r ) · J(r ) − GS (r0 ; r ) · J(r0 )]



@V0

 d r nˆ · − 2 

 r  − r0 · J(r0 ) ; 4|r − r0 |3

(72) (73)

where GS (r0 ; r ) =

1 1 ∇∇ 2  4|r − r0 | k0

(74)

and where nˆ is the normal vector to @V0 . Assuming that V is small, one can invoke the longwave approximation J(r) ≈ J(r0 ) ∀r ∈ V and replace V0 = V . This results in    1 3    d r G0 (r0 ; r ) · J(r ) = M − 2 L · J(r0 ); (75) k0 V  d 3 r  [G0 (r0 ; r ) − GS (r0 ; r )]; (76) M=  L= @V

V

 r  − r0 : d r nˆ · 4|r − r0 |3 2 



(77)

All terms in this expression are regular. Other approximations of the self-term are in use. Lakhtakia and Mulholland [113] discriminate between strong forms of the MoM and the DDA based on the approximation given in Eqs. (75)–(77), and weak forms based on the much cruder approximation  1 d 3 r  G0 (r0 ; r ) · J(r ) = − 2 L · J(r0 ): (78) k0 V The relevance of these approximations for the DDA will be discussed in Section 4.3. 4.2. The method of moments (MoM) A straightforward method for numerically solving Eq. (70) is the method of moments (MoM) [115] which is based    on replacing the volume integral by a discrete sum over 5nite volume elements, − · · · = m Vm · · ·. This is schematically shown in Fig. 4 for the case of a cubical grid. If the volume elements Vm are suOciently small, then one can invoke simplifying assumptions about the variation of (r) and E(r) in each cell. The simplest partitioning scheme uses suOciently small cubical volume elements and assumes that each element is homogeneous (i.e. r is constant throughout each element) and that the long-wavelength approximation can be applied (i.e. E(r) is constant

1814

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

Fig. 3. Fikioris’ method for treating the singular kernel in the VIE: the volume element V is assumed to be small compared to the wave length.

799

Fig. 4. Discretisation of the particle volume in the MoM or the DDA for the example of a cubical grid.

throughout each element) [113]. Thus r (r) ≈ rm

∀r ∈ Vm ; m = 1; : : : ; M;

E(r) ≈ E(rm );

(79)

∀r ∈ Vm ;

(80)

where rm is some pot in Vm , and where M denotes the number of volume cells. Applying the results from Section 4.1.1, the volume integral equation (70) can be transformed into a system of linear equations Einc (rm ) =

M 

Am; k · Ek ;

k=1

Am; k = 1:m; k + (1 − rm ) ·

m = 1; : : : ; M; 

(81)

k02 Mk − Lk

: k = m;

k02 Vm G(rm ; rk )

: k = m:

(82)

that can be solved by using Gaussian elimination or the conjugate gradient method. 4.3. The discrete dipole approximation (DDA) The DDA (also referred to as the coupled dipole method) was originally proposed by de Voe [116] and by Purcell and Penny packer [117]. A review of the popular DDSCAT implementation of Draine and Flatau [118] can be found in Ref. [119]. An excellent review of the formalism and its relation to the MoM can be found in the paper by Lakhtakia and Mulholland [113]. The key issues in the DDA formalism will be reviewed in this section.

1815

800

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

In close analogy to the MoM, the particle volume − is discretised into M small volume cells Vm . Application of this discretisation to the volume-integral equation (69) leads to  M   d 3 r  G(r; r ) · J(r ) = Einc (r) + i! 0 d 3 r  G(r; r ) · J(r ); (83) E(r) − i! 0 Vm

k=1 k =m

Vk

(where r ∈ Vm ). The term on the right hand side M   exc inc d 3 r  G(r; r ) · J(r ) E (r) = E (r) + i! 0 k=1 k =m

Vk

(84)

is interpreted as the 5eld that excites volume cell Vm , because it consists of the incident 5eld and the contributions at r ∈ Vm originating from all other cells Vk , k = m. Invoking the long wave approximation, one can approximate all quantities in Vm by f(r) ≈ f(rm ) = fm where rm is some position vector lying in Vm . Thus Eqs. (83) and (84) in conjunction with Eq. (75) yield   1 exc Em = Em − i! 0 Mm − 2 Lm · Jm ; m = 1; : : : ; M: (85) k0 Further, one obtains from Eq. (65) Jm = i! 0 (1 − r; m )Em :

(86)

Eqs. (85) and (86) in conjunction with Eq. (82) yield Am; m · Em = Emexc

(87)

and thus −1 exc Jm = i! 0 (1 − r; m )Am; m · Em :

(88)

The volume–current density Jm in volume cell Vm can be related to a dipole moment pm =

i V m Jm !

(89)

which becomes with Eq. (88) pm = tm · Emexc ;

(90)

where the polarisability dyad is de5ned as −1 tm = Vm 0 ( r; m − 1)Am; m:

(91)

Thus the e6ect of the exciting 5eld can be interpreted as inducing a dipole moment pm in each discrete volume cell Vm . This is the origin of the term “discrete dipole approximation”.

1816

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

801

Applying the longwave approximation to Eq. (84) yields Emexc = Eminc + i! 0

M 

Vm Gm; k · Jk :

(92)

k=1 k =m

Substitution of Eqs. (89) and (90) into Eq. (92) results in Emexc

=

Eminc

2

+ ! 0

M 

[Gm; k · tk ] · Ekexc

(93)

m = 1; : : : ; M;

(94)

k=1 k =m

or Eminc

=

M 

Qm; k · Ekexc ;

k=1

where Qm; k = 1:m; k − !2 0 [Gm; k · tk ](1 − :m; k ):

(95)

In close analogy to Eqs. (81) and (82) in the MoM, the DDA formalism leads to a system of linear equations that can be inverted numerically by standard techniques, such as Gaussian elimination or the conjugate gradient method. The main di6erence is that the MoM uses as the unknown quantity the total 5elds in the volume cells, whereas the DDA uses the exciting 5elds. As a result, the coeOcient matrices Am; k in the MoM and Qm; k in the DDA are di6erent. However, both methods apply similar numerical procedures for solving the electromagnetic scattering problem. 4.4. Practical issues and applications The computational complexity of the MoM and the DDA critically depends on the method employed for solving the system of linear equations, and on suitable methods for reducing the large number of operations involved in calculating the matrix–vector products. The publically available DDSCAT programme [119] uses a combined fast-Fourier transform conjugate gradient method, which requires that the particle volume is discretised into a periodic lattice. This fast Fourier transform method for performing matrix–vector products can also be used in the MoM, as was done by Lumme and Rahola [120]. An important issue in the DDA is how one relates the relative permittivity r to the polarisability in each cell. The well-known Clausius–Mossotti relation, which Purcell and Penny packer [117] used in their original DDA formulation, is only exact for in5nite cubic lattices in the static limit (! → 0). Various ways of accounting for radiative corrections (! = 0) have therefore been considered by di6erent authors [121–124]. The resulting relations between the polarisability and r and their performances are discussed by Draine and Flatau [118]. Lakhtakia and Mulholland [113] have explained how the various approximations for the polarisability correspond to di6erent ways of treating the singular kernel problem, and how various approximations of the self-term relate to di6erent approximations of the polarisability. This relation can be immediately seen in Eq. (91), in which the polarisability dyad is expressed in terms of the matrix A. The latter is given in

1817

802

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

Eq. (82). It depends on what kind of approximation is used for the self-term, such as for example the strong form [Eq. (75)] or the weak form [Eq. (78)]. The various volume-integral equation implementations that have been proposed and reported [117,118,120,122–127] mainly di6er in the way they treat the self-term, or in what kind of cells they use. The most commonly used cell types are cubical and spherical cells [113]. Volume-integral equation methods have been applied to a large variety of di6erent shapes. For instance, the DDSCAT programme [119] contains geometry implementations for ellipsoids, 5nite circular cylinders, rectangular parallelepipeds, hexagonal prisms, and tetrahedra, as well as for systems of particles, uniaxial particles, and coated particles. New geometries are easy to add. The MoM has, e.g., been applied to irregular particles such as stochastically rough spheres [120]. 4.5. Strengths and drawbacks Volume-integral equation methods have several advantages. They are applicable to arbitrarily shaped, inhomogeneous, anisotropic and optically active particles. The radiation condition is automatically satis5ed, because the free-space dyadic Green’s function satis5es the radiation condition. The computation is con5ned to the volume of the scatterer, which results in fewer unknowns as in many 5nite-di6erence methods. A disadvantage is that the computational accuracy is only improved slowly as the number of cells is increased. Another disadvantage is the need to repeat calculations for each new angle of incidence of the external 5eld. The computational complexity is N ∼O(x9 ) when using Gaussian elimination for solving the system of linear equations, and N ∼O(x3+35 log x) when using conjugate gradient method/fast Fourier transform, where M 5 is the number of iterations in the conjugate gradient method (M being the number of volume cells), and where 0 ¡ 5 ¡ 1. 5. Surface-integral equation methods: The null-3eld method 5.1. The surface-integral equation The starting point for all surface-integral equation methods in general and for the null-5eld method (NFM) in particular is the vector-Green’s identity  d 3 r[X · (∇ × ∇ × Y) − Y · (∇ × ∇ × X)] V

 =

d$nˆ · [Y × ∇ × X − X × ∇ × Y]

(96)

S(V )

which is a general identity that is valid for arbitrary vector 5elds X and Y. From this vector-Green’s identity one can derive the surface-integral equation by substituting X=E; Y=G0 ; V =+ ; S(V )=@, and making use of the Helmholtz equations ∇ × ∇ × E − k02 E = 0;

(97)

∇ × ∇ × G0 − k02 G0 = :3 (r − r ):

(98)

1818

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

803

This yields after reordering factors in the triple-scalar products    i!

Einc (r) + d$(r ) (nˆ+ × H(r )) · G0 (r ; r) + (nˆ+ × E(r )) · (∇ × G0 (r ; r)) c0  =

@

E(r)

: r ∈ + ;

0

: r ∈ − :

(99)

This is the surface-integral equation. The unit vector nˆ+ denotes the outside-pointing normal to @. The surface-integral term represents the scattered 5eld Esca . The lower branch of this equation is called the extended boundary condition or the null-5eld equation. It can be interpreted as an analytic continuation of the incident and scattered 5elds from + into − such that the scattered 5eld exactly cancels the incident 5eld. (Thus the total 5eld in − is equal to the remaining internal 5eld Eint ). Other surface-integral equation methods can be devised by substituting di6erent quantities into the vector-Green’s identity given in Eq. (96). For instance, Kleinman et al. [89] and Ramm [96] substituted an arbitrary element from a set of radiating solution functions to Helmholtz’ equation instead of the free-space Green’s function. However, as shown by Rother et al. [42], the resulting equation is equivalent to the null-5eld equation. 5.2. The null-
     d$(r ) (nˆ+ × H(r )) · G0 (r ; r) + (nˆ+ × E(r )) · (∇ × G0 (r ; r)) c0 @

 =

Esca (r) −E

inc

: r ∈ + ;

(100)

: r ∈ − :

The unknown quantities in Eq. (100) are the scattered 5eld Esca in the surrounding medium + and the 5elds nˆ+ × H|@ and nˆ+ × E|@ on the boundary surface. The NFM proceeds in two steps to determine these unknown quantities. (1) Use the lower branch of Eq. (100), i.e. the null-5eld equation, to determine the unknown boundary-surface 5elds nˆ+ × H|@ and nˆ+ × E|@ in terms of the known incident 5eld Einc . (2) Substitute the results for the surface 5elds into the upper branch of Eq. (100) to determine the unknown scattered 5eld Esca in + . More speci5cally, the incident and scattered 5elds are represented by the expansions given in Eqs. (26) and (27). Likewise, the boundary-surface 5elds are approximated by the

1819

804

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

expansions Ncut  n 

E(kr)|r∈@ =

[cn;(Nm;cutM) Mn;(1)m (kr) + cn;(Nm;cutN) Nn;(1)m (kr)];

(101)

n=0 m=−n

H(kr)|r∈@ =

Ncut  n  n=0 m=−n

[cn;(Nm;cutN) Mn;(1)m (kr) + cn;(Nm;cutM) Nn;(1)m (kr)]:

(102)

In Eqs. (101) and (102) it has been used that the expansion coeOcients of the electric and magnetic 5elds on the boundary surface are not independent. The relation between the expansion coeOcients of the electric and magnetic 5elds on the boundary surface can be derived by applying the vector-Green’s identity to the interior region − , and by using the continuity of the tangential components of the 5elds on the boundary surface [134]. The freespace dyadic Green’s function can be expanded in spherical vector wave functions according to ∞  n  G0 (r ; r) = ik0 (−1)m n=0 m=−n

 ×

Mn;(3)−m (k0 r )Mn;(1)m (k0 r) + Nn;(3)−m (k0 r )Nn;(1)m (k0 r)

: r ¡ r;

Mn;(1)−m (k0 r )Mn;(3)m (k0 r) + Nn;(1)−m (k0 r )Nn;(3)m (k0 r)

: r ¿ r:

(103)

If r ∈ @, the expansion given in Eq. (103) converges uniformly if r lies outside the smallest sphere circumscribing − or inside the largest sphere inscribed into − . The expansions given in Eqs. (26), (27), (101), (102), and (103) are substituted into the surface-integral equation (100). Now one follows the two steps outlined above. First one uses the null-5eld equation to obtain a linear relation between the expansion coeOcients cn;(Nm;cut$) of the boundary-surface 5elds and the expansion cut ) coeOcients an;(Nm; N of the incident 5eld. This resulting linear relation can be written in symbolic matrix–vector form as a = Q · c:

(104)

The elements of the matrix Q are given by Qn; m; M; n
; m
; M = Jn; m; M; n
; m
; N + mr Jn; m; N; n
; m
; M ;

(105)

Qn; m; M; n
; m
; N = Jn; m; N; n
; m
; N + mr Jn; m; M; n
; m
; M ;

(106)

Qn; m; N; n
; m
; M = Jn; m; M; n
; m
; M + mr Jn; m; N; n
; m
; N ;

(107)

Qn; m; N; n
; m
; N = Jn; m; N; n
; m
; M + mr Jn; m; M; n
; m
; N ;

(108)

where Jn; m; M; n
; m
; M = (−1)m



d$nˆ+ · Mn;(3)−m (k0 r) × Mn(1)
; m
(k0 r);

(109)

d$nˆ+ · Nn;(3)−m (k0 r) × Mn(1)
; m
(k0 r);

(110)

@

Jn; m; M; n
; m
; N = (−1)

m



@

1820

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

Jn; m; N; n
; m
; M = (−1)

m

Jn; m; N; n
; m
; N = (−1)m

 @ 

805

d$nˆ+ · Mn;(3)−m (k0 r) × Nn(1)
; m
(k0 r);

(111)

d$nˆ+ · Nn;(3)−m (k0 r) × Nn(1)
; m
(k0 r)

(112)

@

and where mr denotes the refractive index in − . The surface integrals over the boundary surface given in Eqs. (109)–(112) can be evaluated by numerical surface integration. In the second step one uses the upper branch of Eq. (100) to obtain a linear relation between the cut ) expansion coeOcients cn;(Nm;cut$) of the boundary-surface 5elds and the expansion coeOcients pn;(Nm; N of the scattered 5eld, which can be written in symbolic matrix–vector form as p = −RgQ · c:

(113)

The matrix RgQ is obtained by replacing in Eqs. (105)–(108) the elements Jn; m; $; n
; m
; $
by RgJn; m; $; n
; m
; $
, and the elements of the matrix RgJ are obtained by replacing in Eqs. (109)–(112) under the integral the spherical vector wave functions of the third kind Mn;(3)−m (k0 r) and Nn;(3)−m (k0 r) by the corresponding spherical vector wave functions of the 5rst kind Mn;(1)−m (k0 r) and Nn;(1)−m (k0 r), respectively. Inverting Eq. (104) and substitution into Eq. (113) yields the desired expansion coeOcients of the scattered 5eld in terms of those of the incident 5eld p = T · a;

(114)

T = −RgQ · Q−1 :

(115)

where The T matrix in Eq. (115) is computed in the null-5eld method from the Q matrix and the RgQ matrix, the elements of which are computed by numerically evaluating surface integrals over cross products of the spherical vector wave functions. 5.3. NFM: Convergence issues The question of convergence of the NFM solution has been the subject of controversial discussions. In a recent paper, Dallas [135] points out some inconsistencies that can frequently be found in the literature and that contributed to several misunderstandings. For instance, Dallas [135] emphasises the importance of distinguishing between an actual (convergent) in5nite-series representation of the 5elds, and a 5nite sum intended as an approximation of the 5elds. In the present article, the 5elds are consistently represented as a 5nite sum approximation to avoid any misunderstandings. If and how these sums converge (uniform convergence, least-squares convergence) depends on the method employed for determining the coeOcients, as became clear in the discussion of the mathematical justi5cation of the PMM and the GPMM in Section 3.5. The convergence behaviour of the sum approximations of the 5elds can also be di6erent in di6erent spatial regions. Another point to note about the expansions of the 5elds in the NFM is that the question of convergence can easily be confused with convergence issues concerning the expansion of the free-space dyadic Green’s function. The case distinction in the expansion of G0 given in Eq. (103) implies

1821

806

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

sp that Waterman’s NFM solution in the exterior domain + is restricted to a region + outside a circumscribing sphere of the boundary surface @. As noted earlier, the series expansion of G0 sp converges uniformly in + . However, it is incorrect to conclude that the expansion of the scattered sp 5eld converges everywhere in + just because the expansion of G0 does so. The NFM is not based on the assumption that the expansions of the 5elds converge on the boundary surface in the least-squares sense. Therefore the corresponding theorem assuring uniform convergence of the expansion of the scattered 5eld in + cannot be applied. In fact, Waterman [129] tried to avoid invoking any assumptions about the convergence of the expansions of the 5elds on the particle surface, because he doubted that such an expansion would converge at all. Incidentally, Schmidt et al. [136] showed that the NFM expressions for the Q and RgQ matrices can also be derived by substituting the sum approximations of the 5elds into the boundary conditions and applying a projection method, thus assuming that this projection of the sum approximation converges on the boundary surface. However, this approach is clearly not a least-squares approach either. Thus, very little is actually known mathematically about the convergence of the NFM solution. Ramm [96] derived suOciency criteria for the convergence of the NFM solution. But when expanding the 5elds in spherical wave functions, these criteria can only be shown to be satis5ed for spherical boundary surfaces. A major progress has been achieved in the recent work by Dallas [135] who showed that for ellipsoidal boundary surfaces the NFM solution converges in the far 5eld in the least-squares sense. The results obtained by Dallas [135] can be considered the state-of-the-art in our mathematical knowledge about the convergence properties of the NFM solution. On the other hand, it is well-known from numerical experiments for various di6erent shapes and size parameters of the boundary surface that the NFM gives highly accurate far-5eld results. Very little is known about the accuracy and convergence of near-5eld results.

5.4. NFM: Practical issues and applications Perhaps the two most important advantages of the T matrix description of the electromagnetic scattering problem is that one can compute the optical properties of ensembles of randomly oriented particles analytically from the T matrix [20,137,138], and that one can readily exploit geometrical symmetries of particles [139] to drastically reduce CPU-time requirements [140]. However, since these are general advantages of the T matrix formulation that are not restricted to the NFM, these issues will be treated separately and in more depth in Section 8. The computation of the T matrix by means of the NFM involves a numerical matrix inversion, as can be seen in Eq. (115). An important issue in the NFM is how to perform this inversion. Recent studies have shown that the use of quadruple-precision variables [141] and a special LU-decomposition algorithm [142] substantially improve the numerical stability of the NFM and signi5cantly extend its applicability to larger size parameters. To alleviate ill-conditioning problems for highly nonspherical particles, surface-integral equation methods have been devised that are based on using spheroidal vector wave functions instead of spherical vector wave functions [143,144]. The method using multiple expansion origins [98,99], already mentioned in connection with the GPMM [103,104], has also been used in the NFM to extend the applicability range of the method to particles departing signi5cantly from spherical shape. These applications are known as the iterative extended boundary condition method [145,146] or the method of discrete sources [147].

1822

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

807

Another important practical issue concerns the numerical evaluation of the surface integrals in Eqs. (109)–(112) (and analogous expressions for the elements of the matrix RgJ). A numerical quadrature scheme adapted to the particle geometry can reduce the required CPU-time by up to an order of magnitude [140]. The NFM has been applied to several axisymmetric geometries [148], such as 5nite circular cylinders [149], spheroids [32,35,150] and Chebyshev particles [151]. A computer code by Mishchenko [20], described by Mishchenko and Travis [152], for axisymmetric particles is publically available and has been extensively tested and applied [32,35,131,142,149–151,153–155]. The NFM has been extended from homogeneous particles to multilayered [156], clustered [138,156,157] and chiral objects [158], particles on or near a surface [159], composite particles and objects with nonspherical inclusions [160], as well as to di6erent non-axisymmetric geometries, such as triaxial ellipsoids [161,162], cubes [163–165], N -hedral prisms with 3 6 N 6 30 [140], and 6-hedral prisms [166]. 5.5. NFM: Strengths and drawbacks The high accuracy of the NFM makes this method ideally suited for benchmarks computations. The T matrix approximates the full information about a particle’s light scattering and absorption properties at a given wavelength with an accuracy controlled by the truncation parameter Ncut . Thus all quantities of interest, such as the scattering and absorption cross section, the amplitude scattering matrix, and the phase matrix can be expressed in terms of the T matrix. NFM computations do not need to be repeated for each new angle of incidence of the external 5eld. Analytical orientational averaging [20,137] for ensembles of randomly oriented particles is, as mentioned above, one of the most important features of the T matrix formulation in general and of the NFM in particular. For homogeneous symmetric particles, the use of discrete symmetries [139] can speed up NFM computations by several orders of magnitude [140]. However, discrete symmetries could also be exploited in other light scattering methods capable of computing a T matrix. For axisymmetric particles the computational complexity is on the order of N ∼O(x3 ) − O(x4 ). The NFM has been applied to randomly oriented axisymmetric particles with size parameters exceeding 100 [142] and even to circular ice cylinders with size parameters as large as 180 [167], which is unmatched by other light scattering methods based on rigorous theory. On the other hand, for irregular particles the CPU-time and memory requirements of the NFM are signi5cantly higher than for symmetric particles, and the CPU-time requirements for non-axisymmetric particles are comparable to those of the FDTD or the DDA [163] if discrete symmetries are not exploited. The NFM does not lend itself easily to treating inhomogeneous particles, and special implementations are required. The use of spherical vector wave functions for expanding the 5elds results in a decreasing numerical stability when simultaneously considering large size parameters and highly aspherical shapes. Speci5cally adopted implementations are required for circumventing the problems associated with high aspect ratios, such as the iterative extended boundary condition method [146] and the method of discrete sources [147]. Numerical problems are also encountered for particles with high real and/or imaginary parts of the refractive index. Another disadvantage of the NFM is that for non-axisymmetric particles one has to adjust between three and four parameters controlling the convergence of results (the truncation parameter Ncut of the sum over n in the expansions of the 5elds, possibly another truncation parameter Mcut for truncating the sum over m, and the number of polar and azimuthal quadrature points in the numerical evaluation of the surface

1823

808

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

integrals). The large number of convergence-controlling parameters makes it diOcult to implement user-friendly NFM codes for non-axisymmetric particles. 6. Other methods Another method based on rigorous electromagnetic theory is the Fredholm integral equation method [168]. In this approach the internal 5eld in the VIE (70) is subjected to a Fourier transformation into the wavevector space. The resulting Fredholm integral equation is discretised, and the resulting matrix equation is solved numerically. This method allows to eOciently compute scattering by particles of several di6erent refractive indices and for several di6erent directions and polarisations of the incident 5eld. However, this method requires a new implementation for each new particle shape. A time-domain SIE method has been developed [169] for computing transient scattering from dielectric particles. In contrast to the FDTD this method automatically satis5es the radiation condition. Several hybrid methods exist that aim at combining the advantages of two di6erent methods. For instance, Morgan et al. [170] developed a hybrid FEM/NFM method. Yuan et al. [171] devised a combination of the FEM and the MoM. Sheng et al. [77] combined the FEM with SIE methods. The idea of the hybrid methods involving the FEM is to treat complex internal inhomogeneities in the particle by using the FEM, and the “rest” of the problem by the MoM or by SIE methods. Farafonov et al. [172] and Voshchinnikov et al. [173] reported a new light scattering method that combines Farafonov’s formulation of the SVM [13] and the NFM by splitting the 5elds into an axisymmetric and a non-axisymmetric part, introducing Abraham’s potentials for the axisymmetric part, using a superposition of Debye and Hertz potentials for the non-axisymmetric part, and expanding the scalar potentials in spherical wave functions. A system of equations similar to that in the NFM is obtained and solved. 7. Intercomparisons of methods Several intercomparison of di6erent numerical methods for solving the electromagnetic scattering problem are reported in the literature. Most of these reports are primarily concerned with verifying the accuracy of new implementations by checking the results against those obtained with existing and previously tested codes. For instance, Al-Rizzo and Tranquilla [111] compared results for the di6erential scattering cross section obtained with the GPMM and the NFM for oblate spheroids of size parameters up to 20, prolate spheroids of size parameters up to 38, Chebyshev particles of size parameters up to 10, realistic raindrop-shaped particles [174] of size parameters up to 2, and two-layered spheroids of size parameters up to 13. The same authors [104] compared results for the di6erential scattering cross section obtained with the GMT and the NFM for hemispherically capped dielectric cylinders of size parameters up to 11, spheroidal particles with size parameters up to 20, 5nite circular cylinders with a size parameter of 1, and non-convex peanut-shaped particles with size parameters up to 2. Rother and Schmidt [80,81] presented a comparison of DMF and SVM results for the di6erential scattering cross section of spheroids of size parameters 1 and 2. Hovenier et al. [175] compared DDA, NFM, and SVM implementations for prolate and oblate spheroids and for 5nite circular cylinders of an equivalent (equal-volume) sphere size parameters of 5 and a refractive index of 1.5 + 0.1 i. They found good agreement of results for the NFM and SVM computations

1824

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

809

and some deviations for the DDA results. A comparison of the full phase matrix obtained with the SVM and the NFM for spheroidal particles with size parameters up to 10 was presented in Ref. [21]. A comparison of the phase matrix of a cube of a side length l given by k0 l = 10 and a hexahedral prism with hexagonal side length l and prism height h given by k0 l = k0 h = 5 computed with the NFM and the DDA was given in Ref. [140]. The orientation of both the cube and the hexahedral prism with respect to the incident 5eld was 5xed. A similar comparison of NFM and DDA results for an ensemble of randomly oriented hexahedral prisms with k0 l = 4 and k0 h = 2 can be found in Ref. [176]. Lumme and Rahola [120] compared the phase function and the linear polarisation for ensembles of randomly oriented size-shape distributions of completely irregularly shaped stochastically rough spheres computed with the NFM and a VIE code. A critical evaluation of the computational accuracy of di6erent methods and implementations was presented by Voshchinnikov et al. [173], who compared results for the extinction cross section of oblate spheroids of size parameters up to 30 computed with the NFM implementations of Barber and Hill [148] and of Mishchenko [20] and with SVM implementations based on the methods of Asano and Yamamoto [11] and of Farafonov [13]. The main purpose of this comparison was to illustrate the advantages of Farafonov’s method [13] over other methods when considering large, elongated particles. Baran et al. [177] evaluated the relative accuracy of NFM and FDTD computations for monodisperse, randomly oriented hexahedral prisms up to a size of 2kl = kh = 20, where l and h denote, respectively, the side length of the hexagons and the prism height. Computations where performed for refractive indices of 1:30780 + 0:166667 × 10−7 i; 1:29153 + 0:039113i, and 1:27983 + 0:413333i. The authors found for the scalar optical properties that the relative di6erence of results was less than 3% for all refractive indices, and typically less than 1–2% for the cases of low or moderately high imaginary parts of the refractive index. Comparison of the scattering matrix element F11 showed relative di6erences up to 9%, but mainly for scattering angles larger than 150◦ , where the absolute values of the phase function are fairly small. Similar comparisons for the element F12 showed relatively small di6erences between 20◦ and 160◦ . Relative di6erences for the element F22 up to 6.4% where observed, and the di6erences tended to increase with size parameter. The NFM needed between 11 and 35 times less CPU time than the FDTD, and 4 times less memory. In this context it is worth mentioning that in some cases there exist results obtained with methods based on rigorous electromagnetic theory for size parameters that are large enough for allowing a comparison with results obtained with the geometrical optics approximation (GOA). Thus Yang and Liou [49] presented an intercomparison of results obtained with the FDTD and the GOA for 2D scattering by in5nitely long, randomly oriented cylinders with hexagonal cross section with size parameters of k0 l=20 and 40, where l denotes the side length of the hexagon. An intercomparison of NFM and GOA results for the scattering matrix elements F11 and F12 computed for 3D scattering by randomly oriented spheroids and 5nite circular cylinders with size parameters up to 60 was reported by Macke et al. [154,178]. Wielaard et al. [142] compared the full phase matrix computed with the NFM and the GOA for spheroids of size parameters up to 85 and for 5nite circular cylinders of size parameters up to 125. This NFM and GOA intercomparison has recently been extended by Mishchenko et al. [167] to circular ice cylinders with a size parameter of 180. An intercomparison of the CPU- and memory requirements of the DDA, the NFM, and the FDTD has been reported by Wriedt and Comberg [163]. The authors performed computations for cubes of size parameters k0 l up to 15. It was found that the CPU time requirements of all three methods are

1825

810

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

approximately on the same order, whereas the NFM required at least an order of magnitude less storage. Performance intercomparisons of di6erent methods, such as those conducted in Refs. [163,177], are interesting and useful. However, one has to keep in mind that they do not provide an absolute measure for the usefulness of di6erent methods, because the outcome of such comparisons strongly depends on the refractive indices, sizes, and shapes of the particles, and on the speci5c implementations of the di6erent methods. For instance, the NFM code used by Wriedt and Comberg [163] did not make use of the symmetries of a cube (which belongs to the octahedral point-group Oh with order 48). An NFM code that makes full use of a cube’s symmetry can be expected to reduce CPU-time requirements by a factor of 482 ≈ 2300, and memory by almost a factor of 482 (cf. Section 8). The NFM code used by Baran et al. [177] did exploit for the T matrix computations the symmetries of hexahedral prisms (which belong to the dihedral point group D6h with order 24, resulting in a reduction in CPU time by 242 ≈ 600). However, it is unclear if their code also exploited particle symmetries for the analytical orientational averaging procedure (which for hexahedral prisms will reduce computation time for the corresponding subroutines by a factor of 24), and if their code also exploited particle symmetries for reducing memory requirements. On the other hand, the FDTD and the DDA are easier to apply to irregular and inhomogeneous particles. Thus we can expect that an intercomparison of light scattering computations performed for completely irregular and inhomogeneous particles can turn out more favourably for the FDTD and the DDA as compared to the NFM. Moreover, the CPU and memory requirements of di6erent methods can greatly vary with size parameter and refractive index. Thus the results of performance comparisons only apply to the particular range of refractive indices for which the comparisons have been conducted, although they may be useful for indicating certain trends. Finally, a performance-intercomparison is at risk of comparing apples with oranges. Whereas the NFM computes the full information about a particle’s optical properties, the FDTD and the DDA do not. If one is only interested in the scattered 5eld of a particle in a single orientation, than the NFM computes more information than necessary. On the other hand, if one wants to compute the scattered 5eld for many di6erent particle orientations or for an ensemble of randomly oriented particles, then the DDA and the FDTD provide us with too little information, and the computations need to be repeated many times for several di6erent orientations. By contrast, NFM computations only need to be performed once, and the averaging over orientational angles can be done analytically. These aspects are important to consider when interpreting a performance-intercomparison of di6erent methods. Besides comparing computations against independent results, there exist several useful tests for checking the correctness of a calculation. For instance there are numerous properties of the phase matrix [179–183] that can be exploited in testing the correctness of results [184]. These properties have recently been reviewed by Hovenier and van der Mee [185]. Several scattering computations [149,186–188] have used such properties of the phase matrix for automatic testing of results. 8. The T matrix formulation of the electromagnetic scattering problem The traditional description of the scattering problem is based on computing the scattered 5eld for a given incident 5eld. This procedure needs to be repeated for each new angle of incidence or

1826

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

811

each new form of the incident 5eld. By contrast, the T matrix description o6ers a signi5cantly more concise description of a particle’s scattering and absorption properties. The objective in the development of the null-5eld method was to derive an approach suitable for numerically computing the scattered 5eld in the exterior domain + . However, the obtained results are much more comprehensive. The T matrix is a quantity that contains (in approximate form) the full information about a particle’s optical properties at a given wavelength. It depends only on the particle’s size parameter, its shape, its refractive index, and on the particle’s orientation with respect to the coordinate system. The T matrix is independent of the incident 5eld. Thus the T matrix is not just a mere by-product of the NFM, it rather leads to a new and concise matrix formulation of the electromagnetic scattering problem. Historically, this matrix formulation has been developed within the NFM formalism [128–130]. For this reason the NFM has often been referred to as the “T matrix method”. This term may be somewhat misleading, since the NFM is only one among many methods that can be used to compute a T matrix. As mentioned in Section 3.1, a T matrix can be obtained by formulating the SVM in a suitable way [21]. The same principle employed to obtain a T matrix by using the SVM could be applied to the GPMM. The superposition T matrix method [138,156,157] is an approach for computing the T matrix of clusters of particles from the T matrices of the individual cluster components. (This method exploits the translation addition theorem for spherical vector wave functions [189]. A publically available implementation of the superposition T matrix method by Mackowski and Mishchenko exists for clusters of spheres [138]). Rother et al. [42] found a general relation between the surface Green’s function and the T matrix for the exterior-domain problem, which, in principle, allows to compute a T matrix with 5nite-di6erence techniques, for instance. One should therefore strictly discriminate between the null-5eld method and the T matrix formulation of the scattering problem. The fact that the T matrix contains (in approximate form) the full information about a particle’s optical properties leads to a number of major advantages of the T matrix formulation. The two most essential advantages pertain to the treatment of ensembles of randomly oriented particles and to the exploitation of discrete symmetries of the boundary surface. 8.1. Ensembles of randomly oriented particles Once one has determined a particle’s T matrix, it is an easy exercise to compute the scattered 5eld for di6erent incident 5elds from Eq. (114). Each incident 5eld is characterised by its expansion coeOcients contained in the vector a. The scattered 5eld coeOcients contained in the vector p can simply be obtained from a by performing a linear transformation given by the T matrix. Thus, in contrast to methods such as the FDTD or the FEM, numerical computations do not need to be repeated for each new angle of incidence. Conversely, one can also keep the incident 5eld 5xed and rotate the particle. The rotation matrices employed for transforming the T matrix under a rotation [190,191] are the Wigner D-functions that are well-known from the quantum theory of angular momentum [192]. Mishchenko [20] developed a powerful formalism for calculating the phase matrix and the optical cross sections of an ensemble of randomly oriented, axisymmetric particles from the particles’ T matrix by deriving analytical expressions for the orientationally-averaged optical properties. This formalism has been generalised to non-axisymmetric particles [137].

1827

812

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

Fig. 5. Examples of discrete symmetry elements: 5ve-fold rotational symmetry C5 , dihedral symmetry D, reRection plane $h , and inversion centre i.

The derivation of the ensemble-averaged phase matrix of randomly oriented particles makes use of the fact that a particle’s amplitude scattering matrix, and thus a particle’s phase matrix, can be analytically expressed in terms of its T matrix. In conjunction with the transformation of the T matrix under rotations, the integration over orientational angles in the averaging procedure reduces to an integration over Wigner D-functions. Exploiting the Clebsch–Gordan formalism and the orthogonality properties of the Wigner D-functions the integrals can be solved analytically. This results in analytic expressions for the orientationally averaged phase matrix in terms of the T matrix. In a similar way, the orienationally averaged optical cross sections of an ensemble of randomly oriented particles can be expressed analytically in terms of the T matrix. 8.2. Particles with point-group symmetries The second major advantage of the T matrix formulation is that one can readily express discrete symmetries of the boundary surface as symmetry relations of the T matrix [139], which can reduce computation times by several orders of magnitude [140]. Discrete symmetry elements of a scattering particle are coordinate transformations that transform the particle coordinates such that the transformed particle is indistinguishable from the original particle. Examples are indicated in Fig. 5. Discrete symmetries have been extensively used in theoretical chemistry [193], in particular in electronic structure theory and in the classi5cation of molecular vibrations. They have also been applied to potential theory [194,195] and in other surface-integral based formulations of the scalar

1828

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

813

scattering problem [196,197]. In the following, a brief review of the main results from Ref. [139] is given, in which symmetry relations of the T matrix were derived. Let S be a representation of a discrete coordinate transformation in the vector space of the vectors containing the expansion coeOcients of the 5elds. In other words, the expansion coeOcients of the incident and scattered 5elds transform under this coordinate transformation according to a = S · a;

(116)

p = S · p;

(117)

p  = T · a  :

(118)

where Eqs. (114) and (116)–(118) yield T = S−1 · T · S:

(119)

Let’s assume that S represents a coordinate transformation that is a symmetry element of the particle. Then the optical properties are invariant under this transformation. Let’s further assume that the sum approximation of the scattered 5eld given in Eq. (27) has the same symmetry properties under the transformation represented by S as the exact scattered 5eld. Then the T matrix is invariant under this transformation, i.e. T = T, and thus [T; S] = 0;

(120)

where [T; S] = T · S − S · T denotes the commutator of T and S. The close analogy to the description of dynamic symmetries in quantum mechanics is apparent. Invariance (i.e. conservation) ˆ Hˆ ] = 0, where the Hamiltonian Hˆ of an observable Qˆ in a dynamic process is expressed by [Q; describes the dynamics of the system. Analogously, invariance of the optical properties (contained in the T matrix) under a coordinate transformation (represented by the matrix S) is expressed by Eq. (120). Thus the problem of expressing discrete symmetries of scattering particles as explicit symmetry relations of the electromagnetic scattering solution in the T matrix formulation is reduced to 5nding for various symmetry elements representations S in the vector space in which the T matrix operates. This problem can be solved by investigating how the vector wave functions transform under these coordinate transformations [139]. As already mentioned, the derivations in Ref. [139] were based on the assumption that for symmetric particles the sum approximation of the scattered 5eld satis5es the same symmetry relations as the exact scattered 5eld. A derivation of symmetry relations that does not need this assumption has been presented recently [42]. The derivations in Ref. [42] suggest how one can bring the Q matrix; and the RgQ matrix, and the T matrix of non-axisymmetric particles with discrete symmetries into block-diagonal form (where the block matrices correspond to the irreducible representations of the 5nite group), thus reducing computer storage requirements approximately by a factor of the square of the order of the point-group. The set of all symmetry relations of a particle forms a 5nite group, referred to as a point-group [193]. In each point-group one can identify a minimum set of elements from which all other elements can be obtained by combination. Only these elements give rise to independent symmetry relations of the T matrix, as can be seen as follows. Let Q; R, and S represent three symmetry

1829

814

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

elements, i.e. T = Q−1 · T · Q;

(121)

T = R−1 · T · R;

(122)

T = S−1 · T · S;

(123)

and let Q = R · S. Then the symmetry relation (121) can be derived from (122) and (123) by multiplying Eq. (122) from the left with S−1 and from the right with S, and comparing-the resulting equation with Eq. (123). Thus it is only necessary to derive symmetry relations of the T matrix for a relatively small number of symmetry elements in each point-group. It should be noted that the same symmetry relations that hold for the T matrix also hold for the matrices Q; RgQ, and Q−1 . The proof that the matrix Q−1 satis5es the same symmetry relations as the matrix Q becomes almost trivial when reordering the matrix elements according to the irreducible representations of the symmetry group, thus making the matrix block diagonal. This approach is outlined in Ref. [42] for the surface Green’s function description of the scattering problem. The practical implications of exploiting point-group symmetries in the null-5eld method have been systematically investigated in Ref. [140]. For a particle belonging to a point-group of order M0 , the number of Q and RgQ matrix elements that need to be computed can be reduced by a factor of M0 . Also, the area over which the surface-integrals need to be evaluated numerically can be reduced by a factor of M0 . For the computation of the Q and RgQ matrices this results in a reduction in CPU-time on the order of M02 . Numerical tests for rectangular parallelepipeds (M0 = 16) and hexahedral prisms (M0 = 24) could con5rm this [140]. Also, CPU-time requirements for performing the analytical orientational averaging can be reduced by a factor M0 , as demonstrated in Ref. [176]. For future T matrix implementations, a reordering of matrix elements into a block-diagonal structure corresponding to the irreducible representations of a symmetry group can be expected to substantially reduce memory requirements and thus to extend the accessible size parameter regime of T matrix implementations for particles with discrete symmetries. 9. Surface Green’s function (SGF) formulation of the electromagnetic scattering problem The SGF formulation of the electromagnetic scattering problem for separable boundary surfaces is described by Morse and Feshbach [40]. A separable boundary surface is a boundary surface that coincides with a coordinate hypersurface in one of the 11 coordinate systems in which the scalar Helmholtz equation becomes separable. Examples are rectangular parallelepipeds, spheres, spheroids, ellipsoids, or in5nite circular cylinders, which coincide with a coordinate hypersurface in Cartesian, spherical, spheroidal, ellipsoidal, and cylindrical coordinates, respectively. The generalisation of the SGF concept to non-separable boundary surfaces has been introduced by Rother [41]. This generalised SGF formulation o6ers a new, coherent description of the electromagnetic scattering problem that allows us to treat 5nite-di6erence methods and SIE methods under a uni5ed theory. In the following a short review is given of some basic aspects of the SGF formulation. The advantages of having a common framework for treating di6erent light scattering methods is illustrated by a speci5c example.

1830

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

815

The results of Section 3.4, generalised to non-axisymmetric geometries, can be summarised as follows. The exterior solution to the homogeneous Helmholtz equation (∇2r + k 2 )0(r) = 0;

r ∈ +

(124)

subject to the inhomogeneous Dirichlet boundary condition on the boundary surface 0(r) = 0Dirichlet ;

r ∈ @

(125)

and subject to the radiation condition 

@0(r) lim r − ik0(r) = 0 r →∞ @r can be approximated by  0(Ncut ) (r; ; ) =

0



sin  d

 0

(126)

2

(Ncut ) d 0Dirichlet ( ;  )G@ (r; ; ; R( ;  );  ;  );

(127)

where the SGF can be approximated by an expansion in spherical wave functions given by (Ncut ) G@

=

6 Ncut  5   5;6=0 l=−5 l
=−6

(− l )  1 (1) {A@ }− ( ;  )Y5(l) (; ): l; l
; 5; 6 h5 (k0 r)Y6


(128)

1 The expansion coeOcients A− @ of the SGF depend on the geometry of the boundary surface. They can, e.g., be computed with the MoL. The MoL can also be employed to derive a SGF-solution of Helmholtz’ equation that converges on the boundary surface in the least-squares sense [42]. In that case, the solution converges uniformly everywhere in the exterior domain + . However, the SGF can also be derived from the SVM [42]. Rother et al. [42] derived a general T matrix expression in terms of the SGF. From this general expression, Waterman’s surface-integral formulae for the T matrix can be elegantly derived [42]. The SGF is related to the volume Green’s function GV by [42]

G@ (r; r ) =

@GV (r; r ) ; @nˆ−

(129)

where nˆ− denotes the inward pointing normal vector on the boundary surface. The volume Green’s function is a solution to the inhomogeneous Helmholtz equation with a Dirac-delta source in the exterior domain (∇2r + k 2 )GV (r; r ) = :3 (r − r );

r; r ∈ +

(130)

subject to the radiation condition at in5nity, and subject to a homogeneous boundary condition on the boundary surface GV (r; r ) = 0;

r ∈ @:

(131)

A formal solution to the scattering problem given in Eqs. (124)–(126) can be obtained in terms of the volume Green’s function in the following way. Starting from the Green’s identity  [(r)∇2 A(r) − A(r)∇2 (r)]d 3 r (132) +

1831

816

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

  @A(r ) @(r ) (r ) d$(r ) − A(r ) @ n ˆ @ n ˆ − − @

 =

(133)

one substitutes  = 0 and A = GV into Eq. (133). Using Eqs. (124), (125), (130), and (131) leads to  @GV (r ; r) 0(r) = 0Dirichlet (r ) d$(r ): (134) @nˆ− @

From Eq. (134) in conjunction with Eq. (129) one immediately obtains  0(r) = 0Dirichlet (r )G@ (r ; r)d$(r ):

(135)

@

In spherical coordinates, Eq. (135) is identical with Eq. (127). Thus one can see that the SGF formulation of the scattering problem provides a methodical link between the di6erential equation formulation [represented by Eq. (130)] and the surface-integral formulation [represented by Eq. (135)]. This link is provided by Eq. (129). This allows us to treat 5nite-di6erence techniques and SIE methods under the common framework of the SGF formalism. The practical advantage of having such a common formalism can be exempli5ed by considering boundary surfaces with point-group symmetries. The relation between discrete symmetries of the boundary surface and corresponding symmetry properties of the electromagnetic scattering solution in the T matrix formulation has been discussed extensively in Section 8. Such symmetry relations can also be derived for the SGF [42]. Equivalent symmetry relations can be derived for the expansion 1 −1 coeOcients A− @ in Eq. (128) [42]. It turns out that the symmetry relations of the elements of A@ are exactly the same as the corresponding relations of the T; Q, and RgQ matrices [42]. The symmetry relations of the elements of A@ could be used in practice to reduce, e.g., CPU-time requirements in 5nite-di6erence computations. Thus it is in principle possible to exploit discrete symmetries not only in SIE methods, but also in 5nite-di6erence methods. This example illustrates the advantage of having a common formalism for describing 5nite-di6erence and SIE techniques, and how the methodical link between both classes of methods (provided by the SGF) can be used for exploiting certain advantages of one class of methods also in the other class of methods. The SGF concept o6ers the prospect of investigating also other conceptual advantages of SIE methods, such as the ful5lment of the radiation condition and the analytical orientational averaging formalism, with regard to possible applications in 5nite-di6erence methods. A more detailed treatment of the SGF formulation is given by Rother et al. [42]. 10. Concluding remarks In spite of their long history and mathematical simplicity, Maxwell’s equations still pose signi5cant challenges for developing and improving accurate and eOcient numerical solution methods. The overview in this article over some of the most widely used methods has made it clear that no single technique exists (and probably never will exist) that can be applied equally well to all conceivable problems. It is up to the user to decide in each individual case what method is most suitable. The most important issues to consider are the size, the shape, and the dielectric properties of the particle,

1832

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

817

if the scatterer is homogeneous or inhomogeneous, if it possesses certain symmetries, and if one is interested in a single particle in a 5xed orientation or in ensembles of particles with di6erent (or even random) orientations. Some methods, such as the FDTD, FEM, and VIE methods lend themselves easily to treating particles of arbitrary shapes or inhomogeneous particles. Other methods, such as the SVM, NFM and GPMM can more easily account for particle symmetries and are advantageous to use for ensembles of particles in di6erent orientations. For each of these methods there usually exist several implementations of varying degrees of user-friendliness and availability, each having a di6erent applicability range in terms of size parameters, refractive indices and geometries. An excellent starting point for 5nding a suitable electromagnetic scattering code for a given problem are the web-pages maintained by Mishchenko [198], Wriedt [199], and Flatau [200]. In the past decade we have seen major advances in the improvement of numerical methods for solving the electromagnetic scattering problem, but also the development of new formulations, such as the generalised SGF formulation. To further improve the numerical eOciency of existing methods, especially with regard to applications to ensembles of larger particles consisting of many di6erent shapes and in di6erent orientations [120,176,201], will require persistent research e6orts in the future. Besides such practical challenges, numerical techniques for solving the electromagnetic scattering problem still pose unresolved mathematical questions. One example is the question under what conditions convergence of the NFM solution can be guaranteed. Another example is to rigorously prove the convergence behaviour of the MoL for an increasing number of discretisation lines in spherical coordinates. Thus it can be expected that electromagnetic scattering theory will remain a challenging 5eld of active research in the future. Acknowledgements Two anonymous reviewers are gratefully acknowledged for highly constructive comments. The author also wishes to acknowledge Tom Rother for clarifying discussions about the discretised Mie formalism. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

Mie G. BeitrXage zur Optik trXuber Medien, speziell kolloidaler MetallXosungen. Ann Phys 1908;25:377–445. Debye P. Der Lichtdruck auf Kugeln von beliebigem Material. Ann Phys 1909;30:57–136. Stratton JA. Electromagnetic theory. New York: McGraw-Hill, 1941. Born M, Wolf E. Principles of optics. New York: Pergamon Press, 1964. van de Hulst HC. Light scattering by small particles. New York: Wiley, Inc., 1957. Mishchenko MI, Hovenier JW, Travis LD. Light scattering by nonspherical particles. San Diego: Academic Press, 2000. Wriedt T. A review of elastic light scattering theories. Part Part Syst Charact 1998;15:67–74. Flammer C. Spheroidal wave functions. Stanford: Stanford University Press, 1957. Li L-W, Kang X-K, Leong M-S. Spheroidal wave functions in electromagnetic theory. New York: Wiley, 2002. Oguchi T. Scattering properties of oblate raindrops and cross polarization of radio waves due to rain: calculations at 19.3 and 34:8 GHz. J Radio Res Lab Jpn 1973;20:79–118. Asano S, Yamamoto G. Light scattering by a spheroidal particle. Appl Opt 1975;14:29–49. Sinha BP, MacPhie RH. Electromagnetic scattering by prolate spheroids for plane waves with arbitrary polarization and angle of incidence. Radio Sci 1977;12:171–84.

1833

818

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

[13] Farafonov VG. Di6raction of plane electromagnetic wave at a dielectric spheroid. Di6erent Equations 1983;19: 1765–77. [14] Voshchinnikov NV, Farafonov VG. Optical properties of spheroidal particles. Astrophys Space Sci 1993;204: 19–86. [15] Farafonov VG. Optical properties of strongly prolate and oblate spheroidal particles. Opt Spektrosc (USSR) 1990;69:514–8. Z [16] Voshchinnikov NV, Farafonov VG. Numerical treatment of spheroidal wave functions. In: Gustafson BAS, Kolokolova L, Videen G, editors. Electromagnetic and light scattering by nonspherical particles. Adelphi, MD: Army Research Laboratory, 2002. p. 325–8. [17] Chandrasekhar S. Radiative transfer. New York: Dover, 1960. [18] Thomas GE, Stamnes K. Radiative transfer in the atmosphere and ocean. Cambridge: Cambridge Atmospheric and Space Science Series, 2002. [19] Asano S, Sato M. Light scattering by randomly oriented spheroidal particles. Appl Opt 1980;19:962–74. [20] Mishchenko MI. Light scattering by randomly oriented axially symmetric particles. J Opt Soc Am A 1991;8: 871–82. [21] Schulz FM, Stamnes K, Stamnes JJ. Scattering of electromagnetic waves by spheroidal particles: a novel approach exploiting the T-matrix computed in spheroidal coordinates. Appl Opt 1998;37:7875–96. [22] Onaka T. Light scattering by spheroidal grains. Ann Tokyo Astron Observ 1980;18:1–57. [23] Cooray MFR, Ciric IR. Scattering of electromagnetic waves by a coated dielectric spheroid. J Electromagn Waves Appl 1992;6:1491–507. [24] Farafonov VG. Scattering of electromagnetic radiation by two-layer spheroids. Opt Spektrosk 1994;76:79–83. [25] Farafonov VG, Voshchinnikov NV, Somsikov VV. Light scattering by a core-mantle spheroidal particle. Appl Opt 1996;35:5412–26. [26] Voshchinnikov NV. Electromagnetic scattering by homogeneous and coated spheroids: calculations using the separation of variable method. JQSRT 1996;55:627–36. [27] Cooray MFR, Ciric R. Wave scattering by a chiral spheroid. J Opt Soc Am A 1993;10:1197–203. [28] Sinha BP, MacPhie RH. Electromagnetic plane wave scattering by a system of two parallel conducting prolate spheroids. IEEE Trans Antennas Propag 1983;31:294–304. [29] Dalmas J, Deleuil R. Multiple scattering of electromagnetic waves from two in5nitely conducting prolate spheroids which are centered in a plane perpendicular to their axes of revolution. Radio Sci 1985;20:575–81. [30] Cooray MFR, Ciric IR. Scattering of electromagnetic waves by a system of two dielectric spheroids of arbitrary orientation. IEEE Trans Antennas Propag 1991;39:680–4. [31] Nag S, Sinha BP. Electromagnetic plane wave scattering by a system of two uniformly lossy dielectric prolate spheroids in arbitrary orientation. IEEE Trans Antennas Propag 1995;43:322–7. [32] Mishchenko MI, Travis LD, Kahn RA, West RA. Modeling phase functions for dustlike tropospheric aerosols using a shape mixture of randomly oriented polydisperse spheroids. J Geophys Res 1997;102:16,831–47. [33] Schulz FM, Stamnes K, Stamnes JJ. Modeling the radiative transfer properties of media containing particles of moderately and highly elongated shape. Geophys Res Lett 1998;25:4481–4. [34] Schulz FM, Stamnes K, Stamnes JJ. Shape-dependence of the optical properties in size-shape distributions of randomly oriented prolate spheroids, including highly elongated shapes. J Geophys Res 1999;104:9413–21. [35] Zakharova NT, Mishchenko MI. Scattering properties of needle-like and plate-like ice spheroids with moderate size parameters. Appl Opt 2000;39:5052–7. [36] Voshchinnikova NV, Il’in VB, Henning T, Prokopjeva MS. On modelling interstellar polarization. In: Gustafson Z Kolokolova L, Videen G, editors. Electromagnetic and light scattering by nonspherical particles. Adelphi, BAS, MD: Army Research Laboratory, 2002. p. 333–6. [37] Dubkova DN, Voshchinnikov NV, Henning T. Modelling of interstellar extinction with composite dust grains. In: Z Kolokolova L, Videen G, editors. Electromagnetic and light scattering by nonspherical particles. Gustafson BAS, Adelphi, MD: Army Research Laboratory, 2002. p. 61–4. [38] Ciric IR, Cooray FR. Separation of variables for electromagnetic scattering by spheroidal particles. In: Mishchenko MI, Hovenier JW, Travis LD, editors. Light scattering by nonspherical particles. San Diego: Academic Press, 2000. p. 90–130.

1834

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

819

[39] Eide HA, Stamnes JJ, Stamnes K, Schulz FM. New method for computing expansion coeOcients for spheroidal functions. JQSRT 1998;63:191–203. [40] Morse PM, Feshbach H. Methods of theoretical physics. New York: McGraw-Hill Book Company Inc., 1953. [41] Rother T. General aspects of solving Helmholtz’s equation underlying eigenvalue and scattering problems in electromagnetic wave theory. J Electromagn Waves Appl 1999;13:867–88. [42] Rother T, Kahnert M, Doicu A, Wauer J. Surface Green’s function of the Helmholtz equation in spherical coordinates. In: Kong JA, editor. Progress in Electromagnetic Research (PIER), Cambridge, Mass.: EMW Publishing, 2003. in press. [43] TaRove A, editor. Advances in computational electrodynamics: the 5nite-di6erence time-domain method. Boston: Artech House, 1998. [44] Yee SK. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media. IEEE Trans Antennas Propag 1966;14:302–7. [45] Yang P, Liou KN. Finite di6erence time domain method for light scattering by nonspherical and inhomogeneous particles. In: Mishchenko MI, Hovenier JW, Travis LD, editors. Light scattering by nonspherical particles. San Diego: Academic Press, 2000. p. 173–221. [46] Maxwell-Garnett JC. Colours in metal glasses and in metallic 5lms. Philos Trans R Soc A 1904;203:385–420. [47] Umashankar K, TaRove A. A novel method to analyze electromagnetic scattering of complex objects. IEEE Trans Electromagn Compat 1982;24:397–405. [48] Britt CL. Solution of electromagnetic scattering problems using time domain techniques. IEEE Trans Antennas Propag 1989;37:1181–91. [49] Yang P, Liou KN. Light scattering by hexagonal ice crystals: comparison of 5nite-di6erence time domain and geometric optics models. J Opt Soc Am A 1995;12:162–76. [50] Yang P, Liou KN. Finite-di6erence time domain method for light scattering by small ice crystals in three-dimensional space. J Opt Soc Am A 1996;13:2072–85. [51] Mur G. Absorbing boundary condition for the 5nite-di6erence approximation of the time-domain electromagnetic-5eld equations. IEEE Trans Electromagn Compat 1981;23:377–82. [52] Liao Z, Wong HL, Yang B, Yuan Y. A transmitting boundary for transient wave analyses. Scientia Sinica 1984;27:1063–76. [53] Berenger J-P. A perfectly matched layer for the absorption of electromagnetic waves. J Comput Phys 1994;114: 185–200. [54] Berenger J-P. Three-dimensional perfectly matched layer for the absorption of electromagnetic waves. J Comput Phys 1996;127:363–79. [55] Katz DS, Thiele ET, TaRove A. Validation and extension to three dimensions of Berenger PML absorbing boundary condition for FD-TD meshes. IEEE Microwave Guided Wave Lett 1994;4:268–70. [56] Lazzi G, Gandhi OP. On the optimal design of the PML absorbing boundary condition for the FDTD code. IEEE Trans Antennas Propag 1996;45:914–6. [57] Fusco MA. FDTD algorithm in curvilinear coordinates. IEEE Trans Antennas Propag 1990;38:76–89. [58] Fusco MA, Smith MV, Gordon LW. A three-dimensional FDTD algorithm in curvilinear coordinates. IEEE Trans Antennas Propag 1991;39:1463–71. [59] Holland R, Cable VR, Wilson LC. Finite-volume time-domain (FVTD) techniques for EM scattering. IEEE Trans Electromagn Compat 1991;33:281–93. [60] Jurgens TG, TaRove A, Umashankar K, Moore TG. Finite-di6erence time-domain modeling of curved surfaces. IEEE Trans Antennas Propag 1992;40:357–66. [61] Yee SK, Chen JS, Chang AH. Conformal 5nite di6erence time domain (FDTD) with overlapping grids. IEEE Trans Antennas Propag 1992;40:1068–75. [62] Lee JF. Obliquely Cartesian 5nite di6erence time domain algorithm. IEEE Proc H 1993;140:23–7. [63] TaRove A, Brodwin ME. Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell’s equations. IEEE Trans Microwave Theory Tech 1975;23:623–30. [64] Andrew WV, Balanis CA, Tirkas PA, Peng J, Birtcher CR. Finite-di6erence time-domain analysis of HF antennas on helicopter airframes. IEEE Trans Electromagn Compat 1997;39:100–17. [65] Sheen DM, Ali SM, Abouzahra MD, Kong JA. Application of the three-dimensional 5nite-di6erence time-domain method to the analysis of planar microstrip circuits. IEEE Trans Microwave Theory Tech 1990;38:849–57.

1835

820

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

[66] Sullivan DM, Borup DT, Gandhi OP. Use of the 5nite-di6erence time-domain method in calculating EM absorption in human tissues. IEEE Trans Biomed Eng 1987;34:148–57. [67] Videen G, Sun WB, Fu Q. Light scattering from irregular tetrahedral aggregates. Opt Commun 1998;156:5–9. [68] Yang P, Liou KN, Mishchenko MI, Gao B-C. EOcient 5nite-di6erence time-domain scheme for light scattering by dielectric particles: application to aerosols. Appl Opt 2000;39:3727–37. [69] Sun W, Fu Q, Chen Z. Finite-di6erence time domain solution of light scattering by dielectric particles with a perfectly matched layer absorbing boundary condition. Appl Opt 1999;38:3141–51. [70] Sun W, Fu Q. Finite-di6erence time-domain solution of light scattering by dielectric particles with large complex refractive indices. Appl Opt 2000;39:5569–78. [71] Morgan MA, Mei KK. Finite-element computation of scattering by inhomogeneous penetrable bodies of revolution. IEEE Trans Antennas Propag 1979;27:202–14. [72] White RE. An introduction to the 5nite element method with applications to nonlinear problems. New York: Wiley, 1985. [73] Silvester PP, Ferrari RL. Finite elements for electrical engineers. New York: Cambridge Univ Press, 1996. [74] Coccioli R, Itoh T, Pelosi G, Silvester PP. Finite-elements methods in microwaves: a selected bibliography. IEEE Antennas Propag Mag 1996;38:34–48. [75] Mittra R, Ramahi O. Absorbing boundary condition for the direct solution of partial di6erential equations arising in electromagnetic scattering problems. In: Morgan MA, editor. Finite element and 5nite di6erence methods in electromagnetic scattering. New York: Elsevier, 1990. p. 133–73. [76] Volakis JL, Chatterjee A, Kempel LC. Finite element method for electromagnetics. New York: IEEE Press, 1998. [77] Sheng XQ, Jin JM, Song JM, Lu CC, Chu WC. On the formulation of hybrid 5nite-element and boundary-integral methods for 3-D scattering. IEEE Trans Antennas Propag 1998;46:303–11. [78] Pregla R, Pascher W. The method of lines. In: Itoh T, editor. Numerical techniques for microwave and millimeter wave passive structures. New York: Wiley, 1989. p. 381–446. [79] Schiesser WE. The numerical method of lines. New York: Academic Press, 1991. [80] Rother T, Schmidt K. The discretized Mie-formalism for plane wave scattering on dielectric objects with non-separable geometries. JQSRT 1996;55:615–25. [81] Rother T, Schmidt K. The discretized Mie-formalism—a novel algorithm to treat scattering on axisymmetric particles. J Electromagn Waves Appl 1996;10:273–97. [82] Rother T, Schmidt K. The discretized Mie-formalism for electromagnetic scattering. In: Kong JA, editor. Progress in electromagnetic research. Cambridge, MA: EMW Publishing, 1997. p. 91–183. [83] Rother T. Generalization of the separation of variables method for non-spherical scattering on dielectric objects. JQSRT 1998;60:335–53. X [84] Pregla R. About the nature of the method of lines. Arch Elektr Ubertragungstech 1987;41:370–86. [85] Pregla R. Higher order approximation for the di6erence operator in the method of lines. IEEE Microwave and Guided Wave Lett 1995;5:53–5. [86] Oguchi T, Hosoya Y. Scattering properties of oblate raindrops and cross polarization of radio waves due to rain II. Calculations at microwave and millimeter wave regions. J Radio Res Lab Jpn 1974;21:191–259. [87] Al-Rizzo HM, Tranquilla JM. Electromagnetic scattering from dielectrically coated axisymmetric objects using the generalized point-matching technique I. Theoretical formulation. J Comput Phys 1995;119:342–55. [88] Millar RF. The Rayleigh hypothesis and a related least-squares solution to scattering problems for periodic surfaces and other scatterers. Radio Sci 1973;8:785–96. [89] Kleinman RE, Roach GF, StrXom SEG. The null 5eld method and modi5ed Green function. Proc R Soc Lond A 1984;394:121–36. [90] Doicu A, Eremin Y, Wriedt T. Acoustic and electromagnetic scattering analysis using discrete sources. New York: Academic Press, 2000. [91] Sommerfeld A. Partial di6erential equations. New York: Academic Press, 1949. [92] Deriugin LN. Equations for coeOcients of wave reRections from a periodically uneven surface. Dokl Akad Nauk SSSR 1952;87:913–6. [93] Lippmann BA. Note on the theory of gratings. J Opt Soc Amer 1953;43:408. [94] Lysanov IP. Theory of the scattering of waves at periodically uneven surfaces. Sov Phys Acoust 1958;4:1–10.

1836

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

821

[95] Fortuin L. Survey of literature on reRection and scattering of sound waves at the sea surface. J Acoust Soc Am 1970;47:1209–28. [96] Ramm AG. Scattering by obstacles. Dordrecht: Reidel, 1986. [97] Morrison JA, Cross M-J. Scattering of a plane electromagnetic wave by axisymmetric raindrops. Bell Syst Tech J 1974;53:955–1019. [98] Hafner C. BeitrXage zur Berechnung der Ausbreitung elektromagnetischer Wellen in zylindrischen Strukturen mit Hilfe des point-matching Verfahrens. PhD Thesis, Swiss Polytechnical Institute of Technology ZXurich, Switzerland, 1980. [99] Hafner C. The generalized multipole technique for computational electromagnetics. Boston: Artech House, 1991. [100] Kuster N, Ballisti R. MMP method simulation of antenna with scattering objects in the closer-near 5eld. IEEE Trans Magn 1990;26:658–61. [101] Hafner C, Kuster N. Computation of electromagnetic 5elds by multiple multipole method (generalized multipole technique). Radio Sci 1991;26:291–7. [102] Leuchtmann P, Bomholt F. Field modeling with the MMP code. IEEE Trans Electromag Compatbil 1993;35: 170–7. [103] Joo K, Iskander MF. A new procedure of point-matching method for calculating the absorption and scattering of lossy dielectric objects. IEEE Trans Antennas Propag 1990;38:1483–90. [104] Al-Rizzo HM, Tranquilla JM. Electromagnetic wave scattering by highly elongated and geometrically composite objects of large size parameters: The generalized multipole technique. Appl Opt 1995;34:3502–21. [105] Hafner C, Bomholt K. The 3D electrodynamic wave simulator. Chichester: Wiley, 1993. [106] Eremin YA, Sveshinikov AG. The discrete source method for investigating three-dimensional electromagnetic scattering problems. Electromagnetics 1993;13:1–22. [107] Leviatan Y, Baharav Z, Heyman E. Analysis of electromagnetic scattering using arrays of 5ctitious sources. IEEE Trans Antennas Propag 1993;43:1091–8. [108] Kawano M, Ikuno H, Nishimoto M. Numerical analysis of 3-D scattering problems using the Yasuura method. IEICE Trans Electron 1996;E79-C:1358–63. [109] Wriedt T, editor. Generalized multipole techniques for electromagnetic and light scattering. Amsterdam: Elsevier, 1999. [110] Piller NB, Martin OJF. Extension of the generalized multipole technique to three-dimensional anisotropic scatterers. Opt Lett 1998;23:579–81. [111] Al-Rizzo HM, Tranquilla JM. Electromagnetic scattering from dielectrically coated axisymmetric objects using the generalized point-matching technique II. Numerical results and comparison. J Comput Phys 1995;119:356–73. [112] Tai C-T. Dyadic green functions in electromagnetic theory. Piscataway: IEEE Press, 1993. [113] Lakhtakia A, Mulholland GW. On two numerical techniques for light scattering by dielectric agglomerated structures. J Res Natl Inst Stand Technol 1993;98:699–716. [114] Fikioris JG. Electromagnetic 5eld inside a current carrying region. J Math Phys 1965;6:1617–20. [115] Harrington RF. Field computation by moment methods. New York: Macmillan, 1968. [116] DeVoe H. Optical properties of molecular aggregates. I. Classical model of electronic absorption and refraction. J Chem Phys 1964;41:393–400. [117] Purcell EM, Pennypacker CR. Scattering and absorption of light by nonspherical dielectric grains. Astrophys J 1973;186:705–14. [118] Draine BT, Flatau PJ. Discrete-dipole approximation for scattering calculations. J Opt Soc Am A 1994;11:1491–9. [119] Draine BT. The discrete dipole approximation for light scattering by irregular targets. In: Mischenko MI, Hovenier JW, Travis LD, editors. Light scattering by nonspherical particles. San Diego: Academic Press, 2000. p. 131–44. [120] Lumme K, Rahola J. Comparison of light scattering by stochastically rough spheres, best-5t spheroids and spheres. JQSRT 1998;60:439–50. [121] Draine BT. The discrete-dipole approximation and its application to interstellar graphite grains. Astrophys J 1988;333:848–72. [122] Goedecke GH, O’Brian SG. Scattering by irregular inhomogeneous particles via the digitized Green’s function algorithm. Appl Opt 1988;27:2431–8. [123] Dungey CE, Bohren CF. Light scattering by nonspherical particles: a re5nement to the coupled-dipole method. J Opt Soc Am A 1991;8:31–87.

1837

822

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

[124] Draine BT, Goodman JJ. Beyond Clausius-Mossotti: wave propagation on a polarizable point lattice and the discrete dipole approximation. Astrophys J 1993;405:685–97. [125] Iskander MF, Chen HY, Penner JE. Optical scattering and absorption by branched chains of aerosols. Appl Opt 1989;28:3038–91. [126] Hage JI, Greenberg JM. A model for the optical properties of porous grains. Astrophys J 1990;361:251–9. [127] Peltoniemi JI. Variational volume integral equation method for electromagnetic scattering by irregular grains. JQSRT 1996;55:637–47. [128] Waterman PC. Matrix formulation of electromagnetic scattering. Proc IEEE 1965;53:805–12. [129] Waterman PC. Symmetry, unitarity, and geometry in electromagnetic scattering. Phys Rev D 1971;3:825–39. [130] Waterman PC. Matrix methods in potential theory and electromagnetic scattering. J Appl Phys 1979;50:4550–66. [131] Mishchenko MI, Travis LD, Mackowski DW. T -matrix computations of light scattering by nonspherical particles: a review. JQSRT 1996;55:535–75. [132] Mishchenko MI, Travis LD, Macke A. T -matrix method and its applications. In: Mishchenko MI, Hovenier JW, Travis LD, editors. Light scattering by nonspherical particles. San Diego: Academic Press, 2000. p. 147–73. [133] Barber P. Di6erential scattering of electromagnetic waves by homogeneous isotropic dielectric bodies. PhD Thesis, University of California, Los Angeles, 1973. [134] Tsang L, Kong JA, Shin RT. Theory of microwave remote sensing. New York: Wiley, 1995. [135] Dallas AG. On the convergence and numerical stability of the second Waterman scheme for approximation of the acoustic 5eld scattered by a hard object. Technical Report, Dept. of Mathematical Sciences, Univ. of Delaware, No. 2000-7:1–35, 2000. [136] Schmidt K, Rother T, Wauer J. The equivalence of applying the extended boundary condition and the continuity conditions for solving electromagnetic scattering problems. Opt Commun 1998;150:1–4. [137] Khlebtsov NG. Orientational averaging of light-scattering observables in the T-matrix approach. Appl Opt 1992;31:5359–65. [138] Mackowski DW, Mishchenko MI. Calculation of the T matrix and the scattering matrix for ensembles of spheres. J Opt Soc Am A 1996;13:2266–78. [139] Schulz FM, Stamnes K, Stamnes JJ. Point group symmetries in electromagnetic scattering. J Opt Soc Am A 1999;16:853–65. [140] Kahnert FM, Stamnes JJ, Stamnes K. Application of the extended boundary condition method to homogeneous particles with point group symmetries. Appl Opt 2001;40:3110–23. [141] Mishchenko MI, Travis LD. T -matrix computations of light scattering by large spheroidal particles. Opt Commun 1994;109:16–21. [142] Wielaard DJ, Mishchenko MI, Macke A, Carlson BE. Improved T-matrix computations for large, nonabsorbing and weakly absorbing nonspherical particles and comparison with geometrical-optics approximation. Appl Opt 1997;36:4305–13. [143] Hackman RH. The transition matrix for acoustic and elastic wave scattering in prolate spheroidal coordinates. J Acoust Soc Am 1984;75:35–45. [144] Kahnert FM, Stamnes JJ, Stamnes K. Surface-integral formulation for electromagnetic scattering in spheroidal coordinates. JQSRT 2003;77:61–78. [145] Iskander MF, Lakhtakia A, Durney CH. A new procedure for improving the solution stability and extending the frequency range of the EBCM. IEEE Trans Antennas Propag 1983;31:317–24. [146] Iskander MF, Lakhtakia A. Extension of the iterative EBCM to calculate scattering by low-loss or lossless elongated dielectric objects. Appl Opt 1984;23:948–52. [147] Wriedt T, Doicu A. Formulations of the extended boundary condition method for three-dimensional scattering using the method of discrete sources. J Mod Opt 1998;45:199–213. [148] Barber PW, Hill SC. Light scattering by particles: computational methods. Singapore: World Scienti5c, 1990. [149] Mishchenko MI, Travis LD, Macke A. Scattering of light by polydisperse, randomly oriented, 5nite circular cylinders. Appl Opt 1996;35:4927–40. [150] Mishchenko MI, Travis LD. Light scattering by polydispersions of randomly oriented spheroids with sizes comparable to wavelengths of observation. Appl Opt 1994;33:7206–25. [151] Mishchenko MI, Sassen K. Depolarization of lidar returns by small ice crystals: an application to contrails. Geophys Res Lett 1998;25:309–12.

1838

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

823

[152] Mishchenko MI, Travis LD. Capabilities and limitations of a current FORTRAN implementation of the T-matrix method for randomly oriented, rotationally symmetric scatterers. JQSRT 1998;60:309–24. [153] Mishchenko MI. Light scattering by size-shape distributions of randomly oriented axially symmetric particles of a size comparable to a wavelength. Appl Opt 1993;32:4652–65. [154] Macke A, Mishchenko MI, Muinonen K, Carlson BE. Scattering of light by large nonspherical particles: ray-tracing approximation versus T -matrix method. Opt Lett 1995;20:1934–6. [155] Mishchenko MI, Wielaard DJ, Carlson BE. T -matrix computations of zenith-enhanced lidar backscatter from horizontally oriented ice plates. Geophys Res Lett 1997;24:771–4. [156] Peterson B, StrXom S. T matrix formulation of electromagnetic scattering from multilayered scatterers. Phys Rev D 1973;8:3661–78. [157] Peterson B, StrXom S. T matrix for electromagnetic scattering from an arbitrary number of scatterers and representations of E(3)∗ . Phys Rev D 1974;10:2670–84. [158] Lakhtakia A, Varadan VK, Varadan VV. Scattering and absorption characteristics of lossy dielectric, chiral, nonspherical objects. Appl Opt 1985;24:4525–30. [159] Wriedt T, Doicu A. Light scattering from a particle on or near a surface. Opt Commun 1998;152:376–84. [160] Wriedt T, Doicu A. Novel software implementation of the T-matrix method for arbitrary con5gurations of single and clusters of composite nonspherical particles. In: Light scattering by nonspherical particles: Halifax contributions, Adelphi; MD: Army Research Laboratory, 2000. p. 83– 6. [161] Schneider JB, Peden IC. Di6erential cross section of a dielectric ellipsoid by the T -matrix extended boundary condition method. IEEE Trans Antennas Propag 1988;36:1317–21. [162] Schneider J, Brew J, Peden IC. Electromagnetic detection of buried dielectric targets. IEEE Trans Geosci Remote Sens 1991;29:555–62. [163] Wriedt T, Comberg U. Comparison of computational scattering methods. JQSRT 1998;60:411–23. [164] Laitinen H, Lumme K. T -matrix method for general star-shaped particles: 5rst results. JQSRT 1998;60:325–34. [165] Kahnert FM, Stamnes JJ, Stamnes K. Application of the extended boundary condition method to particles with sharp edges: a comparison of two di6erent surface integration approaches. Appl Opt 2001;40:3101–9. [166] Havemann S, Baran AJ. Extension of T matrix to scattering of electromagnetic plane waves by non-axisymmetric dielectric particles: application to hexagonal ice cylinders. JQSRT 2001;70:139–58. [167] Mishchenko MI, Travis LD, Lacis AA. Scattering, absorption, and emission of light by small particles. Cambridge: Cambridge University Press, 2002. [168] Holt AR, Uzunoglu NK, Evans BG. An integral equation solution to the scattering of electromagnetic radiation by dielectric spheroids and ellipsoids. IEEE Trans Antennas Propag 1978;26:706–12. [169] Vechinski DA, Rao SM, Sarkar TK. Transient scattering from three-dimensional arbitrarily shaped dielectric bodies. J Opt Soc Am A 1994;11:1458–70. [170] Morgan MA, Chen CH, Hill SC, Barber PW. Finite element-boundary integral formulation for electromagnetic scattering. Wave Motion 1984;6:91–103. [171] Yuan XC, Lynch DR, Strohbehn JW. Coupling of 5nite-element and moment methods for electromagnetic scattering from inhomogeneous objects. IEEE Trans Antennas Propag 1990;38:386–93. [172] Farafonov VG, Il’in VB, Henning T. A new solution of the light scattering problem for axisymmetric particles. JQSRT 1999;63:205–15. [173] Voshchinnikov NV, Il’in VB, Henning T, Michel B, Farafonov VG. Extinction and polarization of radiation by absorbing spheroids: shape/size e6ects and benchmark results. JQSRT 2000;65:877–93. [174] Pruppacher HR, Pitter RL. A semi-empirical determination of the shape of cloud and rain drops. J Atmos Sci 1971;28:86–94. [175] Hovenier JW, Lumme K, Mishchenko MI, Voshchinnikov NV, Mackowski DW, Rahola J. Computations of scattering matrices of four types of non-spherical particles using diverse methods. JQSRT 1996;55:695–705. [176] Kahnert FM, Stamnes JJ, Stamnes K. Can simple particle shapes be used to model scalar optical properties of an ensemble of wavelength-sized particles with complex shapes? J Opt Soc Am A 2002;19:521–31. [177] Baran AJ, Yang P, Havemann S. Calculation of the single-scattering properties of randomly oriented hexagonal ice columns: a comparison of the T -matrix and the 5nite-di6erence time-domain methods. Appl Opt 2001;40:4376–86. [178] Macke A, Mishchenko MI, Carlson BE, Muinonen K. Scattering of light by large spherical, spheroidal, and circular cylindrical scatterers: geometrical optics approximation versus T matrix method. In: Smith WL, Stamnes K, editors. IRS ’96: current problems in atmospheric radiation. VA: Hampton, 1996.

1839

824

F.M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer 79–80 (2003) 775 – 824

[179] Hovenier JW, van de Hulst HC, van der Mee CVM. Conditions for the elements of the scattering matrix. Astron Astrophys 1986;157:301–10. [180] Hovenier JW, van der Mee CVM. Scattering of polarized light: properties of the elements of the phase matrix. Astron Astrophys 1988;196:287–95. [181] van der Mee CVM, Hovenier JW. Expansion coeOcients in polarized light transfer. Astron Astrophys 1990;228: 559–68. [182] Hovenier JW. Structure of a general pure Mueller matrix. Appl Opt 1994;88:8318–24. [183] Hovenier JW, van der Mee CVM. Basic properties of matrices describing scattering of polarized light in atmospheres and oceans. In: Smith WL, Stamnes K, editors. IRS 96: current problems in atmospheric radiation. VA: Hampton, 1997. p. 797–800. [184] Hovenier JW, van der Mee CVM. Testing scattering matrices: a compendium of recipes. JQSRT 1996;55:649–61. [185] Hovenier JW, van der Mee CVM. Basic relationships for matrices describing scattering by small particles. In: Mishchenko MI, Hovenier JW, Travis LD, editors. Light scattering by nonspherical particles. San Diego: Academic Press, 2000. p. 81–5. [186] Lumme K, Rahola J, Hovenier JW. Light scattering by dense clusters of spheres. Icarus 1997;126:455–69. [187] Hess M, Koelemeijer RBA, Stammes P. Scattering matrices of imperfect hexagonal ice crystals. JQSRT 1998;60: 301–8. [188] Kahnert FM, Stamnes JJ, Stamnes K. On the applicability of models based on simple particle shapes to reproduce Z Kolokolova L, Videen G, editors. the optical properties of small, complex-shaped particles. In: Gustafson BAS, Electromagnetic and light scattering by nonspherical particles. Adelphi, MD: Army Research Laboratory, 2002. p. 151–4. [189] Bruning JH, Lo YT. Multiple scattering of EM waves by spheres. IEEE Trans Antennas Propag 1971;19:378–400. [190] Varadan VV, Varadan VK. Multiple scattering of electromagnetic waves by randomly distributed and oriented dielectric scatterers. Phys Rev D 1980;21:388–94. [191] Tsang L, Kong JA, Shin RT. Radiative transfer theory for active remote sensing of a layer of nonspherical particles. Radio Sci 1984;19:629–42. [192] Varshalovich DA, Moskalev AN, Khersonskii VK. Quantum theory of angular momentum. Singapore: World Scienti5c, 1988. [193] Bishop DM. Group theory and chemistry. Mineola: Dover Publications, 1993. [194] Zakharov EV, Safronov SI, Tarasov RP. Finite-order abelian groups in the numerical analysis of linear boundary-value problems of potential theory. Comput Maths Math Phys 1992;2:34–50. [195] Zakharov EV, Safronov SI, Tarasov RP. Finite group algebras in iterational methods of solving boundary-value problems of potential theory. Comput Maths Math Phys 1993;33:907–17. [196] Zagorodnov IA, Tarasov RP. The problem of scattering from bodies with a noncommutative 5nite group of symmetries and its numerical solution. Comput Maths Math Phys 1996;10:1206–22. [197] Zagorodnov IA, Tarasov RP. Finite groups in numerical solution of electromagnetic scattering problems on non-spherical particles. In: Light scattering by nonspherical particles: Halifax contributions. Adelphi, MD: Army Research Laboratory, 2000. p. 99 –102. [198] http://www.giss.nasa.gov/∼crmim. [199] http://www.t-matrix.de. [200] http://atol.ucsd.edu/∼pRatau. [201] Kahnert FM, Stamnes JJ, Stamnes K. Using simple particle shapes to model the Stokes scattering matrix of ensembles of wavelength-sized particles with complex shapes: possibilities and limitations. JQSRT 2002;74: 167–82.

1840