Discrete dipole approximation for scattering calculations

Vol. 11, No. 4/April 1994/J. Opt. Soc. Am. A 1491 Discrete-dipole approximation for scattering calculations Bruce T. Draine Princeton University Observatory, Princeton, New Jersey 08544-1001 Piotr J. Flatau Scripps Institution of Oceanography, University of California, San Diego, La Jolla, California 92093-0221 Received July 20, 1993; revised manuscript received October 12, 1993; accepted October 18, 1993 The discrete-dipole approximation (DDA) for scattering calculations, including the relationship between the DDA and other methods, is reviewed. Computational considerations, i.e., the use of complex-conjugate gradient algorithms and fast-Fourier-transform methods, are discussed. We test the accuracy of the DDA by using the DDA to compute scattering and absorption by isolated, homogeneous spheres as well as by targets consisting of two contiguous spheres. It is shown that, for dielectric materials (Iml c 2), the DDA permits calculations of scattering and absorption that are accurate to within a few percent. 1. INTRODUCTION The discrete-dipole approximation (DDA) is a flexible and powerful technique for computing scattering and absorp- tion by targets of arbitrary geometry. The development of efficient algorithms and the availability of inexpen- sive computing power together have made the DDA the method of choice for many scattering problems. Bohren and Singham' reviewed recent advances in calculations of backscattering by nonspherical particles, including DDA results. In this paper we review the DDA, with particu- lar attention to recent developments. In Section 2 we briefly summarize the conceptual basis for the DDA. The relation of the DDA to other finite- element methods is discussed. DDA calculations require choices for the locations and the polarizabilities of the point dipoles that are to repre- sent the target. In Section 3 we discuss these choices, with attention to the important question of dipole polar- izabilities. Criteria for the validity of the DDA are also considered. Recent advances in numerical methods now permit the solution of problems involving large values of N, the num- ber of point dipoles. These developments, particularly the use of complex-conjugate gradient (CCG) methods and fast-Fourier-transform (FFT) techniques, are reviewed in Section 4. We illustrate the accuracy of the method in Section 5 by using the DDA to compute scattering by a single sphere and by targets consisting of two spheres in contact and by comparing our results with the exact results for these geometries. Although these are only examples, they make clear that the DDA can be used to compute highly accurate results for targets with dielectric con- stants that are not too large (Iml c 2). A portable FORTRAN implementation of the DDA, the program DDSCAT, has been developed by the authors and is freely available.2 2. WHAT IS THE DISCRETE-DIPOLE APPROXIMATION? A. Approximation Given a target of arbitrary geometry, we seek to calculate its scattering and absorption properties. Exact solutions to Maxwell's equations are known only for special geome- tries such as spheres, spheroids, or infinite cylinders, so approximate methods are in general required. The DDA is one such method. The basic idea of the DDA was introduced in 1964 by DeVoe,3'4 who applied it to study the optical prop- erties of molecular aggregates; retardation effects were not included, so DeVoe's treatment was limited to ag- gregates that were small compared with the wavelength. The DDA, including retardation effects, was proposed in 1973 by Purcell and Pennypacker,5 who used it to study interstellar dust grains. Simply stated, the DDA is an approximation of the continuum target by a finite array of polarizable points. The points acquire dipole moments in response to the local electric field. The dipoles of course interact with one another via their electric fields,5' 6 so the DDA is also sometimes referred to as the coupled dipole approximation.7'8 The theoretical basis for the DDA, in- cluding radiative reaction corrections, is summarized by Draine.6 Nature provides the physical inspiration for the DDA: in 1909 Lorentz showed9 that the dielectric properties of a substance could be directly related to the polarizabil- ities of the individual atoms of which it was composed, with a particularly simple and exact relationship, the Clausius-Mossotti (or Lorentz-Lorenz) relation, when the atoms are located on a cubic lattice. We may expect that, just as a continuum representation of a solid is appropriate on length scales that are large compared with the interatomic spacing, an array of polarizable points can accurately approximate the response of a continuum 0740-3232/94/041491-09$06.00 ©1994 Optical Society of America B. T. Draine and P. J. Flatau

