Numerical solutions of the macroscopic Maxwell equations for scattering by non-spherical particles: A tutorial review

Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Contents lists available at ScienceDirect

Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt

Numerical solutions of the macroscopic Maxwell equations for scattering by non-spherical particles: A tutorial review Michael Kahnert a,b,n a b

Research Department, Swedish Meteorological and Hydrological Institute, Folkborgsvägen 17, SE-601 76 Norrköping, Sweden Department of Earth and Space Science, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden

a r t i c l e i n f o

abstract

Article history: Received 22 September 2015 Received in revised form 29 October 2015 Accepted 29 October 2015

Numerical solution methods for electromagnetic scattering by non-spherical particles comprise a variety of different techniques, which can be traced back to different assumptions and solution strategies applied to the macroscopic Maxwell equations. One can distinguish between time- and frequency-domain methods; further, one can divide numerical techniques into finite-difference methods (which are based on approximating the differential operators), separation-of-variables methods (which are based on expanding the solution in a complete set of functions, thus approximating the fields), and volume integral-equation methods (which are usually solved by discretisation of the target volume and invoking the long-wave approximation in each volume cell). While existing reviews of the topic often tend to have a target audience of program developers and expert users, this tutorial review is intended to accommodate the needs of practitioners as well as novices to the field. The required conciseness is achieved by limiting the presentation to a selection of illustrative methods, and by omitting many technical details that are not essential at a first exposure to the subject. On the other hand, the theoretical basis of numerical methods is explained with little compromises in mathematical rigour; the rationale is that a good grasp of numerical light scattering methods is best achieved by understanding their foundation in Maxwell's theory. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Scattering Numerical methods Nonspherical particles

Contents 1. 2.

3.

n

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Maxwell's equations and the electromagnetic scattering problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1. Microscopic Maxwell equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1. Applications of the microscopic Maxwell equations1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. Macroscopic Maxwell equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3. The link between electromagnetic theory and classical mechanics2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4. Constitutive relations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5. Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6. The frequency-domain macroscopic Maxwell equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical methods for solving the macroscopic Maxwell equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Finite-difference time-domain method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 3 3 4 4 5 5 5 6 7 7

Correspondence address: Research Department, Swedish Meteorological and Hydrological Institute, Folkborgsvägen 17, SE-601 76 Norrköping, Sweden. E-mail address: [email protected]

http://dx.doi.org/10.1016/j.jqsrt.2015.10.029 0022-4073/& 2015 Elsevier Ltd. All rights reserved.

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

2

4. 5.

3.2. The volume integral equation and the discrete dipole approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.3. Separation of variables and point-matching method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.4. The T-matrix formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.5. Waterman's T-matrix method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Comparison of different methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Further reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

1. Introduction Electromagnetic scattering by particles is described by Maxwell's theory of electromagnetism. For spherical particles the scattering problem has been solved over a century ago; this is commonly known as the Lorenz–Mie theory [1]. However, particles in nature often have nonspherical shapes. The optical properties of nonspherical particles usually differ considerably from spherical particles of comparable sizes and chemical composition. For non-spherical particles the electromagnetic scattering problem can only be solved by numerical techniques. The use of such techniques generally places much higher demands on the user than the use of Lorenz–Mie programs. This may have contributed to an overuse of spherical model particles far beyond their narrow range of applicability. This review aims at providing an introduction to the required background knowledge on numerical methods for scattering by nonspherical particles. Numerical methods for solving Maxwell's equations find numerous applications in physics and engineering, e.g. in atmospheric optics, ocean optics, remote sensing, astronomy, biomedical optics, nano- and near-field optics, process engineering, and combustion diagnostics. Opensource computer codes based on various numerical methods are freely available for general download on the internet [2]. Quite many codes are deceptively easy to use, which could mislead the user into thinking that numerical light scattering codes can be deployed as black-box tools that require little understanding of numerical methods. On the other hand, available reviews and monographs on methods for solving Maxwell's equations are often written for seasoned users. As a consequence, they tend to cover a lot more breadth and depth than practitioners and newcomers need to acquire at their first exposure to the subject. It is the aim of this paper to fill a gap in the existing literature by providing a brief introduction to numerical solvers of Maxwell's equations that will be tailored to the needs of novices, students, and practitioners who mainly want to become competent users and learn the main features of numerical light scattering methods, but who do not require to fathom all the technical details of numerical methods. Thus this paper mainly pursues pedagogical aims. However, it is not intended to achieve these aims by lowering the standards of mathematical rigour; the central theme is the development of numerical methods from first principles. Instead, the review will maintain a tight focus on a selected number of key methods. The selection is, on one hand, motivated by the pedagogical need of illustrating different key ideas for solving Maxwell's equations.

Thus, only one method has been selected to illustrate volume integral equation methods, only one finite difference method, etc. On the other hand, the selection is also guided by the popularity of different numerical methods. Fig. 1 shows the relative frequency of numerical methods presented at the Electromagnetic and Light Scattering conference series during 1998–2014. The statistics clearly shows that T-matrix methods and the discrete dipole approximation (DDA) tend to be dominant in the electromagnetic scattering community. However, there is no trend indicating that these two methods are entirely displacing other methods. Quite on the contrary, there are continuing efforts to use, develop, and invent alternative methods, each with their own advantages and niches. Among those other methods that will be omitted in the present review are, e.g., the finite element and boundary element method, the method of lines, the method of moments, Fredholm integral equation methods, or the multiple multipole method. In this review the finite-difference time domain (FDTD) method will be presented as a typical exponent of finitedifference methods. This method also has a great pedagogical value, since it is the most direct method for solving the macroscopic Maxwell equations. The DDA will serve to illustrate volume integral equation methods. The separation of variables method (SVM) will be presented here as a general approach that forms the basis for other approaches, such as the point-matching method and the T-matrix approach. Waterman's extended boundary condition method will be employed to illustrate T-matrix methods. Before proceeding to a discussion of these numerical methods in Section 3, a brief review of Maxwell's equations, boundary conditions, and constitutive relations is Frequency of numerical methods in ELS presentations 1 0.8

DDA FDTD SVM other

0.6 0.4 0.2 0

III

V

VI

VII VIII IX

X

XII XIII XIV

Fig. 1. Relative frequency of numerical methods in papers presented at ELS conferences during 1998–2014. (Data for the fourth and eleventh ELS conference were not available.).

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

3

Fig. 2. Figurative representation of the electromagnetic scattering process: incident field Einc in the absence of any obstacles (left), total field Etot in the presence of an obstacle (centre), and scattered field Esca ¼ Etot  Einc (right).

given in Section 2. Some remarks on comparisons of different methods can be found in Section 4. In view of the pedagogical intent of this paper, literature citations in Sections 2 and 3 have been limited to historical key papers, review articles, and monographs. Suggestions for further reading are given in Section 5.

2. Maxwell's equations and the electromagnetic scattering problem The process of electromagnetic scattering is illustrated in Fig. 2. We consider an electromagnetic field in the absence of any obstructions (left panel), which we call the incident field Einc . When an obstacle is introduced the field is modified (centre panel); the altered field is called the total field Etot . Electromagnetic scattering is the process by which the introduction of the obstacle alters the field. We formally define the scattered field as the difference Esca ¼ Etot  Einc (right panel). We have to discriminate between the direct and the inverse problem. In the direct problem one determines the field Etot (or, equivalently, Esca ) for a given obstacle and incident field. This will be the subject of the present review. In the inverse problem one wishes to obtain information about the obstacle from observations of Etot , and from knowledge of the incident field. Techniques for solving inverse problems very often rely on the ability of solving the direct problem. The direct problem is, therefore, of fundamental importance even in inverse modelling, retrieval theory, and data assimilation (e.g. [3]). The theoretical framework for solving the direct electromagnetic scattering problem is provided by Maxwell's theory of electromagnetism. Maxwell's equations can be formulated in two equivalent forms, the microscopic and the macroscopic Maxwell equations. The main difference between the two forms is how the effect of the electromagnetic field on matter is incorporated into the equations. In the microscopic form, it is expressed as part of the source terms, i.e., the charge and current density. In the macroscopic form, it is expressed as a macroscopic polarisation and magnetisation that enter into the displacement field D and the magnetic field H. A good grasp of Maxwell's equations, especially in its

macroscopic form, is essential for understanding numerical methods for solving light scattering problems. 2.1. Microscopic Maxwell equations The four Maxwell equations in its microscopic form are ϵ0 ∇  Eðr; tÞ ¼ ρðr; tÞ

ð1Þ

∇  Bðr; tÞ ¼ 0;

ð2Þ

1 ∂Eðr; tÞ ∇  Bðr; t Þ  ϵ0 ¼ J ðr; t Þ; μ0 ∂t

ð3Þ

∇  E ðr; t Þ þ

∂Bðr; tÞ ¼0 ∂t

ð4Þ

which are known, respectively, as Gauss's law for the electric field, Gauss's law for the magnetic induction, the Ampère–Maxwell law, and Faraday's induction law. The equations are given in SI-units. Here, ϵ0 and μ0 denote the permittivity and permeability in vacuum, respectively. The sources of the electromagnetic field are the total charge density ρ and the total current density J. Note that the original (magnetostatic) form of Ampère's law is ð1=μ0 Þ∇  B ¼J. However, taking the divergence would yield ∇  J ¼ 0, which is not generally consistent with charge conservation. This was first noted by Maxwell, who introduced the displacement-term  ϵ0 ∂E=∂t into Ampère's law to rectify this problem. Taking the divergence of Eq. (3) yields, in conjunction with Eq. (1) ∇  Jþ

∂ρ ¼ 0: ∂t

ð5Þ

