Convergence of the Discrete Dipole Approximation
Convergence of the discrete dipole approximation. 
I. Theoretical analysis. 
Maxim A. Yurkin 
Faculty of Science, Section Computational Science, of the University of Amsterdam, 
Kruislaan 403, 1098 SJ, Amsterdam, The Netherlands 
Institute of Chemical Kinetics and Combustion, Siberian Branch of the Russian Academy of 
Sciences, Institutskaya 3, Novosibirsk 630090 Russia 
myurkin@science.uva.nl
Valeri P. Maltsev 
Institute of Chemical Kinetics and Combustion, Siberian Branch of the Russian Academy of 
Sciences, Institutskaya 3, Novosibirsk 630090 Russia 
Novosibirsk State University, Pirogova Str. 2, 630090, Novosibirsk, Russia 
Alfons G. Hoekstra 
Faculty of Science, Section Computational Science, of the University of Amsterdam, 
Kruislaan 403, 1098 SJ, Amsterdam, The Netherlands 
alfons@sciene.uva.nl
Abstract 
We performed a rigorous theoretical convergence analysis of the discrete dipole 
approximation (DDA). We prove that errors in any measured quantity are bounded by a sum 
of a linear and quadratic term in the size of a dipole d, when the latter is in the range of DDA 
applicability. Moreover, the linear term is significantly smaller for cubically than for non-
cubically shaped scatterers. Therefore, for small d errors for cubically shaped particles are 
much smaller than for non-cubically shaped. The relative importance of the linear term 
decreases with increasing size, hence convergence of DDA for large enough scatterers is 
quadratic in the common range of d. Extensive numerical simulations were carried out for a 
wide range of d. Finally we discuss a number of new developments in DDA and their 
consequences for convergence. 
Keywords: discrete dipole approximation, non-spherical particle light scattering, convergence 
analysis, accuracy 
OCIS code: 290.5850, 260.2110, 000.4430 
mailto:myurkin@science.uva.nl
mailto:alfons@sciene.uva.nl
1. Introduction 
The discrete dipole approximation (DDA) is a well-known method to solve the light 
scattering problem for arbitrary shaped particles. Since its introduction by Purcell and 
Pennypacker1 it has been improved constantly. The formulation of DDA summarized by 
Draine and Flatau2 more than 10 years ago is still most widely used for many applications,3 
partly due to the publicly available high-quality and user-friendly code DDSCAT.4 Although 
modern improvements of DDA (as discussed in detail in Section 2.F) exist, they are still in the 
research stage because they are not widely used in real applications. 
DDA directly discretizes the volume of the scatterer and hence is applicable to arbitrary 
shaped particles. However, the drawback of this discretization is the extreme computational 
complexity of DDA, although it is significantly decreased by advanced numerical 
techniques.2,5 That is why the usual application strategy for DDA is “single computation”, 
where a discretization is chosen based on available computational resources and some 
empirical estimates of the expected errors.3,4 These error estimates are based on a limited 
number of benchmark calculations3 and hence are external to the light scattering problem 
under investigation. Such error estimates have evident drawbacks, however no better 
alternative is available. Some results of analytical analysis of errors in computational 
electromagnetics are known, e.g. 6,7, however they typically consider the surface integral 
equations. To the best of our knowledge, such analysis has not been done for volume integral 
equations (such as DDA). 
Usually errors in DDA are studied as a function of the size parameter of the scatterer x 
(at a constant or few different total numbers of dipoles N), e.g. 2,8. Only a small number of 
papers directly present errors versus discretization parameter (e.g. d – the size of a single 
dipole).9-17 The range of d typically studied in those papers is limited to a 5 times difference 
between minimum and maximum values, with the exception of two papers11,12 where it is 15 
times. Those plots of errors versus discretization parameter are always used to illustrate the 
performance of a new DDA formulation and compare it with others. No conclusions about the 
convergence properties of DDA, as a function of d, have been made based on these plots. To 
our knowledge, no theoretical analysis of DDA convergence has been performed, but only a 
few limited empirical studies have appeared in the literature. 
In this paper we perform a theoretical analysis of DDA convergence when refining the 
discretization (Section 2). We derive rigorous theoretical bounds on the error in any measured 
quantity for any scatterer. In Section 3 we present extensive numerical results of DDA 
computations for 5 different scatterers using many different discretizations. These results are 
discussed in Section 4 to support conclusions of the theoretical analysis. We formulate the 
conclusions of the paper in Section 5. In a follow-up paper18 (which from now on we refer to 
as Paper 2) the theoretical convergence results are used for an extrapolation technique to 
increase the accuracy of DDA computations.  
2. Theoretical analysis 
In this section we analyze theoretically the errors of DDA computations. We formulate the 
volume integral equation for the internal electric field and its operator counterpart in 
Section 2.A and its discretization in Section 2.B. Section 2.C contains integral and discretized 
formulae for measured quantities that are the final goal of any light scattering simulation. We 
derive the main results in Section 2.D, where we consider errors of the traditional DDA 
formulation2 without shape errors, which are considered separately in Section 2.E. Finally in 
Section 2.F we discuss some recent DDA improvements from the viewpoint of our 
convergence theory. 
A.Integral Equation 
Throughout this paper we assume the )iexp( tω−  time dependence of all fields. The scatterer 
is assumed dielectric but not magnetic (magnetic permittivity 1=μ ), and the electric 
permittivity is assumed isotropic (non-isotropic permittivity will significantly complicate the 
derivations but will not principally change the main conclusion of Section 2 – Eqs. (70) and 
(87)). 
The general form of the integral equation governing the electric field inside the 
dielectric scatterer is the following:19,20 
)()(),(),()()(),(d)()( 00
rErrLrMrErrrGrErE χχ VVr
∂−+′′′′+= ∫ , (1)
where Einc(r), 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 (more general – the volume, which 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 
−=∇∇+=′
)()(ˆˆ),(
kRgRgk IIIrrG , (2)
where I  is the identity dyadic, ck ω=  – free space wave vector, rrR ′−= , R=R , and 
 is a dyadic defined as  (μ, ν are Cartesian components of the vector or 
tensor), and g(R) is the scalar Green’s function 
RR ˆˆ νμμν RRRR =ˆˆ
iexp(
)( = . (3)
M is the following integral associated with the finiteness of the exclusion volume V0
( )∫ ′−′′′′=
)()(),()()(),(d),( s30
rV rErrrGrErrrGrM χχ , (4)
where ),(s rrG ′  is the static limit ( ) of 0→k ),( rrG ′ : 
−−=∇∇=′ 23
11ˆˆ),(
IrrG . (5)
L  is the so-called self-term dyadic: 
rV rL , (6)
where  is an external (as viewed from r) normal to the surface ∂Vn′ˆ 0 at point r'. 
Eq. (1) can be rewritten in operator form as follows 
inc~~~ EEA =⋅ , (7)
where ( 311 )~ CE →=∈ VLH  – functions from V to C3 that have finite L1-norm, 2inc~ H∈E  – 
subspace of H1 containing all functions that satisfy Maxwell equations in free space. A
 is a 
linear operator . Although the Sobolev norm is physically more sound (based on 
the finiteness of energy of the electric field),
21: HH →
6,21 we use the L1-norm. A detailed discussion of 
all assumptions made for the electric field is performed in Section 2.D. 
B.Discretization 
To solve Eq. (1) numerically a discretization is done in the following way.20 Let , 
 for . N denotes the number of subvolumes (dipoles). Assuming  and 
choosing , Eq. 
/0=ji VV I ji ≠ iV∈r
iVV =0 (1) can be rewritten as 
)()(),(),()()(),(d)()( 3inc rErrLrMrErrrGrErE χχ ii
∂−+′′′′+= ∑ ∫
. (8)
The set of Eq. (8) (for all i) is exact. Further one fixed point ri inside each Vi (its center) is 
chosen and  is set. irr =
The usual approximation20 is considering E and χ constant inside each subvolume: 
iiiii V∈==== rrrErErE for)()(,)()( χχχ . (9)
Eq. (8) can then be rewritten as 
( ) iiii
jjjijii V ELMEGEE χχ −++= ∑
inc , (10)
where , )(incinc ii rEE = ),( iii V rLL ∂= , 
( )∫ ′−′′=
iii r ),(),(d
s3 rrGrrGM , (11)
∫ ′′=
ij rV
1 3 rrGG . (12)
A further approximation, which is used in almost all formulations of DDA, is 
),()0( jiij rrGG = . (13)
This assumption is made implicitly by all formulations that start by replacing the scatterer 
with a set of point dipoles, as was done originally by Purcell and Pennypacker.1 For a cubical 
(as well as spherical) cell Vi with ri located at the center of the cell, iL  can be calculated 
analytically yielding22 
=i . (14)
Eq. (10) together with Eqs. (13) and (14) and completely neglecting iM  is equivalent to 
the original DDA by Purcell and Pennypacker (PP).1 The diagonal terms in Eq. (10) are then 
equivalent to the well-known Clausius-Mossotti (CM) polarizability for point dipoles. 
Modifications introduced by other DDA prescriptions are discussed in Section 2.F. 
In matrix notation Eq. (10) reads 
ddd ,incEEA = , (15)
where Ed, Einc,d are elements of  (vectors of size N where each element is a 3D complex 
vector) and 
( )N3C
dA  is a N×N matrix where each element is a 3×3 tensor. d is the size of one 
dipole. In operator notation Eq. (8) (for irr = ) is as follows 
( ) diii ,incinc )(~)(~~ ErErEA == , (16)
We define the discretization error function as 
( ) ( )iddidi ,0)(~~ EArEAh −= , (17)
where E0,d is the exact field at the centers of the dipoles – )(
i rEE = , in contrast to E
d that 
is only an approximation obtained from solution of Eq. (15) (here we neglect the numerical 
error that appears from the solution of Eq. (15) itself, which is acceptable if this error is 
controlled to be much less than other errors). Comparing Eqs. (15) and (17) one can 
immediately obtain the error in internal fields due to discretization δEd: 
( ) ddddd hAEEE 1,0δ −−=−= . (18)
C.Measured quantities 
After having determined the internal electric fields, scattered fields and cross sections can be 
calculated. Scattered fields are obtained by taking the limit ∞→r  of the integral in Eq. (1) 
(see e.g. 23) 
)iexp(
)(sca nFrE
= , (19)
where rrn =  is the unit vector in the scattering direction, and F is the scattering amplitude: 
∑∫ ′′⋅′−′−−=
krnnk )()()iexp(d)ˆˆ(i)( 33 rErnrInF χ . (20)
All other differential scattering properties, such as the amplitude and Mueller scattering 
matrices, and asymmetry parameter >< θcos  can be easily derived from F(n), calculated for 
two incident polarizations.24 We consider an incident polarized plane wave: 
)iexp()( 0inc rkerE ⋅= , (21)
where , a is direction of incidence, and ak k= 10 =e  is assumed. The scattering and 
extinction cross sections (Csca, Cext) are derived from the scattering amplitude:23 
∫ Ω= nFkC , (22)
( )∗⋅= 02ext )(Re
, (23)
where * denote complex conjugation. The expression for absorption cross section (Cabs) 
directly uses the internal fields:23 
( )∑ ∫ ′′′=
abs )()(Imd4 rErχπ , (24)
Since only values of the internal field in the centers of dipoles are known, Eqs. (20) and 
(24) are approximated by (PP) 
∑ ⋅−χ−−=
iii kVnnk )iexp()ˆˆ(i)(
3 nrEInF , (25)
iiiVkC
abs )Im(4 Eχπ . (26)
Corrections to Eq. (26) are discussed in Section 2.F. 
Both Eqs. (20) (for each component) and (24) can be generalized as ( )E~~φ  (a functional 
that is not necessarily linear), which is approximated as: 
( ) ( ) ddd φφφ δ~~ += EE , (27)
where ( )dd Eφ  corresponds to Eqs. (25) or (26) respectively, and the error δφ d consists of  two 
parts: 
( ) ( )[ ] ( ) ( )[ ]ddddddd EEEE φφφφφ −+−= ,0,0~~δ . (28)
The first one comes from discretization (similar to Eq. (17)), and the second from errors in the 
internal fields. 
D.Error analysis 
In this section we perform error analysis for the PP formulation of DDA. Improvements of 
DDA are further discussed in Section 2.F. 
We assume cubical subvolumes with size d. We also assume that the shape of the 
particle is exactly described by these cubical subvolumes (we call this cubically shaped 
scatterer). Moreover, χ is a smooth function inside V (exact assumptions on χ are formulated 
below). An extension of the theory to shapes that do not satisfy these conditions is presented 
in Section 2.E. If there are several regions with different values of χ (smooth inside each 
region), the analysis is still valid but interfaces inside V should be considered the same way as 
the outer boundary of V. We further fix the geometry of the scattering problem and incident 
field. Therefore we will be interested only in variation of discretization (which is 
characterized by the single parameter – d); for reasons that will become clear in the sequel, we 
assume that  (this bound is not limiting since otherwise DDA is generally 
inapplicable
We switch to dimensionless parameters by assuming 1=k , which is equivalent to 
measuring all the distances in units of k1 . The unit of electric field can be chosen arbitrary 
but constant. In all further derivations we will use two sets of constants: γi and ci. γ1-γ13 are 
basic constants that do not depend on the discretization d, but do depend directly upon all 
other problem parameters – size parameter eqkRx =  (Req – volume-equivalent radius), m, 
shape, and incident field – or some of them. On the contrary, c1-c94 are auxiliary values that 
either are numerical constants or can be derived in terms of constants γi. Although the 
dependencies of ci on γi are not explicitly derived in this paper, one can easily obtain them 
following the derivations of this section. That is the main motivation for using such vast 
amount of constants instead of an “order of magnitude” formalism. However, such explicit 
derivation has limited application because, as we will see further, constants in the final result 
depend upon almost all basic constants. Qualitative analysis of these dependencies will be 
performed in the end of this section. It should be noted that the main theoretical results 
concerning DDA convergence (boundedness of errors by a quadratic function, cf. Eq (70)) 
can be formulated and applied without consideration of any constants (which is simpler). 
However our full derivation enables us to make additional conclusions related to the behavior 
of specific error terms. 
The total number of dipoles used to discretize the scatterer is 
−= dN γ . (29)
We assume that the internal field E
 is at least four times differentiable and all these 
derivatives are bounded inside V 
65432 )(,)(,)(,)(,)( γγγγγ τρνμρνμνμμ ≤∂∂∂∂≤∂∂∂≤∂∂≤∂≤ rErErErErE  
for V∈r  and τρνμ ,,,∀ . 
This assumption is acceptable since there are no interfaces inside V, therefore E
 should be a 
smooth function. .  denotes the Euclidian (L2) norm, which is used for all 3D objects: vectors 
and tensors. We use the L1-norm, 
. , for N-dimensional vectors and matrices as well as for 
functions and operators. Eq. (30) immediately implies that . We require that χ 
satisfies Eq. 
~ 1 VL∈E
(30) with constants γ7-γ11. Further we will state an estimate for the norm of 
)(RG  and its derivatives. One can easily obtain from Eq. (2) that for 1>R  )(RG  satisfies 
Eq. (30) (with constants c1-c5), while for 2≤R  
6 )(,)(,)(,)(
−−−− ≤∂∂∂≤∂∂≤∂≤ RcRcRcRc RGRGRGRG ρνμνμμ , 
−≤∂∂∂∂ RcRGτρνμ  for τρνμ ,,,∀ . 
Next we state two auxiliary facts that will be used later. Let Vc be a cube with size d and 
with its center at the origin and f(r) a four times differentiable function inside Vc. Then 
)(max)()(d
fdcffr
d cVV
∈,∫ , (32)
( ) )(max)(
)()(d
d cVV
∂∂∂∂+∇≤−
=∫ . (33)
Eqs. (32) and (33) are the corollary of expanding f into Taylor series. Odd orders of the Taylor 
expansion vanish because of cubical symmetry. 
Our first goal is to estimate 
dh . Starting from Eq. (17) we write  as dih
),()(),(d )0(33 ii
i Vdr
rMPGrPrrGh +
−′′′= ∑ ∫
, (34)
where we have introduced the polarization vector for conciseness 
)()()( rErrP χ= , )( ii rPP = . (35)
It is evident that  also satisfies Eq. )(rP (30) (with constants c13-c17). We start by estimating 
),( iiV rM . Substituting a Taylor expansion of  )(rP
( ) ( )( )∑∑ ∂∂+∂+=
ρρ τρ ),,(
)()()( RrP0P0PRP RRR , (36)
where , into Eq μμ Rr ≤≤
~0 (4) gives 
( ) ( )(∫ ∑∫ ∂∂+−=
ii RRRRV
τρτρ τρ ),,(
)()(d),( 3s3 RrPRGPRGRGrM ) . (37)
The norms of these two terms can be estimated as 
( ) 2183s3 )(d3
)()(d dcRRgR
i ≤=− ∫∫ PIPRGRG , (38)
( ) 21923153 )(d3)),,(~()(d dcRRcRRR
≤≤∂∂ ∫∫ ∑ RGRrPRG
τρτρ τρ . (39)
Eq. (38) follows directly from the definitions in Eqs. (2), (5). To derive Eq. (39) we used 
Eq. (31) and the fact that 23RRR ≤∑
τρ . Finally, Eqs. (37)-(39) lead to 
20),( dcV ii ≤rM . (40)
To estimate the sum in Eq. (34) we consider separately three cases: 1) dipole j lies in a 
complete shell of dipole i (we define the shell below); 2) j lies in a distant shell of dipole i – 
1>−= ijijR rr ; 3) all j that fall between the first two cases (see Fig. 1). We define the first 
shell (S1(i)) of a cubical dipole as a set of dipoles that touch it (including touching in one point 
only). The second shell (S2(i)) is a set of dipoles that touch the outer surface of the first shell, 
and so on. The l-th shell (Sl(i)) is then a set of all dipoles that lie on the boundary of the cube 
with size  and center coinciding with the center of the original dipole. We call a shell 
complete if all its elements lie inside the volume of the scatterer V. A shell is called a distant 
dl )12( +
Kmax K(i) 
(2) (3) 
vacuum scatterer 
Fig. 1. Partition of the scatterer’s volume into three regions relative to dipole i. 
shell if all its elements satisfy , i.e. if its order 1>ijR [ ]dKl 1max => . Let K(i) be the order of 
the first incomplete shell, which is an indicator of how close dipole i is to the surface. We 
demand  to separate cases (1) and (2) described above. All j that fall in the third 
case satisfy  (the exact value of this constant – slightly larger than 
max)( KiK ≤
2<ijR 3  – depends on 
d). The number of dipoles in a shell Sl (which can be incomplete) – ns(l) – can be estimated as 
33 )12()12()( lclllns ≤−−+≤ . (41)
The sum of the error over all dipoles that lie in complete shells is then 
∑ ∑ ∫
)0(33 )(),(d
l iSj
dr PGrPrrG , (42)
Since each shell in Eq. (42) is complete it can be divided into pairs of dipoles that are 
symmetric over the center of the shell (j and –j). For convenience we set . The inner 
sum in Eq. 
0r =i
(42) can then be rewritten as 
( ) ( )∑ ∫
+−′−+′′′
)0(33 )()()(d
dr PPGrPrPrG , (43)
Further we introduce the auxiliary function 
( ) )()()(
)( 0PrPrPru −′−+′=′ , (44)
which satisfies the following inequalities (follows from Eq. (30) for P(r) and Taylor series) 
22 )(,)(,)( crcrc ≤∂∂≤∂≤ rururu νμμ  for νμ,∀ . (45)
Then Eq. (43) is equivalent to 
V l jl j
drdr PGrGuGrurG ∑ ∫∑ ∫
)0(33
)0(33 )(d)()(d , (46)
where . To estimate the first term we apply Eq. )( jj ruu = (32) to the whole function under 
the integral. Using Eqs. (31) and (45) one may obtain 
( ) 325)()(max −
≤′′∂∂ ij
, (47)
and hence 
)0(33 )()(d
V ll j
ldcRdcdr uGrurG , (48)
where we have used Eq. (41) and  for ldRij ≥ )(iSj l∈ . 
It is straightforward to show that 
∑ ∫∑ ∫
′′=′′
3 )(d
iSj ViSj V l jl j
rgrr IrG , (49)
)0( )(
RgIG . (50)
The derivation is based upon Eq. (2) and the equivalence I
 in all sums and integrals 
that satisfy cubical symmetry. Then second part of Eq. (46) is transformed to 
)0(33 )()(d)(d −
+≤−′′≤
−′′ ∑ ∫∑ ∫ ldcldcRgdrgrcdr
V l jl j
PGrG , (51)
where we apply Eq. (33) to derive the second inequality and use the identity  
and the following inequalities 
)()(2 rgrg −=∇
( ) 532131 )(, −− ≤∂∂∂∂≤ RcRgRcRg τρνμ  for τρνμ ,,,∀ . (52)
Substituting Eqs. (48) and (51) into Eq. (42) one can obtain 
( ) 23433
)0(33 )(ln)(),(d diKccdr
l iSj
−′′′∑ ∑ ∫
PGrPrrG , (53)
using the fact that . 1)( ≤diK
We now consider the second part of the sum in Eq. (34) (where ). We first apply 
1>ijR
(32), then use Eq. (30) for P(r) and )(rG , and finally invoke Eq. (29): 
)0(33 )(),(d dcdNcdcdr
ijij j
i ≤≤≤⎟
−′′′ ∑∑ ∫
PGrPrrG . (54)
To analyze the third part of the sum in Eq. (34) we again sum over shells, however since 
they are incomplete we cannot use symmetry considerations. We apply Eq. (33) to the whole 
function under the integral and proceed analogous to the derivation of Eq. (51). Using the 
identity 
)()(2 rGrG −=∇ , (55)
(since we have assumed ) we obtain 1=k
( ) 4372 )()( −= ≤∇ ijRcijRrrPrG , (56)
( ) 738)()(max −
≤′′∂∂∂∂ ij
, (57)
which leads to 
)0(33 )(),(d −−
−′′′∑ ∫ lcdlcdr
PGrPrrG , (58)
and then analogous to Eq. (53): 
)()()(),(d 442
)( )(
)0(33
iKcidKcdr
iKl iSj
−′′′∑ ∑ ∫ PGrPrrG . (59)
Collecting Eqs. (40), (53), (54), (59) we finally obtain 
( ) 24443442141 )(ln)()( diKcciKcidKcdi +++≤ −−h . (60)
Then 
( ) ( )442141
max4443
)(ln −−
+++≤= ∑∑ KcdKcKnNdKcc
d hh , (61)
where n(K) is the number of dipoles whose order of the first incomplete shell is equal to K. It 
is clear that 
NdnKn 12)1()( γ≤≤ , (62)
where γ12 is surface to volume ratio of the scatterer. Finally we obtain 
( )[ ]dcddccNd 46245431 ln +−≤h . (63)
The last term in Eq. (63) is mostly determined by dipoles that lie on the surface (or few 
dipoles deep) because it comes from the K-4 term in Eq. (61) (which rapidly decreases when 
moving from surface). We define surface errors as those associated with the linear term in 
Eq. (63). Our numerical simulation (see Section 0) show that this term is small compared to 
other terms for “typical” values of d, however it is always significant for small enough values 
of d. 
From Eq. (18) we directly obtain 
δ ddd hAE
≤ . (64)
We assume that a bounded solution of Eq. (7) uniquely exists for any , moreover we 
assume that if 
inc~ H∈E
inc =E  then 131
γ≤E . These assumptions are equivalent to the fact that 
1~ −A  exists and is finite (the operator 1
~ −A  is bounded). Because dA  is a discretization of A
one would expect that 
( ) 131
lim γ== −
. (65)
Although Eq. (65) seems intuitively correct, its rigorous prove, even if feasible, lies outside 
the scope of this paper. For an intuitive understanding one may consult the paper by Rahola,25 
where he studied the spectrum of the discretized operator (for scattering by a sphere) and 
showed that it does converge to the spectrum of the integral operator with decreasing d. It 
should however be noted, that convergence of the spectrum only implies the convergence of 
the spectral (L2) norm of the operator and not necessarily the convergence of the L1-norm. 
Therefore Eq. (65) should be considered as an assumption. It implies that there exists a d0 
such that 
for 0dd <  ( ) 47
A , (66)
where c47 is an arbitrary constant larger then γ13 (although d0 depends on its choice). For 
example 1347 2γ=c  should lead to a rather large d0 (a rigorous estimate of d0 does not seem 
feasible). Therefore 
δ dE  satisfies the same constrain as 
dh  (Eq. (63)) but with constants 
c48-c50. 
Next we estimate the errors in the measured quantities and start with the discretization 
error (first part in Eq. (28)). Examining Eqs. (20) and (24) one can see that Eq. (32) may be 
directly applied leading to 
( ) ( ) 252551,0~~ dcdc
dd ≤≤− ∑EE φφ . (67)
The second part in Eq. (28) is estimated as 
( ) ( ) ( ) dcddccdcdc d
55541
,0 lnδδ +−≤≤≤− ∑ EEEE φφ , (68)
where we used Eq. (29). The estimation of the error for Cabs additionally uses the fact 
i c EE δδ 57
c Eδmax
57 . 
By combining Eqs. (67) and (68) we obtain the final result of this section: 
( ) dcddccd 5625558 lnδ +−≤φ . (69)
It is important to remember that the derivation was performed for constant x, m, shape, and 
incident field. There are 13 basic constants (γ1-γ13). γ1 (Eq. (29)) characterizes the total 
volume of the scatterer, hence it depends only on x. γ7-γ11 (Eq. (30) for χ(r)) can be easily 
obtained given the function χ(r), moreover it is completely trivial in the common case of 
homogenous scatterers. γ12 (surface to volume ratio, Eq. (62)) depends on the shape of the 
scatterer and is inversely proportional to x. It is not feasible (except for certain simple shapes) 
to obtain the values of constants γ2-γ6 (Eq. (30)), since it requires an exact solution for the 
internal fields. These constants definitely depend upon all the parameters of the scattering 
problem. Moreover, these dependencies can be rapidly varying, especially near the resonance 
regions. The same is true for γ13 (L1-norm of the inverse of the integral operator, Eq. (65)). 
Finally, there is the important constant d0 that also depends on all the parameters, however 
one may expect it to be large enough (e.g. ) for most of the problems – then its 
variation can be neglected. 
20 ≥d
Before proceeding we introduce the discretization parameter kdmy = . We employ the 
commonly used formula as proposed by Draine,8 however the exact dependence on m is not 
important because all the conclusions are still valid for constant m. Replacing d by y does not 
significantly change the dependence of the constants in Eq. (69) since they all already depend 
on m through the basic constants γ2-γ11, γ13. This leads to 
( ) ycyyccy 6126059 lnδ +−≤φ . (70)
It is not feasible to make any rigorous conclusions about the variation of the constants in 
Eq. (70) with varying parameters because all these constants depend on γ2-γ6, γ13 that in turn 
depend in a complex way upon the parameters of the scattering problem. However we can 
make one conclusion about the general trend of this dependency. 
Following the derivation of the Eq. (70) one can observe that c61 is proportional to γ12, 
while c59 and c60 do not directly depend on it (at least part of the contributions to them are 
independent of γ12). Therefore the general trend will be a decrease of the ratio 5961 cc  with 
increasing x (when all other parameters are fixed). This is a mathematical justification of the 
intuitively evident fact that surface errors are less significant for larger particles. 
In the analysis of the results of the numerical simulations (Section 0) we will neglect the 
variation of the logarithm. Eq. (70) then states that error is bounded by a quadratic function of 
y (for ). However, keep in mind that our derivation does not lead to an optimal error 
estimation, i.e. it overestimates the error and can be improved. For example, the constants γ
0dd ≤
γ6 are usually largest inside a small volume fraction of the scatterer (near the surface or some 
internal resonance regions), while in the rest of the scatterer the internal electric field and its 
derivatives are bounded by significantly smaller constants. However the order of the error is 
estimated correctly, as we will see in the numerical simulations. 
It is important to note that Eq. (70) does not imply that δφ y (which is a signed value) 
actually depends on y as a quadratic function, but we will see later that it is the case for small 
enough y (Section 0, see detailed discussion in Paper 2). Moreover, the coefficients of linear 
and quadratic terms for δφ y may have different signs, which may lead to zero error for non-
zero y (however, this y, if it exists, is unfortunately different for each measured quantity). 
E.Shape errors 
In this section we extend the error analysis as presented in Section 2.D to shapes that cannot 
be described exactly by a set of cubical subvolumes. We perform the discretization the same 
way as in Section 2.B but some of the Vi are not cubical (for Vi ∂∈ , which denotes that 
dipole i lies on the boundary of the volume V). We set ri to be still in the center of the cube 
(circumscribing Vi) not to break the regularity of the lattice. The standard PP prescription uses 
equal volumes ( ) in Eqs. 3dVi = (10), (14), (25), and (26), i.e. the discretization changes the 
shape of the particle a little bit. We will estimate the errors introduced by these boundary 
dipoles. These errors should then be added to those obtained in Section 2.D. We start by 
estimating 
dh . First we consider  for dih Vi ∂∉  
−′′′=
PGrPrrGh )0(33 )(),(d , (71)
which is just a reduction of Eq. (34). For Vi ∂∈   is the same plus the error in the self-term dih
iiiiii
i VVdr
EΙrLrMPGrPrrGh χ
−′′′= ∑ ∫
),(),()(),(d )0(33 . (72)
Let us define 
iij dr
PGrPrrGh )0(33sh )(),(d −′′′= ∫ , (73)
iiiiiiii VV EΙrLrMh χ
⎛ −∂−=
),(),(sh . (74)
We estimate each of the terms in Eq. (73) separately (since there is actually no significant 
cancellation, and the error is of the same order of magnitude as the values themselves) using 
Eq. (30) for P(r) and )(rG  and Eq. (31). This leads to 
h  (75)
To estimate  we assume that the surface of the scatterer is a plane on the scale of the size 
of the dipole. A finite radius of curvature only changes the constants in the following 
expressions. We will prove that 
sh cii ≤h , (76)
therefore we do not need to consider the third term in Eq. (74) (coming from the unity tensor) 
at all, since it is bounded by a constant. 
( ) ( )∫∫ −′′′+′′−′′=
iiii rrV )()(),(d)(),(),(d),(
s3s3 rPrPrrGrPrrGrrGrM . (77)
−′ ic rrThe function in the first integral is always bounded by 
. If  the same is true 
for the second integral and hence 
ii V∈r
dcV ii 66),( ≤rM . (78)
If ii V∉r  we introduce an auxilia  r ′′ry point  that is symmetric to ri over the particle surface 
and apply the identity 
( ) ( ))()()()()()( ii rPrPrPrPrPrP −′′+′′−′=−′  (79)
to the second integral in Eq. (77). Using Taylor expansion of P near  and the fact that  r ′′
irrrr −′≤′′−′  for iV∈′r  one can show that 
∫ ′′+≤ii cdcV ),( s36867r
ir ),(d rrGM , (80)
where the remaining integral can be proven to be equal to ),( iiV rL ∂− . The last prove left 
(see Eqs. (74) and (80)) is to demonstrate that ),( iiV rL ∂  is bounded by a constant. The only 
potential problem may come from the subsurface of iV∂  th  the particle surface 
(because it may be close to r
at is part of
ce is i). This subsurfa med planar. We will calculate the 
integral in Eq. 
(6) over the infinite plane rρrr +=−′ i uch that 0 s =⋅rρ . Then ρρn ±=′  
( ) 223
2d),plane.inf(
mm == ∫irL
+∂V r
r , (81)
which is bounded. The rest of the integral (over the part of the cube surface) is bounded by a 
constant, which is a manifestation of a more general fact that (by its definition) ), ir   ( iVL ∂
does not depend on the size but only on the shape of the volume. Finally we have 
69),( cV ii ≤∂ rL , 
which together with Eqs. (74), (78), and (80) prove Eq. (76). 
Using Eqs. (75) and (76) we obtain 
)ln()( 70
cllnc
+≤+≤ ∑ ∑∑∑ ∑ −hhh 7271
dccNd
, (83)
where we have changed the order of the summation in the double sum and split the 
summation over cubical shells for 
⎝∂∈∂∈∂∈
maxKl ≤  and . Then we have grouped everything 
otal esti
maxKl >
into one sum over boundary dipoles. Eqs. (41) and (62) were used in the last inequality. 
Combining Eqs. (63) and (83) one can obtain the t mate of the 
dh  for any scatterer: 
( ) ( )[ ]ddccddccNd lnln 7273245431 −+−≤h . (84)
Using Eq. (66) we immediately obtain the sam
δ de estimate for E . 
The derivation of the errors in the measured quantities is slightly modified compared to 
 (68) e changed to Section 2.D, by the presence of the shape errors. Eqs. (67) and  ar
( ) ( )~~ dcdcdcdc
,0 +≤+≤− ∑∑
EE φφ , (85)
( ) )( ( ) ( )ddccddccdddd lnln 2,0 −+−≤−Eφ . dcd δ 7655541
53≤ EEφ (86)
The second term in Eq. (85) comes from surface dipoles for which errors are the same order 
as the values themselves. Finally the generalization of Eq. (70) is 
( ) ( )yyccyyccy lnlnδ 797826059 −+−≤φ . (87)
The shape errors “reinforce” the surface errors (the linear term of discretization error), and 
although they both generally decrease with increasing size parameter x one may expect the 
linear term in Eq. (87) to be significant up to higher values of y than in Eq. (70). 
All the derivations in this section can in principle be extended to interfaces inside the 
particle, i.e. when a surface, which cannot be described exactly as a surface of a set of cubes, 
separates two regions where χ(r) varies smoothly. Two parts of the cubical dipole on the 
interface should be considered separately the same way as it was done above. This will 
however not change the main conclusion of this section – Eq.(87) – but only the constants. 
F.Different DDA formulations 
In this section we discuss how different DDA formulations modify the error estimates derived 
in Sections 2.D and 2.E. 
Most of the improvements of PP proposed in the literature are concerned with the self-
term – . They are the Radiative Reaction correction (RR),),( iiV rM
8 the Digitized Green’s 
Function (DGF),23 the formulation by Lakhtakia (LAK),26,27 the a1-term method,28,29 the 
Lattice Dispersion Relation (LDR),30 the formulation by Peltoniemi (PEL),31 and the 
Corrected LDR (CLDR).32 All of them provide an expression for  that is of order d),( iiV rM
(except for RR that is of order d3). For instance LDR is equivalent to 
( )[ ] iii ddSmbmbbV PrM 3223221 i)32(),( +++= , (88)
(remember that we assumed ) where b1=k 1, b2, b3 are numerical constants and S is a 
constant that depends only on the propagation and polarization vectors of the incident field. 
However, none of these formulations can exactly evaluate the integral in Eq. (39), because the 
variation of the electric field is not known beforehand (PEL solves this problem, but only for 
a spherical dipole). Therefore they (hopefully) decrease the constant in Eq. (40), thus 
decreasing the overall error in the measured quantities. However, these formulations are not 
expected to change the order of the error from d2 to some higher order. 
We do not analyze the improvements by Rahmani, Chaumet, and Bryant (RCB)33,34 and 
Surface Corrected LDR (SCLDR),17 as they are limited to certain particle shapes. 
There exist two improvements of the interaction term in PP: Filtered Coupled Dipoles 
(FCD)12 and Integration of Green’s Tensor (IT).35 A rigorous analysis of FCD errors is 
beyond the scope of this paper, but it seems that FCD is not designed to reduce the linear term 
in Eq. (63) that comes from the incomplete (non-symmetric) shells. This is because FCD 
employs sampling theory to improve the accuracy of the overall discretization for regular 
cubical grids. FCD does not improve the accuracy of a single ijG  calculation (approximation 
of an integral over one subvolume). 
IT, which numerically evaluates the integral in Eq. (12), has a more pronounced effect 
on the error estimate. Consider dipole j from l-th shell (incomplete) of dipole i, then 
.),(d),(maxd
)(),(d)(),(d
≤′′′+′∂′′≤
−′′′=−′′′
dlcrrcrrc
jijij
rrGrrG
PrPrrGPGrPrrG
 (89)
Here we have used Eq. (36) and Taylor expansion of Green’s tensor up to the first order. 
Eq. (89) states that the second term in Eq. (58) is completely eliminated and so is the linear 
term in Eqs. (69) and (70) (surface errors). Therefore convergence of DDA with IT for 
cubically shaped scatterers is expected to be purely quadratic (neglecting the logarithm). 
However, for non cubically shaped scatterers the linear term reappears, due to the shape 
errors. Both IT and FCD also modify the self-term, however the effect is basically the same as 
for the other formulations. 
Several papers aimed to reduce shape errors.10,11,36 The first one – Generalized Semi-
Analytical (GSA) method10 – modifies the whole DDA scheme, while the other two propose 
averaging of the susceptibility over the boundary dipoles. We will analyze here Weighted 
Discretization (WD) by Piller11 as probably the most advanced method to reduce shape errors 
available today. 
WD modifies the susceptibility and self-term of the boundary subvolume. We slightly 
modify the definition of the boundary subvolume used in Sections 2.B and 2.E to 
automatically take into account interfaces inside the scatterer. We define Vi to be always 
cubical, but with a possible interface inside. The particle surface, crossing the subvolume Vi, 
is assumed planar and divides the subvolume into two parts: the principal volume  
(containing the center) and the secondary volume  with susceptibilities ,  and 
electric fields ,  respectively. The electric fields are considered constant inside 
each part and related to each other via the boundary condition tensor 
iV ii χχ ≡
ii EE ≡
iT : 
iii ETE =
s . (90)
In WD the susceptibility of the boundary subvolume is replaced by an effective one, defined 
( ) 3ssppe dVV iiiiii TI χχχ += , (91)
which gives the correct total polarization of the cubical dipole. The effective self-term is 
directly evaluated starting from Eq. (4), considering χ and E constant inside each part, 
( ) ( ) iii
rrV ETrrGrrGrrGrrGrM
′−′′+′−′′= ∫∫ ss3ps3
),(),(d),(),(d),( χχ . (92)
Piller evaluated the integrals in Eq. (92) numerically.11 
To take a smooth variation of the electric field and susceptibility into account we define 
 ( r  is defined in Section )(s r ′′= χχ i ′′ 2.E) and iT  is calculated at the surface between ri and 
.  and r ′′ ii PP ≡
iiiiii ETEP
ssss χχ == . Then 
c rrPrP
−≤−′′
min)( 83
s , (93)
where we have assumed that Eq. (30) for χ(r) and E(r) is also valid in . siV
We start estimating errors of WD with  (cf. Eq. shijh (73)) 
( ) ( )∫∫ −′′′+−′′′=
s)0(3p)0(3sh )(),(d)(),(d
jijiij rr PGrPrrGPGrPrrGh , (94)
Using Taylor expansions of  near r)(rP ′ i and r ′′  in  and  correspondingly and Eq. piV
iV (93) 
one may obtain that the main contribution comes from the derivative of Green’s tensor, 
leading to (cf. Eq. (75)) 
h  (95)
iih  is the following (cf. Eq. (74)) 
( ) ( )
( ) ( ) ( )
.),(),(
),(d)(),(d)(),(d
),(),(),(d),(),(d
),(),(
pss3s3p3
ess3ps3
iiiiiii
iiiiiii
PrLErL
PPrrGPrPrrGPrPrrG
ErLPrrGrrGPrrGrrG
PrLrMh
−′′+−′′′+−′′′=
∂−′−′′+′−′′−
 (96)
The first two integrals can be easily shown to be dс85≤  (cf. Eq. (77)) and third one is 
transformed to L  the same way as in Eq. (80), thus 
iiiiiiiiiiii VVVdc ErLPrLPrLh
esspp
sh ),(),(),( χ∂−∂+∂+≤ , (97)
where the second term comes from the fact that averaged PL  is not the same as L  times 
averaged P. This error depends on the geometry of the interface inside Vi and generally is of 
order unity. For example, if the plane interface is described as ε+= izz , taking limit 0→ε  
gives the error ( )zii sp2 PP −π  (using Eq. (81)). Therefore WD does not principally improve 
the error estimate of  given by Eq. shiih (76), although it may significantly decrease the 
constant. On the other hand, since ),( p iiV rL ∂  and ),(
iiV rL ∂  can be (analytically) evaluated 
for a cube intersected by a plane, WD can be further improved to reduce the error in  to 
linear in d, which is a subject of future research. 
Proceeding analogous to the derivation of Eq. (83) one can obtain 
Ndcdccllnc
898887
)( ≤⎟⎟
++≤ ∑ ∑
h . (98)
It can be shown that for the scattering amplitude (Eq. (25)) the error estimate given by 
Eq. (85) can be improved, since WD correctly evaluates the zeroth order of value for the 
boundary dipoles, leading to 
( ) ( ) 291490551,0~~ dcdcdc
dd ≤+≤− ∑∑
EE φφ . (99)
In his original paper11 Piller did not specify the expression that should be used for Cabs. Direct 
application of the susceptibility provided by WD into Eq. (26) does not reduce the order of 
error when compared with the exact Eq. (24) (except when ), since they are not linear 
functions of the electric field. However, if we consider separately  and  (which is 
equivalent to replacing 
0s =iχ
( )( )∗⋅ iiiiV EEeIm χ  by 
2ss2pp )Im()Im( iiiiiii VV ETE χχ + ) the same 
estimate as in Eq. (99) can be derived for Cabs. 
Using Eqs. (98), (99), and the first part of Eq. (86) one can derive the final error 
estimate for WD: 
( ) ycyyccy 9429392 lnδ +−≤φ , (100)
where the constant before the linear term, as compared to Eq. (87), does not contain a 
logarithm and is expected to be significantly smaller, because several factors contributing to it 
are eliminated in WD. Although WD has a potential for improving, it does not seem feasible 
to completely eliminate the linear term in the shape error. The accuracy of evaluation of the 
interaction term over the boundary dipole (cf. Eq. (94)) can be improved by integration of 
Green’s tensor over  and  separately but that would ruin the block-Toeplitz structure of 
the interaction matrix and hinder the FFT-based algorithm for the solution of linear 
0 30 60 90 120 150 180
Scattering angle θ, deg
 cube kD=8
 discretized sphere kD=10
 sphere kD=3 (x10)
 sphere kD=10
 sphere kD=30
Fig. 2. S11(θ ) for all 5 test cases in logarithmic scale. The result for the kD = 3 sphere is multiplied by 
10 for convenience. 
equations.5 Since there is no comparable alternative to FFT nowadays, this method seems 
inapplicable. 
Minor modifications of the expression for Cabs are possible. Draine8 proposed a 
modification of Eq. (26) that was widely used afterwards and which was further modified by 
Chaumet et al.35 However, for many cases these expressions are equivalent and, even when 
they are not, the difference is of order d3, which is neglected in our error analysis. 
3. Numerical simulations 
A.Discrete Dipole Approximation 
The basics of the DDA method were summarized by Draine and Flatau.2 In this paper we use 
the LDR prescription for dipole polarizability,30 which is most widely used nowadays, e.g. in 
the publicly available code DDSCAT 6.1.4 We also employ dipole size correction8 for non-
cubically shaped scatterers to ensure that the cubical approximation of the scatterer has the 
correct volume; this is believed to diminish shape errors, especially for small scatterers.2 We 
use a standard discretization scheme as described in Section 2.E, without any improvements 
for boundary dipoles. It is important to note that all the conclusions are valid for any DDA 
implementation, but with a few changes for specific improvements as discussed in 
Section 2.F. 
Our code – Amsterdam DDA (ADDA) – is capable of running on a cluster of computers 
(parallelizing a single DDA computation), which allows us to use practically unlimited 
number of dipoles, since we are not limited by the memory of a single computer.37,38 We used 
a relative error of residual  as a stopping criterion for the iterative solution of the DDA 
linear system. Tests suggest that the relative error of the measured quantities due to the 
iterative solver is then  (data not shown) and hence can be neglected (total relative 
errors in our simulations are  – see Section 
810−<
710−<
56 1010 −− ÷> 0). More details about our code can 
be found in Paper 2. All DDA simulations were carried out on the Dutch national compute 
cluster LISA.39 
Table 1. Exact values of Qext for the 5 test cases. 
Particle Qext
kD = 8 cube 4.490 
discretized kD = 10 sphere 3.916 
kD = 3 sphere 0.753 
kD = 10 sphere 3.928 
kD = 30 sphere 1.985 
0.1 1
 slope = 0.77
 y = kd·m
 maximum
 θ = 0°
 θ = 45°
 θ = 90°
 θ = 135°
 θ = 180°
 slope = 0.95
Fig. 3. Relative errors of S11 at different angles θ and maximum over all θ versus y for (a) the kD = 8 
cube, (b) the cubical discretization of kD = 10 sphere. A log-log scale is used. A linear fit of maximum 
over θ errors is shown. (m = 1.5). 
B.Results 
We study five test cases: one cube with 8=kD , three spheres with , and a 
particle obtained by a cubical discretization of the 
30,10,3=kD
10=kD  sphere using 16 dipoles per D 
(total 2176 dipoles, x equal to that of a sphere; see detailed description in Paper 2). By D we 
denote the diameter of a sphere or the edge size of a cube. All scatterers are homogenous with 
. Although DDA errors significantly depend on m (see e.g. 5.1=m 14), we limit ourselves to 
one single value and study effects of size and shape of the scatterer. 
The maximum number of dipoles per D (nD) was 256. The values of nD that we used are 
of the form  (p is an integer), except for the discretized sphere, where all np2}7,6,5,4{ ⋅ D are 
0.1 1
0.1 1
0.01 0.1
y = kd·m
 slope 2.29 (c)
 slope 0.91
 maximum
 θ = 0°
 θ = 45°
 θ = 90°
 θ = 135°
 θ = 180°
 slope 1.05
Fig. 4. Same as Fig. 3 but for (a) kD = 3, (b) kD = 10, and (c) kD = 30 spheres. 
multiples of 16 (this is required to exactly describe the shape of the particle composed from a 
number of cubes). The minimum values for nD were 8 for the 3=kD  sphere, 16 for the cube, 
the  sphere, and the discretized sphere, and 40 for the 10=kD 30=kD  sphere. 
All the computations use a direction of incidence parallel to one of the principal axes of 
the cubical dipoles. The scattering plane is parallel to one of the face of the cubical dipoles. In 
this paper we show results only for the extinction efficiency Qext (for incident light polarized 
parallel to one of the principal axes of the cubical dipoles) and phase function S11(θ ) as the 
0.01 0.1 1
 slope = 0.89
 cube kD=8
 discretized sphere kD=10
 sphere kD=3
 sphere kD=10
 sphere kD=30
y=kd·m  
Fig. 5. Relative errors of Qext versus y for all 5 test cases. A log-log scale is used. A linear fit through 5 
finest discretizations of kD = 3 sphere is shown. 
most commonly used in applications. However, the theory applies to any measured quantity. 
For instance, we have also confirmed it for other Mueller matrix elements (data not shown). 
Exact results of S11(θ ) for all 5 test cases are shown in Fig. 2. For spheres this is the 
result of Mie theory (the relative accuracy of the code we used24 is at least ) and for the 
cube and discretized sphere an extrapolation over the 5 finest discretizations (the 
extrapolation technique is presented in Paper 2, together with all details of obtaining these 
results, including their estimated errors). We use such ‘exact’ results because analytical theory 
is unavailable for these shapes and because errors of the best discretization are larger than that 
of the extrapolation. Their use as references for computing real errors (difference between the 
computed and the exact value) of single DDA calculations is justified because all these real 
errors are significantly larger than the errors of the references themselves (see Paper 2; in 
general, real errors obtained this way have an uncertainty of reference error). Exact values of 
610−<
ext for all test cases are presented in Table 1. 
In the following we show the results of DDA convergence. Fig. 3 and Fig. 4 present 
relative errors (absolute values) of S11 at different angles θ and maximum error over all θ 
versus y in log-log scale. In many cases the maximum errors are reached at exact 
backscattering direction, then these two sets of points overlap. Deep minima that happen at 
intermediate values of y for some values of θ (and also sometimes for Qext – Fig. 5) are due to 
the fact that the differences between simulated and reference values change sign near these 
values of y (see Paper 2 for detailed description of this behavior). The solid lines are linear fits 
to all or some points of maximum error. The slopes of these lines are depicted in the figures. 
Fig. 5 shows relative errors of Qext for all 5 studied cases in log-log scale. A linear fit through 
the 5 finest discretizations of the 3=kD  sphere is shown. More results of these numerical 
simulations are presented in Paper 2. 
4. Discussion 
Convergence of DDA for cubically shaped particles (Fig. 3) shows the following trends. All 
curves have linear and quadratic parts (the non-monotonic behavior of errors for some θ are 
also a manifestation of the fact that signed difference can be approximated by a sum of linear 
and quadratic terms that have different signs). The transition between these two regimes 
occurs at different y (which indicates the relative importance of linear and quadratic 
coefficients). While for maximum errors that are close to those of the backscattering direction 
the linear term is significant for larger y, it is much smaller and not significant in the whole 
range of y studied for side scattering ( °= 90θ ). R
for l
esults of DDA convergence for spheres (Fig. 
4) show a different behavior for different sizes. Errors for the small ( 3=kD ) sphere converge 
purely linear (except for small deviation of errors of )90(11 °S   values of y). Similar 
results are obtained for the 10=kD  sphere, but with significant oscillations superimposed 
upon the general trend. Convergence for the large ( 30
=kD ) sphere is quadratic or even faster 
in the range of y studied, also with significant oscillations. 
Comparing Fig. 3 and Fig. 4 (especially Fig. 3(b) and Fig. 4(b) showing results for 
almost the same particles) one can deduce the following differences in DDA convergence for 
cubically and non-cubically shaped scatterers. The linear term for cubically shaped scatterers 
is significantly smaller, resulting in smaller total errors, especially for small y. All these 
conclusions, together with the size dependence of the significance of the linear term in the 
total errors, are in perfect agreement with the theoretical predictions made in Sections 2.D and 
2.E. Errors for non-cubically shaped particles exhibit quasi-random oscillations that are not 
present for cubically shaped particles. This can be explained by the sharp variations of shape 
errors with changing y (discussed in details in Paper 2). Oscillations for the  sphere 
Fig. 4(a)) are very small (but still clearly present), which is due to the small size of the 
particle and hence featurelessness of its light scattering pattern – the surface structure is not 
that important and one may expect rather small shape errors. Results for Qext (Fig. 5) fully 
support the conclusions. Errors of Qext for the large sphere at small values of y are 
unexpectedly smaller than for smaller spheres. This feature requires further study before 
making any firm conclusions, however there is definitely no similar tendency for S11(θ ) (cf. 
Fig. 4). 
We have also studied a  porous cube that was obtained by dividing a cube into 
27 smaller cubes and then removing randomly 9 of them. All the conclusions are the same as 
those reported for the cube, but with slightly higher overall errors (data not shown). 
In this paper we have used a traditional DDA formulation2 for numerical simulations. 
However, as we showed in Section 2.F several modern improvements of DDA (namely IT 
and WD) should significantly change its convergence behavior. IT should completely 
eliminate the linear term for cubically shaped scatterers, which should improve the accuracy 
especially for small y. WD should significantly decrease shape and hence total errors for non-
cubically shaped particles, moreover it should significantly decrease the amplitude of quasi-
random error oscillations because it takes into account the location of the interface inside the 
boundary dipoles. Numerical testing of DDA convergence using IT and WD is a subject of a 
future study. 
5. Conclusion 
To the best of our knowledge, we conducted for the first time a rigorous theoretical 
convergence analysis of DDA. In the range of DDA applicability ( 2<kd ) errors are bounded 
by a sum of a linear and quadratic term in the discretization parameter y; the linear term is 
significantly smaller for cubically than for non-cubically shaped scatterers. Therefore for 
small y errors for cubically shaped particles are much smaller than for non-cubically shaped. 
The relative importance of the linear term decreases with increasing size, hence convergence 
of DDA for large enough scatterers is quadratic in the common range of y. All these 
conclusions were verified by extensive numerical simulations. 
Moreover, these simulations showed that errors are not only bounded by a quadratic 
function (as predicted in Section 2), but actually can be (with good accuracy) described by a 
quadratic function of y. This fact provides a basis for the extrapolation technique presented in 
Paper 2. 
Our theory predicts that modern DDA improvements (namely IT and WD) should 
significantly change the convergence of DDA, however numerical testing of these predictions 
is left for future research. 
Acknowledgements
We thank Gorden Videen and Michiel Min for valuable comments on earlier version of this 
manuscript and anonymous reviewer for helpful suggestions. Our research is supported by the 
NATO Science for Peace program through grant SfP 977976. 
References 
 1.  E. M. Purcell and C. R. Pennypacker, "Scattering and adsorption of light by nonspherical dielectric 
