## Abstract

The fulfillment of the reciprocity by five publicly available scattering programs is investigated for a number of different particles. Reciprocity means that the source and the observation point of a given scattering configuration can be interchanged without changing the result. The programs under consideration are either implementations of T-matrix methods or of the discrete dipole approximation. Similarities and differences concerning their reciprocity behavior are discussed. In particular, it is investigated whether and under which conditions reciprocity tests can be used to evaluate the scattering results obtained by the different programs for the given particles.

© 2012 Optical Society of America

## 1. Introduction

Elastic scattering of electromagnetic waves on single particles represents a basic physical process of great practical importance in such diverse fields as atmospheric and ocean optics, astrophysics, biomedical optics, material sciences, and nano-optics. In the past, a variety of different methods has been developed to deal with the scattering problem of, in particular, nonspherical objects (for an overview see, e.g., [1–11]). In general, they differ in the approaches used, and, consequently, in their capabilities to compute the scattering behavior of various particle classes. Corresponding computer programs are sophisticated, tested, and partly publicly available (see, e.g., [12–14] for a database of numerous programs hosted by the University of Bremen). However, they may lead to slightly different numerical results for a given scattering problem. This can also be the case for various implementations of the same method or even for different versions of the same program. The differences may increase when approaching the limits of the algorithms. On the other hand there are cases where only one single method exists for treating special scatterer types so that no comparative calculations with alternative methods are possible for validating the results obtained. In all these cases it is up to the user to finally judge the accuracy and correctness of the findings. This is, e.g., important in remote sensing applications. Different scattering models can lead to different results in the data processing and finally to different conclusions.

To gain confidence in the results, scattering calculations at the parameter sets of benchmarks can be conducted and compared with each other. Comparisons with the outcome of other reliable programs for scatterers of simple shape and composition is also possible. Both can give some clues how to assess the results of the particles under consideration. Moreover, the tests of the phase and scattering matrices, proposed by Hovenier and van der Mee [15], can be performed. Some necessary conditions, such as that the single scattering albedo has to be equal to or less than 1, can be also checked. Additionally, tests of the energy conservation or the fulfillment of the boundary condition at the particle surface can be taken into account. But successful program runs including passing through some checks can not guarantee the correctness of the results. For instance, the fulfillment of the boundary condition at the scatterer surface and relative convergence of the corresponding differential scattering cross sections (DSCS) in the far field can be two different things. This has been recently illustrated by Rother and Wauer [16]. They conducted a case study on the accuracy behavior of three different T-matrix methods for scalar scattering problems. They found that the fulfillment of the boundary condition by a special T-matrix method within a given accuracy does not necessarily mean that relative convergence of the corresponding DSCS is obtained with the same accuracy. They provided numerical examples in which the convergence was worse. On the other hand, they also showed the contrary behavior for another T-matrix scheme. So the check of the boundary condition seems not to be an appropriate criterion to evaluate far field scattering results.

Additionally, Rother and Wauer [16] investigated the fulfillment of the reciprocity by DSCS. Reciprocity means that the source and the observation point of a given scattering configuration can be interchanged without changing the result (e.g., [17]). It is a consequence of the symmetry of the Green’s function underlying the T-matrix approach (see [10], for example). It always has to be fulfilled by any scattering program within a given numerical accuracy. A reciprocity check requires no detailed knowledge of the numerical method and its implementation. It can be readily done provided that differential scattering quantities are delivered. It is a simple and quick test for both users and program developers of scattering software. One of the conclusions in [16] is that the fulfillment of the reciprocity represents a highly sensitive criterion that can be successfully applied to obtain physically reliable results (see also [1, 2], e.g.). In particular, it could be shown in [16] that the reciprocity is as sensitive as the Barber–Hill criterion [18, 19] for the relative convergence of DSCS. That is, there is a relation between the quality of the scattering behavior in the far field and the fulfillment of the reciprocity. In preliminary investigations, we could also observe this for spheroids by means of the program
`mieschka` [19] as a special T-matrix method implementation. A similar relation has been seen for two touching spheres by using the programs
`scsmfo1b` [20] and
`mstm` [21]. In a recent study [22] of light scattering by fractal aggregates of highly absorbing material computations were performed with the discrete dipole approximation (DDA) and the superposition T-matrix method, using the
`DDSCAT` and
`scsmfo1b` implementations, respectively. The
`DDSCAT` computations were performed in conjunction with the GKDLDR polarizability model (see Section 3.2.2 for details). It was found that the accuracy of the
`DDSCAT` results estimated by the use of the reciprocity condition correlates well with the accuracy estimates based on comparison with the T-matrix results.

To study its practical use, we investigated the fulfillment of the reciprocity condition by existing publicly available scattering programs, considering different particles of relatively simple shapes. In doing so, we focussed on two classes of programs: implementations of T-matrix methods and of the DDA. The results of this study are presented here. In Section 2 the reciprocity conditions, scatterer geometries, and configurations are described. The different programs considered are outlined in Section 3. In Section 4 the numerical results are presented and discussed. Some conclusions are given in Section 5.

## 2. Description of the problem

Consider Fig. 1 for illustration. A plane, linearly polarized incident field

propagates along the positive*z*-axis of a Cartesian coordinate system (

*x*,

*y*,

*z*) and impinges on a particle of, in general, arbitrary size, shape, and composition. Here,

*k*

_{0}= 2

*π*/

*λ*denotes the wave number of the outer free space at a wavelength of

*λ*. The unit vector

**e**

_{0}characterizes the polarization state of the incident field. Let the

*x*–

*z*-plane be the scattering plane in which the differential scattering quantities are measured. Then,

**e**

_{0}=

**e**

*defines a vertically (v) polarized incident field ${\mathbf{E}}_{\text{inc}}^{\text{v}}$. For a horizontally (h) polarized field ${\mathbf{E}}_{\text{inc}}^{\text{h}}$, we have*

_{y}**e**

_{0}=

**e**

*(see, e.g., [18]).*

_{x}**e**

*and*

_{x}**e**

*are unit vectors in the*

_{y}*x*- and

*y*-direction, respectively. Note that time harmonic fields are considered throughout and that their time dependence is always omitted.

The presence of the particle gives rise to a scattered field **E**_{sca}. It is also decomposed into a v- and h-polarized component with respect to the scattering plane. We define
${\mathbf{E}}_{\text{sca}}^{\text{v}}={E}_{\text{sca}}^{\text{v}}{\mathbf{e}}_{\varphi}$ and
${\mathbf{E}}_{\text{sca}}^{\text{h}}={E}_{\text{sca}}^{\text{h}}{\mathbf{e}}_{\theta}$, where **e*** _{ϕ}* and

**e**

*are unit vectors in the*

_{θ}*ϕ*- and

*θ*-direction, respectively (see again [18]).

The relation between the polarized scattered field components (
${E}_{\text{sca}}^{\text{v}}$,
${E}_{\text{sca}}^{\text{h}}$) in the far field approximation and the polarized incident field components (
${E}_{\text{inc}}^{\text{v}}$,
${E}_{\text{inc}}^{\text{h}}$) is governed by the scattering amplitude matrix *ℱ* (e.g., [17, 23]).

Based on the scattering amplitudes *F _{αβ}*, (

*α*,

*β*) = (v, h), all scattering quantities of interest can be computed. In the present study, polarized DSCS (e.g., [18, 23])

