## Abstract

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 (|*m*| ≲ 2), the DDA permits calculations of scattering and absorption that are accurate to within a few percent.

© 1994 Optical Society of America

## 1. INTRODUCTION

The discrete-dipole approximation (DDA) is a flexible and powerful technique for computing scattering and absorption by targets of arbitrary geometry. The development of efficient algorithms and the availability of inexpensive computing power together have made the DDA the method of choice for many scattering problems. Bohren and Singham[1] reviewed recent advances in calculations of backscattering by nonspherical particles, including DDA results. In this paper we review the DDA, with particular 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 represent the target. In Section 3 we discuss these choices, with attention to the important question of dipole polarizabilities. 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 number 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 constants that are not too large (|*m*| ≲ 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 geometries 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 properties of molecular aggregates; retardation effects were not included, so DeVoe’s treatment was limited to aggregates 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, including radiative reaction corrections, is summarized by Draine.[6]

Nature provides the physical inspiration for the DDA: in 1909 Lorentz showed[9] that the dielectric properties of a substance could be directly related to the polarizabilities 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 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 **r*** _{j}* of the dipoles

*j*= 1, …,

*N*) and the dipole polarizabilities

*α*.

_{j}For monochromatic incident waves the self-consistent solution for the oscillating dipole moments **P*** _{j}* 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]*

_{j}With the recognition that the polarizabilities *α _{j}* may be tensors, the DDA can readily be applied to anisotropic materials.[6],[10] The extension of the DDA to treat materials with nonzero magnetic susceptibility is also straightforward,[11],[12] although for most applications magnetic effects are negligible.

#### B. Relation to Other Methods

The DDA is one particular discretization method for solving Maxwell’s equations in the presence of a target. The DDA was anticipated in the engineering literature by the method of moments,[13]–[15] which is a general prescription for discretization of the integral form of Maxwell’s equations. There is a close correspondence between the DDA and discretizations that are based on the digitized Green’s function (DGF) method[16] 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 cubical subvolume and *k* = *ω/c*), it has been shown[11],[12],[21] that the DDA is indeed equivalent to the DGF/VIEF. However, 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 rigorous 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 handling 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 polarizability 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 better job of approximating the boundary geometry. Some researchers have followed this approach[23],[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., ∼ 10^{5}) that are much greater than the largest values (∼ 10^{4}) 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 representation 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 $\widehat{x},\u0177,\u1e91,$.

- 1. Generate a trial lattice that is defined by the lattice spacing
*d*and coordinates (*x*_{0},*y*_{0},*z*_{0}) 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 (*x*_{0},*y*_{0},*z*_{0}) and maximize some goodness-of-fit criterion, subject to the constraint 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*d*^{3}of material centered on the site. In the current implementation of ddscat we do not distinguish between lattice sites near the surface and those in the interior. We therefore rescale the array by requiring that*d*= (*V/N*)^{1/3}so that the volume*Nd*^{3}of the occupied lattice sites is equal to the volume*V*of the original target. - 5. For each occupied site
*j*assign a dipole polarizability*α*._{j}

#### 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 polarizabilities:

where*∊*is the dielectric function of the target material at location

_{j}**r**

*. As is well known,[25] for an infinite cubic lattice the Clausius–Mossotti prescription is exact in the dc limit*

_{j}*kd*→ 0. Draine[6] and Goedecke and O’Brien[16] showed that the polarizabilities should also include a radiative reaction correction of

*O*[(

*kd*)

^{3}]. Goedecke and O’Brien[16] proposed a correction term of

*O*[(

*kd*)

^{2}]; the same correction term was independently proposed by Iskander

*et al.*[20] and by Hage and Greenberg.[21] Unfortunately the derivation of this correction assumed that the electric 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,[27] who arrived at a different

*O*[(

*kd*)

^{2}] correction 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 Goodman[26] considered the following problem: for what polarizability *α*(*ω*) will an infinite lattice of polarizable points have the same dispersion relation as a continuum of refractive index *m*(*ω*)? The lattice dispersion relation (LDR) may be found analytically; in the long-wavelength limit *kd* ≪ 1, the polarizability *α*(*ω*) is given as a series expansion in powers of *kd* and *m*^{2} = *∊*:

*â*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],[20],[21],[27] The LDR prescription for

