Light scattering theories and computer codes

In aerosol science today light scattering simulations are regarded as an indispensable tool to develop new particle characterization techniques or in solving inverse light scattering problems. Light scattering theories and related computational methods have evolved rapidly during the past decade such that scattering computations for wavelength sized nonspherical scatterers can be easily performed. This signi ﬁ cant progress has resulted from rapid advances in computational algorithms developed in this ﬁ eld and from improved computer hardware. In this paper a review of the recent progress of light scattering theories and available computational programs is presented. Short outlines of the various theories is given alongside with informations on is capabilities and restrictions.


Introduction
In aerosol science there is increasing interest not only to determine the size of aerosol particles using optical methods but also the shape or the refractive index of the particles. This can help, for example, to find the source of a particle contamination. To develop new and efficient optical particle characterization methods light scattering simulation methods are needed which can handle a wide variety of particle shapes. In this paper I would like to review the state of the art in this field as it developed over the last decade. I will give a short description of the various theories available and provide information on computation programs.
A number of light scattering theories have been developed, all having their pros and cons. Extensive overviews of available theories have been published by Wriedt [1], Kahnert [2] and more recently by Veronis and Fan [3]. A review of this subject with emphasis on colloid science has been written by Niklasson and Vargas [4]. A paper by Zhao et al. [5] provides an overview of the methods that are currently being used to study the electromagnetics of silver and gold nanoparticles.
Of historic interest might be the critical review by Bouwkamp [6] who presents the progress in classical diffraction theory up 1954. In this paper both scalar and electromagnetic problems are discussed. This report may also serve as an introduction to general diffraction theory.
Interesting reviews on related subjects of nanooptics and metamaterials have recently been published by Myroshnychenko et al. [7] and by Veselago et al. [8].
With the progress of the internet in the 1990s two web sites were erected to provide information on and links to available light scattering codes: Electromagnetic Scattering Programs [9] and SCATTERLIB [10]. With both lists there is an emphasis on open access to computational programs. In the near future this web information will be expanded by a New Scattering Information Portal for the light-scattering community [11]. Most of the mentioned computational programs can be found on either of this two web pages.
For the purpose of this review I will divide the available numerical methods used in electromagnetics into three main categories.
Analytical Unter the heading Differential equation methods two methods are included Finite difference time domain method and Finite Element Method. With these methods the differential equation must be discretized over the entire volume of the scatterer and part of the surrounding medium. Thus both methods need a special boundary condition to truncate the computational domain. Differential equation methods lead to matrices that are highly sparse reducing computer storage.
Integral equation methods are based on the discretization of an integral equations using the Green's function of the surrounding medium. With this methods the radiation conditions for the scattered field is automatically satisfied. The matrix to be solved is a full matrix.
Next, I will shorty review the methods available making reference to review papers where available. For more profound reviews the interested reader is referred to the cited review articles.

