## Abstract

We compare the numerical results obtained by the Finite Element Method (FEM) and the Finite Difference Time Domain Method (FDTD) for near-field spectroscopic studies and intensity map computations. We evaluate their respective efficiencies and we show that an accurate description of the dispersion and of the geometry of the material must be included for a realistic modeling. In particular for the nano-objects, we show that a grid size around Δ*ρ*_{a}
≈ 4*πa*/*λ* (expressed in *λ* units) as well as a Drude-Lorentz’ model of dispersion for FDTD should be used in order to describe more accurately the confinement of the light around the nanostructures (i.e. the high gradients of the electromagnetic field) and to assure the convergence to the physical solution.

© 2005 Optical Society of America

## 1. Introduction

Since the development of near-field optical microscopy and its use to characterize the nano-structures, several numerical methods have been developed and used to simulate the interaction between the electromagnetic field and nanostructures. These methods were reviewed and compared in several papers [1, 2]. A direct comparison between the Finite-Difference Time-Domain (FDTD) method [3] and the Mie theory was also published by Challener et al [4]. In order to take into account a realistic description of the complex geometries and the materials, the Finite-Element Method (FEM) [5, 6, 7] has been introduced to compute near-field optic.

In this paper, we discuss the respective advantages of the FDTD and the FEM and we compare the results to the Mie’s theory in the case of infinite gold cylinders with different radius. We use a total/scattered field formulation and a combination of Drude and Lorentz dispersion models for our FDTD code [8] in contrast to Challener et al [4] who use a scattered field formulation. Finally, we present and compare numerical simulation results for FEM and FDTD on gold nano-square geometries where the geometrical singularities at the edges of the square generate a high gradient of the electric field.

The paper is organized as follows. Section 2 is devoted to the description of the analytical (Mie) and numerical methods (FEM, FDTD) used to compute the total electromagnetic field. In Sec. 3, numerical results obtained by Mie, FEM and FDTD are presented, discussed and compared, before concluding in Sec. 4.

## 2. The methods

In this section, we briefly introduce the analytical method (Mie) and the two numerical methods presented in this paper. The studied diffracting objects are depicted in Fig. 1. They consist of an infinite gold cylinder of radius *a* (see Fig. 1(a)) and an infinite square cylinder of half-length *a* (see Fig. 1(b)) illuminated by a p-polarized incident field with its wave vector **k** along the y-axis.

#### 2.1. Analytical model: Mie’s theory

The Mie’s theory gives, for a time-harmonic electromagnetic field, the analytical solution to the problem of the diffraction of a plane wave by a homogeneous sphere or an infinite cylinder. Although it was first published almost 100 years ago [9], it is still widely used and the subject of intensive research and improvements [10, 11, 12].

For the computation of the electromagnetic field around an infinite cylinder of radius *a*, we define *k* = 2*π*/*λ*, *ρ*_{a}
= *ka* and *ρ*_{r}
= *kr* where *r* is the distance between the point and the center of the cylinder, *ϕ* the polar angle and *m* is the ratio between the optical index of the cylinder and the optical index of the surrounding medium at wavelength *λ*. The total electromagnetic field **E**
_{t} = **E**
_{i} + **E**
_{s} is the sum of the incident field **E**
_{i} and the scattered field Ex defined, in p-polarization, by [13]:

with

where, $j=\sqrt{-1},{Z}_{n}\left({\rho}_{r}\right)={J}_{n}\left({\rho}_{r}\right)$ and *Z*_{n}
(*ρ*_{r}
) = ${H}_{n}^{\left(1\right)}$(*ρ*_{r}
) = *J*_{n}
(*ρ*_{r}
) + *jY*_{n}
(*ρ*_{r}
), for *l* = 1 and *l* = 3, respectively, *J*_{n}
and *Y*_{n}
are the Bessel functions and **e**
_{r} and **e**
_{ϕ} are the unitary vectors. The *D*_{n}
(*z*) are calculated by backward recursion through the following relation:

and the derivative *Z*′_{n} (*x*) is calculated using the well-known relationship for Bessel functions:

#### 2.2. The Finite-Element Method

