Traditional and new simulation techniques for nanoscale optics and photonics a I. Tsukerman*a, F. Čajkoa, A.P. Sokolovb Department of Electrical & Computer Engineering, The University of Akron, OH 44325-3904, USA b Department of Polymer Science, The University of Akron, OH 44325-3909, USA
ABSTRACT

Several classes of computational methods are available for computer simulation of electromagnetic wave propagation and scattering at optical frequencies: Discrete Dipole Approximation, the T-matrix − Extended Boundary Condition methods, the Multiple Multipole Method, Finite Difference (FD) and Finite Element (FE) methods in the time and frequency domain, and others. The paper briefly reviews the relative advantages and disadvantages of these simulation tools and contributes to the development of FD methods. One powerful tool – FE analysis − is applied to optimization of plasmon-enhanced AFM tips in apertureless near-field optical microscopy. Another tool is a new FD calculus of “Flexible Local Approximation MEthods” (FLAME). In this calculus, any desirable local approximations (e.g. scalar and vector spherical harmonics, Bessel functions, plane waves, etc.) are seamlessly incorporated into FD schemes. The notorious ‘staircase’ effect for slanted and curved boundaries on a Cartesian grid is in many cases eliminated – not because the boundary is approximated geometrically on a fine grid but because the solution is approximated algebraically by suitable basis functions. Illustrative examples include problems with plasmon nanoparticles and a photonic crystal with a waveguide bend; FLAME achieves orders of magnitude higher accuracy than the standard FD methods, and even than FEM.

Keywords: wave propagation, computational methods, flexible approximation, photonic crystals, plasmon particles, apertureless near-field microscopy, AFM tips, field enhancement, optimization.

1. INTRODUCTION
The paper presents an overview of computational methods – both traditional and new – with applications in photonics. Simulation examples are given in Section 4 and include plasmon particles and resonances and wave propagation in a photonic crystal. Electromagnetic formulations used throughout the paper are those of classical electrodynamics in the frequency domain, with equivalent material parameters − permittivity and permeability (in general, complex). Time-domain versions of the computational methods do exist but are not discussed here. Although the use of effective material parameters has its limitations, the effective permittivity does adequately represent the key physical effects in many important cases, possibly with adjusted material parameters [17], [32], [35]. Our ultimate goal is to assemble a set of complementary simulation tools for comprehensive analysis and crossvalidation of results. Such cross-validation is critical for several reasons: (i) wave propagation and scattering problems are quite complex, and even more so around the resonance conditions; (ii) our intuition is lacking and often cannot provide even qualitative guidance; (iii) accurate comparison with measurements is extremely difficult on the nanoscale. The traditional methods we use include Finite Element (FE) and Finite Difference (FD) analysis, the Discrete-Dipole Method (with the DDSCAT code by Draine & Flatau [8], [9]), and multipole-multicenter expansions (with or without the T-matrix flavor [25]). In addition to the conventional techniques, one of the authors has developed a new finite difference calculus of “Flexible Local Approximation MEthods” (FLAME), Section 3. This method incorporates accurate approximating functions (such as plane waves, cylindrical or spherical harmonics, etc.) into the difference scheme on simple Cartesian grids. While the computational complexity of FLAME is comparable with that of traditional Finite-Difference (FD) methods, the accuracy can be much higher than the accuracy of not only FD methods but of FE analysis as well.

*

igor@uakron.edu; phone 1-330-972-8041; http://www.ecgf.uakron.edu/~igor/Research/Igor_Research.htm

Plasmonics: Metallic Nanostructures and Their Optical Properties III, edited by Mark I. Stockman, Proceedings of SPIE Vol. 5927 (SPIE, Bellingham, WA, 2005) · 0277-786X/05/$15 · doi: 10.1117/12.615011

Proc. of SPIE 59270L-1