grains," Astrophys. J. 186, 705-714 (1973). 
 2.  B. T. Draine and P. J. Flatau, "Discrete-dipole approximation for scattering calculations," J. Opt. Soc. 
Am. A 11, 1491-1499 (1994). 
 3.  B. T. Draine, "The discrete dipole approximation for light scattering by irregular targets," in Light 
Scattering by Nonspherical Particles, Theory, Measurements, and Applications, M. I. Mishchenko, J. W. 
Hovenier, and L. D. Travis, eds. (Academic Press, New York, 2000), pp. 131-145. 
 4.  B. T. Draine and P. J. Flatau, "User guide for the discrete dipole approximation code DDSCAT 6.1," 
http://xxx.arxiv.org/abs/astro-ph/0409262 (2004). 
 5.  J. J. Goodman, B. T. Draine, and P. J. Flatau, "Application of fast-Fourier-transform techniques to the 
discrete-dipole approximation," Opt. Lett. 16, 1198-1200 (1991). 
 6.  G. C. Hsiao and R. E. Kleinman, "Mathematical foundations for error estimation in numerical solutions of 
integral equations in electromagnetics," IEEE Trans. Ant. Propag. 45, 316-328 (1997). 
 7.  K. F. Warnick and W. C. Chew, "Error analysis of the Moment Method," IEEE Ant. Prop. Mag. 46, 38-53 