In classical electromagnetics, the partial differential equations (PDE) are derived from Maxwell’s equations. The problem can be reduced to the wave equation, or Helmholtz’ equation if a time-harmonic dependance of the form exp(*jωt*) is assumed in the electromagnetic field [14]. In the case of an infinite cylinder along z, the fields have no variation with respect to the Cartesian coordinate (z). The geometry is defined as a function of (x,y)-coordinates and assuming that the illumination field is p-polarized, we compute the *z* component of the magnetic field *H*_{z}
and we deduce the electric field **E**(*x,y*) from Maxwell-Ampére’s equation. The z-component *H*_{z}
of the harmonic magnetic field in a domain Ω verifies:

where the different media are supposed to be non-magnetic, *c* is the velocity of light in vacuum, *ω* the frequency of the harmonic wave and *ε*_{r}
is the dimensionless relative complex permittivity of the considered medium (*ε* = *ε*
_{0}
*ε*_{r}
) which is a function of the spatial coordinates. In this paper, the surrounding medium is vacuum and therefore *ε*_{r}
= *m*
^{2}. Equation (8) is called the homogeneous scalar wave equation for the p-polarization.

To find the physical solution for the electromagnetic problem, a set of boundary conditions on the contour Γ = Γ_{0} + Γ_{1} of the computational domain Ω has to be established. The approximate boundary conditions have been extensively used in problems of wave propagation, radiation and guidance to simulate the material and geometric properties of surfaces [15]. In our approach, we consider the continuity of the tangential components of the electromagnetic field *H*_{z}
(impedance condition). The artificial outer boundary, which limits the region of computation, is characterized by the incoming wave amplitude and the first order Sommerfeld’s (or running waves) boundary conditions:

where *∂*/*∂n* is the outgoing normal derivative operator and *H*_{i}
the illumination field.

To solve the problem, the FEM makes use of a variational formulation (or weak formulation):

where *v* is a test function defined on *L*
^{2} (Ω) (the linear space of the scalar functions *v*, being 2-integrable on Ω). This integral formulation enables us to consider the domain of computations as the union of *N* sub-domains, thus $\Omega =\bigcup _{i=1}^{N}{\Omega}_{i}.$ The problem is solved on each sub-domain and the global solution is then the sum of the solutions in each element. The Eq. (10) is the projection of the solution on a basis of test functions *v*. The solution verifies exactly the PDE on each node for the given boundary conditions. The basis of polynomial functions *v* gives an approximation of the solution into the element. Numerical implementation involves a finite basis of functions and consequently, a linear system is solved by classical numerical LU decomposition or by a bi-conjugate gradient iterative procedure for large systems.

We have developed a home-made code with self-adaptive mesh refinement, which is particularly adapted to the description of complex geometry [6, 7].

#### 2.3. The Finite-Difference Time-Domain method

The Finite-Difference Time-Domain method is a direct space method, first introduced by Yee in 1966 [16], and now widely used for solving electromagnetic problems. Principles and improvements to these method have been published in several books [3, 17,18]. It mainly consists on the resolution of the discrete curl-Maxwell’s equations by iteration over the time.

With a spatial discretization Δ*x* and a temporal discretization Δ*t*, we obtain the standard recursion equations in a non-dispersive and non-magnetic medium for the electric and magnetic fields in a p-polarization:

$${E}_{x}{\mid}_{i,j}^{n+1}={E}_{x}{\mid}_{i,j}^{n}+\frac{\Delta t}{{\epsilon}_{0}{\epsilon}_{i,j}\Delta x}({H}_{z}{\mid}_{i,j}^{\frac{n+1}{2}}-{H}_{z}{\mid}_{i,j-1}^{\frac{n+1}{2}})$$

$${E}_{y}{\mid}_{i,j}^{n+1}={E}_{y}{\mid}_{i,j}^{n}+\frac{\Delta t}{{\epsilon}_{0}{\epsilon}_{i,j}\Delta x}({H}_{z}{\mid}_{i-1,j}^{\frac{n+1}{2}}-{H}_{z}{\mid}_{i,j}^{\frac{n+1}{2}})$$

with ε_{0} and *μ*
_{0} the permittivity and permeability of free space, and *ε*_{i}
the relative permittivity at position *x* = *i*Δ*x*, any field *U* at position *x* = *i*Δ*x* and for the instant *t* = *n*Δ*t* being represented by *U*${|}_{i}^{n}$
. The upper limit for Δ*t* is given by the Current time step ${t}_{c}=\frac{\Delta x}{\left(\sqrt{N}c\right)}$ with *N* the dimension of the problem. In our case we have ${t}_{c}=\frac{\Delta x}{\left(\sqrt{2}c\right)}$. In ref[4], it is stated that for metals or lossy dielectric, a value of Δ*t* ∈ [0.9*t*_{c}
;0.95*t*_{c}
] is generally adequate. Nevertheless, we have verified that for our study, for different values of Δ*x*, Δ*t* = 0.9*t*_{c}
was not sufficient to assure a rapid convergence of the algorithm. Using $\Delta t=\frac{{t}_{c}}{\sqrt{2}}=\frac{\Delta x}{\left(2c\right)}\approx 0.7{t}_{c}$ led to a faster convergence [19]. A way to avoid “stiff” problem would be to not use the recursive accumulator (at the cost of a highermemory usage), but rewrite the FDTD time propagation sheme by using a Runge-Kutta of high order with predictor-corrector (in order to avoid to developp a full implicit time propagation method).

FDTD is well adapted for spectroscopic studies [20], provided that an accurate analytical law of dispersion is used [8]. Indeed, results for a full spectrum can be obtained in a single run of the program. Typical laws used for FDTD simulations are Debye’s, Lorentz’, Drude’s or modified Debye’s dispersion models [20, 21, 22, 23, 24]. At least in principle, any dispersion law could be described in terms of a linear combination of Debye’s and Lorentz’ laws [3, 13]. For gold, Drude’s model only gives an acceptable description of the dispersion function for wavelength above 650 nm. For spectroscopic studies on gold nanostructures with wavelengths between 500 and 1000 nm [25, 26], an additional lorentzian term is needed in order to give a more accurate description of the dispersion function [8]. It can then be written as:

where *ω*_{D}
is the plasma frequency and *γ*_{D}
is the damping coefficient for Drude’s model, Ω_{L} and Γ_{L} stand for the oscillator strength and the spectral width of the Lorentz’ oscillators, respectively and Δ*ε* can be interpreted as a weighting factor. Equations (11) have then to be rewritten in order to take into account the dispersive properties. The two main methods used for
this purpose are the auxiliary differential equation method [27] and the recursive convolution method [17]. The details for the numerical implementation are given in references [3, 17] and reference [8] for the particular case of Drude-Lorentz.

We show in Fig. 2(a-b) the comparison between the real part (a) and the imaginary part (b) of the experimental permittivity and the theoretical dispersion laws calculated with Drude’s and Drude-Lorentz’ models for a gold particle. Figures 2(c-d) show the error of Drude’s (c) and Drude-Lorentz’ (d) models relative to the experimental permittivity as a function of the wavelength *λ* [28]. The relative error on the description of the imaginary part of the permittivities and on the real part for *λ* ≤ 600 nm is not negligible, especially when the Drude’s model is used and can thus induce errors in spectroscopic studies as we show in the next section.

## 3. Comparison of the methods

In this section, we compare results obtained by the different methods on simple examples as shown in Fig. 1. It consists in the computation of the total electromagnetic field around gold infinite circular cylinders and nano-squares (with radius or half-length *a* = 15 nm and *a* = 120 nm) illuminated by an incident p-polarized wave at various wavelength *λ* between 450 and 900 nm. In order to validate the methods, we also compute the spectra at short distances of detection *d*_{y}
(from 17 to 19 nm) from the center of the cylinder along the y-axis where the contribution of the scattered field is supposed to be important. It should be pointed out that the computational procedures are different for the three methods, and they are recapitulated in Table 1. The main advantage of FEM is the use of the experimental values of the relative permittivities of the materials (*ε*_{r}
) and a non-regular non-Cartesian adaptive mesh in order to reproduce de curvature of the object. Therefore, the grid refinement and the computation time are reduced, in contrast to standard FDTD which necessitate the use of a very fine regular or non-regular Cartesian-grid. Moreover, the use of Cartesian-grid in FDTD can produce artificial surface plasmon around the metallic structures for subwavelength cells. In order to reduce the number of cells and thus the computation time, non-regular, non-Cartesian adaptive meshing have also been developped for FDTD [3, 29] in particular cases but is not yet commonly and standardly used. Moreover, in this paper, we focus on the region close to the nano-object where high gradient of the electromagnetic field occurs. In this region, a quasi-regular meshing is necessary for both FEM and FDTD to describe the field variations (as illustrated in Fig. 3). Consequently, the use of a standard regular Cartesian grid for FDTD is not optimal but totally appropriate.

For the intensity map computation at a single wavelength *λ*, the FDTD method with standard PML [3] is applied using a incident sinusoidal plane wave and accurate Drude (or Drude-Lorentz) parameters for this wavelength. For the spectroscopic study, the procedure consists on a full computation on a very fine regular Cartesian-grid with a gaussian incident wave and an accurate long range dispersion law. The electromagnetic field is stored above the structure, time step after time step. Then we apply a Fast Fourier transform (FFT) on the result. This procedure contrasts to the FEM where the full computation is required for each wavelength *λ*. In this case, the use of the FDTD procedure can greatly reduce the computation time, especially when a large number of wavelengths is of interest. Figures 3 show examples of non-regular, non-Cartesian meshes used for FEM (a) and regular Cartesian-grid used for FDTD (b), as well as the ideal boundary of the object. In contrast to the Cartesian regular grid used for FDTD (d), the FEM grid is defined by the nodes on the boundary of the object (c). Consequently, the FEM grid matches closely the geometry of the object. The influence of the mesh on the computation of the fields is studied in the following subsections. We can note that a distinction must be made between: “convergence of the numerical method” and “convergence to the physical solution”. In this paper, we consider that the result of the computation has converged to the physical solution when the relative error between the result and the analytical solution (when it exists) is smaller than a tolerance (e.g. 1%). When no analytical solution exists, this relative error is computed by comparing the solutions obtained with two successive grid refinements.

#### 3.1. Case of the circular cylinder

We compare in Fig. 4 the intensity spectra computed by Mie-FEM (a) with experimental optical index $m=\sqrt{{\epsilon}_{k}},$, Mie-FDTD (b) with Drude’s dispersion law, Mie-FDTD (c) with Drude-Lorentz’ dispersion law and FEM-FDTD (d) for three distances *d*_{y}
= 17, 18 and 19 nm from the center of the particle of radius *a* = 15 nm along the y-axis where the contribution of the scattered field is supposed to be important. In order to avoid the discrepency between the experimental permittivity and the dispersion laws used for FDTD, the range of wavelength *λ* is restricted to [500 nm, 900 nm] and the Drude-Lorentz’ parameters used in the Eq. 12 are given in reference [8] (*ε*
_{∞} = 5.9673, *ω*_{D}
/2*π* = 2113.6 THz, *γ*_{D}
/2*π* = 15.92 THz Ω_{L}/2*π* = 650.07 THz, Ω_{L}/2*π* = 104.86 THz and Δ*ε* = 1.09).

As shown in Fig.4(a), Mie and FEM results are in good agreement for all the sprectrum range. The differences between Mie and FDTD results (see Fig.4(b–c)) are due to both the differences between experimental permittivities and the dispersion law used in the FDTD code, and by the grid refinement around the geometry. This example shows the influence of the dispersion law used for FDTD on the spectroscopic results. In the following, for the spectroscopic studies, the Drude-Lorentz’ model is used as the dispersion law. In contrast to spectroscopic procedure, if an intensity map for a single wavelength is of interest, an accurate Drude’ or Drude-Lorentz model can be used for this single wavelength (see Table 1). In the FDTD case, the number of nodes where the field is computed is large compared to the FEM case (around 4.0 × 10^{6} nodes vs. 2.8 × 10^{4} nodes, respectively). The CPU time on a PC for a full spectrum computation is 10 hours with FDTD. For FEM the CPU time is around 7–8 minutes for each wavelength and the total CPU time for one spectroscopic study is directly related to the number *N* of computed wavelengths (around 10–12 hours CPU time in the case of a wavelength step Δ*λ* = 5 nm).

In Fig. 5, we show the maps of the total field intensity |*E*|^{2} for the gold cylinder of radius *a* = 15 nm illuminated at *λ* = 660 nm for Mie (a) and the two numerical methods FEM (b) and FDTD (c). Mie and FEM computations are performed using the same non-Cartesian grid, which is adapted to describe the curvature of the cylinder and the experimental values for the permitivitty (*ε*_{r}
= -13.6969 + *j*1.0356) at this wavelength *λ*, while FDTD computations are performed using a regular Cartesian-grid (with Drude’s parameters: *ε*
_{∞} = 1.0, *ω*_{D}
/2*π* = 1968.43 THz, *γ*_{D}
/2*π* = 25.18 THz, Ω_{L}/2*π* = 0.0 THz, Ω_{L}/2*π* = 0.0 THz and Δ*ε* = 0.0). The Differences between Mie-FEM (with FEM mesh) and Mie-FDTD (with FDTD-grid) are also presented in Fig. 5(d) and (e), respectively. We can also remark the effect of the artificial surface plasmons, due to the use of a Cartesian-grid sheme in FDTD, which is revealed by higher magnitude of the difference for FDTD than for FEM in the vicitiny of the object. In order to achieve an accurate description of the cylinder geometry and to assure that the result converges to the physical solution (with an relative error ≈ 1%), very small Δ*x* and Δ*y* (Δ*x* = Δ*y* = 0.25 nm, Δ*ρ*_{a}
≈ 4*πa*/*λ*, Δ*t* ≈ 4.2 × 10^{-19} s,*n* ≈ 4.0 × 10^{4}) are necessary.

