On Two Numerical Techniques for Light Scattering by Dielectric Agglomerated Structures

Smoke agglomerates are made of many soot sphcres, and their light scattering response is of interest in fire research. The numerical techniques chiefly used for theoretical scattering studies are the method of moments and the coupled dipole moment. The two methods have been obtained in this tutorial paper directly from the monochromatic Maxwell curl equations and shown to be equivalent. The effects of the finite size of the primary spheres have been numerically delineated.

Null vector 0 Null dyadic a Radius of a spherical subregion a,b,di,di Arbitrary vectors flm (1/2) X maximum linear crosssectional dimension of the subregion V" Aim Dyadic kernel for the MOM equations Akk, Akk Self-terms when MOM is used with spherical subregions b Side of a cubical subregion B{x) Magnetic flux density (weber per square meter) D{x) Electric flux density (coulomb per square meter) Cine Electric field intensity amplitude of an incident plane wave E{x) Electric field intensity (volt per meter)

Incident electric field intensity E"
Electric field intensity at x" E^{x) Scattered electric field intensity F^{u,) Far-zone scattering amplitude in the direction «, 0{x,x') Free space dyadic Green's function

Introduction
Electromagnetic scattering problems involving complicated geometries are treated numerically nowadays. Apart from some low-or highfrequency methods [1, 2] and the T-matrix method [3], implementation of most numerical techniques entails a partitioning of the region occupied by the scatterer into may subregions. This is generally true whether a differential formulation is used or an integro differential formalism.
The method of moments (MOM) [4][5][6], as applied to an inhomogeneous dielectric scatterer, is an approach based on the evaluation of an electric field volume integral equation over the region occupied by the scatterer. This region is partitioned into a number of subregions; the electric field in each subregion is represented by a subregional basis function; and the volume integral equation is converted into a set of simultaneous algebraic equations that are solved using standard procedures. The subregions are generally cubes, although the relevant self-terms are usually evaluated as that of equivoluminal spheres.
Whereas the MOM is an actual field formalism, the coupled dipole method (CDM) is based on the concept of an exciting field. The CDM was formulated intuitively in 1969 by Purcell and Pennypacker [7] for dielectric scatterers, although it had by that time a rich history spanning many decades [8]. Conceived from a microscopic viewpoint, the CDM has only a semi-microscopic basis; indeed. the operational basis for applying the CDM to boundary value problems is totally macroscopic [9]. Both the MOM and the CDM were recently extended to bianisotropic scatterers and their respective weak and strong forms were shown to be equivalent [10].
This tutorial paper contains a complete derivation of the MOM and the CDM for the scattering of time-harmonic electromagnetic waves by an inhomogeneous dielectric object possessing an isotropic permittivity, the starting point of the exercise being the monochromatic Maxwell curl equations. A central topic of the paper is the elucidation of the relationship between the MOM, which is widely used in electrical engineering, and the CDM, which is motivated by concepts in atomic physics [7]. Our treatment is directed towards the researcher who is interested in understanding these methods but who may not be a specialist in electromagnetic field theory.
A current topic of considerable interest is light scattering by agglomerated structures made up of individual spheres arranged in a low-density structure. Examples of such structures include smoke agglomerates formed in fires or internal combustion engines; various materials produced in combustion systems including silica and titanium dioxide; and, perhaps, interstellar dust. The earliest relevant analyses [11,12] of light scattering from smoke were based on the Rayleigh-Debye approximation, in which the field exciting any particular primary sphere is taken to be just the field that is actually incident on the agglomerate. Such a procedure neglects multiple scattering effects, and is acceptable for primary spheres with size parameter ( = radius multiplied by the free space wavenumber) less than about 0.2. The typical size parameters for the smoke agglomerates mentioned above lie in the range 0.1 to about 0.25 at visible frequencies.
Both the CDM [13,14] and the MOM [15,16] have been applied to compute scattering from smoke agglomerates in the recent past, since neither technique suffers from the limitations of the Rayleigh-Debye approximation. While the MOM and the CDM are expected to give the same results for infinitesimally small size parameters [17], the methods-as ordinarily used-do not yield identical results as the primary sphere size increases [18,19]. This is because the CDM has been used chiefly in what may be called its weak form, while it is the strong form of the MOM that is in standard usage [10]. The strong form is valid for larger primary spheres because the effect of the singularity of the free space dyadic Green's function is better estimated therein than in the weak form. This becomes clearer in the following sections.