(2004). 
 8.  B. T. Draine, "The discrete-dipole approximation and its application to interstellar graphite grains," 
Astrophys. J. 333, 848-872 (1988). 
 9.  J. I. Hage, J. M. Greenberg, and R. T. Wang, "Scattering from arbitrarily shaped particles - theory and 
experiment," Appl. Opt. 30, 1141-1152 (1991). 
 10.  F. Rouleau and P. G. Martin, "A new method to calculate the extinction properties of irregularly shaped 
particles," Astrophys. J. 414, 803-814 (1993). 
 11.  N. B. Piller, "Influence of the edge meshes on the accuracy of the coupled-dipole approximation," Opt. 
Lett. 22, 1674-1676 (1997). 
 12.  N. B. Piller and O. J. F. Martin, "Increasing the performance of the coupled-dipole approximation: A 
spectral approach," IEEE Trans. Ant. Propag. 46, 1126-1137 (1998). 
 13.  N. B. Piller, "Coupled-dipole approximation for high permittivity materials," Opt. Comm. 160, 10-14 
(1999). 
 14.  A. G. Hoekstra, J. Rahola, and P. M. A. Sloot, "Accuracy of internal fields in volume integral equation 
simulations of light scattering," Appl. Opt. 37, 8482-8497 (1998). 
 15.  S. D. Druger and B. V. Bronk, "Internal and scattered electric fields in the discrete dipole approximation," 