1492 J. Opt. Soc. Am. A/Vol. 11, No. 4/April 1994 target on length scales that are large compared with the interdipole separation. For a finite array of point dipoles the scattering problem may be solved exactly, so the only approximation that is present in the DDA is the replacement of the continuum target by an array of N-point dipoles. The replacement requires specification of both the geometry (location rj of the dipoles j = 1, ... , N) and the dipole polarizabili- ties aj. For monochromatic incident waves the self-consistent solution for the oscillating dipole moments Pj may be found, as is discussed in Section 4 below; from these P the absorption and scattering cross sections are computed.6 If DDA solutions are obtained for two independent polarizations of the incident wave, then the complete amplitude-scattering matrix can be determined.6 With the recognition that the polarizabilities aj may be tensors, the DDA can readily be applied to anisotropic materials.6" 0 The extension of the DDA to treat materials with nonzero magnetic susceptibility is also straightforward, "" 2 although for most applications magnetic effects are negligible. B. Relation to Other Methods The DDA is one particular discretization method for solv- ing Maxwell's equations in the presence of a target. The DDA was anticipated in the engineering literature by the method of moments,'3 -'5 which is a general prescription for discretization of the integral form of Maxwell's equa- tions. There is a close correspondence between the DDA and discretizations that are based on the digitized Green's function (DGF) method'6 or the volume-integral equation formulation (VIEF).17-22 In the low-frequency limit kd - 0 (where d is the interdipole spacing or the side of a cubi- cal subvolume and k = /c), it has been shown" " 2 2' that the DDA is indeed equivalent to the DGF/VIEF. How- ever, researchers developing the DGF/VIEF approach encountered difficulties in trying to calculate O[(kd)2 ] corrections, whereas the DDA, because the fundamental approximation is simple and physically clear, permits rig- orous derivation of the O[(kd)2] correction (see Subsection 3.B below). On the other hand, the integral form of Maxwell's equations often serves as a convenient point of departure, and there is now a strong interplay between both approaches. The principal limitation of the DDA involves the han- dling of target boundaries, since the geometry of the DDA array has a minimum length scale that is equal to the interdipole spacing d. For targets with large values of the refractive index m, the accuracy of the DDA suffers. In principle we can make the DDA as accurate as we desire by increasing N, the number of dipoles representing the target (thus decreasing d), but very large values of N become computationally prohibitive. 3. SPECIFICATION OF THE DIPOLE ARRAY A. Geometry Each dipole may be thought of as representing the polar- izability of a particular subvolume of target material. If we wished to approximate a certain target geometry (e.g., a sphere) with a finite number of dipoles, then we might consider using some number of closely spaced, weaker dipoles in regions near the target boundaries to do a bet- ter job of approximating the boundary geometry. Some researchers have followed this approach23 24 ; although it clearly offers benefits, there is also a substantial cost: as we show below, FFT methods can greatly accelerate the computations that are required for solving the scattering problem, but only if the dipoles are located on a periodic lattice. With FFT methods we can perform calculations (on a workstation) for numbers N of dipoles (e.g., 105) that are much greater than the largest values (- 104) that could readily be handled without FFT methods. We therefore elect to use FFT methods, accepting the requirement that the dipoles be located on a periodic lattice. We further restrict ourselves to cubic lattices. This option still permits considerable latitude in rep- resentation of the continuum target. For a given target geometry we use the following algorithm to generate the dipole array. We assume that the target orientation is fixed relative to a coordinate system x, 9, i. 1. Generate a trial lattice that is defined by the lattice spacing d and coordinates (xo, yo, zo) of the lattice point nearest the origin. 2. Let the array now be defined as all lattice sites located within the volume V of the continuum target. 3. Try different values of d and (x0, yo, zo) and maxi- mize some goodness-of-fit criterion, subject to the con- straint that the number of points in the dipole array not be so large as to be computationally prohibitive. The goodness-of-fit criteria are basically geometric in nature, as well as somewhat arbitrary, 6 and are not discussed here. 4. We now have a list of occupied sites j = 1, ... , N. We imagine that each of these occupied sites represents a cubic subvolume d3 of material centered on the site. In the current implementation of DDSCAT we do not distin- guish between lattice sites near the surface and those in the interior. We therefore rescale the array by requiring that d = (V/N)" 3 , so that the volume Nd3 of the occu- pied lattice sites is equal to the volume V of the origi- nal target. 5. For each occupied site j assign a dipole polarizabil- ity aj. B. Dipole Polarizabilities There has been some controversy concerning the best method of assigning the dipole polarizabilities. Purcell and Pennypacker 5 used the Clausius-Mossotti polar- izabilities: CM= 3d3 ej-1aj 47r e + 2 (1) where e is the dielectric function of the target material at location rj. As is well known,2 5 for an infinite cubic lattice the Clausius-Mossotti prescription is exact in the dc limit kd - 0. Draine6 and Goedecke and O'Brien 6 showed that the polarizabilities should also include a radiative reaction correction of O[(kd)3]. Goedecke and O'Brien16 proposed a correction term of 0[(kd) 2]; the same correction term was independently proposed by Iskander et al.2 0 and by Hage and Greenberg.2 ' Unfortunately the derivation of this correction assumed that the electric B. T. Draine and P. J. Flatau