*α*(

*ω*) is, by construction, optimal for wave propagation on an infinite lattice, and it is reasonable to assume that it will also be a good choice for finite dipole arrays. Extensive DDA calculations for spheres, comparing different prescriptions for the dipole polarizability, confirm that the LDR prescription appears to be best for |

*m*|

*kd*< 1.[26] The different prescriptions have also been compared in calculations of scattering by two touching spheres[28] (see Fig. 9 below).

#### C. Validity Criteria

There are two obvious criteria for validity of the DDA: (1) |*m*|*kd* ≲ 1 (so that the lattice spacing *d* is small compared 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 *a*_{eff} of a target of volume *V* by *a*_{eff} ≡ (3*V*/4*π*)^{1/3}. The first criterion is then equivalent to

*m*| or scattering problems with large values of

*ka*

_{eff}will require that large numbers of dipoles be used to represent the targets.

Unfortunately the second criterion has not yet been formulated precisely. Draine[6] 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 fraction ∼ *N*^{−1/3} of the total volume that is contributed by the surface monolayer. As a result, when |*m*| ≫ 1 the DDA can seriously overestimate absorption cross sections, even in the dc limit |*m*|*kd* ≪ 1. For materials with |*m*| ≫1 it appears that other techniques (e.g., the method of Rouleau and Martin[29]) 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 *α _{j}*, located at positions

**r**

*. Each dipole has a polarization*

_{j}**P**

*=*

_{j}*α*

_{j}**E**

*, where*

_{j}**E**

*is the electric field at*

_{j}**r**

*that is due to the incident wave*

_{j}**E**

_{inc,}

*=*

_{j}**E**

_{0}exp(

*i*

**k**·

**r**

*−*

_{j}*iωt*) plus the contribution of each of the other

*N*− 1 dipoles:

**A**

_{jk}**P**

*is the electric field at*

_{k}**r**

*that is due to dipole*

_{j}**P**

*. at location*

_{k}**r**

*, including retardation effects. Each element*

_{k}**A**

*is a 3 × 3 matrix:*

_{jk}**1**

_{3}is the 3 × 3 identity matrix. Defining

**A**

*≡*

_{jj}*α*

_{j}^{−1}reduces the scattering problem to finding the polarizations

**P**

*that satisfy a system of 3*

_{j}*N*complex linear equations:

**P**

*, the extinction and absorption cross sections*

_{j}*C*

_{ext}and

*C*

_{abs}may be evaluated[6]:

*C*

_{sca}=

*C*

_{ext}−

*C*

_{abs}. Differential scattering cross sections may also be directly evaluated once the

**P**

*are known.[6] In the far field the scattered electric field is given by*

_{j}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 **E**_{inc} LU decomposition[30] is the method of choice for small problems. It solves problem (2) because once the LU decomposition has been obtained the solution for a new right-hand side **E**_{inc} 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 *N*^{3}.

For homogeneous, rectangular targets the matrix **A** has block-Toeplitz symmetry,[31] and algorithms exist for finding **A**^{−1} in *O*(*N*^{2}) 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**^{−1} 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 3*N* unknown elements, CCG methods in general are only guaranteed to converge in 3*N* iterations. In fact, however, ∼10 − 10^{2} iterations are often found to be sufficient to obtain a solution to high accuracy.[6],[32] The choice of conjugate-gradient variant may influence the convergence rate.[33]–[35]

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 presumably a consequence of the basic symmetries of **A**.

As an iterative technique the CCG requires an initial guess for **P**. The simplest choice is **P** = 0. Attempts to improve on this by use of the scattering-order approximation[36] 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**^{−1} then for each new incident wave **E**_{inc} the solution **P** = **A**^{−1}**E**_{inc} would be obtained from a single matrix–vector multiplication. In principle we can recover eigenvalues and eigenvectors from conjugate- gradient iteration.[37],[38] Some attempts to improve the efficiency of the conjugate-gradient method for multiple incident electromagnetic fields have been reported.[34],[39]–[41] fortran implementations of approximately 15 CCG methods are available from us.[42]

#### C. Fast Fourier Transforms

The computational burden in the CCG method primarily consists of matrix–vector multiplications of the form **A** · **v**. Goodman *et al.*[32] show that the structure of the matrix **A** implies that such multiplications are essentially convolutions, so that FFT methods can be employed to evaluate **A** · **v** in *O* (*N* In *N*) operations rather than in the *O*(*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 8*N _{x}N_{y}N_{z}* points, where

*N*is the number of sites in a rectangular region of the lattice containing all the

_{x}N_{y}N_{z}*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

*N*may be much larger than the actual number of dipoles

_{x}N_{y}N_{z}*N*.

We note that CCG and FFT techniques are routinely used for numerical solutions of electromagnetic problems in engineering.[34],[43],[44]

#### D. Memory and CPU Requirements

If the *N* dipoles were located at arbitrary positions **r*** _{j}*, then the 9

*N*

^{2}elements of

**A**would be nondegenerate. Storage of these elements (with 8 bytes/complex number) would require 72(

*N*/10

^{3})

^{2}Mbytes. By locating the dipoles on a lattice, the elements of

**A**become highly degenerate, since they depend only on the displacement

**r**

*−*

_{j}**r**

*. As a result the memory requirements depend approximately linearly on*

_{i}*N*rather than on

*N*

^{2}. The program ddscat has a total memory requirement of ∼0.58(

*N*/1000) Mbytes, where

_{x}N_{y}N_{z}*N*×

_{x}*N*×

_{y}*N*is the rectangular volume containing all the

_{z}*N*dipoles. Thus for a target fitting into a 32 × 32 × 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 ∼10

^{4}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 _{x}N_{y}N_{z}*/10

^{3}) CPU s; thus one iteration for a

*N*= 17,904 pseudosphere requires ∼100 CPU s. Between 10 and 10

^{2}iterations are typically 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,[45],[46] which exploit efficient procedures for orientational averaging,[47] are competitive for some target geometries as well as recursive

*T*-matrix algorithms currently being developed by Chew and Lu[48] and Chew

*et al.*[49]

## 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 examples 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 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],[51] but only one[26] employed the LDR polarizabilities. The DDA is most accurate for targets with refractive indices *m* near unity, with the accuracy declining (for fixed *ka*_{eff} and *N*) as |*m* − 1| increases. Here we employ the LDR polarizabilities and concentrate on two cases: *m* = 1.33 + 0.01*i*, representative of moderately 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 |*m*| ≫ 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 *ka*_{eff} 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 |**A** · **P** − **E**_{inc} < 10^{−3}|**E**_{inc}|. Figures 1 and 2 show the errors in the absorption and scattering efficiency factors *Q*_{abs} ≡ *C*_{abs}/*πa*^{2} and *Q*_{sca} ≡ *C*_{sca}/*πa*^{2}, as functions of *x* = *πa*^{2}/λ for pseudospheres containing *N* dipoles; the curves are labeled by *N*. Each curve terminates at the point at which |*m*|*kd* ≈ 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.01*i* (*∊*) = 1.769 + 0.0266*i*). It is clear that for a modest refractive index, such as *m*= 1.33 + 0.01*i, Q*_{sca} and *Q*_{abs} may be obtained to ∼ 2% accuracy provided that *N* is large enough that the criterion |*m*|*kd* < 1 is satisfied.

Figure 2 shows results for *m* = 2 + *i*(*∊* = 3 + 4*i*). In this case the larger refractive index leads to larger errors, but the errors in *Q*_{sca} are still quite small. Somewhat surprisingly the largest errors tend to occur in the Rayleigh limit *ka* ≲ 0.5; here *Q*_{sca} is in error by ∼ 3%, and *Q*_{abs} 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 ≳10^{4} dipoles.

### 2. Differential Scattering

Figures 1 and 2 show that the total scattering cross section can be calculated fairly accurately. Differential scattering efficiencies d*Q*/dΩ = *S*_{11}/*π*(*ka*)^{2} for unpolarized light incident upon a sphere are shown in Figs. 3–8 for two refractive indices (*m* = 1.33 + 0.01*i* 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:

*θ*> 150° because d

*Q*/dΩ itself becomes quite small for backscattering; hence small absolute errors in d

*Q*/dΩ correspond to large fractional errors.

From Figs. 1–8 we may draw the following conclusions:

- 1. For pseudospheres with |
*m*| ≲ 3, total scattering and absorption can be computed to accuracies of a few percent if the validity criterion [relation (4)] is satisfied, and*N*≳ 10^{3}, 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.01*i*, fractional errors in d*Q*_{sca}/dΩ are within ∼20% for all scattering directions, with rms errors of ≲ 6%. For *m* = 2 + *i*, fractional errors in d*Q*_{sca}/dΩ are within ∼ 30% for all scattering directions, with rms values of ≲ 11%.

#### B. Two Contiguous Spheres

Recent investigations have studied the general solution to Maxwell’s equations for clusters of spheres.[52],[53] The availability of this exact solution offers a unique opportunity 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, multipole solution with the DDA as a function of *x* = *ka*_{eff}. Results are shown for three different prescriptions for the dipole polarizabilities: CMRR (Clausius–Mossotti plus radiative reaction), DGF–VIEF, and LDR. The refractive index is *m* = 1.33 + 0.01*i*, and the dipoles are placed on a 64 × 32 × 32 cubic lattice. The calculations are for unpolarized light incident at an angle *α* = 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 *Q*_{sca} and *Q*_{abs} 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 (×’s) for *m* = 1.33 + 0.001*i* and *ka*_{eff} = 10. The exact results for *S*_{11} and *S*_{22} are compared with results that were computed with the DDA. The overall agreement is very good, although errors can be seen for scattering angle θ ≳ 50°.

## 6. RECENT RESEARCH

The DDA has recently been applied to a number of diverse problems, including (since 1990) calculations of scattering and absorption by rough and porous particles,[56] interstellar graphite particles,[57],[58] aggregate particles,[50],[59]–[65] inhomogeneous particles,[66] structures on surfaces,[67] microwave scattering by ice crystals,[68] and single scattering properties of ice crystals.[69]

## 7. SUMMARY

The discrete-dipole approximation, or the coupled-dipole method, is a flexible and powerful tool for computing scattering and absorption by arbitrary targets. For targets with |*m*| ≲ 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 ∼ 10^{5}, so that scattering problems with *ka*_{eff} ≲ 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 Atmospheric Radiation Measurement grant DOE DE-FG03- 91ER61198 and Western Regional Center for Global Environmental Change grant UCSD93-1220.

## Figures

## REFERENCES AND NOTES

**1. **C. F. Bohren and S. B. Singham, “Backscattering by non- spherical particles: a review of methods and suggested new approaches,” J. Geophys. Res. **96**, 5269–5277 (1991) [CrossRef] .

**2. **The fortranprogram 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) [CrossRef] .