Mie Theory and its extensions
In a paper within a special edition commemorating Gustav Mie and the Mie theory a short review of Mie theory and its extensions should not be missing. In 1908 Gustav Mie published his famous paper on simulation of the colour effects connected with colloidal Gold particles [12]. In this paper he gave a first outline of how to compute light scattering by small spherical particles.
Cardona and Marx [13] comment, that this paper was almost ignored until about 1945 but its importance rose starting from the 1950s with interest in colloids such that hundred years after publication his paper is still much cited [14] with currently 160 citations a year. The paper is call Dornröschen (Sleeping Beauty) by the researchers at Information Retrieval Services of Max Planck Society [15] because of its late recognition considering the number of citations.
The citations of this paper show that the applications of Mie theory cover a wide range of subjects from climate research, optical particle characterization to interstellar dust in astrophysics.
Understandably prior to the development of electronic computers in the middle of last century there were not many paper written on computing scattering problems using Mie's theory since the computational labour involved in evaluating functions such as Ricatti-Bessel functions was quite extreme.
Even with the rise of the computer it took some time before stable algorithms were developed. Gradually several generally reliable and stable scattering programs were published. Early well known algorithms were published by Giese [16] and Dave [17]. The IBM report by Dave from 1969 was still send out on request in the 1990s. Nowadays a number of efficient algorithms and FOR-TRAN programs are available. A major step was the program MIEV0 written by Wiscombe [19,18], which is based on Lentz's continued-fraction method for the calculation of spherical Bessel functions [20]. The program is well-tested and widely used. It has been demonstrated that Mie's theory can be now successfully applied up to size parameters of 10.000 [21,22,23] and this year Mie theory has even been programmed on a Java enabled Mobile Phone using a Matlab clone [24]. As the scattering of a plane electromagnetic wave by a dielectric sphere is considered a canonical problem, Mie's theory is still used as a standard reference to validate methods intended for more complex scattering problems [25,26].
Apparently As it is easy to consider spherical scatterers there are many extensions of Mie theory covering different aspects. A theory for a coated dielectric sphere was first published by Aden and Kerker [32]. An advanced algorithm is given by Toon and Ackerman [33]. An algorithm for a sphere having two coatings has been presented by Kaiser and Schweiger [34]. These theories have also been extended to spherical particles consisting of multiple layers by Li Kai [35]. Even today improved algorithms are published on this subject [36].
A scattering sphere can also be chiral or optically anisotropic. Theories and programs for such type of scatterers have been published by Bohren in his book for a chiral sphere [38] and by Doicu et al. [37] for anisotropic spheres.
In optical particle characterization scattering by a focused laser beam is quite often of interest. There are different concepts available to handle this problem. A laser beam with Gaussian intensity distribution can be expanded into spherical vector wave functions or into a spectrum of plane waves [39]. The Generalized Lorenz Mie Theory (GLMT) developed by G. Gouesbet and coworkers is based on the first approach. Is has recently been reviewed by Gouesbet [40]. Plane wave expansion is used by Albrecht at al. [41] and it is integrated into the Null-Field Methods with Discrete Sources (NFM-DS) developed by and Doicu and coworkers [37]. Scattering by an aggregate of spheres will also be considered in this sections. There are efficient programs available by Mackowski [42] and Yu-lin Xu [43]. The theory has recently been extended to clusters of rotationally symmetric particles [44] and arbitrary shaped particles [45].
Mie and absorbing surrounding medium has been considered by Sudiarta [46] and by Frisvad et al. [47].
As colour pigments and crystals are anisotropic there is some interest to extend Mie's theory to such kind of scattering particles. Stout et al. [48] established a vector spherical harmonic expansion of the electromagnetic field propagating inside an arbitrary anisotropic medium to solve this problem. A similar problem of an uniaxial anisotropic sphere was considered by Geng et al. [49] to find the coefficients the scattered field. Droplets of liquid crystals may be considered as spherical particle with radial anisotropy. This scattering problem is discussed by Qui at al. [50] using an extension of Mie theory.
Especially with nanosized nobel metal particle of size lower than about 20nm various modifications, extensions and corrections to Mie's original theory are needed to take into account that "sharp" boundary conditions do not hold in the nanoscale. [51]. In his recent survey paper [52] Kreibig lists among others the following supplementary models to the Mie theory, incident wave not plane, non step-like boundary condition, dielectric function dependent on particle size. Applying this extensions help to explain measured absorption spectra of Ag nanoparticles and plasmon polaritons.

Separation of Variables Method
Asano and Yamamoto [53] were apparently were the first who used separation of the vector wave equation in spheroidal coordinates to consider light scattering by homogeneous prolate or oblate spheroidal particles with an arbitrary size. In this way all fields are expanded in terms of spheroidal wavefunctions. This method using spheroidal wavefunctions for field expansion is commonly termed Separation of Variables Method (SVM). Voshchinnikov and Farafonov [54] applied a related approach. Their solution is claimed to be more efficient than the one of Asano and Yamamoto from the computational point of view. The approach has especially an advantage for strongly elongated or flattened particles.
Voshchinnikov's SVM codes based on spheroidal expansion are available from Il'in's website of scattering codes [55]. Scattering by coated spheroids can also computed using SVM [56]. A recursive algorithm for the solution of the problem of scattering by multilayered spheroidal particles has been presented by Farafonov [57]. The method has been reviewed by Farafonov et al. [58]. SVM has also been combined with the Generalized Lorenz-Mie theory to compute light scattering for by a homogeneous spheroid in a shaped laser beam [59].
Suitable Mathematica programs using Spheroidal Wave Functions to compute scattering by spheroidal particles have been published in a book by Le-Wei Li et al. [60]. These includes also programs for coated dielectric spheroids. For more references on this subject using spheroidal wave function I would like to refer the interested reader to this book. In a paper by Schulz [61] a method for computing the T-matrix with the SVM is given.