2. EXISTING METHODS
The Electromagnetic Model. Questions do arise about the applicability of the macroscopic permittivity of metals for nanoparticles. Several physical effects are at play for small particles. Electron scattering at particle boundaries affects the mean free path of conduction electrons, which phenomenologically can be described by an increasing damping constant in the expression for the permittivity as a function of frequency [35], [36]. Furthermore, at the surfaces of silver particles as an example, due to quantum effects the 5s electron density ‘spills out’ into the vacuum, where 5s electronic oscillations are not screened by the 4d electrons [20], [21]. Various physical mechanisms affecting the value of the effective dielectric constant in individual nanoparticles and in particle clusters are discussed in detail in the physical literature [17], [19], [20], [21], [32], [35], [36]. Fortunately, the consensus is that plasmon resonances can in a wide range of applications be described with sufficient accuracy by classical electrodynamic models with effective permittivity. The material parameter can be adjusted for the particle size [32], [35], [36]. The classical electromagnetic model with effective material parameters is used throughout the paper. Analytical Solutions. As an analytical problem, scattering of electromagnetic waves from dielectric objects is quite involved. Closed-form solutions are available only for a few cases (see e.g. [26]): an isotropic homogeneous sphere (the classic Lorenz–Mie-Debye case); concentric core-mantle spheres; concentric multilayered spheres; radially inhomogeneous spheres; a homogeneous infinite circular cylinder; an infinite elliptical cylinder; homogeneous and coremantle spheroids. For objects other than homogeneous spheres or infinite cylinders, the complexity of analytical solutions (if they are available) is so high that the boundary between analytical and numerical methods becomes blurred. At present, further extensions of purely analytical techniques seem unlikely. On the other hand, with the available analytical cases in mind, local analytical approximations to the field are substantially easier to construct than global closed-form solutions. Such local analytical approximations can be incorporated into high-accuracy “Flexible Local Approximation Methods” (FLAME), Section 3. T-matrix methods are very widely used in scattering problems [25], [26]. Mishchenko et al. collected a comprehensive database of references [29] and have developed a T-matrix software package [28]. If a monochromatic wave impinges on a scattering dielectric object of arbitrary shape, both the incident and scattered waves can be expanded into spherical harmonics around the scatterer. If the electromagnetic properties of the scatterer (the permittivity and permeability) are linear, then the expansion coefficients of the scattered wave are linearly related to the coefficients of the incident wave. The matrix governing this linear relationship is called the T- (‘transition’) matrix. For a collection of N scattering particles, the overall field can be sought as a superposition of the individual harmonic expansions around each scatterer. The transformation of vector spherical harmonics centered at one particle to harmonics around another one is accomplished via well-established translation and rotation rules (Theorems) (e.g. [22], [26], [27], [53]). Self-consistency of the multi-centered expansions then leads to a linear system of equations for the expansion coefficients. Since the system matrix is dense, the computational cost may become prohibitively high if the number of scatterers is large. For spherical, spheroidal and other particles that admit a closed-form solution of the wave problem (see above), the Tmatrix can be found analytically. For other shapes, the T-matrix is computed numerically. If the scatterer is homogeneous, the “Extended Boundary Condition Method” (EBCM) (e.g. [1], [26]) is usually the method of choice. EBCM is a combination of integral equations for equivalent surface currents and expansions into vector spherical harmonics [13], [40]. While the T-matrix method is quite suitable for a moderate number of isolated particles and is also very effective for random distributions and orientations of particles (e.g. in atmospheric problems), it is not designed to handle large continuous dielectric regions. It is possible, however, to adapt the method to particles on an infinite substrate at the expense of additional analytical, algorithmic and computational work: plane waves reflected off the substrate are added to the superposition of spherical harmonics scattered from the particles themselves [6], [49]. In the Multiple Multipole Method (MMP), the computational domain is decomposed into homogeneous subdomains, and an appropriate analytical expansion – most frequently, a superposition of multipole expansions as the name suggests – is introduced within each of the subdomains. A system of equations for the expansion coefficients is obtained by collocation of the individual expansions at a set of points on subdomain boundaries. MMP is a common tool in computational electromagnetics and optics. Moreno et al. [30] applied it specifically to the computation of plasmon resonances. They also proposed a heuristic procedure for the choice of the centers of multipole expansions.

Proc. of SPIE 59270L-2