Thus, with Maxwell's displacement term, charge conservation is implicit in the theory. Eqs. (1)–(4) are sometimes referred to as “Maxwell's equations in vacuum”. This terminology may cause confusions, as the charges and currents appearing on the right hand side of the equations are necessarily connected to material particles, such as electrons and protons. Thus, these equations are, quite obviously, not limited to empty space. A term that is, arguably, more suitable is “microscopic Maxwell equations”, hereafter abbreviated by “μME”. This term alludes to the fact that, when interacting with matter, the electromagnetic field induces charge displacements and

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

currents on the microscopic level, which need to be included in the total charge density ρ and the total current density J. 2.1.1. Applications of the microscopic Maxwell equations1 The μME are often the starting point for theoretical investigations. Here I just mention a few examples for the usefulness of these equations without going into the details. Readers who only are interested in practical computations can skip this paragraph and jump directly to Section 2.2. The most important uses of the μME in theoretical physics are in relativistic theories. For instance, one may want to express Maxwell's equations in a form that clearly shows the relativistic nature of the theory, and that can be readily extended to a general relativistic theory of electromagnetism. Or one may want to express the theory in the variational framework of Euler–Lagrange theory. Owing to the fundamental importance of Noether's theorem, this is particularly useful when considering symmetries and gauge invariance, as in the formulation of a quantum field theory of electromagnetism. To see the relativistic nature of the theory, one needs to apply the Lorentz transformation properties of the electromagnetic field. If we demand that the non-existence of magnetic monopoles, expressed in Eq. (2), be valid in all Lorentz frames, then application of the transformation properties of the magnetic induction to Eq. (2) leads to Faraday's law, Eq. (4). Thus Faraday's law follows from Gauss's law for the magnetic induction in conjunction with a symmetry principle, namely, Lorentz-covariance of the laws of nature. In a similar way, one can show that Eqs. (1) and (3) are interlinked by Lorentz covariance. This suggests that the two homogeneous microscopic Maxwell equations as well as the two inhomogeneous microscopic Maxwell equations each should be combined into a single, manifestly Lorentz covariant equation. This can, indeed, be achieved within the language of differential forms applied in 4dimensional space-time. One combines the electric field and magnetic induction into the Faraday 2-form F , one combines the charge and current density into a current 3form J , and one introduces the Maxwell 2-form M ¼ ⋆F , where ⋆ denotes the Hodge operator. Then the microscopic Maxwell equations take the form dF ¼ 0

ð6Þ

dM ¼ J ;

ð7Þ

where “d” denotes the exterior derivative of differential forms. This is, without any doubt, the most elegant, concise, and general formulation of the μME. Both equations are independent of the choice of coordinate system, and they are independent of the metric; hence they are readily extended to a general relativistic theory of electromagnetism. The only place where the metric enters is in the relation M ¼ ⋆F , because the Hodge operator implicitly depends on the metric. More information on Maxwell's equations in the framework of quantum field theory and general relativity can be found in [4] and [5], 1 This section provides additional information that lies somewhat outside the main thread; it can be omitted at a first reading.

respectively. A mathematically stringent proof for the equivalence of Eqs. (6) and (7) with the μME is presented in [6]. A less stringent but more intuitive introduction to μME in the language of differential forms is given in [7]. The μME are most often employed in such theoretical developments (although the equivalent macroscopic Maxwell equations (MME) could be used just as well). In practical applications it can be a formidable task to actually account for the microscopic effects of the fields on individual charges in matter, and to include these effects in the source terms. Even after averaging charge distributions over a characteristic length scale (e.g. 10  10–10  8 m in solids and liquids), the microscopic formulation of Maxwell's equations still proves to be impractical. Such difficulties are avoided by employing the MME, which are discussed in the next section. 2.2. Macroscopic Maxwell equations For practical applications it is a great advantage to include the effect of the fields on matter in the fields themselves, rather than in the source terms. To this end, one can exploit the fact that the sum of all induced electric dipole moments on a microscopic level amounts to a macroscopic polarisation field P. Analogously, induced currents and their related magnetic moments on a microscopic level add up to a macroscopic magnetisation field M. We can formally split the total charge density ρ ¼ ρb þρf into a “bound” charge density ρb that gives rise to the macroscopic polarisation, and a “free” or excess charge density ρf . The relation between the bound charge density and the macroscopic polarisation is given by ρb ¼ ∇  P. We can substitute these relations into Eq. (1) to obtain ∇  ðϵ0 E þ PÞ ¼ ρf :

ð8Þ

Similarly, we split up J ¼ J f þ J b , use J b ¼∇  M þ ∂P=∂t, substitute into Eq. (3), and obtain
 1 ∂ B  M  ð ϵ0 E þ P Þ ¼ J f : ð9Þ ∇ μ0 ∂t We now define the electric displacement D¼ϵ0 E þ P and the magnetic field H ¼ð1=μ0 ÞB  M. We thus obtain Maxwell's equations in the following form: ∇  Dðr; tÞ ¼ ρf ðr; tÞ

ð10Þ

∇  Bðr; tÞ ¼ 0

ð11Þ

∂Dðr; tÞ ¼ J f ðr; t Þ ∂t

ð12Þ

∂Bðr; tÞ ¼ 0: ∂t

ð13Þ

∇  H ðr; t Þ  ∇  Eðr; t Þ þ

These equations are known as the macroscopic Maxwell equations. Although they are mathematically equivalent to Eqs. (1)–(4), these equations are better suited for practical computations, such as solving light scattering problems. The induced charge displacements and currents in matter are no longer included in the source terms on the right hand side; rather, their effects are included in the macroscopic fields D and H on the left hand side. In many

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

applications one deals with non-conducting, electrically neutral scatterers, i.e., J f ¼ 0 and ρf ¼ 0. In such cases the advantages of the MME become particularly obvious; when there are no free currents nor free charges, the MME become homogeneous. Eqs. (10)–(13) will be the foundation of our discussion of numerical light scattering methods in the following sections. Note that taking the divergence of Faraday's law (13) yields ∂∇  B=∂t ¼ 0, hence ∇  B ¼ const. Gauss's law for the magnetic induction (11) tells us that this constant is zero. Thus Eq. (11) can be interpreted as an initial condition for Faraday's law. Similarly, taking the divergence of Ampère's law (12) yields  ∂∇  D=∂t ¼ ∇  J f . Consistency with charge conservation (5) requires that ∇  J f ¼ ∂ρf =∂t, hence ∇  D ρf ¼ const. Again, Gauss's law for the displacement field (10) tells us that this constant is zero; it can therefore be viewed as an initial condition for Ampère's law. 2.3. The link between electromagnetic theory and classical mechanics2 Although not immediately relevant for the rest of this review, I mention here, for the sake of completeness, the equation that provides the link between electromagnetic theory and classical mechanics. While electromagnetism deals with electric and magnetic fields, classical mechanics is concerned with forces that act on bodies, as well as with the effects these forces have on the motion of a body. The force F exerted by the electromagnetic field on a body carrying a charge q and moving with a velocity v is given by the Lorentz force equation F ¼ qðE þ v  BÞ:

ð14Þ

2.4. Constitutive relations Spelled out in components, Eqs. (10)–(13) provide us with eight equations for 12 unknown field components. Our theory is, as yet, underconstrained. To guarantee the uniqueness of the solution, we need additional constraints in the form of constitutive relations among the fields E, D, B, H, and J f . The form of the constitutive relations depends on the medium through which the electromagnetic field propagates. To pick a realistic case, let us consider materials that are linear and isotropic, yet inhomogeneous and temporally dispersive. In that case the constitutive relations have the form Z t dτϵðr; t  τÞEðr; τÞ ð15Þ Dðr; tÞ ¼ 1

Z Bðr; tÞ ¼

t 1

Z J f ðr; tÞ ¼

t 1

dτμðr; t τÞHðr; τÞ

ð16Þ

dτσðr; t  τÞEðr; τÞ:

ð17Þ

2 This section provides additional information that lies somewhat outside the main thread; it can be omitted at a first reading.

5

The conductivity σ is large for conducting materials, while for non-conducting materials one can set J f ¼ 0. The constitutive parameters are often expressed in terms of the relative permittivity ϵr and the relative permeability μr according to ϵ ¼ ϵ0 ϵr

ð18Þ

μ ¼ μ0 μr :

ð19Þ

A standard example that will be used throughout this review is the case of perfect dielectric materials, which are characterised by ϵr a 1, μr ¼ 1, and σ is very small, so that J f ¼ 0. Eqs. (15)–(17) express the fact that in temporally dispersive materials there is a time lag between the fields and the resulting polarisation, magnetisation, and current density. Pulsed waves tend to disperse, i.e., spread out in such materials; thus the term “dispersive”. Note that the integrals in Eqs. (15)–(17) only extend up to t. Thus the quantities on the left hand side do depend on the fields at earlier, but not at later times, which is required by the causality principle. The time lag between the fields and the material response becomes negligible only in the low-frequency limit; in that case Eqs. (15)–(17) simplify to Dðr; tÞ ¼ ϵðrÞEðr; tÞ

ð20Þ

Bðr; tÞ ¼ μðrÞHðr; tÞ

ð21Þ

J f ðr; tÞ ¼ σðrÞEðr; tÞ:

ð22Þ