T-Matrix Method
The T-matrix method or null-field method is one of the most well known light scattering theories to compute scattering by nonspherical particles. A Fortran program for a perfectly electric conducting rotationally symmetric scatterer was published quite early by Waterman [62]. The book by Barber and Hill [63] includes some nice Fortran programs to compute scattering by spheroids. Near field and internal field can also be computed.
Recent reviews of the literature on this method have been published by Mishchenko et al. [64,65]. The advantage of this method is that the T-matrix is computed, which includes the full solution of the scattering problem. This matrix relates the expansion coefficients of the incident field to the expansion coefficients of the scattered field. By using a precomputed T-matrix, scattering by a rotated or a translated particle or orientation averaged scattering can easily be computed. Additionally, multiple scattering or scattering by a particle located near a plane interface can be computed from a stored T-matrix of this particle.
Most well known, tested and widely used are the T-Matrix programs by Mishchenko [66]. The standard T-matrix method is restricted to particles having an aspect ratio not larger than about 1:4 [67]. To solve the stability problems with particles having a larger aspect ratio expansion of the internal field using discrete sources was introduced into the null-field method. This method is reviewed in one of the following sections of this paper.
The Optical tweezers toolbox by Timo Nieminen [68], implemented in Matlab, is based on the T-Matrix method for the computational modelling of optical tweezers. The toolbox is designed for the calculation of optical forces and torques, and can be used for both spherical and nonspherical particles, in both Gaussian and other beams.

Generalized Multipole Technique
The Generalized Multipole Technique (GMT) is a relatively new and fast advancing method which has been developed by different research groups. Ludwig [69] coined the term Generalized Multipole Technique for this spectrum of methods.
The main idea of the Generalized Multipole Technique consists in approximating the solution of the problem by a linear combination of discrete sources. These discrete sources are the fundamental solution of the differential equation of the problem. The introduction of discrete sources method is generally attributed to Kupradze and Aleksidze [70]. The coefficients of the scattered field are found by satisfying the boundary conditions on the surface of the scatterer by combining point matching with a least squares technique.
There are a variety of methods which use "equivalent sources" for field expansion. The "equivalent sources" may be of any type as long as they are solutions of the wave equation. Spherical waves, dipoles, Mie potentials and Gabor functions have been applied for field expansion.
A review of these methods has been published in an edited volume by Wriedt [71] and more recently by Fairweather [72]. A comprehensive presentation of the discrete sources method is given by Doicu et al. [73]. For an analytic foundation of the null field method and the discrete sources method in acoustics and electromagnetics the reader is referred to this book.
The Fortran code of the Multiple Multipole Program (MMP) has been published by Hafner and Bomholt [74]. The book is out of print but an executable is available [75]. Boundary methods for solving Maxwell equations with a focus on the Multiple Multipole Program (MMP) are outlined by Hafner [76]. The 5th order Gaussian beam can also quite easily be integrated into the Multiple Multipole Program [77].
The Discrete Sources Method (DSM) is developed by Yuri Eremin and coworkers [78]. In this method the unknown amplitudes of the discrete sources are determined from the boundary condition in least square sense. DSM is retricted to rotationally symmetric scatterers but the discrete sources can also be positioned in complex plane such that scattering by really flat discs can be computed [79]. A high aspect ratio of 1:50 for elongated particles and 20:1 for flat particles can be achieved [80]. Using DSM scattering by red blood cell can be computed in an efficient way [81]. are positioned on the axis of symmetry of the particle [82]. With oblate rotationally symmetric scatterers it is of advantage to position the discrete sources in a complex plane [83]. In this way, scattering by oblate discs [84], flat Cassini shaped particles [85] has been solved. Scattering by particles having an aspect ratio as large as 100:1 can easily be computed both for extremely oblate or prolate particles. For a full description of the theory we refer the interested reader to the recently published monograph by Doicu et al. [37] which also includes Fortran90 programs on CD-ROM. Recently the method has been extended to compute scattering by biaxially anisotric particles [86]. A sample of a computational result is shown in Fig. 1. The full scope of the Null-Field Method with Discrete Sources has recently been presented in a review paper [87].

