Open Access
Add to BibSonomyAbstract
We present a time-domain formulation of electrodynamics based on the self-consistent derivation of the electromagnetic field in a linear, dispersive, lossy object via the coupled dipole method.
©2008 Optical Society of America
1. Introduction
The interaction of an electromagnetic field with an inhomogeneous, lossy, dispersive three-dimensional object is an important problem in not only optics but many areas of photonics and surface science. There exists a plethora of techniques for solving Maxwell’s equations in a complex environment, in the time-harmonic regime [1], including the finite element method (FEM), the multiple multipole method (MMP), the fast-multipole method (FMM), the method of moments (MoM) and the coupled dipole method (CDM), also called discrete dipole approximation [2, 3]. Both CDM and MoM are based on the volume integral representation of electromagnetic fields. The formal equivalence between the CDM and the MoM is described in Ref. [4].
On the other hand, in the time domain, electromagnetic scattering by an arbitrary three-dimensional structure is usually addressed using the finite difference in time domain (FDTD) method [5, 6, 7]. In the FDTD, the differential form of the time-dependent Maxwell equations is discretized in both space and time. As a result, the entire computational window must be discretized and not just the scatterer. This can lead to a prohibitively large memory requirement. Furthermore, numerical dispersion and/or instabilities may occur when the fields are propagated over large distances (compared to the wavelength). To circumvent these issues, a frequency domain solution of the integral form of Maxwell’s equations can be associated to a temporal Laplace or Fourier transform to study the time evolution of electromagnetic quantities. This leads to an exact frequency-domain equation, which can be solved by one of the standard frequency-domain techniques. We shall use the CDM. This approach has the advantage of restricting the computational domain to the scatterers only and the global error on the computed field is determined mainly by the mesh size over which the integral equation is discretized. Furthermore, because the propagation of the fields at arbitrary distances from the scatterers can be done analytically, once the fields inside the scatterer are known, the fields anywhere else can be computed at a minimal cost. There is, however, a potential hurdle in solving a discretized version of Maxwell’s equations in the frequency domain; if the computation of the fields is to be repeated for a large number of frequencies, the computation time can quickly become prohibitive. This can be avoided by performing the time-harmonic field computation iteratively using a conjugate gradient method with an efficient estimate of the initial field [8]. Note that iterative procedures have been previously used to process efficiently a sequence of CDM calculations but only to account for small changes in the geometry of the scattering problem or the refractive index of the scatterer, e.g. [9, 10], and not to formulate a scattering problem in the time-domain.
2. The coupled-dipole method
To illustrate this approach let us consider a three-dimensional scatterer. In the standard CDM the scatterer is represented by a cubic array of N polarizable subunits [11, 4]. If the subunits are small compared to the wavelength inside the object, the local field at each subunit can be considered as uniform and can be expressed as:
where r i is the position of the i-th subunit, E 0(r i,ω) is the incident field at the position r i, and T⃡ is the free-space field susceptibility tensor [12], i.e., the electric response to an electric dipole. α⃡(r j,ω), the polarizability tensor of subunit j, is expressed as [13]:
where k 0=2π/λ=ω/c is the wave number, I⃡ the identity tensor, and α⃡0(r j,ω) the Clausius-Mossotti polarizability:
In the previous expression d is the lattice spacing of the CDM lattice, and ε⃡(r j,ω) the relative permittivity tensor of the scatterer at subunit j. Once the fields E(r i,ω) are obtained by solving the linear system, Eq. (1), the field scattered by the object at an arbitrary position r is given by:
Notice that if the objects under study are not in free space space, we would only need to change the dyadic tensor T⃡ to account for the environment [14]. The overall approach remains the same.
3. The coupled-dipole method in time domain
Let us now consider the case where the incident field is an electromagnetic pulse with a Gaussian envelop ℱ(t) and a spectrum centered on frequency f 0:
Let F(ω) be the Fourier transform of ℱ(t), the Fourier transform of the incident field 𝓔(r,t), is then F(ω)E 0(r,ω). The linearity of Maxwell’s equations ensures that solving the time harmonic scattering problem with the incident field F(ω)E 0(r,ω), and computing the inverse Fourier transform of the resulting time harmonic fields yields the correct time evolution for the scattered fields. In practice M values are taken in the frequency domain for F(ω) in accordance with the Nyquist-Shannon sampling theorem. Hence the main computational cost stems from having to solve Eq. (1) for a large number of frequencies. This is where an extrapolation procedure [8] can be used to reduce drastically the computation time. This procedure is also known as a marching-on-in frequency method.
In practice, to find the local field at each lattice site at a given frequency, Eq. (1) must be solved. This corresponds to solving a linear system with size 3N×3N. When the number of discretization cells is large the linear system is solved iteratively. We shall use the well known conjugate gradient (CG) method which yields the exact result after 3N iterations. Before we describe our procedure in detail let us rewrite Eq. (1), taken at frequency ω m, in the symbolic form :
where p(ω m) is a vector of length 3N which comprises the dipole moments of the discretization cells, [i.e., p(r i,ω m)=α⃡(r i,ω m)E(r i,ω m)], A is a square matrix which contains the field susceptibility tensors, and E 0 a vector of length 3N which contains the incident field. D is a diagonal matrix which contains the inverse of the polarizabilities. For each iteration, the CG method requires us to compute the matrix vector product A(ω m)p(ω m). This product can be computed efficiently by using the fact that the field-susceptibility tensor actually depends on the relative positions of the source-point and the field-point, rather than on their absolute locations. This means that the matrix-vector product can be cast as a convolution product and computed using fast-Fourier transform techniques [15, 16]. However, so far we have solved the scattering problem for a single frequency whereas Eq. (6) must be solved for M frequencies in order to be able to derive the time evolution of the fields via inverse Fourier transform. Therefore, when the CDM computation is performed for the m th value of the frequency (m=1, …,M), the convergence of the CDM can be accelerated by using an initial guess for the fields (or equivalently the dipole moments). Usually, the incident field is used as the initial guess, but in our case, as we are doing the same CDM calculation for a discrete set of frequencies we can use the solution obtained at frequency m-1 as the initial guess for the computation at frequency m. In fact the initial guess for the CG method can be refined further by using several frequencies:
For the m th frequency, the initial guess is taken as a linear combination of the K previous frequencies, where the coefficients a k are found in minimizing the quantity:
The minimization procedure leads to a linear system of size K×K where the coefficients a k are the unknowns.
4. Results
4.1. Comparison between the CDM and FDTD method
We first start by illustrating the validity of the CDM approach by comparing it to a commercial FDTD code [17]. We consider a homogeneous sphere with relative permittivity ε=2.25 and radius a=120 mm. The parameters of the incident pulse are f 0=2 GHz (λ≈150 mm) and τ=4 ns as presented in Fig. 1. The incident field is a plane wave propagating along z and polarized along x.
For both methods the discretization used is the same, i.e., the size of the subunit is fixed at d=9.6 mm. The number of frequencies used to describe the pulse in the CDM is M=256. Fig. 2 shows the time evolution of the amplitude of the electric field at the center of the sphere as computed by both methods. An excellent agreement is found.
Fig. 2. Component x of the field versus time when the point of observation is located at the center of the sphere. In bold line the CDM is used and in dashed line the FDTD is used.
4.2. Nonhomogeneous, lossy scatterer
We now consider a more complex situation. Consider an inhomogeneous sphere with a radius a=120 mm, a real part of the relative dielectric permittivity profile given by: ε r(ρ)=2.25+2sin2(πρ 2/a 2) where ρ is the distance between the center of the sphere and a point inside the sphere, and an imaginary part for the relative permittivity, ε i=0 or 0.5. The parameters of the incident pulse are the same as previously (f 0=2 GHz (λ≈150 mm) and τ=4 ns with a linear polarization along the x axis.)
We first study the influence of the initialization procedure presented in Sect. 3 on the convergence of the method. The iterative method used is a classical CG method with a tolerance fixed at 10-3. The number of frequencies chosen to discretize the pulse in the frequency domain is M=256. Figure 3(a) shows the number of iterations needed in the CG method to find the internal fields at all the frequencies for different values of K: K=0, …, 5. Notice that K=0 correspond to the use of the Born approximation as the initial guess, and K=1 corresponds to the use of the local field obtained for the frequency i-1 as the initial guess for the frequency i. In the case where ε i=0 the number of iterations is reduced by 20%, compared to the Born approximation, when we use K≥2. Actually, the speed-up factor depends on the frequency. When the computation is done at the frequency of a morphological resonance (which corresponds to a large number of iterations at a given frequency as represented by the bumps in Fig. 3(b) for the dashed line) the number of iterations for the Born approximation (K=0) and the case K=5 are similar. Away from these resonances, however, a value of K > 0 will accelerate the convergence significantly. In the case where ε i=0.5 using K=5 improves the convergence by a factor of 2 compared to K=0 (Born approximation). The increased perfomance of the initial guess approach in this case results from the presence of the term ε i which dampens the morphological resonances [Fig. 3(c)]. Incidentally we emphasize that these issues with morphological resonances can be circumvented by using Laplace transform techniques, however, this topic is beyond the scope of this paper.
We now proceed to compute the time evolution of the fields. Figures 4 shows the component x of the electric field versus time for a point of observation located at the top of the sphere [Fig. 4(a)] and at the center of sphere [Fig. 4(b)]. One can see that when the field is computed close to the surface of the sphere the telltale of a creeping wave can be observed (successive bursts in the amplitude of the field). This particular wave propagates along the surface of the sphere and the 5 ns interval between successive bursts corresponds to the time it takes the wave to travel around the sphere. Logically, when the field is recorded at the center of the sphere no such effect is observed.
4.3. Scatterer with material dispersion
We now consider the scattering of light by a gold particle. The relative permittivity of gold varies with the frequency, therefore dispersion must be accounted for in order to get an accurate answer. With our time-domain formulation of the CDM dispersion is easy to take into account as each frequency forming the spectrum of the incident field is treated separately. Therefore, there is no need to express the optical constant as a function of time. The radius of the sphere is a=60 nm and the frequency of the incident field is f 0=0.57 PHz which corresponds to a wavelength λ0=525 nm. Figure 5 shows the spectrum of the incident field. The relative permittivity of gold is described using a Drude model, improved to take into account the two interband transition (λ≈470 nm and 330 nm) which are present in our spectrum:
Fig. 3. (a) Number of total iterations required to solve the scattering problem for all the frequencies versus K. The computation is done for a sphere with a radius a=120 mm and f 0=2 GHz with ε r=2.25+2sin2(πρ 2/a 2). The line with squares pertains to ε i=0.5i whereas the line with circles pertains to ε i=0. (b) and (c) Number of iterations needed to solve our linear system with the conjugate gradient method at each frequency for K=0 (bold line) and K=5 (dashed line). (b) ε i=0. (c) ε i=0.5i.
Fig. 4. Component x of the electric field versus time. (a) The point of observation is located at the top of the sphere. (b) The point of observation is located at the center of the sphere.
Fig. 5. Spectrum of the pulse used on a gold sphere of radius a=60 nm. The spectrum contains the visible range and is centered at f 0=0.57 PHz (λ0=525 nm).
where the constants ω p, Γ, ε∞ and the functions G 1(ω), G 2(ω) are given in Ref. [18]. We note that the relative permittivity of gold in the spectral range considered can be very large [for λ=1500 nm, ε=-100+10i] and in that case, as shown in Ref. [4], the convergence of the standard formulation of theCDM is impossible to reach. We can, however, use the polarizability derived in Refs. [19, 20] which, for spherical scatterers, improves drastically the convergence of the CDM. Notice that this new form of polarizability gives a tensorial polarizability but as showed in Sect. 2 this is taken into account easily with the CDM. We also use a very fine discretization (d=3 nm) and for each frequency we check that the local fields are computed correctly by comparing the extinction cross section obtained from the CDM to Mie theory. Hence the relative error between Mie theory and the CDM is always less than 0.5% in the range above the 1/e attenuation.
Figure 6 shows the time profile of the x component of the field. The solid line is the incident field, the dashed line is the field at the front of the sphere, and the dashed line with circle the field at the back of the sphere. As the sphere is small compared to the wavelength of illumination there is no “retardation effect” between the field at the top and the bottom of the sphere. However, one can observe that different parts of the pulse are attenuated in a different way by the sphere due to material dispersion.
5. Conclusion
We have extended the CDM to electromagnetic waves that are not time harmonic. Material anisotropy and dispersion are easily taken into account in the method. The main advantage of the method presented here is that, unlike the FDTD, the global error on the computed fields depend mainly on the spatial discretization of the object, since the space between scatterers need not be discretized. While a study of the respective computational performance of the two approaches is beyond the scope of this work, a recent work by Yurkin and co-workers [3] comparing CDM and FDTD in the time-harmonic regime for lossless dielectric scatterers indicates that, quite logically, the performance of each method depends significantly on the scattering configuration. In our time-domain formulation, we anticipate that when extended environments such as substrates and multilayer systems are to be included in the problem, or when scatterers separated by large (compared to the wavelength) distances are present, the CDM will have the distinct advantage over the FDTD of being able to account for the extended environment, or propagate the fields over large distances, using semi-analytical techniques. Similarly, material disperion can be handled using experimental values of the optical constants without the need for a time domain model of the optical constants. However, a quantification of the computational cost of each method for various configurations is left for future studies.
Fig. 6. Solid line: the incident field. Dashed line: field at the top of the sphere. Dashed line with circles: field at the bottom of the sphere.
We emphasize that this work provides the foundation for extending all previous CDM studies to time-varying fields. In particular, because the CDM uses a formulation based of the field susceptibility tensor, it is well suited for the study of electrodynamics problems such as spontaneous emission in complex geometries [16, 21, 22], optical trapping [23, 24], and optical binding [25]. Furthermore, it has been shown previously that the CDM can be used to model the refractive index profile reconstruction of an arbitrary nanostructure from its far-field signature using optical diffraction tomography [26, 27]. The present work gives us a starting point for the modelling of optical diffraction tomography using time-varying field [28].
References and links
1. F. M. Kahnert, “Numerical methods in electromagnetic scattering theory,” J. Quant. Spectrosc. Radiat. Transf. 79–80, 775–824 (2003). [CrossRef]
2. M. A. Yurkin and A. G. Hoekstra, “The discrete dipole approximation: An overview and recent developments,” J. Quant. Spectrosc. Radiat. Transf. 106, 558–589 (2007). [CrossRef]
3. M. A. Yurkin, A. G. Hoekstra, R. S. Brock, and J. Q. Lu, “Systematic comparison of the discrete dipole approximation and the finite difference time domain method for large dielectric scatterers,” Opt. Express 15, 17,902–17,911 (2007). [CrossRef]
4. P. C. Chaumet, A. Sentenac, and A. Rahmani, “Coupled dipole method for scatterers with large permittivity,” Phys. Rev. E 70, 036,606–6 (2004). [CrossRef]
5. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat. 14, 302–307 (1969).
6. A. Taflove and M. E. Brodwin, “Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell’s equations,” IEEE Trans. Microwave Theory Tech. 23, 623”630 (1975). [CrossRef]
7. A. Taflove and S. C. Hagness, “Computational Electrodynamics: The Finite-Difference Time-Domain Method”, 3rd edition, (Artech House Publishers, 2005).
8. Z. Q. Peng and A. G. Tijhuis, “Transient scattering by a lossy dielectric cylinder: marching-on-in-frequency approach,” J. Elect. Waves Appl. 7, 739–763 (1993). [CrossRef]
9. K. Muinonen and E. Zubko, “Optimizing the discrete-dipole approximation for sequences of scatterers with identical shapes but differing sizes or refractive indices,” J. Quant. Spectrosc. Radiat. Transf. 100, 288–294 (2006). [CrossRef]
10. Y. Okada, I. Mann, I. Sano, and S. Mukai, “Acceleration of the iterative solver in the discrete dipole approximation: Application to the orientation variation of irregularly shaped particles,” J. Quant. Spectrosc. Radiat. Transf. 109, 1461–1473 (2008). [CrossRef]
11. E. M. Purcell and C. R. Pennypacker, “Scattering and absorption of light by nonspherical dielectric grains,” Astrophys. J. 186, 705–714 (1973). [CrossRef]
12. J. D. Jackson, Classical Electrodynamics (Wiley, 1975), 2nd ed.
13. P. C. Chaumet, “Comment on “Trapping force, force constant, and potential depths for dielectric spheres in the presence of spherical aberrations,”” Appl. Opt. 43, 1825–1826 (2004). [CrossRef] [PubMed]
14. A. Rahmani, P. C. Chaumet, F. de Fornel, and C. Girard, “Field propagator of a dressed junction: Fluorescence lifetime calculations in a confined geometry,” Phys. Rev. A 56, 3245–3254 (1997). [CrossRef]
15. J. J. Goodman and P. J. Flatau, “Application of fast-fourier-transform techniques to the discrete-dipole approximation,” Opt. Lett. 16, 1198–1200 (2002). [CrossRef]
16. F. Bordas, N. Louvion, S. Callard, P. C. Chaumet, and A. Rahmani, “Coupled dipole method for radiation dynamics in finite photonic crystal structures,” Phys. Rev. E 73, 056,601 (2006). [CrossRef]
17. RSoft Inc., RSoft Fullwave FDTD code, http://www.rsoftdesign.com.
18. P. G. Etchegoin, E. C. Le Ru, and M. Meyer, “An analytic model for the optical properties of gold,” J. Chem. Phys. 125, 164,705–3 (2006). [CrossRef]
19. A. Rahmani, P. C. Chaumet, and G. W. Bryant, “On the Importance of Local-Field Corrections for Polarizable Particles on a Finite Lattice: Application to the Discrete Dipole Approximation,” Astrophys. J. 607, 873–878 (2004). [CrossRef]
20. A. Rahmani, P. C. Chaumet, and G. W. Bryant, “Local-field correction for an interstitial impurity in a crystal,” Opt. Lett. 27, 430–432 (2002). [CrossRef]
21. A. Rahmani, P. C. Chaumet, and F. de Fornel, “Environment-induced modification of spontaneous emission: Single-molecule near-field probe,” Phys. Rev. A. 63023819 (2001). [CrossRef]
22. A. Rahmani and G. W. Bryant, “Spontaneous emission in microcavity electrodynamics,” Phys. Rev. A. 65033817 (2002). [CrossRef]
23. P. C. Chaumet, A. Rahmani, and M. Nieto-Vesperinas, “Optical trapping and manipulation of nano-objects with an apertureless probe,” Phys. Rev. Lett. 88123601 (2002). [CrossRef] [PubMed]
24. A. Rahmani and P. C. Chaumet, “Optical trapping near a photonic crystal,” Opt. Express 14, 6353–6358 (2006). http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-13-6353 [CrossRef] [PubMed]
25. P. C. Chaumet and M. Nieto-Vesperinas, “Optical binding of particles with or without the presence of a flat dielectric surface,” Phys. Rev. B 64, 035422 (2001). [CrossRef]
26. P. C. Chaumet, K. Belkebir, and A. Sentenac, “Three-dimensional sub-wavelength optical imaging using the coupled dipole method,” Phys. Rev. B 69, 245,405–7 (2004). [CrossRef]
27. K. Belkebir, P. C. Chaumet, and A. Sentenac, “Influence of multiple scattering on three-dimensional imaging with optical diffraction tomography,” J. Opt. Soc. Am. A 23, 586–569 (2006). [CrossRef]
28. A. Dubois, J. M. Geffrin, K. Belkebir, and M. Saillard, “Imaging of dielectric cylinders from experimental stepped-frequency data,” Appl. Phys. Lett. 88, 164,104 (2006). [CrossRef]