Volume Integral Equation
As is schematically illustrated in Fig. 1, let all space be divided into two mutually disjoint regions. Kin, and Fe«, that are distinguishable from each other by the occupancy of matter. The region Kit is vacuous; hence, The region Fim is filled with an isotropic, linear, possibly inhomogeneous, dielectric continuum with frequency dependent [exp(-/w/)] constitutive equations where e,{x) is the complex relative permittivity scalar. The square root of €r(x) is the local complex refractive index, Imag [w6oe,(jf)] is the local conductivity, and CO is the circular frequency.
There is no requirement that Kim be a simplyconnected convex region. Sharp corners and cusps should be absent, it being preferable that the boundary surface that separates Knt from K^ be at least once-differentiable everywhere to enable the unambiguous prescription of a unit normal at eveiy point on it. Furthermore, the maximum linear extent of Kint must be bounded so that only the region K", extends out to infinity in all directions.
The monochromatic Maxwell curl equations, in the absence of any externally impressed sources, are given in Kn as  In Eq. (5b), the volume electric current density J(x) = 0 for x G K.", but which implies that/(r) is not merely a mathematically convenient quantity for dielectric scatterers. It is enough that we look for the solution of the differential equations (5a, b) in terms of only the electric field, since the magnetic field everywhere may be obtained from Eq. (5a) if the electric field is known everywhere: /f(jr) = Vx£(x)//w^. On taking the curl of both sides of Eq. (5a), and then substituting for VxH(x) from Eq. (5b) in the result, we get In this fashion a volume electric current density has effectively replaced the dielectric matter occupying the region Vi" [20,. The volume electric current density/(AT) defined by Eq. (5c) is not a fictitious entity in the present context. We must remember that, while the region Vox is vacuous, the region VM is occupied by dielectric matter. The monochromatic polarization current density in this dielectric matter is given by while the conduction current density is given by It is now easy to see that (7) where /to = co(/uoeo)"^ is the free space wavenumber. We take cognizance of the fact that when J{x) is set equal to zero everywhere, Eq. (7) reduces to the celebrated vector Helmholtz equation. As shown in Appendix A, it follows from Eq. (7) that the electric field E{x) is the solution of the volume integral equation or, equivalently, is the free space dyadic Green's function and I is the identity dyadic. The field E;"c(x) is the solution of Eq. (7) when J(x)sO for all x E Fi", + Fen, and is the electric field existing in V-M + Vat if €,(x) = l for all X £ Fini. The standard radiation conditions are satisfied by the right sides of Eqs. (8a, b) as ko\xl^«> [21].

Intermediate Remarks
Equation (8a) is utilized in setting up the CDM, while the integral Eq. (8b) is solved in the MOM, and it is this commonality that partially begets the algorithmic equivalence of the two techniques. In setting up the solutions of Eq. (7) in the form of the integral Eqs. (8a,b), we have brought in two significant topics that need some rumination right now: (i) dyadics and (ii) the free space dyadic Green's function.
Dyadics are as American as apple pie, being the brainchildren of Josiah Willard Gibbs. In 1884, Gibbs circulated a pamphlet introducing the concept and nomenclature of dyadics. Mathematics books with dyadic notation were written every so often earlier in this century, but most mathematicians appear to have eventually discarded dyadics in favor of tensors. In electromagnetics, though, dyadic notation has been used with great profit, the books by van Bladel [22], Fedorov [23] and Chen [24] being immensely popular. A short exposition on dyadics was brought out by Lindell [25] in 1981, and is much recommended to the interested reader.
A dyadic serves as a linear mapping from one vector to another vector; thus, a dyadic D is a mapping from a to 6 given by * = D • a. This property leads to the idea of a dyad that is composed of two vectors, i.e., Dn=didi. It follows that Di2 • a =di(dz • a) and a • D12 = (a • di)d2 are vectors, and Di2Xa =di(d2Xa) and a xDi2=(a y.di)di are dyads, and the oft-used appellation bivectors for dyads appears justified [26]. Because a dyadic can be written as the sum of dyads, the general representation of a dyadic is the sum D = Xkn, Dkmdi4m [27].
The identity dyadic I is such that D • I = D = 1 • D, and the null dyadic 0 is such that D* 0 = 0 = 0-0. The simplest antisymmetric dyadic is « x I, where « is any vector of unit length. Even vector differential operators can be thought of as dyadics; thus, the curl operator is written as V x I, and the divergence operator as V • I, in dyadic notation. It is not possible for us to go through all the wonderful properties of dyadics than can be exploited in theoretical electromagnetics research, so we refer the interested reader to the compendiums in the books by Chen [24], van Bladel [22], and Varadan et al. [28]. The main feature of computational significance is that, as dyads are bivectors, all dyadics used in this paper can be thought of as 3 x 3 matrices.
The second issue concerns the singularity of the dyadic Green's function G(x, x') used in Eqs. (8a,b) and defined in Eq. (8c): the factor exp (//to |jir-x'l)/47rlx-jf'l goes to infinity as lx-x'l-*0. When evaluating G(x, x'), the second derivatives therefore have to be carefully handled. As shown by van Bladel [29], the classical procedure leads one to write where P.V. stands for "principal value." When computing the P.V. integral on the right side of Eq. (9), an infinitesimally small spherical region centered about the point x'=x is excluded from the domain of integration. Computing the P.V. integral is not problematic because [30] G{x, is not singular when the source point x'and the field point x do not coincide; here X=x-x' and ux=XI\X\. Yaghjian [31] has modified the right side of (9) to an expression wherein the excluded infmitesimal region need not be spherical. We will use an alternative approach to evaluate the integral (9), as shown in Appendix 2. In this approach [6,10,32], the region of integration is split into two regions: one region includes the point x'=x in its interior, while the second one is the remainder. Use is made of Gauss' theorem and the Green's function for Poisson's equation to determine this integral. This procedure is attractive as it places very little restrictions on the shape and the size of the region surrounding the point x' =x.

The Method of Moments
The method of moments is a general mathematical procedure for transforming a linear operator equation into a set of simultaneous algebraic equations, and the interested reader is referred to [4], Chap. 1, for a simple introduction. We are, however, confining ourselves here to the application of MOM to electromagnetic scattering problems.
Although the MOM has grown increasingly sophisticated in the last decade [6,33], a simple version [5] suffices for the easy conversion of the integral equation (8b) into a set of simultaneous algebraic equations. We begin by partitioning the scatterer region Fim into simply connected subregions Vm, l£m <M, each bounded by a closed surface 5m, so that Eq. (8b) can be transformed to The main features of this partitioning scheme are as follows: (i) Each subregion Vm is modeled as being homogeneous such that (ii) The maximum linear cross-sectional extent 2flm of Vm is such that flm/A <0.1 and flmler.m"^l/A <0.1; that is, the dimension a" is no more than a tenth of the wavelength A in the exterior region K", as well as of the wavelength in the subregion Vm. (iii) The surface 5^ that bounds the subregion Vm is sufficiently smooth so as to be at least once differentiable, which enables the unambiguous prescription of a normal at any point on Satisfaction of these three conditions, in practice, means that the union SmVm of the subregions is only approximately coincident with the scatterer region VM. Condition (i) can lead to an artificial material discontinuity across the interface of two adjacent subregions, therefore the simple MOM algorithm given in this section works best when. adjacent subregions do not differ widely in their permittivities. Condition (ii) ensures that the spatial variations of the electromagnetic field inside each subregion are small enough so that each subregion can be thought of as a dipole scatterer [34], though different authors use somewhat different upper bounds on the subregional size [35,36]. Condition (iii) is mostly ignored by MOM users, their usual practice being to use cubical subregions. Thus, the adequacy and the accuracy of the MOM-and the CDM-results depend in large measure on the adequacy of the partitioning scheme.
Two ancillary aspects of the partitioning scheme need to be stated for the simple MOM algorithm. First, the incident field must have slow enough spatial variations that it may be considered almost spatially constant over any subregion. Second, not only the permittivity but also the actual field is assumed spatially constant in each subregion. Together, these two assumptions constitute a long wavelength approximation [37], whose consequences are exploited in the MOM as well as in the CDM.
Let the volume of Vm be denoted by Vm = fSIy d'jf, and let Xm denote a distinguished point (such as the centroid) lying inside Vm. On setfingx=jr* in Eq. (11) we get the relation follows from Eq. (A2-6a) with Aim =xt -Xm. In this straightforward and simple version of the MOM, wherein the electric field has been assumed to be piecewise constant over the scatterer region KVK, Eq. (14a) has to be solved for the three cartesian components of all E^ for specified £inc(x).
Bearing in mind that all of our dyadics can be easily thought as 3 x 3 matrices, one can solve Eq. (14a) using a variety of matrix manipulation techniques. The Gauss elimination procedure [38] is simple but places a heavy demand on computer memory. Much less computationally intensive is the conjugate gradient method [39], which is now being heavily employed for MOM calculations [40].
Once all £m have been thus calculated, the scattered electric field in Vc» can be computed as the sum which expression follows from Eq. (8b). Now, from Eq. (A2-6a), we have the asymptote [41] Tend,"^.,",.

The Coupled Dipole Method
The heart of the MOM is Eq. (13b) which involves the electric field Ek that is actually present at Xk. However, in the CDM we have to consider the electric field that excites the subregion Vk, each subregion being explicitly modeled as an electric dipole moment. In order to obtain an expression for this exciting field, we return to Eq. (8a), and rewrite it as after partitioning the scatterer region precisely in the same manner as was done for the MOM solution. We know that Ei"c(x) is the field in the absence of the scatterer and that the quantity is an electric field whose source clearly lies in the subregion V^, mr^k. Thus, the whole of the right side of Eq. (20) is a composite electric field with multiple sources, but none of the sources lies in Vk. This composite electric field
Since the left side of Eq. (20) must be equal to its right side, it follows that

<'^E^,(x)^E(x)-ia>fio Uiy, d^jc' {G(x,*')-/ {x')}\
X e n. (21b) On setting x =Xk in Eq. (21b), and completing the volume integration therein using Eq. (A2-5a), we get the expression We return to Eq. (21a) at this stage, set x =Xk therein, and do the volume integrations; the result is the vector dyadic relation   (28) which relation follows from Eq. (8a), Equation (18) still applies for the scattered electric field in the far zone, but the far-zone scattering amplitude is now given by

Scattering and Absorption
The presence of a scatterer has two observable consequences insofar as the conservation of electromagnetic energy is concerned. A portion of the energy radiated by the source of EiJ^x) is scattered in all directions by the scatterer. Another portion is absorbed by the scatterer and converted into other forms such as heat. All calculations pertaining to the scattered and the absorbed powers at a given circular frequency a are made in terms of power averaged over the time period ITTICO. Because with dn(us) = sinei dft d<fh, as is customary in spherical coordinates, and the asterisk denoting the complex conjugate. After putting Eqs. (18) and (30) in Eq. (31a), we obtain Fs^*(us); (31b) consequently, we are able to ascertain the total time-averaged scattered power as the integral /'sca = (l/2Tio) I d^ I d ft sin 0s