A shortcoming of MMP is that no general systematic procedure for choosing the centers of the multiple multipole expansions is available. The centers are chosen heuristically, which makes it difficult to evaluate and systematically improve the accuracy and convergence characteristics of the method. The Discrete-Dipole Method belongs to the general category of integral equation methods but admits a very simple physical interpretation. Scattering bodies are approximated by a collection of dipoles, each of which is directly related to the local value of the polarization vector. Starting with the volume integral equation for the electric field, one can derive a self-consistent system of equations for the equivalent dipoles [8], [9], [11], [18], [33]. The method has gained popularity in the simulation of plasmon particles (as well as other scattering problems) [15], [23], [42] because of its conceptual simplicity, relative ease of use and the availability of public domain software DDSCAT [8], [9] by Draine & Flatau. DDM has a few disadvantages typical of integral-equation methods. First, the treatment of singularities in DDM is quite involved [18], [33]. Second, the system matrix for the coupled dipoles is dense, and therefore the computational time increases very rapidly with the increasing number of dipoles. If the dipoles are arranged geometrically on a regular grid, the numerical efficiency can be improved by using Fast Fourier Transforms to speed up matrix-vector multiplications in the iterative system solver. However, such a regular arrangement of the sources is not always adequate or convenient in practice. DDM shares one additional disadvantage not with integral-equation methods but rather with finite-difference algorithms. Since the equivalent dipoles reside at the nodes of a Cartesian grid, there is a ‘staircase’ effect in the representation of curved or slanted material boundaries. In DDM simulations (e.g. [10], [23]), there are typically thousands of dipoles in each particle and tens of thousands of dipoles for problems with a few particles on a substrate. As an example, in [23] 11,218 dipoles are used in the particle and 93,911 dipoles in the particle and substrate together, so that the overall system of equations has a dense matrix of dimension ~280,000. The challenge in Finite Difference methods is to overcome the ‘staircase’ artifacts at boundaries not aligned with grid lines. A great variety of approaches to reduce or eliminate the staircase effect in Finite Difference – Time Domain (FDTD) have been proposed [1], [34], [39], [43], [52]. Each case is a trade-off between the simplicity of the original Yee scheme and the ability to represent the interface boundary conditions accurately. On one side of this spectrum lie various adjustments to the Yee scheme: changes in the time-stepping formulas for the magnetic field [5] or heuristic homogenization of material parameters [4], [52]. In some cases, the second order of the FDTD scheme is maintained by including additional geometric parameters [7] or by using partially filled cells [54]. On the other side of the spectrum are Finite Volume – Time Domain methods (FVTD) [2], [43], [51] and the Finite Element Method (FEM). They operate on complex (typically tetrahedral) meshes, which gives higher accuracy, but at a higher algorithmic and computational cost. On tetrahedral meshes, material interfaces are represented much more accurately than on Cartesian grids. (However, adaptive Cartesian grids have also been advocated [50], with cell refinement at the boundaries.) Edge elements [3] with Perfectly Matched Layers [38] at the domain boundary are effective for electromagnetic wave problems in the frequency domain. FLAME schemes described in Section 3 open promising new avenues in this research. These schemes combine geometric simplicity with accuracy of approximation: on the one hand, they can operate on simple Cartesian grids, and on the other hand, employ accurate approximating functions. As Section 3 shows, FLAME can model wave propagation and scattering very efficiently. For example, one standard benchmark problem – scattering from a dielectric cylinder – is solved on a Cartesian grid with a high level of accuracy, and the scheme is of order six, which is hardly possible with any other FD method. In this project, we propose to explore the applications of FLAME to the simulation of plasmon nanoparticles.

3. A NEW FINITE-DIFFERENCE CALCULUS Introduction
Flexible Local Approximation MEthods (FLAME) [44]−[48] is a new finite-difference calculus that seamlessly incorporates any desired local approximating functions into the difference scheme. In particular, FLAME approximates material interfaces algebraically, on geometrically nonconforming simple Cartesian grids, by a suitable set of basis functions, rather than by geometrically conforming meshes. For instance, in the vicinity of cylindrical or spherical dielectric particles, the field can be approximated by cylindrical or spherical harmonics. In many important cases, the

Proc. of SPIE 59270L-3

new approach achieves the same level of accuracy as the Finite Element Method (FEM) with two orders of magnitude fewer unknowns (see [44], [45] and examples below). The setup of FLAME can be outlined as follows. In the computational domain, generate a grid (usually regular Cartesian) and introduce a cover of that domain with subdomains. Each of the subdomains contains a grid stencil and carries a local approximation of the solution as a linear combination of suitable approximating functions {ψα} (e.g. plane waves or spherical harmonics). The approximating set may depend on the a priori knowledge of the local behavior of the solution in the respective subdomain (that is, over the respective grid stencil). In the so-called ‘Trefftz’ version of FLAME, the approximating functions are chosen as local solutions of the underlying differential equation. The key consideration is to provide as high approximation accuracy as possible. With a grid stencil and a set of approximating functions chosen, one generates the matrix comprising the nodal values of the basis functions. The difference scheme is in the null space of this matrix [44], [45]. The consistency error of the FLAME scheme is commensurate with the approximation error by the chosen basis set, and hence the accuracy can be dramatically higher than for standard schemes.