Vol. 11, No. 4/April 1994/J. Opt. Soc. Am. A 1493 field could be taken to be uniform over cubical regions of volume d 3, and this assumption itself introduces errors of O[(kd)2].26 Another approach was taken by Dungey and Bohren,2 7 who arrived at a different O[(kd)2 ] correc- tion by requiring that each point dipole have the same polarizability as a finite sphere of diameter d, but with a modified dielectric constant. To clarify this issue Draine and Goodman2 6 considered the following problem: for what polarizability a(Co) will an infinite lattice of polarizable points have the same dispersion relation as a continuum of refractive index m(co)? The lattice dispersion relation (LDR) may be found analytically; in the long-wavelength limit kd << 1, the polarizability a(co) is given as a series expansion in powers of kd and m2 = : a LDR surface monolayer. As a result, when ml >> 1 the DDA can seriously overestimate absorption cross sections, even in the dc limit Imlkd << 1. For materials with ml >> 1 it appears that other techniques (e.g., the method of Rouleau and Martin29 ) may be superior to the DDA. 4. COMPUTATIONAL CONSIDERATIONS A. Scattering Problem The electromagnetic scattering problem must be solved for the target array of point dipoles (j = 1, ... , N) with polarizabilities aj, located at positions r. Each dipole has a polarization Pj = aEj, where Ej is the electric field at r that is due to the incident wave En = Eo exp(ik r - it) plus the contribution of each of the aCM 1 + (aCM/d3)[(bi + m2b2 + m2b3S)(kd)2 - (2/3)i(kd)3]' (2) bi = -1.891531, b3 = -1.7700004, b2 = 0.1648469, 3 _ 1 (e)2j=1 other N - 1 dipoles: (3) Ej = Eincj - AjkPk, k~j where and are unit vectors defining the incident direction and the polarization state. The O[(kd)2] term in the LDR expansion differs from previously proposed corrections.16'2 0'21'27 The LDR prescription for a(w) is, by construction, optimal for wave propagation on an infi- nite lattice, and it is reasonable to assume that it will also be a good choice for finite dipole arrays. Exten- sive DDA calculations for spheres, comparing different prescriptions for the dipole polarizability, confirm that the LDR prescription appears to be best for mlkd < 1.26 The different prescriptions have also been compared in calculations of scattering by two touching spheres2 8 (see Fig. 9 below). C. Validity Criteria There are two obvious criteria for validity of the DDA: (1) lmlkd ' 1 (so that the lattice spacing d is small com- pared with the wavelength of a plane wave in the target material), and (2) d must be small enough (N must be large enough) to describe the target shape satisfactorily. Define the effective radius aeff of a target of volume V by aeff (3V/4v-)"3 . The first criterion is then equiva- lent to N > (47r/3)lml 3 (kaeff)3 . (4) Thus targets with large values of ml or scattering prob- lems with large values of kaeff will require that large numbers of dipoles be used to represent the targets. Unfortunately the second criterion has not yet been formulated precisely. Draine6 shows that even in the kd - 0 limit the polarizations are too large in the surface monolayer of dipoles in a pseudosphere, and similar errors must also occur for other target shapes. As a result the rate of energy absorption by the dipoles in the surface monolayer is too large, which leads to an error in the overall absorption cross section in proportion to the frac- tion -N- 3 of the total volume that is contributed by the where -AjkPk is the electric field at r that is due to dipole Pk at location rk, including retardation effects. Each element Ajk is a 3 X 3 matrix: = exp(ikrjk)A jk, = ~ 3 1 rjk X Fk2 (jk r jk r13 ) + - (31Pk jk -13) j k, (6) where k coc, rjk rj - rk, rjk (rj - rk)/rjk, and 13 is the 3 X 3 identity matrix. Defining Ajj = aj- re- duces the scattering problem to finding the polarizations Pj that satisfy a system of 3N complex linear equations: N Y AjkPk = Emcj. k=1 (7) Once Eq. (7) has been solved for the unknown polariza- tions Pj, the extinction and absorption cross sections Cext and Cabs may be evaluated6 : 4irk N Cext = Y Im(Ein * j)IEoJ2 _j (8) Cabs = 4Jl- 1 { I [P * (aj-7)*Pj*1 - 2 k3IPi2}. (9) The scattering cross section Csca = Cext - Cabs. Dif- ferential scattering cross sections may also be directly evaluated once the Pj are known.6 In the far field the (5) B. T. Draine and P. J. Flatau