In Eqs. (2) and (3), *ϕ* = 0° or *ϕ* = 180° due to the restriction to the scattering plane, while the scattering angle *θ* runs from 0° to 180° for each *ϕ*. In the further investigations, the following convention is used (see, e.g., [18]). If *ϕ* = 0° then *θ* ∈ [0°, 180°]. If *ϕ* = 180° then *θ* = (360° − −*θ*′) ∈ [180°, 360°] with the angle *θ*′ ∈ [0°, 180°] is chosen. In this way, above relations become independent of *ϕ* while *θ* runs from 0° to 360°.

Now, consider the scattered field of a given particle at an arbitrary but fixed angle *θ* in the scattering plane. In the reciprocal case, the incident field impinges on the particle under this angle *θ*, and the scattered field is observed in negative *z*-direction. At the same time, the polarization states of the fields interchange. The reciprocity condition states that the DSCS have to be equal in both cases. In order to be more general, let (**n**_{inc}, **n**_{sca}) be the incident and scattered field propagation directions and (*α*, *β*) the polarizations with respect to a reference plane spanned by (**n**_{inc}, **n**_{sca}). The reciprocal case is obtained by considering the incident field in the direction −**n**_{sca} and the scattered field in the direction −**n**_{inc}. Then the reciprocity condition for the polarized DSCS is given by (compare, e.g., with [23, 1, 2])

In most scattering programs, however, the particle is rotated rather than the field. Reciprocal configurations are often realized by rotating the scatterer against the fixed incident field and observing the scattered field at an appropriately rotated angle. This is common in T-matrix methods, for instance. In these methods, a body system (*x*_{b}, *y*_{b}, *z*_{b}) is usually introduced which is fixed to the particle and in which its T-matrix is determined. The rotation of the body system with respect to the (laboratory) system (*x*, *y*, *z*) is described by the Eulerian angles (*ϕ*_{p}, *θ*_{p}, *ψ*_{p}) (e.g., [23, 18]). *ϕ*_{p} and *ψ*_{p} represent rotations about different *z*-axes while *θ*_{p} is a rotation about an *y*-axis. One of the advantages of doing so consists in the possibility to utilize symmetries of the particle, if present, within an appropriately chosen body system [24]. The main advantage in the context of reciprocity is, however, the ability to readily compute the scattering behavior of a particle for any fixed orientation by simple rotations of its T-matrix, once the latter is known in the body system (see, e.g., [23, 1, 2]).

In the present study, the following special reciprocal configurations have been considered. *ϕ*_{p} = 0° and *ψ*_{p} = 180° have been chosen throughout. In this way we always remain in the scattering plane. The initial particle orientation is characterized by *θ*_{p} = 0°. In this orientation, the scattered field is mostly observed at an angle *θ* = 90°. The associated reciprocal configuration is given by (*θ*_{p} = 90°, *θ* = 270°) (see Fig. 1). In some special cases, the configurations (*θ*_{p} = 0°, *θ* = 60°) and (*θ*_{p} = 120°, *θ* = 300°) have been also taken into account. Furthermore, only vv- and hh-polarized DSCS have been investigated. Thus the reciprocity condition (4) reduces to

*θ*

_{1}≤ 180°, and

*θ*

_{1}≥ 180°,

*α*= (v, h), omitting the (

*ϕ*

_{p},

*ψ*

_{p})-dependence. The relative reciprocity error

*c*

_{1}= (

*θ*

_{p,1},

*θ*

_{1}) and

*c*

_{2}= (

*θ*

_{p,2},

*θ*

_{2}) represents a measure for its fulfillment.

Now, let us have a look at the scatterer geometries and parameters used in the study. They are summarized in Table 1. Here, *r _{v}*,

*k*

_{0}

*r*and

_{v}*m*are the radius of the volume equivalent sphere, the volume equivalent size parameter, and the complex refractive index, respectively.

*a/b*denotes the aspect ratio of the spheroids with the semi-axes

*a*and

*b*. The aspect ratio of the circular cylinders with a height

*h*and a cross-sectional radius

*r*is

*h*/2

*r*. Furthermore,

*n*and

*ε*are the order and the deformation parameter, respectively, of the Chebyshev particles, the shape of which is given by

To compute the radius *r*_{0} of the underlying sphere from the volume equivalent one, the corresponding expression in [25] has been used. Note that the size parameters have been chosen so that DDA are able to treat them on a desktop PC. The refractive index of (1.6, 0.0005) is that one of a “mean” mineral aerosol (see, e.g., [26]), while (1.3130, 5.889×10^{−10}) is the refractive index of ice at a wavelength of 500 nm [27]. It is obvious that the combination of *k*_{0}*r _{v}* = 15 and

*m*= (1.6, 0.0005) represents the crucial one for the scattering programs in the present study. It characterizes the largest particles, considered, with the highest optical contrast to the background. The deformation parameter

*ε*= 0.05 and orders

*n*= (5, 45) have been taken from [28]. Note furthermore that touching bisphere means a two-sphere cluster (bisphere) with touching components.

The configurations for the spheroids, circular cylinders, and Chebyshev particles are equal to those of the bisphere shown in Figs. 1(a) and 1(b). That is, a plane wave incidence along the rotational axis of these particles is assumed in the initial orientation characterized by *θ*_{p} = 0°. Then, only rotations by *θ*_{p} within the scattering plane are considered to generate the reciprocal configurations. The special case of the cube is illustrated in Figs. 1(c) and 1(d).

## 3. Methods of solution

In this section, the programs considered in the study are outlined. Short descriptions of the underlying methods, the scatterer geometries which can be treated, the convergence strategies, and the most important parameters which influence the accuracy of the results are provided. We focus on practical aspects of using these programs to make our conclusions immediately relevant for their users.

#### 3.1. T-matrix codes

### 3.1.1. Program
`mieschka`

Program
`mieschka` computes integral and differential scattering quantities for rotationally symmetric, homogeneous and isotropic particles in fixed and random orientations with respect to a plane, linearly polarized incident field. It can be applied to particles the size of which is comparable to or smaller than the wavelength of the incoming radiation. A detailed description of
`mieschka` is provided in [19, 10]. Its reliability has been proven in many comparisons with benchmarks and results of other programs as well as in numerous different applications. It is accessible via the German Aerospace Center’s Virtual Laboratory (VL; [19]) for registered users. Additionally, it is delivered in the book by Rother [10].
`mieschka` version 1.0.1 of VL has been used in the present investigations.

Program
`mieschka` is an implementation of a T-matrix method and written in FORTRAN 95. Formulated in spherical coordinates, it involves expansions of the incident, scattered and internal fields (**E**_{inc}, **E**_{sca}, **E**_{int}) in term of vector spherical wave functions **Ψ*** _{τ}_{nl}* and Rg

**Ψ**

*, respectively.*

_{τ}_{nl}Here, *k _{s}* denotes the wave number inside the scatterer. To compute the unknown scattered and internal field expansion coefficients (

*f*,

_{τ}_{nl}*p*), the continuity conditions of the tangential field components across the particle surface are applied [29]. Inserting the field expansions (9), (10), and (11) into the continuity conditions and applying a projection scheme to the resulting equations with certain weighting functions lead to the T-matrix equation

_{τ}_{nl}It relates the vector **f** of the scattered field expansion coefficients to that of the known incident field expansion coefficients **a**. The elements of the T-matrix represent integrals of combinations of vector spherical wave functions and weighting functions over the particle surface. Note that Waterman’s weighting functions [30] have been used in all investigations. If necessary, the internal field expansion coefficients can be calculated in a similar way. Once the scattered field is known, the scattering quantities of interest are calculated.