Test Example: Scattering from a Dielectric Cylinder
As an example, consider a typical 2D problem in the frequency domain for the H-mode (one-component H-field and a transverse electric field) ∇ · ε−1∇H + ω2µH = 0 (1) with the usual notation for the permittivity ε, permeability µ, field H and frequency ω. A similar equation is valid for the E-mode. Several types of boundary conditions are possible. To construct a Trefftz-FLAME scheme, consider the 9-point stencil of 3×3 neighboring nodes. In homogeneous regions, it is natural to choose eight basis functions as plane waves traveling from the central node in the direction of eight other nodes. At the exterior boundary of the domain, the stencil and the basis set are reduced: the stencil nodes outside the domain are ignored, and only outgoing waves are retained in the basis. This is a FLAME version of the PML [45]. In the vicinity of a dielectric scatterer, the Trefftz-FLAME scheme is constructed using cylindrical harmonics as basis functions [44]−[46]. In the vicinity of the scatterer the FLAME basis functions are cylindrical harmonics:
( [ Jn(kairr) + (aH / aJ ) · H n2 ) (kairr)] ejnφ, r > rcyl; Jn(kairr) / aJ · ejnφ, r ≤ rcyl where the coefficients aH, aJ are found from the boundary conditions on the cylinder. The basis set for the 9-point stencil contains 8 functions: the monopole, dipole and all but one quadrupole harmonics. This yields a 6th-order scheme [45].

1.E+00 1.E-01 1.E-02 1.E-03 errors 1.E-04 1.E-05 1.E-06 1.E-07 1.E-08 1.E-09
0.1 mesh size / radius 1

relative error norm without PML Ref. O(h^6)

Fig. 1. Relative error norms for the electric field vs. mesh size. (Mesh size normalized by the radius of the cylinder). Scattering from a dielectric cylinder. FLAME, 9-point scheme.

The numerical results are for the E-mode (one-component electric field along the axis of the cylinder and a twocomponent magnetic field). The 9-point stencil (3×3) is used throughout the domain. For testing and verification

Proc. of SPIE 59270L-4

purposes, the exact analytical value of the electric field on the outer boundary was imposed as a Dirichlet condition. Fig. 1 shows that convergence of the scheme is indeed of 6th order with respect to the mesh size; no standard FD method has comparable performance.

4. SIMULATION EXAMPLES Wave Propagation in a Photonic Crystal
Another application example of the same FLAME scheme is the Fujisawa & Koshiba photonic crystal bend [12], Fig. 2. A waveguide with a bend is formed by eliminating a few dielectric cylindrical rods from a 2D array. For FLAME simulations, the material characteristic of the rods was assumed linear, with the index of refraction equal to 3. For comparison, FE simulations (FEMLABTM) were run with several meshes: from 38,679 d.o.f. (i.e. unknowns) to 154,461 d.o.f., with 2nd order triangles in all cases.

000000000 000000000 000000000 • • •0000 0000O 000 000000000 00000 000

00000 000
Fig. 2. Left: the real part of the E field in the Fujisawa-Koshiba photonic bend (FEMLAB rendition). Electric field parallel to the rods. Center and right: two meshes with comparable accuracy for the Fujisawa-Koshiba photonic bend. An FE mesh with 77,104 second-order elements, 154,461 d.o.f. A 50×50 Trefftz-FLAME Cartesian grid (2,500 d.o.f.), 9-point scheme.

0000

©00

Convergence analysis shows that the FLAME solution on a regular 50×50 grid has the same level of accuracy as the FEM solution with 154,000 unknowns (Fig. 2). Note that for the 50×50 grid there are about 10.5 points per wavelength (ppw) in the air but only 3.5 ppw in the rods, and yet the FLAME results are very accurate because of the special approximation used.

Coupled Cylindrical Plasmon Particles
This set of simulations involves plasmon resonances of two coupled nanoparticles, with the setup proposed by Kottman & Martin [16]. In addition to being a very interesting physical case, this problem can be used for testing and verification of our simulation methods. Two perpendicular directions of propagation of the incident wave (cases (a) and (b) in Fig. 3) give rise to two different field distributions and resonance conditions. Fig. 3 demonstrates the excellent agreement between the FLAME and T-matrix results. The Johnson and Christy data [14] for gold and silver were used in all our simulations.