1494 J. Opt. Soc. Am. A/Vol. 11, No. 4/April 1994 scattered electric field is given by Eca = exp(ikr) exp(-iki rj)(PP - 13 )Pj- (10) The problem is that matrix A is large and full, i.e., a matrix that in general has few zero elements. We would like (1) to solve the linear system of Eq. (7) efficiently and (2) to solve efficiently for multiple cases of the incident plane wave Enc. LU decomposition3 0 is the method of choice for small problems. It solves problem (2) because once the LU decomposition has been obtained the so- lution for a new right-hand side Einc requires just one matrix-vector multiplication. However, it is not feasible or efficient for large problems because of the need to store A and the LU decomposition and because of computer time that is proportional to N3 . For homogeneous, rectangular targets the matrix A has block-Toeplitz symmetry,3' and algorithms exist for finding A-' in O(N2 ) operations rather than in the O(N 3) operations that are required for general matrix inversion. It seems possible that the fundamental symmetries of A may permit efficient algorithms for finding A-' even for targets that are inhomogeneous or nonrectangular. B. Complex-Conjugate Gradient Method Rather than direct methods for solving Eq. (7), CCG methods for finding P iteratively have proven effective and efficient. The particular CCG algorithm that is used in DDSCAT is described by Draine.6 As P has 3N unknown elements, CCG methods in general are only guaranteed to converge in 3N iterations. In fact, how- ever, _ 10_102 iterations are often found to be sufficient to obtain a solution to high accuracy.63 2 The choice of conjugate-gradient variant may influence the conver- gence rate. 33 4 5 When N is large, CCG methods are much faster than are direct methods for finding P. The fact that CCG methods in practice converge relatively rapidly is presum- ably a consequence of the basic symmetries of A. As an iterative technique the CCG requires an ini- tial guess for P. The simplest choice is P = 0. At- tempts to improve on this by use of the scattering-order approximation3 6 as a first estimate for P were found to offer little, if any, advantage over starting with P = 0.6 The iterative method that is described here must be repeated for each incident plane wave [cf. Eq. (7)]. Notice that if we had A-' then for each new incident wave Einc the solution P = A-'Emc would be obtained from a single matrix-vector multiplication. In principle we can recover eigenvalues and eigenvectors from conjugate- gradient iteration.3 7 38 Some attempts to improve the efficiency of the conjugate-gradient method for multiple incident electromagnetic fields have been reported.34 39 ' FORTRAN implementations of approximately 15 CCG methods are available from us.42 C. Fast Fourier Transforms The computational burden in the CCG method primar- ily consists of matrix-vector multiplications of the form A v. Goodman et al.3 2 show that the structure of the matrix A implies that such multiplications are essen- tially convolutions, so that FFT methods can be employed to evaluate A v in O(N ln N) operations rather than in the 0(N 2) operations that are required for general matrix-vector multiplication. Since N is large this is an important calculational breakthrough. As mentioned above, FFT methods require that the dipoles be situated on a periodic lattice, which is most simply taken to be cubic. FFT methods in effect require the computation of three- dimensional FFT's over 8NXNYNZ points, where NXNyNz is the number of sites in a rectangular region of the lattice containing all the N occupied lattice sites. Therefore, in the case of a fractal structure with a large volume-filling factor for vacuum, FFT methods may lose some of their advantage over conventional techniques for evaluating A * v, since NX.NNZ may be much larger than the actual number of dipoles N. We note that CCG and FFT techniques are routinely used for numerical solutions of electromagnetic problems in engineering.3 44 3 ,4 4 D. Memory and CPU Requirements If the N dipoles were located at arbitrary positions rj, then the 9N2 elements of A would be nondegenerate. Storage of these elements (with 8 bytes/complex number) would require 72(N/103 )2 Mbytes. By locating the dipoles on a lattice, the elements of A become highly degenerate, since they depend only on the displacement rj - r,. As a result the memory requirements depend approximately linearly on N rather than on N 2 . The program DDSCAT has a total memory requirement of -0.58(N.NN./1000)'Mbytes, where N. X N, X N is the rectangular volume containing all the N dipoles. Thus for a target fitting into a 32 32 x 32 portion of the lattice (e.g., an N = 17,904 pseudosphere), DDSCAT requires -19 Mbytes of memory. We therefore see that memory requirements begin to be a consideration for targets that are much larger than -104 dipoles. The CPU requirements also are significant for large targets. On a Sun 4/50 (Sparcstation IPX), a single CCG iteration requires 3.0(N.NyN,/10 3 ) CPU s; thus one iteration for a N = 17,904 pseudosphere requires -100 CPU s. Between 10 and 102 iterations are typi- cally necessary to solve for a single incident direction and polarization; thus orientational averaging with the DDA can be time consuming if many target orientations are required. It is for this reason that T-matrix methods,4 5 46 which exploit efficient procedures for ori- entational averaging,4 7 are competitive for some target geometries as well as recursive T-matrix algorithms currently being developed by Chew and Lu4 8 and Chew et al. 4 9 5. ACCURACY OF THE DISCRETE-DIPOLE APPROXIMATION At this time there is no known way to predict precisely the accuracy of DDA calculations. Instead we rely on ex- amples to guide us. Spheres are most convenient since exact solutions are readily available for comparison. Note that the DDA does not give preferential treatment to any geometry (with the possible exception of rectangular targets), so spheres constitute a representative test of B. T. Draine and P. J. Flatau

