Open Access
Add to BibSonomyAbstract
We present a novel numerical scheme for the simulation of the field enhancement by metal nano-particles in the time domain. The algorithm is based on a combination of the finite-difference time-domain method and the pseudo-spectral time-domain method for dispersive materials. The hybrid solver leads to an efficient subgridding algorithm that does not suffer from spurious field spikes as do FDTD schemes. Simulation of the field enhancement by gold particles shows the expected exponential field profile. The enhancement factors are computed for single particles and particle arrays. Due to the geometry conforming mesh the algorithm is stable for long integration times and thus suitable for the simulation of resonance phenomena in coupled nano-particle structures.
©2007 Optical Society of America
1. Introduction
It is well known that small metal particles exhibit resonant behavior at optical wavelengths. The large electromagnetic field enhancement associated with such plasmon resonances is of great importance for applications in nanoscience. In surface enhanced Raman scattering (SERS) [1, 2, 3] field enhancement plays a key role in the amplification of the Raman signal of molecules adsorbed on rough metal particles [4, 5, 6, 7]. In near field optical microscopy, field enhancement can provide a strong, well-defined and localized light source to investigate subwavelength structures. In addition, the investigation of the field enhancement by metallic nano-structures has attracted increasing interest in recent years for use in optical antennas [8]. In order to facilitate experimental studies, it is important to develop numerical tools that are able to predict enhancement factors for various designs and also optimize the numerical efficiency. The finite-difference time-domain method (FDTD) has primarily been used for this purpose [9]. It has been developed and tested extensively and proven to be a useful method for the simulation of metal nano-structures [8, 10]. However, because the FDTD method is carried out on a rectangular grid, structures with curved surfaces have to be approximated with stair-casing algorithms. In the case of field enhancement calculations, the existence of sharp grid corners leads to the development of spurious field spikes that travel around the curved metal surfaces. This problem cannot be removed by decreasing the grid spacing. In fact, a smaller grid spacing creates more field spikes because more sharp corners are present at the stair-case approximation of the surface. The corners causing spikes are a property of the rectangular grid and not of the spatial resolution. This leads to numerical errors in in the near electromagnetic fields which were investigated for conducting surfaces by [11]. As a result, the prediction of the field enhancement factors is only possible in early stages of an FDTD calculation and not possible at all for particle structures that are based on resonance phenomena and therefore require long runtimes. Alternative FDTD schemes make use of non-orthogonal grids in order to improve the surface approximation [12, 13]. Such algorithms suffer however from local reflections and higher global errors than conventional FDTD algorithms. An additional problem lies in the fact that metal nano-particle require a very high spatial resolution to guarantee sufficient accuracy and therefore impose stringent requirements on the simulation hardware.
We therefore propose a different approach to the study of curved metal surfaces in the time domain. We use a subgridding technique, that has been successfully used in the FDTD method [14]. Our algorithm combines the pseudo-spectral time-domain (PSTD) method for local mesh refinement with the FDTD method for the global primary grid. As described by [15, 16], the PSTD method is a suitable algorithm for the approximation of curved surfaces while maintaining high spatial convergence. In order to soften restrictions on the maximum stable time step imposed by the Courant limit, the spectral subgrid is split into multiple domains, which may be of non-rectangular shape. Each domain is projected onto a reference square by means of a conformal mapping technique. Due to the spectral convergence of the PSTD method, only a few grid points are required per domain to obtain high numerical accuracy. In the following we describe how the FDTD and PSTD method can be combined to yield an algorithm that takes advantage of the flexibility of the FDTD method for large simulations and the efficiency of the PSTD method for simulating cylindrical structures. We apply the method to the simulation of field enhancement on metal particles of various materials.
2. The hybrid numerical scheme
The FDTD method has been used extensively for the numerical study of electromagnetic problems. Therefore extensive research has been performed on the implementation of suitable current sources, absorbing boundary conditions and postprocessing options. Because the pseudospectral method has emerged well after the initial suggestion of the FDTD method, less research is available on the topic. However, in the last decade interest in the PSTD method has increased rapidly and the sophistication of the method has improved since its initial proposal. In this article we combine the advantages of both the FDTD and PSTD method to form an algorithm that is able to accurately model geometrical structures with curved surfaces. We apply the hybrid algorithm to the simulation of metal nano-particles. In order to simulate the dispersive properties of metals, we use the concept of auxiliary differential equations, which has been applied successfully in the FDTD algorithm [17]. The dielectric function of metals is approximated with a Drude-Lorentz model, which can be fitted accurately to experimental data. The dielectric function of a metal is given as a multi-pole approximation as follows
In above equation we have split the Drude pole with parameters for the plasma frequency ωP and collision frequency ν into a conductivity term and a Debye pole. The conductivity can be included in Maxwell’s equations without need for additional memory. For each of the N Lorentzian poles we have three parameters, the resonant frequency ωk, the damping constant Γk and a pole amplitude Ak. The dielectric permittivity at infinite frequency is denoted with ε∞. Each pole is linked to a macroscopic polarization in the time domain through an auxiliary differential equation. We then solve Maxwell’s equations with extension for the macroscopic polarization terms for the Drude pole PD and the Lorentzian poles Pk in the following form
(6)
The discretization of the above equations is carried out as suggested in earlier work [17]. The spatial discretization is performed with both the FDTD and the PSTD methods [18]. The above equations are advanced in time on the FDTD grid as described by Taflove [9]. Centered finite-differences yield second-order convergence in both space and time. In the PSTD method, however, the spatial derivatives are approximated in terms of Chebyshev polynomials. Due to their superior convergence properties this method is widely used for the solution of partial differential equations [19]. The grid points in this approach are defined on the Gauss-Lobatto grid, which corresponds to the zeros of the derivative of the Chebyshev polynomials [20]. Evaluation of the spatial derivatives is then achieved by performing a matrix multiplication with a differentiation operator, which is described in detail by [18]. Because all grid points are used for the evaluation of spatial derivatives the method converges exponentially. In order to lessen the computational burden of the matrix multiplications, the spectral domain is broken up into several subdomains on which above equations are updated individually. After each time step patching conditions as described in [15] are applied in order to make sure that information is passed correctly between neighboring domains. The PSTD method is applied to the local simulation of non-rectangular structures and the FDTD method is used to assemble the global solution. Each PSTD subgrid area is embedded in the FDTD parent grid as shown in Fig. 1. On the FDTD parent grid (light red), electric and magnetic fields are staggered in space in order to achieve second-order convergence. The PSTD method (light blue), however, uses collocated electromagnetic fields. Exchange of field information between parent and subgrid is therefore provided by interpolating field values that are placed on the subgrid boundary. Differing from the PSTD method described by Fan et al. [16], we use a staggered time-integration scheme in order to be consistent with the FDTD method, which implies that electric and magnetic fields are not collocated in time. Therefore we interpolate the magnetic fields on the boundary of the PSTD grid, using the magnetic fields of the parent grid, during the electric half time-step. For the electric step, we interpolate electric fields on the FDTD grid boundary using the fields from the PSTD subgrid. The interpolation of boundary values is carried out as described by Zakharian and coworkers in detail in [21]. We use linear interpolation tangential to the subgrid boundary and parabolic interpolation normal to the interface in order to obtain low reflections from the subgrid interface. One cycle of the algorithm is summarized as follows:
1. Update electric fields on the primary grid.
2. Interpolate magnetic fields of the primary grid to yield boundary values for the secondary grids.
3. Update electric fields on the secondary grids.
4. Replace electric boundary fields on the primary grid with averaged fields from the secondary grid.
5. Update magnetic fields on the primary grid.
6. Update magnetic fields on the secondary grids.
Fig. 1. Illustration of the positions of the primary FDTD and secondary PSTD grids for the TE case. The arrows indicate the fields that are used to calculate electric replacement fields for the primary grid or interpolate magnetic boundary fields for the secondary grid.
This cycle is applied until the final integration time is reached. The maximum stable time step that can be used for the integration is given by the Courant limit, which depends on the minimum grid size. Because the grid points for the PSTD method are not placed on a uniform grid, but are rather denser at the boundary of the domain, the time-step decreases quadratically with increasing number of grid points. However, due to the exponential convergence of the PSTD method the accuracy also increases with higher number of grid points. In our experiments we found that roughly 8 points per domain are needed in order to provide sufficient numerical accuracy and also in order to keep the interpolation error along the boundaries small. Therefore the further restriction of the stable time-step does not affect the global time-step significantly and runtimes similar to the FDTD method are observed.
3. Numerical study of metal nano-particles
We apply the hybrid algorithm to the simulation of metallic nano-particles of cylindrical shape. With the FDTD method local spikes emerge due to the discretization of the grid. These spikes occur along the stair-cased boundaries of curved surfaces. Even though the magnitude of the spikes can be very large, the fields stay localized and do not affect the far-field behavior of the particle. However, for the simulation of field enhancement, the far field is required and therefore the FDTD method is only suitable for limited integration times.
3.1. Evaluation of the algorithm
In order to evaluate the algorithm for metallic nano-particles, we compared our method to two commercial products, FIMMWAVE and OmniSim [22]. FIMMWAVE uses eigenmode expansion (EME) to calculate the spectral response of a given numerical setup. The second approach with OmniSim is based on the FDTD method. We consider transmission through an infinite line of metallic particles, which is simulated by applying periodic boundary conditions on the two sides of the computational domain, which are parallel to the direction of propagation, and absorbing boundary conditions in the direction of propagation. We consider a nano-particle of 100 nm diameter. The global grid size is set to 5 nm and fields are propagated for 60 fs. A time-step of 0.0095 fs is used to ensure stability of the algorithm. The subgrid area is composed of nine spectral elements, where 8 grid points are used in each direction of the subdomain. In Fig. 2 we show the calculated results for two material cases. As described in [17], we are able to fit the dielectric response of materials using a Drude-Lorentz model with very high accuracy. We therefore choose the cases of aluminium and nickel, for which material parameter are given in [17] in table 1. In this article the dielectric constant ε∞ is set to 1.0 to simplify the fitting process. As seen in Fig. 2, all methods agree well over a large wavelength range. In particular, we note that the resonance frequency of the particle is identified correctly with all methods. Deviations between our hybrid method and the FDTD are found for longer wavelengths for the simulation of aluminium, which is caused by stair-casing errors and the resulting spikes.
3.2. Field enhancement by a single particle
We now investigate the field enhancement from a single nano-particle. We apply our algorithm to the simulation of a metal cylinder comparing the Drude and Drude-Lorentz model for gold.
Table 1. Fitted material parameters for the the Drude-Lorentz model for Gold, the Drude model of Gold ([23]) and Silver.
Fig. 2. Transmission through an infinite chain of cylindrical nano-particles of 100 nm diameter. Shown are numerical calculations for both Nickel and Aluminium. The simulation was compared to the commercial products FIMMWAVE and OmniSim, using EME and FDTD, respectively.
The parameters for the Drude model were taken from [23] and the parameters for the Drude-Lorentz model were obtained using the fitting method described in [17]. We give all the material parameter values in table 1 as well as the Drude parameters for silver. In Fig. 3 we show the dielectric function of gold and the fitted material models. The Drude model matches the real part of the permittivity ε(λ) in the vicinity of the plasmon resonance of gold well but deviates strongly at longer wavelengths. Furthermore, the imaginary part of ε(λ) is not a good fit to the experimental data. The fitted Drude-Lorentzian model with 3 poles, however, describes both the real (3(a)) and the imaginary part (3(b)) of ε(λ) well for a wide wavelength range and shows a very good fit around 500 nm.
We consider a gold cylinder of 50 nm diameter. The cylinder is embedded in air with a refractive index of 1.0. Details of the geometry used in our calculations are shown in Fig. 4(a). The cylinder is approximated with a square PSTD subgrid with 80 nm long sides, composed of 13 spectral elements. We use eight grid points in each domain in order to achieve good accuracy. The FDTD grid is meshed with a grid spacing of 5 nm, thus the boundary of the subgridding area is composed of 16 FDTD grid points. The electromagnetic fields are integrated using a time-step of 9.5×10-3 fs to make sure that the algorithm is stable. The computational domain is terminated with uniaxial perfectly matched layers in order to absorb reflected fields [9]. The grid used in the calculation is shown in Fig. 4(b). The 13 elements are described in curvilinear coordinates in order to approximate the cylinder accurately. The thin lines in Fig. 4(b) show the location of the Gauss-Lobatto grid points. The simulation results for this setup are given in Fig. 5. We present the calculated enhancement factor for the Ex and the Ez components. The results are shown for calculations using the FDTD-only method (Fig. 5(a) and Fig. 5(c)) and for the hybrid approach (Fig. 5(b)) and Fig. 5(d))). The FDTD-only method is run with a grid resolution of 5 nm for the entire grid, which is equivalent to the global FDTD grid in the hybrid approach. The expected exponential decay from the cylinder boundary is clearly visible only in the hybrid result. The FDTD-only result, however, is greatly distorted by high field singularities along the cylinder boundary. These spikes stay localized but render calculations of field enhancement factors meaningless. As mentioned earlier in the article, the spikes are a result of the rectangular stair-casing approximation of the cylinder. Using a finer grid resolution does not remove the spikes because the corners of the stair-casing approximation will still be present. The variation of field enhancement with cylinder diameter is shown in Fig. 6. We illuminate the gold cylinder with a continuous wave source oscillating at a wavelength of 500 nm, which is close to the plasmon resonance frequency of the cylinder. The excitation source is applied with a half-Gaussian ramp-up, which results in a broadened line-width of the illumination that also covers the plasmon resonance frequency of gold. At 500 nm we also find a very close fit of the Drude-Lorentz model to the experimental data, in particular in the imaginary part of the dielectric function. We measure the enhancement factor for the Ex field component, which is vertical to the propagation direction, and also for the Ez component, which lies in the direction of propagation. The enhancement factor is determined after a total integration time of 3.2 fs and given as the ratio of the local field relative to the reference field. Therefore we need to run two calculations, one without the metal particle in order to obtain the reference field and a second calculation with the particle included. As shown in Fig. 6, the enhancement factor decreases almost linearly with increasing cylinder diameter. For comparison, we also computed the enhancement factor for gold with the dielectric function now modelled with a Drude-Lorentzian model as suggested in [17]. The Drude-Lorentzian model is a better fit to the experimentally determined permittivity function of gold. As shown in Fig. 6, a linear relationship between particle diameter and enhancement factor is also observed. The factor is however lower, because the material loss in the Drude-Lorentz model is higher than in the Drude model as presented in Fig. 3. The Drude model proposed by [23] approximates the real part of the permeability well in the vicinity of the plasmon resonance. However, the imaginary part is significantly lower than the experimental data by a factor of 9.5 and therefore material losses are greatly underestimated.
Fig. 3. The fitted dielectric function of Gold. Shown are the Drude-Lorentz model fit from OMNISIM [22] and the Drude-model fit from Maier [23].
Fig. 4. The geometrical details of the computational domain used for the simulation of cylindrical nano-particles. The spectral element mesh of the PSTD subgrid is shown in the right figure.
Fig. 5. Results from the simulation of a single gold particle embedded in air. Simulation results are shown for both the FDTD method and the new hybrid algorithm. The hybrid approach is virtually spike-free.
Fig. 6. Numerical investigation of the dependency of the field enhancement factor for the Ex component on the cylinder diameter. The factor shows an almost linear behavior.
We also show the enhancement factor for the materials silver and aluminium. Because both materials are modelled with a similar plasma frequency ωP as in the Drude-Lorentz model for gold, the enhancement factor is very similar to that of gold. The linear relationship between enhancement factor and particle diameter is also observed in this case.
3.3. Field enhancement of a particle array
In a more complex setup we consider field enhancement from an array of gold particles. The cylinders are assumed to be of equal diameter of 50 nm. We simulate six cylinders with an interparticle spacing of 100 nm. As in the previous simulation, we keep the time-step at 9.5×10-3 fs. The array is illuminated with a continuous wave source at 500 nm. The numerical results for a time integration time of 4 femto-seconds are shown in Fig. 7. We show the calculation of the field enhancement for the Ex field component. In Fig. 7 we show the magnitude of the electric field. Also presented are results from an equivalent FDTD calculation with a grid resolution of 5 nm. From the above figure it is evident that the FDTD-only result is greatly distorted at the material boundaries due to the sharp corners of the stair-casing approximation. Our hybrid algorithm however shows the smooth exponential increase of the electric fields towards the metal boundary. The hybrid approach is stable for long time integration and does not develop field singularities along the boundaries.
Fig. 7. The magnitude of the electric Ex fields after an integration time of 4 fs. Shown is the enhanced field for an array of six gold cylinders with an inter-particle spacing of 100 nm computed using the FDTD method (a) and the hybrid method (b).
3.4. Performance of the algorithm
In order to test the performance of the algorithm we revisit our first example. We calculate the transmission coefficient for a dielectric cylinder of refractive index 3.477 surrounded by air. The diameter of the cylinder is kept at 50 nm. By measuring the transmission coefficient the stair-casing error is not significant for the FDTD if the fields are sampled far enough from the cylinder. The numerical error of the method is measured as the root mean-square (RMS) error against a high-resolution converged simulation. The calculation parameters from section 3.1 are reused in this section. We compare the results of the hybrid method to the FDTD algorithm with a mesh refinement scheme, which is a feature available in the OMNISIM software. Here one level of refinement means, that the grid spacing in the subgrid area is half the grid spacing of the global grid.
The results of the performance measurements with respect to memory requirement and computation time are shown in table 2. We find that our hybrid method is significantly more efficient when high accuracy is required. For RMS errors below 0.1 % the hybrid method takes only an eighth of the time of the refined FDTD method. This is a result of the exponential convergence of the PSTD algorithm used on the subgrid area in the hybrid method. Therefore fewer grid points have to be used to reach a desired accuracy. This also results in lower memory requirements and therefore less memory throughput through the CPU. In order to achieve a similar numerical error with the FDTD, we are required to use a three-level mesh refined FDTD method. The equivalent FDTD subgrid approach requires storage for roughly twice as much memory as the hybrid approach. Therefore our algorithm shows significant improvements in both computation speed and memory allocation.
4. Conclusion
We have presented a hybrid FDTD-PSTD algorithm for the accurate simulation of field enhancement of metal nano-particles. Our algorithm does not show the field singularities that arise
Table 2. Comparison of the performance of the hybrid algorithm with the FDTD method and a refined FDTD method.
on stair-cased material boundaries with FDTD schemes. Through the introduction of a conformal mapping in the PSTD method, we are able to maintain smooth boundaries and therefore allow for the accurate simulation of field enhancement on metal particles. Because the method does not suffer from field singularities even for long time integration, it is well suited for the simulation of field enhancement in resonant particle structures.
Acknowledgement
W.H.P. Pernice would like to thank Photon Design for funding a DPhil research studentship.
References and links
1. H. Metiu, “Surface enhanced spectroscopy,” Prog. Surf. Sci. 17, 153–320 (1984). [CrossRef]
2. M. Moskovits, “Surface enhanced spectroscopy,” Rev. Mod. Phys. 57, 783–826 (1985). [CrossRef]
3. J. P. Kottmann, O. J. F. Martin, D. R. Smith, and S. Schultz, “Plasmon resonances of silver nanowires with a nonregular cross section,” Phys. Rev. B 64, 235,402 (2001). [CrossRef]
4. K. Kneip, Y. Wang, H. Kneip, L. T. Perelman, I. Itzkan, R. R. Dasari, and M. S. Feld, “Single molecule detection using surface-enhanced Raman scattering,” Phys. Rev. Lett. 78, 1667–1670 (1997). [CrossRef]
5. S. Nie and S. R. Emory, “Probing single molecules and single nanoparticles by surface-enhanced Raman scattering,” Science 275, 1102–1106 (1997). [CrossRef] [PubMed]
6. F. J. Garcia-Vidal and J. B. Pendry, “Collective theory for surface enhanced Raman scattering,” Phys. Rev. Lett. 77, 1163–1166 (1996). [CrossRef] [PubMed]
7. F. J. Garcia-Vidal, J. M. Pitarke, and J. B. Pendry, “Silver-filled carbon nanotubes used as spectroscopic enhancers,” Phys. Rev. B 58, 6783–6786 (1998). [CrossRef]
8. K. B. Crozier, A. Sundaramurthy, G. S. Kino, and C. F. Quate, “Optical antennas: Resonators for local field enhancement,” J. Appl. Phys. 94(7), 463242.
9. A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, Boston, MA, 1995).
10. J. T. Krug, E. J. Sanchez, and X. S. Xie, “Design of near-field optical probes with optimal field enhancement by finite difference time domain,” J. Chem. Phys. 116, 10,895 (2002). [CrossRef]
11. A. C. Cangellaris and D. B. Wright, “Analysis of the numerical error caused by the stair-steppedapproximation of a conducting boundary in FDTD simulations ofelectromagnetic phenomena,” IEEE Trans. Antennas Propag. 39(10), 1518–1525 (1991). [CrossRef]
12. R. Holland, “Finite Difference Solutions of Maxwell’s Equations in Generalized Nonorthogonal Coordinates,” IEEE Trans. Nucl. Sci. NS-30(6), 4589–4591 (1983). [CrossRef]
13. M. A. Fusco, M. V. Smith, and L. W. Gordon, “A Three-Dimensional FDTD Algorithm in Curvilinear Coordinates,” IEEE Trans. Antennas Propag. 39(10), 1463–1471 (1991). [CrossRef]
14. K. M. Krishnaiah and C. J. Railton, “A Stable Subgridding Algorithm and Its Application to Eigenvalue Problems,” IEEE Trans. Microwave Theory Techn. 47, 620–628 (1999). [CrossRef]
15. B. Yang, D. Gottlieb, and J. S. Hesthaven, “Spectral Simulations of ElectromagneticWave Scattering,” J. Comput. Phys. 134, 216–230 (1997). [CrossRef]
16. G.-X. Fan, Q. H. Liu, and J. S. Hesthaven, “Multidomain Pseudospectral Time-Domain Simulations of Scattering by Objects Buried in Lossy Media,” IEEE Trans. Geosci. Remote Sens. 40(6), 1366–1373 (2002).
17. W. Pernice, F. Payne, and D. Gallgher, “A general framework for the finite-difference time-domain simulation of real metals,” IEEE Trans. Antennas Propag.55 (2007). [CrossRef]
18. J. S. Hesthaven, P. G. Dinensen, and J. P. Lynov, “Spectral collocation time-domain modeling of diffractive optical elements,” J. Comput. Phys. 155, 287–306 (1999). [CrossRef]
19. C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods in Fluid Dynamics,, springer series in computational physics ed. (Springer-Verlag, New York, 1987).
20. J. Boyd, Chebyshev and Fourier spectral methods (Dover Publications, Inc, 2001).
21. A. R. Zakharian, M. Brio, and J. V. Moloney, “FDTD based second-order accurate local mesh refinement method for Maxwell’s equations in two space dimensions,” Comm. Math. Sci. 2(3), 497–513 (2004).
22. OmniSim and FIMMWAVE, Photon Design, 34 Leopold Street, Oxford OX4 1TW, UK.
23. S. A. Maier, P. G. Kik, and H. A. Atwater, “Optical pulse propagation in metal nanoparticle chain waveguides,” Phys. Rev. B 67(20), 205,402 (2003).