J. Opt. Soc. Am. B 16, 2239-2246 (1999). 
 16.  Y. L. Xu and B. A. S. Gustafson, "Comparison between multisphere light-scattering calculations: 
Rigorous solution and discrete-dipole approximation," Astrophys. J. 513, 894-909 (1999). 
 17.  M. J. Collinge and B. T. Draine, "Discrete-dipole approximation with polarizabilities that account for both 
finite wavelength and target geometry," J. Opt. Soc. Am. A 21, 2023-2028 (2004). 
 18.  M. A. Yurkin, V. P. Maltsev, and A. G. Hoekstra, "Convergence of the discrete dipole approximation. II. 
An extrapolation technique to increase the accuracy," J. Opt. Soc. Am. A 23, 2592-2601 (2006). 
 19.  A. Lakhtakia, "Strong and weak forms of the method of moments and the coupled dipole method for 
scattering of time-harmonic electromagnetic-fields," Int. J. Mod. Phys. C 3, 583-603 (1992). 
 20.  F. M. Kahnert, "Numerical methods in electromagnetic scattering theory," J. Quant. Spectrosc. Radiat. 
Transf. 79, 775-824 (2003). 
 21.  C. P. Davis and K. F. Warnick, "On the physical interpretation of the Sobolev norm in error estimation," 
Appl. Comp. ElectroMagn. Soc. J. 20, 144-150 (2005). 
http://xxx.arxiv.org/abs/astro-ph/0409262
 22.  A. D. Yanghjian, "Electric dyadic Green's function in the source region," IEEE Proc. 68, 248-263 (1980). 
 23.  G. H. Goedecke and S. G. O'Brien, "Scattering by irregular inhomogeneous particles via the digitized 
Green's function algorithm," Appl. Opt. 27, 2431-2438 (1988). 
 24.  C. F. Bohren and D. R. Huffman, Absorption and scattering of Light by Small Particles, (Wiley, New 
York, 1983). 
 25.  J. Rahola, "On the eigenvalues of the volume integral operator of electromagnetic scattering," SIAM J. 
Sci. Comp. 21, 1740-1754 (2000). 
 26.  A. Lakhtakia and G. W. Mulholland, "On 2 numerical techniques for light-scattering by dielectric 
agglomerated structures," J. Res. Nat. Inst. Stand. Technol. 98, 699-716 (1993). 
 27.  J. I. Hage and J. M. Greenberg, "A model for the optical-properties of porous grains," Astrophys. J. 361, 
251-259 (1990). 
 28.  C. E. Dungey and C. F. Bohren, "Light-scattering by nonspherical particles - a refinement to the coupled-
dipole method," J. Opt. Soc. Am. A 8, 81-87 (1991). 
 29.  H. Okamoto, "Light scattering by clusters: the a1-term method," Opt. Rev. 2, 407-412 (1995). 
 30.  B. T. Draine and J. J. Goodman, "Beyond clausius-mossotti - wave-propagation on a polarizable point 
lattice and the discrete dipole approximation," Astrophys. J. 405, 685-697 (1993). 
 31.  J. I. Peltoniemi, "Variational volume integral equation method for electromagnetic scattering by irregular 
grains," J. Quant. Spectrosc. Radiat. Transf. 55, 637-647 (1996). 
 32.  D. Gutkowicz-Krusin and B. T. Draine, "Propagation of electromagnetic waves on a rectangular lattice of 
polarizable points," http://xxx.arxiv.org/abs/astro-ph/0403082 (2004). 
 33.  A. Rahmani, P. C. Chaumet, and G. W. Bryant, "Coupled dipole method with an exact long-wavelength 
limit and improved accuracy at finite frequencies," Opt. Lett. 27, 2118-2120 (2002). 
 34.  A. Rahmani, P. C. Chaumet, and G. W. Bryant, "On the importance of local-field corrections for 
polarizable particles on a finite lattice: Application to the discrete dipole approximation," Astrophys. J. 
607, 873-878 (2004). 
 35.  P. C. Chaumet, A. Sentenac, and A. Rahmani, "Coupled dipole method for scatterers with large 
permittivity," Phys. Rev. E 70, 036606 (2004). 
 36.  K. F. Evans and G. L. Stephens, "Microwave radiative-transfer through clouds composed of realistically 
shaped ice crystals .1. Single scattering properties," J. Atmos. Sci. 52, 2041-2057 (1995). 
 37.  A. G. Hoekstra, M. D. Grimminck, and P. M. A. Sloot, "Large scale simulations of elastic light scattering 