In any practical computation, the series expansions in (9), (10), and (11) have to be truncated at a finite *n*_{cut}, leading to a finite sized T-matrix. Moreover, the sum over *l* can also be truncated at a certain *l*_{cut} ≤ *n*_{cut}, depending on the particle shape. The parameters (*n*_{cut}, *l*_{cut}) are automatically determined in the standard mode of
`mieschka`, which was used in the computations. The procedure is similar to that by Barber and Hill [18] and has been described in detail in [19, 10]. It is however important to note that (*n*_{cut}, *l*_{cut}) are fixed in such a way that the hh- and vv-polarized DSCS are convergent within a relative error of 5% at 80% of the scattering angles. In doing so, *n*_{cut} is increased in steps of three, beginning from a defined starting value, while *l*_{cut} is increased in steps of one. The accuracy in computing the surface integrals of the T-matrix elements is also automatically controlled. Due to the rotational symmetry of the scatterers, only *θ*-integrals remain. They are numerically solved by a Gauss-Kronrod quadrature with an automatic step-size determination. The latter is realized in such a way that a relative integration error of 0.1% is achieved, resulting in a special number of integration points *n*_{int}.

### 3.1.2. Programs
`scsmfo1b` and
`mstm`

Mackowski, Fuller, and Mishchenko [20] have developed codes for the calculation of the scattering matrix and cross sections of neighboring, non–intersecting spheres. These codes, written in FORTRAN 77, are available via ftp since the late 1990s. In particular, the program
`scsmfo1b` is designed to calculate the entire 2-D scattering matrix and cross sections for large-scale sphere clusters in fixed orientations of the cluster with respect to a plane, linearly-polarized incident field. It is based on the work in [31, 32] and represents a superposition extension of the Mie theory to multiple spheres.

In
`scsmfo1b`, each of the individual spheres is described by an expansion in terms of outgoing vector spherical wave functions analogous to Eq. (9), centered about the origin of the sphere. The application of the continuity conditions at the surface of each sphere results in a system of interaction equations for the scattered field expansion coefficients. These equations are solved iteratively for a given cluster configuration. The scattered field of the entire cluster is described by a single vector spherical harmonic expansion (9), centered about the cluster origin, by translating the individual sphere expansions to the common cluster origin. Based on this field, above mentioned scattering quantities are computed. Consequently, the iteration strategy to obtain the solution and its numerical accuracy are governed by the parameters (*itermax*, *eps*, *meth*, *qeps*1, *qeps*2) which are fixed in the input file *scsmfo.inp*. Here, *itermax* denotes the maximum number of iterations in the solution method, *eps* the relative residual error tolerance of the solution, *meth* the solution method, *qeps*1 the error tolerance for determining the single-sphere harmonic order truncation, and *qeps*2 the error tolerance for determining the cluster expansion truncation limit. Note that the original parameter values, provided by the program authors, have not been changed in our computations for the bispheres. That is, they were (*itermax*, *eps*, *meth*, *qeps*1, *qeps*2) = (400, 1 × 10^{−6}, 0, 1 × 10^{−5}, 1 × 10^{−8}), throughout. Note furthermore that above theoretical description can be converted into a T-matrix formulation.

Recently, a revised and modernized version of above codes, called
`mstm`, has been released [21]. It is written in FORTRAN 90 and combines above codes for fixed and random orientation into a single package. It is able to treat systems of several thousand spheres. The source code and supporting documentation are available for free download at [33]. Note that the program version 2.1, released on 30 March 2011, has been used in the present investigations.

Several changes have been made during the revision. In particular, target T-matrices for each sphere and a cluster T-matrix, which treats the ensemble of spheres as a single particle, are derived. One of the main steps in a program run is the iterative solution of the system of interaction equations for the sphere-target T-matrices. So, the key parameters influencing the convergence behavior of the program
`mstm` are (*mie_epsilon*, *translation_epsilon*, *solution_epsilon*, *max_number_iterations*). Here, *mie_epsilon* characterizes the convergence criterion for determining the number of orders in the Mie expansions for each sphere, *translation_epsilon* the convergence criterion for estimating the maximum order of the cluster T-matrix, *solution_epsilon* the error criterion for the solution of the interaction equations, and *max_number_iterations* the maximum number of iterations used in the biconjugate gradient scheme for a particular right hand side. The following parameter values have been chosen when conducted the computations for the bispheres: (*mie_epsilon*, *translation_epsilon*, *solution_epsilon*, *max_number_iterations*) = (1 × 10^{−4}, 1 × 10^{−3}, 1 × 10^{−10}, 2000). Note that they are given in the
`mstm` 2.1 manual.

#### 3.2. DDA codes

The DDA is based on replacing a scatterer by a set of small volume elements (dipoles) which are characterized by their polarizabilities *α*. In case of an isotropic material, the latter – generally a tensor – is expressed as

*V*

_{d}is the volume of the dipole,

*χ*= (

*ε*− 1)/4

*π*and

*ε*are the susceptibility and the electric permittivity of the medium at the dipole center, and

**Ī**is the identity tensor.

*α*

^{CM}denotes the Clausius-Mossotti (CM) polarizability which is equivalent to assume

**M̄**= 0.

The term **M̄** is associated with the finiteness of the dipole which can be derived from different approximate assumptions. The most commonly used polarizability formulation is the lattice dispersion relation [34] (LDR; sometimes also abbreviated as LATTDR)

*m*is the refractive index at the dipole center,

*d*the dipole size (edge length of a cubical dipole), and

*S*a factor depending on the geometry of the (linearly polarized) incident wave. The latter is characterized by two unit vectors: the propagation direction

**a**and the polarization direction

**e**

^{0}of the incident electric field. The summation in the expression for

*S*is performed over three vector components. Later a minor flaw in the LDR derivation was found and corrected [35]. This corrected LDR (CLDR; also referred to as Gutkowicz-Krusin and Draine LDR (GKDLDR)) is independent on the incident polarization but leads to a diagonal polarizability tensor instead of scalar

*δ*is the Kronecker symbol.

_{μν}There exist a number of other polarizability formulations [36, 37], but the LDR and the CLDR are unique in their dependence on the incident direction. As an example of alternative formulation, we consider the radiative reaction correction (RRC, [38])

The core part of DDA consists in the solution of a system of linear equations. One of the equivalent forms of the latter is [36]

**P**

*are the unknown dipole polarizations, ${\mathbf{E}}_{i}^{\text{inc}}$ is the incident electric field at the dipole center,*

_{i}**Ḡ**

*is the interaction term, and the indices*

_{ij}*i*and

*j*enumerate the dipoles. For the interaction term, the Greens tensor for point dipoles is used:

**r**

*is the radius-vector of the dipole center,*

_{i}**R**=

**r**

*−*

_{j}**r**

*,*

_{i}*R*=

*|*

**R**

*|*, and

*R̂R̂*=

*R*. Although finite-size corrections for

_{μ}R_{ν}**Ḡ**

*are possible [36], we do not discuss them here. We want to refer to, e.g, [37] for further details.*

_{ij}### 3.2.1. Program
`ADDA`

`ADDA` is an open-source implementation of the DDA which is capable of running on a cluster of computers parallelizing a single DDA computation [36]. The solution of Eq. (19) is performed by an iterative method [36] until the relative residual (norm of the difference between the left-and right-hand-sides of Eq. (19) divided by the norm of
${\mathbf{E}}_{i}^{\text{inc}}$) is less than the specified threshold *ε*_{it}. The typical value for the latter is 10^{−5} which is much smaller than the typical total error of an
`ADDA` simulation, measured by the errors in scattering quantities. Hence, a particular value of *ε*_{it} is not relevant in many cases, but there is an important exception. While the finite size of dipoles is accounted for by a polarizability formulation, Eq. (19) is a direct implication of the Maxwell equations for a set of point dipoles with specified polarizabilities. Hence, an exact solution of Eq. (19) exactly satisfies the Maxwell equations (and all its implications), albeit not for the original scatterer but for a set of point dipoles. In this respect, finite residual (and hence the value of *ε*_{it}) is the only cause of inaccuracy for any tests based on satisfaction of Maxwell equations.