Eqs. (20) and (21) also apply to empty space, with ϵ ¼ ϵ0 and μ ¼ μ0 . A more detailed discussion of constitutive relations for different types of materials, including anisotropic, biisotropic, and nonlinear media can be found in [8]. There one can also find a concise discussion of the conditions that the constitutive relations need to fulfil in order to guarantee that the MME become well-posed. Further, it is explained how the constitutive relations are constrained by the requirement of Lorentz-covariance. 2.5. Boundary conditions Maxwell's equations provide us with a theory of the behaviour of electromagnetic fields propagating through continuous media, where all fields, sources, and constitutive parameters are continuously differentiable functions of space and time. However, at the interfaces between different media the constitutive relations or the source terms may make discontinuous jumps. In such case one needs additional conditions that relate the fields in the two spatial domains across the adjoining surface. Following [9], I take the point of view that the boundary conditions are postulated rather than derivable. Thus they are stated here as facts based on our experimental knowledge. For a discussion of how to motivate the boundary conditions by use of Maxwell's integral equations the reader is referred to [8]. However, as pointed out in [9], it is important to understand that such an approach does by no means constitute a rigorous “proof” of the boundary conditions. The Maxwell integral equations are derived from

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

the corresponding differential equations. A sufficient condition for the existence of the integrals is the continuity of the integrands. However, without any prior knowledge of the boundary conditions we do not know if the integrands are continuous across the boundary, hence if the integrals even exist. It is, then, not entirely consistent to use the integral form of Maxwell's equations for “deriving” the boundary conditions. Let n^ denote the surface normal to the boundary surface in question (where n^ will, in general, depend on both position and time). Then the boundary conditions for the tangential components of the electric and magnetic fields are n^  ðE 2  E1 Þ ¼ 0

ð23Þ

n^  ðH 2  H 1 Þ ¼ J s

ð24Þ

where the subscripts 1 and 2 label the fields in the two spatial regions on either side of the boundary, and where J s represents a surface current on the boundary. For a nonconducting surface Eq. (24) becomes homogeneous, i.e., the tangential component of H is continuous across the adjoining surface. There are also boundary conditions that state that the normal component of the B-field is continuous, while the normal component of the D-field makes a jump equal to the surface charge density. However, only the two conditions given in Eqs. (23) and (24) are required for the uniqueness of the solution. 2.6. The frequency-domain macroscopic Maxwell equations One frequently invokes the assumption that all timedependencies are harmonic. In such case one can represent, e.g., the electric field as   Eðr; tÞ ¼ Re EðrÞ expð  iωtÞ ð25Þ with a complex field amplitude EðrÞ. Similar representations are assumed for D, B, H, J f , and ρf . Note that the notational distinction between the time-dependent real quantities represented by italic letters, and the complex time-independent quantities represented by upright letters. Substitution of Eq. (25) and analogous relations for the other fields and source terms into Eqs. (10)–(13) yields ∇  DðrÞ ¼ ρf ðrÞ

ð26Þ

∇  BðrÞ ¼ 0

ð27Þ

∇  HðrÞ þiωDðrÞ ¼ Jf ðrÞ

ð28Þ

∇  EðrÞ  iωBðrÞ ¼ 0:

ð29Þ

These are the frequency-domain MME. The assumption of a harmonic time-dependency is considerably more general than it may seem at a first glance. Firstly, nonmonochromatic fields can be represented by a Fouriersuperposition of monochromatic, i.e., time-harmonic fields. Secondly, even non-coherent fields can be described by Eqs. (26)–(29), provided that the coherence time, i.e., the time scale on which the random phase fluctuations occur, is much longer than the oscillation period T ¼ 2π=ω of the fields.

The frequency-domain version of the continuity equation (5) becomes ∇  Jf  iωρf ¼ 0:

ð30Þ

To work in the frequency-domain has a number of advantages over the time-domain. For instance, if we substitute Eq. (25), and analogous relation for the other fields, into Eqs. (15)–(17), then the convolution integrals in the constitutive relations become simple products in the frequency domain: DðrÞ ¼ ϵðr; ωÞEðrÞ

ð31Þ

BðrÞ ¼ μðr; ωÞHðrÞ

ð32Þ

Jf ðrÞ ¼ σðr; ωÞEðrÞ;

ð33Þ

where Z ϵðr; ωÞ ¼

1

0

dt ϵðr; tÞ expðiωtÞ;

ð34Þ

and similarly for μ and σ. Thus the time-dispersive realvalued permittivity ϵ, permeability μ, and conductivity σ become frequency-dispersive, complex valued quantities ϵ, μ, and σ, respectively. (Note, again, the notational distinction between italic and upright letters!) Hence, time dispersion and frequency dispersion are two sides of the same coin. Another characteristic of the frequency domain is the following. Since we no longer deal with a time-dependent problem (all time-derivatives have been eliminated), Gauss's laws for the displacement field and for the magnetic induction no longer provide any initial conditions for solving the two Maxwell-curl equations. This can be seen by taking the divergence of Eq. (29), which directly yields Eq. (27). Recall that in the time domain we had obtained ∇  B ¼ const. Thus, Eq. (27) no longer provides any extra information; it is now implicit in Eq. (29). Similarly, taking the divergence of Eq. (28) in conjunction with Eq. (30) yields ∇  D ρf ¼ 0, while in the time-domain we had obtained ∇  D ρf ¼ const. Hence in the frequency domain Eq. (26) no longer provides us with an initial condition for Ampère's law. (However, we can only “derive” Eq. (26) from Eq. (28) if we assume Eq. (30). Alternatively, we need Eq. (26) in conjunction with Eq. (28) to derive charge conservation, Eq. (30).) Eqs. (31) and (33) inserted into Eq. (28) yield ∇  H þ iωϵE ¼ σE. Thus, with the definition σðr; ωÞ εðr; ωÞ ¼ ϵðr; ωÞ þ i ω

ð35Þ

one obtains ∇  HðrÞ þiωεðr; ωÞEðrÞ ¼ 0:

ð36Þ

For non-conducting materials, one simply has ε ¼ ϵ. Similarly, Eq. (29) in conjunction with Eq. (32) yields ∇  EðrÞ iωμðr; ωÞHðrÞ ¼ 0:

ð37Þ

Eqs. (36) and (37) will be the starting points for developing frequency-domain numerical methods for frequencydispersive media.

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

7

3. Numerical methods for solving the macroscopic Maxwell equations

and likewise for Hðr; tÞ. Substituting this into Eq. (40) and making the substitution t τ ¼ t 0 yields

3.1. Finite-difference time-domain method

∂ ∇  HðrÞ expð  iωt Þ EðrÞ expð  iωt Þ ∂t

For didactic reasons, I start this overview with the finite-difference time-domain method (FDTD), because its main concepts are quite easy to grasp. The method has been originally proposed by Yee [10]. The idea is to start from the MME in the time domain, Eqs. (10)–(13), and to discretise both the spatial and the temporal coordinates. However, we should be careful to clearly state which assumptions we make in the formulation of the method. We assume that there are no free charges (ρf ¼ 0), and that the material is a perfect dielectric, i.e., it is nonconducting (J f ¼ 0) and non-magnetic (μ ¼ μ0 ). Besides that, it would be tempting to assume the constitutive relation given in (20); then Eqs. (12) and (13) would take the rather simple form ∂Eðr; tÞ ∇  H ðr; t Þ  ϵðrÞ ¼0 ∂t ∇  Eðr; t Þ þμ0

∂Hðr; tÞ ¼ 0: ∂t

ð38Þ ð39Þ

This would be straightforward to discretise. However, as noted in Section 2.4, the constitutive relation given in Eq. (20) is only valid in the low-frequency limit. This is too restrictive for many important applications. On the other hand, if we employ the more general time-dispersive constitutive relation (15), then Ampère's law takes the form of an integro-differential equation Z ∂ t ∇  H ðr; t Þ  dτ ϵðr; t  τÞEðr; τÞ ¼ 0: ð40Þ ∂t  1 This is non-local in time and could be cumbersome to integrate by numerical time-iteration. We wish to find a way to formulate Ampère's law in the time domain in such a way that it is local in time, but without making too restrictive assumptions about the constitutive relations. This can be achieved by restricting the time-dependency of the fields. We assume that the fields be time-harmonic, i.e. Eðr; tÞ ¼ EðrÞ expð iωtÞ;

ð41Þ

Z

1 0

dt 0 ϵðr; t 0 Þ expðiωt 0 Þ ¼ 0:

ð42Þ Thus, in conjunction with Eq. (34), we obtain ∇  HðrÞ expð  iωtÞ þiωϵðr; ωÞEðrÞ expð iωtÞ ¼ 0:

ð43Þ

We want to return now to the time domain. To this end, we explicitly write the complex permittivity in terms of its real and imaginary parts, ϵ ¼ ϵR þ iϵI , and use once more iωEðrÞexpð iωtÞ ¼ ð∂=∂tÞEðrÞexpð  iωtÞ. This yields ∇  H ðr; t Þ  ωϵI ðr; ωÞEðr; t Þ ϵR ðr; ωÞ

∂Eðr; tÞ ¼ 0: ∂t

ð44Þ

This form of Ampère's law in conjunction with Faraday's law given in Eq. (39) can be used as starting points for the FDTD for perfect dielectric, dispersive materials. The price we pay for avoiding the difficulties we encountered in Eq. (40) is that we have limited ourselves to time-harmonic fields. Thus, if we want to apply Eqs. (39) and (44) to more general cases (e.g. laser pulses as incident fields), then we would need to solve these equations for a number of frequencies and obtain the general solution by performing a Fourier superposition. The discretisation of the equations in the FDTD is straightforward. First one fixes a computational domain that contains the target particle, and one selects a coordinate system, e.g., Cartesian coordinates x; y; z. Next one discretises the coordinates by specifying a grid in the chosen coordinate system, e.g. ðx; y; zÞ ¼ ðiΔx; jΔy; kΔzÞ. This is illustrated in Fig. 3 (centre). Similarly, one discretises time, t ¼ nΔt. Thus all field components are characterised by specifying their values at the discrete spacetime points, e.g. El ðiΔx; jΔy; kΔz; nΔtÞ ¼ El ði; j; k; nÞ;