by a fast discrete dipole approximation," Int. J. Mod. Phys. C 9, 87-102 (1998). 
 38.  M. A. Yurkin, K. A. Semyanov, P. A. Tarasov, A. V. Chernyshev, A. G. Hoekstra, and V. P. Maltsev, 
"Experimental and theoretical study of light scattering by individual mature red blood cells with scanning 
flow cytometry and discrete dipole approximation," Appl. Opt. 44, 5249-5256 (2005). 
 39.  "Description of the national compute cluster Lisa," http://www.sara.nl/userinfo/lisa/description/ (2005). 
http://xxx.arxiv.org/abs/astro-ph/0403082
http://xxx.arxiv.org/abs/astro-ph/0409262
  /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 performed a rigorous theoretical convergence analysis of the discrete
dipole approximation (DDA). We prove that errors in any measured quantity are
bounded by a sum of a linear and quadratic term in the size of a dipole d, when
the latter is in the range of DDA applicability. Moreover, the linear term is
significantly smaller for cubically than for non-cubically shaped scatterers.
Therefore, for small d errors for cubically shaped particles are much smaller
than for non-cubically shaped. The relative importance of the linear term
decreases with increasing size, hence convergence of DDA for large enough
scatterers is quadratic in the common range of d. Extensive numerical
simulations were carried out for a wide range of d. Finally we discuss a number
of new developments in DDA and their consequences for convergence.