A standard DDA simulation is performed for a single value of *d*. Then the accuracy of the result obtained is largely unknown. There exist a number of benchmark studies, reviewed in [37], but their applicability for a specific application is often arguable. Therefore, the only feasible approach to estimate the DDA accuracy, when no reference result is available, is to perform DDA simulations with different *d*. The accuracy can then be estimated from a variation of the results with decreasing *d*. For example, simple differences between successive values of *d* can be used, similar to the convergence criterion of the T-matrix method (see Section 3.1.1). A more rigorous approach is based on the extrapolation of the result to zero value of *d*. A specific procedure is described in detail in [39]. Although it is, to a large extent, an empirical technique, its consistency (in particular, reliability of the error estimate) was studied on a number of cases in the original paper [39], and it was successfully applied in [40–42]. In the current contribution, we use the extrapolation technique for a few specific cases to obtain accurate
`ADDA` results together with estimate of its error (uncertainty).

`ADDA` 1.0 was used for almost all simulations. However, the ability to handle Cheby-shev particles was implemented specially for this manuscript. Hence,
`ADDA` 1.1b1 was taken for these shapes. All internal parameters of
`ADDA` were set to default values, unless noted otherwise. In particular, LDR polarizability was chosen, *ε*_{it} = 10^{−5}, and dipole size *d* = max(*λ*/(10|*m*|), *D _{x}*/16), where

*D*is the size of the particle along the

_{x}*x*

_{b}-axis in the body system.

`ADDA`uses double precision for all calculations.

### 3.2.2. Program
`DDSCAT`

`DDSCAT` is a publicly available, open-source DDA code written in FORTRAN 90 for computing electromagnetic scattering and absorption by objects with arbitrary geometries and, in general, inhomogeneous and anisotropic dielectric properties. The objects can be finite particles or 1D or 2D periodic structures. Details about
`DDSCAT` can be found in [43] and at [44] or [45]. The user guide is provided in [46]. Note that
`DDSCAT` version 7.1 has been used in our study.

The accuracy of
`DDSCAT` computations for finite targets is mainly controlled by the dipole spacing *d* and the error tolerance *TOL*. Depending on the chosen target geometry, the user specifies the size of the target and the number of dipoles along different target axes. For example, for a rectangular prism with side lengths *a*, *b*, and *c*, the user specifies *a/d*, *b/d*, and *c/d*, as well as the volume-equivalent radius *a*_{eff}, so that
$4\pi {a}_{\text{eff}}^{3}/3=abc$. For computing differential scattering properties, it is recommended to choose *d* such that |*m*|*k*_{0}*d* < 0.5, where *m* denotes here the maximum of the (in general varying) refractive index of the target. The error tolerance determines the termination of the conjugate gradient iteration. The iteration scheme continues until the linear equation system is solved with a fractional error *TOL*. Most sample input files in
`DDSCAT` version 7.1 use *TOL* = 10^{−5}, which is also the setting we used in our calculations. All calculations are performed in single precision. The output files give differential scattering properties up to the third digit after the decimal point, which is the overall accuracy with which we present
`DDSCAT` results in this paper.

There are presently 24 classes of geometries for finite targets and 11 classes of periodic targets implemented, for which the code automatically generates an array of dipoles with user-specified target dimensions and dipole spacing. In addition, the code can read in user-generated input files of dipole arrays, as well as files with spectral dielectric properties for different materials, thus extending the applicability of the code to practically arbitrary shapes and arbitrary inhomogeneous and anisotropic compositions.

In the input parameter file, the user must specify the prescription of the dipole polarizabilities. The two available options are LATTDR and GKDLDR (see Section 3.2 for details). The user manual [46] recommends the use of GKDLDR; but it is noted that LATTDR in many cases gives results of similar quality. In all
`DDSCAT` calculations presented in this paper, we used GKDLDR (CLDR).

## 4. Computational results and discussion

#### 4.1. T-matrix results

### 4.1.1. Program
`mieschka`

Table 2 shows the reciprocity errors *e*_{hh} and *e*_{vv}, defined by Eq. (7), for the reciprocal configurations (*θ*_{p} = 0°, *θ* = 90°) and (*θ*_{p} = 90°, *θ* = 270°) of different particles. In most cases, the errors are less than 5%, which is below the required accuracy for determining DSCS within
`mieschka` (see Section 3.1.1). The only exception is *e*_{vv} for cyl_2 of 15%. Fig. 2 shows the corresponding differential scattering behavior of this particle. It is seen that the DSCS values at *θ* = 90° for *θ*_{p} = 0° (solid line) and at *θ* = 270° for *θ*_{p} = 90° (dashed line) are located near deep down spikes. Relative convergence of DSCS is, however, not tried to be achieved in such regions according to the convergence strategy of
`mieschka`. So the poor fulfillment of the reciprocity condition (5) is not surprising in this particular case. But decreasing the relative error within the Barber–Hill criterion of 5% to 1% while keeping the integration error fixed has also an indirect effect on this. We have *e*_{hh} = 4.3% and *e*_{vv} = 15% with the accuracy of 5%. In the case of 1% accuracy, *e*_{hh} = 0.4% and *e*_{vv} = 1.8% are obtained. Note that a refinement of the numerical surface integration does not significantly influence these results. To circumvent such situations as for cyl_2, more appropriate reciprocal configurations should be taken into account.

It is furthermore seen that the reciprocity is, in general, better fulfilled for the spheroids and Chebychev particles than for the circular cylinders by
`mieschka`. The reciprocity error ranges between 3.2 × 10^{−4}% and 3.1% for the former particles whereas it is within 1.2 × 10^{−1}% – 4.3% for the cylinders (except *e*_{vv} for cyl_2, of course). The latter error range is, however, still below the required accuracy of 5%. The main difference between the former and the latter particles consists in the fact that finite circular cylinders exhibit 90°-edges, while spheroids and Chebychev particles have smooth, continuously differentiable boundary surfaces. Based on above investigations for cyl_2, one can conclude that the differences in the fulfillment of the reciprocity condition originate from the inappropriate field expansions (9) – (11) in terms of vector spherical wave functions especially around the cylindrical edges.

Now let us turn to the question whether the reciprocity fulfillment represents an appropriate criterion to evaluate the scattering results of
`mieschka`. The accuracy of a
`mieschka` run, i.e. the relative convergence of DSCS, is controlled by the relative error of the Barber–Hill criterion (see Section 3.1.1). Therefore, investigations concerning the relation between this error and the resulting reciprocity error have been conducted. Fig. 3(a) provides an example of these investigations for the scatterer sphd_p_2 in the two reciprocal configurations (*θ*_{p} = 0°, *θ* = 90°) and (*θ*_{p} = 90°, *θ* = 270°). First of all it can be observed that the reciprocity errors are always smaller than the required relative errors. It is furthermore seen that the reciprocity errors decrease and remain constant, respectively, with decreasing relative error. In the former case a decreased relative error leads to larger values of the convergence parameters (*n*_{cut}, *l*_{cut}, *n*_{int}). The resulting higher order approximation of the scattered field leads, in turn, to a better reciprocity fulfillment as long as numerical instabilities do not arise. On the other hand, a given set of (*n*_{cut}, *l*_{cut}, *n*_{int}) can meet different accuracy requirements of a
`mieschka` run (see Table 3). Equal sets of convergence parameters yield, of course, the same DSCS and, consequently, the same reciprocity errors. This represents one reason for the plateaus of constant *e _{αα}* values observed in Fig. 3(a). Equal DSCS within a given numerical accuracy result also if

