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.

 

Fig. 1. Geometry of the study for (a) the infinite circular cylinder along the z-axis and (b) the infinite square cylinder.

Download Full Size | PDF

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 (ε = ε ₀ ε_r ) which is a function of the spatial coordinates. In this paper, the surrounding medium is vacuum and therefore ε_r = m ². 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 Γ = Γ₀ + Γ₁ 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 ² (Ω) (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:

with ε₀ and μ ₀ 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.9t_c ;0.95t_c ] is generally adequate. Nevertheless, we have verified that for our study, for different values of Δx, Δt = 0.9t_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.7t_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:

 

Fig. 2. Comparison of the Real (a), Imaginary (b) part of the permittivity of gold, calculated with Drude’s and Drude-Lorentz’ models and the corresponding relative error with the experimental permittivity (Drude’s (c) and Drude-Lorentz’ (d) models).

Download Full Size | PDF

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.

Table 1. Computational procedure for the three methods.

View Table

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.

 

Fig. 3. Examples of non-Cartesian mesh used for FEM (a) and the Cartesian-grid used for FDTD (b). Zoom around the curve domain for FEM (c) and FDTD (d).

Download Full Size | PDF

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).

 

Fig. 4. Comparison of the spectra of intensities computed by Mie and FEM (a), Mie and FDTD with Drude’s model (b), Mie and FDTD with Drude-Lorentz’ model (c) and computed by FEM and FDTD with Drude-Lorentz’ model (d) for different distances d_y from the center of the particle along the y-axis for radius a = 15 nm. The vertical dashed line in (b), (c) and (d) shows the limit of FDTD computations.

Download Full Size | PDF

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⁶ nodes vs. 2.8 × 10⁴ 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).

 

Fig. 5. Maps of the total electric field intensity |E|² in the xy-plane for a gold cylinder of radius a = 15 nm for a p-polarized illumination at λ = 660 nm computed by the Mie’s theory (a), FEM (b), FDTD (c) and the absolute difference maps (Mie-FEM—) (d) and (Mie-FDTD) (e).

Download Full Size | PDF

In Fig. 5, we show the maps of the total field intensity |E|² 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 + j1.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⁴) 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).

 

Fig. 6. Comparison between (a) the intensity and (b) the relative errors as functions of the distance d_x from the center of the nano-object (a = 15 nm), along the x-axis computed by Mie, FEM and FDTD for different grid sizes.

Download Full Size | PDF

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)|² where the field is computed, the number of nodes is 4.0 × 10⁶ (around 2–3 hours CPU time) for FDTD and 2.8 × 10⁴ (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|² for the gold cylinder illuminated at λ = 550 nm and computed, with experimental values for the permitivitty (ε_r = -5.9473 + j2.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⁷ for FDTD (around 3–4 hours CPU time, Δx = Δy = 2 nm, Δρ_a ≈ 4π/λ) and 6.5 × 10⁴ 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.

 

Fig. 7. Maps of the total electric field intensity |E|² for a gold cylinder of radius a = 120 nm illuminated in a p-polarization source at λ = 550 nm computed by the Mie’s theory (a), FEM (b), FDTD (c) and the absolute difference maps (Me-FEM) (d) and (Mie-FDTD) (e).

Download Full Size | PDF

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⁵) 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⁷ – 10⁸ nodes vs. 6.5 – 10.1 × 10⁴ 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⁷ nodes, 9–10 hours and 6.5 × 10⁴, 22 hours (N × 15 minutes), respectively.

 

Fig. 8. Comparison of the spectrum of the total field intensities |E|² computed by FEM and FDTD for different distances d_y from the center of the particle along the y-axis for the half-length a = 15 nm.

Download Full Size | PDF

We show in Fig. 9 and Fig. 10 the maps of the total field intensity |E|² 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⁷ (4 hours CPU time, Δx = Δy = 0.25 nm, Δρ_a ≈ 4πa/λ, n ≈ 6.0 × 10⁴) for FDTD and 6.5 × 10⁴ (15 minutes CPU time) for FEM. For a = 120 nm (Fig. 10), the node numbers are 10⁸ (5 hours CPU time, Δx = Δy = 2.0 nm Δρ_a ≈ 4πa/λ, n ≈ 1.0 × 10⁵) for FDTD and 1.1 × 10⁵ (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.

 

Fig. 9. Maps of the total electric field intensity |E|² in the xy-plane of a gold nano-square of half-length a = 15 nm for a p-polarized illumination at λ = 660 nm computed by the FEM (a) and FDTD (b).

Download Full Size | PDF

 

Fig. 10. Maps of the total electric field intensity |E|² in the xy-plane of a gold nano-square of half-length a = 120 nm for a p-polarized illumination at λ = 550 nm computed by the FEM (a) and FDTD (b).

Download Full Size | PDF

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.