Plasmon-Enhanced Optical Tips
While FLAME schemes are quite effective in the situations described above, FEM is indispensable for complex geometries in a great variety of problems. Here FEM is employed for analysis and optimization of plasmon enhancement of optical tips. The ultimate practical goal is to achieve the highest and most robust field amplification around the apex of the plasmon-enhanced AFM tip in apertureless near-field optical microscopy.

Proc. of SPIE 59270L-5

The FEMLABTM software package is integrated with MatlabTM and gives the user access to the FE structures. FEMLAB commands can be incorporated into Matlab scripts, in particular for additional postprocessing and optimization. We have taken advantage of this capability and optimized the field enhancement with respect to geometric and physical parameters, for several setups.
The electric field along the central line betw een silver cylinders. Diameter 50 nm. Separation 10 nm.
T-m atrix: λ = 340nm

30 Electric field enhancement

E

0

25 20 15 10 5 0 -50 -30 -10 x [nm] 10

T-m atrix: λ = 350 nm T-m atrix: λ = 380 nm FLA M E: λ = 340nm FLA M E: λ = 350nm FLA M E: λ = 380nm

kinc

Case (a), Kottman & Martin [16]

30

50

The electric field along the central line between particles. Diameter 50 nm. Separation 10 nm. Wavelength 350 nm.

10
FLAME n=250 L=250

Electric field enhancement

E0 kinc Case (b), Kottman & Martin [16]

8 6 4 2 0 -50 -30 -10 y [nm]

T matrix method

10

30

50

Fig. 3. Two silver cylinders. Left: setup due to Kottman & Martin for two directions of illumination. Right: the magnitude of the electric field along the dashed line; comparison of FLAME and T-matrix results.

As an example, some representative results of such optimization are shown below. The setup, due to Martin et al. [24], includes a semi-spheroidal gold particle attached to a silicon AFM tip (Fig. 4, 5). The optimization parameters
TM

FEMLAB is a registered trademark of Comsol, Inc. (www.comsol.com). MATLAB is a registered trademark of The MathWorks, Inc.

Proc. of SPIE 59270L-6

were the axes of the particle and the wavelength. The objective was to maximize the average field over a reference disk 20 nm in diameter located 1 nm below the plasmon particle. Optimization was performed for several values of the permittivity of the tip. The height of the tip was fixed at 100 nm. Under these conditions, the optimal ratio of the axes of the spheroid equals 5.8 and the average field enhancement over the reference disk is about 90.
Filt

Edit opti995DHp9i9t69dfl/59bd9H J at! zrn
36

HI

—

A:0 flId IgyO (I _____________— j
40

J.i J I
Mt9:E0

y [nm]

HAT—ti 1]i
15039
170. 30)

0

0

4

B

ID

IGRID IEQUAL I

Fig. 4. A FEMLAB simulation of a conical AFM tip enhanced with a semi-spheroidal plasmon particle. Electrostatic approximation. Amplitude of electric field in space and the decay of field along the y-axis (upper plot; −10 nm < y < 0) and the x-axis (lower plot; 0 < x < 10 nm).
Gold semi-ellipsoidal plasmon particle (2b=20 nm) on a conical AFM tip (d=34 nm, h=300 nm) eps_core=15, a/b=3.2, θ=45º eps_core=10, a/b=3.2, θ=45º eps_core=15, a/b=5.8, θ=45º eps_core=15, a/b=5.8, θ=30º

AFM tip θ 2b θ h

160 140 120 | electric field | 100 80 60 40 20

plasmon particle d a

0 450

500

550 600 650 wavelength [nm]

700

750

Fig. 5. Plasmon-enhanced conical AFM tip. Incident field vertically polarized. Left: simulated geometry of the conical tip and gold plasmon particle. Right: comparison of the electric field amplitude in the electrostatic approximation for different parameters.

The level of plasmon enhancement depends strongly on the material of the tip. This is illustrated by Fig. 6, where the amplification of the field is plotted for different materials with fixed geometrical parameters of the model. Finally, the

Proc. of SPIE 59270L-7

dependence of the field enhancement on the ratio of the particle axes for a tungsten tip is shown in Fig. 7. The dielectric characteristics of tungsten measured by Ordal et al. [31] were used.
Cylindrical tip (height 100nm, radius 17nm) with semi-spheroidal apex (a = 40 nm, b = 8 nm). Not optimized.
140 Electric field enhancement 120 100 80 60 40 20 0 450 tip-tungsten + apex-Au tip-silicon + apex-Au tip-silicon + apex-Ag tip-tungsten + apex-Ag

500

550