*l*

_{cut}varies only slightly with varying relative error at oblique incidence while (

*n*

_{cut},

*n*

_{int}) do not change. In this case,

*e*becomes also constant for different relative errors (compare Fig. 3(a) and Table 3). So,

_{αα}*e*alters with varying relative error only if

_{αα}*n*

_{cut}changes at the same time. The relation between

*e*and

_{αα}*n*

_{cut}is illustrated in Fig. 3(b). Note that this behavior has also been found above for cyl_2. A decrease of the relative error from 5% to 1% led to an increase of

*n*

_{cut}from 36 to 42 while (

*l*

_{cut},

*n*

_{int}) remained constant. As a result, the reciprocity error decreased. So the relative error of the Barber–Hill criterion determines the relative convergence of DSCS and the reciprocity fulfillment in the same way via (

*n*

_{cut},

*l*

_{cut},

*n*

_{int}). The values of

*e*represent, therefore, a measure for the DSCS accuracy.

_{αα}### 4.1.2. Programs
`scsmfo1b` and
`mstm`

Table 4 summarizes the reciprocity errors *e*_{hh} and *e*_{vv} for the reciprocal configurations (*θ*_{p} = 0°, *θ* = 90°) and (*θ*_{p} = 90°, *θ* = 270°), obtained with
`scsmfo1b` and
`mstm`, respectively, for the bispheres. The fulfillment of the reciprocity condition by both programs is comparable with each other when using the convergence parameters given in Section 3.1.2. The orders of magnitude of their reciprocity errors are also comparable with those of the program
`mieschka` for, e.g., the spheroids (see Table 2).

Since both programs use, in particular, scattered field expansions the question arises whether there is any relation between a corresponding truncation parameter and the reciprocity error as found above for the program
`mieschka`. That is, does the reciprocity fulfillment tell us anything about the degree of the scattered field approximation and the resulting DSCS ? Fig. 4 provides an example of the corresponding investigations conducted with the newer program
`mstm`. It shows the relative reciprocity errors *e*_{vv} and *e*_{hh} for bi_sph_2 in the two reciprocal configurations (*θ*_{p} = 0°, *θ* = 90°) and (*θ*_{p} = 90°, *θ* = 270°) at different values of *translation_epsilon*. As noted in Section 3.1.2, this parameter determines the maximum order of the cluster T-matrix and, in this way, the number of the scattered field expansion terms. It has been increased in steps of a factor 10 in Fig. 4. The other parameters, influencing the convergence behavior of
`mstm`, have been fixed to the values listed in Section 3.1.2. In this way, the quality to treat each individual sphere and to solve the interaction equations remained unchanged. The curves of *e*_{vv} and *e*_{hh} versus *translation_epsilon* exhibit some variations. Nevertheless, the tendency of decreasing reciprocal errors with a decreased *translation_epsilon* is seen in the range between *translation_epsilon* = 1 to 1 × 10^{−7}. A further decrease of this parameter beyond 1 × 10^{−7} seems to result in constant reciprocal errors. Since smaller *translation_epsilon* values yield, in general, higher order approximations of the scattered field and of DSCS, the reciprocal error can be regarded as an indicator for the DSCS accuracy, provided that no numerical instabilities occur.

#### 4.2. DDA results

`ADDA` and
`DDSCAT` have not only a common theoretical basis, described in Section 3.2, and use partly the same polarizability models (LDR/LATTDR and CLDR/GKDLDR). Their numerical results are also equivalent. Test computations by
`ADDA` with CLDR and
`DDSCAT` with GKDLDR for sphd_p_1 in the two orientations *θ*_{p} = 0° and *θ*_{p} = 90° have provided the same DSCS within all digits. In doing so, the same dipole size have been used. This conclusion has been also confirmed in [47]. Therefore, each of both implementations can be regarded to be a part of a single comprehensive DDA program with various selectable polarizability models. Correspondingly, we will summarize the results obtained by
`ADDA` and
`DDSCAT` in a single section and refer only to the special program, used, as necessary.

Let us have a look at the
`ADDA`–LDR results first. Table 2 shows the reciprocity errors *e*_{hh} and *e*_{vv} for the configurations (*θ*_{p} = 0°, *θ* = 90°) and (*θ*_{p} = 90°, *θ* = 270°). It should be noted that results for Chebyshev particles by a DDA implementation are presented here, to the best of our knowledge, for the first time. It is seen that the reciprocity errors range between about *e*_{hh} = 6.3 × 10^{−6}% for sphd_p_4 and *e*_{vv} = 1.2 × 10^{−1}% for cyl_2. Only *e*_{vv} for the cubes is significantly larger being within about 2.1% – 160%. A comparison for the spheroids, cylinders, and Chebyshev particles shows that the reciprocity is in most cases significantly better fulfilled than by
`mieschka` with its accuracy of 5%. Only *e*_{hh} for sphd_p_2 and cheb_45_3, and *e*_{vv} for sphd_o_2 are slightly larger. The reciprocity errors for the bispheres are also smaller than those of
`scsmfo1b` and
`mstm`. The question is, why DDA with LDR yields, in most cases, a much better reciprocity fulfillment than above T-matrix programs. And what about *e*_{vv} for the cubes?

An answer for these questions has been hinted in Section 3.2.1 using a notion of a set of point dipoles. Here we expand it in a more rigorous manner. Consider the equation system (19) or its initial integral equation in [37], for instance. One of the main ingredients in these relations is the free space dyadic Green’s function **Ḡ** (**r*** _{i}*,

**r**

*). It obeys the symmetry relation*

_{j}**Ḡ**(

**r**

*,*

_{i}**r**

*) =*

_{j}**Ḡ**(

**r**

*,*

_{j}**r**

*). As was shown by Rother [10], this and the symmetry relation for the Green’s dyadic, belonging to the given boundary value problem, lead to a corresponding symmetry relation for the matrix elements of the dyadic interaction operator (or the T-matrix). This, in turn, results in the reciprocity property of the scatterer considered. So, above symmetry relation for*

_{i}**Ḡ**(

**r**

*,*

_{i}**r**

*) represents one condition for fulfilling reciprocity. As long as the other quantities in (19) such as the polarizabilities do not extinguish the symmetry effect of*

_{j}**Ḡ**(

**r**

*,*

_{i}**r**

*), i.e. their values are the same for reciprocal configurations, the solution of the equation system meets a priori the reciprocity condition. The numerical accuracy of its fulfillment is, then, solely determined by the solution scheme applied to the equation system.*

_{j}In all cases, except the cubes in the vv-polarization, we have just this situation. It is a direct consequence of the LDR polarizability model. Remember that the formula for the LDR polarizability contains the factor *S* (see Eq. (15)). This factor depends on the propagation direction **a** and polarization direction **e**^{0} of the incident electric field in the body system (see Eq. (16)). If either **a** or **e**^{0} is along a coordinate axis, then *S* = 0. As a consequence, the symmetry effect of **Ḡ** (**r*** _{i}*,

**r**

*) is not extinguished so that the reciprocity condition is fulfilled a priori. In the*

_{j}`ADDA`simulations, the cube is always oriented along the coordinate axes in the body system. Hence, the scattering plane has to be rotated along the z-axis to match the scattering configuration shown in Figs. 1(c) and 1(d). Therefore, we have