ð45Þ

where El denotes the lth component of the field vector E, and where i; j; k; n are cardinal numbers. All derivatives are numerically approximated by finite difference quotients, e.g. ∂El ðx; y; z; tÞ El ði; j; k; n þ 1=2Þ El ði; j; k; n  1=2Þ  ; ∂t Δt

ð46Þ

Fig. 3. Scattering particle (left); discretised computational domain in finite-difference methods, such as FDTD (centre); discretised particle volume in volume-integral equation methods, such as DDA (right).

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8

and similarly for the spatial derivatives. For instance, the x-component of Eq. (39) becomes     H x i þ 1=2; j; k; n þ1=2 ¼ H x iþ 1=2; j; k; n 1=2 1 1 Ez ðiþ 1=2; j þ 1=2; k; nÞ  Δt μ0 Δy

Ez ði þ 1=2; j 1=2; k; nÞ 1 Ey ði þ 1=2; j; k þ 1=2; nÞ Δz 

Ey ðiþ 1=2; j; k  1=2; nÞ ;



ð47Þ

and similarly for Eq. (44). We see that this provides an algorithm for solving Maxwell's equations by a timeiterative process. The component Hx at time-step ðn þ 1=2Þ is computed in terms of the electric and magnetic fields at earlier time steps. To use the algorithm, one has to provide an initial field. Among the main advantages of the FDTD is its conceptual simplicity; it is based on directly discretising the time-domain MME. Also, it is straightforward to implement the method. However, there are several technical issues that need to be addressed to make this algorithm work in practice. The most important practical problem arises from the fact that the computational domain is finite; thus one has to make sure that there are no spurious reflections at the boundary back into the computational domain. One of the most successful methods to eliminate such reflections is the construction of an absorbing multi-layered medium at the domain boundary, which is known as a uniaxial perfectly matched layer. Other technical questions pertain to the choice of the discrete spatial mesh, the method for computing an average permittivity per grid cell, the required size of Δx, Δy, Δz relative to the wavelength of light, the size of Δt and its relation to the spatial increments, and how to post-process the near-field time-domain results to obtain a far-field frequency-domain result, from which one can compute the single-scattering and absorption optical properties. A recent review of the FDTD and its applications to particles and surfaces can be found in [11]. That review provides a detailed discussion of the theoretical basis and of important technical issues related to the FDTD. Some of the more subtle practical details, including numerical stability problems and the representation of dispersive materials in FDTD computations, are discussed in [12]. 3.2. The volume integral equation and the discrete dipole approximation We now move from the time- to the frequency domain, and we begin our discussion with one of the most general methods for solving the frequency-domain MME. Our starting points are Eqs. (36) and (37). To pick a specific case we consider a perfect dielectric particle, i.e., the magnetic permeability is constant, μ ¼ μ0 , and the conductivity σ is zero, which implies εðr; ωÞ ¼ ϵðr; ωÞ ¼ ϵ0 ϵr ðr; ωÞ — see Eqs. (18) and (35). In such case, if we take the curl of Eq. (37), then the second term will be ∇  ðμðr; ωÞHðrÞ ¼ μ0 ∇  HðrÞ, for which we can substitute

Eq. (36). This yields ∇  ∇  EðrÞ  ω2 μ0 ϵ0 ϵr ðr; ωÞEðrÞ ¼ 0:

ð48Þ

This is the wave equation for the electric field in the frequency domain. It is a homogeneous equation with nonconstant coefficients (on account of the positiondependent factor ϵr ðr; ωÞ). We would like to circumvent the difficulties related to solving differential equations with non-constant coefficients. The trick is to trade this problem for an inhomogeneous equation with constant coefficients. This is achieved as follows. Let us formally introduce jðrÞ ¼  iωϵ0 μ0 ½ϵr ðrÞ  1EðrÞ; ð49Þ pffiffiffiffiffiffiffiffiffiffi as well as k0 ¼ ω ϵ0 μ0 . Then Eq. (48) can be written in the form 2

∇  ∇  EðrÞ  k0 EðrÞ ¼ iωjðrÞ:

ð50Þ

Now the left hand side of this equation only has constant coefficients. The inhomogeneity on the right hand side behaves just like a source term. Note that ϵr 1 ¼ 0 in empty space. Thus if we consider scattering by a dielectric particle surrounded by empty space then the source term is non-zero only inside the particle. Linear inhomogeneous equations can be solved by Green's function method. Green's function G is a solution to the inhomogeneous equation with a point source, i.e. ∇  ∇  Gðr; r0 Þ k0 Gðr; r0 Þ ¼ 1δðr r0 Þ; 2

ð51Þ

where the operator ∇ acts upon the first argument r. The solution to this equation is given by " # 1 expðik0 jr  r0 jÞ 0 ð52Þ Gðr; r Þ ¼ 1 þ 2 ∇∇ 4πjr r0 j k 0

(see, e.g., [13]). A general source jðrÞ can be thought of as being composed of a continuous distribution of point sources. Thus, owing to the linearity of the wave equation, once we know the solution of Eq. (51) for a point source, the solution of Eq. (50) for a general source can be obtained by R 3 linear superposition, i.e. EðrÞ ¼ iω V int d r 0 Gðr ; r0 Þ  jðr0 Þ. The integration only extends over the volume V int of the scattering target, since in the surrounding medium the source term j vanishes. We can add to this a solution to the corresponding homogeneous equation, which is given by the incident field Einc . Thus the general solution to the inhomogeneous wave Eq. (50) is given by Z 3 d r 0 Gðr; r0 Þ  jðr0 Þ; ð53Þ EðrÞ ¼ Einc ðrÞ þiω V int

where r A R3 , and the integration variable r0 varies over all points inside the particle volume. (To double-check the correctness of this equation, one can apply the operator 2 ð∇  ∇  1  k0 Þ, use the fact that Einc solves the homogeneous wave equation, and make use of Eqs. (51) and (50)). Eq. (53) is known as the volume-integral equation (VIE). As an illustration of numerical VIE methods, let us consider the discrete dipole approximation (DDA), originally proposed in [14]. The main idea is to discretise the particle volume into M small volume cells V 1 ; …; V M . This can be done in different ways, but one of the most

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

common ones is to use a cubical grid. This is illustrated in Fig. 3 (right). Since the integration domain in the VIE only extends over the inside of the particle, the computations will not involve any points in the surrounding medium, as was the case in the FDTD (see Fig. 3, centre). By discretising the particle volume the VIE will be transformed into an algebraic system of linear equations, which can be solved with standard numerical techniques. It is instructive to understand, at least, the broad outline of this derivation. We are interested in the field inside volume cell m; Eq. (53) can be written as Z 3 EðrÞ  iω d r 0 Gðr; r0 Þ  jðr0 Þ Vm

¼ Einc ðrÞ þ iω

M Z X k ¼ 1 kam

Vk

d r 0 Gðr; r0 Þ  jðr0 Þ; 3

ð54Þ

where we now consider r A V m . The reason for isolating the integral over Vm from the rest of the volume integral is because for r ¼ r0 Green's function becomes singular. This singularity needs to be treated with care. A more detailed discussion of this problem can be found in [15]. The expression on the right hand side is known as the exciting field Eexc , because it can be interpreted as the field that is exciting volume cell m. It is composed of the incident field plus the field contributions originating from all other volume cells inside the target. Up to this point we have not made any approximations. We now invoke the assumption that the wavelength of the field is large compared to the size of the volume cells, so that the amplitude variation of the fields over each volume cell can be neglected. In other words, in each cell k the field amplitude EðrÞ is assumed to have a constant value Ek inc (and similarly for Eexc k , Ek , jk , etc.). Physically, this means that all charges inside a volume cell oscillate in phase, just like a single dipole. Thus we effectively describe the target as an array of dipoles with a dipole spacing equal to the dimension of the grid cells (thus the term “discrete dipole approximation”). With the long-wave (or discrete dipole) approximation Eq. (54) can now be written as ^ Eexc m ¼ Em  iωV m G m;m  jm inc Eexc m ¼ Em þ iω

M X

V k Gm;k  jk :

ð55Þ ð56Þ

k ¼ 1 kam