We compare in Fig. 6 the results obtained with FEM and FDTD for (a) the intensity and (b) the error relatively to Mie’s results, as a function of the distance *d*_{x}
from the center of the nano-object (*a* = 15 nm) along the x-axis for different grid sizes. We can remark, for a convergence of the numerical results to the physical solution (with an relative error around 1%), the FDTD requires a finer cell size than FEM, especially to describe accurately the high gradient close to the object (see Fig. 6(a) for Δ*x* = Δ*y* ≥ 0.5 nm).

It should be noted that in the case of dielectric objects, this kind of problem has been solved using the effective permittivities in order to smooth the transition between the inner and the outer parts of the structure [30, 31]. Nevertheless, in the case of dispersive materials, such a method would require the determination of all of the parameters of all the laws of dispersion involved for the description of the intermediate zone. Surface impedance boundary conditions (SIBC) could be a solution, but the numerical stability criteria for the time-step have still to be established [3]. For the map intensity |*E*(*x,y*)|^{2} where the field is computed, the number of nodes is 4.0 × 10^{6} (around 2–3 hours CPU time) for FDTD and 2.8 × 10^{4} (8 minutes CPU time) for FEM.

It is important to note that the time required for the FDTD computation does not only depend on the number of nodes but also on the way the computation is done. For a spectroscopic study, the computational process has to be continued until the full incident wave (gaussian shape) vanishes from the computation window, and its computational time depends on the spectral width of the wave. For the computation of an intensity map, the incident wave is a simple plane wave and a real-time Fourier transform is performed in order to compute the amplitude of the field in the computational window. In both cases, the time can strongly vary, depending on the resonant behavior of the diffracting structure [32]. Therefore the values of CPU time, given in this paper, are indicative values instead of absolute values. The computation time of FEM strongly depends on the method used to inverse the matrix obtained from the variational sheme. In this paper, we use an optimized LU decomposition for sparse matrix [33] or, for larger matrix, a iterative biconjugate stabilized gradient with hermitian preconditioning.