**a**= (−0.7, 0.7, 0) and

**e**

^{0}= (0.7, 0.7, 0), and hence

*S*= 0.5, for the cubes at a v-polarized field incidence through the edge in the configuration (

*θ*

_{p}= 90°,

*θ*= 270°). This configuration leads to a different polarizability of the dipoles so that reciprocity is not expected to be satisfied as accurately as in the other cases. Note that a similar situation can be found for the spheroids sphd_p_2 and sphd_o_2 in the configurations (

*θ*

_{p}= 0°,

*θ*= 60°) and (

*θ*

_{p}= 120°,

*θ*= 300°), for instance. Here,

*e*

_{vv}= 1.3 × 10

^{−3}% and 1.6×10

^{−2}%, respectively, whereas

*e*

_{hh}= 22% and 52%, respectively. Again, we have

*S*≠ 0 in the second configuration (

*θ*

_{p}= 120°,

*θ*= 300°), but now for the h-polarized incident field.

In order to illustrate this explanation, consider cube_1 and change the polarizability model from LDR to RRC. In doing so, all other parameters are fixed, i.e. we have 16 dipoles per dipole edge and *ε*_{it} = 10^{−5}. The RRC model is independent of the incident field, i.e. of both **a** and **e**^{0} (see Section 3.2). Therefore it is expected that the *e*_{vv} value decreases to the order of *e*_{hh} in Table 2. And indeed, *e*_{vv} drops down from about 3.6% to 1.1 × 10^{−3}% while *e*_{hh} = 1.2 × 10^{−3}% stays nearly on the same order as with LDR (see second column in Table 5 and compare with Table 2). A refinement of the discretization from 16 to 32 dipoles per dipole edge changes the values of the hh- and vv-polarized DSCS, but has only a minor effect on the reciprocity error (compare second and third column of Table 5). The latter stays on the order of about 10^{−3}% as for 16 dipoles. A decrease of the threshold of the iterative solver down to 10^{−10} has only a small effect on the DSCS values themselves, but reduces the reciprocity errors down to the order of 10^{−9}% (compare second and fourth column of Table 5). So the accuracy of the reciprocity fulfillment at RRC is mainly determined by the threshold *ε*_{it} of the iterative solver, as stated above. This is, however, not the main factor which determines the accuracy of a DDA simulation itself. It is first of all affected by the discretization, i.e., by the number of dipoles. The discretization, in turn, does not show a relation to the reciprocity accuracy within RRC.

That is, the reciprocity does not correlate with the DSCS quality if the former is a priori fulfilled either by a specially chosen polarizability model, which is independent of the incident field as RRC, or by the scattering configuration, as above for most particles at LDR. And small values of the reciprocity error do not necessarily correspond to accurate DSCS. Fig. 5 provides an example for this. It shows a comparison of the hh-polarized DSCS for the scatterer sphd_o_2 in the orientations *θ*_{p} = 0° and *θ*_{p} = 90° obtained by
`ADDA`–LDR and
`mieschka`. Remember that the corresponding reciprocity errors *e*_{hh} are 1.3 × 10^{−3}% for
`ADDA`–LDR and 8.0×10^{−3}% for
`mieschka`. They are both on the same order. The
`mieschka` DSCS can be regarded as the reference curves since its relatively small reciprocity error correlates with a correspondingly convergent differential scattering behavior as shown in Section 4.1.1. Later we will come back again to this example and will see yet another indication that this assumption is justified. In Fig. 5, larger differences are observed in the side- and backward scattering directions for both orientations. The latter is crucial in LIDAR applications, for instance. That is, the
`ADDA`–LDR reciprocity error is much smaller than the error in DSCS, compared to
`mieschka`. Note that Gasteiger et al. [48] have recently applied reciprocity tests to their
`ADDA` results. They do not give any details, but the indirect conclusion is that the reciprocity error is also smaller than the DSCS error. Note furthermore that we have also found cases where a good DSCS agreement with
`mieschka` has been obtained.

So the cases where the reciprocity is not automatically fulfilled are more interesting. They offer the possibility to investigate whether some DSCS accuracy information can be drawn out of the reciprocity error. In doing so, let us have a look at the CLDR polarizability. It depends on the propagation direction **a** of the incident field as LDR (see Eq. (17)). Table 2 shows the reciprocity errors obtained by
`DDSCAT` with GKDLDR (CLDR). The reciprocity condition (5) is well fulfilled for all scatterers at the smaller size parameter *k*_{0}*r _{v}* of 3. The same holds also for the particles with

*k*

_{0}

*r*= 15 and the smaller refractive index

_{v}*m*= (1.313, 5.889 × 10

^{−10}). The reciprocity errors range between 0% (within the given output accuracy of

`DDSCAT`up to the third digit after the decimal point) and about 5.5% for these cases. More critical are, however, the cases of scatterers having

*k*

_{0}

*r*= 15 and

_{v}*m*= (1.6, 0.0005). Here,

*e*

_{vv}is about 8.4% for cyl_2, 44% for cube_2, and 48% for sphd_p_2.

*e*

_{hh}ranges between about 18% for cyl_2 and 280% for sphd_o_2. Comparisons between the results of

`DDSCAT`–GKDLDR and

`mieschka`for the latter particle have given larger differences in the scattering angle range of about 60° – 120°, consequently in 240° – 300°, and in the backscattering region between 170° – 180° for the orientation

*θ*

_{p}= 0°. In contrast to this, a relatively good agreement has been found at

*θ*

_{p}= 90°. It can therefore be assumed that this large reciprocity error is caused by problems of

`DDSCAT`–GKDLDR in the former orientation.

A comparison of the DDA results in Table 2 shows that the use of LDR gives, in most cases considered, smaller reciprocity errors than CLDR (GKDLDR). This indicates that the reciprocity is not automatically fulfilled and, consequently, that there is a relation between the reciprocity error and DSCS for CLDR at the given particle configurations. Indeed, investigations by use of
`ADDA`–CLDR have shown a correlation between the reciprocity error and the number of dipoles used in the computations. Fig. 6 provides an example for the scatterer sphd_o_1 in the two configurations (*θ*_{p} = 0°, *θ* = 90°) and (*θ*_{p} = 90°, *θ* = 270°). It is seen that *e*_{vv} and *e*_{hh} decrease with an increasing dipole number. As already mentioned above, an increase of the dipole number should in turn lead to improved scattering cross sections. This is demonstrated by Fig. 7(a). It shows a comparison of the hh-polarized DSCS for this particle in the orientation *θ*_{p} = 0° obtained by using
`ADDA`–CLDR at different discretizations and by
`mieschka`. The
`ADDA`–CLDR results approach that of
`mieschka` when refining the discretization from 32 to 128 dipoles in *x*-direction. A further improvement can be achieved by applying the extrapolation technique described in Section 3.2.1 (see blue curve in Fig. 7(a)). So we have here a case where the reciprocity error correlates with the DSCS accuracy.

One can further think about determining the dipole number from reciprocity investigations in such cases. But additional studies are needed to gain reliable and more generally valid criteria to do so, due to the involved interplay of the different polarizability models and scattering configurations. Moreover, there are certain fundamental limitations of this approach. In particular, the DDA error for any scattering quantity (e.g. that for DSCS) can be divided into so-called shape and discretization parts [39]. The former is defined as the difference between exact results for real (smooth) scatterer and that for cubically–discretized one (i.e. a union of cubical dipoles). The result of the (standard) DDA, and hence the reciprocity error, is identical for both shapes. Therefore, the reciprocity error can not directly correlate with shape error, but only with discretization one. And the interrelation between these two parts of DDA errors is not trivial. While for a given scatterer both errors correlate with dipole size, and hence with each other, the quantitative relationship largely depends on a particular scatterer. For instance, total error equals the discretization one for a cube, but order-of-magnitude difference is easily possible for a sphere [39].