Vol. 11, No. 4/April 1994/J. Opt. Soc. Am. A 1495 the accuracy of the DDA. Clusters of spheres, for which exact solutions can now also be computed, are also useful. Several earlier studies have examined the accuracy of the DDA,5,6,23,27 50,5' but only one26 employed the LDR polarizabilities. The DDA is most accurate for targets with refractive indices m near unity, with the accuracy declining (for fixed kaeff and N) as Im - 11 increases. Here we employ the LDR polarizabilities and concentrate on two cases: m = 1.33 + 0.01i, representative of moder- ately polarizable material, and m = 2 + i, representative of a strongly polarizable, strongly absorbing material in the visible or the ultraviolet. As is well known, very large refractive indices (Imi >> 1) occur for conducting materials in the infrared; the DDA is not particularly well suited for this regime, and we do not consider any such examples here. We will examine the dependence of the accuracy on both kaeff and N. A. Spheres 1. Total Cross Sections For a given incident wave the scattering problem is solved by CCG iteration until a specified level of accuracy is achieved. All the results shown here are converged until IA P - EincI < lO-3lE mc. Figures 1 and 2 show the errors in the absorption and scattering efficiency fac- tors Qabs 3 Cabs/7ra 2 and Qsca 3 Csca/ra2 , as functions of x = 2ra/A, for pseudospheres containing N dipoles; the curves are labeled by N. Each curve terminates at the point at which Imikd - 1. The pseudospheres are obviously not precisely spherically symmetric; the DDA calculations in Figs. 1-8 are for radiation incident along the (1, 1, 1) direction. Figure 1 shows results for m = 1.33 + 0.01i (e = 1.769 + 0.0266i). It is clear that for a modest refractive index, such as m = 1.33 + 001i, Qsca and Qabs may be obtained to 2% accuracy provided that N is large enough that the criterion Imlkd < 1 is satisfied. Figure 2 shows results for m = 2 + i (e = 3 + 4i). In this case the larger refractive index leads to larger errors, but the errors in Qsca are still quite small. Somewhat surprisingly the largest errors tend to occur in the Rayleigh limit ka ' 0.5; here Qsca is in error by - 3%, and Qabs is in error by 5% for N = 1064. For fixed ka the errors decline as N is increased. It is apparent that even for a refractive index m = 2 + i, which is moderately large and highly absorptive, scattering and absorption cross sections can be calculated to accuracies of - 1% and - 4%, respectively, for pseudospheres containing - 104 dipoles. 2. Differential Scattering Figures 1 and 2 show that the total scattering cross section can be calculated fairly accurately. Differential scattering efficiencies dQ/dfl = S, 1/r(ka) 2 for unpolar- ized light incident upon a sphere are shown in Figs. 3-8 for two refractive indices (m = 1.33 + 0.01i and m = 2 + i) and ka = 3, 5, and 7. The DDA results are shown for different numbers of dipoles N. The lower panel in Figs. 3-8 displays the fractional error: eo dQ/dfl(DDA) - dQ/dfl(Mie) (11) error = dQ/dfl (Mie) (1 The fractional error becomes quite large for scattering angle 0 > 1500 because dQ/dfl itself becomes quite small for backscattering; hence small absolute errors in dQ/dfl correspond to large fractional errors. From Figs. 1-8 we may draw the following conclusions: 1. For pseudospheres with Iml ' 3, total scattering and absorption can be computed to accuracies of a few Q) a(50. C.o 0.01 01 -3 CY~2 1 3 0 40 0 1 2 3 4 5 6 7 8 9 10 11 12 13 x=ka Fig. 1. Scattering and absorption for a sphere with refractive index m = 1.33 + 0.01i. The upper panel shows the exact values of the scattering and absorption efficiency factors Qsca and Qabs, obtained with Mie theory, as functions of ka, where k = 2X/A. The middle and lower panels show fractional errors in Qsca and Qabs, obtained from DDA pseudospheres, labeled by the number N of dipoles in the pseudosphere. a 5 0.1 a 0.01 IR 4 aI 0 0._ O V -2 1 6 a 4 .5 02 a) 0 L 0 1 2 3 4 5 6 7 x=ka ame as Fig. 1, but for the refractive a 9 10 B. T. Draine and P. J. Flatau Fig. 2. index m= 2 + i.

1496 J. Opt. Soc. Am. A/Vol. 11, No. 4/April 1994 B. T. Draine and P. J. Flatau a I-0.l 0.01 15 10 = 5 .5 -5 °-1 0 -15 -20 0 20 40 60 80 0(.)100 120 140 160 180 Fig. 3. Differential scattering cross section for pseudospheres with m = 1.33 + 0.01i and ka = 3. The curves are labeled by N, the number of dipoles in the pseudosphere. 10 0 C. 0.01 F N ._ '0 0) 20 15 10 5 0 -5 Mie theory N= 17,904 =7664 N=2320 ~~~~_ \ ~~~~~~~~~~Qsc./47r m= 1.33+0.Oli 27ra/X=5 rms, N=17,904: 1.8% 7664: 3.1% 2320: 3.4% ,jj,'~, I' I. 0 20 40 60 80 5(I)100 120 140 Fig. 4. Same as Fig. 3 but for ka = 5. 160 180 percent if the validity criterion [relation (4)] is satisfied, and N ' 103, so that the surface geometry is reasonably well approximated. 2. For pseudospheres, if the validity criterion [relation (4)] is satisfied, differential scattering cross sections can be calculated to reasonable accuracy. For m = 1.33 + 0.01i, fractional errors in dQsca/dl are within 20% for all scattering directions, with rms errors of G 6%. For m = 2 + i, fractional errors in dQ4Ca/dQ are within - 30% for all scattering directions, with rms values of S 11%. B. Two Contiguous Spheres Recent investigations have studied the general solution to Maxwell's equations for clusters of spheres.52 5 3 The availability of this exact solution offers a unique opportu- nity to test the accuracy of the DDA when it is applied to nonspherical targets. 28 54 Figure 9 presents a comparison of a two-sphere, multi- pole solution with the DDA as a function of x = kaeff- Results are shown for three different prescriptions for the dipole polarizabilities: CMRR (Clausius-Mossotti plus radiative reaction), DGF-VIEF, and LDR. The refrac- tive index is m = 1.33 + 0.01i, and the dipoles are placed on a 64 X 32 X 32 cubic lattice. The calculations are for 10 c: 1'0 I C. 'X 0.1 0.01 15 ig 10 C' 0 5 ._ 0 - 5 0 20 40 60 Fig. 5. Same C-0 ' 0.1 '0 0.01 cr \ 5 .531 0 '--l -20 0 20 40 80 (0)l0° as Fig. 3 but 60 80 0(.100 120 Fig. 6. Same as Fig. 3 but for m = 2 + i. Mie theory N=59,728 N= 17,904 N=7664 -- N=2320 Q.../4Tr \ m=2+i - 2nra/X=3 rmNs N=59,728: 4.5% 17,904: 6.6% 7664: 8.1% i 2320: 10.7% ~~ \ / 120 140 160 180 for ka = 7. 140 160 180 ..........................

B. T. Draine and P. J. Flatau Vol. 11, No. 4/April 1994/J. Opt. Soc. Am. A 1497 1 V I 0.1 0.01 5 3U o-10 -15 _s Mie theory N=59728 N=17,904 N=7664 -m=2+i\ 2Ta/ =5 /\/ 17,904: 7.9% 7664: 8.3% 0 0 0 0 0 100 |120 140,160,180,,,1, 0 20 40 60 80 or()100 120 140 160 180 Fig. 7. Same as Fig. 3, but for m = 2 + i and a = 5. 10 1c a. 0.1 0.01 5 '0c -5 C.- 1 0'0 -1 h-20 q) -25 -30 0 20 40 60 80 0(.)100 120 140 160 180 Fig. 8. Same as Fig. 3, but for m = 2 + i and ka = 7. 6. RECENT RESEARCH The DDA has recently been applied to a number of diverse problems, including (since 1990) calculations of scattering 0.1 0.01 0.001 2 c 1 I c0 3- -1 0 4) -2 6 " 4 o 2 0 -2 2 6 a 10 2 4 6 8 10 x=ka. 1 4 Fig. 9. (a) Qsca CCa/raeff and Qabs a Cabs/7aeff (calculated by modal analysis) for two contiguous spheres of refractive index m = 1.33 + 0.01i, for radiation incident at an angle of 30° relative to the axis passing through the center of the two spheres. (b) Fractional error in computed value of Qsca and (c) Qabs for a DDA representation containing N = 2 X 17256 dipoles. Results are shown for three different prescriptions for the dipole polarizabilities: LDR, DGF/VIEF, and CMRR as functions of x = kaeff (source, Ref. 28). 0 20 40 60 80 100 120 140 10 102 = 10' )2 0.1 10-2 103 102 unpolarized light incident at an angle a = 30° relative to the axis passing through the centers of the two spheres. As is noted in Section 3, the LDR method for specifying the polarizabilities results in very accurate values for Qsca and Qabs and is clearly preferable to the CMRR and DGF/VIEF prescriptions. In Fig. 10 we compare the exact solution for scattering by two touching spheres (solid curve) with the solution derived from the DDA (X's) for m = 1.33 + 0.001i and kaeff = 10. The exact results for S11 and S22 are com- pared with results that were computed with the DDA. The overall agreement is very good, although errors can be seen for scattering angle 0 ' 500. S3 2)2 10l 0.1 10-2 160 180 0 20 40 60 80 100 120 140 160 180 O() Fig. 10. Scattering by two contiguous spheres with refractive index m = 1.33 + 0.001i for x = kaeff = 10. Radiation is in- cident along the axis passing through the centers of the two spheres. Results are shown for two elements (S1u and S 22 ) of the 4 x 4 scattering matrix.55 S and S22 computed with the DDA are compared with exact results (source, Ref. 28). -,X,