Now, we consider the infinite cylinder of radius axis *a* = 120 nm illuminated by an incident p-polarized field. In Fig. 7, we show the maps of the total field intensity |*E*|^{2} for the gold cylinder illuminated at *λ* = 550 nm and computed, with experimental values for the permitivitty (*ε*_{r}
= -5.9473 + *j*2.0904) and the related optical index $m=\sqrt{{\epsilon}_{k}},$, by Mie (a), FEM (b), FDTD (c) (with Drude’ parameters: *ε*
_{∞} = 1.0, *ω*_{D}
/2*π* = 1837.29 THz, *γ*_{D}
/2*π* = 104.16 THz, Ω_{L}/2*π* = 0.0 THz, Γ_{L}/2*π* = 0.0 THz and Δ*ε* = 0.0) and the differences between Mie-FEM (d) and Mie-FDTD (e). In this case, the quadrupolar term in Mie’s series is predominant and a high gradient of the electric field is expected at four positions around the cylinder. As in the previous simulation, we can remark the influence of the edges of the grid used in FDTD on the computational results, especially in the region of strong field confinement. This effect could decrease the accuracy of the computation of the field enhancement. The number of nodes where the field is computed is 10^{7} for FDTD (around 3–4 hours CPU time, Δ*x* = Δ*y* = 2 nm, Δ*ρ*_{a}
≈ 4*π*/*λ*) and 6.5 × 10^{4} with FEM (15 minutes CPU time). The dissymetry in Fig. 7(d) is related to the non regularity of the FEM mesh. The improvement of the mesh generation is a challenge in FEM community and is currently in progress.

#### 3.2. Case of the gold nano-square

In this last subsection, we consider the case of a square nano-object for which Yee’s cells seem to be better adapted. Strong confinement of the field is observed at the edges of the square section. Like in the previous examples and in order to assure a convergence of the numerical result to the physical solution (relative error around 1%), the mesh refinement is crucial and very small Δ*x*, Δ*y*, (Δ*ρ*_{a}
≈ 4*πa*/*λ*, *n* ≈ 1.3 × 10^{5}) are necessary for the FDTD code. Indeed, for greater values of Δ*x* and Δ*y*, the intensity at the corners is larger than the values obtained with FEM. The number of nodes for FDTD where the field must be computed, is larger than for FEM (around 10^{7} – 10^{8} nodes vs. 6.5 – 10.1 × 10^{4} nodes, respectively)