Finite Difference Time Domain
The most common electromagnetic modelling technique for scattering simulations is the Finite Difference Time Domain method (FDTD). The Finite Difference Time Domain (FDTD) method is an electromagnetic modelling technique which is becoming more popular in computing light scattering by nonspherical or inhomogeneous particles. The Finite Difference Time Domain Algorithm has been introduced in 1966 by Yee [88].
For the application of FDTD the computational space and time is divided into a grid of discrete points and then the derivatives of the Maxwell equations are approximated by finite differences. With the FDTD method an entire volume including the scatterer is discretized into a cubical grid. The basic element of discretization is the Yee cell [88]. A difference approximation is applied to evaluate the space and time derivatives of the field. The fields are computed using a time marching scheme until a steady state solution is obtained. Not only the scattering particle has to be gridded but also some part of the surrounding medium. To truncate the computational domain the FDTD needs the use of absorbing boundary conditions that require a layer of grid cells all around the computational domain. The absorbing boundary conditions ensure that the wave is not reflected at the open boundary of the discretized volume.
A near-to far-field transformation is needed to compute the scattered farfield from the near-field values of the computational domain. FDTD can be applied to any kind of problem involving complex materials with arbitrary geometry, and it continues to be an active field of research [89].
Recently, it has been demonstrated that a GPU (Graphics Processing Unit) can be a cost effective co-processor [90]. Currently almost all vendors of FDTD programs offer an acceleration card based on a GPU.
There are about 27 commercial FDTD software vendors [91]. An other list of commercial programs is given in a review paper by Veselago et al. [8]. A survey of the Finite-Difference Time Domain literature has been published by Shlager and Schneider [92].
The Transmission Line Matrix (TLM) method is related to the FDTD method. It has similar capabilities. With TLM the nodes of the cubic grid are connected by virtual transmission lines. [98]. A number of commercial TLM packages is also listed in the review by Veselago et al. [8].

Finite Element Method
In light scattering applications the Finite Element Method (FEM) is not that well known. But you will nevertheless find some paper on applying this method. FEM is based on discretion of the total simulation volume including the scattering particle into a set of finite elements (tetrahedra or triangular prisms) which are best suited to describe the geometry of arbitrary shaped scattering particle. Then the unknown functions are represented by simple interpolation functions defined on each element with unknown coefficients. This functions provide an approximation to the solution within a finite element. Thus a system of equations is formulated which has to be solved.
To limit the region of computation like in FDTD the computational space has to be terminated using some artificial absorbing boundary. To solve the problem, the FEM makes use of a variational formulation [99]. The resulting matrix to be solved is sparse which is considered an advantage of this method [100]. A review of the finite element method for light scattering computations has been published by Violakis et al. [101].
Currently one of the most popular software packages implementing a 3D finite element method is Comsol Multiphysics. COMSOL Multiphysics (formerly FEMLAB) is a finite element analysis and solver software package for various physics and engineering applications [102]. Using this program the polarization dependence of greatly enhanced electromagnetic fields, in coupled gold nanoparticle-nanowire systems has recently been investigated [103]. The enhancement is strongest when the incident light is polarized across the junction of the particle wire system.

Method of Moments
In electrical engineering the Boundary Element Method (BEM) is commonly call Method of Moments (MoM). By MoM the surface integral equations are solved by discretization and reducing them to a system of linear equations. The development of the method is commonly attributed to Harringon [104]. The equation solved by MoM generally has the form of an electric field integral equation (EFIE), a magnetic field integral equation (MFIE) or combined field integral equation (CFIE) [105].
In this method the scatterer is replaced by a number of interacting equivalent electric or magnetic surface currents. Thus just the boundary surface has to be discretizing. A full matrix is generated which has to be solved to find the surface currents. From the equivalent currents on the surface of the scatterer all farfield quantities can be computed. Currently there is much progress with this method using various acceleration method such as the Fast Multipole Method (FMM). The fast multipole method (FMM) accelerates the matrix-vector product in the conjugate gradient (CG) method used to solve the matrix equation iteratively [106].
A personal history of the development of the Method of Moments as applied to electromagnetics is presented by Harrington [107].
Open FMM is a free collection of electromagnetic software for scattering at very large dielectric objects [108] based on an implementation by Fostier [109]. A FMM program is able to handle and solve very large linear systems. Thus scattering by objects as large as 250 λ in diameter has been solved using a parallel computer [110]. A number of commercial companies offering 3D MoM simulation package are listed in appendix F of Davidson book [111].

