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

 

Fig. 1 Configurations of the bisphere with touching components ((a) and (b)) and of the cube ((c) and (d)) where the incident field E_inc propagates along the positive z-axis. The y-axis points to the reader. (a) The rotational axis is oriented along the positive z-axis (θ_p = 0°). The scattered field E_sca is observed in positive x-direction (θ = 90°). (b) The rotational axis is oriented along the positive x-axis (θ_p = 90°). The scattered field E_sca is observed in negative x-direction (θ = 270°). (c) The central axis, which stands perpendicularly on the top and bottom faces of the cube while going through its center, is oriented along the positive z-axis (θ_p = 0°). The scattered field E_sca is observed in positive x-direction (θ = 90°) out of a cube edge. (d) The central axis is oriented along the positive x-axis (θ_p = 90°). The scattered field E_sca is observed in negative x-direction (θ = 270°).

Download Full Size | PDF

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

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₀r_v and 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/2r. Furthermore, n and ε are the order and the deformation parameter, respectively, of the Chebyshev particles, the shape of which is given by

Table 1. Scatterer geometries considered in the investigations.

View Table | View all tables in this article

To compute the radius r₀ 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⁻¹⁰) is the refractive index of ice at a wavelength of 500 nm [27]. It is obvious that the combination of k₀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Ψ_τnl, respectively.

Here, k_s denotes the wave number inside the scatterer. To compute the unknown scattered and internal field expansion coefficients (f_τnl, p_τnl), 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

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, qeps1, qeps2) 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, qeps1 the error tolerance for determining the single-sphere harmonic order truncation, and qeps2 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, qeps1, qeps2) = (400, 1 × 10⁻⁶, 0, 1 × 10⁻⁵, 1 × 10⁻⁸), 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⁻⁴, 1 × 10⁻³, 1 × 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

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)

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]

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⁻⁵ 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⁻⁵, and dipole size d = max(λ/(10|m|), D_x/16), where D_x is the size of the particle along the 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₀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⁻⁵, 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.

 

Fig. 2 vv-polarized DSCS of the scatterer cyl_2 in the orientations θ_p = 0° (black line) and θ_p = 90° (red line) calculated with mieschka. The corresponding reciprocal scattering angles θ = 90° and θ = 270° are marked by vertical lines.

Download Full Size | PDF

Table 2. Reciprocity errors e_hh and e_vv, defined by Eq. (7), for the reciprocal configurations (θ_p = 0°, θ = 90°) and (θ_p = 90°, θ = 270°) using mieschka, ADDA with LDR polarizability, and DDSCATwith GKDLDR (CLDR) polarizability.

View Table | View all tables in this article

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⁻⁴% and 3.1% for the former particles whereas it is within 1.2 × 10⁻¹% – 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.

 

Fig. 3 Reciprocity errors e_vv and e_hh according to Eq. (7) for sphd_p_2 in the two reciprocal configurations (θ_p = 0°, θ = 90°) and (θ_p = 90°, θ = 270°) using mieschka (a) versus the relative error within the chosen Barber–Hill convergence criterion, and (b) versus n_cut. (n_cut, l_cut, n_int) are given by Table 3.

Download Full Size | PDF

Table 3. Convergence parameters (n_cut, l_cut, n_int) determined automatically in the standard mode of mieschka for sphd_p_2 in the orientation θ_p = 90° at different input values of the relative error within the chosen Barber–Hill convergence criterion. Note that the relative integration error has been fixed to 0.1%.

View Table | View all tables in this article

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

Table 4. Reciprocity errors e_hh and e_vv, defined by Eq. (7), for the reciprocal configurations (θ_p = 0°, θ = 90°) and (θ_p = 90°, θ = 270°) using scsmfo1b and mstm for the bispheres.

View Table | View all tables in this article

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⁻⁷. A further decrease of this parameter beyond 1 × 10⁻⁷ 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.

 

Fig. 4 Reciprocity errors e_vv and e_hh according to Eq. (7) for bi_sph_2 in the two reciprocal configurations (θ_p = 0°, θ = 90°) and (θ_p = 90°, θ = 270°) using mstm versus translation_epsilon.

Download Full Size | PDF

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⁻⁶% for sphd_p_4 and e_vv = 1.2 × 10⁻¹% 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_j). It obeys the symmetry relation Ḡ (r_i, r_j) = Ḡ (r_j, r_i). 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 Ḡ (r_i, r_j) 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 Ḡ (r_i, r_j), 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.

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⁰ of the incident electric field in the body system (see Eq. (16)). If either a or e⁰ is along a coordinate axis, then S = 0. As a consequence, the symmetry effect of Ḡ (r_i, r_j) is not extinguished so that the reciprocity condition is fulfilled a priori. In the 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.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⁻³% and 1.6×10⁻²%, 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⁻⁵. The RRC model is independent of the incident field, i.e. of both a and e⁰ (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⁻³% while e_hh = 1.2 × 10⁻³% 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⁻³% as for 16 dipoles. A decrease of the threshold of the iterative solver down to 10⁻¹⁰ has only a small effect on the DSCS values themselves, but reduces the reciprocity errors down to the order of 10⁻⁹% (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.

Table 5. Polarized DSCS and reciprocity errors e_hh and e_vv, defined by Eq. (7), for the reciprocal configurations c₁ = (θ_p = 0°, θ = 90°) and c₂ = (θ_p = 90°, θ = 270°) using ADDA with RRC polarizability for cube 1 at different program parameters.

View Table | View all tables in this article

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⁻³% for ADDA–LDR and 8.0×10⁻³% 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.

 

Fig. 5 hh-polarized DSCS of the scatterer sphd_o_2 in the orientations θ_p = 0° and θ_p = 90° calculated using ADDA with LDR polarizability (black line) and mieschka (red line).

Download Full Size | PDF

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₀r_v of 3. The same holds also for the particles with k₀r_v = 15 and the smaller refractive index m = (1.313, 5.889 × 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₀r_v = 15 and 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.

 

Fig. 6 Reciprocity errors e_vv and e_hh according to Eq. (7) for sphd_o_1 in the two reciprocal configurations (θ_p = 0°, θ = 90°) and (θ_p = 90°, θ = 270°) using ADDA with CLDR polarizability versus discretization level.

Download Full Size | PDF

 

Fig. 7 hh-polarized DSCS of the scatterer (a) sphd_o_1 and (b) sphd_o_2 in the orientation θ_p = 0° calculated using ADDA at different discretizations and mieschka (green line). (a) ADDA with CLDR polarizability (black line – 32 dipoles in x-direction, red line – 128 dipoles in x-direction, blue line – extrapolated) (b) ADDA with LDR polarizability (black line – 88 dipoles in x-direction, red line – 256 dipoles in x-direction, blue line – extrapolated)

Download Full Size | PDF

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⁻⁴% 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⁻²% 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⁻⁶% 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.