In Fig. 8 we show the comparison of the intensity spectrum computed by FEM and FDTD for different distances *d*_{y}
= 17, 18 and 19 nm from the center of the particle along the y-axis for the half-length *a* = 15 nm. The main differences in the spectra between FEM and FDTD are located for wavelengths below 550 nm, i.e. in the region where the relative error of the dispersion law used in FDTD for the real part of the complex permittivity *ε*_{r}
is the larger (see Fig. 2(d)). The number of nodes used for the computation and the CPU times for FDTD and FEM for a full spectrum are 10^{7} nodes, 9–10 hours and 6.5 × 10^{4}, 22 hours (*N* × 15 minutes), respectively.

We show in Fig. 9 and Fig. 10 the maps of the total field intensity |*E*|^{2} computed by FEM (a) and FDTD (b) for the gold nano-squares of half-length *a* = 15 nm illuminated at *λ* = 660 nm and of half-length *a* = 120 nm illuminated at *λ* = 550 nm, respectively. For *a* = 15 nm (Fig. 9), the node numbers are 10^{7} (4 hours CPU time, Δ*x* = Δ*y* = 0.25 nm, Δ*ρ*_{a}
≈ 4*πa*/*λ*, *n* ≈ 6.0 × 10^{4}) for FDTD and 6.5 × 10^{4} (15 minutes CPU time) for FEM. For *a* = 120 nm (Fig. 10), the node numbers are 10^{8} (5 hours CPU time, Δ*x* = Δ*y* = 2.0 nm Δ*ρ*_{a}
≈ 4*πa*/*λ*, *n* ≈ 1.0 × 10^{5}) for FDTD and 1.1 × 10^{5} (30 minutes CPU time) for FEM. The maximum values of the colorbars are limited in order to compare the global intensity around the nano-square geometry. The values of the peak intensities at the corners of the object are directly related to the grid size refinement used in the FDTD method and reveal the difficulty to describe a strong gradient of the field near a material geometric syngularity. The value of the intensity computed with FDTD is higher than with FEM due to the low convergence of the FDTD by decreasing the cell size.

## 4. Conclusion

In this paper, we have presented two numerical methods (FEM and FDTD) to compute the total electromagnetic field in context of near-field spectroscopy. The computation results are in good agreement with the analytical ones. We have shown that, in order to achieve realistic spectroscopic studies, the models must take into account both the geometric and dispersive characteristics of the complex materials. In particular, a sufficient grid size refinement must be considere in order to assure convergence and stabilization of the solution, particularely in the regions where a strong confinement of the light around the nanostructure occurs. We have shown that FEM is more adapted to describe accurately the complex geometries and to represent the maps of the fields. On the other hand, FDTD is more adapted to spectroscopic studies when a large number of wavelengths is of interest and/or when a non-regular and/or non-Cartesian mesh refinement is implemented.

## References and links

**
1
. **
D.
Barchiesi
,
C.
Girard
,
O.J.F.
Martin
,
D.
Van Labeke
, and
D.
Courjon
, “
Computing the optical near-field distributions around complex subwavelength surface structures: A comparative study of different methods
,”
Phys. Rev. E
**
54
**
,
4285
–
4292
(
1996
). [CrossRef]

**
2
. **
B.
Guizal
,
D.
Barchiesi
, and
D.
Felbacq
, “
Electromagnetic beam diffraction by a finite lamellar structure
,”
J. Opt. Soc. Am. A
**
20
**
,
2274
–
2280
(
2003
). [CrossRef]