**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) [CrossRef] .

**5. **E. M. Purcell and C. R. Pennypacker, “Scattering and absorption of light by nonspherical dielectric grains,” Astrophys. J. **186**, 705–714 (1973) [CrossRef] .

**6. **B. T. Draine, “The discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J. **333**, 848–872 (1988) [CrossRef] .

**7. **S. B. Singham and G. C. Salzman, “Evaluation of the scattering matrix of an arbitrary particle using the coupled dipole approximation,” J. Chem. Phys. **84**, 2658–2667 (1986) [CrossRef] .

**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) [CrossRef] [PubMed] .

**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) [CrossRef] .

**11. **A. Lakhtakia, “General theory of the Purcell-Pennypacker scattering approach and its extension to bianisotropic scatterers,” Astrophys. J. **394**, 494–499 (1992) [CrossRef] .

**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 arbitrary cross-section shape,” IEEE Trans. Antennas Propag. **AP-13**, 334–343 (1965) [CrossRef] .

**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) [CrossRef] .

**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) [CrossRef] [PubMed] .

**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) [CrossRef] .

**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 electromagnetics,” IEEE Trans. Antennas Propag. **36**, 1281–1305 (1988) [CrossRef] .

**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) [CrossRef] [PubMed] .

**21. **J. I. Hage and J. M. Greenberg, “A model for the optical properties of porous grains,” Astrophys. J. **361**, 251–259 (1990) [CrossRef] .