1498 J. Opt. Soc. Am. A/Vol. 11, No. 4/April 1994 and absorption by rough and porous particles,5 6 inter- stellar graphite particles,5 75 8 aggregate particles,5 0'5 9 -6 5 inhomogeneous particles, 6 6 structures on surfaces, 6 7 mi- crowave scattering by ice crystals,68 and single scattering properties of ice crystals. 6 9 7. SUMMARY The discrete-dipole approximation, or the coupled-dipole method, is a flexible and powerful tool for computing scat- tering and absorption by arbitrary targets. For targets with Iml S 2 (with m being the complex refractive index), scattering and absorption cross sections can be evaluated to accuracies of a few percent provided that the number N of dipoles satisfies relation (4). The use of CCG and FFT algorithms permits calculations for N as large as - 105, so that scattering problems with kaff ' 10 can be studied with scientific workstations. ACKNOWLEDGMENTS B. T. Draine's research has been supported in part by National Science Foundation grant AST-9017082. P. J. Flatau has been supported in part by Atmo- spheric Radiation Measurement grant DOE DE-FG03- 91ER61198 and Western Regional Center for Global Environmental Change grant UCSD93-1220. REFERENCES AND NOTES 1. C. F. Bohren and S. B. Singham, "Backscattering by non- spherical particles: a review of methods antsuggested new approaches," J. Geophys. Res. 96, 5269-5277 (1991). 2. The FORTRAN program DDSCAT.4b is available from the authors. Direct queries to the Internet address draine astro.princeton.edu or pflatau@macao.ucsd.edu. 3. H. DeVoe, "Optical properties of molecular aggregates. I. Classical model of electronic absorption and refraction," J. Chem. Phys. 41, 393-400 (1964). 4. H. DeVoe, "Optical properties of molecular aggregates. II. Classical theory of the refraction, absorption, and optical activity of solutions and crystals," J. Chem. Phys. 43, 3199-3208 (1965). 5. E. M. Purcell and C. R. Pennypacker, "Scattering and ab- sorption of light by nonspherical dielectric grains," Astro- phys. J. 186, 705-714 (1973). 6. B. T. Draine, "The discrete-dipole approximation and its application to interstellar graphite grains," Astrophys. J. 333, 848-872 (1988). 7. S. B. Singham and G. C. Salzman, "Evaluation of the scatter- ing matrix of an arbitrary particle using the coupled dipole approximation," J. Chem. Phys. 84, 2658-2667 (1986). 8. S. B. Singham and C. F. Bohren, "Light scattering by an arbitrary particle: a physical reformulation of the coupled- dipoles method," Opt. Lett. 12, 10-12 (1987). 9. H. A. Lorentz, Theory of Electrons (Teubner, Leipzig, 1909). 10. E. L. Wright, "The ultraviolet extinction from interstellar graphitic onions," Nature (London) 366, 227-228 (1988). 11. A. Lakhtakia, "General theory of the Purcell-Pennypacker scattering approach and its extension to bianisotropic scat- terers," Astrophys. J. 394, 494-499 (1992). 12. A. Lakhtakia, "Strong and weak forms of the method of moments and the coupled dipole method for scattering of time-harmonic electromagnetic fields," Int. J. Mod. Phys. C3, 583-603 (1992). 13. J. H. Richmond, "Scattering by a dielectric cylinder of arbi- trary cross-section shape," IEEE Trans. Antennas Propag. AP-13, 334-343 (1965). 14. R. F. Harrington, Field Computation by Moment Methods (Macmillan, New York, 1968). 15. R. Harrington, "Origin and development of the method of moments for field computation," IEEE Antennas Propag. Mag. 32(3), 31-35 (1990). 16. G. H. Goedecke and S. G. O'Brien, "Scattering by irregular inhomogeneous particles via the digitized Green's function algorithm," Appl. Opt. 27, 2431-2438 (1988). 17. D. E. Livesay and K. Chen, "Electromagnetic fields induced inside arbitrarily shaped biological bodies," IEEE Trans. Microwave Theory Tech. MTT-22, 1273-1280 (1974). 18. A. W. Glison, "Recent advances in frequency domain techniques for electromagnetic scattering problems," IEEE Trans. Antennas Propag. 25, 2867-2871 (1989). 19. E. K. Miller, "A selective survey of computational elec- tromagnetics," IEEE Trans. Antennas Propag. 36, 1281- 1305 (1988). 20. M. F. Iskander, H. Y. Chen, and J. E. Penner, "Optical scattering and absorption by branched chains of aerosols," Appl. Opt. 28, 3083-3091 (1989). 21. J. I. Hage and J. M. Greenberg, "A model for the opti- cal properties of porous grains," Astrophys. J. 361, 251- 259 (1990). 22. E. H. Newman and K. Kingsley, "An introduction to the method of moments," J. Comput. Phys. 68, 1-18 (1991). 23. C. Bourrely, P. Chiappetta, T. Lemaire, and B. Torr6sani, "Multidipole formulation of the coupled dipole method for electromagnetic scattering by an arbitrary particle," J. Opt. Soc. Am. A 9, 1336-1340 (1992). 24. J. M. Perrin and J. P. Sivan, "Light scattering by dust grains: effects of the state of the surface on the validity of the discrete dipole approximation," C. R. Acad. Sci. Paris Ser. II 316, 47-53 (1993). 25. J. D. Jackson, Classical Electromagnetism (Wiley, New York, 1975). 26. B. T. Draine and J. Goodman, "Beyond Clausius-Mossotti: wave propagation on a polarizable point lattice and the discrete dipole approximation," Astrophys. J. 405, 685- 697 (1993). 27. C. E. Dungey and C. F. Bohren, "Light scattering by nonspherical particles: a refinement to the coupled-dipole method," J. Opt. Soc. Am. A 8, 81-87 (1991). 28. P. J. Flatau, K. A. Fuller, and D. W. Mackowski, "Scattering by two spheres in contact: comparisons between discrete dipole approximation and modal analysis," Appl. Opt. 32, 3302-3305 (1993). 29. F. Rouleau and P. G. Martin, "A new method to calculate the extinction properties of irregularly haped particles," Astrophys. J. 414, 803-814 (1993). 30. W. Hager, Applied Numerical Linear Algebra (Prentice-Hall, Englewood Cliffs, N.J., 1988). 31. P. J. Flatau, G. L. Stephens, and B. T. Draine, "Light scattering by rectangular solids in the discrete-dipole ap- proximation: a new algorithm exploiting the block-Toeplitz structure," J. Opt. Soc. Am. A 7, 593-600 (1990). 32. J. J. Goodman, B. T. Draine, and P. J. Flatau, "Application of fast-Fourier-transform techniques to the discrete-dipole approximation," Opt. Lett. 16, 1198-1200 (1991). 33. T. K. Sarkar, X. Yang, and E. Arvas, "A limited survey of various conjugate gradient methods for complex matrix equations arising in electromagnetic wave interactions," Wave Motion 10, 527-546 (1988). 34. A. F. Peterson, S. L. Ray, C. H. Chan, and R. Mittra, "Nu- merical implementation of the conjugate gradient method and the CG-FFT for electromagnetic scattering," in Applica- tion of Conjugate Gradient Method to Electromagnetics and Signal Processing, T. K. Sarkar, ed. (Elsevier, New York, 1991), Chap. 5. 35. R. W. Freund and N. M. Nachtigal, "QMR: a quasi-minimal residual method for non-Hermitian linear systems," Numer. Math. 60, 315-339 (1991). 36. S. B. Singham and C. F. Bohren, "Light scattering by an arbitrary particle: the scattering order formulation of the coupled-dipole method," J. Opt. Soc. Am. A 5, 1867- 1872 (1988). 37. L. Knockaert, "A note on the relationship between the con- jugate gradient method and polynomials orthogonal over the spectrum of a linear operator," IEEE Trans. Antennas Propag. AP-35, 1089-1091 (1987). B. T. Draine and P. J. Flatau