Introduction 
The discrete dipole approximation (DDA) is a well-known method to solve the light 
scattering problem for arbitrary shaped particles. Since its introduction by Purcell and 
Pennypacker1 it has been improved constantly. The formulation of DDA summarized by 
Draine and Flatau2 more than 10 years ago is still most widely used for many applications,3 
partly due to the publicly available high-quality and user-friendly code DDSCAT.4 Although 
modern improvements of DDA (as discussed in detail in Section 2.F) exist, they are still in the 
research stage because they are not widely used in real applications. 
DDA directly discretizes the volume of the scatterer and hence is applicable to arbitrary 
shaped particles. However, the drawback of this discretization is the extreme computational 
complexity of DDA, although it is significantly decreased by advanced numerical 
techniques.2,5 That is why the usual application strategy for DDA is “single computation”, 
where a discretization is chosen based on available computational resources and some 
empirical estimates of the expected errors.3,4 These error estimates are based on a limited 
number of benchmark calculations3 and hence are external to the light scattering problem 
under investigation. Such error estimates have evident drawbacks, however no better 
alternative is available. Some results of analytical analysis of errors in computational 
electromagnetics are known, e.g. 6,7, however they typically consider the surface integral 
equations. To the best of our knowledge, such analysis has not been done for volume integral 
equations (such as DDA). 
Usually errors in DDA are studied as a function of the size parameter of the scatterer x 
(at a constant or few different total numbers of dipoles N), e.g. 2,8. Only a small number of 
papers directly present errors versus discretization parameter (e.g. d – the size of a single 
dipole).9-17 The range of d typically studied in those papers is limited to a 5 times difference 
between minimum and maximum values, with the exception of two papers11,12 where it is 15 
times. Those plots of errors versus discretization parameter are always used to illustrate the 
performance of a new DDA formulation and compare it with others. No conclusions about the 
convergence properties of DDA, as a function of d, have been made based on these plots. To 
our knowledge, no theoretical analysis of DDA convergence has been performed, but only a 
few limited empirical studies have appeared in the literature. 
In this paper we perform a theoretical analysis of DDA convergence when refining the 
discretization (Section 2). We derive rigorous theoretical bounds on the error in any measured 
quantity for any scatterer. In Section 3 we present extensive numerical results of DDA 
computations for 5 different scatterers using many different discretizations. These results are 
discussed in Section 4 to support conclusions of the theoretical analysis. We formulate the 
conclusions of the paper in Section 5. In a follow-up paper18 (which from now on we refer to 
as Paper 2) the theoretical convergence results are used for an extrapolation technique to 
increase the accuracy of DDA computations.  
2. Theoretical analysis 
In this section we analyze theoretically the errors of DDA computations. We formulate the 
volume integral equation for the internal electric field and its operator counterpart in 
Section 2.A and its discretization in Section 2.B. Section 2.C contains integral and discretized 
formulae for measured quantities that are the final goal of any light scattering simulation. We 
derive the main results in Section 2.D, where we consider errors of the traditional DDA 
formulation2 without shape errors, which are considered separately in Section 2.E. Finally in 
Section 2.F we discuss some recent DDA improvements from the viewpoint of our 
convergence theory. 
A.Integral Equation 
Throughout this paper we assume the )iexp( tω−  time dependence of all fields. The scatterer 
is assumed dielectric but not magnetic (magnetic permittivity 1=μ ), and the electric 
permittivity is assumed isotropic (non-isotropic permittivity will significantly complicate the 
derivations but will not principally change the main conclusion of Section 2 – Eqs. (70) and 
(87)). 
The general form of the integral equation governing the electric field inside the 
dielectric scatterer is the following:19,20 
)()(),(),()()(),(d)()( 00
rErrLrMrErrrGrErE χχ VVr
∂−+′′′′+= ∫ , (1)
where Einc(r), 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 (more general – the volume, which 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 
−=∇∇+=′
)()(ˆˆ),(
kRgRgk IIIrrG , (2)
where I  is the identity dyadic, ck ω=  – free space wave vector, rrR ′−= , R=R , and 
 is a dyadic defined as  (μ, ν are Cartesian components of the vector or 
tensor), and g(R) is the scalar Green’s function 
RR ˆˆ νμμν RRRR =ˆˆ
iexp(
)( = . (3)
M is the following integral associated with the finiteness of the exclusion volume V0
( )∫ ′−′′′′=
)()(),()()(),(d),( s30
rV rErrrGrErrrGrM χχ , (4)
where ),(s rrG ′  is the static limit ( ) of 0→k ),( rrG ′ : 
−−=∇∇=′ 23
11ˆˆ),(
IrrG . (5)
L  is the so-called self-term dyadic: 
rV rL , (6)
where  is an external (as viewed from r) normal to the surface ∂Vn′ˆ 0 at point r'. 
Eq. (1) can be rewritten in operator form as follows 
inc~~~ EEA =⋅ , (7)
where ( 311 )~ CE →=∈ VLH  – functions from V to C3 that have finite L1-norm, 2inc~ H∈E  – 
subspace of H1 containing all functions that satisfy Maxwell equations in free space. A
 is a 
linear operator . Although the Sobolev norm is physically more sound (based on 
the finiteness of energy of the electric field),
21: HH →
6,21 we use the L1-norm. A detailed discussion of 
all assumptions made for the electric field is performed in Section 2.D. 
B.Discretization 
To solve Eq. (1) numerically a discretization is done in the following way.20 Let , 
 for . N denotes the number of subvolumes (dipoles). Assuming  and 
choosing , Eq. 
/0=ji VV I ji ≠ iV∈r
iVV =0 (1) can be rewritten as 
)()(),(),()()(),(d)()( 3inc rErrLrMrErrrGrErE χχ ii
∂−+′′′′+= ∑ ∫
. (8)
The set of Eq. (8) (for all i) is exact. Further one fixed point ri inside each Vi (its center) is 
chosen and  is set. irr =
The usual approximation20 is considering E and χ constant inside each subvolume: 
iiiii V∈==== rrrErErE for)()(,)()( χχχ . (9)
Eq. (8) can then be rewritten as 
( ) iiii
jjjijii V ELMEGEE χχ −++= ∑
inc , (10)
where , )(incinc ii rEE = ),( iii V rLL ∂= , 
( )∫ ′−′′=
iii r ),(),(d
s3 rrGrrGM , (11)
∫ ′′=
ij rV
1 3 rrGG . (12)
A further approximation, which is used in almost all formulations of DDA, is 
),()0( jiij rrGG = . (13)
This assumption is made implicitly by all formulations that start by replacing the scatterer 
with a set of point dipoles, as was done originally by Purcell and Pennypacker.1 For a cubical 
(as well as spherical) cell Vi with ri located at the center of the cell, iL  can be calculated 
analytically yielding22 
=i . (14)
Eq. (10) together with Eqs. (13) and (14) and completely neglecting iM  is equivalent to 
the original DDA by Purcell and Pennypacker (PP).1 The diagonal terms in Eq. (10) are then 
equivalent to the well-known Clausius-Mossotti (CM) polarizability for point dipoles. 
Modifications introduced by other DDA prescriptions are discussed in Section 2.F. 
In matrix notation Eq. (10) reads 
ddd ,incEEA = , (15)
where Ed, Einc,d are elements of  (vectors of size N where each element is a 3D complex 
vector) and 
( )N3C
dA  is a N×N matrix where each element is a 3×3 tensor. d is the size of one 
dipole. In operator notation Eq. (8) (for irr = ) is as follows 
( ) diii ,incinc )(~)(~~ ErErEA == , (16)
We define the discretization error function as 
( ) ( )iddidi ,0)(~~ EArEAh −= , (17)
where E0,d is the exact field at the centers of the dipoles – )(
i rEE = , in contrast to E
d that 
is only an approximation obtained from solution of Eq. (15) (here we neglect the numerical 
error that appears from the solution of Eq. (15) itself, which is acceptable if this error is 
controlled to be much less than other errors). Comparing Eqs. (15) and (17) one can 
immediately obtain the error in internal fields due to discretization δEd: 
( ) ddddd hAEEE 1,0δ −−=−= . (18)
C.Measured quantities 
After having determined the internal electric fields, scattered fields and cross sections can be 
calculated. Scattered fields are obtained by taking the limit ∞→r  of the integral in Eq. (1) 
(see e.g. 23) 
)iexp(
)(sca nFrE
= , (19)
where rrn =  is the unit vector in the scattering direction, and F is the scattering amplitude: 
∑∫ ′′⋅′−′−−=
krnnk )()()iexp(d)ˆˆ(i)( 33 rErnrInF χ . (20)
All other differential scattering properties, such as the amplitude and Mueller scattering 
matrices, and asymmetry parameter >< θcos  can be easily derived from F(n), calculated for 
two incident polarizations.24 We consider an incident polarized plane wave: 
)iexp()( 0inc rkerE ⋅= , (21)
where , a is direction of incidence, and ak k= 10 =e  is assumed. The scattering and 
extinction cross sections (Csca, Cext) are derived from the scattering amplitude:23 
∫ Ω= nFkC , (22)
( )∗⋅= 02ext )(Re
, (23)
where * denote complex conjugation. The expression for absorption cross section (Cabs) 
directly uses the internal fields:23 
( )∑ ∫ ′′′=
abs )()(Imd4 rErχπ , (24)
Since only values of the internal field in the centers of dipoles are known, Eqs. (20) and 
(24) are approximated by (PP) 
∑ ⋅−χ−−=
iii kVnnk )iexp()ˆˆ(i)(
3 nrEInF , (25)
iiiVkC
abs )Im(4 Eχπ . (26)
Corrections to Eq. (26) are discussed in Section 2.F. 
Both Eqs. (20) (for each component) and (24) can be generalized as ( )E~~φ  (a functional 
that is not necessarily linear), which is approximated as: 
( ) ( ) ddd φφφ δ~~ += EE , (27)
where ( )dd Eφ  corresponds to Eqs. (25) or (26) respectively, and the error δφ d consists of  two 
parts: 
( ) ( )[ ] ( ) ( )[ ]ddddddd EEEE φφφφφ −+−= ,0,0~~δ . (28)
The first one comes from discretization (similar to Eq. (17)), and the second from errors in the 
internal fields. 
D.Error analysis 
In this section we perform error analysis for the PP formulation of DDA. Improvements of 
DDA are further discussed in Section 2.F. 
We assume cubical subvolumes with size d. We also assume that the shape of the 
particle is exactly described by these cubical subvolumes (we call this cubically shaped 
scatterer). Moreover, χ is a smooth function inside V (exact assumptions on χ are formulated 
below). An extension of the theory to shapes that do not satisfy these conditions is presented 
in Section 2.E. If there are several regions with different values of χ (smooth inside each 
region), the analysis is still valid but interfaces inside V should be considered the same way as 
the outer boundary of V. We further fix the geometry of the scattering problem and incident 
field. Therefore we will be interested only in variation of discretization (which is 
characterized by the single parameter – d); for reasons that will become clear in the sequel, we 
assume that  (this bound is not limiting since otherwise DDA is generally 
inapplicable
We switch to dimensionless parameters by assuming 1=k , which is equivalent to 
measuring all the distances in units of k1 . The unit of electric field can be chosen arbitrary 
but constant. In all further derivations we will use two sets of constants: γi and ci. γ1-γ13 are 
basic constants that do not depend on the discretization d, but do depend directly upon all 
other problem parameters – size parameter eqkRx =  (Req – volume-equivalent radius), m, 
shape, and incident field – or some of them. On the contrary, c1-c94 are auxiliary values that 
either are numerical constants or can be derived in terms of constants γi. Although the 
dependencies of ci on γi are not explicitly derived in this paper, one can easily obtain them 
following the derivations of this section. That is the main motivation for using such vast 
amount of constants instead of an “order of magnitude” formalism. However, such explicit 
derivation has limited application because, as we will see further, constants in the final result 
depend upon almost all basic constants. Qualitative analysis of these dependencies will be 
performed in the end of this section. It should be noted that the main theoretical results 
concerning DDA convergence (boundedness of errors by a quadratic function, cf. Eq (70)) 
can be formulated and applied without consideration of any constants (which is simpler). 
However our full derivation enables us to make additional conclusions related to the behavior 
of specific error terms. 
The total number of dipoles used to discretize the scatterer is 
−= dN γ . (29)
We assume that the internal field E
 is at least four times differentiable and all these 
derivatives are bounded inside V 
65432 )(,)(,)(,)(,)( γγγγγ τρνμρνμνμμ ≤∂∂∂∂≤∂∂∂≤∂∂≤∂≤ rErErErErE  
for V∈r  and τρνμ ,,,∀ . 
This assumption is acceptable since there are no interfaces inside V, therefore E
 should be a 
smooth function. .  denotes the Euclidian (L2) norm, which is used for all 3D objects: vectors 
and tensors. We use the L1-norm, 
. , for N-dimensional vectors and matrices as well as for 
functions and operators. Eq. (30) immediately implies that . We require that χ 
satisfies Eq. 
~ 1 VL∈E
(30) with constants γ7-γ11. Further we will state an estimate for the norm of 
)(RG  and its derivatives. One can easily obtain from Eq. (2) that for 1>R  )(RG  satisfies 
Eq. (30) (with constants c1-c5), while for 2≤R  
6 )(,)(,)(,)(
−−−− ≤∂∂∂≤∂∂≤∂≤ RcRcRcRc RGRGRGRG ρνμνμμ , 
−≤∂∂∂∂ RcRGτρνμ  for τρνμ ,,,∀ . 
Next we state two auxiliary facts that will be used later. Let Vc be a cube with size d and 
with its center at the origin and f(r) a four times differentiable function inside Vc. Then 
)(max)()(d
fdcffr
d cVV
∈,∫ , (32)
( ) )(max)(
)()(d
d cVV
∂∂∂∂+∇≤−
=∫ . (33)
Eqs. (32) and (33) are the corollary of expanding f into Taylor series. Odd orders of the Taylor 
expansion vanish because of cubical symmetry. 
Our first goal is to estimate 
dh . Starting from Eq. (17) we write  as dih
),()(),(d )0(33 ii
i Vdr
rMPGrPrrGh +
−′′′= ∑ ∫
, (34)
where we have introduced the polarization vector for conciseness 
)()()( rErrP χ= , )( ii rPP = . (35)
It is evident that  also satisfies Eq. )(rP (30) (with constants c13-c17). We start by estimating 
),( iiV rM . Substituting a Taylor expansion of  )(rP
( ) ( )( )∑∑ ∂∂+∂+=
ρρ τρ ),,(
)()()( RrP0P0PRP RRR , (36)
where , into Eq μμ Rr ≤≤
~0 (4) gives 
( ) ( )(∫ ∑∫ ∂∂+−=
ii RRRRV
τρτρ τρ ),,(
)()(d),( 3s3 RrPRGPRGRGrM ) . (37)
The norms of these two terms can be estimated as 
( ) 2183s3 )(d3
)()(d dcRRgR
i ≤=− ∫∫ PIPRGRG , (38)
( ) 21923153 )(d3)),,(~()(d dcRRcRRR
≤≤∂∂ ∫∫ ∑ RGRrPRG
τρτρ τρ . (39)
Eq. (38) follows directly from the definitions in Eqs. (2), (5). To derive Eq. (39) we used 
Eq. (31) and the fact that 23RRR ≤∑
τρ . Finally, Eqs. (37)-(39) lead to 
20),( dcV ii ≤rM . (40)
To estimate the sum in Eq. (34) we consider separately three cases: 1) dipole j lies in a 
complete shell of dipole i (we define the shell below); 2) j lies in a distant shell of dipole i – 
1>−= ijijR rr ; 3) all j that fall between the first two cases (see Fig. 1). We define the first 
shell (S1(i)) of a cubical dipole as a set of dipoles that touch it (including touching in one point 
only). The second shell (S2(i)) is a set of dipoles that touch the outer surface of the first shell, 
and so on. The l-th shell (Sl(i)) is then a set of all dipoles that lie on the boundary of the cube 
with size  and center coinciding with the center of the original dipole. We call a shell 
complete if all its elements lie inside the volume of the scatterer V. A shell is called a distant 
dl )12( +
Kmax K(i) 
(2) (3) 
vacuum scatterer 
Fig. 1. Partition of the scatterer’s volume into three regions relative to dipole i. 
shell if all its elements satisfy , i.e. if its order 1>ijR [ ]dKl 1max => . Let K(i) be the order of 
the first incomplete shell, which is an indicator of how close dipole i is to the surface. We 
demand  to separate cases (1) and (2) described above. All j that fall in the third 
case satisfy  (the exact value of this constant – slightly larger than 
max)( KiK ≤
2<ijR 3  – depends on 
d). The number of dipoles in a shell Sl (which can be incomplete) – ns(l) – can be estimated as 
33 )12()12()( lclllns ≤−−+≤ . (41)
The sum of the error over all dipoles that lie in complete shells is then 
∑ ∑ ∫
)0(33 )(),(d
l iSj
dr PGrPrrG , (42)
Since each shell in Eq. (42) is complete it can be divided into pairs of dipoles that are 
symmetric over the center of the shell (j and –j). For convenience we set . The inner 
sum in Eq. 
0r =i
(42) can then be rewritten as 
( ) ( )∑ ∫
+−′−+′′′
)0(33 )()()(d
dr PPGrPrPrG , (43)
Further we introduce the auxiliary function 
( ) )()()(
)( 0PrPrPru −′−+′=′ , (44)
which satisfies the following inequalities (follows from Eq. (30) for P(r) and Taylor series) 
22 )(,)(,)( crcrc ≤∂∂≤∂≤ rururu νμμ  for νμ,∀ . (45)
Then Eq. (43) is equivalent to 
V l jl j
drdr PGrGuGrurG ∑ ∫∑ ∫
)0(33
)0(33 )(d)()(d , (46)
where . To estimate the first term we apply Eq. )( jj ruu = (32) to the whole function under 
the integral. Using Eqs. (31) and (45) one may obtain 
( ) 325)()(max −
≤′′∂∂ ij
, (47)
and hence 
)0(33 )()(d
V ll j
ldcRdcdr uGrurG , (48)
where we have used Eq. (41) and  for ldRij ≥ )(iSj l∈ . 
It is straightforward to show that 
∑ ∫∑ ∫
′′=′′
3 )(d
iSj ViSj V l jl j
rgrr IrG , (49)
)0( )(
RgIG . (50)
The derivation is based upon Eq. (2) and the equivalence I
 in all sums and integrals 
that satisfy cubical symmetry. Then second part of Eq. (46) is transformed to 
)0(33 )()(d)(d −
+≤−′′≤
−′′ ∑ ∫∑ ∫ ldcldcRgdrgrcdr
V l jl j
PGrG , (51)
where we apply Eq. (33) to derive the second inequality and use the identity  
and the following inequalities 
)()(2 rgrg −=∇
( ) 532131 )(, −− ≤∂∂∂∂≤ RcRgRcRg τρνμ  for τρνμ ,,,∀ . (52)
Substituting Eqs. (48) and (51) into Eq. (42) one can obtain 
( ) 23433
)0(33 )(ln)(),(d diKccdr
l iSj
−′′′∑ ∑ ∫
PGrPrrG , (53)
using the fact that . 1)( ≤diK
We now consider the second part of the sum in Eq. (34) (where ). We first apply 
1>ijR
(32), then use Eq. (30) for P(r) and )(rG , and finally invoke Eq. (29): 
)0(33 )(),(d dcdNcdcdr
ijij j
i ≤≤≤⎟
−′′′ ∑∑ ∫
PGrPrrG . (54)
To analyze the third part of the sum in Eq. (34) we again sum over shells, however since 
they are incomplete we cannot use symmetry considerations. We apply Eq. (33) to the whole 
function under the integral and proceed analogous to the derivation of Eq. (51). Using the 
identity 
)()(2 rGrG −=∇ , (55)
(since we have assumed ) we obtain 1=k
( ) 4372 )()( −= ≤∇ ijRcijRrrPrG , (56)
( ) 738)()(max −
≤′′∂∂∂∂ ij
, (57)
which leads to 
)0(33 )(),(d −−
−′′′∑ ∫ lcdlcdr
PGrPrrG , (58)
and then analogous to Eq. (53): 
)()()(),(d 442
)( )(
)0(33
iKcidKcdr
iKl iSj
−′′′∑ ∑ ∫ PGrPrrG . (59)
Collecting Eqs. (40), (53), (54), (59) we finally obtain 
( ) 24443442141 )(ln)()( diKcciKcidKcdi +++≤ −−h . (60)
Then 
( ) ( )442141
max4443
)(ln −−
+++≤= ∑∑ KcdKcKnNdKcc
d hh , (61)
where n(K) is the number of dipoles whose order of the first incomplete shell is equal to K. It 
is clear that 
NdnKn 12)1()( γ≤≤ , (62)
where γ12 is surface to volume ratio of the scatterer. Finally we obtain 
( )[ ]dcddccNd 46245431 ln +−≤h . (63)
The last term in Eq. (63) is mostly determined by dipoles that lie on the surface (or few 
dipoles deep) because it comes from the K-4 term in Eq. (61) (which rapidly decreases when 
moving from surface). We define surface errors as those associated with the linear term in 
Eq. (63). Our numerical simulation (see Section 0) show that this term is small compared to 
other terms for “typical” values of d, however it is always significant for small enough values 
of d. 
From Eq. (18) we directly obtain 
δ ddd hAE
≤ . (64)
We assume that a bounded solution of Eq. (7) uniquely exists for any , moreover we 
assume that if 
inc~ H∈E
inc =E  then 131
γ≤E . These assumptions are equivalent to the fact that 
1~ −A  exists and is finite (the operator 1
~ −A  is bounded). Because dA  is a discretization of A
one would expect that 
( ) 131
lim γ== −
. (65)
Although Eq. (65) seems intuitively correct, its rigorous prove, even if feasible, lies outside 
the scope of this paper. For an intuitive understanding one may consult the paper by Rahola,25 
where he studied the spectrum of the discretized operator (for scattering by a sphere) and 
showed that it does converge to the spectrum of the integral operator with decreasing d. It 
should however be noted, that convergence of the spectrum only implies the convergence of 
the spectral (L2) norm of the operator and not necessarily the convergence of the L1-norm. 
Therefore Eq. (65) should be considered as an assumption. It implies that there exists a d0 
such that 
for 0dd <  ( ) 47
A , (66)
where c47 is an arbitrary constant larger then γ13 (although d0 depends on its choice). For 
example 1347 2γ=c  should lead to a rather large d0 (a rigorous estimate of d0 does not seem 
feasible). Therefore 
δ dE  satisfies the same constrain as 
dh  (Eq. (63)) but with constants 
c48-c50. 
Next we estimate the errors in the measured quantities and start with the discretization 
error (first part in Eq. (28)). Examining Eqs. (20) and (24) one can see that Eq. (32) may be 
directly applied leading to 
( ) ( ) 252551,0~~ dcdc
dd ≤≤− ∑EE φφ . (67)
The second part in Eq. (28) is estimated as 
( ) ( ) ( ) dcddccdcdc d
55541
,0 lnδδ +−≤≤≤− ∑ EEEE φφ , (68)
where we used Eq. (29). The estimation of the error for Cabs additionally uses the fact 
i c EE δδ 57
c Eδmax
57 . 
By combining Eqs. (67) and (68) we obtain the final result of this section: 
( ) dcddccd 5625558 lnδ +−≤φ . (69)
It is important to remember that the derivation was performed for constant x, m, shape, and 
incident field. There are 13 basic constants (γ1-γ13). γ1 (Eq. (29)) characterizes the total 
volume of the scatterer, hence it depends only on x. γ7-γ11 (Eq. (30) for χ(r)) can be easily 
obtained given the function χ(r), moreover it is completely trivial in the common case of 
homogenous scatterers. γ12 (surface to volume ratio, Eq. (62)) depends on the shape of the 
scatterer and is inversely proportional to x. It is not feasible (except for certain simple shapes) 
to obtain the values of constants γ2-γ6 (Eq. (30)), since it requires an exact solution for the 
internal fields. These constants definitely depend upon all the parameters of the scattering 
problem. Moreover, these dependencies can be rapidly varying, especially near the resonance 
regions. The same is true for γ13 (L1-norm of the inverse of the integral operator, Eq. (65)). 
Finally, there is the important constant d0 that also depends on all the parameters, however 
one may expect it to be large enough (e.g. ) for most of the problems – then its 
variation can be neglected. 
20 ≥d
Before proceeding we introduce the discretization parameter kdmy = . We employ the 
commonly used formula as proposed by Draine,8 however the exact dependence on m is not 
important because all the conclusions are still valid for constant m. Replacing d by y does not 
significantly change the dependence of the constants in Eq. (69) since they all already depend 
on m through the basic constants γ2-γ11, γ13. This leads to 
( ) ycyyccy 6126059 lnδ +−≤φ . (70)
It is not feasible to make any rigorous conclusions about the variation of the constants in 
Eq. (70) with varying parameters because all these constants depend on γ2-γ6, γ13 that in turn 
depend in a complex way upon the parameters of the scattering problem. However we can 
make one conclusion about the general trend of this dependency. 
Following the derivation of the Eq. (70) one can observe that c61 is proportional to γ12, 
while c59 and c60 do not directly depend on it (at least part of the contributions to them are 
independent of γ12). Therefore the general trend will be a decrease of the ratio 5961 cc  with 
increasing x (when all other parameters are fixed). This is a mathematical justification of the 
intuitively evident fact that surface errors are less significant for larger particles. 
In the analysis of the results of the numerical simulations (Section 0) we will neglect the 
variation of the logarithm. Eq. (70) then states that error is bounded by a quadratic function of 
y (for ). However, keep in mind that our derivation does not lead to an optimal error 
estimation, i.e. it overestimates the error and can be improved. For example, the constants γ
0dd ≤
γ6 are usually largest inside a small volume fraction of the scatterer (near the surface or some 
internal resonance regions), while in the rest of the scatterer the internal electric field and its 
derivatives are bounded by significantly smaller constants. However the order of the error is 
estimated correctly, as we will see in the numerical simulations. 
It is important to note that Eq. (70) does not imply that δφ y (which is a signed value) 
actually depends on y as a quadratic function, but we will see later that it is the case for small 
enough y (Section 0, see detailed discussion in Paper 2). Moreover, the coefficients of linear 
and quadratic terms for δφ y may have different signs, which may lead to zero error for non-
zero y (however, this y, if it exists, is unfortunately different for each measured quantity). 
E.Shape errors 
In this section we extend the error analysis as presented in Section 2.D to shapes that cannot 
be described exactly by a set of cubical subvolumes. We perform the discretization the same 
way as in Section 2.B but some of the Vi are not cubical (for Vi ∂∈ , which denotes that 
dipole i lies on the boundary of the volume V). We set ri to be still in the center of the cube 
(circumscribing Vi) not to break the regularity of the lattice. The standard PP prescription uses 
equal volumes ( ) in Eqs. 3dVi = (10), (14), (25), and (26), i.e. the discretization changes the 
shape of the particle a little bit. We will estimate the errors introduced by these boundary 
dipoles. These errors should then be added to those obtained in Section 2.D. We start by 
estimating 
dh . First we consider  for dih Vi ∂∉  
−′′′=
PGrPrrGh )0(33 )(),(d , (71)
which is just a reduction of Eq. (34). For Vi ∂∈   is the same plus the error in the self-term dih
iiiiii
i VVdr
EΙrLrMPGrPrrGh χ
−′′′= ∑ ∫
),(),()(),(d )0(33 . (72)
Let us define 
iij dr
PGrPrrGh )0(33sh )(),(d −′′′= ∫ , (73)
iiiiiiii VV EΙrLrMh χ
⎛ −∂−=
),(),(sh . (74)
We estimate each of the terms in Eq. (73) separately (since there is actually no significant 
cancellation, and the error is of the same order of magnitude as the values themselves) using 
Eq. (30) for P(r) and )(rG  and Eq. (31). This leads to 
h  (75)
To estimate  we assume that the surface of the scatterer is a plane on the scale of the size 
of the dipole. A finite radius of curvature only changes the constants in the following 
expressions. We will prove that 
sh cii ≤h , (76)
therefore we do not need to consider the third term in Eq. (74) (coming from the unity tensor) 
at all, since it is bounded by a constant. 
( ) ( )∫∫ −′′′+′′−′′=
iiii rrV )()(),(d)(),(),(d),(
s3s3 rPrPrrGrPrrGrrGrM . (77)
−′ ic rrThe function in the first integral is always bounded by 
. If  the same is true 
for the second integral and hence 
ii V∈r
dcV ii 66),( ≤rM . (78)
If ii V∉r  we introduce an auxilia  r ′′ry point  that is symmetric to ri over the particle surface 
and apply the identity 
( ) ( ))()()()()()( ii rPrPrPrPrPrP −′′+′′−′=−′  (79)
to the second integral in Eq. (77). Using Taylor expansion of P near  and the fact that  r ′′
irrrr −′≤′′−′  for iV∈′r  one can show that 
∫ ′′+≤ii cdcV ),( s36867r
ir ),(d rrGM , (80)
where the remaining integral can be proven to be equal to ),( iiV rL ∂− . The last prove left 
(see Eqs. (74) and (80)) is to demonstrate that ),( iiV rL ∂  is bounded by a constant. The only 
potential problem may come from the subsurface of iV∂  th  the particle surface 
(because it may be close to r
at is part of
ce is i). This subsurfa med planar. We will calculate the 
integral in Eq. 
(6) over the infinite plane rρrr +=−′ i uch that 0 s =⋅rρ . Then ρρn ±=′  
( ) 223
2d),plane.inf(
mm == ∫irL
+∂V r
r , (81)
which is bounded. The rest of the integral (over the part of the cube surface) is bounded by a 
constant, which is a manifestation of a more general fact that (by its definition) ), ir   ( iVL ∂
does not depend on the size but only on the shape of the volume. Finally we have 
69),( cV ii ≤∂ rL , 
which together with Eqs. (74), (78), and (80) prove Eq. (76). 
Using Eqs. (75) and (76) we obtain 
)ln()( 70
cllnc
+≤+≤ ∑ ∑∑∑ ∑ −hhh 7271
dccNd
, (83)
where we have changed the order of the summation in the double sum and split the 
summation over cubical shells for 
⎝∂∈∂∈∂∈
maxKl ≤  and . Then we have grouped everything 
otal esti
maxKl >
into one sum over boundary dipoles. Eqs. (41) and (62) were used in the last inequality. 
Combining Eqs. (63) and (83) one can obtain the t mate of the 
dh  for any scatterer: 
( ) ( )[ ]ddccddccNd lnln 7273245431 −+−≤h . (84)
Using Eq. (66) we immediately obtain the sam
δ de estimate for E . 
The derivation of the errors in the measured quantities is slightly modified compared to 
 (68) e changed to Section 2.D, by the presence of the shape errors. Eqs. (67) and  ar
( ) ( )~~ dcdcdcdc
,0 +≤+≤− ∑∑
EE φφ , (85)
( ) )( ( ) ( )ddccddccdddd lnln 2,0 −+−≤−Eφ . dcd δ 7655541
53≤ EEφ (86)
The second term in Eq. (85) comes from surface dipoles for which errors are the same order 
as the values themselves. Finally the generalization of Eq. (70) is 
( ) ( )yyccyyccy lnlnδ 797826059 −+−≤φ . (87)
The shape errors “reinforce” the surface errors (the linear term of discretization error), and 
although they both generally decrease with increasing size parameter x one may expect the 
linear term in Eq. (87) to be significant up to higher values of y than in Eq. (70). 
All the derivations in this section can in principle be extended to interfaces inside the 
particle, i.e. when a surface, which cannot be described exactly as a surface of a set of cubes, 
separates two regions where χ(r) varies smoothly. Two parts of the cubical dipole on the 
interface should be considered separately the same way as it was done above. This will 
however not change the main conclusion of this section – Eq.(87) – but only the constants. 
F.Different DDA formulations 
In this section we discuss how different DDA formulations modify the error estimates derived 
in Sections 2.D and 2.E. 
Most of the improvements of PP proposed in the literature are concerned with the self-
term – . They are the Radiative Reaction correction (RR),),( iiV rM
8 the Digitized Green’s 
Function (DGF),23 the formulation by Lakhtakia (LAK),26,27 the a1-term method,28,29 the 
Lattice Dispersion Relation (LDR),30 the formulation by Peltoniemi (PEL),31 and the 
Corrected LDR (CLDR).32 All of them provide an expression for  that is of order d),( iiV rM
(except for RR that is of order d3). For instance LDR is equivalent to 
( )[ ] iii ddSmbmbbV PrM 3223221 i)32(),( +++= , (88)
(remember that we assumed ) where b1=k 1, b2, b3 are numerical constants and S is a 
constant that depends only on the propagation and polarization vectors of the incident field. 
However, none of these formulations can exactly evaluate the integral in Eq. (39), because the 
variation of the electric field is not known beforehand (PEL solves this problem, but only for 
a spherical dipole). Therefore they (hopefully) decrease the constant in Eq. (40), thus 
decreasing the overall error in the measured quantities. However, these formulations are not 
expected to change the order of the error from d2 to some higher order. 
We do not analyze the improvements by Rahmani, Chaumet, and Bryant (RCB)33,34 and 
Surface Corrected LDR (SCLDR),17 as they are limited to certain particle shapes. 
There exist two improvements of the interaction term in PP: Filtered Coupled Dipoles 
(FCD)12 and Integration of Green’s Tensor (IT).35 A rigorous analysis of FCD errors is 
beyond the scope of this paper, but it seems that FCD is not designed to reduce the linear term 
in Eq. (63) that comes from the incomplete (non-symmetric) shells. This is because FCD 
employs sampling theory to improve the accuracy of the overall discretization for regular 
cubical grids. FCD does not improve the accuracy of a single ijG  calculation (approximation 
of an integral over one subvolume). 
IT, which numerically evaluates the integral in Eq. (12), has a more pronounced effect 
on the error estimate. Consider dipole j from l-th shell (incomplete) of dipole i, then 
.),(d),(maxd
)(),(d)(),(d
≤′′′+′∂′′≤
−′′′=−′′′
dlcrrcrrc
jijij
rrGrrG
PrPrrGPGrPrrG
 (89)
Here we have used Eq. (36) and Taylor expansion of Green’s tensor up to the first order. 
Eq. (89) states that the second term in Eq. (58) is completely eliminated and so is the linear 
term in Eqs. (69) and (70) (surface errors). Therefore convergence of DDA with IT for 
cubically shaped scatterers is expected to be purely quadratic (neglecting the logarithm). 
However, for non cubically shaped scatterers the linear term reappears, due to the shape 
errors. Both IT and FCD also modify the self-term, however the effect is basically the same as 
for the other formulations. 
Several papers aimed to reduce shape errors.10,11,36 The first one – Generalized Semi-
Analytical (GSA) method10 – modifies the whole DDA scheme, while the other two propose 
averaging of the susceptibility over the boundary dipoles. We will analyze here Weighted 
Discretization (WD) by Piller11 as probably the most advanced method to reduce shape errors 
available today. 
WD modifies the susceptibility and self-term of the boundary subvolume. We slightly 
modify the definition of the boundary subvolume used in Sections 2.B and 2.E to 
automatically take into account interfaces inside the scatterer. We define Vi to be always 
cubical, but with a possible interface inside. The particle surface, crossing the subvolume Vi, 
is assumed planar and divides the subvolume into two parts: the principal volume  
(containing the center) and the secondary volume  with susceptibilities ,  and 
electric fields ,  respectively. The electric fields are considered constant inside 
each part and related to each other via the boundary condition tensor 
iV ii χχ ≡
ii EE ≡
iT : 
iii ETE =
s . (90)
In WD the susceptibility of the boundary subvolume is replaced by an effective one, defined 
( ) 3ssppe dVV iiiiii TI χχχ += , (91)
which gives the correct total polarization of the cubical dipole. The effective self-term is 
directly evaluated starting from Eq. (4), considering χ and E constant inside each part, 
( ) ( ) iii
rrV ETrrGrrGrrGrrGrM
′−′′+′−′′= ∫∫ ss3ps3
),(),(d),(),(d),( χχ . (92)
Piller evaluated the integrals in Eq. (92) numerically.11 
To take a smooth variation of the electric field and susceptibility into account we define 
 ( r  is defined in Section )(s r ′′= χχ i ′′ 2.E) and iT  is calculated at the surface between ri and 
.  and r ′′ ii PP ≡
iiiiii ETEP
ssss χχ == . Then 
c rrPrP
−≤−′′
min)( 83
s , (93)
where we have assumed that Eq. (30) for χ(r) and E(r) is also valid in . siV
We start estimating errors of WD with  (cf. Eq. shijh (73)) 
( ) ( )∫∫ −′′′+−′′′=
s)0(3p)0(3sh )(),(d)(),(d
jijiij rr PGrPrrGPGrPrrGh , (94)
Using Taylor expansions of  near r)(rP ′ i and r ′′  in  and  correspondingly and Eq. piV
iV (93) 
one may obtain that the main contribution comes from the derivative of Green’s tensor, 
leading to (cf. Eq. (75)) 
h  (95)
iih  is the following (cf. Eq. (74)) 
( ) ( )
( ) ( ) ( )
.),(),(
),(d)(),(d)(),(d
),(),(),(d),(),(d
),(),(
pss3s3p3
ess3ps3
iiiiiii
iiiiiii
PrLErL
PPrrGPrPrrGPrPrrG
ErLPrrGrrGPrrGrrG
PrLrMh
−′′+−′′′+−′′′=
∂−′−′′+′−′′−
 (96)
The first two integrals can be easily shown to be dс85≤  (cf. Eq. (77)) and third one is 
transformed to L  the same way as in Eq. (80), thus 
iiiiiiiiiiii VVVdc ErLPrLPrLh
esspp
sh ),(),(),( χ∂−∂+∂+≤ , (97)
where the second term comes from the fact that averaged PL  is not the same as L  times 
averaged P. This error depends on the geometry of the interface inside Vi and generally is of 
order unity. For example, if the plane interface is described as ε+= izz , taking limit 0→ε  
gives the error ( )zii sp2 PP −π  (using Eq. (81)). Therefore WD does not principally improve 
the error estimate of  given by Eq. shiih (76), although it may significantly decrease the 
constant. On the other hand, since ),( p iiV rL ∂  and ),(
iiV rL ∂  can be (analytically) evaluated 
for a cube intersected by a plane, WD can be further improved to reduce the error in  to 
linear in d, which is a subject of future research. 
Proceeding analogous to the derivation of Eq. (83) one can obtain 
Ndcdccllnc
898887
)( ≤⎟⎟
++≤ ∑ ∑
h . (98)
It can be shown that for the scattering amplitude (Eq. (25)) the error estimate given by 
Eq. (85) can be improved, since WD correctly evaluates the zeroth order of value for the 
boundary dipoles, leading to 
( ) ( ) 291490551,0~~ dcdcdc
dd ≤+≤− ∑∑
EE φφ . (99)
In his original paper11 Piller did not specify the expression that should be used for Cabs. Direct 
application of the susceptibility provided by WD into Eq. (26) does not reduce the order of 
error when compared with the exact Eq. (24) (except when ), since they are not linear 
functions of the electric field. However, if we consider separately  and  (which is 
equivalent to replacing 
0s =iχ
( )( )∗⋅ iiiiV EEeIm χ  by 
2ss2pp )Im()Im( iiiiiii VV ETE χχ + ) the same 
estimate as in Eq. (99) can be derived for Cabs. 
Using Eqs. (98), (99), and the first part of Eq. (86) one can derive the final error 
estimate for WD: 
( ) ycyyccy 9429392 lnδ +−≤φ , (100)
where the constant before the linear term, as compared to Eq. (87), does not contain a 
logarithm and is expected to be significantly smaller, because several factors contributing to it 
are eliminated in WD. Although WD has a potential for improving, it does not seem feasible 
to completely eliminate the linear term in the shape error. The accuracy of evaluation of the 
interaction term over the boundary dipole (cf. Eq. (94)) can be improved by integration of 
Green’s tensor over  and  separately but that would ruin the block-Toeplitz structure of 
the interaction matrix and hinder the FFT-based algorithm for the solution of linear 
0 30 60 90 120 150 180
Scattering angle θ, deg
 cube kD=8
 discretized sphere kD=10
 sphere kD=3 (x10)
 sphere kD=10
 sphere kD=30
Fig. 2. S11(θ ) for all 5 test cases in logarithmic scale. The result for the kD = 3 sphere is multiplied by 
10 for convenience. 
equations.5 Since there is no comparable alternative to FFT nowadays, this method seems 
inapplicable. 
Minor modifications of the expression for Cabs are possible. Draine8 proposed a 
modification of Eq. (26) that was widely used afterwards and which was further modified by 
Chaumet et al.35 However, for many cases these expressions are equivalent and, even when 
they are not, the difference is of order d3, which is neglected in our error analysis. 
3. Numerical simulations 
A.Discrete Dipole Approximation 
The basics of the DDA method were summarized by Draine and Flatau.2 In this paper we use 
the LDR prescription for dipole polarizability,30 which is most widely used nowadays, e.g. in 
the publicly available code DDSCAT 6.1.4 We also employ dipole size correction8 for non-
cubically shaped scatterers to ensure that the cubical approximation of the scatterer has the 
correct volume; this is believed to diminish shape errors, especially for small scatterers.2 We 
use a standard discretization scheme as described in Section 2.E, without any improvements 
for boundary dipoles. It is important to note that all the conclusions are valid for any DDA 
implementation, but with a few changes for specific improvements as discussed in 
Section 2.F. 
Our code – Amsterdam DDA (ADDA) – is capable of running on a cluster of computers 
(parallelizing a single DDA computation), which allows us to use practically unlimited 
number of dipoles, since we are not limited by the memory of a single computer.37,38 We used 
a relative error of residual  as a stopping criterion for the iterative solution of the DDA 
linear system. Tests suggest that the relative error of the measured quantities due to the 
iterative solver is then  (data not shown) and hence can be neglected (total relative 
errors in our simulations are  – see Section 
810−<
710−<
56 1010 −− ÷> 0). More details about our code can 
be found in Paper 2. All DDA simulations were carried out on the Dutch national compute 
cluster LISA.39 
Table 1. Exact values of Qext for the 5 test cases. 
Particle Qext
kD = 8 cube 4.490 
discretized kD = 10 sphere 3.916 
kD = 3 sphere 0.753 
kD = 10 sphere 3.928 
kD = 30 sphere 1.985 
0.1 1
 slope = 0.77
 y = kd·m
 maximum
 θ = 0°
 θ = 45°
 θ = 90°
 θ = 135°
 θ = 180°
 slope = 0.95
Fig. 3. Relative errors of S11 at different angles θ and maximum over all θ versus y for (a) the kD = 8 
cube, (b) the cubical discretization of kD = 10 sphere. A log-log scale is used. A linear fit of maximum 
over θ errors is shown. (m = 1.5). 
B.Results 
We study five test cases: one cube with 8=kD , three spheres with , and a 
particle obtained by a cubical discretization of the 
30,10,3=kD
10=kD  sphere using 16 dipoles per D 
(total 2176 dipoles, x equal to that of a sphere; see detailed description in Paper 2). By D we 
denote the diameter of a sphere or the edge size of a cube. All scatterers are homogenous with 
. Although DDA errors significantly depend on m (see e.g. 5.1=m 14), we limit ourselves to 
one single value and study effects of size and shape of the scatterer. 
The maximum number of dipoles per D (nD) was 256. The values of nD that we used are 
of the form  (p is an integer), except for the discretized sphere, where all np2}7,6,5,4{ ⋅ D are 
0.1 1
0.1 1
0.01 0.1
y = kd·m
 slope 2.29 (c)
 slope 0.91
 maximum
 θ = 0°
 θ = 45°
 θ = 90°
 θ = 135°
 θ = 180°
 slope 1.05
Fig. 4. Same as Fig. 3 but for (a) kD = 3, (b) kD = 10, and (c) kD = 30 spheres. 
multiples of 16 (this is required to exactly describe the shape of the particle composed from a 
number of cubes). The minimum values for nD were 8 for the 3=kD  sphere, 16 for the cube, 
the  sphere, and the discretized sphere, and 40 for the 10=kD 30=kD  sphere. 
All the computations use a direction of incidence parallel to one of the principal axes of 
the cubical dipoles. The scattering plane is parallel to one of the face of the cubical dipoles. In 
this paper we show results only for the extinction efficiency Qext (for incident light polarized 
parallel to one of the principal axes of the cubical dipoles) and phase function S11(θ ) as the 
0.01 0.1 1
 slope = 0.89
 cube kD=8
 discretized sphere kD=10
 sphere kD=3
 sphere kD=10
 sphere kD=30
y=kd·m  
Fig. 5. Relative errors of Qext versus y for all 5 test cases. A log-log scale is used. A linear fit through 5 
finest discretizations of kD = 3 sphere is shown. 
most commonly used in applications. However, the theory applies to any measured quantity. 
For instance, we have also confirmed it for other Mueller matrix elements (data not shown). 
Exact results of S11(θ ) for all 5 test cases are shown in Fig. 2. For spheres this is the 
result of Mie theory (the relative accuracy of the code we used24 is at least ) and for the 
cube and discretized sphere an extrapolation over the 5 finest discretizations (the 
extrapolation technique is presented in Paper 2, together with all details of obtaining these 
results, including their estimated errors). We use such ‘exact’ results because analytical theory 
is unavailable for these shapes and because errors of the best discretization are larger than that 
of the extrapolation. Their use as references for computing real errors (difference between the 
computed and the exact value) of single DDA calculations is justified because all these real 
errors are significantly larger than the errors of the references themselves (see Paper 2; in 
general, real errors obtained this way have an uncertainty of reference error). Exact values of 
610−<
ext for all test cases are presented in Table 1. 
In the following we show the results of DDA convergence. Fig. 3 and Fig. 4 present 
relative errors (absolute values) of S11 at different angles θ and maximum error over all θ 
versus y in log-log scale. In many cases the maximum errors are reached at exact 
backscattering direction, then these two sets of points overlap. Deep minima that happen at 
intermediate values of y for some values of θ (and also sometimes for Qext – Fig. 5) are due to 
the fact that the differences between simulated and reference values change sign near these 
values of y (see Paper 2 for detailed description of this behavior). The solid lines are linear fits 
to all or some points of maximum error. The slopes of these lines are depicted in the figures. 
Fig. 5 shows relative errors of Qext for all 5 studied cases in log-log scale. A linear fit through 
the 5 finest discretizations of the 3=kD  sphere is shown. More results of these numerical 
simulations are presented in Paper 2. 
4. Discussion 
Convergence of DDA for cubically shaped particles (Fig. 3) shows the following trends. All 
curves have linear and quadratic parts (the non-monotonic behavior of errors for some θ are 
also a manifestation of the fact that signed difference can be approximated by a sum of linear 
and quadratic terms that have different signs). The transition between these two regimes 
occurs at different y (which indicates the relative importance of linear and quadratic 
coefficients). While for maximum errors that are close to those of the backscattering direction 
the linear term is significant for larger y, it is much smaller and not significant in the whole 
range of y studied for side scattering ( °= 90θ ). R
for l
esults of DDA convergence for spheres (Fig. 
4) show a different behavior for different sizes. Errors for the small ( 3=kD ) sphere converge 
purely linear (except for small deviation of errors of )90(11 °S   values of y). Similar 
results are obtained for the 10=kD  sphere, but with significant oscillations superimposed 
upon the general trend. Convergence for the large ( 30
=kD ) sphere is quadratic or even faster 
in the range of y studied, also with significant oscillations. 
Comparing Fig. 3 and Fig. 4 (especially Fig. 3(b) and Fig. 4(b) showing results for 
almost the same particles) one can deduce the following differences in DDA convergence for 
cubically and non-cubically shaped scatterers. The linear term for cubically shaped scatterers 
is significantly smaller, resulting in smaller total errors, especially for small y. All these 
conclusions, together with the size dependence of the significance of the linear term in the 
total errors, are in perfect agreement with the theoretical predictions made in Sections 2.D and 
2.E. Errors for non-cubically shaped particles exhibit quasi-random oscillations that are not 
present for cubically shaped particles. This can be explained by the sharp variations of shape 
errors with changing y (discussed in details in Paper 2). Oscillations for the  sphere 
Fig. 4(a)) are very small (but still clearly present), which is due to the small size of the 
particle and hence featurelessness of its light scattering pattern – the surface structure is not 
that important and one may expect rather small shape errors. Results for Qext (Fig. 5) fully 
support the conclusions. Errors of Qext for the large sphere at small values of y are 
unexpectedly smaller than for smaller spheres. This feature requires further study before 
making any firm conclusions, however there is definitely no similar tendency for S11(θ ) (cf. 
Fig. 4). 
We have also studied a  porous cube that was obtained by dividing a cube into 
27 smaller cubes and then removing randomly 9 of them. All the conclusions are the same as 
those reported for the cube, but with slightly higher overall errors (data not shown). 
In this paper we have used a traditional DDA formulation2 for numerical simulations. 
However, as we showed in Section 2.F several modern improvements of DDA (namely IT 
and WD) should significantly change its convergence behavior. IT should completely 
eliminate the linear term for cubically shaped scatterers, which should improve the accuracy 
especially for small y. WD should significantly decrease shape and hence total errors for non-
cubically shaped particles, moreover it should significantly decrease the amplitude of quasi-
random error oscillations because it takes into account the location of the interface inside the 
boundary dipoles. Numerical testing of DDA convergence using IT and WD is a subject of a 
future study. 
5. Conclusion 
To the best of our knowledge, we conducted for the first time a rigorous theoretical 
convergence analysis of DDA. In the range of DDA applicability ( 2<kd ) errors are bounded 
by a sum of a linear and quadratic term in the discretization parameter y; the linear term is 
significantly smaller for cubically than for non-cubically shaped scatterers. Therefore for 
small y errors for cubically shaped particles are much smaller than for non-cubically shaped. 
The relative importance of the linear term decreases with increasing size, hence convergence 
of DDA for large enough scatterers is quadratic in the common range of y. All these 
conclusions were verified by extensive numerical simulations. 
Moreover, these simulations showed that errors are not only bounded by a quadratic 
function (as predicted in Section 2), but actually can be (with good accuracy) described by a 
quadratic function of y. This fact provides a basis for the extrapolation technique presented in 
Paper 2. 
Our theory predicts that modern DDA improvements (namely IT and WD) should 
significantly change the convergence of DDA, however numerical testing of these predictions 
is left for future research. 
Acknowledgements
We thank Gorden Videen and Michiel Min for valuable comments on earlier version of this 
manuscript and anonymous reviewer for helpful suggestions. Our research is supported by the 
NATO Science for Peace program through grant SfP 977976. 
References 
 1.  E. M. Purcell and C. R. Pennypacker, "Scattering and adsorption of light by nonspherical dielectric 