Volume Integral Equation Method
The Volume Integral Equation (VIE) may be solved by two related methods to compute scattering by an arbitrarily shaped or inhomogeneous scatterer [112]. Most well known and apparently most used is the Discrete Dipole Approximation (DDA) proposed by Purcell and Pennypacker [113]. In DDA an arbitrarily shaped particle is treated as a three-dimensional assembly of dipoles on a cubic grid. Each dipole cell is assigned a complex polarizability which can be com- puted from the complex refractive index of the bulk material and the number of dipoles in a unit volume. An extensive review of the method has recently been published by Yurkin and Hoekstra [114]. A Wiki on this method covering all aspects from theory to program has been prepared by Draine and Flatau [115].
A number of programs based on this related methods is freely available. This includes DDSCAT [115], ADDA [116], CG-FFT by Manuel F. Catedra, jscat by Jouni Peltoniemi and LIBGREEN by Nicolas Piller [117]. For the FORTRAN sources programs please have a look at the web site Electromagnetic Scattering Programs [9].
Also with this method there is a project to accelerate light scattering simulation by reconfigurable computing which intends to implement DDA in an FPGA (field-programmable gate array) chip [118]. A new release of the discretedipole approximation scattering code DDSCAT recently became available [115]. DDSCAT 7.0 contains a number of new features such as postprocessing to calculate the fields within or near a scattering particle. With this version absorption and scattering by targets containing arbitrarily-oriented anisotropic materials can be calculated.
To generate particle shape input data consisting of positions of the dipoles together with the index of refraction for DDA programes from 3d-model file formats based on a surface triangulation such as Wavefront .obj the "point in polyhedron" algorithm can be used [119]. Fig. 2 gives an example of a triangulated Fresnel-2 shape and Fig. 3 presents the corresponting DDA input data as a collection of spheres representing dipoles.

Validation
In development of computational programs in the field of light scattering validations is an important aspect. Mostly, results obtained from different computer programs are compared in the literature.
Examples using Discrete Sources Method (DSM) [120] for validation of the NFM-DS have been published for the flat discs [85], concave Cassini oval based shapes [84], long fibres [121].
Results for a dielectric cube comparing NFM-DS with DDA, FDTD [122] and VIEM [123] have been published by Wriedt and Comberg [124]. NFM-DS, DDA and MMP [74] have been compared for a two particle multiple scattering problem by Comberg and Wriedt [125]. NFM-DS, MMP, DSM, DDA and Finite Integration Technique (FIT) [126] have all been applied to compute scattering by a red blood cell [127]. Especially DSM proved to be an efficient method to compute scattering by the erythocyte. Fig. 4 shows the cross sections of two erythocytes of different shapes used for the sample results computed using DSM pesenented next. Fig. 5 presents the scattering diagram for incident ppolarisation and Fig. 6 for the corresponding s-polarisation. In each case the incident plane wave propagates along the axis of symmetry of the scattering erythocyte.
Recently, several computer implementations of the DDA method have been compared in terms of their accuracy, speed and usability. Additionally, the accuracy of the implementations have been studied by comparing results against results from Mie, T-matrix or cluster T-Matrix codes [128].
For validation microwave analog experiments can also be used. Sabouroux et al. [129] measured amplitude and phase of the scattered field of a cluster of spheres and the measurements compared with theoretical predictions obtained from a multiple-scattering T-matrix method and discrete dipole approximation (DDA). Excellent agreement demonstrates the validities of both the experiment and the models.
Grosges et al. [25] compared numerical scattering results in context of nearfield spectroscopy for FEM and FDTD on gold nano structures where the geometrical singularities at the edges of the square generate a high gradient of the electric near field. A good agreement between the different methods was observed.
Some numerical benchmarks results computed using various commercial and other software have been published by IEEE Germany [130]. For example problem 1 is a sphere with complex valued permittivity in a field of a plane wave  [132]. DDA demonstrated to be faster for smaller refractive indices, while the FDTD poved to be more efficient for larger values.

Conclusion
In this paper an overview of the progress in developing light scattering programs suitable for application in optical particle characterization has been given. A fundamental description of each technique has been presented and suitable references have been provided such that the reader can find more detailed information on the various methods and sources of computer programs.
The state-of-the-art in numerical light scattering modeling is progressing rapidly especially with fast advancing research fields such as nanophotonics and nearfield optics. With almost all concepts there was much progress in the recent years.
The recent advance in computer hardware, and the development of fast algorithms with reduced computational demands and memory requirements, have made an exact numerical solution of the problem of scattering from large scat-