(The caret on top of G^ m;m indicates that, owing to the singularity of Green's function, this term is not simply given by Gðrm ; rm Þ; as mentioned before, this integral needs some care, but it can be shown to be regular.) The discretised form of Eq. (49) now reads jm ¼  iωϵ0 μ0 ½ϵr m 1Em :

ð57Þ

Substitution of the last equation into Eq. (55) yields 2 ^ Eexc m ¼ Am;m  Em , where Am;m ¼ 1 ω μ0 ϵ0 ðϵr m 1ÞV m G m;m . 1 Thus Eq. (57) yields jm ¼  iωϵ0 μ0 ðϵr m  1ÞAm;m  Eexc m . Substitution of this equation into Eq. (56) leads to inc 2 Eexc m ¼ Em þ ω ϵ0 μ0

M X k ¼ 1 kam

1 V k Gm;k  ðϵr k 1ÞAk;k  Eexc k ;

ð58Þ

9

or Einc m ¼

M X

Q m;k  Eexc k ;

ð59Þ

k¼1

where 1 Q m;k ¼ 1δm;k  k0 ð1  δm;k ÞV k Gm;k  ðϵr k  1ÞAk;k : 2

ð60Þ

Eq. (59) is a system of linear equations that can be solved by standard numerical techniques for the unknown exciting field in terms of the known incident field. Clearly, the solution of this system of equations depends on the incident field. Thus, if one wants to solve light scattering problems for different incident fields (e.g. fields incident from different directions or, equivalently, for particles in different orientations relative to the incident field), then the solution procedure needs to be repeated for each new incident field. Alternatively, one needs to solve the system of equations simultaneously for many incident fields by increasing the size of the system of equations. In either case this results in an increase of computation time. What we eventually want is the field outside the target in the surrounding medium. To this end we must use the VIE (53) with r lying in the outside medium, and with the source term jðr0 Þ expressed through the exciting field. The current density jm in cell m can be related to the induced dipole moment pm via pm ¼ ði=ωÞV m jm . Also, the induced dipole moment is linearly related to the exciting field, i.e. pm ¼ tm  Eexc m , where the matrix tm is the polarisability matrix in the mth cell. Thus Eq. (53) becomes EðrÞ ¼ Einc ðrÞ þ ω2

M X

Gðr; rm Þ  tm  Eexc m :

ð61Þ

m¼1

To summarise, the DDA amounts to numerically solving the system of linear equations given in Eq. (59) for the in each grid cell k, and to use unknown exciting field Eexc k this in Eq. (61) to compute the field EðrÞ outside the target. In practice one needs to employ a polarisability model that relates the target's dielectric properties given by ϵr m in each volume cell to the corresponding polarisability tensor tm . Note that the numerical algorithm formulated here does not make any restricting assumptions about the particle shape or the variation of the dielectric properties inside the particle. Thus the method can be applied to arbitrary morphologies and to dielectrically inhomogeneous targets. By decreasing the dipole spacing one can increase the numerical accuracy of the method (at the cost of increasing computation time). The algorithm turns out to be highly stable. However, assuming an over-simplified polarisability model may introduce errors. A detailed review of the DDA, its relation to the method of moments, and a discussion of different polarisability models, as well as how such models are related to the treatment of the singularity of the Green's function, can be found in [15]. A recent review with a thorough treatment of the theory and a discussion of important technical issues of the DDA is given in [16].

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

10

3.3. Separation of variables and point-matching method Let us start again from Eqs. (36) and (37), and consider a perfect dielectric target with no excess charges, i.e., μðr; ωÞ ¼ μ0 , σðr; ωÞ ¼ 0, and ρf ðrÞ ¼ 0. In addition, we assume that the particle is homogeneous (or, at least, piecewise homogeneous), so that we can set εðr; ωÞ ¼ ϵ0 ϵr ðωÞ. We then apply the curl operator to each of Eqs. (36) and (37). Since the magnetic permeability is not dependent of position, it can be pulled out of the curl operation, and we arrive, as we saw before, at the wave equation (48), where ϵr now depends only on ω. Analogously, we can apply the curl operator to Eq. (37). Since the permittivity is now assumed to be position-independent, it can be pulled out of the curl, and we can substitute Eq. (36), which yields ∇  ∇  HðrÞ  ω2 μ0 ϵ0 ϵr ðωÞHðrÞ ¼ 0:

ð62Þ

We can further use the identity ∇  ∇  A ¼∇ð∇  AÞ  ∇2 A, which is valid for arbitrary smooth vector fields A. Since we have no excess charges, and since the magnetic and dielectric properties are position-independent, Eqs. (26), (27), (31) and (32) imply ∇E¼0

ð63Þ

∇  H ¼ 0;

ð64Þ

whence ∇  ∇  E ¼  ∇ E,pand likewise for H. With this in ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi conjunction with kðωÞ ¼ ω ϵ0 ϵr ðωÞμ0 , the wave equations become h i 2 ∇2 þ k ðωÞ EðrÞ ¼ 0 ð65Þ 2

h

i 2 ∇2 þ k ðωÞ HðrÞ ¼ 0:

ð66Þ

This is known as the Helmholtz equation for the complex amplitudes of the electric and magnetic fields. Note that in the surrounding medium ϵr ¼ 1, hence k ¼ k0 . The separation of variables method (SVM) is a general approach for solving linear partial differential equations, such as the Helmholtz equation. The key idea is to find a complete set of functions, where each function in the set solves the equation in question. A general solution of the equation can then be written as a linear superposition of the solution functions. In ourh case, one i first considers the scalar Helmholtz 2 equation ∇2 þk ψðrÞ ¼ 0, where ψ is a scalar wave function. Suppose we want to express the position vector in spherical coordinates ðr; θ; ϕÞ. Then we make a separation ansatz, i.e., we assume that the wavefunction can be written as a product of three functions ψðr; θ; ϕÞ ¼ RðrÞΘðθÞΦðϕÞ:

Helmholtz equation is given by ðsÞ ψ n;m ðr; θ; ϕÞ p znðsÞ ðkrÞP ðmÞ n ð cos θÞ expðimϕÞ; n A N0 ; m ¼ n;  n þ 1; …; n;

znð1Þ

znð3Þ

ð1Þ ¼ hn

ð68Þ zð4Þ n

ð2Þ

¼ jn , ¼ nn , ¼ jn þ inn , and ¼ hn where ¼ jn  inn are, respectively, spherical Bessel, Neumann, and are Hankel functions of the first and second kind, P ðmÞ n associated Legendre functions, and expðimϕÞ ¼ cos ðmϕÞ þ i sin ðmϕÞ are complex harmonic functions. The proportionality sign is used here to emphasise that the functions can be multiplied by a suitable normalisation factor. The index n is called the degree, m is called the order of the wavefunction. The superscript s is the kind of the wavefunction. From the scalar wavefunctions one can construct different sets of vector wavefunctions. The most common ones are ðsÞ ðsÞ ðr; θ; ϕÞ p ∇ψ n;m ðr; θ; ϕÞ ψ n;m;0

ð69Þ

 ðsÞ ðsÞ ðr; θ; ϕÞ p ∇  rψ n;m ðr; θ; ϕÞ ψ n;m;1

ð70Þ

 ðsÞ ðr; θ; ϕÞ p ∇  ∇  rψ ðsÞ ψ n;m;2 n;m ðr; θ; ϕÞ :

ð71Þ

The index q ¼ 0; 1; 2 is known as the mode of the wavefunction. It can be shown that the vector spherical wave ðsÞ of all degrees, orders, modes, and kinds functions ψ n;m;q satisfy the vector Helmholtz equation in spherical coordinates. Owing to the linearity of Eq. (65) a general solution of Helmholtz' equation can be written as a linear superposition of vector wave functions. One makes the following ansatz for the incident field Einc , the scattered field Esca , and the internal field Eint inside the target: Einc ðk0 r; θ; ϕÞ ¼

Ncut n 2 X X X n ¼ 0 m ¼ n q ¼ 1

Esca ðk0 r; θ; ϕÞ ¼

Ncut n 2 X X X n ¼ 0 m ¼ n q ¼ 1

Eint ðkr; θ; ϕÞ ¼

Ncut n 2 X X X n ¼ 0 m ¼ n q ¼ 1

ð1Þ an;m;q ψ n;m;q ðk0 r; θ; ϕÞ

ð72Þ

ð3Þ pn;m;q ψ n;m;q ðk0 r; θ; ϕÞ

ð73Þ

cn;m;q ψ ð1Þ n;m;q ðkr; θ; ϕÞ:

ð74Þ

For the magnetic field, one makes an analogous ansatz. Here it is important to note the following facts.

 The sum over the degree n extends, formally, to infinity.

ð67Þ

Substitution into Helmholtz' equation (where one needs to express the Laplacian ∇2 in spherical coordinates) yields three ordinary differential equations, one for the radial, one for the polar, and one for the azimuthal functions. These turn out to be the spherical Bessel, associated Legendre, and harmonic differential equations, respectively. Thus, the set of solution function of the scalar

zð2Þ n



In numerical computations it needs to be truncated at some finite value Ncut that will depend on the required accuracy. Thus, in this approach one approximates the solution function, rather than the differential operator in Helmholtz’ equation. By contrast, the main idea in the FDTD and, indeed, in any finite-difference method is to approximate the differential operator (as in Eq. (46)). In the sum over q the mode q ¼0 has not been included in the ansatz. This is because Eqs. (63) and (64) require the fields to be divergence-free. Since the divergence of a curl is identically zero, it is manifest in Eqs. (69)–(71) that the wavefunctions for the modes q¼1,2 are divergence-free, but the mode q ¼0 is not. (Incidentally, note

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎









that a Fourier transformation of Eq. (63) from position space to wavevector space (i.e., from r to k0 ) would transform Eq. (63) into k0  E ¼ 0, which says that the field is transverse to the direction of propagation. Thus the requirement that, in the absence of excess charges, the fields be divergence-free is equivalent to the requirement that the fields be transverse. Therefore, the mode q ¼0 is also known as the longitudinal mode. It only contributes to the electric field if ρf a 0, in which case the electric field is not divergence-free.) Note that the wavenumber in free space k0 appears in the argument of the incident and scattered fields, while the wavenumber inside the scattering target k is found in the argument of the internal field. Note further that the incident and the internal field are expanded in wavefunctions of the first kind. These are the only kind of wavefunctions that are regular at the origin. Owing to the singularity of the Neumann functions nn, all other kinds of wavefunctions are not regular at the origin. The scattered field is expanded in wavefunctions of the third kind. These kind of wavefunctions represent outgoing wavefunctions that satisfy the boundary condition at infinity, known as the radiation condition, which is given in spherical coordinates by
 r 1 ð3Þ ∇  ψ ð3Þ ðr; θ; ϕÞ ¼ O  ψ n;m;q : ð75Þ n;m;q ðr; θ; ϕÞ  ik0 jrj r This simply says that the wavefunctions converge sufficiently fast to zero at large distances from the target particle. The expansions of the incident, scattered, and internal magnetic fields follow the same pattern as that of the electric field. However, note that these expansions do not contain any new expansion coefficients. Rather, the electric and magnetic fields are related via Eq. (37). It turns out that in the expansion of the magnetic field the role of the modes q ¼1 and q ¼2 are exchanged as compared to the corresponding expansions of the electric fields. Both the FDTD and the DDA require the use of a computational grid — see Fig. 3, centre and right. As a result the geometry of the particle is approximated by a more or less jagged shape. No such approximation of the particle geometry is invoked in the SVM approach.

To summarise, the expansions in Eqs. (72)–(74) together with similar expressions for the magnetic field represent a general ansatz for a divergence-free solution of the vector Helmholtz equation for which the incident and internal fields are free of singularities, and for which the scattered field satisfies the radiation condition at infinity. A detailed explanation of the separation of variables ansatz applied to light scattering problems, as well as a more complete explanation of the properties of vector spherical wave functions can be found in [13]. To solve a specific scattering problem one needs to know the incident field, i.e., one has to specify the coefficients an;m;q in Eq. (72). One then uses the boundary conditions (23) and (24) with J s ¼ 0 for perfect dielectric materials. Expressed in terms of the incident, scattered,

11

and internal fields, the boundary conditions are n^  ðEinc þEsca  Eint Þjboundary ¼ 0

ð76Þ

n^  ðHinc þ Hsca  Hint Þjboundary ¼ 0

ð77Þ

Substitution of the expansion given in Eqs. (72)–(74) yields a relation for the unknown coefficients pn;m;q and cn;m;q of the scattered and internal fields in terms of those of the incident field. One way to formulate a numerical method for non-spherical geometries is to enforce the boundary conditions in a finite number of matching points on the particle surface. If the number of matching points is equal to the number of unknown expansion coefficients, then one obtains a system of linear equations that can be solved by standard numerical techniques. This method is known as the ordinary point-matching method, which was first proposed in [17]. The expansion of the fields in this method often has poor convergence properties for reasons that are discussed in more detail in [18]. A numerically much more stable method is obtained by using more than twice as many matching points as unknowns. The resulting over-determined system of equations can be solved with the least-squares approach. This method is known as the generalised point-matching method [19]. 3.4. The T-matrix formulation The boundary conditions, Eqs. (76) and (77), are linear equations. Thus, substitution of the linear superposition of vector wavefunctions, Eqs. (72)–(74), will yield linear relations among the expansion coefficients, which in obvious matrix-vector notation take the form a¼Q c

ð78Þ

p ¼  RgQ  c

ð79Þ

p ¼ T  a:

ð80Þ

The somewhat unusual notation for the matrix RgQ as well as the minus sign in the second equation are a matter of convention. What we actually want is the T-matrix, which gives us the expansion coefficients p of the scattered field in term of those of the known incident field a. Combining these equations, we obtain T ¼  RgQ  Q  1 :

ð81Þ

Thus, once we know the matrices Q and RgQ , we can compute the T-matrix. The T-matrix provides us with a powerful formulation of the scattering problem. In all numerical methods we encountered thus far one computes for a given particle and incident field the scattered field. For any new incident field (or particle orientation) the procedure needs to be repeated. By contrast, the T-matrix is independent of the direction of the incident field. It only depends on the particle's size, dielectric properties, shape, and on the wavelength of light. It contains the complete information on the particle's spectral scattering and absorption properties. This has several major advantages over other formulations of the scattering problem.

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

12

 Optical properties that characterise scattering and







absorption in the far-field, such as the optical cross sections, the polarimetric differential scattering cross section, the Mueller matrix, etc., can be expressed analytically in terms of the T-matrix. If the T-matrix has been computed for one specific orientation of the target particle, one can compute the T-matrix for other particle orientations (equivalently, for other directions of the incident field) by use of analytic expressions that make use of suitable representations of the rotation group. Closed-form expressions for ensemble-averaged optical properties of randomly oriented particles can be obtained from the T-matrix. These expressions yield a highly efficient and numerically accurate algorithm for computing ensemble-averaged optical properties. Geometrical symmetries of a particle manifest themselves as symmetry relations (known as commutation relations) of the T-matrix. Exploiting these symmetries in numerical algorithms can reduce computation time by several orders of magnitude. This is particularly effective for axisymmetric particles.

The T-matrix formulation was originally introduced by Waterman [20]. Recent reviews of Waterman's approach and the orientational averaging method can be found in [21,22]. In [22] one can also find detailed explanations of how to compute optical properties, such as the Mueller matrix and optical cross sections, from the T-matrix. Green's function perspective on the T-matrix formulation, as well as a review of the use of symmetries in T-matrix computations is given in [13]. 3.5. Waterman's T-matrix method At this point we have little more than a T-matrix formulation of the scattering problem; but we still need an actual method for computing it. It turns out that there are many different methods that can be used for computing a T-matrix. The original method was proposed by Waterman [20]; this approach is also referred to as the extended boundary condition method or the null-field method. Later, the superposition T-matrix method [23] has been developed for solving light scattering problems by clusters of spheres. Also, the SVM in spheroidal coordinates [24], the DDA [25], the point-matching method [26], and the invariant embedding method [27] have been suitably formulated for computing a T-matrix with these methods. Here, Waterman's method will serve as an illustration of a numerical T-matrix method. It has been selected here mainly because of its historical significance, but also because this method is still one of the most widely used T-matrix methods. The derivation of the method starts from the vector Green's identity Z 3 d r 0 ½X  ð∇  ∇  YÞ  Y  ð∇  ∇  XÞ V I dσ 0 n^  ½Y  ∇  X X  ∇  Y ; ¼ ð82Þ ∂V

which is valid for any vector fields X, Y that are sufficiently

well-behaved in a bounded domain V  R3 with boundary ∂V. n^ denotes a unit normal vector on ∂V pointing out of V, and dσ 0 denotes the area element. In our case, we divide three-dimensional space into the region V ext surrounding the scattering target, and the region V int inside the particle. We let V ¼ V ext , so that ∂V will be the union of two sets, the boundary at infinity, ∂R3 , and the boundary of the scattering target, ∂V int . We further substitute X ¼ E and Y ¼ G. The double curl terms simplify owing to the fact that E and G satisfy wave equations, ∇  ∇  E k0 E ¼ 0

2

ð83Þ

∇  ∇  G  k0 G ¼ 1δðr0  rÞ:

ð84Þ

2

This yields Z h i 3 2 2 d r 0 Eðr0 Þ  ðk0 Gðr0 ; rÞ þ1δðr0  rÞÞ  Gðr0 ; rÞ  k0 Eðr0 Þ V ext

¼

I ∂R3 [ ∂V int

dσ 0 n^  ½Gðr0 ; rÞ  ∇  Eðr0 Þ  Eðr0 Þ  ∇  Gðr0 ; rÞ:

ð85Þ The first and third term under the volume integral cancel. The remaining volume integral over the Dirac-delta function yields the electric field EðrÞ if r lies in the integration domain V ext , otherwise 0. The first term under the surface integral can be simplified by use of Faraday's law in the form of Eq. (37). One obtains I  I dσ 0 n^  ½iωμG  H  E  ∇  G þ ∂R3 ∂V int ( EðrÞ : r A V ext : ð86Þ ¼ 0 : r A Vint Thus we can formally think of the field in the surrounding medium to be composed of two terms, which are given by the two surface integrals over ∂R3 and ∂V int , respectively. The integral over the boundary at infinity represents that part of the field that does not vanish at infinity, namely, the incident field Einc . By subtracting the incident field from both sides of the equation, we arrive at I 

dσ 0 iωμðn^ þ  HÞ  Gþ ðn^ þ  EÞ  ð∇  GÞ ∂V int ( sca E ðrÞ : r A V ext : ð87Þ ¼ Einc : r A V int Here we replaced n^ ¼  n^ þ , where n^ points out of V ext , hence into V int , while n^ þ points in the opposite direction. We further used the fact that all box products can be cyclically permuted, i.e. a  ðb  cÞ ¼b  ðc  aÞ¼c  ða  bÞ. At this point we should pause and inspect the last equation, because the main idea of how to solve the scattering problem in Waterman's method now becomes apparent. We want to compute the scattered field, given by the upper branch of Eq. (87). But to use this equation we need to know the electric and magnetic fields on the particle boundary ∂V int . They can be obtained from the lower branch of Eq. (87) in terms of the known incident field. However, what we actually wanted is not the scattered field for a specific incident field, but the T-matrix, which is independent of the incident field. The strategy is to expand

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

each of the fields in Eq. (87) in vector wave functions. We also need to employ a corresponding expansion of Green's function in Eq. (52), which turns out to be Gðr0 ; rÞ ¼ ik0

N cut n 2 X X X n ¼ 0 m ¼ n q ¼ 1

ð  1Þm

8 ð1Þ < ψ ð3Þ n;  m;q ðk0 rÞψ n;m;q ðk0 rÞ :

r o r0

ð3Þ : ψ ð1Þ n;  m;q ðk0 rÞψ n;m;q ðk0 rÞ :

r 4 r0

:

ð88Þ We substitute Eqs. (72)–(74) for the electric field, and corresponding expansions for the magnetic field, into Eq. (87); the surface fields in the integrand are represented by the internal field inside the particle. Thus each branch of Eq. (87) will give us a linear relation among expansion coefficients of the fields. The lower branch will give us the coefficients a of the incident field in terms of those of the internal field c. As we can see in Eq. (78) this is just the matrix Q. The upper branch of Eq. (87) gives us the coefficients p of the scattered field in terms of those of the internal field c. Inspection of Eq. (79) reveals that this just gives us the matrix RgQ . From the matrices Q and RgQ we finally obtain the T-matrix by use of Eq. (81). With a bit of imagination we see that substitution of the field and Green's function expansions into Eq. (87) will result in expressions for the Q- and RgQ-matrices that involve surface integrals over vector products of vector wave functions, i.e. I ð1Þ ðsÞ dσ 0 n^ þ  ψ n; ð89Þ  m;q ðk0 rÞ  ψ n0 ;m0 ;q0 ðkrÞ; ∂V int

where the kinds s¼1 and s¼3 are involved in the expressions for the matrices RgQ and Q, respectively. The exact formulas are somewhat messy, and the details are not particularly enlightening. In fact, to understand the main numerical properties of the method we already have everything in place. Let us inspect Eqs. (72)–(74), (81) and (89).

 We saw that the expansion over the degree n in Eqs.





(72)–(74) is, in fact, an infinite series. Thus the T-matrix is formally a matrix of infinite dimension that needs to be truncated at a finite expansion order N cut . In practice, this truncation index depends mainly on the size parameter k0 r of the particle, and to some extent on the refractive index and the particle shape. The larger the particle in relation to the wavelength, the higher Ncut , and the more matrix elements need to be computed. Hence the CPU time increases with particle size, usually on the order of CPU ðk0 rÞl where l typically ranges from 3 to 5.5, depending on the presence or absence of particle symmetries. (Increasing CPU time requirements with increasing size is, of course, a general feature of all numerical light scattering methods.) Eq. (81) reveals that the computation of the T-matrix requires a numerical inversion of the Q-matrix. This can become an ill-conditioned problem and introduce numerical inaccuracies. The integrands in Eq. (89) can be rapidly oscillating functions. This does not pose a numerical problem for small oscillation amplitudes. However, the vector wave functions of the third kind, j¼3, which contain spherical Neumann functions (see Eq. (68) and accompanying

13

text), can give rise to integrands with amplitudes that are many orders of magnitude larger than the actual value of the integral. In such cases the numerical evaluation of the integrals will run into severe loss-ofprecision problems. Storage of the large integrand values will consume much of the floating point precision, and summation of these large positive and negative values in the numerical integration process will leave little or no precision for the evaluation of the integral. Within the limits set by these numerical issues, Waterman's T-matrix method is among the most accurate and fast methods for solving far-field light scattering problems. It is particularly efficient for axisymmetric particles, but it can also be applied to non-axisymmetric shapes. A detailed derivation and complete list of all resulting surface integral expressions can be found in [28,22]. Specific implementations of the method and illustrative examples are discussed in [22,13]. Loss-of-precision problems are investigated in [29]. Alternative derivations of Waterman's T-matrix equations have been reported [30], which allow for more flexibility in the choice of functions that appear in the integrand of the surface integrals [13]. This can be important for near-field applications, where a choice of weighting functions different from that in Waterman's method can substantially improve the convergence of near-field results. Further refinements that build on Waterman's method exist. One such approach is the null-field method with discrete sources [31], which is numerically more stable, e.g., for particles that strongly depart from spherical shape. Another extension is known as the shape-matrix method [32], which efficiently solves the light scattering problem for particles of equal shapes and different sizes. A T-matrix reference data base [33] provides a comprehensive listing of articles on T-matrix theory and applications. Updates of the database are published periodically.

4. Comparison of different methods Comparisons of numerical methods for solving the MME should be interpreted with caution. First of all, the numerical performance can be strongly dependent on the individual implementation of the method. Some codes may also be tailored to specific applications. Run-time comparison among different methods may also be misleading when done for only a few isolated cases. This is because the size-scaling of CPU-time may vary considerably among different methods. It may also depend on numerical details of the implementation that are not directly related to the MME-solution method as such. For instance, in those methods that require a solution of a system of linear equations the size-scaling of CPU-time can be much more favourable when using the conjugantgradient rather than other methods, such as Gaussian elimination. Further, there may be factors other than the particle's size parameter, such as the refractive index or the shape, that can impact and bias a run-time comparison. For instance, a cubical particle may be a “home

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

14

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Table 1 Comparison of FDTD, DDA, SVM, and T-matrix methods. x denotes the size parameter. Numerical aspects Computational complexity Arbitrary morphologies Analytic orientation averaging Particle symmetries

FDTD

DDA

SVM

Oðx4 Þ þ – –

Oðx3 þ 3α log xÞ þ – –

a

Oðx3 Þ–Oðx4 Þ – ðþÞ d þ

T-matrix b

Oðx3 Þ  Oðx5:5 Þ – þ þ

c

a If N denotes the number of volume cells, then Nα is the number of iterations in the conjugate-gradient method for solving the system of equations in the DDA. b This is valid for spheroidal particles. c This is strongly dependent on the particle's symmetries, and to what extent symmetries are exploited in the implementation of the T-matrix method. d In most SVM implementations analytical orientation averaging is not possible. This only becomes possible when reformulating the method according to the suggestions in [24]. These ideas have also been applied to the point-matching method [26].

match” for the DDA with a rectangular grid, but a challenge to a T-matrix code using spherical coordinates. On the other hand, a particle with small-scale surface roughness may require a very fine discretisation in a DDA code, hence a high computation time, while a T-matrix code could be much faster. When performing method-intercomparisons, it is important to realise that different methods do not compute the same amount of information. This complicates method-intercomparisons. All methods that compute a Tmatrix deliver the complete information on a particle's monochromatic optical properties. Hence optical properties of particles in different or in random orientations can be computed rapidly. By contrast, methods that do not compute a T-matrix only yield information limited to one or a few selected orientations. Hence any run-time intercomparison of such methods, e.g., T-matrix and DDA, should be taken with a grain of salt. With these reservations in mind, Table 1 makes an attempt to summarise some of the main features of the methods that have been discussed in this review. As one can see, the main advantage of the FDTD and DDA is that they can be applied to particles with arbitrary morphologies. Thus these methods are most popular in fundamental studies that make use of morphologically complex particles resembling real particles encountered in nature. Their main limitation is the high computation time once the size parameter exceeds 20. However, with the availability of open-access DDA codes that can be run on parallel machines, and with the development of more powerful parallel computers, the range of particle size parameters that can be covered with DDA computations has steadily increased over the past 2 decades. Nevertheless, even by modern standards the computation time required for extending DDA or FDTD computations to size parameters that overlap with the geometric optics (GO) regime would be very high. Thus it is currently impractical to apply rigorous light scattering methods to arbitrary morphologies at intermediate size parameters in between wavelengthscale particles and the GO regime. The T-matrix method would, without the use of particle symmetries, be limited to a similar size range. However, this method lends itself easily to exploiting symmetries,

which drastically reduces CPU-time requirements. This extends not only the accessible range of size parameters, but it also allows us to model optical properties of ensembles of particle with many different sizes and shapes. The price one pays for this is that one is limited to simplified shapes with a more or less high degree of symmetry. T-matrix codes are very popular in applications in which the emphasis is on mimicking ensembles of a large number of particles. This can be relevant in theoretical investigations in which light-scattering computations are to be coupled to radiative transfer studies. A typical example would be studies of the radiative impact of tropospheric aerosol populations. Another important application would be the generation of look-up tables for use in either climate models or remote sensing retrieval systems. T-matrix calculations for mildly aspherical axisymmetric particles have been performed for large size parameters that overlap with the validity range of the GO approximation. For non-axisymmetric or strongly aspherical particles there do exist T-matrix methods that have been tailored to specific geometries or morphological properties, such as elongated spheroids or particles with smallscale surface roughness. Some of such custom-made methods can considerably extend the range of accessible size parameters, but they are limited to certain model geometries that are often characterised by a high degree of symmetry. Another important factor in method-intercomparisons is numerical stability. For instance, the DDA is often highly stable and mostly limited by computation time, while the T-matrix method may encounter numerical instabilities near the limit of its range of validity [34]. Convergence problems can also be encountered in the SVM. The ordinary point-matching method is often quite unstable, while the generalised point-matching method is much more robust. Finally, there are qualitative factors that can be important when choosing a numerical method. One is the availability of open-source computer codes, another is the degree of user-friendliness of a method and its implementations. The developers of DDA and T-matrix codes were among the first who made their computer codes freely available in the 1980s and early 1990s. They also

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

invested considerable efforts into making their codes userfriendly. This may have contributed to the fact that these methods have established themselves in the community as the most popular methods, as we saw in Fig. 1.

5. Further reading As advertised in the introduction, this review is meant to provide a first step to understanding light scattering theory and numerical methods. The derivations, although mathematically rigorous, were trimmed to the essential parts that are needed to understand the basic ideas in various numerical approaches. Also, the presentation was limited to an illustrative selection of numerical methods. This section will offer some guidance in taking the second step in studying numerical methods. It provides suggestions for review articles and monographs, which will lead further into the literature on more specialised topics. Among the classics in the field are the books by van de Hulst [35] and by Bohren and Huffman [36]. They offer insightful introductions to light scattering, Mie theory, approximate methods, and optical properties of bulk matter. A more recent treatment of Mie theory is presented in the book by Wriedt and Hergert [37]. A comprehensive treatment of electromagnetic scattering theory can be found in the book by Mishchenko et al. [22]; it contains a detailed introduction to the T-matrix method, as well as a short review of other numerical methods. An introduction to electromagnetic scattering theory for single particles and groups of particles, as well as a microphysical derivation of radiative transfer theory is given in the recent book by Mishchenko [9]. The monograph by Rother and Kahnert [13] discusses light scattering as a boundary value problem from a Green's function, separation-of-variables, and T-matrix perspective; also, a review of the use of group theory for exploiting particle symmetries in T-matrix computations is included. A quite exhaustive presentation of the finite-difference time-domain method, including technical aspects as well as applications to particles, periodic structures, circuits, antennas, and photonics, can be found in the book by Taflove and Hagness [38]. More concise reviews can be found in the recent articles by Sun et al. [11] and Cole et al. [12]. To the best of my knowledge, no monographs dedicated to the discrete dipole approximation have been published yet. An excellent starting point for learning more about this method is the review article by Lakhtakia and Mulholland [15], and the more recent review by Yurkin and Hoekstra [16]. The review by Lakhtakia and Mulholland also contains a discussion of the method of moments and its relation to the DDA. Numerous other light scattering methods and recent developments have not been treated in this brief introduction. Among the most important omissions is the superposition T-matrix method, which is based on an entirely different numerical scheme for computing the T-matrix. A detailed description of the method is given by Mackowski [39]. It has been applied extensively to light scattering by aggregated particles. More recently, this method has made it possible to study large groups of

15

particles that mimic macroscopic particulate and heterogeneous media. Thus it is now within reach of our computational capabilities to employ numerically exact solvers of the macroscopic Maxwell equations for modelling light propagation through discrete random media. An up-todate numerical algorithm and its numerical capabilities are discussed in [40]. Another important method is the boundary element method, which is particularly popular in engineering sciences. A review of the method is given by Buffa and Hiptmair [41]. Light scattering by laser beams of various shapes is treated in the generalised Lorenz–Mie theory, to which the book by Gouesbet and Gréhan is dedicated [42]. The null-field method with discrete sources and its numerical implementation is the subject of the monograph by Doicu et al. [43]. Numerical methods for light scattering by aggregated particles are discussed in the book by Borghese et al. [44]. A more complete overview of different numerical light scattering methods, their theoretical foundations, and their merits and drawbacks can be found in the book edited by Mishchenko et al. [45] and in the review by Kahnert [18,46]. A good resource for obtaining up-to-date review articles on all aspects of light scattering theory, radiative transfer, inverse problems, and applications is Light Scattering Reviews, a series of Springer books published annually and edited by A. Kokhanovsky.

Acknowledgements This work was funded by the Swedish Research Council (contract 621-2011-3346).

References [1] Mie G. Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen. Ann Phys 1908;25:377–445. [2] ScattPort. 〈http://www.scattport.org.〉. [3] Rodgers CD. In: Inverse methods for atmospheric sounding. Singapore: World Scientific; 2000. [4] Greiner W, Reinhardt J. In: Quantum electrodynamics. Berlin: Springer; 2008. [5] Weinberg S. In: Gravitation and cosmology. New York: Wiley; 1972. [6] Lee JM. In: Manifolds and differential geometry. Providence: American Mathematical Society; 2009. [7] Misner CW, Thorne KS, Wheeler JA. In: Gravitation. New York: W. H. Freeman; 1973. [8] Rothwell EJ, Cloud MJ. In: Electromagnetics. Boca Raton: CRC Press; 2001. [9] Mishchenko MI. Electromagnetic scattering by particles and particle groups. Cambridge: Cambridge University Press; 2014. [10] Yee SK. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Trans Antennas Propag 1966;14:302–7. [11] Sun W, Videen G, Fu Q, Tanev S, Lin B, Hu Y, Liu Z, Huang J. Finitedifference time-domain solution of light scattering by arbitrarily shaped particles and surfaces. In: Kokhanovsky A, editor. Light scattering reviews, vol. 6. Berlin: Springer; 2011. p. 75–113. [12] Cole JB, Okada N, Banerjee S. Advances in finite-difference timedomain calculation methods. In: Kokhanovsky A, editor. Light scattering reviews, vol. 6. Berlin: Springer; 2011. p. 115–175. [13] Rother T, Kahnert M. In: Electromagnetic wave scattering on nonspherical particles: basic methodology and simulations. 2nd editionBerlin: Springer; 2014 〈http://dx.doi.org/10.1007/978-3-64236745-8〉.

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i

16

M. Kahnert / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

[14] Purcell EM, Pennypacker CR. Scattering and absorption of light by nonspherical dielectric grains. Astrophys J 1973;186:705–14. [15] Lakhtakia A, Mulholland GW. On two numerical techniques for light scattering by dielectric agglomerated structures. J Res Natl Inst Stand Technol 1993;98:699–716. [16] Yurkin MA, Hoekstra AG. The discrete dipole approximation: an overview and recent developments. J Quant Spectrosc Radiat Transfer 2007;106:558–89. [17] Oguchi T. Scattering properties of oblate raindrops and cross polarization of radio waves due to rain: calculations at 19.3 and 34.8 GHz. J Radio Res Lab Jpn 1973;20:79–118. [18] Kahnert FM. Numerical methods in electromagnetic scattering theory. J Quant Spectrosc Radiat Transfer 2003;79–80:775–824. [19] Oguchi T, Hosoya Y. Scattering properties of oblate raindrops and cross polarization of radio waves due to rain II. Calculations at microwave and millimeter wave regions. J Radio Res Lab Jpn 1974;21:191–259. [20] Waterman PC. Matrix formulation of electromagnetic scattering. In: Proceedings of IEEE 1965; 53: 805–12. [21] Mishchenko MI, Travis LD, Mackowski DW. T-matrix computations of light scattering by nonspherical particles: a review. J Quant Spectrosc Radiat Transfer 1996;55:535–75. [22] Mishchenko MI, Travis LD, Lacis AA. In: Scattering, absorption, and emission of light by small particles.Cambridge: Cambridge University Press; 2002. [23] Mackowski DW, Mishchenko MI. Calculation of the T matrix and the scattering matrix for ensembles of spheres. J Opt Soc Am A 1996;13: 2266–78. [24] Schulz FM, Stamnes K, Stamnes JJ. Scattering of electromagnetic waves by spheroidal particles: a novel approach exploiting the Tmatrix computed in spheroidal coordinates. Appl Opt 1998;37: 7875–96. [25] Mackowski DW. Discrete dipole moment method for calculation of the T matrix for nonspherical particles. J Opt Soc Am A 2002;19: 881–93. [26] Nieminen TA, Rubinsztein-Dunlop H, Heckenberg NR. Calculation of the T-matrix: general considerations and application of the pointmatching method. J Quant Spectrosc Radiat Transfer 2003;79–80: 1019–29. [27] Johnson BR. Invariant imbedding T matrix approach to electromagnetic scattering. Appl Opt 1988;27:4861–73. [28] Tsang L, Kong JA, Shin RT. In: Theory of microwave remote sensing. New York: John Wiley & Sons; 1985. [29] Somerville WRC, Auguié B, Ru ECL. Severe loss of precision in calculations of the T-matrix integrals. J Quant Spectrosc Radiat Transfer 2012;113:524–35. [30] Schmidt K, Rother T, Wauer J. The equivalence of applying the extended boundary condition and the continuity conditions for

[31]

[32]

[33]

[34]

[35] [36] [37] [38] [39] [40]

[41]

[42] [43]

[44]

[45] [46]

solving electromagnetic scattering problems. Opt Commun 1998;150:1–4. Wriedt T. Studies of light scattering by complex particles using the null-field method with discrete sources. In: Kokhanovsky A, editor, Light scattering reviews, vol. 2. Berlin: Springer; 2007. p. 269–94. Petrov D, Synelnyk E, Shkuvatov Y, Videen G. TheT-matrix technique for calculations of scattering properties of ensembles of randomly oriented particles with different size. J Quant Spectrosc Radiat Transfer 2006;102:85–110. Mishchenko MI, Videen G, Babenko VA, Khlebtsov NG, Wriedt T. T-matrix theory of electromagnetic scattering by particles and its applications: a comprehensive reference database. J Quant Spectrosc Radiat Transfer 2004;88:357–406. Yurkin MA, Kahnert M. Light scattering by a cube: accuracy limits of the discrete dipole approximation and the T-matrix method. J Quant Spectrosc Radiat Transfer 2013;123:176–83. van de Hulst HC. In: Light scattering by small particles. New York: John Wiley & Sons, Inc.; 1957. Bohren CF, Huffman DR. In: Absorption and scattering of light by small particles. Weinheim: Wiley-VCH; 1983. Wriedt T, Hergert W. In: The Mie theory: basics and applications. Berlin: Springer; 2012. Taflove A, Hagness SC. In: Computational electrodynamics—the finite-difference time-domain method. Boston: Artech House; 2005. Mackowski DW. Calculation of total cross sections of multiplesphere clusters. J Opt Soc Am A 1994;11:2851–61. Mackowski DW, Mishchenko MI. A multiple sphere T-matrix Fortran code for use on parallel computer clusters. J Quant Spectrosc Radiat Transfer 2011;112:2182–92. Buffa A, Hiptmair R. Galerkin boundary element methods for electromagnetic scattering. In: Ainsworth M, Davies P, Duncan D, Rynne B, Martin P, editors. Topics in computational wave propagation, Lecture notes in computational science and engineering, vol. 31. Berlin, Heidelberg: Springer; 2003. p. 83–124. Gouesbet G, Gréhan G. In: Generalized Lorenz–Mie theories. Berlin: Springer; 2011. Doicu A, Wriedt T, Eremin YA. Light scattering by systems of particles.Null-field method with discrete sources: theory and programs. Berlin: Springer; 2006. Borghese F, Denti P, Saija R. Scattering from model nonspherical particles.Theory and applications to environmental physics. Berlin: Springer; 2010. Mishchenko MI, Hovenier JW, Travis LD. In: Light scattering by nonspherical particles. San Diego: Academic Press; 2000. Kahnert M. Electromagnetic scattering by nonspherical particles: recent advances. J Quant Spectrosc Radiat Transfer 2010;111: 1788–90.

Please cite this article as: Kahnert M. Numerical solutions of the macroscopic Maxwell equations for scattering by nonspherical particles: A tutorial review. J Quant Spectrosc Radiat Transfer (2015), http://dx.doi.org/10.1016/j. jqsrt.2015.10.029i