Anyway, reliable DDA results can always be obtained by increasing the dipole number and/or applying the extrapolation technique. To demonstrate this, let us come back to the above example of the scatterer sphd_o_2 with LDR of Fig. 5. Fig. 7(b) shows its hh-polarized DSCS for the orientation *θ*_{p} = 0° obtained by using the programs
`ADDA`–LDR at different discretizations and
`mieschka`. It is seen that an increase of the dipole number from 88, originally used in Fig. 5, to 256 leads to a convergence of the
`ADDA` curve towards that of
`mieschka`. The differential scattering behavior differs only in the peak at about *θ* = 165° and in the backscattering region. The blue curve represents the extrapolated result. Its agreement with
`mieschka` is comparable to the former 256 grid curve. Larger differences are seen only in the peak around *θ* = 105° and again in the backscattering region. Note that the convergence of the DDA results towards that of
`mieschka` is another indication to justify the use of the latter DSCS as reference, assumed above.

## 5. Summary and conclusions

In this paper, the fulfillment of the reciprocity by different T-matrix and DDA implementations for various sets of scatterers has been investigated. In most cases, the reciprocal configurations (*θ*_{p} = 0°, *θ* = 90°) and (*θ*_{p} = 90°, *θ* = 270°) have been considered. Although we focused on specific DDA and T-matrix implementations, we believe that our study should be seen in a broader context; the results help to establish the reciprocity condition as a useful accuracy test in light scattering computations.

The reciprocity error (7) obtained by the T-matrix program
`mieschka` ranges between about 3 × 10^{−4}% and nearly 15%. In general, the reciprocity is better satisfied for particles having a smooth, continuously differentiable boundary surface such as spheroids and Chebyshev particles than for scatterers with edges like finite circular cylinders. Except for *e*_{vv} of cyl_2 with the 15%, the reciprocity error is always below the required accuracy of 5%, chosen within the convergence procedure of
`mieschka`. The error of 15% results from the configurations considered, leading to investigations near deep down spikes of DSCS. Since it is hard for most scattering programs to achieve accurate DSCS values in such regions, other configurations should be taken into account to circumvent these difficulties. This becomes more important for particles with larger size parameters and/or refractive indices due to their complex differential scattering behavior.

Furthermore, a correlation between the reciprocity error, on one hand, and the relative error of the Barber–Hill convergence criterion and the series truncation parameter *n*_{cut}, respectively, has been found. Since the latter two parameters determine the relative convergence and consequently the accuracy of DSCS, reciprocity checks represent an appropriate tool to evaluate the far field scattering behavior obtained by
`mieschka` for a given particle.

The reciprocity errors obtained by the programs
`scsmfo1b` and
`mstm` for the bispheres are within the range of about 1.4 × 10^{−2}% to 5.8%. Similar to
`mieschka`, a relation between the reciprocity error and the convergence parameter *translation epsilon* of the newer program
`mstm` has been found. The latter parameter determines the maximum order of the cluster T-matrix and, in this way, the number of the scattered field expansion terms. Therefore, the reciprocal error can be regarded as an indicator for the DSCS accuracy of
`mstm`.

Only very little has been mathematically proven concerning the convergence of T-matrix approaches for nonspherical particles. Nevertheless, above results seem to indicate a general behavior of these methods. An increase of the number of field expansion terms, either directly or indirectly via corresponding convergence parameters, leads to an improvement of the relative convergence of DSCS, provided that no numerical instabilities arise. At the same time, the reciprocity error decreases. Thus the latter seems always to represent a measure for the DSCS accuracy of these methods.

The variation of the reciprocity error obtained by DDA for different polarizability models is essentially larger than that of above T-matrix programs. The error ranges between about 6.3 × 10^{−6}% and 280%. One reason for such small reciprocity errors is the fact that the reciprocity can be fulfilled a priori. These cases occur either if the employed polarizability model is independent of the incident field (such as RRC) or it becomes independent of the incident field for specific scattering configurations (as sometimes with LDR, e.g.). The numerical accuracy of the reciprocity fulfillment is, then, solely determined by the solution scheme applied to the DDA equation system. So, blindly applying reciprocity tests to DDA can lead to false conclusions.

The reciprocity criterion is only applicable in conjunction with polarizability models that are dependent on the incident field. The CLDR model has met this requirement for all the particles and reciprocal configurations considered in this paper. Investigations by means of
`ADDA`–CLDR have revealed a correlation between the reciprocity error and the number of dipoles used in the computations. It could further be shown that the increase of the dipole number leads, in turn, to improved scattering cross sections. Consequently, reciprocity checks can potentially provide information about the DSCS accuracy of DDA.

One of the main aim of this study was to identify programs and situations in which the reciprocity fulfillment can be used as an criterion to evaluate scattering results. Additional investigations for a multitude of different reciprocal configurations and particle parameters are now needed to derive ready-to-use criteria to quantitatively estimate the relative error of DSCS on the basis of reciprocity tests. Concerning DDA, this can become intricately due to the involved interplay of the different polarizability models and scattering configurations. But if available, both users and program developers of scattering software can benefit from them. They can assist a user in deciding whether the program at hand is appropriate for his needs. On the other hand, they can help a developer to find errors during coding or reveal merits and demerits of implementations.

## Acknowledgments

We would like to thank B. T. Draine, P. J. Flatau, D. W. Mackowski, K. Fuller, and M. Mishchenko for making their codes available and for allowing the publication of results obtained with them. Furthermore, we would like to thank T. Rother for helpful discussions. M. Yurkin acknowledges support of the program of the Russian Government “Research and educational personnel of innovative Russia” (Contract Nos. P422, P1039, and 14.740.11.0921), by grant from Russian Government 11.G34.31.0034, and by Russian Foundation of Basic Research (grant 12-04-00737-a). M. Kahnert acknowledges funding from the Swedish Research Council (Vetenskapsrådet) under projects 621-2008-4387 and 621-2011-3346.

## References and links

**1. **M. I. Mishchenko, J. W. Hovenier, and L. D. Travis, eds., *Light Scattering by Nonspherical Particles: Theory, Measurements, and Applications* (Academic, 2000).

**2. **M. I. Mishchenko, L. D. Travis, and A. A. Lacis, *Scattering, Absorption, and Emission of Light by Small Particles* (Cambridge University Press, 2002).

**3. **F. Borghese, P. Denti, and R. Saija, *Scattering from Model Nonspherical Particles: Theory and Applications to Environmental Physics*(Springer, 2003).

**4. **M. Kahnert, “Numerical methods in electromagnetic scattering theory,” J. Quant. Spectrosc. Radiat. Transfer **79–80**, 775–824 (2003). [CrossRef]

**5. **A. Doicu, T. Wriedt, and Y. A. Eremin, *Light Scattering by Systems of Particles. Null-Field Method with Discrete Sources: Theory and Programs* (Springer, 2006).

**6. **A. A. Kokhanovsky, ed., *Light Scattering Reviews*, Vol. 1 (Springer, 2006). [CrossRef]

**7. **A. A. Kokhanovsky, ed., *Light Scattering Reviews*, Vol. 2 (Springer, 2007). [CrossRef]

**8. **A. A. Kokhanovsky, ed., *Light Scattering Reviews*, Vol. 3 (Springer, 2008). [CrossRef]