Vol. 11, No. 4/April 1994/J. Opt. Soc. Am. A 1499 38. A. F. Peterson, C. F. Smith, and R. Mittra, "Eigenvalues of the moment-method matrix and their effect on the con- vergence of the conjugate gradient algorithm," IEEE Trans. Antennas Propag. 36, 1177-1179 (1988). 39. C. F. Smith, A. F. Peterson, and R. Mittra, "A conjugate gradient algorithm for the treatment of multiple incident electromagnetic fields," IEEE Trans. Antennas Propag. 37, 1490-1493 (1989). 40. P. Joly, "Resolution de systemes lineaires avec plusieurs members par la methode du gradient conjugu6," Tech. Rep. R-91012 (Publications du Laboratoire d'Analyse Numerique, Universit6 Pierre et Marie Curie, Paris, 1991). 41. V. Simoncini and E. Gallopoulos, "An iterative method for nonsymmetric systems with multiple right-hand sides," Tech. Rep. 1242 (Center for Supercomputing Research and Development, University of Illinois at Urbana-Champaign, Champaign, Ill., 1992). 42. P. J. Flatau, T. Schneider, and F. Evans, CCG- PAK-FORTRAN. Conjugate gradient package for solv- ing complex matrix equations (1993). Available from pflatau~ucsd.edu. 43. D. T. Borup and 0. P. Gandhi, "Calculation of high- resolution SAR distributions in biological bodies using the FFT algorithm and conjugate gradient method," IEEE Trans. Microwave Theory Tech. MTT-33, 417-419 (1985). 44. T. K. Sarkar, E. Arvas, and S. M. Rao, "Application of the fast Fourier transform and the conjugate gradient method for efficient solution of electromagnetic scattering from both electrically large and small conducting bodies," Electromag- netics 5, 99-122 (1985). 45. P. Barber and C. Yeh, "Scattering of electromagnetic waves by arbitrarily shaped dielectric bodies," Appl. Opt. 14, 2864-2872 (1975). 46. D.-S. Wang and P. W. Barber, "Scattering by inhomogeneous nonspherical objects," Appl. Opt. 18, 1190-1197 (1979). 47. M. I. Mishchenko, "Light scattering by randomly oriented axially symmetric particles," J. Opt. Soc. Am. A 8, 871- 882 (1991). 48. W. C. Chew and C.-C. Lu, "NEPAL-an algorithm for solving the volume integral equation," Microwave Opt. Tech. Lett. 6, 185-188 (1993). 49. W. C. Chew, Y. M. Wang, and L. Gurel, "Recursive algorithm for wave-scattering using windowed addition theorem," J. Electromagn. Waves Appl. 6, 1537-1560 (1992). 50. J. C. Ku and K.-H. Shim, "A comparison of solutions for light scattering and absorption by agglomerated or arbitrarily- ,shaped particles," J. Quant. Spectrosc. Radiat. Transfer 47, 201-220 (1992). 51. J. C. Ku, "Comparisons of coupled-dipole solutions and dipole refractive indices for light scattering and absorption by arbitrarily shaped or agglomerated particles," J. Opt. Soc. Am. A 10, 336-342 (1993). 52. K. A. Fuller, "Optical resonances and two-sphere systems," Appl. Opt. 30, 4716-4731 (1991). 53. D. W. Mackowski, "Analysis of radiative scattering for mul- tiple sphere configurations," Proc. R. Soc. London Ser. A 433, 599-614 (1991). 54. G. W. Kattawar and T. J. Humphreys, "Electromagnetic scattering from two identical pseudospheres," in Light Scat- tering by Irregularly Shaped Particles, D. W. Schuerman, ed. (Plenum, New York, 1980), pp. 177-190. 55. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, New York, 1983). 56. J. M. Perrin and J. P. Sivan, "Scattering and polarisation of light by rough and porous interstellar grains," Astron. Astrophys. 247, 497-504 (1991). 57. J. M. Perrin and J. P. Sivan, "Porosity and impurities within interstellar grains. Is the ultraviolet bump still explained by carbonaceous material?" Astron. Astrophys. 228, 238-245 (1990). 58. B. T. Draine and S. Malhotra, "On graphite and the 2175 A extinction profile," Astrophys. J. 414, 632-645 (1993). 59. M. F. Iskander, H. Y. Chen, and J. E. Penner, "Resonance op- tical absorption by fractal agglomerates of smoke aerosols," Atmos. Environ. 25A, 2563-2569 (1991). 60. R. A. West, "Optical properties of aggregate particles whose outer diameter is comparable to the wavelength," Appl. Opt. 30, 5316-5324 (1991). 61. R. A. West and P. H. Smith, "Evidence for aggregate par- ticles in the atmospheres of Titan and Jupiter," Icarus 90, 330-333 (1991). 62. T. Kozasa, J. Blum, and T. Mukai, "Optical properties of dust aggregates. I. Wavelength dependence," Astron. Astrophys. 263, 423-432 (1992). 63. T. Mukai, H. Ishimoto, T. Kozasa, J. Blum, and J. M. Green- berg, "Radiation pressure forces of fluffy porous grains," Astron. Astrophys. 262, 315-320 (1992). 64. M. J. Wolff, G. C. Clayton, P. G. Martin, and R. E. Schulte- Ladbeck, "Modeling composite and fluffy grains: the effects of porosity," Astrophys. J. (to be published). 65. S. B. Singham and C. F. Bohren, "Scattering of unpolarized and polarized light by particle aggregates of different size and fractal dimension," Langmuir 9, 1431-1435 (1993). 66. J. M. Perrin and P. L. Lamy, "On the validity of effective- medium theories in the case of light extinction by inhomo- geneous dust particles," Astrophys. J. 364, 146-151 (1990). 67. M. A. Taubenblatt and T. K. Tran, "Calculation of light scattering from particles and structures on a surface by the coupled-dipole method," J. Opt. Soc. Am. A 10, 912- 919 (1993). 68. K. F. Evans and J. Vivekanandan, "Multiparameter radar and microwave radiative transfer modeling of nonspheri- cal atmospheric ice particles," IEEE Trans. Geosci. Remote Sensing 28, 423-437 (1990). 69. P. J. Flatau, "Scattering by irregular particles in anomalous diffraction and discrete dipole approximations," Atmos. Sci. Paper 517 (Department of Atmospheric Science, Colorado State University, Fort Collins, Colo., 1992). B. T. Draine and P. J. Flatau