600 wavelength [nm]

650

700

750

Fig. 6. Enhancement of cylindrical tip with semi-spheroidal plasmon particle for different materials as a function of the wavelength.

Cylindrical tip (height 100nm, radius 17nm) with semi-spheroidal apex (a=40nm). Wavelength 657 nm. Not optimized. 35 30 25 20 15 10 5 0 2 4 6 8 10 12 14 semiaxes ratio a/b 16 18 20 tip-tungsten + apex-Ag tip-tungsten + apex-Au

Fig. 7. Enhancement of cylindrical tungsten tip with semi-spheroidal plasmon particle vs. semi-axes ratio for gold and silver.

Electric field enhancement

5. CONCLUSION
Traditional and new computational techniques for frequency-domain simulations in photonics have been reviewed and applied to electromagnetic resonances in plasmon particles and to wave propagation in photonic crystals. The most

Proc. of SPIE 59270L-8

powerful of the conventional techniques is the FE method that is suitable for analysis and optimization of geometrically and physically complex structures such as plasmon-enhanced AFM tips. Optimization of field amplification has been carried out for a semi-spheroidal gold particle at the apex of the tip (the setup due to Martin et al.), with the maximum field enhancement around 100. The “Flexible Local Approximation” method (FLAME) provides seamless and simple integration of local analytical approximations of the solution into finite-difference schemes. When such local approximations are available (e.g. via plane waves or cylindrical harmonics), the method is both simple and very accurate. In this paper, FLAME has been applied to plasmon particles and to wave propagation in a photonic crystal bend. In all examples, the accuracy of the new method is dramatically higher than that of the standard FD or FE analysis. References
[1] P. Barber, C. Yeh, Scattering of electromagnetic waves by arbitrarily shaped dielectric bodies, Applied Optics, vol. 14, Issue 12, pp. 2864–2872, 1975. [2] Dirk Baumann, Christophe Fumeaux, Pascal Leuchtmann, and Rüdiger Vahldieck, Finite-volume time-domain (FVTD) modelling of a broadband double-ridged horn antenna, Int. J. Numer. Model. 17, pp. 285–298, 2004. [3] A. Bossavit, Computational Electromagnetism, Academic Press: San Diego, 1998. [4] S. Dey, R. Mittra, A conformal Finite Difference Time Domain technique for modeling cylindrical dielectric resonators, IEEE Trans. MTT, vol. 47, No. 9, pp. 1737–1739, 1999. [5] S. Dey, R. Mittra, and S. Chebolu, A technique for implementing the FDTD algorithm on a nonorthogonal grid, Microwave and Optical Technology Letters, Vol. 14, No. 4, pp. 213–215, 1997. [6] A. Doicu, Yu. A. Eremin, T. Wriedt, Convergence of the T-matrix method for light scattering from a particle on or near a surface, Optics Communications, vol. 159, pp. 266–277, 1999. [7] A. Ditkowski, K. Dridi, and J. S. Hesthaven, Convergent Cartesian grid methods for Maxwell’s equations in complex geometries, Journal of Computational Physics, 170, pp. 39–80, 2001 [8] B. Draine and P. Flatau, Discrete-dipole approximation for scattering calculations, J. Opt. Soc. Am. A, vol. 11, pp. 1491−1499, 1994. [9] B.T. Draine and Piotr J. Flatau, User Guide to the Discrete Dipole Approximation Code DDSCAT.6.0, 2003. http://arxiv.org/abs/astro-ph/0300969, http://www.astro.princeton.edu/~draine/DDSCAT.6.0.html [10] Nordin Félidj, Jean Aubard, and Georges Lévi, Discrete dipole approximation for ultraviolet–visible extinction spectra simulation of silver and gold colloids, Journal of Chemical Physics, vol. 111, No. 3, 1999 [11] Piotr J. Flatau, Improvements in the discrete-dipole approximation method of computing scattering and absorption, Optics Letters, vol. 22, No. 16, 1997. [12] T. Fujisawa, M. Koshiba, Time-domain beam propagation method for nonlinear optical propagation analysis and its application to photonic crystal circuits, J. of Lightwave Technology, vol. 22, No. 2, pp. 684–691, 2004. [13] Roger F. Harrington, Time-Harmonic Electromagnetic Fields, Wiley–IEEE Press; 2nd edition, 2001. [14] P. B. Johnson and R. W. Christy, Optical constants of the noble metals, Phys. Rev. B, vol. 6, pp. 4370–4379, 1972. [15] K. Lancea Kelly, Anne A. Lazarides, George C. Schatz, Computational electromagnetics of metal nanoparticles and their aggregates, Computing in Science and Eng., vol. 3, Issue: 4, pp. 67–73, 2001. [16] Jörg P. Kottmann, Olivier J. F. Martin, Retardation-induced plasmon resonances in coupled nanoparticles, Optics Letters, vol. 26, Issue 14, pp. 1096–1098, 2001. [17] Uwe Kreibig, Michael Vollmer, Optical Properties of Metal Clusters, Springer, 1995. [18] A. Lakhtakia, G. Mulholland, On two numerical techniques for light scattering by dielectric agglomerated structures, J. of Research of the Nat. Inst. of Standards and Tech., Vol. 98, No. 6, 699–716, 1993. http://fire.nist.gov/bfrlpubs/fire93/art023.html. [19] J. Lermé, B. Palpant, E. Cottancin, M. Pellarin, B. Prével, J. L. Vialle and M. Broyer, Quantum extension of Mie's theory in the dipolar approximation, Phys. Rev. B, vol. 60, pp. 16151–16156, 1999. [20] A. Liebsch, Surface-plasmon dispersion and size dependence of Mie resonance: Silver versus simple metals, Phys. Rev. B, vol. 48, pp. 11317–11328, 1993. [21] A. Liebsch, Surface plasmon dispersion of Ag, Physical Review Letters, vol. 71, Issue: 1, pp. 145–148, 1993. [22] Daniel W. Mackowski, Analysis of radiative scattering for multiple sphere configurations, Proc. Royal Soc. London A, vol. 433, pp. 599–614, 1991. [23] Michelle Duval Malinsky, K. Lance Kelly, George C. Schatz, and Richard P. Van Duyne, Nanosphere lithography: effect of substrate on the localized surface plasmon resonance spectrum of silver nanoparticles, J. Phys. Chem. B, vol. 105, pp. 23432350, 2001. [24] Yves C. Martin, Hendrik F. Hamann, and H. Kumar Wickramasinghe, Strength of the electric field in apertureless near-field optical microscopy, J. Appl. Phys., vol. 89, No. 10, pp. 5774–5778, 2001.