**22. **E. H. Newman and K. Kingsley, “An introduction to the method of moments,” J. Comput. Phys. **68**, 1–18 (1991) [CrossRef] .

**23. **C. Bourrely, P. Chiappetta, T. Lemaire, and B. Torrésani, “Multidipole formulation of the coupled dipole method for electromagnetic scattering by an arbitrary particle,” J. Opt. Soc. Am. A **9**, 1336–1340 (1992) [CrossRef] .

**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) [CrossRef] .

**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) [CrossRef] .

**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) [CrossRef] [PubMed] .

**29. **F. Rouleau and P. G. Martin, “A new method to calculate the extinction properties of irregularly shaped particles,” Astrophys. J. **414**, 803–814 (1993) [CrossRef] .

**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 approximation: a new algorithm exploiting the block-Toeplitz structure,” J. Opt. Soc. Am. A **7**, 593–600 (1990) [CrossRef] .

**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) [CrossRef] [PubMed] .

**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) [CrossRef] .

**34. **A. F. Peterson, S. L. Ray, C. H. Chan, and R. Mittra, “Numerical implementation of the conjugate gradient method and the CG-FFT for electromagnetic scattering,” in *Application 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) [CrossRef] .

**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) [CrossRef] [PubMed] .

**37. **L. Knockaert, “A note on the relationship between the conjugate gradient method and polynomials orthogonal over the spectrum of a linear operator,” IEEE Trans. Antennas Propag. **AP-35**, 1089–1091 (1987) [CrossRef] .

