The discrete dipole approximation: an overview and recent 
developments 
M.A. Yurkina,b,∗ and A.G. Hoekstraa
a Section Computational Science, Faculty of Science, University of Amsterdam, Kruislaan 
403, 1098 SJ, Amsterdam, The Netherlands 
b Institute of Chemical Kinetics and Combustion, Siberian Branch of the Russian Academy of 
Sciences, Institutskaya Str. 3, 630090, Novosibirsk, Russia 
Abstract 
We present a review of the discrete dipole approximation (DDA), which is a general method 
to simulate light scattering by arbitrarily shaped particles. We put the method in historical 
context and discuss recent developments, taking the viewpoint of a general framework based 
on the integral equations for the electric field. We review both the theory of the DDA and its 
numerical aspects, the latter being of critical importance for any practical application of the 
method. Finally, the position of the DDA among other methods of light scattering simulation 
is shown and possible future developments are discussed. 
Keywords: discrete dipole approximation, review, light scattering simulation 
                                                 
∗ Corresponding author: Tel.: +31-20-525-7562; fax: +31-20-525-7490. 
  E-mail address: myurkin@science.uva.nl 
Contents 
1 Introduction ......................................................................................................................... 2 
2 General framework.............................................................................................................. 3 
3 Various DDA models .......................................................................................................... 7 
3.1 Theoretical base of the DDA........................................................................................ 7 
3.2 Accuracy of DDA simulations ................................................................................... 13 
3.3 The DDA for clusters of spheres................................................................................ 16 
3.4 Modifications and extensions of the DDA................................................................. 18 
4 Numerical considerations.................................................................................................. 19 
4.1 Direct vs. iterative methods........................................................................................ 19 
4.2 Scattering order formulation ...................................................................................... 22 
4.3 Block-Toeplitz ........................................................................................................... 23 
4.4 FFT............................................................................................................................. 24 
4.5 Fast multipole method................................................................................................ 24 
4.6 Orientation averaging and repeated calculations ....................................................... 25 
5 Comparison of the DDA to other methods ....................................................................... 27 
6 Concluding remarks .......................................................................................................... 28 
Acknowledgements .................................................................................................................. 28 
Appendix. Description of used acronyms and symbols ........................................................... 28 
References ................................................................................................................................ 31 
1 Introduction 
The discrete dipole approximation (DDA) is a general method to compute scattering and 
absorption of electromagnetic waves by particles of arbitrary geometry and composition. 
Initially the DDA was proposed by Purcell and Pennypacker (PP) [1], who replaced the 
scatterer by a set of point dipoles. These dipoles interact with each other and the incident 
field, giving rise to a system of linear equations, which is solved to obtain dipole 
polarizations. All the measured scattering quantities can be obtained from these polarizations. 
The DDA was further developed by Draine and coworkers [2-5], who popularized the method 
by developing a publicly available computer code DDSCAT [6]. Later it was shown that the 
DDA also can be derived from the integral equation for the electric field, which is discretized 
by dividing the scatterer into small cubical subvolumes. This derivation was apparently first 
performed by Goedecke and O'Brien [7] and further developed by others (see, for instance, 
[8-11]). It is important to note that the final equations, produced by both lines of derivation of 
the DDA are essentially the same. The only difference is that derivations based on the integral 
equations give more mathematical insight into the approximation, thus pointing at ways to 
improve the method, while the model based on point dipoles is physically clearer. 
The DDA is called the coupled dipole method or approximation by some researchers 
[12,13]. There are also other methods, such as the volume integral equation formulation [14] 
and the digitized Green’s function (DGF) [7], which were developed completely 
independently from PP. However, later they were shown to be equivalent to DDA [8,15]. In 
this review we will use the term DDA to refer to all such methods, since we describe them in 
terms of one general framework. However, it is difficult to separate unambiguously the DDA 
from other similar methods, based on the volume integral equations for the electromagnetic 
fields, such as a broad range of method of moments (MoM) with different bases and testing 
functions [16-19]. In our opinion, one fundamental aspect of the DDA is that the solution for 
the “physically meaningful” internal fields or their direct derivatives, e.g. polarization, plays 
an integral role in the process. In other words, any DDA formulation can be interpreted as 
replacing a scatterer by a set of interacting dipoles; this is further discussed in Section 2. An 
example of method that is not considered DDA is the MoM with higher-order hierarchical 
Legendre basis functions [17]. 
The DDA is a popular method in the light-scattering community and it has been 
reviewed by several authors. An extensive review by Draine and Flatau [4] covers almost all 
DDA developments up to 1994. A more recent review by Draine [5] mainly concerns 
applications and numerical considerations. DDA theory was discussed together with other 
methods for light scattering simulations in reviews by Wriedt [20], Chiappetta and Torresani 
[21], and Kahnert [15] and in books by Mishchenko et al. [22] and Tsang et al. [23]. Jones 
[24] placed the DDA in context of different methods with respect to particle characterization. 
However, many important DDA developments since 1994 are not mentioned in any of these 
manuscripts. Those that are mentioned are usually considered as side-steps, and are not placed 
into a general framework. Moreover, to the best of our knowledge numerical aspects of the 
DDA have never been reviewed extensively – each paper discusses only a few particular 
aspects. In this review we try to fill this gap. 
A general framework is developed in Section 2 to ease the further discussion of different 
DDA models. This framework is based on the integral equation because it allows a uniform 
description of all the DDA development. However, connection to a physically clearer model 
of point dipoles is discussed throughout the section. The sources of errors in the DDA 
formulation are also discussed there. 
In Section 3 the physical principles of the DDA are reviewed and results of different 
models are compared. In Subsection 3.1 different improvements of polarizabilities and 
interaction terms are reviewed from a theoretical point of view. Different expressions for Cabs 
also are discussed. Comparison of simulation results using different formulations is given in 
Subsection 3.2. Subsection 3.3 covers the special case of a cluster of spheres that allows 
particular improvements and simplifications. In Section 3.4 different significant modifications 
are reviewed, which do not fall completely into the general framework described in Section 2. 
Enhancements of the DDA for some special purposes also are discussed. 
Different numerical aspects of the DDA are reviewed in Section 4. These are concerned 
primarily with solving very large systems of linear equations (Subsection 4.1). Subsection 4.2 
describes the simplest iterative procedure to solve DDA linear system, which has a clear 
physical meaning. The special structure of the DDA interaction matrix for a rectangular grid 
and its application to decrease computational costs are described in subsections 4.3 and 4.4 
respectively. General methods to accelerate calculations, which do not require a rectangular 
grid, are discussed in Subsection 4.5. Subsection 4.6 covers special techniques to increase the 
efficiency of repeated calculations (e.g. in orientation averaging). 
A numerical comparison of the DDA with other methods is reviewed in Section 5; its 
strong and weak points are discussed. Section 6 concludes the review and discusses future 
development of the DDA. 
2 General framework 
The )iexp( tω−  time dependence of all fields is assumed throughout this review. The scatterer 
is assumed dielectric but not magnetic (magnetic permittivity 1=μ ). The electric permittivity 
is assumed isotropic to simplify the derivations; however, extension to arbitrary dielectric 
tensors is straightforward.1
The general form of the integral equation governing the electric field inside the 
dielectric scatterer is the following [8,15]: 
)()(),(),()()(),(d)()( 00
rErrLrMrErrrGrErE χχ VVr
∂−+′′′′+= ∫ , (1)
                                                 
1 In most formulae scalar values can be replaced directly by tensors, but there are exceptions. Extensions of 
DDA to optically anisotropic scatterers are discussed in Section 3.4. 
where Einc(r) and E(r) are the incident and total electric field at location r; 
πεχ 4)1)(()( −= rr  is the susceptibility of the medium at point r (ε(r) – relative 
permittivity). V is the volume of the particle, i.e., the volume that contains all points where the 
susceptibility is not zero. V0 is a smaller volume such that , VV ⊂0 00 \ VV ∂∈r . ),( rrG ′  is 
the free space dyadic Green’s function, defined as 
−=∇∇+=′ 222
i1ˆˆ)iexp()iexp(ˆˆ),(
k IIIrrG , (2)
where I  is the identity dyadic, ck ω=  is the free space wave vector, rrR ′−= , R=R , 
and  is a dyadic defined as  (μ and ν  are Cartesian components of the vector 
or tensor). M is the following integral associated with the finiteness of the exclusion volume 
RR ˆˆ νμμν RRRR =ˆˆ
( )∫ ′−′′′′=
)()(),()()(),(d),( s30
rV rErrrGrErrrGrM χχ , (3)
where ),(s rrG ′  is the static limit ( ) of 0→k ),( rrG ′ : 
−−=∇∇=′
11ˆˆ),(
IrrG . (4)
L  is the so-called self-term dyadic: 
rV rL , (5)
where  is an external normal to the surface ∂Vn′ˆ 0 at point r'. L  is always a real, symmetric 
dyadic with trace equal to 4π [25]. It is important to note that L  does not depend on the size 
of the volume V0, but only on its shape (and location of the point r inside it). On the contrary, 
M does depend on the size of the volume, moreover it approaches zero when the size of the 
volume decreases [8] (if both χ(r) and E(r) are continuous inside V0). 
When deriving Eq. (1) the singularity of the Green’s function has been treated explicitly, 
therefore it is preferable to the commonly used formulation [8,15]: 
∫ ′′′′+=
r )()(),(d)()( 3inc rErrrGrErE χ . (6)
Moreover, Yanghjian noted [25] that there exist several methods for treating the singularity in 
Eq. (6) leading to different results. He also proved that the derivation of Eq. (6) is false in the 
vicinity of the singularity of ),( rrG ′ . Hence it can be considered correct only if the 
singularity is then treated in a way similar to that of Lakhtakia [8], resulting in the correct 
Eq. (1). 
Discretization of Eq. (1) is done in the following way [15]. Let , U
= /0=ji VV I  
for ji ≠ ; N denotes number of subvolumes.2 Although the formulation is applicable to any 
set of subvolumes Vi, in most applications standard (equal) cells are used. Then the shape of 
the scatterer cannot always be described exactly by such standard cells. Hence, the 
discretization may be only approximately correct. Assuming iV∈r  and choosing iVV =0 , 
Eq. (1) can be rewritten as 
)()(),(),()()(),(d)()( 3inc rErrLrMrErrrGrErE χχ ii
∂−+′′′′+= ∑ ∫
. (7)
                                                 
2 In the framework of the DDA we usually call a subvolume a dipole. 
The set of Eq. (7) (for all i) is exact. Further, one fixed point ri inside each Vi (its center) is 
chosen and  is set. In many cases the following assumptions can be made: irr =
)()()()(),(d3 jjijj
rErGrErrrG χχ =′′′′∫ , (8)
)()(),( iiiiiV rErMrM χ= , (9)
which state that integrals in Eq. (7) linearly depend upon the values of χ and E at point ri. 
Eq. (7) can then be rewritten as 
( ) iiii
jjjijii V ELMEGEE χχ −++= ∑
inc , (10)
where , , )( ii rEE = )(
incinc
ii rEE = )( ii rχχ = , ),( iii V rLL ∂= . 
The usual approximation [15] is to consider E and χ constant inside each subvolume: 
iii V∈== rrErE for)(,)( χχ , (11)
which automatically implies Eqs. (8), (9) and 
( )∫ ′−′′=
iii r ),(),(d
s3)0( rrGrrGM , (12)
∫ ′′=
ij rV
1 3)0( rrGG . (13)
Superscript (0) denotes approximate values of the dyadics. A further approximation, which is 
used in almost all formulations of the DDA, including e.g. [8], is  
),()0( jiij rrGG = . (14)
This assumption is made implicitly by all formulations that start by replacing the scatterer 
with a set of point dipoles. It is important to note that Eq. (10) and derivations resulting from 
it require weaker assumptions (Eqs. (8), (9)) than imposed by Eq. (11) and, moreover, 
Eq. (14). It is possible to formulate the DDA based on Eq. (10), e.g. the Peltoniemi 
formulation [26] that is described in Section 3.1. We postulate Eq. (10) as a distinctive feature 
of the DDA, i.e. a method is called the DDA if and only if its main equation is equivalent to 
Eq. (10) with any Vi, χi, iM , iL , and ijG . 
Kahnert [15] distinguished the DDA from the MoM by the fact that the MoM solves 
directly Eq. (10) for unknown Ei, while the DDA seeks not the total, but the exciting electric 
fields 
( )( ) selfexc iiiiiii EEEMLIE −=−+= χ , (15)
( ) iiiii ELME χ−=self , (16)
where  is the field induced by the subvolume on itself. Eq. selfiE (10) is then equivalent to  
jjijii
excexcinc EαGEE , (17)
where iα  is the polarizability tensor defined as 
( )( ) 1−−+= iiiiii V χχ MLIα . (18)
However, an alternative formulation of the DDA exists [4] seeking a solution for unknown 
polarizations Pi: 
iiiiii V EEαP χ==
exc , (19)
jijiii PGPαE
1inc . (20)
It is important to note that Pi, defined by Eq. (19), is only an approximation to the polarization 
of the subvolume Vi. This approximation is exact only under the assumption of Eq. (11), while 
the formulation itself does not require it. The formulation, using Eq. (20), can be thought as 
an intermediary between the DDA and the MoM as classified by Kahnert [15], therefore 
revealing complete equivalence of these two formulations. The special structure of the matrix 
ijG  makes Eq. (20) preferable over Eqs. (10), (17) to find a numerical solution. This is 
discussed in Section 4. 
Lakhtakia [8] classified strong and weak forms of the DDA as those accounting for or 
neglecting iM  respectively. The weak form approaches the strong form when the size of the 
cell decreases, because iM  approaches zero. For a cubical cell Vi and with ri located at the 
center of the cell, iL  can be calculated analytically yielding [25] 
=i . (21)
Using Eq. (18), this results in the well-known Clausius-Mossotti (CM) polarizability (used 
originally by Purcell and Pennypacker [1]) for the weak form of the DDA: 
ii d ε
α IIα , (22)
where )( ii rεε = , and d is the size of the cubical cell. 
After the internal electric fields are determined, the scattered fields and cross sections 
can be calculated. The scattered fields are obtained by taking the limit ∞→r  of the integral 
in Eq. (1) (see e.g. [7]): 
)iexp(
)(sca nFrE
= , (23)
where rrn =  is the unit vector in the scattering direction, and F is the scattering amplitude: 
∑∫ ′′⋅′−′−−=
krnnk )()()iexp(d)ˆˆ(i)( 33 rErnrInF χ . (24)
All other differential scattering properties, such as amplitude and Mueller scattering matrices, 
and asymmetry parameter >< θcos  can be derived from F(n), calculated for two incident 
polarizations [27]. Radiation forces also can be calculated [28-30]. Consider an incident 
polarized plane wave3
)iexp()( 0inc rkerE ⋅= , (25)
where , a is the incident direction, and ak k= 10 =e . The scattering cross section Csca is [27] 
2sca )(d
C . (26)
Absorption and extinction cross sections (Cabs, Cext) are derived [7,14] directly from the 
internal fields: 
( )∑ ∫ ′′′=
abs )()(Imd4 rErχπ , (27)
[ ]( ) ( )∗∗ ⋅=′⋅′′′= ∑∫ 02inc3ext )(Re4)()()(Imd4 eaFrErEr krkC i Vi
χπ , (28)
where * denotes a complex conjugate. Conservation of energy necessitates that 
absextsca CCC −= . (29)
However, as was noted by Draine [2], use of Eq. (29) for evaluation of Csca can lead to larger 
errors than Eq. (26), especially when . scaabs CC >>
The easiest way to express Eqs. (24) and (27) in terms of the internal fields in the 
subvolumes centers is to assume Eq. (11), yielding 
∑ ∫ ⋅′−′−−=
krnnk )iexp(d)ˆˆ(i)( 33)0( nrEInF χ , (30)
                                                 
3 DDA can be used for any incident wave, e.g. Gaussian beams [31]; however, we do not discuss this here. 
∑∑ ∗==
iii kVkC )Im(4)Im(4
abs EPE πχπ . (31)
Further approximation of Eq. (30), leaving only the lowest order expansion of the exponent 
around ri, leads to 
∑ ⋅−−−=
ii knnk )iexp()ˆˆ(i)(
3)0( nrPInF , (32)
which together with Eq. (28), leads to 
( )∑ ∗⋅=
inc)0(
ext Im4 EPπ . (33)
Eqs. (32) and (33) are identical to those used by Purcell and Pennypacker [1] and then by 
Draine [2], while expressions for Cabs (compared to Eq. (31)) are slightly different. These 
differences are discussed in Subsection 3.1. Unfortunately, many researchers do not specify 
explicitly how the scattering quantities are obtained from the computed internal fields or 
polarizations. Those who do usually use Draine’s prescription (Eqs. (26), (32), (33), and 
(35)). 
Errors of the formulation can be classified as associated with the finite cell size d 
(discretization errors), and with approximating the particle shape with a set of standard cells, 
e.g. cubical (shape errors). Discretization errors result from considering E constant inside 
each cell and the approximate evaluation of iM  and ijG . Shape errors also can be considered 
as resulting from the assumption of constant χ and E inside bordering cells, which is false 
since the edge of the particle crosses these cells. On the other hand, shape errors can be 
viewed as a difference of the results for the exact particle shape and for that comprised of the 
set of standard cells. Both errors approach zero when ∞→N , while the geometry of the 
scatterer and parameters of the incident field are fixed. However, the same does not apply if 
 while N is fixed, i.e. the DDA is not exact in the long-wavelength limit. Moreover, 
both errors are sensitive to the size of the scatterer in the resonance region (see discussion in 
Subsection 
3.2). The behavior of these errors was studied by Yurkin et al. [32]. 
3 Various DDA models 
3.1 Theoretical base of the DDA 
Since the original manuscript by Purcell and Pennypacker [1], many attempts have been made 
to improve the DDA. The first stage (1988-1993) of these improvements was reviewed by 
Draine and Flatau [4]. It has been noted [2] that Eq. (22) does not satisfy energy conservation, 
and results obtained using this formulation do not satisfy the optical theorem. Based on the 
well-known [33] “radiative reaction” (RR) electric field, a correction to the polarizability for a 
finite dipole was added [2]: 
i)32(1 α
= . (34)
Draine [2] also proposed the following expression for the absorption cross section: 
[ ]∑ ⋅−⋅=
iiii kkC
*3*exc)0(
abs )32()Im(4 PPEPπ , (35)
derived from Eq. (29) applied to a single point dipole. The PP formulation uses Eq. (35) 
without the second part. It can be verified that Eq. (35) results in zero absorption for any 
scatterer if the polarizability is of the following form: 
IAα 31 i)32( kii −=
− , Hii AA = , (36)
where H denotes the conjugate transpose of a tensor. For real refractive index m, RR and all 
other expressions specified below result in α  satisfying Eq. (36), which makes Eq. (35) 
clearly favorable over e.g. the PP formulation. It must be noted however that the original PP 
formulation, where CM polarizability was used, also results in zero absorption for real m. 
The correction in Eq. (34) is ( )3)(O kd . Several other corrections of ( )2)(O kd  have been 
proposed. The first one was proposed by Goedecke and O’Brien [7] and independently in two 
other manuscripts [34,35]. They started from Eqs. (10)-(12) and used the following 
simplifying fact for a cubical cell (also valid for spherical cells), resulting from symmetry: 
)(d IRRf
RRf , (37)
where the origin is in the center of the cube. Eq. (37) is valid for any f(R) that has a singularity 
of less than third order for , i.e. the integrals on both sides are defined. They obtained 0→R
32)0( )iexp(d
Rki IM . (38)
By expanding exp(ikR) in Taylor series one can obtain 
( )⎟⎟
++= ∫ 423
2)0( Oi
ki IM . (39)
The remaining integral was evaluated by approximating the cube by a volume-equivalent 
sphere, resulting in 
( )( )432DGF1)0( )(O))i(32()( kdkdkdbi ++= IM , (40)
611992.1)34( 31DGF1 ≈= πb . (41)
An exact evaluation, obtained without expanding the exponent, of Eq. (38) for the equivolume 
sphere with radius 31)43( πda =  was performed by Livensay and Chen [36] and 
implemented into the DGF formulation of the DDA by Hage and Greenberg [14,35] and later 
Lakhtakia [37]: 
[ ]1)iexp()i1()38()0( −−= kakai IM π . (42)
In terms of the first two orders of expansion, this yields an identical result as Eq. (40). Finally 
the polarizability is obtained as 
( )( )32DGF13CM
)i()32()(1 kdkdbd +−
α . (43)
We denote the method based on Eq. (42) as LAK. Differences between LAK and DGF should 
be noticeable only for large values of kd. 
Dungey and Bohren [38], using results by Doyle [39], proposed the following treatment 
of the polarizability. First, each cubic cell is replaced by the inscribed sphere that is called a 
dipolar subunit with a higher relative electric permittivity εs as determined by the Maxwell-
Garnett effective medium theory [27]: 
f , (44)
where 6π=f  is the volume filling factor. Other effective medium theories also may be used 
[40]. Next, the dipole moment of the equivalent sphere is determined using the Mie theory, 
and the polarizability is defined as [39] 
i= , (45)
where α1 is the electric dipole coefficient from the Mie theory (see e.g. [41]): 
)()()()(
)()()()(
sssssss
sssssss
1 xmxxxmm
xmxxxmm
α , (46)
where ψ, ξ are Riccati-Bessel functions; 2s kdx =  and ss ε=m  are the size parameter and 
the relative refractive index of the equivalent sphere. We denote this formulation for the 
polarizability as the a1-term method (note that this terminology was introduced later [42]). It 
has the particular property that 1constCMM ≠→αα  when , contrary to all other 
polarization prescription, for which this ratio approaches 1. It should be noted that the Mie 
theory is based on the assumption that the external electric field is a plane wave. In most 
applications of the DDA this is true for the incident electric field, but not for the field created 
by other subvolumes. Therefore the a
1-term method is expected to be correct only for very 
small cell size. Hence it is not clear whether this method has advantages even compared to 
CM. On the other hand, this method may be more justified for clusters of small spheres, 
where each sphere can be considered as a dipole (see Subsection 3.3). 
Draine and Goodman [3] pointed out that considering electric fields constant for 
evaluating integrals over a cell introduces errors of order ( )( )2O kd . This represents a problem 
for many polarizability corrections, based on integral equations. Draine and Goodman 
approached this problem from a different angle. They determined the optimal polarizability in 
the sense that an infinite lattice of point dipoles with such polarizability would lead to the 
same propagation of plane waves4 as in a medium with a given refractive index. This 
polarizability was called LDR (Lattice Dispersion Relation) and is, as expected, CM plus 
high-order corrections. These corrections in turn depend on the direction of propagation a and 
the polarization of the incident field e0: 
( ) ( )[ ]322LDR32LDR2LDR13CM
)i()32()(1 kdkdSmbmbbd +++−
α , (47)
8915316.1LDR1 ≈b , , , 1648469.0
2 −≈b 7700004.1
3 ≈b (48)
( )∑=
20eaS . (49)
We use a reverse sign convention in the denominator of Eq. (47) and the LDR coefficients as 
compared to the original paper [3]. 
Recently it has been shown [43] that the LDR derivation is not completely accurate, 
since the resulting dipole moment does not satisfy the transversality condition, for which a 
correction was proposed. This corrected LDR (CLDR) differs principally in the fact that the 
polarizability tensor can not be made isotropic but only diagonal [43], though not dependent 
on the incident polarization: 
( ) ( )[ ]3222LDR32LDR2LDR13CM
)i()32()(1 kdkdambmbbd +++−
α . (50)
Another flaw of LDR is that it is evidently not correct for dipoles near the particle surface. 
However, it is not clear how to evaluate the effect of these mistreated surface dipoles on the 
overall results, e.g. on the scattering cross section. 
Further improvement of the DDA was initiated by Peltoniemi [26] (PEL) who showed 
that the term M(Vi) in Eq. (7) can be evaluated exactly up to the third order of kd by 
expanding the term )()( rEr ′′χ  under the integral in a Taylor series over the point irr =′ , 
yielding 
( ) ( ,)(O3i3)iexp(d
)iexp(
EkdERRRRkRRk
ERkRRk
REMVM
ντρτρνμ
+∂∂−+−
∂−++=
                                                
 (51)
where χ, E and their derivatives are all considered at the point ri. Eq. (51) is correct up to the 
third order of kd since the third term in the Taylor series vanishes because of symmetry. For 
spherical Vi of radius a, the integrals can be evaluated exactly [26] in a way similar to 
4 with certain direction of propagation and polarization state. 
obtaining Eq. (42), but only terms of less than fourth order of kd are significant, which results 
( EkaakakaVi χχχχπ 42232 )(O)(10
)( +⎥
⎛ ⋅∇∇−∇−⎟
⎛ += EEEM ). (52)
If χ is constant inside the cell then the Maxwell equations state that 
EE 222 km−=∇ , 0=⋅∇ E . (53)
Hence Eq. (9) is valid up to the third order of ka and 
( )[ ]322 )(i)32()()101(1)34( kakami ++= IM π . (54)
Piller and Martin [44] proposed using sampling theory to evaluate the integrals in 
Eq. (1). The electric field and the susceptibility is sampled: 
∑ −′=′′
iiih )()()()()(
r rErrrrEr χχ , (55)
where hr(r) is the impulse response function of an antialiasing filter defined as 
)cos()sin(
qrqrqr
=r , (56)
where dq π2= . Eq. (1) is then transformed to Eq. (10) with the so-called filtered Green’s 
function, defined as 
∫ −′′′=
r3 )(),(d
ij hrV
rrrrGG . (57)
Eq. (57) can be viewed as a generalization of Eq. (13). The latter is obtained if a pulse 
function is considered instead of hr. The integral in Eq. (57) is evaluated analytically [44], 
taking V0 to be infinitesimally small. The filtered Green’s function does not have a singularity 
when , therefore ji rr = iiii V GM = . It was shown that the Fourier spectrum of E(r) lies on a 
sphere with radius m(r)k, if m is constant in the vicinity of r. Therefore at least two sampling 
points per wavelength in the scatterer are required. The susceptibility is also filtered, either by 
a mean value filter or a more complicated one, e.g. a Hanning window. This approach is 
called FCD (filtered coupled dipoles), and a computer code library for evaluation of filtered 
Green’s function is available [45]. 
Chaumet et al. [11] proposed direct integration of the Green’s tensor (IT) in Eqs. (12), 
(13). A Weyl expansion of the Green’s tensor is performed, transforming it to a form allowing 
efficient numerical computation of the self-term ( LM − ). They also proposed a correction to 
the second term in Draine’s expression for Cabs (Eq. (35)). Extension of their results to a non-
isotropic self-term is 
( ) ( )[ ]∑ −⋅+⋅=
iiiiiii VkC /)(ImIm4
***exc)0(
abs PLMPEPπ
The corrected second term is based on radiation energy of a finite dipole [11]: , in 
contrast to a point dipole used in the derivation of Eq. 
)Im( self ∗⋅ ii PE
(35). One can see that Eqs. (58) and 
(31) are equivalent. Moreover, both of them are equivalent to Eq. (35) if and only if 
IAM iii Vk
3i)32(+= , Hii AA = . (59)
This condition is similar, but not equivalent, to Eq. (36) and is always satisfied for RR, DGF, 
and LAK. Other polarizability prescriptions satisfy Eq. (59) for real m, then both Eqs. (58) 
and (35) result in zero absorption. 
Rahmani, Chaumet, and Bryant [46] proposed a new method (RCB) to determine 
polarizability based on the known solution of the electrostatic problem for the same scatterer. 
In the static limit the electric field at any point is linearly related to the incident field 
)()()( 01 rErCrE −= . (60)
Substituting Eq. (60) into Eq. (20) with the static Green’s tensor, one can obtain the 
polarizability, which would give an exact solution in the static limit, as 
1RCB −= iiii V Λχα , (61)
ii CCrrGCΛ
1),( −
∑+= χ , (62)
where )( ii rCC = . This static polarizability then replaces the CM polarizability, and the RR 
(Eq. (34)) is applied to it [46] to obtain the final polarizability for DDA simulations. It was 
later shown that RCB polarizabilities differ significantly from CM only for dipoles closer than 
2d to the interface [47]. 
In their next manuscript [48] Rahmani et al. stated that the previous derivation is correct 
only if the tensor C  is constant inside the particle (e.g. for ellipsoids), since otherwise the 
polarizability tensor obtained from Eq. (61) is generally not symmetric, which is physically 
impossible in the static case. This shows that a particle with a non-constant C  is not 
equivalent to any set of physical point dipoles even in the static regime. However, it is 
equivalent to a set of non-physical dipoles with an asymmetric polarizability. Therefore, the 
polarization defined by Eq. (61) formally can be used, by itself or with RR, even when C  is 
not constant. 
Collinge and Draine [47] empirically combined the RCB prescription with CLDR to get 
the surface-corrected LDR (SCLDR): 
( )( ) 13RCBRCBSCLDR −−= BαIαα d , (63)
where B  is the correction matrix (analogous to Eq. (50)): 
( )[ ]3222LDR32LDR2LDR1 )i()32()( kdkdambmbbB +++= μμνμν δ . (64)
All methods based on the paper by Rahmani et al. [46] are initially limited to very specific 
shapes of the scatterer (ellipsoids, infinite slabs and cylinders). Expansion of its applicability 
to other shapes is debatable [48] and would anyway require a preliminary solution of the 
electrostatic problem for the same shape, which is generally not trivial. 
All DDA formulations are schematically depicted in Fig. 1, which also shows 
interrelations between them. Some formulations can be compared unambiguously in terms of 
theoretical soundness: one is an improvement of the other, i.e. it employs fewer 
approximations. Such formulations are depicted in the same column on Fig. 1, while others 
cannot be compared directly with each other; they give rise to different columns. Comparison 
between formulations from different columns can and has been made almost exclusively 
empirically by comparing the accuracy of the simulation results (see Subsection 3.2). 
All the above techniques are aimed at reducing discretization errors; only a few aim at 
reducing shape errors. Some of them employ adaptive discretization (different dipole sizes) to 
better describe the shape of the scatterer (see Subsection 3.4). Another approach is to average 
susceptibility in boundary subvolumes. The simplest averaging using the Lorentz-Lorenz 
mixing rule was proposed by Evans and Stephens [49] for the case of the boundary between 
the scatterer and its surrounding medium 
3434 e
, (65)
where  is the effective susceptibility, and f is the volume fraction of the subvolume actually 
occupied by scatterer. 
A more advanced averaging, called the weighted discretization (WD), was proposed by 
Piller [13]. It modifies the susceptibility and self-term of the boundary subvolume.5 The 
particle surface, crossing the subvolume Vi, is assumed linear and divides the subvolume into 
two parts: the principal  that contains the center and a secondary  with susceptibilities piV
                                                 
5 any subvolume that has non-zero intersection with both the scatterer and the outer medium. All such 
subvolumes are accounted for. 
 Integral Eq. (1)  
discretization 
(no assumptions)
Eq. (7) 
General 
formulation of 
DDA – Eq. (20) 
Eqs. (8), (9) 
DGF, LAK 
Eqs. (11) 
Eq. (14)
(weak form)
CLDR 
a1 term 
SCLDR FCD 
sampling with 
antialiasing 
filter 
removing 
antialiasing 
filter 
improving polarizability starting 
from dipole formulation 
complies 
Eq. (14) 
simplifies to 
Fig. 1. Scheme of interrelation between the different DDA models discussed in Section 3.1. Arrows 
down correspond to assumptions employed. Vertical position of the method qualitatively corresponds 
to its accuracy (higher = better), however methods in different columns cannot be compared directly. 
iχ ,  and electric fields ,  respectively. Electric fields are considered constant 
inside each part and related to each other via a boundary condition tensor 
iχ ii EE ≡
iT : 
iii ETE =
s . (66)
Then the total polarization of the subvolume can be evaluated as follows: 
iiiiiiiii
i VVVr
EEErErP essspp3 )()(d χχχχ =+=′′′= ∫ , (67)
( ) iiiiiii VVV TI ssppe χχχ += . (68)
The susceptibility of the boundary subvolume is replaced by an effective one. 
The effective self-term is evaluated directly starting from Eq. (3), considering χ and E 
constant inside each part: 
( ) ( ) ii
rr TrrGrrGrrGrrGM ss3ps3ee
),(),(d),(),(d χχχ ∫∫ ′−′′+′−′′= . (69)
Piller [13] evaluated the integrals in Eq. (69) numerically. The final equations are the same as 
Eq. (20), where polarizabilities are obtained from Eq. (18) using effective susceptibilities and 
self-terms for boundary subvolumes. Hence, WD does not modify the general numerical 
scheme. 
Currently, there are no rigorous theoretical reasons for preferring one formulation over 
others. However, theoretical analyses of DDA convergence when refining discretization 
recently conducted by Yurkin et al. [32], showed that IT and WD significantly improve the 
convergence of shape and discretization errors, respectively. Experimental verification of 
these theoretical conclusions is still to be performed. 
Table 1. Accuracy of different DDA formulations for a sphere.a
Value Method x a/d y m Error, % Ref.
Cext a1-term 1÷2 2÷4
c 0.65 
0.85 
1.33+0.05i 
1.7+0.1i 
CSec, S11 LAK 9 
0.44 
0.42 
0.51 
1.05 
1.33+0.01i 
2.5+1.4i 
0.05, 37 
0.5, 35 
4, 15 
Csca, Cabs DGF ≤3.2c 16 ≤1 4+3i 5, 10÷30 [3]
CSec LDR ≤8c 16 ≤0.5, ≤0.1 m-1≤1 1, 2 
Csca 
LDR ≤7c 16 ≤1 
≤0.5, ≤1 
2+i 1.5 
3, 4 
CSec ≤16c
25 ≤1 1.6+0.0008i 
2.5+0.02i 
LDR [51]
[4]eCSec 
LDR any any ≤1 |m|≤3 5 
20÷30 
Csca LDR ≤10c 16 kd≤0.63 0.69 
0.41 
0.29 
[148]
S11 LDR ≤10c 24 kd≤0.42 0.69 
0.41 
Cext,  RMS11S
3.2 Accuracy of DDA simulations 
Over the years many results on the accuracy of DDA simulations have been published. It is, 
however, generally hard to systematically compare the relevant manuscripts because they all 
use different independent parameters, such as the size parameter x, refractive index m, or 
discretization, as a function of which the error is measured. We will describe discretization by 
the parameter kdmy =  or . The former is used wherever possible; however, 
in some cases a description of results is more straightforward in terms of y
kdmy )Re(Re =
Re. Accuracy results 
LDR 20÷160 
20÷130 
20÷60 
20÷30 
20÷30 
32÷256c
40÷256c
48÷128c
56÷80c
64÷88c
0.61÷0.65d
0.56÷0.64d
0.58÷0.65d
0.57÷0.60d
0.56÷0.62d
0.62d
1.05 
0.04, 38 
0.4, 23 
1, 59 
4.4, 56 
5.7, 105 
2.0, 86 
[113]
Ψ FCD π, 2π 2.8, 5.6c 1.7 1.5 1 [44]
WD-FCD 0.5÷3.2c
1.5÷3.8c
0.9÷1.5c
yRe=0.63 |m|<7b
|m|<2.5b
|m|<4b
[10]Ψ 
IT ≤5.2c CSec 
≤2.1c
≤1.1c
8 ≤1 1.5+0.3i 
3.5+1.4i 
7.1+0.7i 
Cabs 
CSec RCB-RR ≤8.2c 
≤7.5c 
≤5.9c 
≤3.4c 
≤1.3c
16 ≤1 1.8+0.4i 
1.9+i 
2.5+i 
2.5+4i 
7.4+9.4i 
CSec SCLDR 
SCLDR 
≤7.2c 
≤1.5c 
≤1.5c
12 ≤0.8 1.33+0.1i 
5+4i 
5+4i 
a All errors are relative. CSec denotes the maximum error over all cross sections, S11 and  correspond to 
maximum and root mean square error over the range of scattering angles, Ψ is the normalized mean error of the 
far-field electric fields [44]. In some cases two errors are shown in one cell separated by a coma. They 
correspond to two values of one of the parameters in the same row. 
b approximate description of the range. 
c this value is determined by other values in the same row. 
d this value is slightly different for different size parameters. 
e this corresponds to the “rule of thumb” for spheres. 
for scattering by a sphere are summarized in Table 1. All manuscripts on this subject can be 
divided into two classes: those that fix x and vary N (or equivalently, the number of dipoles 
per sphere radius a/d) with y, and those that fix a/d and vary the size parameter with y. The 
former is easier to interpret; the latter is easier to simulate. To facilitate comparison between 
different methods we provide both x and a/d, however one of them is dependent on the other. 
Some additional information on these results follows below. 
Draine and Goodman [3] compared RR, DGF, and LDR for cross sections of a sphere 
with . DGF is generally more accurate than RR. For 16/ =da 1|1| ≤−m  LDR gives superior 
or comparable results to DGF, for i2+=m  LDR and DGF are comparable, and for 
 DGF is preferable over LDR. In the review of LDR DDA, Draine and Flatau [4] 
summarized that for  cross sections can be evaluated to accuracies of a few percent 
provided . In that case differential cross sections have satisfactory accuracy: relative 
errors up to 20-30%, but only where the absolute value of the differential cross sections is 
small. For spheres, such results are obtained even for 
i34 +=m
2|| ≤m
3|| ≤m . Comparison of CLDR to LDR 
[43] only results in minor differences. Generally CLDR results in slightly better accuracy for 
Csca, but worse for Cabs. 
Piller and Martin [44] compared FCD to LAK by studying the dependence of the mean 
relative error of the far-field electric fields (Ψ) on y for spheres with π=x , 2π and . It 
was shown that FCD (with a Hanning window filter for the electric permittivity ε) is roughly 
3 times more accurate than LAK in the range 
5.1=m
5.27.0 ≤≤ y  and gives similar accuracies for 
 (for larger spheres). Comparison of WD to traditional methods [13] was performed 
for spheres with 
4.0≤y
π=x , 2π and 32.1=m , i7.01.2 + . LAK was used to determine 
polarizabilities. For  in the range 32.1=m 3.14.0 ≤≤ y  overall accuracy was only slightly 
improved, but error peaks for certain values of y were smoothed out. For i7.01.2 +=m  
accuracy was improved 4-5 times over the whole range 3.1≤y . Piller also showed [10] that a 
combination of WD and FCD gives even better results. Generally FCD decreases the negative 
effects of Re(ε) on accuracy and WD those of Im(ε). 
Rahmani et al. [48] showed that RCB was clearly superior to CM in calculating cross 
sections for fixed  and m from 16/ =da i4.08.1 +  to i4.94.7 +  in the range . Two 
corrections (LDR and RR) over the static case were compared, and they gave similar overall 
results. Improvement of overall accuracy compared with CM was 2-5 times in all cases 
studied. For a thin slab, it was shown [46,48] that the internal fields calculated using RCB 
differ from those by CM mostly near the interfaces, where RCB yields much smaller errors, 
almost the same as far from interfaces. 
Collinge and Draine [47] compared LDR, RCB, and SCLDR in calculations of cross 
sections of spheres with . It was shown that for 12/ =da i01.033.1 +=m , LDR and SCLDR 
are superior in the range , while for 8.0≤y i45 +=m , SCLDR and RCB are superior. 
Convergence of cross sections for spheres and ellipsoids for increasing N with fixed x and 
different m (from  to i01.033.1 + i45+ ) also was studied. SCLDR showed the most stable 
results for all cases, being the most or close to the most accurate one; however, for ellipsoids 
with large Im(m) RCB gave significantly more accurate results for Csca, especially for larger y. 
Performance of the DDA for more complex shapes also was studied by different 
authors. Flatau et al. [50] compared DDA simulations for a bisphere with an exact solution 
from a multipole expansion. For i01.033.1 +=m , 16/ =da , and 8.0≤y , LDR was several 
times more accurate than DGF and resulted in errors of less than 0.5% for both Csca and Cabs. 
Xu and Gustafson [51] made a similar but much more extended study of LDR. For 
, i008.06.1 +=m 25/ =da , and , errors in C4.0≤y ext, Cabs, and θcos  are within 10%. For 
, errors in the angular dependence of S81.0=y 11 are up to 20% while S12 and S21 were 
completely wrong. For , errors in cross sections exceed 10% for . i02.05.2 +=m 3.0≥y
Errors in the angular dependencies of the Mueller matrix elements are within 10-20% for 
 and increase rapidly with increasing y. For a fixed 3.0=y 3=x  and , errors 
i004.06.1 +=m
ext, Cabs, and >< θcos  decrease from 10% to 1% while y decrease from 1 to 0.2. For 
, the angular dependence of S33.0=y 11 is in good agreement with the rigorous solution, 
while S12 and S21 differ significantly for certain orientations of the bisphere. 
Hage and Greenberg [14] compared LAK to experimental results obtained from 
microwave experiments on porous cubes. Using i005.0362.1 +=m , 64.0=y  and , 
they obtained a difference of less than 40% with the experimental results of angular scattering 
patterns, except for deep minima. Light scattering of cubes, tiles, and cylinders with similar 
parameters also was studied and comparable differences between experiment and theory were 
obtained. Theoretical errors were estimated to be less than 10%, except for deep minima. 
5504=N
Iskander et al. [34] conducted a limited test of LAK for small elongated spheroids, 
comparing the results to those obtained using an iterative extended boundary condition 
method. Using , calculations were performed for aspect ratios up to 20 with maximum 
size parameter of the long axis being 10 and 0.5 for 
i01.033.1 +=m  and  
respectively. Errors in scattering cross section were 21% and 11%, respectively. Ku [52] 
compared LAK with CM and the a
i28.076.1 +
1-term for different shapes, but his conclusions are based 
on a large parameter y (up to 2), and are therefore suspicious and not further discussed here. 
Andersen et al. [53] studied the performance of the DDA for Rayleigh-sized clusters of 
a few spheres (most DDA formulations are then equivalent to CM). Several constituent 
materials were tested, all with high refractive indices in the studied region. It was shown that 
the DDA failed to converge using the fixed computational resources for very high (up to 13.0) 
and very low (down to 0.12) Re(m); up to 30 dipoles were used per diameter of a single 
sphere. 
It can be concluded that particles with more complex shapes than spheres are more 
difficult to model with the DDA, leading to larger errors for the same m and y. This effect can 
be explained in general by the increase of surface to volume ratio and hence larger fraction of 
boundary subvolumes [32]. Another possible reason is complex regions, e.g. contact between 
two particles in a cluster, where rapid variation of the electric field deteriorates the overall 
accuracy. There is, however, a notable exception from this general tendency. Shapes, which 
can be modeled exactly by a set of cubical dipoles, e.g. a cube, can be simulated using the 
DDA much more accurately than spheres, especially for small y [32]. 
Draine and Flatau [4] have introduced a “rule of thumb” for discretization: use 10 
dipoles per wavelength in the medium (i.e. either y or yRe equal to 0.63, depending on the 
interpretation). Though it is widely used, the accuracy of the results, when using such 
discretization, is hard to deduce a priori. Draine and Flatau themselves derived an estimate of 
the error based on a set of test simulations. This estimate is described above and mentioned in 
Table 1; it is usually cited as a “few percent accuracy in cross sections.” However, it may 
significantly over- or under-estimate the error, especially for large size parameters. Moreover, 
it does not completely account for the dependence on m, even in the stated range of its 
application ( ), since DDA accuracy deteriorates rapidly with increasing m (see 2|| ≤m Table 
1). Still, the rule of thumb is good first guess for many applications. 
Most studies of DDA accuracy are limited to integral scattering quantities and, at most, 
the angular dependence of S11. In only a few manuscripts are other scattering quantities 
studied. For instance, Singham [54] simulated the angular dependence of Mueller matrix 
element S34 for spheres and less compact particles, using CM polarizability. It was shown that 
an accurate simulation of this element requires smaller values of y than for S11. For 55.1=x  
and  a calculation of S33.1=m 11 was accurate already for 8.0=y , while  was 
required for S
2.0≤y
34. It was also reported that for less compact objects like discs and rods, the 
required y was larger, 0.4 and 0.55 respectively, because of the smaller interaction between 
the dipoles. However, Hoekstra and Sloot argued [55] that this effect is mostly caused by the 
pronounced S34 sensitivity to surface roughness, which is significant for smaller size if y is 
fixed. They showed that for  and 7.10=x 05.1=m , very high accuracy is achieved with 
 because of the larger number of dipoles used. 66.0=y
Internal fields are an intermediate result in the DDA. They cannot be directly compared 
to the experimental results; however, all measured scattering quantities are derived from 
them. Therefore, a study of their accuracy can reveal greater understanding of the nature of 
DDA errors. Hoekstra et al. [56] performed such a study for LAK polarizability. Three 
spheres were examined with , 9, 5 and 9=x 05.1=m , i01.033.1 + ,  respectively. 
Values of y were 0.44, 0.42, and 0.51 respectively. The most significant errors in the 
amplitude of the internal field were localized at the boundary of the spheres with maximum 
relative errors of 3.4%, 19%, and 120% respectively. Errors in S
i4.15.2 +
12, S33, S34 were significant 
only for the third sphere. It was shown that for a given yRe these errors rapidly increase with m 
but only slightly depend upon x in the range from 1 to 10. Moreover, the DDA is capable of 
reproducing resonances of Mie theory, although their positions are slightly shifted (less than 
1% in m). 
Druger and Bronk [57] studied the accuracy of the internal fields for single and coated 
spheres. They used 5.1=x , , and CM polarizability. Errors in the internal fields were 
localized at the interfaces, with average errors larger than 30% for a single sphere with 
 and , and less than 7% for a single and concentric sphere with  and 
. The core of the concentric sphere has 
8.1≤m
8.1=m 17.0=y 3.1=m
08.0=y 1.1=m  and its diameter is half the total 
diameter. The angular dependence of the absolute values of S1 and S2 had significant errors in 
the side- and backscattering. It can be concluded that shape errors contribute mostly to the 
internal fields near the boundary, and increase with m. 
All the literature discussing DDA accuracy shows errors as a function of input 
parameters and discretization, which is the most straightforward way. The only exception so 
far is the rule of thumb, which is too general and approximate to be applied in many particular 
cases. A more useful way to present errors is to fix the desired accuracy for certain input 
parameters and find the discretization that results in such accuracy. Such an analysis can be 
applied directly to practical calculations and can be used to derive rigorous estimates of DDA 
computational requirements [58]. 
In a number of manuscripts the origin of errors in the DDA was examined to try to 
separate and compare shape and discretization errors [49,59-62]; however, no definite 
conclusions were reached. The uncertainty was due to the indirect methods used that have 
inherent interpretation problems. Recently, Yurkin et al. [63] proposed a direct method to 
separate shape and discretization errors, which can be used to study their fundamental 
properties. This method also can be applied to study the performance of different formulations 
aimed at decreasing shape errors, e.g. WD. For example, it has been shown that the maximum 
errors of S11(θ) for a sphere with  and 5=x 5.1=m , discretized using 16 dipoles per diameter 
( ), are mostly due to shape errors. However the same is not true for all measured 
quantities. In another manuscript [32] it was suggested that the discretization error should 
decrease more rapidly with decreasing y than shape errors. However, it is still hard to deduce 
a priori the importance of shape errors for a certain scatterer and y; hence, further systematic 
quantitative study is required. 
93.0=y
3.3 The DDA for clusters of spheres 
There are two main peculiarities when the DDA is applied to clusters of spheres. First, such 
particles are generally less compact, yielding smaller interactions between dipoles. This leads 
to a smaller condition number of the DDA interaction matrix and hence faster convergence of 
the iterative solver (see Section 4.1). Second, when the constituent spheres are small 
compared to the wavelength, each sphere can be modeled as one spherical subvolume, 
yielding some theoretical simplifications. 
A general theory exists [64] based on the Mie theory (generalized multiparticle Mie 
solution (GMM) [65]) that allows for highly accurate simulations of clusters of spheres. 
However, when many small spheres are used one wants to minimize the number of unknowns 
in the linear system. Direct reduction of the GMM to the lowest order (using only the first 
order expansion coefficients) leads to DDA + CM [64]. Improving accuracy in the GMM is 
done by accounting for higher multipole moments, while the DDA introduces higher order 
corrections to the coefficients of the linear system. It is not clear how the accuracy of these 
two methods compare with each other; however, the former should lead to a formulation 
similar to a coupled multipole method (Subsection 3.4) with a larger number of unknowns. 
DDA-based methods (starting usually with the integral equations introduced in Section 2) 
should be successful in making the formulation more accurate without increasing the number 
of unknowns, which is the goal for large clusters of small spheres. Moreover, the DDA may 
employ fast algorithms for solving the linear system. In this setting, the fast multipole method 
(FMM) (see Subsection 4.5) seems most promising. 
It should be noted, however, that a cluster having a small size parameter (i.e. in the 
electrostatic approximation) does not imply that all expansion coefficients, except the first 
one, are negligible. This is because the size of the constituent particles is also very small and 
the fields inside them are far from constant, especially when the spheres are located close to 
each other and have large refractive indices [66]. Therefore, the DDA does have some 
principal difficulties of calculating scattering by clusters of spheres. Mackowski [67], for 
instance, found that for some systems composed of spheres much smaller than the 
wavelength, up to 10 expansion terms were necessary to achieve convergence. In studies of 
osculating spheres, Ngo et al. [68] proved that the GMM could be chaotic and were able to 
calculate Lyapunov exponents, and that the slow convergence for the touching spheres was 
the result of the system lying in an attractor region. A recent paper by Markel et al. [69] 
presented computationally efficient modifications of the GMM in the static limit and 
demonstrated the insufficiency of the DDA to compute scattering properties of fractal 
aggregates accurately. However, Kim et al. [70] showed that the DDA is satisfactory in 
calculating the static polarizability of dielectric nanoclusters, especially of clusters with a 
large number of constituents. 
The development of DDA-based methods for calculating light scattering by clusters of 
small spheres was started by Jones [71,72], who developed a method similar to CM. Iskander 
et al. [34] used a method equivalent to LAK to calculate scattering of chained aerosol 
clusters. This subject was further investigated by Kosaza [73,74]. Lou and Charalampopoulos 
[75] (LC) further improved the calculations of the interaction term and scattering quantities. 
Starting from an integral equation for the internal field equivalent to Eq. (1), they assumed 
Eq. (11). After that the integrals in Eqs. (12) and (13) over spherical subvolumes can be 
evaluated analytically. The result for the interaction term is the following: 
),()()0( jiij ka rrGG η= , (70)
where a correction function η is defined as 
)O()101(1
cossin
3)( 42
x +−=
=η . (71)
Eq. (30) also is evaluated analytically, yielding 
∑ ⋅−−−=
ii knnkak )iexp()ˆˆ)((i)(
3)0( nrPInF η , (72)
( )∑ ∗⋅=
iikakC
inc)0(
ext Im)(4 EPηπ . (73)
The following expression for Cabs is stated without derivation: 
iikakC )Im()(4
abs EPηπ . (74)
Markel et al. [76] applied the DDA to fractal clusters of spheres, and studied their 
optical properties. However, they have not fixed the polarizability of a single dipole but rather 
treated it as a variable, calculating the dependence of a cluster’s optical characteristics upon it. 
Pustovit et al. [77] argued that the DDA is inaccurate for touching spheres. They developed a 
hybrid of the DDA and the GMM, which considers only pair interactions between spheres (as 
the DDA) but, when calculating them, accounts for higher multipole terms. This formulation 
can be considered as the one providing a more accurate evaluation of the interaction term 
(Eq. (13)), and hence similar to LC. 
LC was compared to DGF and LAK in a Csca computation of a cluster of 10 particles for 
 and . Differences between DGF and LAK are less than 1% (as 
expected), while the difference between LC and LAK increases quadratically with ka, 
reaching 10% for 
i7.07.1 +=m 5.005.0 ≤≤ ka
5.0=ka . However, as no exact (e.g. GMM) solution is presented, the 
accuracy of each individual method is not clear. 
Okamoto [42] tested the a1-term method for clusters of up to 3 touching spheres. No 
effective medium is needed in this case, making the method sounder. It was shown that the a1-
term is clearly superior to LDR in cross-sections calculations, when each sphere is treated as a 
single dipole. Errors of the a1-term are less than 10% for 2.1≤y  when . For 
three collinear touching spheres the errors are 30% and 40% for  and 2.8 when 
 and  respectively. However, errors do not seem to diminish 
significantly for small y (results are presented only down to 
i01.033.1 +=m
9.1≤y
i01.033.1 +=m i2 +
2.0=y ). Therefore, the a1-term 
seems suitable for obtaining quick crude estimations of cross sections. 
In the sequel of this subsection we mention several applications of the DDA to 
scattering from clusters of spheres. It was applied to describe the scattering by astrophysical 
dust aggregates [78,79] using the a1-term method. Hull et al. [80] applied CM DDA to Diesel 
soot particles. LC was applied [81] to the computation of light scattering by randomly 
branched chain aggregates. Lumme and Rahola [40] studied scattering properties of clusters 
of large spheres (each modeled by a set of dipoles) with the a1-term method considering 
astrophysical applications. Hage and Greenberg [35] studied scattering by porous particles, 
which were modeled as clusters of cubical cells making their method equivalent to standard 
LAK. Recently the DDA with LDR was used [82] to model scattering by porous dust grains 
and compare them to approximate theories, e.g. effective medium theories. It also was used to 
study light scattering by fractal aggregates [83], especially its dependence on the internal 
structure [84]. 
3.4 Modifications and extensions of the DDA 
Bourrely et al. [85] proposed to use small d to minimize surface roughness, but larger dipoles 
inside the particle. Starting with small dipoles with CM polarizability, one dipole is combined 
with 6 adjacent ones (if they all have the same polarizability) producing a dipole, located at 
the same point but with a 7 times larger polarizability. This operation is repeated while 
possible. Interaction terms are considered in their simplest form (Eq. (14)). This method 
allows the decrease of the shape errors with only a minor increase in the number of dipoles. 
The authors showed that this method is more than two times more accurate than CM for some 
test cases. 
Rouleau and Martin [86] proposed a generalized semi-analytical method. A dynamic 
grid is used to evaluate the integral in Eq. (1). First, a static grid is built inside the particle. 
Then each point on the static grid is used as an origin of a spherical coordinate system, and 
the particle is approximated by an ensemble of volume elements in these spherical 
coordinates. As usual, the polarization inside each subvolume is assumed constant, but 
Eq. (13) can be evaluated analytically in spherical coordinates. Polarization inside a 
subvolume is obtained by interpolation of its values at the points of the static grid. In addition, 
adaptive gridding is employed, where smaller subvolumes are used at the boundary of the 
particle. 
Mulholland et al. [87] proposed a coupled electric and magnetic dipole method 
(CEMD), where a magnetic dipole is considered at each subvolume together with an electric 
dipole. Polarizabilities are derived from the a1 and b1 terms of the Mie theory. CEMD 
requires two times more variables in the linear system, since the electric and magnetic fields 
are interconnected. Lemaire [88] went further and developed the coupled multipole method, 
considering also the electric quadrupole. Addition of the electric quadrupole can be 
considered as a more accurate evaluation of the interaction term in Eq. (13), as compared to 
Eq. (14). It results in even better accuracy than CEMD, but at the expense of additional 
computation time. The major disadvantage of all these four methods is that the matrix of the 
system of linear equations does not seem to have any special form, suitable for faster 
algorithms (see Section 4). Therefore computational costs are much larger compared to 
regular methods, thus limiting their practical use. In what follows, several DDA extensions 
are mentioned without further discussion. 
The theoretical basis for application of the DDA to optically anisotropic particles was 
summarized by Lakhtakia [89]. Loiko and Molochko [90] applied the DDA to study light 
scattering by liquid-crystal spherical droplets. Smith and Stokes [91] used the DDA to 
calculate the Faraday effect for nanoparticles. Researchers in the electrical engineering 
community applied MoM (in a variation that is equivalent to the DDA) to anisotropic 
scatterers [92,93]. 
Rectangular parallelepipeds can be used as subvolumes in the DDA [11,23,43]. This 
allows an accurate description of light scattering by particles with large aspect ratios, using 
fewer dipoles and is also compatible with FFT techniques (Subsection 4.4). 
Khlebtsov [94] proposed a simplification of the DDA, based on the assumption that all 
polarizations are parallel to the incident electric field. The number of variables is thus reduced 
three times, however at a cost of accuracy. Moreover, depolarization is completely ignored. 
Markel [95] analytically solved the DDA equations for scattering by an infinite one-
dimensional periodic dipole array. This approach is similar to the one used in obtaining the 
LDR formulation for dipole polarizability [3]. 
Chaumet et al. [96] generalized the DDA to periodic structures, and further to defects in 
a periodic grating on a surface [97]. The idea of using the complex Green’s tensor in the 
standard DDA formulation was summarized by Martin [98]. 
Yang et al. [99] used the DDA to calculate surface electromagnetic fields and determine 
Raman intensities for small metal particles of arbitrary shape.  
Lemaire and Bassrei [100] showed that the shape of an object can be reconstructed from 
the measured angle dependence of scattered intensities. This procedure can be thought of as 
an inversion of the dependence between dipole polarizabilities and scattering. This 
dependence is taken from the DDA. A similar idea is used in recent manuscripts on optical 
tomography [101-103]. 
Zubko et al. [104] modified the Green’s tensor used in the DDA to study the 
backscattering of debris particles. They showed that the far-field part of the Green’s tensor is 
responsible for both the backscattering brightness surge and the negative polarization branch. 
4 Numerical considerations 
In this section the numerical aspects of the DDA are discussed. One should keep in mind, 
however, that final simulation times depend not only on the chosen numerical methods but 
also on the particular implementation. Recently, Penttila et al. [105] have compared four 
different computer programs for the DDA. These are based on almost identical numerical 
methods: the Krylov-subspace iterative method (Section 4.1) combined with a FFT 
acceleration of the matrix-vector product (Section 4.4). However, simulation times may differ 
by several factors. Optimizations of computer codes are not further discussed in this review. 
4.1 Direct vs. iterative methods 
There are two general types of methods to solve linear systems of equations , where x 
is an unknown vector and A and y are known matrix and vector, respectively: direct and 
yAx =
iterative [106]. Direct methods give results in a fixed number of steps, while the number of 
iterations required in iterative methods is generally not known a priori. The most usual 
example of a direct method is LU decomposition, which allows quick solving for multiple y 
once the decomposition is performed. Iterative methods are usually faster, less memory 
consuming and numerically more stable. However, iterative methods cannot be considered 
superior over direct, since they strongly depend on the problem to solve [107]. 
For a general n×n matrix (in DDA Nn 3= ) computation time of LU decomposition is 
O(n3) and storage requirements O(n2), while computation time for one iteration is O(n2) [107]. 
Iterative methods for a general matrix converge in O(n) iterations, although some of them 
may not converge at all. However, in many cases satisfactory accuracy can be obtained after a 
much smaller number of iterations. In these cases, iterative methods can provide significant 
increases in speed, especially for large n. Most iterative methods access the matrix A only 
through matrix-vector multiplication (sometimes also with the transposed matrix), which 
allows the construction of special routines for calculation of these products. Such routines 
may decrease memory requirements, since it is no longer necessary to store the entire matrix, 
especially for matrices of special form (see Subsection 4.3). A special structure of the matrix 
may also allow acceleration of the matrix-vector product from O(n2) to O(nlnn) (see 
subsections 4.4, 4.5). However, the same applies to direct methods (see Subsection 4.3). 
Throughout DDA history, mostly iterative methods were employed (however see 
Subsection 4.6). At first, they were used to accelerate computations [1], but they also allowed 
larger numbers of dipoles to be simulated [6,108], since storage of the entire matrix is 
prohibitive for direct methods. The most widely used iterative methods in the DDA are 
Krylov-space methods, such as [107] conjugate gradient (CG), CG applied to the Normalized 
equation with minimization of Residual norm (CGNR), Bi-CG, Bi-CG stabilized (Bi-
CGSTAB), CG squared (CGS), generalized minimal residual (GMRES), quasi-minimal 
residual (QMR), transpose free QMR (TFQMR), and generalized product-type methods based 
on Bi-CG (GPBi-CG) [109]. 
An important part of the iterative solver is preconditioning, which effectively decreases 
the condition number of the matrix A and therefore speeds up convergence. However, this 
requires additional computational time during both initialization and each iteration. 
Preconditioning of the initial system can be summarized as [107] 
yMxMAMM 12
21 )( =
− , (75)
where M1 and M2 are left and right preconditioners, respectively. Preconditioners should 
either allow fast inversion or be integrated into the iteration process. The simplest 
preconditioner of the first type is the Jacobi (point), which is just the diagonal part of matrix 
A. An example of the second type of preconditioner is the Neumann polynomial 
preconditioner of order l: 
)( AIM . (76)
QMR and Bi-CG can be made to employ the complex symmetric (CS) property of the 
DDA interaction matrix to halve the number of matrix-vector multiplications [110] (and thus 
computational time). Lumme and Rahola [40] were the first to apply QMR(CS) to the DDA 
and compared it with CGNR. They used m from i1.06.1 +  to i43+ , and x from 1.3 to 13.5, 
corresponding to N from 136 to 20336. For all cases studied QMR(CS) was 2-4 times faster 
than CGNR. 
Rahola [9] further studied QMR(CS) and compared it to CGNR, Bi-CG(CS), Bi-
CGSTAB, CGS, GMRES (full and with different memory length). For a “typical small 
problem” (parameters were not specified, unfortunately) the convergence of different methods 
was tested and QMR(CS) along with Bi-CG(CS) showed the best results. Although full 
GMRES was able to converge in fewer iterations, GMRES with as much as 40 memory 
lengths was slower than QMR(CS). 
Flatau [111] reviewed the use of iterative algorithms in the DDA and tested many of 
them, together with several preconditioners. He calculated scattering of a homogenous sphere 
with  and m from 1.33 up to 1.0=x i0001.05+ , 1=x  and m from 1.33 up to  and 
. Left (L) and right (R) Jacobi-, and first-order Neumann polynomial 
preconditioners were tested. Unfortunately the number of dipoles N was not specified, which 
hampers comparison with other studies. For small particles CG(L) was superior for all 
refractive indices studied. CG and CG(R) showed similar results, while CGNR(L) and Bi-
CGSTAB(L) were about 4 times slower. For 
i33.1 +
i0001.03+
1=x  Bi-CGSTAB(L) was superior while Bi-
CGSTAB,(R) and CGS,(L),(R) were slightly worse. TFQMR (both with and without Jacobi 
preconditioner) was 3-4 times slower. The first-order Neumann preconditioner showed 
unsatisfactory results. It was concluded that Bi-CGSTAB(L) is the most satisfactory choice 
for the DDA, and that method is the default one used in the DDSCAT program [6]. 
Recently Fan et al. [112] have compared GMRES, QMR(CS), Bi-CGSTAB, GPBi-CG, 
and Bi-CG(CS). They tested them on wavelength-sized scatterers (x up to 10) with m up to 
, and concluded that GMRES with memory depth 30 was the fastest, although it 
required four times more memory than the other methods. However, only the times of the 
matrix-vector product was compared, while other parts of the iteration may also take 
significant time, especially for GMRES(30). Choosing from less memory-consuming 
methods, QMR(CS) and Bi-CG(CS) showed a better convergence rate than Bi-CGSTAB and 
GPBi-CG, especially when 
i2.05.4 +
2>m . Moreover, the authors pointed out some flaws in the 
comparison by Flatau [111], making his conclusions insufficient. 
Yurkin et al. [113] employed QMR(CS), Bi-CG(CS), and Bi-CGSTAB to simulate light 
scattering by spheres with x up to 160 and 40 for 05.1=m  and 2, respectively. It was shown 
that convergence of the iterative methods becomes very slow with increasing x and m (up to 
105 iterations are required), and none of them is clearly preferable to the others. Moreover, 
there seems to be no systematic dependency of the choice of the best iterative solver on x and 
m; however, the difference in computational time was less than a factor of two, except for the 
largest x and m studied. 
Rahola [114] showed that the spectrum of the integral scattering operator for any 
homogenous scatterer is a line in the complex plane going from 1 to m2, except for a small 
amount of points, which corresponds to refractive indices that cause resonances for the 
specific shape. The spectrum of A is similar, since this matrix is obtained in the DDA by 
discretization of the integral operator (see also [9]). Assuming that the spectrum of A exactly 
lies on the specified line, it was shown that an estimate for the optimal reduction factor6 γ can 
be given as 
Eq. (77) is an approximation valid for small particle sizes, where no, or only few, resonances 
are present. However, in all cases the spectrum of A resembles the spectrum of the linear 
operator, which is defined by shape, size and refractive index of the scatterer. Therefore, the 
spectrum, and thus convergence, should not depend significantly on the discretization. This 
fact was confirmed empirically in other manuscripts [9,63]. 
Budko and Samokhin [115] generalized Rahola’s results to arbitrary inhomogeneous 
and anisotropic scatterers. They described a region in the complex plane that contains the 
whole spectrum of the integral scattering operator. This region depends only on the values of 
m inside the scatterer and does not depend on x. They showed that for purely real m or for m 
with very small imaginary part this region may come close to the origin, therefore the 
spectrum may contain very small eigenvalues for particles larger than the wavelength. This 
                                                 
6 Norm of the residual is decreased by this factor every iteration. 
may explain the extremely slow convergence of the iterative solver for real m and large x, 
which was recently obtained in numerical simulations [113]. Based on the analysis of the 
spectrum of the integral scattering operator for particles much smaller than the wavelength, 
Budko et al. [116] proposed an efficient iteration method for this particular case. 
It can be concluded that there are several modern iterative methods (QMR(CS), Bi-
CG(CS), and Bi-CGSTAB) that have proved to be efficient when applied to the DDA. 
However, none of them can be claimed superior to the others, and one should test them for 
particular light-scattering problems. Moreover, except for the simplest cases, preconditioning 
of the DDA interaction matrix is almost not studied, while there is a need for it for large x and 
m, since then all methods converge extremely slowly or even diverge. It seems to us that the 
next major numerical advance in the DDA will be achieved by developing a powerful 
preconditioner for the DDA matrix. 
A large number of dipoles requires large computational power and, hence, parallel 
computers are commonly used, e.g. [108,113]. Parallel efficiency is not discussed here, but 
for iterative solvers, it is generally close to 1 [117]. However, this is not true for all 
preconditioners [107], and hence heavy preconditioners requiring large computational time in 
combination with a parallel DDA implementation should be employed with caution. 
4.2 Scattering order formulation 
The Rayleigh-Debye-Gans (RDG) approximation [27] consists in considering E(r) equal to 
Einc(r). F(n) is then obtained directly from Eq. (24). Generalization of the RDG approach is 
obtained by iteratively solving the integral equation (1), which can be rewritten as 
)()()( inc rΛErErE += , (78)
where Λ is a linear integral operator describing the scatterer. The iterative scheme is readily 
obtained by inserting the current (l-th) iteration of the electric field E(l)(r) into the right side of 
Eq. (78) and calculating the next iteration in the left side: 
)()()( )(inc)1( rΛErErE ll +=+ . (79)
The starting value is taken the same as in RDG, , and the general formula for 
the solution is the following: 
)()( inc)0( rErE =
inc )()(
l rEΛrE , (80)
which is a direct implementation of the well-known Neumann series: 
lΛΛI , (81)
where I is the unitary operator. A necessary and sufficient condition for Neumann-series 
convergence is 
1<Λ . (82)
Physical sense of this iterative method lies in successive calculations of interaction 
between different parts of the scatterer. The zeroth approximation (or RDG) accounts for no 
interaction; the first approximation considers the influence of scattering of each dipole on the 
others once, and so on. Eq. (82) states that the interaction inside the scatterer should be small, 
but not as small as required for the applicability of RDG ( 1<<Λ ). In scattering problems, 
especially in quantum physics, Eq. (80) is called the Born expansion. 
Although theoretically clear, the Born expansion is not directly applicable [118], since 
each successive iteration requires analytical evaluation of multidimensional integrals with 
rising complexity, which quickly becomes unfeasible even for the simplest scatterers. The 
latest result is probably that of Acquista [118], who evaluated the Born expansion for a 
homogenous sphere up to second order. Therefore, realistic application of the Born expansion 
does require discretization of the integral operator, which is naturally done in the DDA. 
A scattering order formulation (SOF) of the DDA was developed independently by 
Chiappetta [119] and Singham and Bohren [12,120] by applying the Neumann series to 
Eq. (17). Λ is then a matrix defined as jijij αGΛ = , where each element is a dyadic, which 
can be expressed as a 3×3 matrix. An explicit check of Eq. (82) for a certain scatterer is not 
feasible numerically, however de Hoop [121] derived a sufficient condition for scalar waves: 
1)(max)(2 20 <r
χπ kR , (83)
where R0 is the radius of the smallest sphere circumscribing the scatterer. Although not 
directly applicable to light scattering, Eq. (83) can be used as an estimate. 
The range of size parameter and refractive index where SOF converges is limited [120]. 
Moreover, even when SOF converges, more advanced iterative methods converge faster (see 
Subsection 4.1). However, SOF has clear physical sense and can be used to study the 
importance of multiple scattering. 
4.3 Block-Toeplitz 
A square matrix A is called Toeplitz if jiij aA −= , i.e. matrix elements on any line parallel to 
the main diagonal are the same [106]. In a block-Toeplitz (BT) matrix (of order K) elements 
ai are not numbers, but square matrices themselves: 
. (84)
A 2-level BT matrix has BT matrices as components ai. Proceeding recursively a multilevel 
BT (MBT) matrix for any number of levels is defined. 
Let us consider a rectangular lattice nx×ny×nz, numbered in the following way 
zzyyzxxzy ininninnni +−+−= )1()1( , (85)
where  indicates the position of the element along the axes. Let us also define 
the vector index . Then one can verify that the interaction matrix in Eq. 
},...,1{ μμ ni ∈
),,( zyx iii=i (20), 
defined by Eq. (13), satisfies the following: 
jiGGG −′== jiij . (86)
This equation alone can be used to greatly reduce the storage requirements of iterative 
methods by use of indirect addressing. Further improvement is to note that Eq. (86) defines a 
symmetric 3-level BT matrix (orders of subsequent levels – nx, ny, nz) whose smallest blocks 
are 3×3 matrices (dyadics) ijG . 
A rectangular lattice is not much of a restriction, since any scatterer can be embedded in 
an appropriate rectangular grid. However, additional “empty” dipoles should be introduced to 
build up the grid up to the full parallelepiped. Moreover, position and size of the dipoles 
cannot be chosen arbitrarily to better describe the shape of the scatterer. This is especially 
problematic for highly porous particles or clusters of particles, where the monomer has a size 
comparable to a single dipole. For all other cases these restrictions are minor compared to the 
large increase in computational speed, imposed by the BT-structure of the interaction matrix. 
A matrix-vector multiplication can be transformed to a convolution, which is computed using 
a fast Fourier transform (FFT) technique in O(nln(n)) operations (see Subsection 4.4). Note 
however, that alternative techniques exist that do not require a regular grid (see Subsection 
4.5). 
The BT-structure also permits acceleration of direct methods. Flatau et al. [122] used an 
algorithm for inversion of symmetric BT-matrices. It has complexity )(O 3 xnn  and storage 
requirements )(O 2 xnn , since only 2 block columns of the inverse matrix need to be stored. 
In this case the x-axis is oriented along the longest particle dimension. Recently Flatau [123] 
studied the special case of 1D DDA where all dipoles are located on a straight line and 
equally spaced, in which systems of equations for different components can be separated. The 
interaction matrix for each component is symmetric Toeplitz, and a modern fast algorithm can 
be applied for its inversion. This method requires preliminarily solving linear equations for 
two right sides (e.g., by some iterative technique); then multiplication of the inverse matrix by 
any vector (i.e., a solution of the linear system for any right part) requires only O(nln(n)) 
operations. However, Flatau pointed out a strict limitation for all methods for fast calculation 
of the inverse of the interaction matrix: they are applicable only when polarizabilities of all 
dipoles are the same, since otherwise the first term on the right side of Eq. (20) ruins the BT 
structure on the diagonal of the interaction matrix. Therefore, they are currently limited to 
homogenous rectangular scatterers. Fortunately, it is not a problem for matrix-vector 
multiplication, since the diagonal term can be evaluated independently and added to the final 
result. 
4.4 FFT 
Goodman et al. [124] showed that multiplication of the interaction matrix for a rectangular 
lattice (see Subsection 4.3) by a vector can be transformed into a discrete convolution 
jii PGPGPGy ′′=′== ∑∑∑
)2,2,2(
)1,1,1(
)1,1,1(1
zyxzyx nnnnnnN
jij , (87)
where iG′  is defined by Eq. (86) (and 0=′0G ) for μμ ni ≤ and 
⎧ ≤≤∀
otherwise,0
1:, μμμ njj
P . (88)
Both G ′  and  are then regarded as periodic in each dimension μ with period 2nP′ μ. A 
discrete convolution can be transformed with a FFT to an element-wise product of two 
vectors, which is easily computed. It requires evaluation of a direct and inverse FFT for each 
matrix-vector product. Each of them is a 3D FFT of order 2nx×2ny×2nz. This operation is done 
for each of the 3 Cartesian components of P′  and preliminary calculations is performed for 6 
independent tensor components of G ′ . 
A slightly different method can be devised based on the paper by Barrowes et al. [125], 
who developed an algorithm for multiplication of any MBT by a vector. The multiplication is 
brought down to a 1D convolution that is evaluated by two 1D FFTs of order 
)12)(12)(12( −−− zyx nnn . Flatau [123] proposed an algorithm of matrix-vector 
multiplication for BT interaction matrix (e.g. 1D DDA), which requires twice as many FFTs 
as the standard algorithm, but of order n instead of 2n. Although Flatau stated that an 
extension of this algorithm to the general 3D case is straightforward, it is at least not trivial 
and probably its complexity will scale the same as standard methods. 
4.5 Fast multipole method 
The fast multipole method (FMM) was developed by Greengard and Rokhlin [126] for 
efficient evaluation of the potential and force fields in N-body simulations where all pairwise 
interactions of N particles are computed. The FMM is based on truncated potential expansions 
[127]. It is also called a hierarchical tree method because particles are grouped together in a 
hierarchical way, and the interaction between single particles and this hierarchy of particle 
groups is calculated [128]. However, some researchers distinguish between single- and 
multilevel FMM [129,130]; only the latter is truly hierarchical. The FMM naturally fits the 
DDA, since the matrix-vector multiplication is actually computing the total field on each 
single dipole due to all other dipoles, as was noted by Hoekstra and Sloot [128]. The 
computational complexity of the FMM (see below) is similar to FFT-based methods (see 
Subsection 4.4), but it does not require any regularity of the grid, thus making it applicable to 
any scatterer. The drawback is that the FMM is conceptually more complex, making it much 
harder to code. Nonetheless, the FMM was implemented in the DDA by Rahola [9,127]. 
Error analysis is critical for the FMM, since the acceleration is obtained by using 
approximations, in contrast to exact FFT-based methods. Approximation parameters are 
chosen to keep an error, calculated according to some estimate, in certain bounds. The more 
exact the error estimate is, the less computations are required; thus, the faster the whole 
algorithm. Therefore, algorithm complexity is directly connected to error analysis [131]. 
Koc and Chew [129] described the application of multilevel FMM to the DDA. They 
used semi-empirical formulae to determine the number of terms in multipole series, and 
obtained O(N) complexity. However, rigorous, close to exact, error analysis is still lacking for 
the FMM applied to the DDA. It will allow obtaining a real algorithm complexity with 
guaranteed accuracy. Such an analysis has been conducted for 2D acoustic scattering [130], 
and for light scattering formulated in terms of surface integrals [131]. In both cases the FMM 
was proven to have an asymptotic complexity O(Nln2(N)). Application of the FMM to 
surface-integrals formulation of light scattering was reviewed by Dembart and Yip [132]. 
Another problem of implementing the FMM is that it is completely dependent upon the 
exact form of the interaction potential ijG . All manuscripts mentioned above deal with 
interaction between point dipoles, i.e. Eq. (14). If a more complex expression for ijG  is used 
(e.g. IT), most of the FMM should be developed anew. This makes integration of the FMM 
and the DDA a formidable problem. 
The FMM is a promising method to calculate light scattering by particles that cannot be 
mapped effectively on a rectangular grid; however, there is still space for improving its theory 
to make it more robust and guarantee certain accuracy. 
The FMM is not the only hierarchical tree method available. For instance, a very 
intuitively simple method was proposed by Barnes and Hut [133,134]. Multipole expansions 
over the center of mass in gravitational computations are used, contrary to geometrical center 
in the FMM. It automatically eliminates the second term in the multipole expansion, and 
allows fast evaluation of monopole terms. Though this method is much simpler and clearer 
than the FMM, it has very little control over the errors that can be studied almost exclusively 
empirically. It can be applied to the DDA without significant increase in the total 
computational errors.7
An alternative approach was proposed by Ding and Tsang [135]. They studied scattering 
from trees and used a sparse matrix iterative approach. The interaction matrix is divided into a 
strong part, which accounts for interaction between nearby dipoles, and a complement weak 
part: . The strong part is sparse and therefore allows quick solution of the linear 
system. The weak part is a small correction that is accounted for iteratively: 
ws AAA +=
yxA =)0(s , . )(w)1(s ll xAyxA −=+ (89)
The authors demonstrate potential of this approach for some test cases. 
4.6 Orientation averaging and repeated calculations 
In many physical applications, one is interested in optical properties of an ensemble of 
randomly oriented particles. When the concentration of particles is small, multiple scattering 
is negligible and the optical properties are obtained by averaging single-particle scattering 
over different particle orientations. More general problems, where particles are not identical 
or multiple scattering is significant, are not considered here. 
Orientation averaging of any scattering property can be described as the integral over 
the Euler’s orientation angles (including a probability distribution function if necessary), 
                                                 
7 Hoekstra AG, unpublished results 
which is brought down to a sum by appropriate quadrature. The problem therefore consists in 
calculation of some scattering property for a set of different orientations of the same particle. 
The easiest way is to calculate it by solving sequentially and independently each problem 
from the set. However, the large size of this set calls for some means of reducing the 
calculations. This is especially relevant when the particle is asymmetric; hence, its optical 
properties are sensitive to particle orientation. Let us further assume for clarity that we are 
interested in the scattering matrix at a certain scattering angle. All the discussion for other 
scattering properties is analogous or even simpler. 
Singham et al. [136] noted that the set of problems described above is physically 
equivalent to a fixed orientation of the particle and different incident and scattering directions. 
The latter are determined by transformation of the laboratory reference frame to the reference 
frame associated with the particle. The amplitude scattering matrix, and hence the Mueller 
matrix, also is transformed along with the reference frame (see e.g. [137] for transformation 
formulae). There are two immediate advantages of using a fixed particle orientation. First, A 
is kept constant (see though discussion below), and therefore the construction of A is done 
only once. Second, the amplitude matrix for any scattering angle is quickly obtained after the 
linear system is solved (for two incident polarizations). Hence, integration over one Euler 
angle is relatively fast. 
The constancy of A can be exploited to further reduce the time of orientation-averaging. 
If  or its LU decomposition is obtained [75,136], a single solution for any right part y can 
be obtained in n
2 operations – the same or less time than required for one iteration using 
general iterative methods (see Subsection 4.1). Moreover, Singham et al. [136] and McClain 
and Ghoul [138] independently proposed an analytical way of averaging the scattering matrix 
at any scattering angle, which requires O(n2) operations once  is known. Khlebtsov [139] 
extended this technique to averaging of extinction and absorption cross sections. 
However, by employing special properties of the matrix A in the DDA allows 
computing matrix-vector products in O(nln(n)) operations (see subsections 4.4, 4.5). Although 
some acceleration of direct methods also can be performed (see Subsection 4.3), they are still 
O(n2) or slower. For large n, iterative methods (assuming that they converge in much less than 
n iterations) are clearly preferable, even if many quadrature points are used. Moreover, large n 
is unattainable by direct methods because of storage requirements. Another improvement 
could be using a heavy preconditioner, which has large initialization cost and greatly 
increases convergence rate. Initialization cost is then justified because it is computed only 
once. Possible candidates are incomplete factorization preconditioners [107]. 
Above it was stated that A is constant for a fixed-orientation particle. However, modern 
DDA formulations (e.g. LDR) take into account the direction of light incidence. Hence A 
depends upon this direction, but only weakly through ( )( )2O kd  corrections. This complicates 
the techniques described above, however probably they still may be used together with some 
special methods to correct for small changes in A on every step. Such methods have not been 
developed as yet. 
Another possibility to perform orientation averaging is to first compute the T-matrix of 
the particle, which then allows analytical averaging [140]. The T-matrix formalism is based 
on the multipole expansion which is truncated at some order N0. Although N0 is hard to 
deduce a priori, usually it is several times x [141,142]. The number of rows in the T-matrix 
equals . The simplest way to evaluate the T-matrix based on the DDA is to solve 
for every incident spherical wave (i.e. for each row of the T-matrix) independently [141]. 
Then the above discussion about optimizing this repeated calculation is relevant. Using 
iterative techniques with N
)2(2 00 +NN
iter number of iterations, computation time is 
( ) ( )[ ]20iter20 O)ln(O NNNNN + , where the first term in the sum is the time for solving the linear 
system, and the second one is the actual computation of the values in the row of the T-matrix. 
A new method to obtain the T-matrix from the DDA interaction matrix was proposed by 
Mackowski [141]. This requires two summations with computational time ( ))ln(O 20 NNN  and 
( )NN 40O . Mackowski showed that for 5=x  his method is an order of magnitude faster than 
the straightforward one. 
Recently Muinonen and Zubko [143] have proposed a way to optimize ensemble 
averaging of DDA results over different sizes and refractive indices. It is based on calculating 
a “good guess” for the initial vector in the iterative solver using results of the calculations 
with similar parameters. Similar ideas can be used to optimize simulation of a set of slightly 
different shapes or orientational averaging. 
Use of repeated calculations to increase the accuracy of DDA simulations was proposed 
recently by Yurkin et al. [63]. Several independent simulations with different discretization 
parameter were performed and results were extrapolated to the infinite discretization giving 
better accuracy than those of a single DDA simulation. 
5 Comparison of the DDA to other methods 
Hovenier et al. [144] compared the DDA, the extended boundary condition method (EBCM), 
and the separation of variables method (SVM) for calculations of scattering by spheroids, 
finite cylinders and bispheres. Parameters of the problems were as follows: , 
equivolume size parameter , 
i01.05.1 +=m
5=x 6.0=y . The angular dependencies of scattering matrix 
elements were calculated. The EBCM and SVM seemed to achieve an exact solution, and the 
DDA showed little errors, except for backscattering angles, where they were up to 10-20%. 
Wriedt and Comberg [145] compared the DDA, EBCM, and finite difference time 
domain (FDTD) method for a cube with 33.1=m , 1.5 and 9.2=x , 4.9, 9.7. For  and 
4.9. The DDA and EBCM achieved good accuracy in calculation of scattering intensity angle 
dependence; the DDA was 2-5 times faster, but consumed 8-16 times more memory (y was in 
the range 0.3-0.5). The FDTD had similar computational requirements as the DDA but was 
less accurate. For  the DDA was the only one to achieve little errors within the given 
computational resources. 
9.2=x
7.9=x
Comberg and Wriedt [146] compared the DDA, GMM (see Subsection 3.3) and the 
generalized multipole technique (GMT) for clusters of a few spheres. A single sphere had x in 
the range 4–20 and 33.1=m , 1.5. All the methods managed to achieve good accuracy, but the 
GMM was one order of magnitude (and for large x even several orders) faster than the other 
two. The DDA and GMT also were used to compute scattering by a cluster of two oblate 
spheroids with  and . The DDA was less accurate and consumed 4 times more 
memory, but was 6 times faster than the GMT. 
5=x 33.1=m
Wriedt et al. [147] compared the DDA, FDTD, GMT, and discrete sources method 
(DSM) for the calculation of light scattering by a red blood cell (RBC) with  and 
. Accuracy of all methods was similar. The DDA and GMT showed similar 
calculation times; they were 7 times faster than the FDTD and 12 times slower than the DSM. 
It should be noted that the latter employed the axisymmetric property of RBC. 
06.1=m
Recently Yurkin et al. [58] systematically compared the DDA and the FDTD for 
spheres with m from 1.02 to 2 and x from 10 to 100, depending on m. It was shown that 
numerical performance of the DDA is much more sensitive to the refractive index than that of 
the FDTD. Therefore, the DDA is preferable for small m, the FDTD for larger m. Cleary, the 
crossover point is not well defined and will depend on the details of the problem at hand as 
well as on the particular implementations of both methods. 
The main advantage of the DDA is that it is one of the most general methods, having a 
very broad range of applicability, limited only by available computational power. The reverse 
of this advantage is that it has almost no means to use the symmetry of the scatterer. Thus the 
DDA is not able to compete with the EBCM for homogenous axisymmetric scatterers. For 
homogenous non-axisymmetric scatterers the DDA is competitive with the EBCM for single-
particle orientation, but the latter allows much faster orientation averaging. The EBCM has 
little applicability to inhomogeneous scatterers, where the DDA can be applied without any 
changes. Comparison between the FDTD and the DDA suggests that the DDA is more 
suitable for small m. It also should be noted that the FDTD is even more general, being easily 
applicable to non-harmonic incident electric fields. Moreover, simulation of one pulse 
incident wave with the FDTD gives the solution for a complete spectrum of incident harmonic 
plane waves, but with a limitation on accuracy. 
6 Concluding remarks 
The DDA has been reviewed using a general framework based on the integral equation for the 
electric field. Although mainstream DDA algorithms as used in several production computer 
programs, has not changed significantly since 1994, many different improvements have been 
proposed since that time. Some of them do improve the accuracy or numerical performance of 
the DDA; however, they still wait for a wide acceptance. It seems that a critical mass of new 
improvements is building up, hopefully resulting in a next breakthrough in the field of the 
DDA. 
In our opinion, future major improvements in the DDA computer implementations will 
be connected with one of the following: 
1) Decreasing shape errors by implementing WD or similar techniques. 
2) Improving polarizability and interaction terms by techniques that are still to be 
developed similar to IT and PEL. 
3) Studying different preconditioners for the DDA interaction matrix, either trying some of 
the known types or developing one considering the special structure of the matrix. 
Item (1) should improve the overall accuracy of the DDA, especially for cases where shape 
errors are dominant, item (2) should expand the DDA applicability region to higher refractive 
indices, and item (3) should boost overall performance, especially for large size parameters 
and/or refractive indices. 
Acknowledgements 
We thank Dan Mackowski for clarifying discussion on the simulations of scattering by 
clusters of spheres and Gorden Videen for critically reading the manuscript and for valuable 
discussions. Our research is supported by Siberian Branch of the Russian Academy of 
Sciences through the grant 2006-03. 
Appendix. Description of used acronyms and symbols 
See Tables A1 and A2. 
Table A1. Acronyms in alphabetical order. 
Acronym Description Sectiona
(L) left Jacobi preconditioner 4.1
I right Jacobi preconditioner 4.1
a1-term (M) dipole term in the Mie theory 3.1
Bi-CGSTAB Bi-CG stabilized 4.1
BT block-Toeplitz 4.3
CEMD coupled electric and magnetic dipole 3.4
CG conjugate gradient 4.1
CGNR CG applied to normalized equation with minimization of residual norm 4.1
CGS CG squared 4.1
CS complex symmetric 4.1
CLDR corrected LDR 3.1
CM Clausius-Mossotti 2
DDA discrete dipole approximation 1
DGF digitized Green’s function 1
DSM discrete sources method 5
EBCM extended boundary condition method 5
FCD filtered coupled dipoles 3.1
FDTD finite difference time domain 5
FFT fast Fourier transform 4.4
FMM fast multipole method 4.5
GMM generalized multiparticle Mie solution 3.3
GMRES generalized minimal residual 4.1
GMT generalized multipole technique 5
GPBi-CG generalized product-type methods based on Bi-CG 4.1
IT integration of Green’s tensor 3.1
LAK Lakhtakia 3.1
LC Lou and Charalampopoulos 3.3
LDR lattice dispersion relation 3.1
MBT multilevel BT 4.3
MoM method of moments 3.1
PEL Peltoniemi 3.1
PP Purcell and Pennypacker 1
QMR quasi-minimal residual 4.1
RBC red blood cell 5
RCB Rahmani, Chaumet, and Bryant 3.1
RDG Rayleigh-Debye-Gans approximation 4.2
RR radiative reaction correction 3.1
SCLDR surface-corrected LDR 3.1
SOF scattering order formulation 4.2
SVM separation of variables method 5
TFQMR transpose free QMR 4.1
WD weighted discretization 3.1
a where it is explained or first appears (if no explanation is given). 
Table A2. Symbols used, Latin and Greek letters in alphabetical order.a
Symbols Description Section
(0) superscript: approximate value (usually under constant field assumption) 2
(n) superscript: after n-th iteration 4.2
>< θcos  asymmetry parameter 2
* superscript: complex conjugate 2
A a matrix 3.1
a kk  2
a radius of (equivalent) sphere 3.1
B correction matrix in SCLDR 3.1, Eq. (64)
b1 – b3 numerical coefficients in polarization prescriptions 3.1
C  tensor of electrostatic solution 3.1, Eq. (60)
Csca, Cabs, Cext scattering, absorption, extinction cross section 2
c speed of light in vacuum 2
d size of a cubical cell 2
E, Einc, Eexc, Eself, Esca (total) electric field, incident, exciting, self-induced, scattered 2
e0 polarization vector of the incident wave, 10 =e  2
e superscript: effective 3.1
F scattering amplitude 2
f a function; 
volume filling factor 
G  free space dyadic Green’s function (tensor) 2
2sG  G  in static limit 
ijG  interaction term 
H superscript: conjugate transpose 3.1
hr impulse response function of a filter 3.1
2, 4.1I , I identity dyadic (tensor), operator (matrix) 
i, j subscript: vector indices 4.3
i imaginary unity 2
i, j subscript: number of the dipole 2
K order of a BT matrix 4.3
k free space wave vector 2
L  self-term dyadic 2
M integral associated with finiteness of V0; 
preconditioner 
M  dyadic associated with M 2
m refractive index (relative) 3.1
N total number of dipoles 2
n rr  2
n size of a matrix 4.1
n′ˆ  external normal to the surface 2
nx, ny, nz sizes of the rectangular lattice 4.3
P polarization 2
p superscript: principal 3.1
q dπ2  3.1
R rr ′−  2
R0 radius of the smallest sphere circumscribing the scatterer 4.2
2r, r′ radius-vectors 
S LDR coefficient dependent on incident polarization 3.1, Eq. (49)
Si amplitude matrix element 3.2
Sij Mueller matrix element 3.2
s superscript: secondary; 
strong; 
subscript: equivalent spherical dipole 
T  boundary condition tensor 3.1, Eq. (66) 
t time 2
V volume of the scatterer 2
V0 exclusion volume 2
w superscript: weak 4.5
x unknown vector 4.1
Table A2 (continued) 
Symbols Description Section
x size parameter of scatterer 3.2
x, y, z Cartesian coordinates 4.3
y a known vector (right side of a linear system) 4.1
y kdm ||  3.2
yRe kdm)Re(  3.2
2α, α  polarizability, tensor 
4.1γ optimal reduction factor 
3.1δ Kronecker symbol 
2ε electric permittivity (relative) 
3.3η correction function 
Λ intermediate tensor in RCB method 3.1, Eq. (62)
4.2Λ linear integral operator, its matrix 
μ, ν, ρ, τ, … sub-, superscript: Cartesian components of vectors (tensors) 2
ξ, ψ Riccati-Bessel functions 3.1
2χ electric susceptibility 
3.2Ψ mean relative error of far-field electric field 
2Ω solid angle 
2ω circular frequency of the harmonic electric field 
a common sub- and superscripts are given on their own. For all vectors – the same symbol but in italic (instead of 
bold) denotes Euclidian norm of the vector (except unitary vectors). 
References 
 [1]  Purcell EM, Pennypacker CR. Scattering and adsorption of light by nonspherical dielectric grains. 
Astrophys J 1973;186:705-714. 
 [2]  Draine BT. The discrete-dipole approximation and its application to interstellar graphite grains. 
Astrophys J 1988;333:848-872. 
 [3]  Draine BT, Goodman JJ. Beyond clausius-mossotti - wave-propagation on a polarizable point lattice and 
the discrete dipole approximation. Astrophys J 1993;405:685-697. 
 [4]  Draine BT, Flatau PJ. Discrete-dipole approximation for scattering calculations. J Opt Soc Am A 
1994;11:1491-1499. 
 [5]  Draine BT. The discrete dipole approximation for light scattering by irregular targets. In: Mishchenko 
MI, Hovenier, JW, Travis, LD, editors. Light Scattering by Nonspherical Particles, Theory, 
Measurements, and Applications. New York: Academic Press, 2000. p. 131-145. 
 [6]  Draine BT, Flatau PJ. User guide for the discrete dipole approximation code DDSCAT 6.1. 
http://xxx.arxiv.org/abs/astro-ph/0409262, 2004. 
 [7]  Goedecke GH, O'Brien SG. Scattering by irregular inhomogeneous particles via the digitized Green's 
function algorithm. Appl Opt 1988;27:2431-2438. 
 [8]  Lakhtakia A. 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 C 1992;3:583-603. 
 [9]  Rahola J. Solution of dense systems of linear equations in the discrete-dipole approximation. SIAM J Sci 
Comp 1996;17:78-89. 
 [10]  Piller NB. Coupled-dipole approximation for high permittivity materials. Opt Comm 1999;160:10-14. 
 [11]  Chaumet PC, Sentenac A, Rahmani A. Coupled dipole method for scatterers with large permittivity. 
Phys Rev E 2004;70:036606. 
 [12]  Singham SB, Bohren CF. Light scattering by an arbitrary particle: a physical reformulation of the 
coupled dipole method. Opt Lett 1987;12:10-12. 
 [13]  Piller NB. Influence of the edge meshes on the accuracy of the coupled-dipole approximation. Opt Lett 
1997;22:1674-1676. 
 [14]  Hage JI, Greenberg JM, Wang RT. Scattering from arbitrarily shaped particles - theory and experiment. 
Appl Opt 1991;30:1141-1152. 
 [15]  Kahnert FM. Numerical methods in electromagnetic scattering theory. J Quant Spectrosc Radiat Transf 
2003;79:775-824. 
 [16]  Peterson AW, Ray SL, Mittra R. Computational Methods of Electromagnetic Scattering. IEEE Press, 
1998. 
 [17]  Kim OS, Meincke P, Breinbjerg O, Jorgensen E. Method of moments solution of volume integral 
equations using higher-order hierarchical Legendre basis functions. Radio Science 2004;39. 
 [18]  Lu CC. A fast algorithm based on volume integral equation for analysis of arbitrarily shaped dielectric 
radomes. IEEE Trans Ant Propag 2003;51:606-612. 
http://xxx.arxiv.org/abs/astro-ph/0409262,
 [19]  Ivakhnenko V, Eremin Y. Light scattering by needle-type and disk-type particles. J Quant Spectrosc 
Radiat Transf 2006;100:165-172. 
 [20]  Wriedt T. A review of elastic light scattering theories. Part Part Sys Charact 1998;15:67-74. 
 [21]  Chiappetta P, Torresani B. Some approximate methods for computing electromagnetic fields scattered by 
complex objects. Meas Sci Technol 1998;9:171-182. 
 [22]  Mishchenko MI, Travis LD, Lacis AA. Scattering, Absorption, and Emission of Light by Small Particles. 
Cambridge: Cambridge University Press, 2002. 
 [23]  Tsang L, Kong JA, Ding KH, Ao CO. Scattering of Electromagnetic Waves: Numerical Simulations. 
New York: Wiley, 2001. 
 [24]  Jones AR. Light scattering for particle characterization. Prog Ener Comb Sci 1999;25:1-53. 
 [25]  Yanghjian AD. Electric dyadic Green's function in the source region. IEEE Proc 1980;68:248-263. 
 [26]  Peltoniemi JI. Variational volume integral equation method for electromagnetic scattering by irregular 
grains. J Quant Spectrosc Radiat Transf 1996;55:637-647. 
 [27]  Bohren CF, Huffman DR. Absorption and scattering of Light by Small Particles. New York: Wiley, 
1983. 
 [28]  Draine BT, Weingartner JC. Radiative torques on interstellar grains .1. Superthermal spin-up. Astrophys 
J 1996;470:551-565. 
 [29]  Hoekstra AG, Frijlink M, Waters LBFM, Sloot PMA. Radiation forces in the discrete-dipole 
approximation. J Opt Soc Am A 2001;18:1944-1953. 
 [30]  Chaumet PC, Rahmani A, Sentenac A, Bryant GW. Efficient computation of optical forces with the 
coupled dipole method. Phys Rev E 2005;72:046708. 
 [31]  Hoekstra AG. Computer simulations of elastic light scattering. PhD thesis. University of Amsterdam, 
Amsterdam, 1994. 
 [32]  Yurkin MA, Maltsev VP, Hoekstra AG. Convergence of the discrete dipole approximation. I. Theoretical 
analysis. J Opt Soc Am A 2006;23:2578-2591. 
 [33]  Jackson JD. Classical Electrodynamics. New York: Wiley, 1975. 
 [34]  Iskander MF, Chen HY, Penner JE. Optical-scattering and absorption by branched chains of aerosols. 
Appl Opt 1989;28:3083-3091. 
 [35]  Hage JI, Greenberg JM. A model for the optical-properties of porous grains. Astrophys J 1990;361:251-
259. 
 [36]  Livesay DE, Chen KM. Electromagnetic fields induced inside arbitrarily shaped biological bodies. IEEE 
Trans Microw Theory Tech 1974;22:1273-1280. 
 [37]  Lakhtakia A, Mulholland GW. On 2 numerical techniques for light-scattering by dielectric agglomerated 
structures. J Res Nat Inst Stand Technol 1993;98:699-716. 
 [38]  Dungey CE, Bohren CF. Light-scattering by nonspherical particles - a refinement to the coupled-dipole 
method. J Opt Soc Am A 1991;8:81-87. 
 [39]  Doyle WT. Optical properties of a suspension of metal spheres. Phys Rev B 1989;39:9852-9858. 
 [40]  Lumme K, Rahola J. Light-scattering by porous dust particles in the discrete-dipole approximation. 
Astrophys J 1994;425:653-667. 
 [41]  van de Hulst HC. Light Scattering by Small Particles. New York: Dover, 1981. 
 [42]  Okamoto H. Light scattering by clusters: the a1-term method. Opt Rev 1995;2:407-412. 
 [43]  Gutkowicz-Krusin D, Draine BT. Propagation of electromagnetic waves on a rectangular lattice of 
polarizable points. http://xxx.arxiv.org/abs/astro-ph/0403082, 2004. 
 [44]  Piller NB, Martin OJF. Increasing the performance of the coupled-dipole approximation: A spectral 
approach. IEEE Trans Ant Propag 1998;46:1126-1137. 
 [45]  Gay-Balmaz P, Martin OJF. A library for computing the filtered and non-filtered 3D Green's tensor 
associated with infinite homogeneous space and surfaces. Comp Phys Comm 2002;144:111-120. 
 [46]  Rahmani A, Chaumet PC, Bryant GW. Coupled dipole method with an exact long-wavelength limit and 
improved accuracy at finite frequencies. Opt Lett 2002;27:2118-2120. 
 [47]  Collinge MJ, Draine BT. Discrete-dipole approximation with polarizabilities that account for both finite 
wavelength and target geometry. J Opt Soc Am A 2004;21:2023-2028. 
 [48]  Rahmani A, Chaumet PC, Bryant GW. On the importance of local-field corrections for polarizable 
particles on a finite lattice: Application to the discrete dipole approximation. Astrophys J 2004;607:873-
878. 
 [49]  Evans KF, Stephens GL. Microwave radiative-transfer through clouds composed of realistically shaped 
ice crystals .1. Single scattering properties. J Atmos Sci 1995;52:2041-2057. 
 [50]  Flatau PJ, Fuller KA, Mackowski DW. Scattering by 2 spheres in contact - comparisons between 
discrete-dipole approximation and modal-analysis. Appl Opt 1993;32:3302-3305. 
 [51]  Xu YL, Gustafson BAS. Comparison between multisphere light-scattering calculations: Rigorous 
solution and discrete-dipole approximation. Astrophys J 1999;513:894-909. 
 [52]  Ku JC. Comparisons of coupled-dipole solutions and dipole refractive-indexes for light-scattering and 
absorption by arbitrarily shaped or agglomerated particles. J Opt Soc Am A 1993;10:336-342. 
http://xxx.arxiv.org/abs/astro-ph/0403082,
 [53]  Andersen AC, Mutschke H, Posch T, Min M, Tamanai A. Infrared extinction by homogeneous particle 
aggregates of SiC, FeO and SiO2: Comparison of different theoretical approaches. J Quant Spectrosc 
Radiat Transf 2006;100:4-15. 
 [54]  Singham SB. Theoretical factors in modeling polarized light scattering by arbitrary particles. Appl Opt 
1989;28:5058-5064. 
 [55]  Hoekstra AG, Sloot PMA. Dipolar unit size in coupled-dipole calculations of the scattering matrix-
elements. Opt Lett 1993;18:1211-1213. 
 [56]  Hoekstra AG, Rahola J, Sloot PMA. Accuracy of internal fields in volume integral equation simulations 
of light scattering. Appl Opt 1998;37:8482-8497. 
 [57]  Druger SD, Bronk BV. Internal and scattered electric fields in the discrete dipole approximation. J Opt 
Soc Am B 1999;16:2239-2246. 
 [58]  Yurkin MA, Brock RS, Lu JQ, Hoekstra AG. Systematic comparison of the discrete dipole 
approximation and the finite difference time domain method. (in preparation) 
 [59]  Okamoto H, Macke A, Quante M, Raschke E. Modeling of backscattering by non-spherical ice particles 
for the interpretation of cloud radar signals at 94 GHz. An error analysis. Contrib Atmos Phys 
1995;68:319-334. 
 [60]  Liu CL, Illingworth AJ. Error analysis of backscatter from discrete dipole approximation for different ice 
particle shapes. Atmos Res 1997;44:231-241. 
 [61]  Lemke H, Okamoto H, Quante M. Comment on error analysis of backscatter from discrete dipole 
approximation for different ice particle shapes [ Liu, C.-L., Illingworth, A.J., 1997, Atmos. Res. 44, 231-
241.]. Atmos Res 1998;49:189-197. 
 [62]  Liu CL, Illingworth AJ. Reply to comment by Lemke, Okamoto and Quante on 'Error analysis of 
backscatter from discrete dipole approximation for different ice particle shapes'. Atmos Res 1999;50:1-2. 
 [63]  Yurkin MA, Maltsev VP, Hoekstra AG. Convergence of the discrete dipole approximation. II. An 
extrapolation technique to increase the accuracy. J Opt Soc Am A 2006;23:2592-2601. 
 [64]  Fuller KA, Mackowski DW. Electromagnetic scattering by compounded spherical particles. In: 
Mishchenko MI, Hovenier, JW, Travis, LD, editors. Light Scattering by Nonspherical Particles, Theory, 
Measurements, and Applications. New York: Academic Press, 2000. p. 223-272. 
 [65]  Xu YL. Scattering Mueller matrix of an ensemble of variously shaped small particles. J Opt Soc Am A 
2003;20:2093-2105. 
 [66]  Mackowski DW. Electrostatics analysis of radiative absorption by sphere clusters in the rayleigh limit - 
application to soot particles. Appl Opt 1995;34:3535-3545. 
 [67]  Mackowski DW. Calculation of total cross-sections of multiple-sphere clusters. J Opt Soc Am A 
1994;11:2851-2861. 
 [68]  Ngo D, Videen G, Dalling R. Chaotic light scattering from a system of osculating, conducting spheres. 
Physics Letters A 1997;227:197-202. 
 [69]  Markel VA, Pustovit VN, Karpov SV, Obuschenko AV, Gerasimov VS, Isaev IL. Electromagnetic 
density of states and absorption of radiation by aggregates of nanospheres with multipole interactions. 
Phys Rev B 2004;70:054202. 
 [70]  Kim HY, Sofo JO, Velegol D, Cole MW, Mukhopadhyay G. Static polarizabilities of dielectric 
nanoclusters. Phys Rev A 2005;72:053201. 
 [71]  Jones AR. Electromagnetic wave scattering by assemblies of particles in the Rayleigh approximation. 
Proc R Soc London A 1979;366:111-127. 
 [72]  Jones AR. Scattering efficiency factors for agglomerates for small spheres. J Phys D 1979;12:1661-1672. 
 [73]  Kozasa T, Blum J, Mukai T. Optical-properties of dust aggregates .1. Wavelength dependence. Astron 
Astrophys 1992;263:423-432. 
 [74]  Kozasa T, Blum J, Okamoto H, Mukai T. Optical-properties of dust aggregates .2. Angular-dependence 
of scattered-light. Astron Astrophys 1993;276:278-288. 
 [75]  Lou WJ, Charalampopoulos TT. On the electromagnetic scattering and absorption of agglomerated small 
spherical-particles. J Phys D 1994;27:2258-2270. 
 [76]  Markel VA, Shalaev VM, Stechel EB, Kim W, Armstrong RL. Small-particle composites .1. Linear 
optical properties. Phys Rev B 1996;53:2425-2436. 
 [77]  Pustovit VN, Sotelo JA, Niklasson GA. Coupled multipolar interactions in small-particle metallic 
clusters. J Opt Soc Am A 2002;19:513-518. 
 [78]  Lumme K, Rahola J, Hovenier JW. Light scattering by dense clusters of spheres. Icarus 1997;126:455-
469. 
 [79]  Kimura H, Mann I. Light scattering by large clusters of dipoles as an analog for cometary dust 
aggregates. J Quant Spectrosc Radiat Transf 2004;89:155-164. 
 [80]  Hull P, Shepherd I, Hunt A. Modeling light scattering from Diesel soot particles. Appl Opt 
2004;43:3433-3441. 
 [81]  Venizelos DT, Lou WJ, Charalampopoulos TT. Development of an algorithm for the calculation of the 
scattering properties of agglomerates. Appl Opt 1996;35:542-548. 
 [82]  Voshchinnikov NV, Il'in VB, Henning T. Modelling the optical properties of composite and porous 
interstellar grains. Astron Astrophys 2005;429:371-381. 
 [83]  Kohler M, Kimura H, Mann I. Applicability of the discrete-dipole approximation to light-scattering 
simulations of large cosmic dust aggregates. Astron Astrophys 2006;448:395-399. 
 [84]  Zubko E, Petrov D, Shkuratov Y, Videen G. Discrete dipole approximation simulations of scattering by 
particles with hierarchical structure. Appl Opt 2005;44:6479-6485. 
 [85]  Bourrely C, Chiappetta P, Lemaire TJ, Torresani B. Multidipole formulation of the coupled dipole 
method for electromagnetic scattering by an arbitrary particle. J Opt Soc Am A 1992;9:1336-1340. 
 [86]  Rouleau F, Martin PG. A new method to calculate the extinction properties of irregularly shaped 
particles. Astrophys J 1993;414:803-814. 
 [87]  Mulholland GW, Bohren CF, Fuller KA. Light-scattering by agglomerates - coupled electric and 
magnetic dipole method. Langmuir 1994;10:2533-2546. 
 [88]  Lemaire TJ. Coupled-multipole formulation for the treatment of electromagnetic scattering by a small 
dielectric particle of arbitrary shape. J Opt Soc Am A 1997;14:470-474. 
 [89]  Lakhtakia A. General-theory of the purcell-pennypacker scattering approach and its extension to 
bianisotropic scatterers. Astrophys J 1992;394:494-499. 
 [90]  Loiko VA, Molochko VI. Polymer dispersed liquid crystal droplets: Methods of calculation of optical 
characteristics. Liq Crys 1998;25:603-612. 
 [91]  Smith DA, Stokes KL. Discrete dipole approximation for magneto-optical scattering calculations. Opt 
Expr 2006;14:5746-5754. 
 [92]  Su CC. Electromagnetic scattering by a dielectric body with arbitrary inhomogeneity and anisotropy. 
IEEE Trans Ant Propag 1989;37:384-389. 
 [93]  Chen RS, Fan ZH, Yung EKN. Analysis of electromagnetic scattering of three-dimensional dielectric 
bodies using Krylov subspace FFT iterative methods. Microwave Opt Tech Lett 2003;39:261-267. 
 [94]  Khlebtsov NG. An approximate method for calculating scattering and absorption of light by fractal 
aggregates. Opt Spec 2000;88:594-601. 
 [95]  Markel VA. Coupled-dipole approach to scattering of light from a one-dimensional periodic dipole 
structure. J Mod Opt 1993;40:2281-2291. 
 [96]  Chaumet PC, Rahmani A, Bryant GW. Generalization of the coupled dipole method to periodic 
structures. Phys Rev B 2003;67:165404. 
 [97]  Chaumet PC, Sentenac A. Numerical simulations of the electromagnetic field scattered by defects in a 
double-periodic structure. Phys Rev B 2005;72:205437. 
 [98]  Martin OJF. Efficient scattering calculations in complex backgrounds. AEU-Int J Electr Comm 
2004;58:93-99. 
 [99]  Yang WH, Schatz GC, Vanduyne RP. Discrete dipole approximation for calculating extinction and 
raman intensities for small particles with arbitrary shapes. J Chem Phys 1995;103:869-875. 
 [100]  Lemaire TJ, Bassrei A. Three-dimensional reconstruction of dielectric objects by the coupled-dipole 
method. Appl Opt 2000;39:1272-1278. 
 [101]  Belkebir K, Chaumet PC, Sentenac A. Superresolution in total internal reflection tomography. J Opt Soc 
Am A 2005;22:1889-1897. 
 [102]  Chaumet PC, Belkebir K, Sentenac A. Three-dimensional subwavelength optical imaging using the 
coupled dipole method. Phys Rev B 2004;69:245405. 
 [103]  Chaumet PC, Belkebir K, Lencrerot R. Three-dimensional optical imaging in layered media. Opt Expr 
2006;14:3415-3426. 
 [104]  Zubko E, Shkuratov Y, Videen G. Discrete-dipole analysis of backscatter features of agglomerated 
debris particles comparable in size with wavelength. J Quant Spectrosc Radiat Transf 2006;100:483-488. 
 [105]  Penttila A, Zubko E, Lumme K, Muinonen K, Yurkin MA, Draine BT, Rahola J, Hoekstra AG, 
Shkuratov Y. Comparison between discrete dipole implementations and exact techniques. J Quant 
Spectrosc Radiat Transf 2007, doi:10.1016/j.jqsrt.2007.01.26. 
 [106]  Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical Recipes in C. The Art of Scientific 
Computing. New York: Cambridge University Press, 1990. 
 [107]  Barrett R, Berry M, Chan TF, Demmel J, Donato J, Dongarra J, Eijkhout V, Pozo R, Romine C, van der 
Vorst HA. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. SIAM, 
1994. 
 [108]  Hoekstra AG, Grimminck MD, Sloot PMA. Large scale simulations of elastic light scattering by a fast 
discrete dipole approximation. Int J Mod Phys C 1998;9:87-102. 
 [109]  Zhang SL. GPBi-CG: Generalized product-type methods based on Bi-CG for solving nonsymmetric 
linear systems. SIAM J Sci Comp 1997;18:537-551. 
 [110]  Freund RW. Conjugate gradient-type methods for linear-systems with complex symmetrical coefficient 
matrices. SIAM J Sci Stat Comp 1992;13:425-448. 
 [111]  Flatau PJ. Improvements in the discrete-dipole approximation method of computing scattering and 
absorption. Opt Lett 1997;22:1205-1207. 
 [112]  Fan ZH, Wang DX, Chen RS, Yung EKN. The application of iterative solvers in discrete dipole 
approximation method for computing electromagnetic scattering. Microwave Opt Tech Lett 
2006;48:1741-1746. 
 [113]  Yurkin MA, Maltsev VP, Hoekstra AG. The discrete dipole approximation for simulation of light 
scattering by particles much larger than the wavelength. J Quant Spectrosc Radiat Transf 2007, 
doi:10.1016/j.jqsrt.2007.01.33. 
 [114]  Rahola J. On the eigenvalues of the volume integral operator of electromagnetic scattering. SIAM J Sci 
Comp 2000;21:1740-1754. 
 [115]  Budko NV, Samokhin AB. Spectrum of the volume integral operator of electromagnetic scattering. 
SIAM J Sci Comp 2006;28:682-700. 
 [116]  Budko NV, Samokhin AB, Samokhin AA. A generalized overrelaxation method for solving singular 
volume integral equations in low-frequency scattering problems. Differ Eq 2005;41:1262-1266. 
 [117]  Hoekstra AG, Sloot PMA. Coupled dipole simulations of elastic light scattering on parallel systems. Int J 
Mod Phys C 1995;6:663-679. 
 [118]  Acquista C. Light scattering by tenuous particles: a generalization of the Rayleigh-Gans-Rocard 
approach. Appl Opt 1976;15:2932-2936. 
 [119]  Chiappetta P. Multiple scattering approach to light scattering by arbitrarily shaped particles. J Phys A 
1980;13:2101-2108. 
 [120]  Singham SB, Bohren CF. Light-scattering by an arbitrary particle - the scattering-order formulation of 
the coupled-dipole method. J Opt Soc Am A 1988;5:1867-1872. 
 [121]  de Hoop AT. Convergence criterion for the time-domain iterative Born approximation to scattering by an 
inhomogeneous, dispersive object. J Opt Soc Am A 1991;8:1256-1260. 
 [122]  Flatau PJ, Stephens GL, Draine BT. Light-scattering by rectangular solids in the discrete-dipole 
approximation - a new algorithm exploiting the block-Toeplitz structure. J Opt Soc Am A 1990;7:593-
600. 
 [123]  Flatau PJ. Fast solvers for one dimensional light scattering in the discrete dipole approximation. Opt 
Expr 2004;12:3149-3155. 
 [124]  Goodman JJ, Draine BT, Flatau PJ. Application of fast-Fourier-transform techniques to the discrete-
dipole approximation. Opt Lett 1991;16:1198-1200. 
 [125]  Barrowes BE, Teixeira FL, Kong JA. Fast algorithm for matrix-vector multiply of asymmetric multilevel 
block-Toeplitz matrices in 3-D scattering. Microwave Opt Tech Lett 2001;31:28-32. 
 [126]  Greengard L, Rokhlin V. A fast algorithm for particle simulations. J Comp Phys 1987;73:325-348. 
 [127]  Rahola J. Diagonal forms of the translation operators in the fast multipole algorithm for scattering 
problems. BIT 1996;36:333-358. 
 [128]  Hoekstra AG, Sloot PMA. New computational techniques to simulate light-scattering from arbitrary 
particles. Part Part Sys Charact 1994;11:189-193. 
 [129]  Koc S, Chew WC. Multilevel fast multipole algorithm for the discrete dipole approximation. J Electrom 
Wav Applic 2001;15:1447-1468. 
 [130]  Amini S, Profit ATJ. Multi-level fast multipole solution of the scattering problem. Engin Anal Bound 
Elem 2003;27:547-564. 
 [131]  Darve E. The fast multipole method I: error analysis and asymptotic complexity. SIAM J Num Anal 
2000;38:98-128. 
 [132]  Dembart B, Yip E. The accuracy of fast multipole methods for Maxwell's equations. IEEE Comp Sci 
Engin 1998;5:48-56. 
 [133]  Barnes JE, Hut P. A hierarchical O(N log N) force-calculation algorithm. Nature 1986;324:446-449. 
 [134]  Barnes JE, Hut P. Error analysis of a tree code. Astrophys J Suppl 1989;70:389-417. 
 [135]  Ding KH, Tsang L. A sparse matrix iterative approach for modeling tree scattering. Microwave Opt Tech 
Lett 2003;38:198-202. 
 [136]  Singham MK, Singham SB, Salzman GC. The scattering matrix for randomly oriented particles. J Chem 
Phys 1986;85:3807-3815. 
 [137]  Mishchenko MI. Calculation of the amplitude matrix for a nonspherical particle in a fixed orientation. 
Appl Opt 2000;39:1026-1031. 
 [138]  McClain WM, Ghoul WA. Elastic light scattering by randomly oriented macromolecules: Computation 
of the complete set of observables. J Chem Phys 1986;84:6609-6622. 
 [139]  Khlebtsov NG. Orientational averaging of integrated cross sections in the discrete dipole method. Opt 
Spec 2001;90:408-415. 
 [140]  Mishchenko MI, Travis LD, Mackowski DW. T-matrix computations of light scattering by nonspherical 
particles: A review. J Quant Spectrosc Radiat Transf 1996;55:535-575. 
 [141]  Mackowski DW. Discrete dipole moment method for calculation of the T matrix for nonspherical 
particles. J Opt Soc Am A 2002;19:881-893. 
 [142]  Mishchenko MI. Light-scattering by size shape distributions of randomly oriented axially-symmetrical 
particles of a size comparable to a wavelength. Appl Opt 1993;32:4652-4666. 
 [143]  Muinonen K, Zubko E. Optimizing the discrete-dipole approximation for sequences of scatterers with 
identical shapes but differing sizes or refractive indices. J Quant Spectrosc Radiat Transf 2006;100:288-
294. 
 [144]  Hovenier JW, Lumme K, Mishchenko MI, Voshchinnikov NV, Mackowski DW, Rahola J. 
Computations of scattering matrices of four types of non-spherical particles using diverse methods. J 
Quant Spectrosc Radiat Transf 1996;55:695-705. 
 [145]  Wriedt T, Comberg U. Comparison of computational scattering methods. J Quant Spectrosc Radiat 
Transf 1998;60:411-423. 
 [146]  Comberg U, Wriedt T. Comparison of scattering calculations for aggregated particles based on different 
models. J Quant Spectrosc Radiat Transf 1999;63:149-162. 
 [147]  Wriedt T, Hellmers J, Eremina E, Schuh R. Light scattering by single erythrocyte: Comparison of 
different methods. J Quant Spectrosc Radiat Transf 2006;100:444-456. 
 [148]  Laczik Z. Discrete-dipole-approximation-based light-scattering calculations for particles with a real 
refractive index smaller than unity. Appl Opt 1996;35:3736-3745. 
	1 Introduction
	2 General framework
	3 Various DDA models
	3.1 Theoretical base of the DDA
	3.2 Accuracy of DDA simulations
	3.3 The DDA for clusters of spheres
	3.4 Modifications and extensions of the DDA
	4 Numerical considerations
	4.1 Direct vs. iterative methods
	4.2 Scattering order formulation
	4.3 Block-Toeplitz
	4.4 FFT
	4.5 Fast multipole method
	4.6 Orientation averaging and repeated calculations
	5 Comparison of the DDA to other methods
	6 Concluding remarks
	Acknowledgements
	Appendix. Description of used acronyms and symbols
	References
  /ASCII85EncodePages false
  /AllowTransparency false
  /AutoPositionEPSFiles true
  /AutoRotatePages /None
  /Binding /Left
  /CalGrayProfile (Dot Gain 20%)
  /CalRGBProfile (sRGB IEC61966-2.1)
  /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2)
  /sRGBProfile (sRGB IEC61966-2.1)
  /CannotEmbedFontPolicy /Error
  /CompatibilityLevel 1.4
  /CompressObjects /Tags
  /CompressPages true
  /ConvertImagesToIndexed true
  /PassThroughJPEGImages true
  /CreateJDFFile false
  /CreateJobTicket false
  /DefaultRenderingIntent /Default
  /DetectBlends true
  /DetectCurves 0.0000
  /ColorConversionStrategy /CMYK
  /DoThumbnails false
  /EmbedAllFonts true
  /EmbedOpenType false
  /ParseICCProfilesInComments true
  /EmbedJobOptions true
  /DSCReportingLevel 0
  /EmitDSCWarnings false
  /EndPage -1
  /ImageMemory 1048576
  /LockDistillerParams false
  /MaxSubsetPct 100
  /Optimize true
  /OPM 1
  /ParseDSCComments true
  /ParseDSCCommentsForDocInfo true
  /PreserveCopyPage true
  /PreserveDICMYKValues true
  /PreserveEPSInfo true
  /PreserveFlatness true
  /PreserveHalftoneInfo false
  /PreserveOPIComments true
  /PreserveOverprintSettings true
  /StartPage 1
  /SubsetFonts true
  /TransferFunctionInfo /Apply
  /UCRandBGInfo /Preserve
  /UsePrologue false
  /ColorSettingsFile ()
  /AlwaysEmbed [ true
  /NeverEmbed [ true
  /AntiAliasColorImages false
  /CropColorImages true
  /ColorImageMinResolution 300
  /ColorImageMinResolutionPolicy /OK
  /DownsampleColorImages true
  /ColorImageDownsampleType /Bicubic
  /ColorImageResolution 300
  /ColorImageDepth -1
  /ColorImageMinDownsampleDepth 1
  /ColorImageDownsampleThreshold 1.50000
  /EncodeColorImages true
  /ColorImageFilter /DCTEncode
  /AutoFilterColorImages true
  /ColorImageAutoFilterStrategy /JPEG
  /ColorACSImageDict <<
    /QFactor 0.15
    /HSamples [1 1 1 1] /VSamples [1 1 1 1]
  /ColorImageDict <<
    /QFactor 0.15
    /HSamples [1 1 1 1] /VSamples [1 1 1 1]
  /JPEG2000ColorACSImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 30
  /JPEG2000ColorImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 30
  /AntiAliasGrayImages false
  /CropGrayImages true
  /GrayImageMinResolution 300
  /GrayImageMinResolutionPolicy /OK
  /DownsampleGrayImages true
  /GrayImageDownsampleType /Bicubic
  /GrayImageResolution 300
  /GrayImageDepth -1
  /GrayImageMinDownsampleDepth 2
  /GrayImageDownsampleThreshold 1.50000
  /EncodeGrayImages true
  /GrayImageFilter /DCTEncode
  /AutoFilterGrayImages true
  /GrayImageAutoFilterStrategy /JPEG
  /GrayACSImageDict <<
    /QFactor 0.15
    /HSamples [1 1 1 1] /VSamples [1 1 1 1]
  /GrayImageDict <<
    /QFactor 0.15
    /HSamples [1 1 1 1] /VSamples [1 1 1 1]
  /JPEG2000GrayACSImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 30
  /JPEG2000GrayImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 30
  /AntiAliasMonoImages false
  /CropMonoImages true
  /MonoImageMinResolution 1200
  /MonoImageMinResolutionPolicy /OK
  /DownsampleMonoImages true
  /MonoImageDownsampleType /Bicubic
  /MonoImageResolution 1200
  /MonoImageDepth -1
  /MonoImageDownsampleThreshold 1.50000
  /EncodeMonoImages true
  /MonoImageFilter /CCITTFaxEncode
  /MonoImageDict <<
    /K -1
  /AllowPSXObjects false
  /CheckCompliance [
    /None
  /PDFX1aCheck false
  /PDFX3Check false
  /PDFXCompliantPDFOnly false
  /PDFXNoTrimBoxError true
  /PDFXTrimBoxToMediaBoxOffset [
    0.00000
    0.00000
    0.00000
    0.00000
  /PDFXSetBleedBoxToMediaBox true
  /PDFXBleedBoxToTrimBoxOffset [
    0.00000
    0.00000
    0.00000
    0.00000
  /PDFXOutputIntentProfile ()
  /PDFXOutputConditionIdentifier ()
  /PDFXOutputCondition ()
  /PDFXRegistryName ()
  /PDFXTrapped /False
  /Description <<
    /CHS <FEFF4f7f75288fd94e9b8bbe5b9a521b5efa7684002000410064006f006200650020005000440046002065876863900275284e8e9ad88d2891cf76845370524d53705237300260a853ef4ee54f7f75280020004100630072006f0062006100740020548c002000410064006f00620065002000520065006100640065007200200035002e003000204ee553ca66f49ad87248672c676562535f00521b5efa768400200050004400460020658768633002>
    /CHT <FEFF4f7f752890194e9b8a2d7f6e5efa7acb7684002000410064006f006200650020005000440046002065874ef69069752865bc9ad854c18cea76845370524d5370523786557406300260a853ef4ee54f7f75280020004100630072006f0062006100740020548c002000410064006f00620065002000520065006100640065007200200035002e003000204ee553ca66f49ad87248672c4f86958b555f5df25efa7acb76840020005000440046002065874ef63002>
    /DAN <FEFF004200720075006700200069006e0064007300740069006c006c0069006e006700650072006e0065002000740069006c0020006100740020006f007000720065007400740065002000410064006f006200650020005000440046002d0064006f006b0075006d0065006e007400650072002c0020006400650072002000620065006400730074002000650067006e006500720020007300690067002000740069006c002000700072006500700072006500730073002d007500640073006b007200690076006e0069006e00670020006100660020006800f8006a0020006b00760061006c0069007400650074002e0020004400650020006f007000720065007400740065006400650020005000440046002d0064006f006b0075006d0065006e0074006500720020006b0061006e002000e50062006e00650073002000690020004100630072006f00620061007400200065006c006c006500720020004100630072006f006200610074002000520065006100640065007200200035002e00300020006f00670020006e0079006500720065002e>
    /DEU <FEFF00560065007200770065006e00640065006e0020005300690065002000640069006500730065002000450069006e007300740065006c006c0075006e00670065006e0020007a0075006d002000450072007300740065006c006c0065006e00200076006f006e002000410064006f006200650020005000440046002d0044006f006b0075006d0065006e00740065006e002c00200076006f006e002000640065006e0065006e002000530069006500200068006f006300680077006500720074006900670065002000500072006500700072006500730073002d0044007200750063006b0065002000650072007a0065007500670065006e0020006d00f60063006800740065006e002e002000450072007300740065006c006c007400650020005000440046002d0044006f006b0075006d0065006e007400650020006b00f6006e006e0065006e0020006d006900740020004100630072006f00620061007400200075006e0064002000410064006f00620065002000520065006100640065007200200035002e00300020006f0064006500720020006800f600680065007200200067006500f600660066006e00650074002000770065007200640065006e002e>
    /ESP <FEFF005500740069006c0069006300650020006500730074006100200063006f006e0066006900670075007200610063006900f3006e0020007000610072006100200063007200650061007200200064006f00630075006d0065006e0074006f00730020005000440046002000640065002000410064006f0062006500200061006400650063007500610064006f00730020007000610072006100200069006d0070007200650073006900f3006e0020007000720065002d0065006400690074006f007200690061006c00200064006500200061006c00740061002000630061006c0069006400610064002e002000530065002000700075006500640065006e00200061006200720069007200200064006f00630075006d0065006e0074006f00730020005000440046002000630072006500610064006f007300200063006f006e0020004100630072006f006200610074002c002000410064006f00620065002000520065006100640065007200200035002e003000200079002000760065007200730069006f006e0065007300200070006f00730074006500720069006f007200650073002e>
    /FRA <FEFF005500740069006c006900730065007a00200063006500730020006f007000740069006f006e00730020006100660069006e00200064006500200063007200e900650072002000640065007300200064006f00630075006d0065006e00740073002000410064006f00620065002000500044004600200070006f0075007200200075006e00650020007100750061006c0069007400e90020006400270069006d007000720065007300730069006f006e00200070007200e9007000720065007300730065002e0020004c0065007300200064006f00630075006d0065006e00740073002000500044004600200063007200e900e90073002000700065007500760065006e0074002000ea0074007200650020006f007500760065007200740073002000640061006e00730020004100630072006f006200610074002c002000610069006e00730069002000710075002700410064006f00620065002000520065006100640065007200200035002e0030002000650074002000760065007200730069006f006e007300200075006c007400e90072006900650075007200650073002e>
    /ITA <FEFF005500740069006c0069007a007a006100720065002000710075006500730074006500200069006d0070006f007300740061007a0069006f006e00690020007000650072002000630072006500610072006500200064006f00630075006d0065006e00740069002000410064006f00620065002000500044004600200070006900f900200061006400610074007400690020006100200075006e00610020007000720065007300740061006d0070006100200064006900200061006c007400610020007100750061006c0069007400e0002e0020004900200064006f00630075006d0065006e007400690020005000440046002000630072006500610074006900200070006f00730073006f006e006f0020006500730073006500720065002000610070006500720074006900200063006f006e0020004100630072006f00620061007400200065002000410064006f00620065002000520065006100640065007200200035002e003000200065002000760065007200730069006f006e006900200073007500630063006500730073006900760065002e>
    /JPN <FEFF9ad854c18cea306a30d730ea30d730ec30b951fa529b7528002000410064006f0062006500200050004400460020658766f8306e4f5c6210306b4f7f75283057307e305930023053306e8a2d5b9a30674f5c62103055308c305f0020005000440046002030d530a130a430eb306f3001004100630072006f0062006100740020304a30883073002000410064006f00620065002000520065006100640065007200200035002e003000204ee5964d3067958b304f30533068304c3067304d307e305930023053306e8a2d5b9a306b306f30d530a930f330c8306e57cb30818fbc307f304c5fc59808306730593002>
    /KOR <FEFFc7740020c124c815c7440020c0acc6a9d558c5ec0020ace0d488c9c80020c2dcd5d80020c778c1c4c5d00020ac00c7a50020c801d569d55c002000410064006f0062006500200050004400460020bb38c11cb97c0020c791c131d569b2c8b2e4002e0020c774b807ac8c0020c791c131b41c00200050004400460020bb38c11cb2940020004100630072006f0062006100740020bc0f002000410064006f00620065002000520065006100640065007200200035002e00300020c774c0c1c5d0c11c0020c5f40020c2180020c788c2b5b2c8b2e4002e>
    /NLD (Gebruik deze instellingen om Adobe PDF-documenten te maken die zijn geoptimaliseerd voor prepress-afdrukken van hoge kwaliteit. De gemaakte PDF-documenten kunnen worden geopend met Acrobat en Adobe Reader 5.0 en hoger.)
    /NOR <FEFF004200720075006b00200064006900730073006500200069006e006e007300740069006c006c0069006e00670065006e0065002000740069006c002000e50020006f0070007000720065007400740065002000410064006f006200650020005000440046002d0064006f006b0075006d0065006e00740065007200200073006f006d00200065007200200062006500730074002000650067006e0065007400200066006f00720020006600f80072007400720079006b006b0073007500740073006b00720069006600740020006100760020006800f800790020006b00760061006c0069007400650074002e0020005000440046002d0064006f006b0075006d0065006e00740065006e00650020006b0061006e002000e50070006e00650073002000690020004100630072006f00620061007400200065006c006c00650072002000410064006f00620065002000520065006100640065007200200035002e003000200065006c006c00650072002000730065006e006500720065002e>
    /PTB <FEFF005500740069006c0069007a006500200065007300730061007300200063006f006e00660069006700750072006100e700f50065007300200064006500200066006f0072006d00610020006100200063007200690061007200200064006f00630075006d0065006e0074006f0073002000410064006f0062006500200050004400460020006d00610069007300200061006400650071007500610064006f00730020007000610072006100200070007200e9002d0069006d0070007200650073007300f50065007300200064006500200061006c007400610020007100750061006c00690064006100640065002e0020004f007300200064006f00630075006d0065006e0074006f00730020005000440046002000630072006900610064006f007300200070006f00640065006d0020007300650072002000610062006500720074006f007300200063006f006d0020006f0020004100630072006f006200610074002000650020006f002000410064006f00620065002000520065006100640065007200200035002e0030002000650020007600650072007300f50065007300200070006f00730074006500720069006f007200650073002e>
    /SUO <FEFF004b00e40079007400e40020006e00e40069007400e4002000610073006500740075006b007300690061002c0020006b0075006e0020006c0075006f00740020006c00e400680069006e006e00e4002000760061006100740069007600610061006e0020007000610069006e006100740075006b00730065006e002000760061006c006d0069007300740065006c00750074007900f6006800f6006e00200073006f00700069007600690061002000410064006f0062006500200050004400460020002d0064006f006b0075006d0065006e007400740065006a0061002e0020004c0075006f0064007500740020005000440046002d0064006f006b0075006d0065006e00740069007400200076006f0069006400610061006e0020006100760061007400610020004100630072006f0062006100740069006c006c00610020006a0061002000410064006f00620065002000520065006100640065007200200035002e0030003a006c006c00610020006a006100200075007500640065006d006d0069006c006c0061002e>
    /SVE <FEFF0041006e007600e4006e00640020006400650020006800e4007200200069006e0073007400e4006c006c006e0069006e006700610072006e00610020006f006d002000640075002000760069006c006c00200073006b006100700061002000410064006f006200650020005000440046002d0064006f006b0075006d0065006e007400200073006f006d002000e400720020006c00e4006d0070006c0069006700610020006600f60072002000700072006500700072006500730073002d007500740073006b00720069006600740020006d006500640020006800f600670020006b00760061006c0069007400650074002e002000200053006b006100700061006400650020005000440046002d0064006f006b0075006d0065006e00740020006b0061006e002000f600700070006e00610073002000690020004100630072006f0062006100740020006f00630068002000410064006f00620065002000520065006100640065007200200035002e00300020006f00630068002000730065006e006100720065002e>
    /ENU (Use these settings to create Adobe PDF documents best suited for high-quality prepress printing.  Created PDF documents can be opened with Acrobat and Adobe Reader 5.0 and later.)
  /Namespace [
    (Adobe)
    (Common)
    (1.0)
  /OtherNamespaces [
    <<
      /AsReaderSpreads false
      /CropImagesToFrames true
      /ErrorControl /WarnAndContinue
      /FlattenerIgnoreSpreadOverrides false
      /IncludeGuidesGrids false
      /IncludeNonPrinting false
      /IncludeSlug false
      /Namespace [
        (Adobe)
        (InDesign)
        (4.0)
      ]
      /OmitPlacedBitmaps false
      /OmitPlacedEPS false
      /OmitPlacedPDF false
      /SimulateOverprint /Legacy
    >>
    <<
      /AddBleedMarks false
      /AddColorBars false
      /AddCropMarks false
      /AddPageInfo false
      /AddRegMarks false
      /ConvertColors /ConvertToCMYK
      /DestinationProfileName ()
      /DestinationProfileSelector /DocumentCMYK
      /Downsample16BitImages true
      /FlattenerPreset <<
        /PresetSelector /MediumResolution
      >>
      /FormElements false
      /GenerateStructure false
      /IncludeBookmarks false
      /IncludeHyperlinks false
      /IncludeInteractive false
      /IncludeLayers false
      /IncludeProfiles false
      /MultimediaHandling /UseObjectSettings
      /Namespace [
        (Adobe)
        (CreativeSuite)
        (2.0)
      ]
      /PDFXOutputIntentProfileSelector /DocumentCMYK
      /PreserveEditing true
      /UntaggedCMYKHandling /LeaveUntagged
      /UntaggedRGBHandling /UseDocumentProfile
      /UseDocumentBleed false
    >>
>> setdistillerparams
  /HWResolution [2400 2400]
  /PageSize [612.000 792.000]
>> setpagedevice

  We present a review of the discrete dipole approximation (DDA), which is a
general method to simulate light scattering by arbitrarily shaped particles. We
put the method in historical context and discuss recent developments, taking
the viewpoint of a general framework based on the integral equations for the
electric field. We review both the theory of the DDA and its numerical aspects,
the latter being of critical importance for any practical application of the
method. Finally, the position of the DDA among other methods of light
scattering simulation is shown and possible future developments are discussed.

Introduction ......................................................................................................................... 2 
2 General framework.............................................................................................................. 3 
3 Various DDA models .......................................................................................................... 7 
3.1 Theoretical base of the DDA........................................................................................ 7 
3.2 Accuracy of DDA simulations ................................................................................... 13 
3.3 The DDA for clusters of spheres................................................................................ 16 
3.4 Modifications and extensions of the DDA................................................................. 18 
4 Numerical considerations.................................................................................................. 19 
4.1 Direct vs. iterative methods........................................................................................ 19 
4.2 Scattering order formulation ...................................................................................... 22 
4.3 Block-Toeplitz ........................................................................................................... 23 
4.4 FFT............................................................................................................................. 24 
4.5 Fast multipole method................................................................................................ 24 
4.6 Orientation averaging and repeated calculations ....................................................... 25 
5 Comparison of the DDA to other methods ....................................................................... 27 
6 Concluding remarks .......................................................................................................... 28 
Acknowledgements .................................................................................................................. 28 
Appendix. Description of used acronyms and symbols ........................................................... 28 
References ................................................................................................................................ 31 
1 Introduction 
The discrete dipole approximation (DDA) is a general method to compute scattering and 
absorption of electromagnetic waves by particles of arbitrary geometry and composition. 
Initially the DDA was proposed by Purcell and Pennypacker (PP) [1], who replaced the 
scatterer by a set of point dipoles. These dipoles interact with each other and the incident 
field, giving rise to a system of linear equations, which is solved to obtain dipole 
polarizations. All the measured scattering quantities can be obtained from these polarizations. 
The DDA was further developed by Draine and coworkers [2-5], who popularized the method 
by developing a publicly available computer code DDSCAT [6]. Later it was shown that the 
DDA also can be derived from the integral equation for the electric field, which is discretized 
by dividing the scatterer into small cubical subvolumes. This derivation was apparently first 
performed by Goedecke and O'Brien [7] and further developed by others (see, for instance, 
[8-11]). It is important to note that the final equations, produced by both lines of derivation of 
the DDA are essentially the same. The only difference is that derivations based on the integral 
equations give more mathematical insight into the approximation, thus pointing at ways to 
improve the method, while the model based on point dipoles is physically clearer. 
The DDA is called the coupled dipole method or approximation by some researchers 
[12,13]. There are also other methods, such as the volume integral equation formulation [14] 
and the digitized Green’s function (DGF) [7], which were developed completely 
independently from PP. However, later they were shown to be equivalent to DDA [8,15]. In 
this review we will use the term DDA to refer to all such methods, since we describe them in 
terms of one general framework. However, it is difficult to separate unambiguously the DDA 
from other similar methods, based on the volume integral equations for the electromagnetic 
fields, such as a broad range of method of moments (MoM) with different bases and testing 
functions [16-19]. In our opinion, one fundamental aspect of the DDA is that the solution for 
the “physically meaningful” internal fields or their direct derivatives, e.g. polarization, plays 
an integral role in the process. In other words, any DDA formulation can be interpreted as 
replacing a scatterer by a set of interacting dipoles; this is further discussed in Section 2. An 
example of method that is not considered DDA is the MoM with higher-order hierarchical 
Legendre basis functions [17]. 
The DDA is a popular method in the light-scattering community and it has been 
reviewed by several authors. An extensive review by Draine and Flatau [4] covers almost all 
DDA developments up to 1994. A more recent review by Draine [5] mainly concerns 
applications and numerical considerations. DDA theory was discussed together with other 
methods for light scattering simulations in reviews by Wriedt [20], Chiappetta and Torresani 
[21], and Kahnert [15] and in books by Mishchenko et al. [22] and Tsang et al. [23]. Jones 
[24] placed the DDA in context of different methods with respect to particle characterization. 
However, many important DDA developments since 1994 are not mentioned in any of these 
manuscripts. Those that are mentioned are usually considered as side-steps, and are not placed 
into a general framework. Moreover, to the best of our knowledge numerical aspects of the 
DDA have never been reviewed extensively – each paper discusses only a few particular 
aspects. In this review we try to fill this gap. 
A general framework is developed in Section 2 to ease the further discussion of different 
DDA models. This framework is based on the integral equation because it allows a uniform 
description of all the DDA development. However, connection to a physically clearer model 
of point dipoles is discussed throughout the section. The sources of errors in the DDA 
formulation are also discussed there. 
In Section 3 the physical principles of the DDA are reviewed and results of different 
models are compared. In Subsection 3.1 different improvements of polarizabilities and 
interaction terms are reviewed from a theoretical point of view. Different expressions for Cabs 
also are discussed. Comparison of simulation results using different formulations is given in 
Subsection 3.2. Subsection 3.3 covers the special case of a cluster of spheres that allows 
particular improvements and simplifications. In Section 3.4 different significant modifications 
are reviewed, which do not fall completely into the general framework described in Section 2. 
Enhancements of the DDA for some special purposes also are discussed. 
Different numerical aspects of the DDA are reviewed in Section 4. These are concerned 
primarily with solving very large systems of linear equations (Subsection 4.1). Subsection 4.2 
describes the simplest iterative procedure to solve DDA linear system, which has a clear 
physical meaning. The special structure of the DDA interaction matrix for a rectangular grid 
and its application to decrease computational costs are described in subsections 4.3 and 4.4 
respectively. General methods to accelerate calculations, which do not require a rectangular 
grid, are discussed in Subsection 4.5. Subsection 4.6 covers special techniques to increase the 
efficiency of repeated calculations (e.g. in orientation averaging). 
A numerical comparison of the DDA with other methods is reviewed in Section 5; its 
strong and weak points are discussed. Section 6 concludes the review and discusses future 
development of the DDA. 
2 General framework 
The )iexp( tω−  time dependence of all fields is assumed throughout this review. The scatterer 
is assumed dielectric but not magnetic (magnetic permittivity 1=μ ). The electric permittivity 
is assumed isotropic to simplify the derivations; however, extension to arbitrary dielectric 
tensors is straightforward.1
The general form of the integral equation governing the electric field inside the 
dielectric scatterer is the following [8,15]: 
)()(),(),()()(),(d)()( 00
rErrLrMrErrrGrErE χχ VVr
∂−+′′′′+= ∫ , (1)
                                                 
1 In most formulae scalar values can be replaced directly by tensors, but there are exceptions. Extensions of 
DDA to optically anisotropic scatterers are discussed in Section 3.4. 
where Einc(r) and E(r) are the incident and total electric field at location r; 
πεχ 4)1)(()( −= rr  is the susceptibility of the medium at point r (ε(r) – relative 
permittivity). V is the volume of the particle, i.e., the volume that contains all points where the 
susceptibility is not zero. V0 is a smaller volume such that , VV ⊂0 00 \ VV ∂∈r . ),( rrG ′  is 
the free space dyadic Green’s function, defined as 
−=∇∇+=′ 222
i1ˆˆ)iexp()iexp(ˆˆ),(
k IIIrrG , (2)
where I  is the identity dyadic, ck ω=  is the free space wave vector, rrR ′−= , R=R , 
and  is a dyadic defined as  (μ and ν  are Cartesian components of the vector 
or tensor). M is the following integral associated with the finiteness of the exclusion volume 
RR ˆˆ νμμν RRRR =ˆˆ
( )∫ ′−′′′′=
)()(),()()(),(d),( s30
rV rErrrGrErrrGrM χχ , (3)
where ),(s rrG ′  is the static limit ( ) of 0→k ),( rrG ′ : 
−−=∇∇=′
11ˆˆ),(
IrrG . (4)
L  is the so-called self-term dyadic: 
rV rL , (5)
where  is an external normal to the surface ∂Vn′ˆ 0 at point r'. L  is always a real, symmetric 
dyadic with trace equal to 4π [25]. It is important to note that L  does not depend on the size 
of the volume V0, but only on its shape (and location of the point r inside it). On the contrary, 
M does depend on the size of the volume, moreover it approaches zero when the size of the 
volume decreases [8] (if both χ(r) and E(r) are continuous inside V0). 
When deriving Eq. (1) the singularity of the Green’s function has been treated explicitly, 
therefore it is preferable to the commonly used formulation [8,15]: 
∫ ′′′′+=
r )()(),(d)()( 3inc rErrrGrErE χ . (6)
Moreover, Yanghjian noted [25] that there exist several methods for treating the singularity in 
Eq. (6) leading to different results. He also proved that the derivation of Eq. (6) is false in the 
vicinity of the singularity of ),( rrG ′ . Hence it can be considered correct only if the 
singularity is then treated in a way similar to that of Lakhtakia [8], resulting in the correct 
Eq. (1). 
Discretization of Eq. (1) is done in the following way [15]. Let , U
= /0=ji VV I  
for ji ≠ ; N denotes number of subvolumes.2 Although the formulation is applicable to any 
set of subvolumes Vi, in most applications standard (equal) cells are used. Then the shape of 
the scatterer cannot always be described exactly by such standard cells. Hence, the 
discretization may be only approximately correct. Assuming iV∈r  and choosing iVV =0 , 
Eq. (1) can be rewritten as 
)()(),(),()()(),(d)()( 3inc rErrLrMrErrrGrErE χχ ii
∂−+′′′′+= ∑ ∫
. (7)
                                                 
2 In the framework of the DDA we usually call a subvolume a dipole. 
The set of Eq. (7) (for all i) is exact. Further, one fixed point ri inside each Vi (its center) is 
chosen and  is set. In many cases the following assumptions can be made: irr =
)()()()(),(d3 jjijj
rErGrErrrG χχ =′′′′∫ , (8)
)()(),( iiiiiV rErMrM χ= , (9)
which state that integrals in Eq. (7) linearly depend upon the values of χ and E at point ri. 
Eq. (7) can then be rewritten as 
( ) iiii
jjjijii V ELMEGEE χχ −++= ∑
inc , (10)
where , , )( ii rEE = )(
incinc
ii rEE = )( ii rχχ = , ),( iii V rLL ∂= . 
The usual approximation [15] is to consider E and χ constant inside each subvolume: 
iii V∈== rrErE for)(,)( χχ , (11)
which automatically implies Eqs. (8), (9) and 
( )∫ ′−′′=
iii r ),(),(d
s3)0( rrGrrGM , (12)
∫ ′′=
ij rV
1 3)0( rrGG . (13)
Superscript (0) denotes approximate values of the dyadics. A further approximation, which is 
used in almost all formulations of the DDA, including e.g. [8], is  
),()0( jiij rrGG = . (14)
This assumption is made implicitly by all formulations that start by replacing the scatterer 
with a set of point dipoles. It is important to note that Eq. (10) and derivations resulting from 
it require weaker assumptions (Eqs. (8), (9)) than imposed by Eq. (11) and, moreover, 
Eq. (14). It is possible to formulate the DDA based on Eq. (10), e.g. the Peltoniemi 
formulation [26] that is described in Section 3.1. We postulate Eq. (10) as a distinctive feature 
of the DDA, i.e. a method is called the DDA if and only if its main equation is equivalent to 
Eq. (10) with any Vi, χi, iM , iL , and ijG . 
Kahnert [15] distinguished the DDA from the MoM by the fact that the MoM solves 
directly Eq. (10) for unknown Ei, while the DDA seeks not the total, but the exciting electric 
fields 
( )( ) selfexc iiiiiii EEEMLIE −=−+= χ , (15)
( ) iiiii ELME χ−=self , (16)
where  is the field induced by the subvolume on itself. Eq. selfiE (10) is then equivalent to  
jjijii
excexcinc EαGEE , (17)
where iα  is the polarizability tensor defined as 
( )( ) 1−−+= iiiiii V χχ MLIα . (18)
However, an alternative formulation of the DDA exists [4] seeking a solution for unknown 
polarizations Pi: 
iiiiii V EEαP χ==
exc , (19)
jijiii PGPαE
1inc . (20)
It is important to note that Pi, defined by Eq. (19), is only an approximation to the polarization 
of the subvolume Vi. This approximation is exact only under the assumption of Eq. (11), while 
the formulation itself does not require it. The formulation, using Eq. (20), can be thought as 
an intermediary between the DDA and the MoM as classified by Kahnert [15], therefore 
revealing complete equivalence of these two formulations. The special structure of the matrix 
ijG  makes Eq. (20) preferable over Eqs. (10), (17) to find a numerical solution. This is 
discussed in Section 4. 
Lakhtakia [8] classified strong and weak forms of the DDA as those accounting for or 
neglecting iM  respectively. The weak form approaches the strong form when the size of the 
cell decreases, because iM  approaches zero. For a cubical cell Vi and with ri located at the 
center of the cell, iL  can be calculated analytically yielding [25] 
=i . (21)
Using Eq. (18), this results in the well-known Clausius-Mossotti (CM) polarizability (used 
originally by Purcell and Pennypacker [1]) for the weak form of the DDA: 
ii d ε
α IIα , (22)
where )( ii rεε = , and d is the size of the cubical cell. 
After the internal electric fields are determined, the scattered fields and cross sections 
can be calculated. The scattered fields are obtained by taking the limit ∞→r  of the integral 
in Eq. (1) (see e.g. [7]): 
)iexp(
)(sca nFrE
= , (23)
where rrn =  is the unit vector in the scattering direction, and F is the scattering amplitude: 
∑∫ ′′⋅′−′−−=
krnnk )()()iexp(d)ˆˆ(i)( 33 rErnrInF χ . (24)
All other differential scattering properties, such as amplitude and Mueller scattering matrices, 
and asymmetry parameter >< θcos  can be derived from F(n), calculated for two incident 
polarizations [27]. Radiation forces also can be calculated [28-30]. Consider an incident 
polarized plane wave3
)iexp()( 0inc rkerE ⋅= , (25)
where , a is the incident direction, and ak k= 10 =e . The scattering cross section Csca is [27] 
2sca )(d
C . (26)
Absorption and extinction cross sections (Cabs, Cext) are derived [7,14] directly from the 
internal fields: 
( )∑ ∫ ′′′=
abs )()(Imd4 rErχπ , (27)
[ ]( ) ( )∗∗ ⋅=′⋅′′′= ∑∫ 02inc3ext )(Re4)()()(Imd4 eaFrErEr krkC i Vi
χπ , (28)
where * denotes a complex conjugate. Conservation of energy necessitates that 
absextsca CCC −= . (29)
However, as was noted by Draine [2], use of Eq. (29) for evaluation of Csca can lead to larger 
errors than Eq. (26), especially when . scaabs CC >>
The easiest way to express Eqs. (24) and (27) in terms of the internal fields in the 
subvolumes centers is to assume Eq. (11), yielding 
∑ ∫ ⋅′−′−−=
krnnk )iexp(d)ˆˆ(i)( 33)0( nrEInF χ , (30)
                                                 
3 DDA can be used for any incident wave, e.g. Gaussian beams [31]; however, we do not discuss this here. 
∑∑ ∗==
iii kVkC )Im(4)Im(4
abs EPE πχπ . (31)
Further approximation of Eq. (30), leaving only the lowest order expansion of the exponent 
around ri, leads to 
∑ ⋅−−−=
ii knnk )iexp()ˆˆ(i)(
3)0( nrPInF , (32)
which together with Eq. (28), leads to 
( )∑ ∗⋅=
inc)0(
ext Im4 EPπ . (33)
Eqs. (32) and (33) are identical to those used by Purcell and Pennypacker [1] and then by 
Draine [2], while expressions for Cabs (compared to Eq. (31)) are slightly different. These 
differences are discussed in Subsection 3.1. Unfortunately, many researchers do not specify 
explicitly how the scattering quantities are obtained from the computed internal fields or 
polarizations. Those who do usually use Draine’s prescription (Eqs. (26), (32), (33), and 
(35)). 
Errors of the formulation can be classified as associated with the finite cell size d 
(discretization errors), and with approximating the particle shape with a set of standard cells, 
e.g. cubical (shape errors). Discretization errors result from considering E constant inside 
each cell and the approximate evaluation of iM  and ijG . Shape errors also can be considered 
as resulting from the assumption of constant χ and E inside bordering cells, which is false 
since the edge of the particle crosses these cells. On the other hand, shape errors can be 
viewed as a difference of the results for the exact particle shape and for that comprised of the 
set of standard cells. Both errors approach zero when ∞→N , while the geometry of the 
scatterer and parameters of the incident field are fixed. However, the same does not apply if 
 while N is fixed, i.e. the DDA is not exact in the long-wavelength limit. Moreover, 
both errors are sensitive to the size of the scatterer in the resonance region (see discussion in 
Subsection 
3.2). The behavior of these errors was studied by Yurkin et al. [32]. 
3 Various DDA models 
3.1 Theoretical base of the DDA 
Since the original manuscript by Purcell and Pennypacker [1], many attempts have been made 
to improve the DDA. The first stage (1988-1993) of these improvements was reviewed by 
Draine and Flatau [4]. It has been noted [2] that Eq. (22) does not satisfy energy conservation, 
and results obtained using this formulation do not satisfy the optical theorem. Based on the 
well-known [33] “radiative reaction” (RR) electric field, a correction to the polarizability for a 
finite dipole was added [2]: 
i)32(1 α
= . (34)
Draine [2] also proposed the following expression for the absorption cross section: 
[ ]∑ ⋅−⋅=
iiii kkC
*3*exc)0(
abs )32()Im(4 PPEPπ , (35)
derived from Eq. (29) applied to a single point dipole. The PP formulation uses Eq. (35) 
without the second part. It can be verified that Eq. (35) results in zero absorption for any 
scatterer if the polarizability is of the following form: 
IAα 31 i)32( kii −=
− , Hii AA = , (36)
where H denotes the conjugate transpose of a tensor. For real refractive index m, RR and all 
other expressions specified below result in α  satisfying Eq. (36), which makes Eq. (35) 
clearly favorable over e.g. the PP formulation. It must be noted however that the original PP 
formulation, where CM polarizability was used, also results in zero absorption for real m. 
The correction in Eq. (34) is ( )3)(O kd . Several other corrections of ( )2)(O kd  have been 
proposed. The first one was proposed by Goedecke and O’Brien [7] and independently in two 
other manuscripts [34,35]. They started from Eqs. (10)-(12) and used the following 
simplifying fact for a cubical cell (also valid for spherical cells), resulting from symmetry: 
)(d IRRf
RRf , (37)
where the origin is in the center of the cube. Eq. (37) is valid for any f(R) that has a singularity 
of less than third order for , i.e. the integrals on both sides are defined. They obtained 0→R
32)0( )iexp(d
Rki IM . (38)
By expanding exp(ikR) in Taylor series one can obtain 
( )⎟⎟
++= ∫ 423
2)0( Oi
ki IM . (39)
The remaining integral was evaluated by approximating the cube by a volume-equivalent 
sphere, resulting in 
( )( )432DGF1)0( )(O))i(32()( kdkdkdbi ++= IM , (40)
611992.1)34( 31DGF1 ≈= πb . (41)
An exact evaluation, obtained without expanding the exponent, of Eq. (38) for the equivolume 
sphere with radius 31)43( πda =  was performed by Livensay and Chen [36] and 
implemented into the DGF formulation of the DDA by Hage and Greenberg [14,35] and later 
Lakhtakia [37]: 
[ ]1)iexp()i1()38()0( −−= kakai IM π . (42)
In terms of the first two orders of expansion, this yields an identical result as Eq. (40). Finally 
the polarizability is obtained as 
( )( )32DGF13CM
)i()32()(1 kdkdbd +−
α . (43)
We denote the method based on Eq. (42) as LAK. Differences between LAK and DGF should 
be noticeable only for large values of kd. 
Dungey and Bohren [38], using results by Doyle [39], proposed the following treatment 
of the polarizability. First, each cubic cell is replaced by the inscribed sphere that is called a 
dipolar subunit with a higher relative electric permittivity εs as determined by the Maxwell-
Garnett effective medium theory [27]: 
f , (44)
where 6π=f  is the volume filling factor. Other effective medium theories also may be used 
[40]. Next, the dipole moment of the equivalent sphere is determined using the Mie theory, 
and the polarizability is defined as [39] 
i= , (45)
where α1 is the electric dipole coefficient from the Mie theory (see e.g. [41]): 
)()()()(
)()()()(
sssssss
sssssss
1 xmxxxmm
xmxxxmm
α , (46)
where ψ, ξ are Riccati-Bessel functions; 2s kdx =  and ss ε=m  are the size parameter and 
the relative refractive index of the equivalent sphere. We denote this formulation for the 
polarizability as the a1-term method (note that this terminology was introduced later [42]). It 
has the particular property that 1constCMM ≠→αα  when , contrary to all other 
polarization prescription, for which this ratio approaches 1. It should be noted that the Mie 
theory is based on the assumption that the external electric field is a plane wave. In most 
applications of the DDA this is true for the incident electric field, but not for the field created 
by other subvolumes. Therefore the a
1-term method is expected to be correct only for very 
small cell size. Hence it is not clear whether this method has advantages even compared to 
CM. On the other hand, this method may be more justified for clusters of small spheres, 
where each sphere can be considered as a dipole (see Subsection 3.3). 
Draine and Goodman [3] pointed out that considering electric fields constant for 
evaluating integrals over a cell introduces errors of order ( )( )2O kd . This represents a problem 
for many polarizability corrections, based on integral equations. Draine and Goodman 
approached this problem from a different angle. They determined the optimal polarizability in 
the sense that an infinite lattice of point dipoles with such polarizability would lead to the 
same propagation of plane waves4 as in a medium with a given refractive index. This 
polarizability was called LDR (Lattice Dispersion Relation) and is, as expected, CM plus 
high-order corrections. These corrections in turn depend on the direction of propagation a and 
the polarization of the incident field e0: 
( ) ( )[ ]322LDR32LDR2LDR13CM
)i()32()(1 kdkdSmbmbbd +++−
α , (47)
8915316.1LDR1 ≈b , , , 1648469.0
2 −≈b 7700004.1
3 ≈b (48)
( )∑=
20eaS . (49)
We use a reverse sign convention in the denominator of Eq. (47) and the LDR coefficients as 
compared to the original paper [3]. 
Recently it has been shown [43] that the LDR derivation is not completely accurate, 
since the resulting dipole moment does not satisfy the transversality condition, for which a 
correction was proposed. This corrected LDR (CLDR) differs principally in the fact that the 
polarizability tensor can not be made isotropic but only diagonal [43], though not dependent 
on the incident polarization: 
( ) ( )[ ]3222LDR32LDR2LDR13CM
)i()32()(1 kdkdambmbbd +++−
α . (50)
Another flaw of LDR is that it is evidently not correct for dipoles near the particle surface. 
However, it is not clear how to evaluate the effect of these mistreated surface dipoles on the 
overall results, e.g. on the scattering cross section. 
Further improvement of the DDA was initiated by Peltoniemi [26] (PEL) who showed 
that the term M(Vi) in Eq. (7) can be evaluated exactly up to the third order of kd by 
expanding the term )()( rEr ′′χ  under the integral in a Taylor series over the point irr =′ , 
yielding 
( ) ( ,)(O3i3)iexp(d
)iexp(
EkdERRRRkRRk
ERkRRk
REMVM
ντρτρνμ
+∂∂−+−
∂−++=
                                                
 (51)
where χ, E and their derivatives are all considered at the point ri. Eq. (51) is correct up to the 
third order of kd since the third term in the Taylor series vanishes because of symmetry. For 
spherical Vi of radius a, the integrals can be evaluated exactly [26] in a way similar to 
4 with certain direction of propagation and polarization state. 
obtaining Eq. (42), but only terms of less than fourth order of kd are significant, which results 
( EkaakakaVi χχχχπ 42232 )(O)(10
)( +⎥
⎛ ⋅∇∇−∇−⎟
⎛ += EEEM ). (52)
If χ is constant inside the cell then the Maxwell equations state that 
EE 222 km−=∇ , 0=⋅∇ E . (53)
Hence Eq. (9) is valid up to the third order of ka and 
( )[ ]322 )(i)32()()101(1)34( kakami ++= IM π . (54)
Piller and Martin [44] proposed using sampling theory to evaluate the integrals in 
Eq. (1). The electric field and the susceptibility is sampled: 
∑ −′=′′
iiih )()()()()(
r rErrrrEr χχ , (55)
where hr(r) is the impulse response function of an antialiasing filter defined as 
)cos()sin(
qrqrqr
=r , (56)
where dq π2= . Eq. (1) is then transformed to Eq. (10) with the so-called filtered Green’s 
function, defined as 
∫ −′′′=
r3 )(),(d
ij hrV
rrrrGG . (57)
Eq. (57) can be viewed as a generalization of Eq. (13). The latter is obtained if a pulse 
function is considered instead of hr. The integral in Eq. (57) is evaluated analytically [44], 
taking V0 to be infinitesimally small. The filtered Green’s function does not have a singularity 
when , therefore ji rr = iiii V GM = . It was shown that the Fourier spectrum of E(r) lies on a 
sphere with radius m(r)k, if m is constant in the vicinity of r. Therefore at least two sampling 
points per wavelength in the scatterer are required. The susceptibility is also filtered, either by 
a mean value filter or a more complicated one, e.g. a Hanning window. This approach is 
called FCD (filtered coupled dipoles), and a computer code library for evaluation of filtered 
Green’s function is available [45]. 
Chaumet et al. [11] proposed direct integration of the Green’s tensor (IT) in Eqs. (12), 
(13). A Weyl expansion of the Green’s tensor is performed, transforming it to a form allowing 
efficient numerical computation of the self-term ( LM − ). They also proposed a correction to 
the second term in Draine’s expression for Cabs (Eq. (35)). Extension of their results to a non-
isotropic self-term is 
( ) ( )[ ]∑ −⋅+⋅=
iiiiiii VkC /)(ImIm4
***exc)0(
abs PLMPEPπ
The corrected second term is based on radiation energy of a finite dipole [11]: , in 
contrast to a point dipole used in the derivation of Eq. 
)Im( self ∗⋅ ii PE
(35). One can see that Eqs. (58) and 
(31) are equivalent. Moreover, both of them are equivalent to Eq. (35) if and only if 
IAM iii Vk
3i)32(+= , Hii AA = . (59)
This condition is similar, but not equivalent, to Eq. (36) and is always satisfied for RR, DGF, 
and LAK. Other polarizability prescriptions satisfy Eq. (59) for real m, then both Eqs. (58) 
and (35) result in zero absorption. 
Rahmani, Chaumet, and Bryant [46] proposed a new method (RCB) to determine 
polarizability based on the known solution of the electrostatic problem for the same scatterer. 
In the static limit the electric field at any point is linearly related to the incident field 
)()()( 01 rErCrE −= . (60)
Substituting Eq. (60) into Eq. (20) with the static Green’s tensor, one can obtain the 
polarizability, which would give an exact solution in the static limit, as 
1RCB −= iiii V Λχα , (61)
ii CCrrGCΛ
1),( −
∑+= χ , (62)
where )( ii rCC = . This static polarizability then replaces the CM polarizability, and the RR 
(Eq. (34)) is applied to it [46] to obtain the final polarizability for DDA simulations. It was 
later shown that RCB polarizabilities differ significantly from CM only for dipoles closer than 
2d to the interface [47]. 
In their next manuscript [48] Rahmani et al. stated that the previous derivation is correct 
only if the tensor C  is constant inside the particle (e.g. for ellipsoids), since otherwise the 
polarizability tensor obtained from Eq. (61) is generally not symmetric, which is physically 
impossible in the static case. This shows that a particle with a non-constant C  is not 
equivalent to any set of physical point dipoles even in the static regime. However, it is 
equivalent to a set of non-physical dipoles with an asymmetric polarizability. Therefore, the 
polarization defined by Eq. (61) formally can be used, by itself or with RR, even when C  is 
not constant. 
Collinge and Draine [47] empirically combined the RCB prescription with CLDR to get 
the surface-corrected LDR (SCLDR): 
( )( ) 13RCBRCBSCLDR −−= BαIαα d , (63)
where B  is the correction matrix (analogous to Eq. (50)): 
( )[ ]3222LDR32LDR2LDR1 )i()32()( kdkdambmbbB +++= μμνμν δ . (64)
All methods based on the paper by Rahmani et al. [46] are initially limited to very specific 
shapes of the scatterer (ellipsoids, infinite slabs and cylinders). Expansion of its applicability 
to other shapes is debatable [48] and would anyway require a preliminary solution of the 
electrostatic problem for the same shape, which is generally not trivial. 
All DDA formulations are schematically depicted in Fig. 1, which also shows 
interrelations between them. Some formulations can be compared unambiguously in terms of 
theoretical soundness: one is an improvement of the other, i.e. it employs fewer 
approximations. Such formulations are depicted in the same column on Fig. 1, while others 
cannot be compared directly with each other; they give rise to different columns. Comparison 
between formulations from different columns can and has been made almost exclusively 
empirically by comparing the accuracy of the simulation results (see Subsection 3.2). 
All the above techniques are aimed at reducing discretization errors; only a few aim at 
reducing shape errors. Some of them employ adaptive discretization (different dipole sizes) to 
better describe the shape of the scatterer (see Subsection 3.4). Another approach is to average 
susceptibility in boundary subvolumes. The simplest averaging using the Lorentz-Lorenz 
mixing rule was proposed by Evans and Stephens [49] for the case of the boundary between 
the scatterer and its surrounding medium 
3434 e
, (65)
where  is the effective susceptibility, and f is the volume fraction of the subvolume actually 
occupied by scatterer. 
A more advanced averaging, called the weighted discretization (WD), was proposed by 
Piller [13]. It modifies the susceptibility and self-term of the boundary subvolume.5 The 
particle surface, crossing the subvolume Vi, is assumed linear and divides the subvolume into 
two parts: the principal  that contains the center and a secondary  with susceptibilities piV
                                                 
5 any subvolume that has non-zero intersection with both the scatterer and the outer medium. All such 
subvolumes are accounted for. 
 Integral Eq. (1)  
discretization 
(no assumptions)
Eq. (7) 
General 
formulation of 
DDA – Eq. (20) 
Eqs. (8), (9) 
DGF, LAK 
Eqs. (11) 
Eq. (14)
(weak form)
CLDR 
a1 term 
SCLDR FCD 
sampling with 
antialiasing 
filter 
removing 
antialiasing 
filter 
improving polarizability starting 
from dipole formulation 
complies 
Eq. (14) 
simplifies to 
Fig. 1. Scheme of interrelation between the different DDA models discussed in Section 3.1. Arrows 
down correspond to assumptions employed. Vertical position of the method qualitatively corresponds 
to its accuracy (higher = better), however methods in different columns cannot be compared directly. 
iχ ,  and electric fields ,  respectively. Electric fields are considered constant 
inside each part and related to each other via a boundary condition tensor 
iχ ii EE ≡
iT : 
iii ETE =
s . (66)
Then the total polarization of the subvolume can be evaluated as follows: 
iiiiiiiii
i VVVr
EEErErP essspp3 )()(d χχχχ =+=′′′= ∫ , (67)
( ) iiiiiii VVV TI ssppe χχχ += . (68)
The susceptibility of the boundary subvolume is replaced by an effective one. 
The effective self-term is evaluated directly starting from Eq. (3), considering χ and E 
constant inside each part: 
( ) ( ) ii
rr TrrGrrGrrGrrGM ss3ps3ee
),(),(d),(),(d χχχ ∫∫ ′−′′+′−′′= . (69)
Piller [13] evaluated the integrals in Eq. (69) numerically. The final equations are the same as 
Eq. (20), where polarizabilities are obtained from Eq. (18) using effective susceptibilities and 
self-terms for boundary subvolumes. Hence, WD does not modify the general numerical 
scheme. 
Currently, there are no rigorous theoretical reasons for preferring one formulation over 
others. However, theoretical analyses of DDA convergence when refining discretization 
recently conducted by Yurkin et al. [32], showed that IT and WD significantly improve the 
convergence of shape and discretization errors, respectively. Experimental verification of 
these theoretical conclusions is still to be performed. 
Table 1. Accuracy of different DDA formulations for a sphere.a
Value Method x a/d y m Error, % Ref.
Cext a1-term 1÷2 2÷4
c 0.65 
0.85 
1.33+0.05i 
1.7+0.1i 
CSec, S11 LAK 9 
0.44 
0.42 
0.51 
1.05 
1.33+0.01i 
2.5+1.4i 
0.05, 37 
0.5, 35 
4, 15 
Csca, Cabs DGF ≤3.2c 16 ≤1 4+3i 5, 10÷30 [3]
CSec LDR ≤8c 16 ≤0.5, ≤0.1 m-1≤1 1, 2 
Csca 
LDR ≤7c 16 ≤1 
≤0.5, ≤1 
2+i 1.5 
3, 4 
CSec ≤16c
25 ≤1 1.6+0.0008i 
2.5+0.02i 
LDR [51]
[4]eCSec 
LDR any any ≤1 |m|≤3 5 
20÷30 
Csca LDR ≤10c 16 kd≤0.63 0.69 
0.41 
0.29 
[148]
S11 LDR ≤10c 24 kd≤0.42 0.69 
0.41 
Cext,  RMS11S
3.2 Accuracy of DDA simulations 
Over the years many results on the accuracy of DDA simulations have been published. It is, 
however, generally hard to systematically compare the relevant manuscripts because they all 
use different independent parameters, such as the size parameter x, refractive index m, or 
discretization, as a function of which the error is measured. We will describe discretization by 
the parameter kdmy =  or . The former is used wherever possible; however, 
in some cases a description of results is more straightforward in terms of y
kdmy )Re(Re =
Re. Accuracy results 
LDR 20÷160 
20÷130 
20÷60 
20÷30 
20÷30 
32÷256c
40÷256c
48÷128c
56÷80c
64÷88c
0.61÷0.65d
0.56÷0.64d
0.58÷0.65d
0.57÷0.60d
0.56÷0.62d
0.62d
1.05 
0.04, 38 
0.4, 23 
1, 59 
4.4, 56 
5.7, 105 
2.0, 86 
[113]
Ψ FCD π, 2π 2.8, 5.6c 1.7 1.5 1 [44]
WD-FCD 0.5÷3.2c
1.5÷3.8c
0.9÷1.5c
yRe=0.63 |m|<7b
|m|<2.5b
|m|<4b
[10]Ψ 
IT ≤5.2c CSec 
≤2.1c
≤1.1c
8 ≤1 1.5+0.3i 
3.5+1.4i 
7.1+0.7i 
Cabs 
CSec RCB-RR ≤8.2c 
≤7.5c 
≤5.9c 
≤3.4c 
≤1.3c
16 ≤1 1.8+0.4i 
1.9+i 
2.5+i 
2.5+4i 
7.4+9.4i 
CSec SCLDR 
SCLDR 
≤7.2c 
≤1.5c 
≤1.5c
12 ≤0.8 1.33+0.1i 
5+4i 
5+4i 
a All errors are relative. CSec denotes the maximum error over all cross sections, S11 and  correspond to 
maximum and root mean square error over the range of scattering angles, Ψ is the normalized mean error of the 
far-field electric fields [44]. In some cases two errors are shown in one cell separated by a coma. They 
correspond to two values of one of the parameters in the same row. 
b approximate description of the range. 
c this value is determined by other values in the same row. 
d this value is slightly different for different size parameters. 
e this corresponds to the “rule of thumb” for spheres. 
for scattering by a sphere are summarized in Table 1. All manuscripts on this subject can be 
divided into two classes: those that fix x and vary N (or equivalently, the number of dipoles 
per sphere radius a/d) with y, and those that fix a/d and vary the size parameter with y. The 
former is easier to interpret; the latter is easier to simulate. To facilitate comparison between 
different methods we provide both x and a/d, however one of them is dependent on the other. 
Some additional information on these results follows below. 
Draine and Goodman [3] compared RR, DGF, and LDR for cross sections of a sphere 
with . DGF is generally more accurate than RR. For 16/ =da 1|1| ≤−m  LDR gives superior 
or comparable results to DGF, for i2+=m  LDR and DGF are comparable, and for 
 DGF is preferable over LDR. In the review of LDR DDA, Draine and Flatau [4] 
summarized that for  cross sections can be evaluated to accuracies of a few percent 
provided . In that case differential cross sections have satisfactory accuracy: relative 
errors up to 20-30%, but only where the absolute value of the differential cross sections is 
small. For spheres, such results are obtained even for 
i34 +=m
2|| ≤m
3|| ≤m . Comparison of CLDR to LDR 
[43] only results in minor differences. Generally CLDR results in slightly better accuracy for 
Csca, but worse for Cabs. 
Piller and Martin [44] compared FCD to LAK by studying the dependence of the mean 
relative error of the far-field electric fields (Ψ) on y for spheres with π=x , 2π and . It 
was shown that FCD (with a Hanning window filter for the electric permittivity ε) is roughly 
3 times more accurate than LAK in the range 
5.1=m
5.27.0 ≤≤ y  and gives similar accuracies for 
 (for larger spheres). Comparison of WD to traditional methods [13] was performed 
for spheres with 
4.0≤y
π=x , 2π and 32.1=m , i7.01.2 + . LAK was used to determine 
polarizabilities. For  in the range 32.1=m 3.14.0 ≤≤ y  overall accuracy was only slightly 
improved, but error peaks for certain values of y were smoothed out. For i7.01.2 +=m  
accuracy was improved 4-5 times over the whole range 3.1≤y . Piller also showed [10] that a 
combination of WD and FCD gives even better results. Generally FCD decreases the negative 
effects of Re(ε) on accuracy and WD those of Im(ε). 
Rahmani et al. [48] showed that RCB was clearly superior to CM in calculating cross 
sections for fixed  and m from 16/ =da i4.08.1 +  to i4.94.7 +  in the range . Two 
corrections (LDR and RR) over the static case were compared, and they gave similar overall 
results. Improvement of overall accuracy compared with CM was 2-5 times in all cases 
studied. For a thin slab, it was shown [46,48] that the internal fields calculated using RCB 
differ from those by CM mostly near the interfaces, where RCB yields much smaller errors, 
almost the same as far from interfaces. 
Collinge and Draine [47] compared LDR, RCB, and SCLDR in calculations of cross 
sections of spheres with . It was shown that for 12/ =da i01.033.1 +=m , LDR and SCLDR 
are superior in the range , while for 8.0≤y i45 +=m , SCLDR and RCB are superior. 
Convergence of cross sections for spheres and ellipsoids for increasing N with fixed x and 
different m (from  to i01.033.1 + i45+ ) also was studied. SCLDR showed the most stable 
results for all cases, being the most or close to the most accurate one; however, for ellipsoids 
with large Im(m) RCB gave significantly more accurate results for Csca, especially for larger y. 
Performance of the DDA for more complex shapes also was studied by different 
authors. Flatau et al. [50] compared DDA simulations for a bisphere with an exact solution 
from a multipole expansion. For i01.033.1 +=m , 16/ =da , and 8.0≤y , LDR was several 
times more accurate than DGF and resulted in errors of less than 0.5% for both Csca and Cabs. 
Xu and Gustafson [51] made a similar but much more extended study of LDR. For 
, i008.06.1 +=m 25/ =da , and , errors in C4.0≤y ext, Cabs, and θcos  are within 10%. For 
, errors in the angular dependence of S81.0=y 11 are up to 20% while S12 and S21 were 
completely wrong. For , errors in cross sections exceed 10% for . i02.05.2 +=m 3.0≥y
Errors in the angular dependencies of the Mueller matrix elements are within 10-20% for 
 and increase rapidly with increasing y. For a fixed 3.0=y 3=x  and , errors 
i004.06.1 +=m
ext, Cabs, and >< θcos  decrease from 10% to 1% while y decrease from 1 to 0.2. For 
, the angular dependence of S33.0=y 11 is in good agreement with the rigorous solution, 
while S12 and S21 differ significantly for certain orientations of the bisphere. 
Hage and Greenberg [14] compared LAK to experimental results obtained from 
microwave experiments on porous cubes. Using i005.0362.1 +=m , 64.0=y  and , 
they obtained a difference of less than 40% with the experimental results of angular scattering 
patterns, except for deep minima. Light scattering of cubes, tiles, and cylinders with similar 
parameters also was studied and comparable differences between experiment and theory were 
obtained. Theoretical errors were estimated to be less than 10%, except for deep minima. 
5504=N
Iskander et al. [34] conducted a limited test of LAK for small elongated spheroids, 
comparing the results to those obtained using an iterative extended boundary condition 
method. Using , calculations were performed for aspect ratios up to 20 with maximum 
size parameter of the long axis being 10 and 0.5 for 
i01.033.1 +=m  and  
respectively. Errors in scattering cross section were 21% and 11%, respectively. Ku [52] 
compared LAK with CM and the a
i28.076.1 +
1-term for different shapes, but his conclusions are based 
on a large parameter y (up to 2), and are therefore suspicious and not further discussed here. 
Andersen et al. [53] studied the performance of the DDA for Rayleigh-sized clusters of 
a few spheres (most DDA formulations are then equivalent to CM). Several constituent 
materials were tested, all with high refractive indices in the studied region. It was shown that 
the DDA failed to converge using the fixed computational resources for very high (up to 13.0) 
and very low (down to 0.12) Re(m); up to 30 dipoles were used per diameter of a single 
sphere. 
It can be concluded that particles with more complex shapes than spheres are more 
difficult to model with the DDA, leading to larger errors for the same m and y. This effect can 
be explained in general by the increase of surface to volume ratio and hence larger fraction of 
boundary subvolumes [32]. Another possible reason is complex regions, e.g. contact between 
two particles in a cluster, where rapid variation of the electric field deteriorates the overall 
accuracy. There is, however, a notable exception from this general tendency. Shapes, which 
can be modeled exactly by a set of cubical dipoles, e.g. a cube, can be simulated using the 
DDA much more accurately than spheres, especially for small y [32]. 
Draine and Flatau [4] have introduced a “rule of thumb” for discretization: use 10 
dipoles per wavelength in the medium (i.e. either y or yRe equal to 0.63, depending on the 
interpretation). Though it is widely used, the accuracy of the results, when using such 
discretization, is hard to deduce a priori. Draine and Flatau themselves derived an estimate of 
the error based on a set of test simulations. This estimate is described above and mentioned in 
Table 1; it is usually cited as a “few percent accuracy in cross sections.” However, it may 
significantly over- or under-estimate the error, especially for large size parameters. Moreover, 
it does not completely account for the dependence on m, even in the stated range of its 
application ( ), since DDA accuracy deteriorates rapidly with increasing m (see 2|| ≤m Table 
1). Still, the rule of thumb is good first guess for many applications. 
Most studies of DDA accuracy are limited to integral scattering quantities and, at most, 
the angular dependence of S11. In only a few manuscripts are other scattering quantities 
studied. For instance, Singham [54] simulated the angular dependence of Mueller matrix 
element S34 for spheres and less compact particles, using CM polarizability. It was shown that 
an accurate simulation of this element requires smaller values of y than for S11. For 55.1=x  
and  a calculation of S33.1=m 11 was accurate already for 8.0=y , while  was 
required for S
2.0≤y
34. It was also reported that for less compact objects like discs and rods, the 
required y was larger, 0.4 and 0.55 respectively, because of the smaller interaction between 
the dipoles. However, Hoekstra and Sloot argued [55] that this effect is mostly caused by the 
pronounced S34 sensitivity to surface roughness, which is significant for smaller size if y is 
fixed. They showed that for  and 7.10=x 05.1=m , very high accuracy is achieved with 
 because of the larger number of dipoles used. 66.0=y
Internal fields are an intermediate result in the DDA. They cannot be directly compared 
to the experimental results; however, all measured scattering quantities are derived from 
them. Therefore, a study of their accuracy can reveal greater understanding of the nature of 
DDA errors. Hoekstra et al. [56] performed such a study for LAK polarizability. Three 
spheres were examined with , 9, 5 and 9=x 05.1=m , i01.033.1 + ,  respectively. 
Values of y were 0.44, 0.42, and 0.51 respectively. The most significant errors in the 
amplitude of the internal field were localized at the boundary of the spheres with maximum 
relative errors of 3.4%, 19%, and 120% respectively. Errors in S
i4.15.2 +
12, S33, S34 were significant 
only for the third sphere. It was shown that for a given yRe these errors rapidly increase with m 
but only slightly depend upon x in the range from 1 to 10. Moreover, the DDA is capable of 
reproducing resonances of Mie theory, although their positions are slightly shifted (less than 
1% in m). 
Druger and Bronk [57] studied the accuracy of the internal fields for single and coated 
spheres. They used 5.1=x , , and CM polarizability. Errors in the internal fields were 
localized at the interfaces, with average errors larger than 30% for a single sphere with 
 and , and less than 7% for a single and concentric sphere with  and 
. The core of the concentric sphere has 
8.1≤m
8.1=m 17.0=y 3.1=m
08.0=y 1.1=m  and its diameter is half the total 
diameter. The angular dependence of the absolute values of S1 and S2 had significant errors in 
the side- and backscattering. It can be concluded that shape errors contribute mostly to the 
internal fields near the boundary, and increase with m. 
All the literature discussing DDA accuracy shows errors as a function of input 
parameters and discretization, which is the most straightforward way. The only exception so 
far is the rule of thumb, which is too general and approximate to be applied in many particular 
cases. A more useful way to present errors is to fix the desired accuracy for certain input 
parameters and find the discretization that results in such accuracy. Such an analysis can be 
applied directly to practical calculations and can be used to derive rigorous estimates of DDA 
computational requirements [58]. 
In a number of manuscripts the origin of errors in the DDA was examined to try to 
separate and compare shape and discretization errors [49,59-62]; however, no definite 
conclusions were reached. The uncertainty was due to the indirect methods used that have 
inherent interpretation problems. Recently, Yurkin et al. [63] proposed a direct method to 
separate shape and discretization errors, which can be used to study their fundamental 
properties. This method also can be applied to study the performance of different formulations 
aimed at decreasing shape errors, e.g. WD. For example, it has been shown that the maximum 
errors of S11(θ) for a sphere with  and 5=x 5.1=m , discretized using 16 dipoles per diameter 
( ), are mostly due to shape errors. However the same is not true for all measured 
quantities. In another manuscript [32] it was suggested that the discretization error should 
decrease more rapidly with decreasing y than shape errors. However, it is still hard to deduce 
a priori the importance of shape errors for a certain scatterer and y; hence, further systematic 
quantitative study is required. 
93.0=y
3.3 The DDA for clusters of spheres 
There are two main peculiarities when the DDA is applied to clusters of spheres. First, such 
particles are generally less compact, yielding smaller interactions between dipoles. This leads 
to a smaller condition number of the DDA interaction matrix and hence faster convergence of 
the iterative solver (see Section 4.1). Second, when the constituent spheres are small 
compared to the wavelength, each sphere can be modeled as one spherical subvolume, 
yielding some theoretical simplifications. 
A general theory exists [64] based on the Mie theory (generalized multiparticle Mie 
solution (GMM) [65]) that allows for highly accurate simulations of clusters of spheres. 
However, when many small spheres are used one wants to minimize the number of unknowns 
in the linear system. Direct reduction of the GMM to the lowest order (using only the first 
order expansion coefficients) leads to DDA + CM [64]. Improving accuracy in the GMM is 
done by accounting for higher multipole moments, while the DDA introduces higher order 
corrections to the coefficients of the linear system. It is not clear how the accuracy of these 
two methods compare with each other; however, the former should lead to a formulation 
similar to a coupled multipole method (Subsection 3.4) with a larger number of unknowns. 
DDA-based methods (starting usually with the integral equations introduced in Section 2) 
should be successful in making the formulation more accurate without increasing the number 
of unknowns, which is the goal for large clusters of small spheres. Moreover, the DDA may 
employ fast algorithms for solving the linear system. In this setting, the fast multipole method 
(FMM) (see Subsection 4.5) seems most promising. 
It should be noted, however, that a cluster having a small size parameter (i.e. in the 
electrostatic approximation) does not imply that all expansion coefficients, except the first 
one, are negligible. This is because the size of the constituent particles is also very small and 
the fields inside them are far from constant, especially when the spheres are located close to 
each other and have large refractive indices [66]. Therefore, the DDA does have some 
principal difficulties of calculating scattering by clusters of spheres. Mackowski [67], for 
instance, found that for some systems composed of spheres much smaller than the 
wavelength, up to 10 expansion terms were necessary to achieve convergence. In studies of 
osculating spheres, Ngo et al. [68] proved that the GMM could be chaotic and were able to 
calculate Lyapunov exponents, and that the slow convergence for the touching spheres was 
the result of the system lying in an attractor region. A recent paper by Markel et al. [69] 
presented computationally efficient modifications of the GMM in the static limit and 
demonstrated the insufficiency of the DDA to compute scattering properties of fractal 
aggregates accurately. However, Kim et al. [70] showed that the DDA is satisfactory in 
calculating the static polarizability of dielectric nanoclusters, especially of clusters with a 
large number of constituents. 
The development of DDA-based methods for calculating light scattering by clusters of 
small spheres was started by Jones [71,72], who developed a method similar to CM. Iskander 
et al. [34] used a method equivalent to LAK to calculate scattering of chained aerosol 
clusters. This subject was further investigated by Kosaza [73,74]. Lou and Charalampopoulos 
[75] (LC) further improved the calculations of the interaction term and scattering quantities. 
Starting from an integral equation for the internal field equivalent to Eq. (1), they assumed 
Eq. (11). After that the integrals in Eqs. (12) and (13) over spherical subvolumes can be 
evaluated analytically. The result for the interaction term is the following: 
),()()0( jiij ka rrGG η= , (70)
where a correction function η is defined as 
)O()101(1
cossin
3)( 42
x +−=
=η . (71)
Eq. (30) also is evaluated analytically, yielding 
∑ ⋅−−−=
ii knnkak )iexp()ˆˆ)((i)(
3)0( nrPInF η , (72)
( )∑ ∗⋅=
iikakC
inc)0(
ext Im)(4 EPηπ . (73)
The following expression for Cabs is stated without derivation: 
iikakC )Im()(4
abs EPηπ . (74)
Markel et al. [76] applied the DDA to fractal clusters of spheres, and studied their 
optical properties. However, they have not fixed the polarizability of a single dipole but rather 
treated it as a variable, calculating the dependence of a cluster’s optical characteristics upon it. 
Pustovit et al. [77] argued that the DDA is inaccurate for touching spheres. They developed a 
hybrid of the DDA and the GMM, which considers only pair interactions between spheres (as 
the DDA) but, when calculating them, accounts for higher multipole terms. This formulation 
can be considered as the one providing a more accurate evaluation of the interaction term 
(Eq. (13)), and hence similar to LC. 
LC was compared to DGF and LAK in a Csca computation of a cluster of 10 particles for 
 and . Differences between DGF and LAK are less than 1% (as 
expected), while the difference between LC and LAK increases quadratically with ka, 
reaching 10% for 
i7.07.1 +=m 5.005.0 ≤≤ ka
5.0=ka . However, as no exact (e.g. GMM) solution is presented, the 
accuracy of each individual method is not clear. 
Okamoto [42] tested the a1-term method for clusters of up to 3 touching spheres. No 
effective medium is needed in this case, making the method sounder. It was shown that the a1-
term is clearly superior to LDR in cross-sections calculations, when each sphere is treated as a 
single dipole. Errors of the a1-term are less than 10% for 2.1≤y  when . For 
three collinear touching spheres the errors are 30% and 40% for  and 2.8 when 
 and  respectively. However, errors do not seem to diminish 
significantly for small y (results are presented only down to 
i01.033.1 +=m
9.1≤y
i01.033.1 +=m i2 +
2.0=y ). Therefore, the a1-term 
seems suitable for obtaining quick crude estimations of cross sections. 
In the sequel of this subsection we mention several applications of the DDA to 
scattering from clusters of spheres. It was applied to describe the scattering by astrophysical 
dust aggregates [78,79] using the a1-term method. Hull et al. [80] applied CM DDA to Diesel 
soot particles. LC was applied [81] to the computation of light scattering by randomly 
branched chain aggregates. Lumme and Rahola [40] studied scattering properties of clusters 
of large spheres (each modeled by a set of dipoles) with the a1-term method considering 
astrophysical applications. Hage and Greenberg [35] studied scattering by porous particles, 
which were modeled as clusters of cubical cells making their method equivalent to standard 
LAK. Recently the DDA with LDR was used [82] to model scattering by porous dust grains 
and compare them to approximate theories, e.g. effective medium theories. It also was used to 
study light scattering by fractal aggregates [83], especially its dependence on the internal 
structure [84]. 
3.4 Modifications and extensions of the DDA 
Bourrely et al. [85] proposed to use small d to minimize surface roughness, but larger dipoles 
inside the particle. Starting with small dipoles with CM polarizability, one dipole is combined 
with 6 adjacent ones (if they all have the same polarizability) producing a dipole, located at 
the same point but with a 7 times larger polarizability. This operation is repeated while 
possible. Interaction terms are considered in their simplest form (Eq. (14)). This method 
allows the decrease of the shape errors with only a minor increase in the number of dipoles. 
The authors showed that this method is more than two times more accurate than CM for some 
test cases. 
Rouleau and Martin [86] proposed a generalized semi-analytical method. A dynamic 
grid is used to evaluate the integral in Eq. (1). First, a static grid is built inside the particle. 
Then each point on the static grid is used as an origin of a spherical coordinate system, and 
the particle is approximated by an ensemble of volume elements in these spherical 
coordinates. As usual, the polarization inside each subvolume is assumed constant, but 
Eq. (13) can be evaluated analytically in spherical coordinates. Polarization inside a 
subvolume is obtained by interpolation of its values at the points of the static grid. In addition, 
adaptive gridding is employed, where smaller subvolumes are used at the boundary of the 
particle. 
Mulholland et al. [87] proposed a coupled electric and magnetic dipole method 
(CEMD), where a magnetic dipole is considered at each subvolume together with an electric 
dipole. Polarizabilities are derived from the a1 and b1 terms of the Mie theory. CEMD 
requires two times more variables in the linear system, since the electric and magnetic fields 
are interconnected. Lemaire [88] went further and developed the coupled multipole method, 
considering also the electric quadrupole. Addition of the electric quadrupole can be 
considered as a more accurate evaluation of the interaction term in Eq. (13), as compared to 
Eq. (14). It results in even better accuracy than CEMD, but at the expense of additional 
computation time. The major disadvantage of all these four methods is that the matrix of the 
system of linear equations does not seem to have any special form, suitable for faster 
algorithms (see Section 4). Therefore computational costs are much larger compared to 
regular methods, thus limiting their practical use. In what follows, several DDA extensions 
are mentioned without further discussion. 
The theoretical basis for application of the DDA to optically anisotropic particles was 
summarized by Lakhtakia [89]. Loiko and Molochko [90] applied the DDA to study light 
scattering by liquid-crystal spherical droplets. Smith and Stokes [91] used the DDA to 
calculate the Faraday effect for nanoparticles. Researchers in the electrical engineering 
community applied MoM (in a variation that is equivalent to the DDA) to anisotropic 
scatterers [92,93]. 
Rectangular parallelepipeds can be used as subvolumes in the DDA [11,23,43]. This 
allows an accurate description of light scattering by particles with large aspect ratios, using 
fewer dipoles and is also compatible with FFT techniques (Subsection 4.4). 
Khlebtsov [94] proposed a simplification of the DDA, based on the assumption that all 
polarizations are parallel to the incident electric field. The number of variables is thus reduced 
three times, however at a cost of accuracy. Moreover, depolarization is completely ignored. 
Markel [95] analytically solved the DDA equations for scattering by an infinite one-
dimensional periodic dipole array. This approach is similar to the one used in obtaining the 
LDR formulation for dipole polarizability [3]. 
Chaumet et al. [96] generalized the DDA to periodic structures, and further to defects in 
a periodic grating on a surface [97]. The idea of using the complex Green’s tensor in the 
standard DDA formulation was summarized by Martin [98]. 
Yang et al. [99] used the DDA to calculate surface electromagnetic fields and determine 
Raman intensities for small metal particles of arbitrary shape.  
Lemaire and Bassrei [100] showed that the shape of an object can be reconstructed from 
the measured angle dependence of scattered intensities. This procedure can be thought of as 
an inversion of the dependence between dipole polarizabilities and scattering. This 
dependence is taken from the DDA. A similar idea is used in recent manuscripts on optical 
tomography [101-103]. 
Zubko et al. [104] modified the Green’s tensor used in the DDA to study the 
backscattering of debris particles. They showed that the far-field part of the Green’s tensor is 
responsible for both the backscattering brightness surge and the negative polarization branch. 
4 Numerical considerations 
In this section the numerical aspects of the DDA are discussed. One should keep in mind, 
however, that final simulation times depend not only on the chosen numerical methods but 
also on the particular implementation. Recently, Penttila et al. [105] have compared four 
different computer programs for the DDA. These are based on almost identical numerical 
methods: the Krylov-subspace iterative method (Section 4.1) combined with a FFT 
acceleration of the matrix-vector product (Section 4.4). However, simulation times may differ 
by several factors. Optimizations of computer codes are not further discussed in this review. 
4.1 Direct vs. iterative methods 
There are two general types of methods to solve linear systems of equations , where x 
is an unknown vector and A and y are known matrix and vector, respectively: direct and 
yAx =
iterative [106]. Direct methods give results in a fixed number of steps, while the number of 
iterations required in iterative methods is generally not known a priori. The most usual 
example of a direct method is LU decomposition, which allows quick solving for multiple y 
once the decomposition is performed. Iterative methods are usually faster, less memory 
consuming and numerically more stable. However, iterative methods cannot be considered 
superior over direct, since they strongly depend on the problem to solve [107]. 
For a general n×n matrix (in DDA Nn 3= ) computation time of LU decomposition is 
O(n3) and storage requirements O(n2), while computation time for one iteration is O(n2) [107]. 
Iterative methods for a general matrix converge in O(n) iterations, although some of them 
may not converge at all. However, in many cases satisfactory accuracy can be obtained after a 
much smaller number of iterations. In these cases, iterative methods can provide significant 
increases in speed, especially for large n. Most iterative methods access the matrix A only 
through matrix-vector multiplication (sometimes also with the transposed matrix), which 
allows the construction of special routines for calculation of these products. Such routines 
may decrease memory requirements, since it is no longer necessary to store the entire matrix, 
especially for matrices of special form (see Subsection 4.3). A special structure of the matrix 
may also allow acceleration of the matrix-vector product from O(n2) to O(nlnn) (see 
subsections 4.4, 4.5). However, the same applies to direct methods (see Subsection 4.3). 
Throughout DDA history, mostly iterative methods were employed (however see 
Subsection 4.6). At first, they were used to accelerate computations [1], but they also allowed 
larger numbers of dipoles to be simulated [6,108], since storage of the entire matrix is 
prohibitive for direct methods. The most widely used iterative methods in the DDA are 
Krylov-space methods, such as [107] conjugate gradient (CG), CG applied to the Normalized 
equation with minimization of Residual norm (CGNR), Bi-CG, Bi-CG stabilized (Bi-
CGSTAB), CG squared (CGS), generalized minimal residual (GMRES), quasi-minimal 
residual (QMR), transpose free QMR (TFQMR), and generalized product-type methods based 
on Bi-CG (GPBi-CG) [109]. 
An important part of the iterative solver is preconditioning, which effectively decreases 
the condition number of the matrix A and therefore speeds up convergence. However, this 
requires additional computational time during both initialization and each iteration. 
Preconditioning of the initial system can be summarized as [107] 
yMxMAMM 12
21 )( =
− , (75)
where M1 and M2 are left and right preconditioners, respectively. Preconditioners should 
either allow fast inversion or be integrated into the iteration process. The simplest 
preconditioner of the first type is the Jacobi (point), which is just the diagonal part of matrix 
A. An example of the second type of preconditioner is the Neumann polynomial 
preconditioner of order l: 
)( AIM . (76)
QMR and Bi-CG can be made to employ the complex symmetric (CS) property of the 
DDA interaction matrix to halve the number of matrix-vector multiplications [110] (and thus 
computational time). Lumme and Rahola [40] were the first to apply QMR(CS) to the DDA 
and compared it with CGNR. They used m from i1.06.1 +  to i43+ , and x from 1.3 to 13.5, 
corresponding to N from 136 to 20336. For all cases studied QMR(CS) was 2-4 times faster 
than CGNR. 
Rahola [9] further studied QMR(CS) and compared it to CGNR, Bi-CG(CS), Bi-
CGSTAB, CGS, GMRES (full and with different memory length). For a “typical small 
problem” (parameters were not specified, unfortunately) the convergence of different methods 
was tested and QMR(CS) along with Bi-CG(CS) showed the best results. Although full 
GMRES was able to converge in fewer iterations, GMRES with as much as 40 memory 
lengths was slower than QMR(CS). 
Flatau [111] reviewed the use of iterative algorithms in the DDA and tested many of 
them, together with several preconditioners. He calculated scattering of a homogenous sphere 
with  and m from 1.33 up to 1.0=x i0001.05+ , 1=x  and m from 1.33 up to  and 
. Left (L) and right (R) Jacobi-, and first-order Neumann polynomial 
preconditioners were tested. Unfortunately the number of dipoles N was not specified, which 
hampers comparison with other studies. For small particles CG(L) was superior for all 
refractive indices studied. CG and CG(R) showed similar results, while CGNR(L) and Bi-
CGSTAB(L) were about 4 times slower. For 
i33.1 +
i0001.03+
1=x  Bi-CGSTAB(L) was superior while Bi-
CGSTAB,(R) and CGS,(L),(R) were slightly worse. TFQMR (both with and without Jacobi 
preconditioner) was 3-4 times slower. The first-order Neumann preconditioner showed 
unsatisfactory results. It was concluded that Bi-CGSTAB(L) is the most satisfactory choice 
for the DDA, and that method is the default one used in the DDSCAT program [6]. 
Recently Fan et al. [112] have compared GMRES, QMR(CS), Bi-CGSTAB, GPBi-CG, 
and Bi-CG(CS). They tested them on wavelength-sized scatterers (x up to 10) with m up to 
, and concluded that GMRES with memory depth 30 was the fastest, although it 
required four times more memory than the other methods. However, only the times of the 
matrix-vector product was compared, while other parts of the iteration may also take 
significant time, especially for GMRES(30). Choosing from less memory-consuming 
methods, QMR(CS) and Bi-CG(CS) showed a better convergence rate than Bi-CGSTAB and 
GPBi-CG, especially when 
i2.05.4 +
2>m . Moreover, the authors pointed out some flaws in the 
comparison by Flatau [111], making his conclusions insufficient. 
Yurkin et al. [113] employed QMR(CS), Bi-CG(CS), and Bi-CGSTAB to simulate light 
scattering by spheres with x up to 160 and 40 for 05.1=m  and 2, respectively. It was shown 
that convergence of the iterative methods becomes very slow with increasing x and m (up to 
105 iterations are required), and none of them is clearly preferable to the others. Moreover, 
there seems to be no systematic dependency of the choice of the best iterative solver on x and 
m; however, the difference in computational time was less than a factor of two, except for the 
largest x and m studied. 
Rahola [114] showed that the spectrum of the integral scattering operator for any 
homogenous scatterer is a line in the complex plane going from 1 to m2, except for a small 
amount of points, which corresponds to refractive indices that cause resonances for the 
specific shape. The spectrum of A is similar, since this matrix is obtained in the DDA by 
discretization of the integral operator (see also [9]). Assuming that the spectrum of A exactly 
lies on the specified line, it was shown that an estimate for the optimal reduction factor6 γ can 
be given as 
Eq. (77) is an approximation valid for small particle sizes, where no, or only few, resonances 
are present. However, in all cases the spectrum of A resembles the spectrum of the linear 
operator, which is defined by shape, size and refractive index of the scatterer. Therefore, the 
spectrum, and thus convergence, should not depend significantly on the discretization. This 
fact was confirmed empirically in other manuscripts [9,63]. 
Budko and Samokhin [115] generalized Rahola’s results to arbitrary inhomogeneous 
and anisotropic scatterers. They described a region in the complex plane that contains the 
whole spectrum of the integral scattering operator. This region depends only on the values of 
m inside the scatterer and does not depend on x. They showed that for purely real m or for m 
with very small imaginary part this region may come close to the origin, therefore the 
spectrum may contain very small eigenvalues for particles larger than the wavelength. This 
                                                 
6 Norm of the residual is decreased by this factor every iteration. 
may explain the extremely slow convergence of the iterative solver for real m and large x, 
which was recently obtained in numerical simulations [113]. Based on the analysis of the 
spectrum of the integral scattering operator for particles much smaller than the wavelength, 
Budko et al. [116] proposed an efficient iteration method for this particular case. 
It can be concluded that there are several modern iterative methods (QMR(CS), Bi-
CG(CS), and Bi-CGSTAB) that have proved to be efficient when applied to the DDA. 
However, none of them can be claimed superior to the others, and one should test them for 
particular light-scattering problems. Moreover, except for the simplest cases, preconditioning 
of the DDA interaction matrix is almost not studied, while there is a need for it for large x and 
m, since then all methods converge extremely slowly or even diverge. It seems to us that the 
next major numerical advance in the DDA will be achieved by developing a powerful 
preconditioner for the DDA matrix. 
A large number of dipoles requires large computational power and, hence, parallel 
computers are commonly used, e.g. [108,113]. Parallel efficiency is not discussed here, but 
for iterative solvers, it is generally close to 1 [117]. However, this is not true for all 
preconditioners [107], and hence heavy preconditioners requiring large computational time in 
combination with a parallel DDA implementation should be employed with caution. 
4.2 Scattering order formulation 
The Rayleigh-Debye-Gans (RDG) approximation [27] consists in considering E(r) equal to 
Einc(r). F(n) is then obtained directly from Eq. (24). Generalization of the RDG approach is 
obtained by iteratively solving the integral equation (1), which can be rewritten as 
)()()( inc rΛErErE += , (78)
where Λ is a linear integral operator describing the scatterer. The iterative scheme is readily 
obtained by inserting the current (l-th) iteration of the electric field E(l)(r) into the right side of 
Eq. (78) and calculating the next iteration in the left side: 
)()()( )(inc)1( rΛErErE ll +=+ . (79)
The starting value is taken the same as in RDG, , and the general formula for 
the solution is the following: 
)()( inc)0( rErE =
inc )()(
l rEΛrE , (80)
which is a direct implementation of the well-known Neumann series: 
lΛΛI , (81)
where I is the unitary operator. A necessary and sufficient condition for Neumann-series 
convergence is 
1<Λ . (82)
Physical sense of this iterative method lies in successive calculations of interaction 
between different parts of the scatterer. The zeroth approximation (or RDG) accounts for no 
interaction; the first approximation considers the influence of scattering of each dipole on the 
others once, and so on. Eq. (82) states that the interaction inside the scatterer should be small, 
but not as small as required for the applicability of RDG ( 1<<Λ ). In scattering problems, 
especially in quantum physics, Eq. (80) is called the Born expansion. 
Although theoretically clear, the Born expansion is not directly applicable [118], since 
each successive iteration requires analytical evaluation of multidimensional integrals with 
rising complexity, which quickly becomes unfeasible even for the simplest scatterers. The 
latest result is probably that of Acquista [118], who evaluated the Born expansion for a 
homogenous sphere up to second order. Therefore, realistic application of the Born expansion 
does require discretization of the integral operator, which is naturally done in the DDA. 
A scattering order formulation (SOF) of the DDA was developed independently by 
Chiappetta [119] and Singham and Bohren [12,120] by applying the Neumann series to 
Eq. (17). Λ is then a matrix defined as jijij αGΛ = , where each element is a dyadic, which 
can be expressed as a 3×3 matrix. An explicit check of Eq. (82) for a certain scatterer is not 
feasible numerically, however de Hoop [121] derived a sufficient condition for scalar waves: 
1)(max)(2 20 <r
χπ kR , (83)
where R0 is the radius of the smallest sphere circumscribing the scatterer. Although not 
directly applicable to light scattering, Eq. (83) can be used as an estimate. 
The range of size parameter and refractive index where SOF converges is limited [120]. 
Moreover, even when SOF converges, more advanced iterative methods converge faster (see 
Subsection 4.1). However, SOF has clear physical sense and can be used to study the 
importance of multiple scattering. 
4.3 Block-Toeplitz 
A square matrix A is called Toeplitz if jiij aA −= , i.e. matrix elements on any line parallel to 
the main diagonal are the same [106]. In a block-Toeplitz (BT) matrix (of order K) elements 
ai are not numbers, but square matrices themselves: 
. (84)
A 2-level BT matrix has BT matrices as components ai. Proceeding recursively a multilevel 
BT (MBT) matrix for any number of levels is defined. 
Let us consider a rectangular lattice nx×ny×nz, numbered in the following way 
zzyyzxxzy ininninnni +−+−= )1()1( , (85)
where  indicates the position of the element along the axes. Let us also define 
the vector index . Then one can verify that the interaction matrix in Eq. 
},...,1{ μμ ni ∈
),,( zyx iii=i (20), 
defined by Eq. (13), satisfies the following: 
jiGGG −′== jiij . (86)
This equation alone can be used to greatly reduce the storage requirements of iterative 
methods by use of indirect addressing. Further improvement is to note that Eq. (86) defines a 
symmetric 3-level BT matrix (orders of subsequent levels – nx, ny, nz) whose smallest blocks 
are 3×3 matrices (dyadics) ijG . 
A rectangular lattice is not much of a restriction, since any scatterer can be embedded in 
an appropriate rectangular grid. However, additional “empty” dipoles should be introduced to 
build up the grid up to the full parallelepiped. Moreover, position and size of the dipoles 
cannot be chosen arbitrarily to better describe the shape of the scatterer. This is especially 
problematic for highly porous particles or clusters of particles, where the monomer has a size 
comparable to a single dipole. For all other cases these restrictions are minor compared to the 
large increase in computational speed, imposed by the BT-structure of the interaction matrix. 
A matrix-vector multiplication can be transformed to a convolution, which is computed using 
a fast Fourier transform (FFT) technique in O(nln(n)) operations (see Subsection 4.4). Note 
however, that alternative techniques exist that do not require a regular grid (see Subsection 
4.5). 
The BT-structure also permits acceleration of direct methods. Flatau et al. [122] used an 
algorithm for inversion of symmetric BT-matrices. It has complexity )(O 3 xnn  and storage 
requirements )(O 2 xnn , since only 2 block columns of the inverse matrix need to be stored. 
In this case the x-axis is oriented along the longest particle dimension. Recently Flatau [123] 
studied the special case of 1D DDA where all dipoles are located on a straight line and 
equally spaced, in which systems of equations for different components can be separated. The 
interaction matrix for each component is symmetric Toeplitz, and a modern fast algorithm can 
be applied for its inversion. This method requires preliminarily solving linear equations for 
two right sides (e.g., by some iterative technique); then multiplication of the inverse matrix by 
any vector (i.e., a solution of the linear system for any right part) requires only O(nln(n)) 
operations. However, Flatau pointed out a strict limitation for all methods for fast calculation 
of the inverse of the interaction matrix: they are applicable only when polarizabilities of all 
dipoles are the same, since otherwise the first term on the right side of Eq. (20) ruins the BT 
structure on the diagonal of the interaction matrix. Therefore, they are currently limited to 
homogenous rectangular scatterers. Fortunately, it is not a problem for matrix-vector 
multiplication, since the diagonal term can be evaluated independently and added to the final 
result. 
4.4 FFT 
Goodman et al. [124] showed that multiplication of the interaction matrix for a rectangular 
lattice (see Subsection 4.3) by a vector can be transformed into a discrete convolution 
jii PGPGPGy ′′=′== ∑∑∑
)2,2,2(
)1,1,1(
)1,1,1(1
zyxzyx nnnnnnN
jij , (87)
where iG′  is defined by Eq. (86) (and 0=′0G ) for μμ ni ≤ and 
⎧ ≤≤∀
otherwise,0
1:, μμμ njj
P . (88)
Both G ′  and  are then regarded as periodic in each dimension μ with period 2nP′ μ. A 
discrete convolution can be transformed with a FFT to an element-wise product of two 
vectors, which is easily computed. It requires evaluation of a direct and inverse FFT for each 
matrix-vector product. Each of them is a 3D FFT of order 2nx×2ny×2nz. This operation is done 
for each of the 3 Cartesian components of P′  and preliminary calculations is performed for 6 
independent tensor components of G ′ . 
A slightly different method can be devised based on the paper by Barrowes et al. [125], 
who developed an algorithm for multiplication of any MBT by a vector. The multiplication is 
brought down to a 1D convolution that is evaluated by two 1D FFTs of order 
)12)(12)(12( −−− zyx nnn . Flatau [123] proposed an algorithm of matrix-vector 
multiplication for BT interaction matrix (e.g. 1D DDA), which requires twice as many FFTs 
as the standard algorithm, but of order n instead of 2n. Although Flatau stated that an 
extension of this algorithm to the general 3D case is straightforward, it is at least not trivial 
and probably its complexity will scale the same as standard methods. 
4.5 Fast multipole method 
The fast multipole method (FMM) was developed by Greengard and Rokhlin [126] for 
efficient evaluation of the potential and force fields in N-body simulations where all pairwise 
interactions of N particles are computed. The FMM is based on truncated potential expansions 
[127]. It is also called a hierarchical tree method because particles are grouped together in a 
hierarchical way, and the interaction between single particles and this hierarchy of particle 
groups is calculated [128]. However, some researchers distinguish between single- and 
multilevel FMM [129,130]; only the latter is truly hierarchical. The FMM naturally fits the 
DDA, since the matrix-vector multiplication is actually computing the total field on each 
single dipole due to all other dipoles, as was noted by Hoekstra and Sloot [128]. The 
computational complexity of the FMM (see below) is similar to FFT-based methods (see 
Subsection 4.4), but it does not require any regularity of the grid, thus making it applicable to 
any scatterer. The drawback is that the FMM is conceptually more complex, making it much 
harder to code. Nonetheless, the FMM was implemented in the DDA by Rahola [9,127]. 
Error analysis is critical for the FMM, since the acceleration is obtained by using 
approximations, in contrast to exact FFT-based methods. Approximation parameters are 
chosen to keep an error, calculated according to some estimate, in certain bounds. The more 
exact the error estimate is, the less computations are required; thus, the faster the whole 
algorithm. Therefore, algorithm complexity is directly connected to error analysis [131]. 
Koc and Chew [129] described the application of multilevel FMM to the DDA. They 
used semi-empirical formulae to determine the number of terms in multipole series, and 
obtained O(N) complexity. However, rigorous, close to exact, error analysis is still lacking for 
the FMM applied to the DDA. It will allow obtaining a real algorithm complexity with 
guaranteed accuracy. Such an analysis has been conducted for 2D acoustic scattering [130], 
and for light scattering formulated in terms of surface integrals [131]. In both cases the FMM 
was proven to have an asymptotic complexity O(Nln2(N)). Application of the FMM to 
surface-integrals formulation of light scattering was reviewed by Dembart and Yip [132]. 
Another problem of implementing the FMM is that it is completely dependent upon the 
exact form of the interaction potential ijG . All manuscripts mentioned above deal with 
interaction between point dipoles, i.e. Eq. (14). If a more complex expression for ijG  is used 
(e.g. IT), most of the FMM should be developed anew. This makes integration of the FMM 
and the DDA a formidable problem. 
The FMM is a promising method to calculate light scattering by particles that cannot be 
mapped effectively on a rectangular grid; however, there is still space for improving its theory 
to make it more robust and guarantee certain accuracy. 
The FMM is not the only hierarchical tree method available. For instance, a very 
intuitively simple method was proposed by Barnes and Hut [133,134]. Multipole expansions 
over the center of mass in gravitational computations are used, contrary to geometrical center 
in the FMM. It automatically eliminates the second term in the multipole expansion, and 
allows fast evaluation of monopole terms. Though this method is much simpler and clearer 
than the FMM, it has very little control over the errors that can be studied almost exclusively 
empirically. It can be applied to the DDA without significant increase in the total 
computational errors.7
An alternative approach was proposed by Ding and Tsang [135]. They studied scattering 
from trees and used a sparse matrix iterative approach. The interaction matrix is divided into a 
strong part, which accounts for interaction between nearby dipoles, and a complement weak 
part: . The strong part is sparse and therefore allows quick solution of the linear 
system. The weak part is a small correction that is accounted for iteratively: 
ws AAA +=
yxA =)0(s , . )(w)1(s ll xAyxA −=+ (89)
The authors demonstrate potential of this approach for some test cases. 
4.6 Orientation averaging and repeated calculations 
In many physical applications, one is interested in optical properties of an ensemble of 
randomly oriented particles. When the concentration of particles is small, multiple scattering 
is negligible and the optical properties are obtained by averaging single-particle scattering 
over different particle orientations. More general problems, where particles are not identical 
or multiple scattering is significant, are not considered here. 
Orientation averaging of any scattering property can be described as the integral over 
the Euler’s orientation angles (including a probability distribution function if necessary), 
                                                 
7 Hoekstra AG, unpublished results 
which is brought down to a sum by appropriate quadrature. The problem therefore consists in 
calculation of some scattering property for a set of different orientations of the same particle. 
The easiest way is to calculate it by solving sequentially and independently each problem 
from the set. However, the large size of this set calls for some means of reducing the 
calculations. This is especially relevant when the particle is asymmetric; hence, its optical 
properties are sensitive to particle orientation. Let us further assume for clarity that we are 
interested in the scattering matrix at a certain scattering angle. All the discussion for other 
scattering properties is analogous or even simpler. 
Singham et al. [136] noted that the set of problems described above is physically 
equivalent to a fixed orientation of the particle and different incident and scattering directions. 
The latter are determined by transformation of the laboratory reference frame to the reference 
frame associated with the particle. The amplitude scattering matrix, and hence the Mueller 
matrix, also is transformed along with the reference frame (see e.g. [137] for transformation 
formulae). There are two immediate advantages of using a fixed particle orientation. First, A 
is kept constant (see though discussion below), and therefore the construction of A is done 
only once. Second, the amplitude matrix for any scattering angle is quickly obtained after the 
linear system is solved (for two incident polarizations). Hence, integration over one Euler 
angle is relatively fast. 
The constancy of A can be exploited to further reduce the time of orientation-averaging. 
If  or its LU decomposition is obtained [75,136], a single solution for any right part y can 
be obtained in n
2 operations – the same or less time than required for one iteration using 
general iterative methods (see Subsection 4.1). Moreover, Singham et al. [136] and McClain 
and Ghoul [138] independently proposed an analytical way of averaging the scattering matrix 
at any scattering angle, which requires O(n2) operations once  is known. Khlebtsov [139] 
extended this technique to averaging of extinction and absorption cross sections. 
However, by employing special properties of the matrix A in the DDA allows 
computing matrix-vector products in O(nln(n)) operations (see subsections 4.4, 4.5). Although 
some acceleration of direct methods also can be performed (see Subsection 4.3), they are still 
O(n2) or slower. For large n, iterative methods (assuming that they converge in much less than 
n iterations) are clearly preferable, even if many quadrature points are used. Moreover, large n 
is unattainable by direct methods because of storage requirements. Another improvement 
could be using a heavy preconditioner, which has large initialization cost and greatly 
increases convergence rate. Initialization cost is then justified because it is computed only 
once. Possible candidates are incomplete factorization preconditioners [107]. 
Above it was stated that A is constant for a fixed-orientation particle. However, modern 
DDA formulations (e.g. LDR) take into account the direction of light incidence. Hence A 
depends upon this direction, but only weakly through ( )( )2O kd  corrections. This complicates 
the techniques described above, however probably they still may be used together with some 
special methods to correct for small changes in A on every step. Such methods have not been 
developed as yet. 
Another possibility to perform orientation averaging is to first compute the T-matrix of 
the particle, which then allows analytical averaging [140]. The T-matrix formalism is based 
on the multipole expansion which is truncated at some order N0. Although N0 is hard to 
deduce a priori, usually it is several times x [141,142]. The number of rows in the T-matrix 
equals . The simplest way to evaluate the T-matrix based on the DDA is to solve 
for every incident spherical wave (i.e. for each row of the T-matrix) independently [141]. 
Then the above discussion about optimizing this repeated calculation is relevant. Using 
iterative techniques with N
)2(2 00 +NN
iter number of iterations, computation time is 
( ) ( )[ ]20iter20 O)ln(O NNNNN + , where the first term in the sum is the time for solving the linear 
system, and the second one is the actual computation of the values in the row of the T-matrix. 
A new method to obtain the T-matrix from the DDA interaction matrix was proposed by 
Mackowski [141]. This requires two summations with computational time ( ))ln(O 20 NNN  and 
( )NN 40O . Mackowski showed that for 5=x  his method is an order of magnitude faster than 
the straightforward one. 
Recently Muinonen and Zubko [143] have proposed a way to optimize ensemble 
averaging of DDA results over different sizes and refractive indices. It is based on calculating 
a “good guess” for the initial vector in the iterative solver using results of the calculations 
with similar parameters. Similar ideas can be used to optimize simulation of a set of slightly 
different shapes or orientational averaging. 
Use of repeated calculations to increase the accuracy of DDA simulations was proposed 
recently by Yurkin et al. [63]. Several independent simulations with different discretization 
parameter were performed and results were extrapolated to the infinite discretization giving 
better accuracy than those of a single DDA simulation. 
5 Comparison of the DDA to other methods 
Hovenier et al. [144] compared the DDA, the extended boundary condition method (EBCM), 
and the separation of variables method (SVM) for calculations of scattering by spheroids, 
finite cylinders and bispheres. Parameters of the problems were as follows: , 
equivolume size parameter , 
i01.05.1 +=m
5=x 6.0=y . The angular dependencies of scattering matrix 
elements were calculated. The EBCM and SVM seemed to achieve an exact solution, and the 
DDA showed little errors, except for backscattering angles, where they were up to 10-20%. 
Wriedt and Comberg [145] compared the DDA, EBCM, and finite difference time 
domain (FDTD) method for a cube with 33.1=m , 1.5 and 9.2=x , 4.9, 9.7. For  and 
4.9. The DDA and EBCM achieved good accuracy in calculation of scattering intensity angle 
dependence; the DDA was 2-5 times faster, but consumed 8-16 times more memory (y was in 
the range 0.3-0.5). The FDTD had similar computational requirements as the DDA but was 
less accurate. For  the DDA was the only one to achieve little errors within the given 
computational resources. 
9.2=x
7.9=x
Comberg and Wriedt [146] compared the DDA, GMM (see Subsection 3.3) and the 
generalized multipole technique (GMT) for clusters of a few spheres. A single sphere had x in 
the range 4–20 and 33.1=m , 1.5. All the methods managed to achieve good accuracy, but the 
GMM was one order of magnitude (and for large x even several orders) faster than the other 
two. The DDA and GMT also were used to compute scattering by a cluster of two oblate 
spheroids with  and . The DDA was less accurate and consumed 4 times more 
memory, but was 6 times faster than the GMT. 
5=x 33.1=m
Wriedt et al. [147] compared the DDA, FDTD, GMT, and discrete sources method 
(DSM) for the calculation of light scattering by a red blood cell (RBC) with  and 
. Accuracy of all methods was similar. The DDA and GMT showed similar 
calculation times; they were 7 times faster than the FDTD and 12 times slower than the DSM. 
It should be noted that the latter employed the axisymmetric property of RBC. 
06.1=m
Recently Yurkin et al. [58] systematically compared the DDA and the FDTD for 
spheres with m from 1.02 to 2 and x from 10 to 100, depending on m. It was shown that 
numerical performance of the DDA is much more sensitive to the refractive index than that of 
the FDTD. Therefore, the DDA is preferable for small m, the FDTD for larger m. Cleary, the 
crossover point is not well defined and will depend on the details of the problem at hand as 
well as on the particular implementations of both methods. 
The main advantage of the DDA is that it is one of the most general methods, having a 
very broad range of applicability, limited only by available computational power. The reverse 
of this advantage is that it has almost no means to use the symmetry of the scatterer. Thus the 
DDA is not able to compete with the EBCM for homogenous axisymmetric scatterers. For 
homogenous non-axisymmetric scatterers the DDA is competitive with the EBCM for single-
particle orientation, but the latter allows much faster orientation averaging. The EBCM has 
little applicability to inhomogeneous scatterers, where the DDA can be applied without any 
changes. Comparison between the FDTD and the DDA suggests that the DDA is more 
suitable for small m. It also should be noted that the FDTD is even more general, being easily 
applicable to non-harmonic incident electric fields. Moreover, simulation of one pulse 
incident wave with the FDTD gives the solution for a complete spectrum of incident harmonic 
plane waves, but with a limitation on accuracy. 
6 Concluding remarks 
The DDA has been reviewed using a general framework based on the integral equation for the 
electric field. Although mainstream DDA algorithms as used in several production computer 
programs, has not changed significantly since 1994, many different improvements have been 
proposed since that time. Some of them do improve the accuracy or numerical performance of 
the DDA; however, they still wait for a wide acceptance. It seems that a critical mass of new 
improvements is building up, hopefully resulting in a next breakthrough in the field of the 
DDA. 
In our opinion, future major improvements in the DDA computer implementations will 
be connected with one of the following: 
1) Decreasing shape errors by implementing WD or similar techniques. 
2) Improving polarizability and interaction terms by techniques that are still to be 
developed similar to IT and PEL. 
3) Studying different preconditioners for the DDA interaction matrix, either trying some of 
the known types or developing one considering the special structure of the matrix. 
Item (1) should improve the overall accuracy of the DDA, especially for cases where shape 
errors are dominant, item (2) should expand the DDA applicability region to higher refractive 
indices, and item (3) should boost overall performance, especially for large size parameters 
and/or refractive indices. 
Acknowledgements 
We thank Dan Mackowski for clarifying discussion on the simulations of scattering by 
clusters of spheres and Gorden Videen for critically reading the manuscript and for valuable 
discussions. Our research is supported by Siberian Branch of the Russian Academy of 
Sciences through the grant 2006-03. 
Appendix. Description of used acronyms and symbols 
See Tables A1 and A2. 
Table A1. Acronyms in alphabetical order. 
Acronym Description Sectiona
(L) left Jacobi preconditioner 4.1
I right Jacobi preconditioner 4.1
a1-term (M) dipole term in the Mie theory 3.1
Bi-CGSTAB Bi-CG stabilized 4.1
BT block-Toeplitz 4.3
CEMD coupled electric and magnetic dipole 3.4
CG conjugate gradient 4.1
CGNR CG applied to normalized equation with minimization of residual norm 4.1
CGS CG squared 4.1
CS complex symmetric 4.1
CLDR corrected LDR 3.1
CM Clausius-Mossotti 2
DDA discrete dipole approximation 1
DGF digitized Green’s function 1
DSM discrete sources method 5
EBCM extended boundary condition method 5
FCD filtered coupled dipoles 3.1
FDTD finite difference time domain 5
FFT fast Fourier transform 4.4
FMM fast multipole method 4.5
GMM generalized multiparticle Mie solution 3.3
GMRES generalized minimal residual 4.1
GMT generalized multipole technique 5
GPBi-CG generalized product-type methods based on Bi-CG 4.1
IT integration of Green’s tensor 3.1
LAK Lakhtakia 3.1
LC Lou and Charalampopoulos 3.3
LDR lattice dispersion relation 3.1
MBT multilevel BT 4.3
MoM method of moments 3.1
PEL Peltoniemi 3.1
PP Purcell and Pennypacker 1
QMR quasi-minimal residual 4.1
RBC red blood cell 5
RCB Rahmani, Chaumet, and Bryant 3.1
RDG Rayleigh-Debye-Gans approximation 4.2
RR radiative reaction correction 3.1
SCLDR surface-corrected LDR 3.1
SOF scattering order formulation 4.2
SVM separation of variables method 5
TFQMR transpose free QMR 4.1
WD weighted discretization 3.1
a where it is explained or first appears (if no explanation is given). 
Table A2. Symbols used, Latin and Greek letters in alphabetical order.a
Symbols Description Section
(0) superscript: approximate value (usually under constant field assumption) 2
(n) superscript: after n-th iteration 4.2
>< θcos  asymmetry parameter 2
* superscript: complex conjugate 2
A a matrix 3.1
a kk  2
a radius of (equivalent) sphere 3.1
B correction matrix in SCLDR 3.1, Eq. (64)
b1 – b3 numerical coefficients in polarization prescriptions 3.1
C  tensor of electrostatic solution 3.1, Eq. (60)
Csca, Cabs, Cext scattering, absorption, extinction cross section 2
c speed of light in vacuum 2
d size of a cubical cell 2
E, Einc, Eexc, Eself, Esca (total) electric field, incident, exciting, self-induced, scattered 2
e0 polarization vector of the incident wave, 10 =e  2
e superscript: effective 3.1
F scattering amplitude 2
f a function; 
volume filling factor 
G  free space dyadic Green’s function (tensor) 2
2sG  G  in static limit 
ijG  interaction term 
H superscript: conjugate transpose 3.1
hr impulse response function of a filter 3.1
2, 4.1I , I identity dyadic (tensor), operator (matrix) 
i, j subscript: vector indices 4.3
i imaginary unity 2
i, j subscript: number of the dipole 2
K order of a BT matrix 4.3
k free space wave vector 2
L  self-term dyadic 2
M integral associated with finiteness of V0; 
preconditioner 
M  dyadic associated with M 2
m refractive index (relative) 3.1
N total number of dipoles 2
n rr  2
n size of a matrix 4.1
n′ˆ  external normal to the surface 2
nx, ny, nz sizes of the rectangular lattice 4.3
P polarization 2
p superscript: principal 3.1
q dπ2  3.1
R rr ′−  2
R0 radius of the smallest sphere circumscribing the scatterer 4.2
2r, r′ radius-vectors 
S LDR coefficient dependent on incident polarization 3.1, Eq. (49)
Si amplitude matrix element 3.2
Sij Mueller matrix element 3.2
s superscript: secondary; 
strong; 
subscript: equivalent spherical dipole 
T  boundary condition tensor 3.1, Eq. (66) 
t time 2
V volume of the scatterer 2
V0 exclusion volume 2
w superscript: weak 4.5
x unknown vector 4.1
Table A2 (continued) 
Symbols Description Section
x size parameter of scatterer 3.2
x, y, z Cartesian coordinates 4.3
y a known vector (right side of a linear system) 4.1
y kdm ||  3.2
yRe kdm)Re(  3.2
2α, α  polarizability, tensor 
4.1γ optimal reduction factor 
3.1δ Kronecker symbol 
2ε electric permittivity (relative) 
3.3η correction function 
Λ intermediate tensor in RCB method 3.1, Eq. (62)
4.2Λ linear integral operator, its matrix 
μ, ν, ρ, τ, … sub-, superscript: Cartesian components of vectors (tensors) 2
ξ, ψ Riccati-Bessel functions 3.1
2χ electric susceptibility 
3.2Ψ mean relative error of far-field electric field 
2Ω solid angle 
2ω circular frequency of the harmonic electric field 
a common sub- and superscripts are given on their own. For all vectors – the same symbol but in italic (instead of 
bold) denotes Euclidian norm of the vector (except unitary vectors). 
References 
 [1]  Purcell EM, Pennypacker CR. Scattering and adsorption of light by nonspherical dielectric grains. 
Astrophys J 1973;186:705-714. 
 [2]  Draine BT. The discrete-dipole approximation and its application to interstellar graphite grains. 
Astrophys J 1988;333:848-872. 
 [3]  Draine BT, Goodman JJ. Beyond clausius-mossotti - wave-propagation on a polarizable point lattice and 
the discrete dipole approximation. Astrophys J 1993;405:685-697. 
 [4]  Draine BT, Flatau PJ. Discrete-dipole approximation for scattering calculations. J Opt Soc Am A 
1994;11:1491-1499. 
 [5]  Draine BT. The discrete dipole approximation for light scattering by irregular targets. In: Mishchenko 
MI, Hovenier, JW, Travis, LD, editors. Light Scattering by Nonspherical Particles, Theory, 
Measurements, and Applications. New York: Academic Press, 2000. p. 131-145. 
 [6]  Draine BT, Flatau PJ. User guide for the discrete dipole approximation code DDSCAT 6.1. 
http://xxx.arxiv.org/abs/astro-ph/0409262, 2004. 
 [7]  Goedecke GH, O'Brien SG. Scattering by irregular inhomogeneous particles via the digitized Green's 
function algorithm. Appl Opt 1988;27:2431-2438. 
 [8]  Lakhtakia A. 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 C 1992;3:583-603. 
 [9]  Rahola J. Solution of dense systems of linear equations in the discrete-dipole approximation. SIAM J Sci 
Comp 1996;17:78-89. 
 [10]  Piller NB. Coupled-dipole approximation for high permittivity materials. Opt Comm 1999;160:10-14. 
 [11]  Chaumet PC, Sentenac A, Rahmani A. Coupled dipole method for scatterers with large permittivity. 
Phys Rev E 2004;70:036606. 
 [12]  Singham SB, Bohren CF. Light scattering by an arbitrary particle: a physical reformulation of the 
coupled dipole method. Opt Lett 1987;12:10-12. 
 [13]  Piller NB. Influence of the edge meshes on the accuracy of the coupled-dipole approximation. Opt Lett 
1997;22:1674-1676. 
 [14]  Hage JI, Greenberg JM, Wang RT. Scattering from arbitrarily shaped particles - theory and experiment. 
Appl Opt 1991;30:1141-1152. 
 [15]  Kahnert FM. Numerical methods in electromagnetic scattering theory. J Quant Spectrosc Radiat Transf 
2003;79:775-824. 
 [16]  Peterson AW, Ray SL, Mittra R. Computational Methods of Electromagnetic Scattering. IEEE Press, 
1998. 
 [17]  Kim OS, Meincke P, Breinbjerg O, Jorgensen E. Method of moments solution of volume integral 
equations using higher-order hierarchical Legendre basis functions. Radio Science 2004;39. 
 [18]  Lu CC. A fast algorithm based on volume integral equation for analysis of arbitrarily shaped dielectric 
radomes. IEEE Trans Ant Propag 2003;51:606-612. 
http://xxx.arxiv.org/abs/astro-ph/0409262,
 [19]  Ivakhnenko V, Eremin Y. Light scattering by needle-type and disk-type particles. J Quant Spectrosc 
Radiat Transf 2006;100:165-172. 
 [20]  Wriedt T. A review of elastic light scattering theories. Part Part Sys Charact 1998;15:67-74. 
 [21]  Chiappetta P, Torresani B. Some approximate methods for computing electromagnetic fields scattered by 
complex objects. Meas Sci Technol 1998;9:171-182. 
 [22]  Mishchenko MI, Travis LD, Lacis AA. Scattering, Absorption, and Emission of Light by Small Particles. 
Cambridge: Cambridge University Press, 2002. 
 [23]  Tsang L, Kong JA, Ding KH, Ao CO. Scattering of Electromagnetic Waves: Numerical Simulations. 
New York: Wiley, 2001. 
 [24]  Jones AR. Light scattering for particle characterization. Prog Ener Comb Sci 1999;25:1-53. 
 [25]  Yanghjian AD. Electric dyadic Green's function in the source region. IEEE Proc 1980;68:248-263. 
 [26]  Peltoniemi JI. Variational volume integral equation method for electromagnetic scattering by irregular 
grains. J Quant Spectrosc Radiat Transf 1996;55:637-647. 
 [27]  Bohren CF, Huffman DR. Absorption and scattering of Light by Small Particles. New York: Wiley, 
1983. 
 [28]  Draine BT, Weingartner JC. Radiative torques on interstellar grains .1. Superthermal spin-up. Astrophys 
J 1996;470:551-565. 
 [29]  Hoekstra AG, Frijlink M, Waters LBFM, Sloot PMA. Radiation forces in the discrete-dipole 
approximation. J Opt Soc Am A 2001;18:1944-1953. 
 [30]  Chaumet PC, Rahmani A, Sentenac A, Bryant GW. Efficient computation of optical forces with the 
coupled dipole method. Phys Rev E 2005;72:046708. 
 [31]  Hoekstra AG. Computer simulations of elastic light scattering. PhD thesis. University of Amsterdam, 
Amsterdam, 1994. 
 [32]  Yurkin MA, Maltsev VP, Hoekstra AG. Convergence of the discrete dipole approximation. I. Theoretical 
analysis. J Opt Soc Am A 2006;23:2578-2591. 
 [33]  Jackson JD. Classical Electrodynamics. New York: Wiley, 1975. 
 [34]  Iskander MF, Chen HY, Penner JE. Optical-scattering and absorption by branched chains of aerosols. 
Appl Opt 1989;28:3083-3091. 
 [35]  Hage JI, Greenberg JM. A model for the optical-properties of porous grains. Astrophys J 1990;361:251-
259. 
 [36]  Livesay DE, Chen KM. Electromagnetic fields induced inside arbitrarily shaped biological bodies. IEEE 
Trans Microw Theory Tech 1974;22:1273-1280. 
 [37]  Lakhtakia A, Mulholland GW. On 2 numerical techniques for light-scattering by dielectric agglomerated 
structures. J Res Nat Inst Stand Technol 1993;98:699-716. 
 [38]  Dungey CE, Bohren CF. Light-scattering by nonspherical particles - a refinement to the coupled-dipole 
method. J Opt Soc Am A 1991;8:81-87. 
 [39]  Doyle WT. Optical properties of a suspension of metal spheres. Phys Rev B 1989;39:9852-9858. 
 [40]  Lumme K, Rahola J. Light-scattering by porous dust particles in the discrete-dipole approximation. 
Astrophys J 1994;425:653-667. 
 [41]  van de Hulst HC. Light Scattering by Small Particles. New York: Dover, 1981. 
 [42]  Okamoto H. Light scattering by clusters: the a1-term method. Opt Rev 1995;2:407-412. 
 [43]  Gutkowicz-Krusin D, Draine BT. Propagation of electromagnetic waves on a rectangular lattice of 
polarizable points. http://xxx.arxiv.org/abs/astro-ph/0403082, 2004. 
 [44]  Piller NB, Martin OJF. Increasing the performance of the coupled-dipole approximation: A spectral 
approach. IEEE Trans Ant Propag 1998;46:1126-1137. 
 [45]  Gay-Balmaz P, Martin OJF. A library for computing the filtered and non-filtered 3D Green's tensor 
associated with infinite homogeneous space and surfaces. Comp Phys Comm 2002;144:111-120. 
 [46]  Rahmani A, Chaumet PC, Bryant GW. Coupled dipole method with an exact long-wavelength limit and 
improved accuracy at finite frequencies. Opt Lett 2002;27:2118-2120. 
 [47]  Collinge MJ, Draine BT. Discrete-dipole approximation with polarizabilities that account for both finite 
wavelength and target geometry. J Opt Soc Am A 2004;21:2023-2028. 
 [48]  Rahmani A, Chaumet PC, Bryant GW. On the importance of local-field corrections for polarizable 
particles on a finite lattice: Application to the discrete dipole approximation. Astrophys J 2004;607:873-
878. 
 [49]  Evans KF, Stephens GL. Microwave radiative-transfer through clouds composed of realistically shaped 
ice crystals .1. Single scattering properties. J Atmos Sci 1995;52:2041-2057. 
 [50]  Flatau PJ, Fuller KA, Mackowski DW. Scattering by 2 spheres in contact - comparisons between 
discrete-dipole approximation and modal-analysis. Appl Opt 1993;32:3302-3305. 
 [51]  Xu YL, Gustafson BAS. Comparison between multisphere light-scattering calculations: Rigorous 
solution and discrete-dipole approximation. Astrophys J 1999;513:894-909. 
 [52]  Ku JC. Comparisons of coupled-dipole solutions and dipole refractive-indexes for light-scattering and 
absorption by arbitrarily shaped or agglomerated particles. J Opt Soc Am A 1993;10:336-342. 
http://xxx.arxiv.org/abs/astro-ph/0403082,
 [53]  Andersen AC, Mutschke H, Posch T, Min M, Tamanai A. Infrared extinction by homogeneous particle 
aggregates of SiC, FeO and SiO2: Comparison of different theoretical approaches. J Quant Spectrosc 
Radiat Transf 2006;100:4-15. 
 [54]  Singham SB. Theoretical factors in modeling polarized light scattering by arbitrary particles. Appl Opt 
1989;28:5058-5064. 
 [55]  Hoekstra AG, Sloot PMA. Dipolar unit size in coupled-dipole calculations of the scattering matrix-
elements. Opt Lett 1993;18:1211-1213. 
 [56]  Hoekstra AG, Rahola J, Sloot PMA. Accuracy of internal fields in volume integral equation simulations 
of light scattering. Appl Opt 1998;37:8482-8497. 
 [57]  Druger SD, Bronk BV. Internal and scattered electric fields in the discrete dipole approximation. J Opt 
Soc Am B 1999;16:2239-2246. 
 [58]  Yurkin MA, Brock RS, Lu JQ, Hoekstra AG. Systematic comparison of the discrete dipole 
approximation and the finite difference time domain method. (in preparation) 
 [59]  Okamoto H, Macke A, Quante M, Raschke E. Modeling of backscattering by non-spherical ice particles 
for the interpretation of cloud radar signals at 94 GHz. An error analysis. Contrib Atmos Phys 
1995;68:319-334. 
 [60]  Liu CL, Illingworth AJ. Error analysis of backscatter from discrete dipole approximation for different ice 
particle shapes. Atmos Res 1997;44:231-241. 
 [61]  Lemke H, Okamoto H, Quante M. Comment on error analysis of backscatter from discrete dipole 
approximation for different ice particle shapes [ Liu, C.-L., Illingworth, A.J., 1997, Atmos. Res. 44, 231-
241.]. Atmos Res 1998;49:189-197. 
 [62]  Liu CL, Illingworth AJ. Reply to comment by Lemke, Okamoto and Quante on 'Error analysis of 
backscatter from discrete dipole approximation for different ice particle shapes'. Atmos Res 1999;50:1-2. 
 [63]  Yurkin MA, Maltsev VP, Hoekstra AG. Convergence of the discrete dipole approximation. II. An 
extrapolation technique to increase the accuracy. J Opt Soc Am A 2006;23:2592-2601. 
 [64]  Fuller KA, Mackowski DW. Electromagnetic scattering by compounded spherical particles. In: 
Mishchenko MI, Hovenier, JW, Travis, LD, editors. Light Scattering by Nonspherical Particles, Theory, 
Measurements, and Applications. New York: Academic Press, 2000. p. 223-272. 
 [65]  Xu YL. Scattering Mueller matrix of an ensemble of variously shaped small particles. J Opt Soc Am A 
2003;20:2093-2105. 
 [66]  Mackowski DW. Electrostatics analysis of radiative absorption by sphere clusters in the rayleigh limit - 
application to soot particles. Appl Opt 1995;34:3535-3545. 
 [67]  Mackowski DW. Calculation of total cross-sections of multiple-sphere clusters. J Opt Soc Am A 
1994;11:2851-2861. 
 [68]  Ngo D, Videen G, Dalling R. Chaotic light scattering from a system of osculating, conducting spheres. 
Physics Letters A 1997;227:197-202. 
 [69]  Markel VA, Pustovit VN, Karpov SV, Obuschenko AV, Gerasimov VS, Isaev IL. Electromagnetic 
density of states and absorption of radiation by aggregates of nanospheres with multipole interactions. 
Phys Rev B 2004;70:054202. 
 [70]  Kim HY, Sofo JO, Velegol D, Cole MW, Mukhopadhyay G. Static polarizabilities of dielectric 
nanoclusters. Phys Rev A 2005;72:053201. 
 [71]  Jones AR. Electromagnetic wave scattering by assemblies of particles in the Rayleigh approximation. 
Proc R Soc London A 1979;366:111-127. 
 [72]  Jones AR. Scattering efficiency factors for agglomerates for small spheres. J Phys D 1979;12:1661-1672. 
 [73]  Kozasa T, Blum J, Mukai T. Optical-properties of dust aggregates .1. Wavelength dependence. Astron 
Astrophys 1992;263:423-432. 
 [74]  Kozasa T, Blum J, Okamoto H, Mukai T. Optical-properties of dust aggregates .2. Angular-dependence 
of scattered-light. Astron Astrophys 1993;276:278-288. 
 [75]  Lou WJ, Charalampopoulos TT. On the electromagnetic scattering and absorption of agglomerated small 
spherical-particles. J Phys D 1994;27:2258-2270. 
 [76]  Markel VA, Shalaev VM, Stechel EB, Kim W, Armstrong RL. Small-particle composites .1. Linear 
optical properties. Phys Rev B 1996;53:2425-2436. 
 [77]  Pustovit VN, Sotelo JA, Niklasson GA. Coupled multipolar interactions in small-particle metallic 
clusters. J Opt Soc Am A 2002;19:513-518. 
 [78]  Lumme K, Rahola J, Hovenier JW. Light scattering by dense clusters of spheres. Icarus 1997;126:455-
469. 
 [79]  Kimura H, Mann I. Light scattering by large clusters of dipoles as an analog for cometary dust 
aggregates. J Quant Spectrosc Radiat Transf 2004;89:155-164. 
 [80]  Hull P, Shepherd I, Hunt A. Modeling light scattering from Diesel soot particles. Appl Opt 
2004;43:3433-3441. 
 [81]  Venizelos DT, Lou WJ, Charalampopoulos TT. Development of an algorithm for the calculation of the 
scattering properties of agglomerates. Appl Opt 1996;35:542-548. 
 [82]  Voshchinnikov NV, Il'in VB, Henning T. Modelling the optical properties of composite and porous 
interstellar grains. Astron Astrophys 2005;429:371-381. 
 [83]  Kohler M, Kimura H, Mann I. Applicability of the discrete-dipole approximation to light-scattering 
simulations of large cosmic dust aggregates. Astron Astrophys 2006;448:395-399. 
 [84]  Zubko E, Petrov D, Shkuratov Y, Videen G. Discrete dipole approximation simulations of scattering by 
particles with hierarchical structure. Appl Opt 2005;44:6479-6485. 
 [85]  Bourrely C, Chiappetta P, Lemaire TJ, Torresani B. Multidipole formulation of the coupled dipole 
method for electromagnetic scattering by an arbitrary particle. J Opt Soc Am A 1992;9:1336-1340. 
 [86]  Rouleau F, Martin PG. A new method to calculate the extinction properties of irregularly shaped 
particles. Astrophys J 1993;414:803-814. 
 [87]  Mulholland GW, Bohren CF, Fuller KA. Light-scattering by agglomerates - coupled electric and 
magnetic dipole method. Langmuir 1994;10:2533-2546. 
 [88]  Lemaire TJ. Coupled-multipole formulation for the treatment of electromagnetic scattering by a small 
dielectric particle of arbitrary shape. J Opt Soc Am A 1997;14:470-474. 
 [89]  Lakhtakia A. General-theory of the purcell-pennypacker scattering approach and its extension to 
bianisotropic scatterers. Astrophys J 1992;394:494-499. 
 [90]  Loiko VA, Molochko VI. Polymer dispersed liquid crystal droplets: Methods of calculation of optical 
characteristics. Liq Crys 1998;25:603-612. 
 [91]  Smith DA, Stokes KL. Discrete dipole approximation for magneto-optical scattering calculations. Opt 
Expr 2006;14:5746-5754. 
 [92]  Su CC. Electromagnetic scattering by a dielectric body with arbitrary inhomogeneity and anisotropy. 
IEEE Trans Ant Propag 1989;37:384-389. 
 [93]  Chen RS, Fan ZH, Yung EKN. Analysis of electromagnetic scattering of three-dimensional dielectric 
bodies using Krylov subspace FFT iterative methods. Microwave Opt Tech Lett 2003;39:261-267. 
 [94]  Khlebtsov NG. An approximate method for calculating scattering and absorption of light by fractal 
aggregates. Opt Spec 2000;88:594-601. 
 [95]  Markel VA. Coupled-dipole approach to scattering of light from a one-dimensional periodic dipole 
structure. J Mod Opt 1993;40:2281-2291. 
 [96]  Chaumet PC, Rahmani A, Bryant GW. Generalization of the coupled dipole method to periodic 
structures. Phys Rev B 2003;67:165404. 
 [97]  Chaumet PC, Sentenac A. Numerical simulations of the electromagnetic field scattered by defects in a 
double-periodic structure. Phys Rev B 2005;72:205437. 
 [98]  Martin OJF. Efficient scattering calculations in complex backgrounds. AEU-Int J Electr Comm 
2004;58:93-99. 
 [99]  Yang WH, Schatz GC, Vanduyne RP. Discrete dipole approximation for calculating extinction and 
raman intensities for small particles with arbitrary shapes. J Chem Phys 1995;103:869-875. 
 [100]  Lemaire TJ, Bassrei A. Three-dimensional reconstruction of dielectric objects by the coupled-dipole 
method. Appl Opt 2000;39:1272-1278. 
 [101]  Belkebir K, Chaumet PC, Sentenac A. Superresolution in total internal reflection tomography. J Opt Soc 
Am A 2005;22:1889-1897. 
 [102]  Chaumet PC, Belkebir K, Sentenac A. Three-dimensional subwavelength optical imaging using the 
coupled dipole method. Phys Rev B 2004;69:245405. 
 [103]  Chaumet PC, Belkebir K, Lencrerot R. Three-dimensional optical imaging in layered media. Opt Expr 
2006;14:3415-3426. 
 [104]  Zubko E, Shkuratov Y, Videen G. Discrete-dipole analysis of backscatter features of agglomerated 
debris particles comparable in size with wavelength. J Quant Spectrosc Radiat Transf 2006;100:483-488. 
 [105]  Penttila A, Zubko E, Lumme K, Muinonen K, Yurkin MA, Draine BT, Rahola J, Hoekstra AG, 
Shkuratov Y. Comparison between discrete dipole implementations and exact techniques. J Quant 
Spectrosc Radiat Transf 2007, doi:10.1016/j.jqsrt.2007.01.26. 
 [106]  Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical Recipes in C. The Art of Scientific 
Computing. New York: Cambridge University Press, 1990. 
 [107]  Barrett R, Berry M, Chan TF, Demmel J, Donato J, Dongarra J, Eijkhout V, Pozo R, Romine C, van der 
Vorst HA. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. SIAM, 
1994. 
 [108]  Hoekstra AG, Grimminck MD, Sloot PMA. Large scale simulations of elastic light scattering by a fast 
discrete dipole approximation. Int J Mod Phys C 1998;9:87-102. 
 [109]  Zhang SL. GPBi-CG: Generalized product-type methods based on Bi-CG for solving nonsymmetric 
linear systems. SIAM J Sci Comp 1997;18:537-551. 
 [110]  Freund RW. Conjugate gradient-type methods for linear-systems with complex symmetrical coefficient 
matrices. SIAM J Sci Stat Comp 1992;13:425-448. 
 [111]  Flatau PJ. Improvements in the discrete-dipole approximation method of computing scattering and 
absorption. Opt Lett 1997;22:1205-1207. 
 [112]  Fan ZH, Wang DX, Chen RS, Yung EKN. The application of iterative solvers in discrete dipole 
approximation method for computing electromagnetic scattering. Microwave Opt Tech Lett 
2006;48:1741-1746. 
 [113]  Yurkin MA, Maltsev VP, Hoekstra AG. The discrete dipole approximation for simulation of light 
scattering by particles much larger than the wavelength. J Quant Spectrosc Radiat Transf 2007, 
doi:10.1016/j.jqsrt.2007.01.33. 
 [114]  Rahola J. On the eigenvalues of the volume integral operator of electromagnetic scattering. SIAM J Sci 
Comp 2000;21:1740-1754. 
 [115]  Budko NV, Samokhin AB. Spectrum of the volume integral operator of electromagnetic scattering. 
SIAM J Sci Comp 2006;28:682-700. 
 [116]  Budko NV, Samokhin AB, Samokhin AA. A generalized overrelaxation method for solving singular 
volume integral equations in low-frequency scattering problems. Differ Eq 2005;41:1262-1266. 
 [117]  Hoekstra AG, Sloot PMA. Coupled dipole simulations of elastic light scattering on parallel systems. Int J 
Mod Phys C 1995;6:663-679. 
 [118]  Acquista C. Light scattering by tenuous particles: a generalization of the Rayleigh-Gans-Rocard 
approach. Appl Opt 1976;15:2932-2936. 
 [119]  Chiappetta P. Multiple scattering approach to light scattering by arbitrarily shaped particles. J Phys A 
1980;13:2101-2108. 
 [120]  Singham SB, Bohren CF. Light-scattering by an arbitrary particle - the scattering-order formulation of 
the coupled-dipole method. J Opt Soc Am A 1988;5:1867-1872. 
 [121]  de Hoop AT. Convergence criterion for the time-domain iterative Born approximation to scattering by an 
inhomogeneous, dispersive object. J Opt Soc Am A 1991;8:1256-1260. 
 [122]  Flatau PJ, Stephens GL, Draine BT. Light-scattering by rectangular solids in the discrete-dipole 
approximation - a new algorithm exploiting the block-Toeplitz structure. J Opt Soc Am A 1990;7:593-
600. 
 [123]  Flatau PJ. Fast solvers for one dimensional light scattering in the discrete dipole approximation. Opt 
Expr 2004;12:3149-3155. 
 [124]  Goodman JJ, Draine BT, Flatau PJ. Application of fast-Fourier-transform techniques to the discrete-
dipole approximation. Opt Lett 1991;16:1198-1200. 
 [125]  Barrowes BE, Teixeira FL, Kong JA. Fast algorithm for matrix-vector multiply of asymmetric multilevel 
block-Toeplitz matrices in 3-D scattering. Microwave Opt Tech Lett 2001;31:28-32. 
 [126]  Greengard L, Rokhlin V. A fast algorithm for particle simulations. J Comp Phys 1987;73:325-348. 
 [127]  Rahola J. Diagonal forms of the translation operators in the fast multipole algorithm for scattering 
problems. BIT 1996;36:333-358. 
 [128]  Hoekstra AG, Sloot PMA. New computational techniques to simulate light-scattering from arbitrary 
particles. Part Part Sys Charact 1994;11:189-193. 
 [129]  Koc S, Chew WC. Multilevel fast multipole algorithm for the discrete dipole approximation. J Electrom 
Wav Applic 2001;15:1447-1468. 
 [130]  Amini S, Profit ATJ. Multi-level fast multipole solution of the scattering problem. Engin Anal Bound 
Elem 2003;27:547-564. 
 [131]  Darve E. The fast multipole method I: error analysis and asymptotic complexity. SIAM J Num Anal 
2000;38:98-128. 
 [132]  Dembart B, Yip E. The accuracy of fast multipole methods for Maxwell's equations. IEEE Comp Sci 
Engin 1998;5:48-56. 
 [133]  Barnes JE, Hut P. A hierarchical O(N log N) force-calculation algorithm. Nature 1986;324:446-449. 
 [134]  Barnes JE, Hut P. Error analysis of a tree code. Astrophys J Suppl 1989;70:389-417. 
 [135]  Ding KH, Tsang L. A sparse matrix iterative approach for modeling tree scattering. Microwave Opt Tech 
Lett 2003;38:198-202. 
 [136]  Singham MK, Singham SB, Salzman GC. The scattering matrix for randomly oriented particles. J Chem 
Phys 1986;85:3807-3815. 
 [137]  Mishchenko MI. Calculation of the amplitude matrix for a nonspherical particle in a fixed orientation. 
Appl Opt 2000;39:1026-1031. 
 [138]  McClain WM, Ghoul WA. Elastic light scattering by randomly oriented macromolecules: Computation 
of the complete set of observables. J Chem Phys 1986;84:6609-6622. 
 [139]  Khlebtsov NG. Orientational averaging of integrated cross sections in the discrete dipole method. Opt 
Spec 2001;90:408-415. 
 [140]  Mishchenko MI, Travis LD, Mackowski DW. T-matrix computations of light scattering by nonspherical 
particles: A review. J Quant Spectrosc Radiat Transf 1996;55:535-575. 
 [141]  Mackowski DW. Discrete dipole moment method for calculation of the T matrix for nonspherical 
particles. J Opt Soc Am A 2002;19:881-893. 
 [142]  Mishchenko MI. Light-scattering by size shape distributions of randomly oriented axially-symmetrical 
particles of a size comparable to a wavelength. Appl Opt 1993;32:4652-4666. 
 [143]  Muinonen K, Zubko E. Optimizing the discrete-dipole approximation for sequences of scatterers with 
identical shapes but differing sizes or refractive indices. J Quant Spectrosc Radiat Transf 2006;100:288-
294. 
 [144]  Hovenier JW, Lumme K, Mishchenko MI, Voshchinnikov NV, Mackowski DW, Rahola J. 
Computations of scattering matrices of four types of non-spherical particles using diverse methods. J 
Quant Spectrosc Radiat Transf 1996;55:695-705. 
 [145]  Wriedt T, Comberg U. Comparison of computational scattering methods. J Quant Spectrosc Radiat 
Transf 1998;60:411-423. 
 [146]  Comberg U, Wriedt T. Comparison of scattering calculations for aggregated particles based on different 
models. J Quant Spectrosc Radiat Transf 1999;63:149-162. 
 [147]  Wriedt T, Hellmers J, Eremina E, Schuh R. Light scattering by single erythrocyte: Comparison of 
different methods. J Quant Spectrosc Radiat Transf 2006;100:444-456. 
 [148]  Laczik Z. Discrete-dipole-approximation-based light-scattering calculations for particles with a real 
refractive index smaller than unity. Appl Opt 1996;35:3736-3745. 
	1 Introduction
	2 General framework
	3 Various DDA models
	3.1 Theoretical base of the DDA
	3.2 Accuracy of DDA simulations
	3.3 The DDA for clusters of spheres
	3.4 Modifications and extensions of the DDA
	4 Numerical considerations
	4.1 Direct vs. iterative methods
	4.2 Scattering order formulation
	4.3 Block-Toeplitz
	4.4 FFT
	4.5 Fast multipole method
	4.6 Orientation averaging and repeated calculations
	5 Comparison of the DDA to other methods
	6 Concluding remarks
	Acknowledgements
	Appendix. Description of used acronyms and symbols
	References
  /ASCII85EncodePages false
  /AllowTransparency false
  /AutoPositionEPSFiles true
  /AutoRotatePages /None
  /Binding /Left
  /CalGrayProfile (Dot Gain 20%)
  /CalRGBProfile (sRGB IEC61966-2.1)
  /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2)
  /sRGBProfile (sRGB IEC61966-2.1)
  /CannotEmbedFontPolicy /Error
  /CompatibilityLevel 1.4
  /CompressObjects /Tags
  /CompressPages true
  /ConvertImagesToIndexed true
  /PassThroughJPEGImages true
  /CreateJDFFile false
  /CreateJobTicket false
  /DefaultRenderingIntent /Default
  /DetectBlends true
  /DetectCurves 0.0000
  /ColorConversionStrategy /CMYK
  /DoThumbnails false
  /EmbedAllFonts true
  /EmbedOpenType false
  /ParseICCProfilesInComments true
  /EmbedJobOptions true
  /DSCReportingLevel 0
  /EmitDSCWarnings false
  /EndPage -1
  /ImageMemory 1048576
  /LockDistillerParams false
  /MaxSubsetPct 100
  /Optimize true
  /OPM 1
  /ParseDSCComments true
  /ParseDSCCommentsForDocInfo true
  /PreserveCopyPage true
  /PreserveDICMYKValues true
  /PreserveEPSInfo true
  /PreserveFlatness true
  /PreserveHalftoneInfo false
  /PreserveOPIComments true
  /PreserveOverprintSettings true
  /StartPage 1
  /SubsetFonts true
  /TransferFunctionInfo /Apply
  /UCRandBGInfo /Preserve
  /UsePrologue false
  /ColorSettingsFile ()
  /AlwaysEmbed [ true
  /NeverEmbed [ true
  /AntiAliasColorImages false
  /CropColorImages true
  /ColorImageMinResolution 300
  /ColorImageMinResolutionPolicy /OK
  /DownsampleColorImages true
  /ColorImageDownsampleType /Bicubic
  /ColorImageResolution 300
  /ColorImageDepth -1
  /ColorImageMinDownsampleDepth 1
  /ColorImageDownsampleThreshold 1.50000
  /EncodeColorImages true
  /ColorImageFilter /DCTEncode
  /AutoFilterColorImages true
  /ColorImageAutoFilterStrategy /JPEG
  /ColorACSImageDict <<
    /QFactor 0.15
    /HSamples [1 1 1 1] /VSamples [1 1 1 1]
  /ColorImageDict <<
    /QFactor 0.15
    /HSamples [1 1 1 1] /VSamples [1 1 1 1]
  /JPEG2000ColorACSImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 30
  /JPEG2000ColorImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 30
  /AntiAliasGrayImages false
  /CropGrayImages true
  /GrayImageMinResolution 300
  /GrayImageMinResolutionPolicy /OK
  /DownsampleGrayImages true
  /GrayImageDownsampleType /Bicubic
  /GrayImageResolution 300
  /GrayImageDepth -1
  /GrayImageMinDownsampleDepth 2
  /GrayImageDownsampleThreshold 1.50000
  /EncodeGrayImages true
  /GrayImageFilter /DCTEncode
  /AutoFilterGrayImages true
  /GrayImageAutoFilterStrategy /JPEG
  /GrayACSImageDict <<
    /QFactor 0.15
    /HSamples [1 1 1 1] /VSamples [1 1 1 1]
  /GrayImageDict <<
    /QFactor 0.15
    /HSamples [1 1 1 1] /VSamples [1 1 1 1]
  /JPEG2000GrayACSImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 30
  /JPEG2000GrayImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 30
  /AntiAliasMonoImages false
  /CropMonoImages true
  /MonoImageMinResolution 1200
  /MonoImageMinResolutionPolicy /OK
  /DownsampleMonoImages true
  /MonoImageDownsampleType /Bicubic
  /MonoImageResolution 1200
  /MonoImageDepth -1
  /MonoImageDownsampleThreshold 1.50000
  /EncodeMonoImages true
  /MonoImageFilter /CCITTFaxEncode
  /MonoImageDict <<
    /K -1
  /AllowPSXObjects false
  /CheckCompliance [
    /None
  /PDFX1aCheck false
  /PDFX3Check false
  /PDFXCompliantPDFOnly false
  /PDFXNoTrimBoxError true
  /PDFXTrimBoxToMediaBoxOffset [
    0.00000
    0.00000
    0.00000
    0.00000
  /PDFXSetBleedBoxToMediaBox true
  /PDFXBleedBoxToTrimBoxOffset [
    0.00000
    0.00000
    0.00000
    0.00000
  /PDFXOutputIntentProfile ()
  /PDFXOutputConditionIdentifier ()
  /PDFXOutputCondition ()
  /PDFXRegistryName ()
  /PDFXTrapped /False
  /Description <<
    /CHS <FEFF4f7f75288fd94e9b8bbe5b9a521b5efa7684002000410064006f006200650020005000440046002065876863900275284e8e9ad88d2891cf76845370524d53705237300260a853ef4ee54f7f75280020004100630072006f0062006100740020548c002000410064006f00620065002000520065006100640065007200200035002e003000204ee553ca66f49ad87248672c676562535f00521b5efa768400200050004400460020658768633002>
    /CHT <FEFF4f7f752890194e9b8a2d7f6e5efa7acb7684002000410064006f006200650020005000440046002065874ef69069752865bc9ad854c18cea76845370524d5370523786557406300260a853ef4ee54f7f75280020004100630072006f0062006100740020548c002000410064006f00620065002000520065006100640065007200200035002e003000204ee553ca66f49ad87248672c4f86958b555f5df25efa7acb76840020005000440046002065874ef63002>
    /DAN <FEFF004200720075006700200069006e0064007300740069006c006c0069006e006700650072006e0065002000740069006c0020006100740020006f007000720065007400740065002000410064006f006200650020005000440046002d0064006f006b0075006d0065006e007400650072002c0020006400650072002000620065006400730074002000650067006e006500720020007300690067002000740069006c002000700072006500700072006500730073002d007500640073006b007200690076006e0069006e00670020006100660020006800f8006a0020006b00760061006c0069007400650074002e0020004400650020006f007000720065007400740065006400650020005000440046002d0064006f006b0075006d0065006e0074006500720020006b0061006e002000e50062006e00650073002000690020004100630072006f00620061007400200065006c006c006500720020004100630072006f006200610074002000520065006100640065007200200035002e00300020006f00670020006e0079006500720065002e>
    /DEU <FEFF00560065007200770065006e00640065006e0020005300690065002000640069006500730065002000450069006e007300740065006c006c0075006e00670065006e0020007a0075006d002000450072007300740065006c006c0065006e00200076006f006e002000410064006f006200650020005000440046002d0044006f006b0075006d0065006e00740065006e002c00200076006f006e002000640065006e0065006e002000530069006500200068006f006300680077006500720074006900670065002000500072006500700072006500730073002d0044007200750063006b0065002000650072007a0065007500670065006e0020006d00f60063006800740065006e002e002000450072007300740065006c006c007400650020005000440046002d0044006f006b0075006d0065006e007400650020006b00f6006e006e0065006e0020006d006900740020004100630072006f00620061007400200075006e0064002000410064006f00620065002000520065006100640065007200200035002e00300020006f0064006500720020006800f600680065007200200067006500f600660066006e00650074002000770065007200640065006e002e>
    /ESP <FEFF005500740069006c0069006300650020006500730074006100200063006f006e0066006900670075007200610063006900f3006e0020007000610072006100200063007200650061007200200064006f00630075006d0065006e0074006f00730020005000440046002000640065002000410064006f0062006500200061006400650063007500610064006f00730020007000610072006100200069006d0070007200650073006900f3006e0020007000720065002d0065006400690074006f007200690061006c00200064006500200061006c00740061002000630061006c0069006400610064002e002000530065002000700075006500640065006e00200061006200720069007200200064006f00630075006d0065006e0074006f00730020005000440046002000630072006500610064006f007300200063006f006e0020004100630072006f006200610074002c002000410064006f00620065002000520065006100640065007200200035002e003000200079002000760065007200730069006f006e0065007300200070006f00730074006500720069006f007200650073002e>
    /FRA <FEFF005500740069006c006900730065007a00200063006500730020006f007000740069006f006e00730020006100660069006e00200064006500200063007200e900650072002000640065007300200064006f00630075006d0065006e00740073002000410064006f00620065002000500044004600200070006f0075007200200075006e00650020007100750061006c0069007400e90020006400270069006d007000720065007300730069006f006e00200070007200e9007000720065007300730065002e0020004c0065007300200064006f00630075006d0065006e00740073002000500044004600200063007200e900e90073002000700065007500760065006e0074002000ea0074007200650020006f007500760065007200740073002000640061006e00730020004100630072006f006200610074002c002000610069006e00730069002000710075002700410064006f00620065002000520065006100640065007200200035002e0030002000650074002000760065007200730069006f006e007300200075006c007400e90072006900650075007200650073002e>
    /ITA <FEFF005500740069006c0069007a007a006100720065002000710075006500730074006500200069006d0070006f007300740061007a0069006f006e00690020007000650072002000630072006500610072006500200064006f00630075006d0065006e00740069002000410064006f00620065002000500044004600200070006900f900200061006400610074007400690020006100200075006e00610020007000720065007300740061006d0070006100200064006900200061006c007400610020007100750061006c0069007400e0002e0020004900200064006f00630075006d0065006e007400690020005000440046002000630072006500610074006900200070006f00730073006f006e006f0020006500730073006500720065002000610070006500720074006900200063006f006e0020004100630072006f00620061007400200065002000410064006f00620065002000520065006100640065007200200035002e003000200065002000760065007200730069006f006e006900200073007500630063006500730073006900760065002e>
    /JPN <FEFF9ad854c18cea306a30d730ea30d730ec30b951fa529b7528002000410064006f0062006500200050004400460020658766f8306e4f5c6210306b4f7f75283057307e305930023053306e8a2d5b9a30674f5c62103055308c305f0020005000440046002030d530a130a430eb306f3001004100630072006f0062006100740020304a30883073002000410064006f00620065002000520065006100640065007200200035002e003000204ee5964d3067958b304f30533068304c3067304d307e305930023053306e8a2d5b9a306b306f30d530a930f330c8306e57cb30818fbc307f304c5fc59808306730593002>
    /KOR <FEFFc7740020c124c815c7440020c0acc6a9d558c5ec0020ace0d488c9c80020c2dcd5d80020c778c1c4c5d00020ac00c7a50020c801d569d55c002000410064006f0062006500200050004400460020bb38c11cb97c0020c791c131d569b2c8b2e4002e0020c774b807ac8c0020c791c131b41c00200050004400460020bb38c11cb2940020004100630072006f0062006100740020bc0f002000410064006f00620065002000520065006100640065007200200035002e00300020c774c0c1c5d0c11c0020c5f40020c2180020c788c2b5b2c8b2e4002e>
    /NLD (Gebruik deze instellingen om Adobe PDF-documenten te maken die zijn geoptimaliseerd voor prepress-afdrukken van hoge kwaliteit. De gemaakte PDF-documenten kunnen worden geopend met Acrobat en Adobe Reader 5.0 en hoger.)
    /NOR <FEFF004200720075006b00200064006900730073006500200069006e006e007300740069006c006c0069006e00670065006e0065002000740069006c002000e50020006f0070007000720065007400740065002000410064006f006200650020005000440046002d0064006f006b0075006d0065006e00740065007200200073006f006d00200065007200200062006500730074002000650067006e0065007400200066006f00720020006600f80072007400720079006b006b0073007500740073006b00720069006600740020006100760020006800f800790020006b00760061006c0069007400650074002e0020005000440046002d0064006f006b0075006d0065006e00740065006e00650020006b0061006e002000e50070006e00650073002000690020004100630072006f00620061007400200065006c006c00650072002000410064006f00620065002000520065006100640065007200200035002e003000200065006c006c00650072002000730065006e006500720065002e>
    /PTB <FEFF005500740069006c0069007a006500200065007300730061007300200063006f006e00660069006700750072006100e700f50065007300200064006500200066006f0072006d00610020006100200063007200690061007200200064006f00630075006d0065006e0074006f0073002000410064006f0062006500200050004400460020006d00610069007300200061006400650071007500610064006f00730020007000610072006100200070007200e9002d0069006d0070007200650073007300f50065007300200064006500200061006c007400610020007100750061006c00690064006100640065002e0020004f007300200064006f00630075006d0065006e0074006f00730020005000440046002000630072006900610064006f007300200070006f00640065006d0020007300650072002000610062006500720074006f007300200063006f006d0020006f0020004100630072006f006200610074002000650020006f002000410064006f00620065002000520065006100640065007200200035002e0030002000650020007600650072007300f50065007300200070006f00730074006500720069006f007200650073002e>
    /SUO <FEFF004b00e40079007400e40020006e00e40069007400e4002000610073006500740075006b007300690061002c0020006b0075006e0020006c0075006f00740020006c00e400680069006e006e00e4002000760061006100740069007600610061006e0020007000610069006e006100740075006b00730065006e002000760061006c006d0069007300740065006c00750074007900f6006800f6006e00200073006f00700069007600690061002000410064006f0062006500200050004400460020002d0064006f006b0075006d0065006e007400740065006a0061002e0020004c0075006f0064007500740020005000440046002d0064006f006b0075006d0065006e00740069007400200076006f0069006400610061006e0020006100760061007400610020004100630072006f0062006100740069006c006c00610020006a0061002000410064006f00620065002000520065006100640065007200200035002e0030003a006c006c00610020006a006100200075007500640065006d006d0069006c006c0061002e>
    /SVE <FEFF0041006e007600e4006e00640020006400650020006800e4007200200069006e0073007400e4006c006c006e0069006e006700610072006e00610020006f006d002000640075002000760069006c006c00200073006b006100700061002000410064006f006200650020005000440046002d0064006f006b0075006d0065006e007400200073006f006d002000e400720020006c00e4006d0070006c0069006700610020006600f60072002000700072006500700072006500730073002d007500740073006b00720069006600740020006d006500640020006800f600670020006b00760061006c0069007400650074002e002000200053006b006100700061006400650020005000440046002d0064006f006b0075006d0065006e00740020006b0061006e002000f600700070006e00610073002000690020004100630072006f0062006100740020006f00630068002000410064006f00620065002000520065006100640065007200200035002e00300020006f00630068002000730065006e006100720065002e>
    /ENU (Use these settings to create Adobe PDF documents best suited for high-quality prepress printing.  Created PDF documents can be opened with Acrobat and Adobe Reader 5.0 and later.)
  /Namespace [
    (Adobe)
    (Common)
    (1.0)
  /OtherNamespaces [
    <<
      /AsReaderSpreads false
      /CropImagesToFrames true
      /ErrorControl /WarnAndContinue
      /FlattenerIgnoreSpreadOverrides false
      /IncludeGuidesGrids false
      /IncludeNonPrinting false
      /IncludeSlug false
      /Namespace [
        (Adobe)
        (InDesign)
        (4.0)
      ]
      /OmitPlacedBitmaps false
      /OmitPlacedEPS false
      /OmitPlacedPDF false
      /SimulateOverprint /Legacy
    >>
    <<
      /AddBleedMarks false
      /AddColorBars false
      /AddCropMarks false
      /AddPageInfo false
      /AddRegMarks false
      /ConvertColors /ConvertToCMYK
      /DestinationProfileName ()
      /DestinationProfileSelector /DocumentCMYK
      /Downsample16BitImages true
      /FlattenerPreset <<
        /PresetSelector /MediumResolution
      >>
      /FormElements false
      /GenerateStructure false
      /IncludeBookmarks false
      /IncludeHyperlinks false
      /IncludeInteractive false
      /IncludeLayers false
      /IncludeProfiles false
      /MultimediaHandling /UseObjectSettings
      /Namespace [
        (Adobe)
        (CreativeSuite)
        (2.0)
      ]
      /PDFXOutputIntentProfileSelector /DocumentCMYK
      /PreserveEditing true
      /UntaggedCMYKHandling /LeaveUntagged
      /UntaggedRGBHandling /UseDocumentProfile
      /UseDocumentBleed false
    >>
>> setdistillerparams
  /HWResolution [2400 2400]
  /PageSize [612.000 792.000]
>> setpagedevice