**
3
. **
A.
Taflove
and
S.C.
Hagness
,
*
Computational Electrodynamics: The Finite-Difference Time-Domain Method
*
(
Artech House, Norwood
,
2000
).

**
4
. **
W.A.
Challener
,
I.K.
Sendur
, and
C.
Peng
, “
Scattered field formulation of finite difference time domain for a focused light beam in dense media with lossy material
,”
Opt. Express
**
11
**
,
3160
–
3170
(
2003
),
http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-23-3160
[CrossRef] [PubMed]

**
5
. **
R.
Fikri
,
D.
Barchiesi
,
F.
H’Dhili
,
R.
Bachelot
,
A.
Vial
, and
P.
Royer
, “
Modeling recent experiments of aperture-less near-field optical microscopy using 2D finite element method
,”
Opt. Commun.
**
221
**
,
13
–
22
(
2003
). [CrossRef]

**
6
. **
R.
Fikri
,
T.
Grosges
, and
D.
Barchiesi
, “
Apertureless scanning near-field optical microscopy : On the need of the tip vibration modelling
,”
Opt. Lett.
**
28
**
,
2147
–
2149
(
2003
). [CrossRef] [PubMed]

**
7
. **
R.
Fikri
,
T.
Grosges
, and
D.
Barchiesi
, “
Apertureless scanning near-field optical microscopy: Numerical modeling of the lock-in detection
,”
Opt. Commun.
**
232
**
,
15
–
23
(
2004
). [CrossRef]

**
8
. **
A.
Vial
,
A.S.
Grimault
,
D.
Macías
,
D.
Barchiesi
, and
M. Lamy de la
Chapelle
, “
Improved analytical fit of gold dispersion: application to the modelling of extinction spectra with the FDTD method
,”
Phys. Rev. B
**
71
**
,
085416
–
085422
(
2005
). [CrossRef]

**
9
. **
G.
Mie
, “
Beiträge zur Optik trìber Medien, speziell kolloidaler Metallösungen
,”
Ann. Phys.
**
25
**
,
377
–
445
(
1908
). [CrossRef]

**
10
. **
C.
Gréhan
,
G.
Gouesbet
, and
F.
Guilloteau
, “
Comparison of the diffraction theory and the generalized lorenz-mie theory for a sphere arbitrarily located into a laser beam
,”
Opt. Commun.
**
90
**
,
1
–
6
(
1992
). [CrossRef]

**
11
. **
H.
Du
, “
Mie-scattering calculation
,”
Appl. Opt.
**
43
**
,
1951
–
1956
(
2004
). [CrossRef] [PubMed]

**
12
. **
H.
Xu
, “
Calculation of the near field of aggregates of arbitrary spheres
,”
J. Opt. Soc. Am. A
**
21
**
,
804
–
809
(
2004
). [CrossRef]

**
13
. **
C.F.
Bohren
and
D.R.
Huffman
,
*
Absorption and scattering of light by small particles
*
(
John Wiley and Sons, New York
,
1983
).

**
14
. **
M.
Born
and
E.
Wolf
,
*
Principle of Optics
*
(
Pergamon Press, Oxford
,
1993
).

**
15
. **
J.
Jin
,
*
The Finite Element Method in Electromagnetics
*
(
John Wiley and Sons, New York
,
1993
).

**
16
. **
K.S.
Yee
, “
Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media
,”
IEEE Trans. Antennas Propag.
**
16
**
,
302
–
307
(
1966
).

**
17
. **
K.
Kunz
and
R.
Luebbers
,
*
The Finite Difference Time Domain Method for Electromagnetics
*
(
CRC Press
,
1993
).

**
18
. **
A.
Taflove
,
*
Advances in Computational Electrodynamics, the Finite-Difference Time-Domain Method
*
(
Artech House, Norwood
,
1998
).

**
19
. **
W.M.
Saj
, “
FDTD simulations of 2D plasmon waveguide on silver nanorods in hexagonal lattice
,”
Opt. Express
**
13
**
,
4818
–
4827
(
2005
),
http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-13-4818
[CrossRef] [PubMed]

**
20
. **
M.C.
Beard
and
C.A.
Schmuttenmaer
, “
Using the finite-difference time-domain pulse propagation method to simulate time-resolved the experiments
,”
J. Chem. Phys.
**
114
**
,
2903
–
2909
(
2001
). [CrossRef]