F^.(us)-F^*(us).
(31c) Unless the scatterer material is intrinsically lossless, there is absorption of electromagnetic energy in Fint. Because the scatterer is simply a dielectric object here, the time-averaged power absorbed in Fint may be computed as the volume integral from Poynting's theorem for monochromatic fields [24,Chap. 2]. After using Eq. (2a) on the right side of this relation, we get the result /'ab, = Real [(/w/2) (32b) Finally, after using Eq. (12) as well as the longwavelength approximation made heretofore, we are able to reduce this expression to the sum Quite often, one is interested in the extinction of the plane wave where fine carries the units of volts per meter and Mine is a dimensionless unit vector such that emc"«inc = 0. In this case, the total time-averaged power extinguished by the presence of matter in V\M can be estimated using the forward scattering amplitude as [22,Chap. 8]: Substitution of Eqs. (29a) and (36a) in this relation gives us whence, for use with the MOM.

Spherical and Cubical Subregions
Although spheroidal and ellipsoidal subregions have been used [43,44], it is commonplace in literature to have cubical or spherical subregions. CDM users are more comfortable with spherical subregions, while MOM users are fond of cubical ones. Cubes and spheres have the same dyadic L, and in many MOM codes the M dyadic of a cube is estimated as that of an equi-voluminal sphere; see Appendix C. Without any particular loss of generality therefore, we take the subregions to be spherical and of identical radii in the remainder of this paper.
Let the subregion Vk be the sphere of radius a with its center at jc*. As a result, the volume Both Avk and Akk should be called self-terms; instead, it seems only Akk has been accorded that honor in some MOM papers, e.g., [16][17][18][19]45].
The CDM self-term is somewhat obscure, being hidden as the polarizability dyadic tt of Eq. (24b). In the present instance, the polarizability dyadic reduces to a scalar because the subregion is spherical; hence, where and Tk = ak/(1+Akk/Akk), Ok^Aira^ eo(er.* -l)/(er,* + 2) (42c) is the Mossotti-Clausius polarizability [46,47] of the electrically small dielectric sphere. Let A^ofl < < 1 in Eq. (41c) and Au be evaluated correct to order fco^a^. Then, and we observe that the {2iko^/3)ak/4Treo term in the denominator of the right side of Eq. (43) is the radiative reaction term used by Draine [48] in his CDM formulation. More often, Akk is evaluated correct only to the first order in k(fl, leading to Tk = ak [7], and thereby giving rise to the semimicroscopic flavor of this numerical approach [9].

Strong and Weak Forms
The various approximations that can be made for the self-terms lead us to the strong and the weak forms of the MOM and the CDM [10]. In the strong forms (S-CDM & S-MOM), Eq. (A2-5a) is used to estimate the self-term {m=k)'\x\ Eqs. (13a) and (21b). In the Weak forms (W-CDM & W-MOM), Eq. (A2-7) is used in place of Eq. (A2-5a) for this estimation. For isotropic dielectric scatterers, the W-CDM is exemplified by Purcell and Pennypacker [7] using spherical subregions, and the S-CDM has become recently available [17] for the same problem. The S-MOM has been used for isotropic dielectric scatterers for many years [6,33], but the idea of W-MOM is of more recent provenance [17].
The W-MOM corresponds exactly to the W-CDM, while the S-MOM corresponds exactly to the S-CDM. When all subregional volumetric capacities are very small, the S-MOM/S-CDM effectively transmutes into the W-MOMAV-CDM. Generally stated, therefore, it can be surmised that the scatterer region V\m must be broken up into a larger number of subregions when the W-MOM/ W-CDM is used than if the S-MOM/S-CDM is used. Comparison of S-MOM results with the W-CDM results, with identical partitioning of the scatterer region, [e.g., 16,18,19], in some instances may be like comparing apples with oranges. The difference between S-CDM/S-MOM and W-CDM/ W-MOM is solely due to the inclusion or the exclusion of the dyadic M* in the expression for A«, therefore computational time is marginally increased and the memory requirement negligibly, when one shifts from W-CDM/W-MOM to S-CDM/S-MOM.

Spherical Subregions and Finite-Size Effects
An assessment of the MOM and the CDM for isotropic dielectric scatterers with spherical subregions is now in order. To facilitate such a comparison, we reiterate that indicates that the weak forms become increasingly inadequate as the size parameter of the subregion increases. Equations (46a) and (46b) tell us that strong forms incorporate finite-size effects meaningfully, while the weak forms do so trivially. This conclusion is reinforced by Figs. 2-4 that show plots of r/a versus the normalized radius k^ of an electrically small sphere for 6,= 1.5, 2.5 and 4.0. In these figures we have ensured that the maximum value of Arofller"^ is about 0.5. We observe-not surprisingly -that T is complex while a is purely real for these input parameters. Furthermore, in all three figures the ratio Ir/al lies between 1.02 and 1.03 when A:uflle,"^s0.5; while T/a = l for k<fi =0, of course. Similar conclusions can be drawn from Fig. 5, wherein the calculations of r/a have been reported for a lossy dielectric sphere with e, = 3.75-t-0.25/. These figures highlight the understanding on finitesize effects drawn from Eqs. (46a,b), and it follows that coarser partitioning of Ki", may be acceptable when S-CDM/S-MOM is used than when W-CDM/ W-MOM is used. These thoughts can be validated by careful examination of the graphs recently published by Ku [19].  ng. 5. Same as Fig. 2, but t, = 3.75 + 0.25i.
The strong forms discussed above are not the only way of introducing finite-size effects, there being at least two more CDM algorithms available for that purpose. As alluded to earlier, Draine [48] replaced the Mossotti-Clausius polarizability by [D-CDM] (47) However, TD is fully contained in T, as we have seen in Eq. (43); and TD approaches a as the size parameter kaa becomes vanishingly small.
The second source for the incorporation of finite-size effects has been the Lorentz-Mie-Debye analysis for scattering by a dielectric sphere [46]. The electromagnetic field scattered by any object can be described in terms of a multipole expansion [49]. With inspiration from Doyle [50], Dungey and Bohren [51] introduced finite-size effects in the CDM by using the lowest order Mie coefficientcorresponding to the electric dipole term of the mutipole expansion -for the polarizability in place of the Mossotti-Clausius polarizability a. Thus,
Draine [48] and Dungey and Bohren [51] concluded from their numerical investigations that D-CDM and DB-CDM, respectively, generally provide scattering results superior to those from W-CDM. This does not come as a surprise because the self-terms in W-CDM (or W-MOM) are estimated with the least accuracy. On the other hand, although it is difficult to provide general enough conclusions for the adequacy of either D-CDM or DB-CDM vis-a-vis that of the S-CDM/ S-MOM, it is safe to state that any claims of superiority-based purely on the estimation of some gross parameter, such as the total scattering crosssection-are debatable. Indeed, the only good way of deciding on the superiority of S-CDM, D-CDM or DB-CDM is by making calculations for the specific problem under consideration: The scatterer region should be subdivided into different numbers of subregions, and computed data from the various algorithms for different partitioning schemes compared for the property of interest [52][53][54].
The refractive index of soot at visible frequencies has been measured by a variety of experimental techniques, and has been found to be dependent on the source materials [55]. We conclude with calculations of r/a and TDBIU for er = 2.40 + 2.38j, corresponding to a complex refractive index of 1.7 + 0.7/, in Fig. 6. It is not difficult to gather from this figure that conclusions drawn in the previous paragraphs of this section apply to light scattering by soot agglomerates.

Appendix A. The Volume Integral
Equation for the Electric Field that depends on ihe forcing function, J(x) being the forcing function here, and the complementary function that holds when the forcing function is The free space Green's function G(x,x') satisfies identically zero everywhere. The right side of Eq. the dyadic differential equation

VxVxG(x,x')-ko^G(x,x') = \S(x-x'),
where S(x-x') is the Dirac delta function. Since J(x') is not a function of x, it follows that

Appendix B. The Self-Integral
The integral under consideration is given as

(Al-6)
The solution of every linear differential equation can be divided into two parts: the particular solution As shown by Fikioris [56], we can transform Eq. (A2-1) to is an auxiliary dyadic related to the Green's function for Poisson's equation [10]; Vo is an exclusionary region bounded by the surface So, as shown in Fig. 7; u" is the unit normal to 5o, pointing away from Vo, at the point x'^Sa. The exclusionary region Vo should be small but not necessarily infinitesimal, and it must be wholly contained within V. Moreover, there is no requirement that So be a miniature copy of S. The use of the long-wavelength approximation b{x)sb(xo) for all x £ Freduces Eq. (A2-2) to a(xo)=nSy_y"d\' [G(xo,x')'b(xo)] +