**9. **A. A. Kokhanovsky, ed., *Light Scattering Reviews*, Vol. 4 (Springer, 2009). [CrossRef]

**10. **T. Rother, *Electromagnetic Wave Scattering on Nonspherical Particles: Basic Methodology and Simulations* (Springer, 2009). [CrossRef]

**11. **M. Kahnert, “Electromagnetic scattering by nonspherical particles: recent advances,” J. Quant. Spectrosc. Radiat. Transfer **111**, 1788–1790 (2010). [CrossRef]

**12. **T. Wriedt and J. Hellmers, “New scattering information portal for the light–scattering community,” J. Quant. Spectrosc. Radiat. Transfer **109**, 1536–1542 (2008). [CrossRef]

**13. **J. Hellmers and T. Wriedt, “New approaches for a light scattering internet information portal and categorization schemes for light scattering software,” J. Quant. Spectrosc. Radiat. Transfer **110**, 1511–1517 (2009). [CrossRef]

**15. **J. W. Hovenier and C. V. M. van der Mee, “Basic relationships for matrices describing scattering by small particles,” in *Light Scattering by Nonspherical Particles: Theory, Measurements, and Applications*, M. I. Mishchenko, J. W. Hovenier, and L. D. Travis, eds. (Academic, 2000), pp. 61–85. [CrossRef]

**16. **T. Rother and J. Wauer, “Case study about the accuracy behavior of three different T-matrix methods,” Appl. Opt. **49**, 5746–5756 (2010). [CrossRef] [PubMed]

**17. **H. C. van de Hulst, *Light Scattering by Small Particles* (John Wiley & Sons, 1957).

**18. **P. Barber and S. C. Hill, *Light Scattering by Particles: Computational Methods* (World Scientific, 1990).

**19. **J. Wauer, K. Schmidt, T. Rother, T. Ernst, and M. Hess, “Two software tools for the plane–wave scattering on nonspherical particles in the German Aerospace Center’s virtual laboratory,” Appl. Opt. **43**, 6371–6379 (2004). [CrossRef] [PubMed]

**20. **D. W. Mackowski, K. Fuller, and M. I. Mishchenko, “Codes for calculation of scattering by clusters of spheres,” ftp://ftp.eng.auburn.edu/pub/dmckwski/scatcodes/index.html.

**21. **D. W. Mackowski and M. I. Mishchenko, “A multiple sphere T-matrix Fortran code for use on parallel computer clusters,” J. Quant. Spectrosc. Radiat. Transfer **112**, 2182–2192 (2011). [CrossRef]

**22. **M. Kahnert, T. Nousiainen, H. Lindqvist, and M. Ebert, “Optical properties of light absorbing carbon aggregates mixed with sulphate: assessment of different model geometries for climate forcing calculations,” Opt. Express20, 10042–10058 (2012), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-9-10042. [CrossRef] [PubMed]

**23. **L. Tsang, J. A. Kong, and R. T. Shin, *Theory of Microwave Remote Sensing* (John Wiley & Sons, 1985).

**24. **M. Kahnert, “Irreducible representations of finite groups in the T-matrix formulation of the electromagnetic scattering problem”, J. Opt. Soc. Am. A **22**, 1187–1199 (2005). [CrossRef]

**25. **A. Mugnai and W. J. Wiscombe, “Scattering of radiation by moderately nonspherical particles,” J. Atmos. Sci. **37**, 1291–1307 (1980). [CrossRef]

**26. **H. Volten, O. Muñoz, J. W. Hovenier, J. F. de Haan, W. Vassen, W. J. van der Zande, and L. B. F. M. Waters, “WWW scattering matrix database for small mineral particles at 441.6 and 632.8 nm,” J. Quant. Spectrosc. Radiat. Transfer **90**, 191–206 (2005). [CrossRef]

**27. **S. G. Warren and R. E. Brandt, “Optical constants of ice from the ultraviolet to the microwave: a revised compilation,” J. Geophys. Res. **113**, D14220 (2008). [CrossRef]

**28. **T. Rother, K. Schmidt, J. Wauer, V. Shcherbakov, and J.-F. Gayet, “Light scattering on Chebyshev particles of higher order,” Appl. Opt. **45**, 6030–6037 (2006). [CrossRef] [PubMed]

**29. **T. Rother, “Generalization of the separation of variables method for nonspherical scattering on dielectric objects,” J. Quant. Spectrosc. Radiat. Transfer **60**, 335–353 (1998). [CrossRef]

**30. **P. C. Waterman, “Symmetry, unitarity, and geometry in electromagnetic scattering,” Phys. Rev. D **3**, 825–839 (1971). [CrossRef]

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

**32. **D. W. Mackowski, “Calculation of total cross sections of multiple sphere clusters,” J. Opt. Soc. Am. A **11**, 2851–2861 (1994). [CrossRef]

**33. **The MSTM package is available at http://www.eng.auburn.edu/users/dmckwski/scatcodes.

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

**35. **D. Gutkowicz-Krusin and B. T. Draine, “Propagation of electromagnetic waves on a rectangular lattice of polarizable points,” (2004), http://arxiv.org/abs/astro-ph/0403082.

**36. **M. A. Yurkin and A. G. Hoekstra, “The discrete–dipole–approximation code ADDA: capabilities and known limitations,” J. Quant. Spectrosc. Radiat. Transfer **112**, 2234–2247 (2011). [CrossRef]

**37. **M. A. Yurkin and A. G. Hoekstra, “The discrete dipole approximation: an overview and recent developments,” J. Quant. Spectrosc. Radiat. Transfer **106**, 558–589 (2007). [CrossRef]

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

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

**40. **M. A. Yurkin, D. de Kanter, and A. G. Hoekstra, “Accuracy of the discrete dipole approximation for simulation of optical properties of gold nanoparticles,” J. Nanophoton. **4**, 041585 (2010). [CrossRef]

**41. **M. A. Yurkin, A. G. Hoekstra, R. S. Brock, and J. Q. Lu, “Systematic comparison of the discrete dipole approximation and the finite difference time domain method for large dielectric scatterers,” Opt. Express15, 17902–17911 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-26-17902. [CrossRef] [PubMed]

**42. **M. A. Yurkin, M. Min, and A. G. Hoekstra, “Application of the discrete dipole approximation to very large refractive indices: filtered coupled dipoles revived,” Phys. Rev. E **82**, 036703 (2010). [CrossRef]

**43. **B. T. Draine and P. J. Flatau, “Discrete–dipole approximation for scattering calculations,” J. Opt. Soc. Am. A **11**, 1491–1499 (1994). [CrossRef]

**44. **http://www.astro.princeton.edu/draine/DDSCAT.html

**45. **http://ddscat.wikidot.com/start

**46. **B. T. Draine and P. J. Flatau, “User guide to the discrete dipole approximation code DDSCAT 7.1,” (2010), http://arXiv.org/abs/1002.1505v1.

**47. **X. A. Penttila, E. Zubko, K. Lumme, K. Muinonen, M. A. Yurkin, B. T. Draine, J. Rahola, A. G. Hoekstra, and Y. Shkuratov,“Comparison between discrete dipole implementations and exact techniques,” J. Quant. Spectrosc. Radiat. Transfer **106**, 417–436 (2007). [CrossRef]

**48. **J. Gasteiger, M. Wiegner, S. Groß, V. Freudenthaler, C. Toledano, M. Tesche, and K. Kandler, “Modelling lidar–relevant optical properties of complex mineral dust aerosols,” Tellus B **63**, 725–741 (2011). [CrossRef]