grains," Astrophys. J. 186, 705-714 (1973). 
 2.  B. T. Draine and P. J. Flatau, "Discrete-dipole approximation for scattering calculations," J. Opt. Soc. 
Am. A 11, 1491-1499 (1994). 
 3.  B. T. Draine, "The discrete dipole approximation for light scattering by irregular targets," in Light 
Scattering by Nonspherical Particles, Theory, Measurements, and Applications, M. I. Mishchenko, J. W. 
Hovenier, and L. D. Travis, eds. (Academic Press, New York, 2000), pp. 131-145. 
 4.  B. T. Draine and P. J. Flatau, "User guide for the discrete dipole approximation code DDSCAT 6.1," 
http://xxx.arxiv.org/abs/astro-ph/0409262 (2004). 
 5.  J. J. Goodman, B. T. Draine, and P. J. Flatau, "Application of fast-Fourier-transform techniques to the 
discrete-dipole approximation," Opt. Lett. 16, 1198-1200 (1991). 
 6.  G. C. Hsiao and R. E. Kleinman, "Mathematical foundations for error estimation in numerical solutions of 
integral equations in electromagnetics," IEEE Trans. Ant. Propag. 45, 316-328 (1997). 
 7.  K. F. Warnick and W. C. Chew, "Error analysis of the Moment Method," IEEE Ant. Prop. Mag. 46, 38-53 