**38. **A. F. Peterson, C. F. Smith, and R. Mittra, “Eigenvalues of the moment-method matrix and their effect on the convergence of the conjugate gradient algorithm,” IEEE Trans. Antennas Propag. **36**, 1177–1179 (1988) [CrossRef] .

**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) [CrossRef] .

**40. **P. Joly, “Résolution de systèmes linéaires avec plusieurs members par la méthode du gradient conjugué,” Tech. Rep. R-91012 (Publications du Laboratoire d’Analyse Numérique, Université 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, III., 1992).

**42. **P. J. Flatau, T. Schneider, and F. Evans, ccg-pak—fortran. Conjugate gradient package for solving complex matrix equations (1993).Available from pflatau@ucsd.edu.

**43. **D. T. Borup and O. 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) [CrossRef] .

**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,” Electromagnetics **5**, 99–122 (1985) [CrossRef] .

**45. **P. Barber and C. Yeh, “Scattering of electromagnetic waves by arbitrarily shaped dielectric bodies,” Appl. Opt. **14**, 2864–2872 (1975) [CrossRef] [PubMed] .

**46. **D.-S. Wang and P. W. Barber, “Scattering by inhomogeneous nonspherical objects,” Appl. Opt. **18**, 1190–1197 (1979) [CrossRef] [PubMed] .

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

**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) [CrossRef] .

**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) [CrossRef] .

**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) [CrossRef] .

**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) [CrossRef] .

**52. **K. A. Fuller, “Optical resonances and two-sphere systems,” Appl. Opt. **30**, 4716–4731 (1991) [CrossRef] [PubMed] .

**53. **D. W. Mackowski, “Analysis of radiative scattering for multiple sphere configurations,” Proc. R. Soc. London Ser. A **433**, 599–614 (1991) [CrossRef] .

**54. **G. W. Kattawar and T. J. Humphreys, “Electromagnetic scattering from two identical pseudospheres,” in *Light Scattering by Irregularly Shaped Particles*, D. W. Schuerman, ed. (Plenum, New York, 1980), pp. 177–190 [CrossRef] .

**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 Å extinction profile,” Astrophys. J. **414**, 632–645 (1993) [CrossRef] .

**59. **M. F. Iskander, H. Y. Chen, and J. E. Penner, “Resonance optical 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) [CrossRef] [PubMed] .

**61. **R. A. West and P. H. Smith, “Evidence for aggregate particles in the atmospheres of Titan and Jupiter,” Icarus **90**, 330–333 (1991) [CrossRef] .

**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. Greenberg, “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) [CrossRef] .

**66. **J. M. Perrin and P. L. Lamy, “On the validity of effective- medium theories in the case of light extinction by inhomogeneous dust particles,” Astrophys. J. **364**, 146–151 (1990) [CrossRef] .

**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) [CrossRef] .

**68. **K. F. Evans and J. Vivekanandan, “Multiparameter radar and microwave radiative transfer modeling of nonspherical atmospheric ice particles,” IEEE Trans. Geosci. Remote Sensing **28**, 423–437 (1990) [CrossRef] .

**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).