Proc. of SPIE 59270L-9

[25] M.I. Mishchenko, L.D. Travis, and D.W. Mackowski, T-matrix computations of light scattering by nonspherical particles: A review, J. Quant. Spectrosc. Radiat. Transfer, vol. 55, pp. 535–575, 1996. [26] M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles, Cambridge University Press, 2002. [27] D.W. Mackowski and M.I. Mishchenko, Calculation of the T matrix and the scattering matrix for ensembles of spheres, J. Optic. Soc. Amer. A, vol. 13, pp. 2266–2277, 1996. [28] M.I. Mishchenko and L.D. Travis, Capabilities and limitations of a current FORTRAN implementation of the T-matrix method for randomly oriented, rotationally symmetric scatterers. J. Quant. Spectrosc. Radiat. Transfer, vol. 60, pp. 309–324. 1998. [29] M.I. Mishchenko, G. Videen, V.A. Babenko, N.G. Khlebtsov, and T. Wriedt. T-matrix theory of electromagnetic scattering by particles and its applications: A comprehensive reference database. J. Quant. Spectrosc. Radiat. Transfer, vol. 88, pp. 357-406, 2004. doi:10.1016/j.jqsrt.2004.05.002. http://pubs.giss.nasa.gov/abstracts/2004/MishchenkoVideenBKW.html [30] Esteban Moreno, Daniel Erni, Christian Hafner, and Rüdiger Vahldieck, Multiple multipole method with automatic multipole setting applied to the simulation of surface plasmons in metallic nanostructures, Opt. Soc. Am. A, vol. 19, No. 1, pp. 101-111, 2002. [31] M. A. Ordal et al., Optical properties of the metals Al, Co, Cu, Au, Fe, Pb, Ni, Pd, Pt, Ag, Ti, and W in the infrared and far infrared, Applied Optics, vol. 22, No. 7, pp. 1099–1119, 1983. [32] B. Palpant, B. Prével, J. Lermé, E. Cottancin, M. Pellarin, M. Treilleux, A. Perez, J. L. Vialle, and M. Broyer, Optical properties of gold clusters in the size range 2-4 nm, Physical Review B, vol. 57, No. 3, pp. 1963–1970, 1998. [33] J. Peltoniemi, Electromagnetic scattering by irregular grains using variational volume integral equation method, J. Quant. Spectrosc. Radiat. Transfer, vol. 55, No. 5, pp. 637-647, 1996. [34] A.F. Peterson, S.L. Ray, R. Mittra, Computational Methods for Electromagnetics, Oxford University Press, 1998. [35] M. Quinten, Optical constants of gold and silver clusters in the spectral range between 1.5 eV and 4.5 eV, Zeitschrift für Physik D, vol. 101, pp. 211−217, 1996. [36] Michael Quinten, Optical effects associated with aggregates of clusters, Journal of Cluster Science, vol. 10, No. 2, pp. 319–358, 1999. [37] Kazuaki Sakoda, Optical Properties of Photonic Crystals, Springer-Verlag, 2001. [38] Z.S. Sacks, D.M. Kingsland, R. Lee, J.-F. Lee, A perfectly matched anysotropic absorber for use as an absorbing boundary condition, IEEE Trans. on Antennas and Propag. vol. 43, No.12, pp.1460-1463, 1993. [39] R. Schuhmann and T. Weiland, Recent advances in finite integration technique for high frequency applications, in: Springer Series: Mathematics in Industry, vol. 4, W.H.A. Schilders, Jan W. ter Maten, S.H.M.J. Houben, (Eds.), pp. 46–57, 2004. ISBN: 3-540-21372-4. [40] Julius A. Stratton, Electromagnetic Theory, McGraw-Hill Book Company, Inc.: New York, 1941. [41] Llorenç Serra, Angel Rubio, Optical response of Ag clusters, Zeitschrift für Physik D Atoms, Molecules and Clusters, vol. 40, Issue: 1-4, pp. 262 – 264, 1997. [42] K.-H. Su, Q.-H. Wei, X. Zhang, J. J. Mock, D. R. Smith, and S. Schultz, Interparticle coupling effects on plasmon resonances of nanogold particles, Nanoletters, vol. 3, No. 8, pp. 1087−1090, 2003. [43] A. Taflove, S.C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House Publishers, 2000. [44] Igor Tsukerman, Difference schemes with flexible local approximation, J. Comp. Phys., to appear. [45] Igor Tsukerman, Electromagnetic applications of a new finite-difference calculus, IEEE Trans. Magn., July 2005 (to appear). [46] Igor Tsukerman, (1) A new computational method for plasmon resonances of nanoparticles and for wave propagation. (2) Electrostatic finite-difference analysis without the ‘staircase’ effect. IEEE/ACES International Conference on Wireless Comm. and Applied Computational Electromagnetics, Honolulu, April 2005. [47] Igor Tsukerman, Flexible local approximation method for electro- and magnetostatics, IEEE Trans. Magn., vol. 40, No. 2, pp. 941– 944, 2004. [48] Achim Basermann, Igor Tsukerman, Parallel Generalized Finite Element Method for Magnetic Multiparticle Problems, in Springer Series: Lecture Notes in Computational Science and Engineering, LNCS 3402, M. Daydé et al. (Eds.), pp. 325–339, 2005. http://www.springerlink.com/index/10.1007/11403937_26. ISBN: 3-540-21372-4. [49] T. Wriedt and A. Doicu, T-Matrix Method for Light Scattering from a Particle on or Near an Infinite Surface, in: F. Moreno and F. González (Eds.), Springer Lecture Notes in Physics 534, pp. 113−132, Springer-Verlag, 2000. [50] Z.J. Wang, A.J. Przekwas, Yen Liu, A FV-TD electromagnetic solver using adaptive Cartesian grids, Computer Physics Communications, 148, pp. 17–29, 2002. [51] K.S. Yee, J.S. Chen, The finite difference time domain (FDTD) and the finite volume (FVTD) methods for solving Maxwell's equations, IEEE Trans. Antennas & Propag. vol. 45, No. 3, pp. 354–363, 1997. [52] W. Yu and R. Mittra, A conformal finite difference time domain technique for modeling curved dielectric surfaces, IEEE Microwave Wireless Comp. Lett., vol. 11, pp. 25–27, 2001. [53] Y.-l. Xu, Electromagnetic scattering by an aggregate of spheres, Appl. Opt., vol. 34, pp. 4573–4588, 1995. [54] I.A. Zagorodnov, R. Schuhmann and T. Weiland, A uniformly stable conformal FDTD-method in Cartesian grids, Int. J. Numer. Model. 16, pp. 127–141, 2003.

Proc. of SPIE 59270L-10