(2004). 
 8.  B. T. Draine, "The discrete-dipole approximation and its application to interstellar graphite grains," 
Astrophys. J. 333, 848-872 (1988). 
 9.  J. I. Hage, J. M. Greenberg, and R. T. Wang, "Scattering from arbitrarily shaped particles - theory and 
experiment," Appl. Opt. 30, 1141-1152 (1991). 
 10.  F. Rouleau and P. G. Martin, "A new method to calculate the extinction properties of irregularly shaped 
particles," Astrophys. J. 414, 803-814 (1993). 
 11.  N. B. Piller, "Influence of the edge meshes on the accuracy of the coupled-dipole approximation," Opt. 
Lett. 22, 1674-1676 (1997). 
 12.  N. B. Piller and O. J. F. Martin, "Increasing the performance of the coupled-dipole approximation: A 
spectral approach," IEEE Trans. Ant. Propag. 46, 1126-1137 (1998). 
 13.  N. B. Piller, "Coupled-dipole approximation for high permittivity materials," Opt. Comm. 160, 10-14 
(1999). 
 14.  A. G. Hoekstra, J. Rahola, and P. M. A. Sloot, "Accuracy of internal fields in volume integral equation 
simulations of light scattering," Appl. Opt. 37, 8482-8497 (1998). 
 15.  S. D. Druger and B. V. Bronk, "Internal and scattered electric fields in the discrete dipole approximation," 
J. Opt. Soc. Am. B 16, 2239-2246 (1999). 
 16.  Y. L. Xu and B. A. S. Gustafson, "Comparison between multisphere light-scattering calculations: 
Rigorous solution and discrete-dipole approximation," Astrophys. J. 513, 894-909 (1999). 
 17.  M. J. Collinge and B. T. Draine, "Discrete-dipole approximation with polarizabilities that account for both 
finite wavelength and target geometry," J. Opt. Soc. Am. A 21, 2023-2028 (2004). 
 18.  M. A. Yurkin, V. P. Maltsev, and A. G. Hoekstra, "Convergence of the discrete dipole approximation. II. 
An extrapolation technique to increase the accuracy," J. Opt. Soc. Am. A 23, 2592-2601 (2006). 
 19.  A. Lakhtakia, "Strong and weak forms of the method of moments and the coupled dipole method for 
scattering of time-harmonic electromagnetic-fields," Int. J. Mod. Phys. C 3, 583-603 (1992). 
 20.  F. M. Kahnert, "Numerical methods in electromagnetic scattering theory," J. Quant. Spectrosc. Radiat. 
Transf. 79, 775-824 (2003). 
 21.  C. P. Davis and K. F. Warnick, "On the physical interpretation of the Sobolev norm in error estimation," 
Appl. Comp. ElectroMagn. Soc. J. 20, 144-150 (2005). 
http://xxx.arxiv.org/abs/astro-ph/0409262
 22.  A. D. Yanghjian, "Electric dyadic Green's function in the source region," IEEE Proc. 68, 248-263 (1980). 
 23.  G. H. Goedecke and S. G. O'Brien, "Scattering by irregular inhomogeneous particles via the digitized 
Green's function algorithm," Appl. Opt. 27, 2431-2438 (1988). 
 24.  C. F. Bohren and D. R. Huffman, Absorption and scattering of Light by Small Particles, (Wiley, New 
York, 1983). 
 25.  J. Rahola, "On the eigenvalues of the volume integral operator of electromagnetic scattering," SIAM J. 
Sci. Comp. 21, 1740-1754 (2000). 
 26.  A. Lakhtakia and G. W. Mulholland, "On 2 numerical techniques for light-scattering by dielectric 
agglomerated structures," J. Res. Nat. Inst. Stand. Technol. 98, 699-716 (1993). 
 27.  J. I. Hage and J. M. Greenberg, "A model for the optical-properties of porous grains," Astrophys. J. 361, 
251-259 (1990). 
 28.  C. E. Dungey and C. F. Bohren, "Light-scattering by nonspherical particles - a refinement to the coupled-
dipole method," J. Opt. Soc. Am. A 8, 81-87 (1991). 
 29.  H. Okamoto, "Light scattering by clusters: the a1-term method," Opt. Rev. 2, 407-412 (1995). 
 30.  B. T. Draine and J. J. Goodman, "Beyond clausius-mossotti - wave-propagation on a polarizable point 
lattice and the discrete dipole approximation," Astrophys. J. 405, 685-697 (1993). 
 31.  J. I. Peltoniemi, "Variational volume integral equation method for electromagnetic scattering by irregular 
grains," J. Quant. Spectrosc. Radiat. Transf. 55, 637-647 (1996). 
 32.  D. Gutkowicz-Krusin and B. T. Draine, "Propagation of electromagnetic waves on a rectangular lattice of 
polarizable points," http://xxx.arxiv.org/abs/astro-ph/0403082 (2004). 
 33.  A. Rahmani, P. C. Chaumet, and G. W. Bryant, "Coupled dipole method with an exact long-wavelength 
limit and improved accuracy at finite frequencies," Opt. Lett. 27, 2118-2120 (2002). 
 34.  A. Rahmani, P. C. Chaumet, and G. W. Bryant, "On the importance of local-field corrections for 
polarizable particles on a finite lattice: Application to the discrete dipole approximation," Astrophys. J. 
607, 873-878 (2004). 
 35.  P. C. Chaumet, A. Sentenac, and A. Rahmani, "Coupled dipole method for scatterers with large 
permittivity," Phys. Rev. E 70, 036606 (2004). 
 36.  K. F. Evans and G. L. Stephens, "Microwave radiative-transfer through clouds composed of realistically 
shaped ice crystals .1. Single scattering properties," J. Atmos. Sci. 52, 2041-2057 (1995). 
 37.  A. G. Hoekstra, M. D. Grimminck, and P. M. A. Sloot, "Large scale simulations of elastic light scattering 
by a fast discrete dipole approximation," Int. J. Mod. Phys. C 9, 87-102 (1998). 
 38.  M. A. Yurkin, K. A. Semyanov, P. A. Tarasov, A. V. Chernyshev, A. G. Hoekstra, and V. P. Maltsev, 
"Experimental and theoretical study of light scattering by individual mature red blood cells with scanning 
flow cytometry and discrete dipole approximation," Appl. Opt. 44, 5249-5256 (2005). 
 39.  "Description of the national compute cluster Lisa," http://www.sara.nl/userinfo/lisa/description/ (2005). 
http://xxx.arxiv.org/abs/astro-ph/0403082
http://xxx.arxiv.org/abs/astro-ph/0409262
  /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