**
21
. **
F.L.
Teixeira
,
W.C.
Chew
,
M.
Straka
,
M.L.
Oristaglio
, and
T.
Wang
, “
Finite-difference time-domain simulation of ground penetrating radar on dispersive, inhomogeneous, and conductive soils
,”
IEEE Trans. Geosci. Remote Sens.
**
36
**
,
1928
–
1937
(
1998
). [CrossRef]

**
22
. **
S.K.
Gray
and
T.
Kupka
, “
Propagation of light in metallic nanowire arrays: Finite-difference time-domain studies of silver cylinders
,”
Phys. Rev. B
**
68
**
,
045415
–
045425
(
2003
). [CrossRef]

**
23
. **
M.
Futamata
,
Y.
Maruyama
, and
M.
Ishikawa
, “
Local electric field and scattering cross section of Ag nanoparticles under surface plasmon resonance by finite difference time domain method
,”
J. Phys. Chem. B
**
107
**
,
7607
–
7617
(
2003
). [CrossRef]

**
24
. **
J.T.
Krug
II
,
E.J.
Sanchez
, and
X.S.
Xie
, “
Design of near-field optical probes with optimal field enhancement by finite difference time domain electromagnetic simulation
,”
J. Chem. Phys.
**
116
**
,
10895
–
10901
(
2002
). [CrossRef]

**
25
. **
N.
Félidj
,
J.
Aubard
,
G.
Lévi
,
J.R.
Krenn
,
M.
Salerno
,
G.
Schider
,
B.
Lamprecht
,
A.
Leitner
, and
F.R.
Aussenegg
, “
Controlling the optical response of regular arrays of gold particles for surface-enhanced Raman scattering
,”
Phys. Rev. B
**
65
**
,
075419
–
075427
(
2002
). [CrossRef]

**
26
. **
J.
Grand
,
S.
Kostcheev
,
J.L.
Bijeon
,
M. Lamy de la
Chapelle
,
P.M.
Adam
,
A.
Rumyantseva
,
G.
Lérondel
, and
P.
Royer
, “
Optimization of SERS-active substrates for near-field raman spectroscopy
,”
Syn. Metals
**
139
**
,
621
–
624
(
2003
). [CrossRef]

**
27
. **
T.O.
Körner
and
W.
Fichtner
, “
Auxiliary differential equation: efficient implementation in the finite-difference time-domain method
,”
Opt. Lett.
**
22
**
,
1586
–
1588
(
1997
). [CrossRef]

**
28
. **
P.
Johnson
and
R.
Christy
, “
Optical constants of the noble metals
,”
Phys. Rev.
**
6
**
,
4370
–
4379
(
1972
).

**
29
. **
T.
Laroche
,
F.I.
Baida
, and
D.
Van Labeke
, “
Three-dimensional time-difference time-domain study of enhanced second harmonic generation at the end of a apertureless scanning near-field optical microscope metal tip
,”
J. Opt. Soc. Am. B
**
22
**
,
1045
–
1051
(
2005
). [CrossRef]

**
30
. **
S.
Dey
and
R.
Mittra
, “
A conformal finite-difference time-domain technique for modeling cylindrical dielectric resonators
,”
IEEE Trans. Microwave Theory Tech.
**
47
**
,
1737
–
1739
(
1999
). [CrossRef]

**
31
. **
W.H.
Yu
and
R.
Mittra
, “
A conformal finite difference time domain technique for modeling curved dielectric surfaces
,”
IEEE Microw. Wirel. Compon. Lett.
**
11
**
,
25
–
27
(
2001
). [CrossRef]

**
32
. **
C.
Ropers
,
D.J.
Park
,
G.
Stibenz
,
G.
Steinmeyer
,
J.
Kim
,
D.S.
Kim
, and
C.
Lienau
, “
Femtosecond Light Transmission and Subradiant Damping in Plasmonic Crystals
,”
Phys. Rev. Lett.
**
94
**
,
113901
–4 (
2005
). [CrossRef] [PubMed]

**
33
. **
T.A.
Davis
and
I.S.
Duff
, “
A combined unifrontal multifrontal method for unsymmetric sparse matrices
,”
ACM T. Math Software
**
25
**
,
1
–
20
(
1999
). [CrossRef]