Open Access
Add to BibSonomyAbstract
We demonstrate that simulating plasmonic nanostructures by means of curved elements (CEs) significantly increases the accuracy and computation speed not only in the linear but also in the nonlinear regime. We implemented CEs within the discontinuous Galerkin (DG) method and, as an example of a nonlinear effect, investigated second-harmonic generation (SHG) at a silver nanoparticle. The second-harmonic response of the material is simulated by an extended Lorentz model (ELM). In the linear regime the CEs are ≈ 9 times faster than ordinary elements for the same accuracy, provide a much better convergence and show fewer unphysical field artifacts. For DG-SHG calculations CEs are almost indispensable to obtain physically reasonable results at all. Additionally, their boundary approximation has to be of the same order as their polynomial degree to achieve artifact-free field distributions. In return, the use of such CEs with the DG method pays off more than evidently, since the additional computation time is only 1%.
©2011 Optical Society of America
1. Introduction
Plamonic nanostructures are well known for their striking optical properties allowing, e.g., for robust biolabels [1], local field scatterers [2] and surface-enhanced Raman scattering [3, 4]. Recently, also their nonlinear optical properties, such as four-wave mixing [5], second/third-harmonic generation (SHG/THG) [6], or interactions with molecules [7] have gained considerable interest. However, to obtain fundamental insight into these effects, a proper modeling is necessary that (i) works in the time domain, (ii) can properly describe the geometry, (iii) is reasonably fast and (iv) introduces no unphysical artifacts.
In the linear regime, analytical models exist for only a few basic plasmonic problems (e.g. for a sphere [8]) and, hence, in most cases numerical calculations have to be conducted to simulate realistic structures. Due to the small feature sizes, the dissipative properties of metals, and large field enhancements, many numerical problems arise in plasmonics; especially curved geometries such as spheres, nanorods or corners of larger objects are difficult to model. Recently, we showed for linear 2D scattering problems that curved elements (CEs) implemented within the discontinuous Galerkin (DG) method can master these difficulties and improve nano-optical simulations significantly [9]. Niegemann et al. further supported this statement by investigations on a 3D spherical cavity [10]. Nevertheless, the influence of CEs on linear scattering problems and on nonlinear optical calculations has not been reported so far.
In this paper we investigate how CEs implemented within the DG method can be used to simulate linear and nonlinear plasmonic problems. As a realistic test case we describe scattering and second-harmonic generation (SHG) from a spherical silver nanoparticle of 80 nm diameter (see Fig. 1).
Fig. 1 Schematic drawing of the numerical setup used for calculating the SHG response from a 3D sphere: The scattered field / total field (SF/TF) contour is used to inject a short laser pulse from the left exciting the 80-nm Ag nanoparticle. The resulting scattered light is then monitored and absorbed by a perfectly matched layer (PML, parameters as in [9]).
The outline of this paper is as follows: Firstly, the basic features of the DG method will be recalled, and higher-order boundary approximations of a sphere using CEs will be introduced. Secondly, the scattering cross section and the linear near field of the spherical particle will be calculated in order to determine the performance of CEs for various orders of boundary and polynomial approximations in the linear regime. Thirdly, an advanced material model will be introduced and the SHG of the sphere will be calculated. The resulting near-field distributions are investigated with respect to convergence and compared to the linear results with regard to accuracy. Finally, the impact of the nonlinearities on the numerical effort will be discussed.
2. Discontinuous Galerkin method and curved elements
The discontinuous Galerkin (DG) method is a relatively young approach in the field of optics. It was first proposed by Hesthaven and Warburton [11], then extended by Lu to meet nano-optical requirements [12] and recently applied to plasmonic systems [13, 14]. Nevertheless, it already outperforms finite-difference time-domain calculations for some complex systems [15]. A detailed coverage of DG can be found in [16]. The most relevant features of the DG method are summarized in the following.
The DG method is a volume method, i.e., the problem is solved by discretizing the relevant space into basic elements in which the field can be easily expanded into higher-order basis functions such as polynomials. Generally these elements can be of any arbitrary shape, which allows for describing the problem in a body-conform way, i.e., approximating any geometry to the best fit with these elements. However, triangles (2D) or tetrahedrons (3D) are used most commonly as elements, similarly as in the finite-element method (FEM). The unique feature of DG is that the electric field no longer needs to be continuous across the boundary between two neighboring cells. In order to still achieve a coupling between cells, a numerical flux is introduced as a penalty term on the basis of a conservation law, i.e., the net flux flowing into or out of each element through its sides/faces, is minimized. The numerical flux is obtained by integration over the different sides/faces of the spatial elements and is one scalar number per coordinate. Therefore, in contrast to FEM, higher-order basis functions (cf. Fig. 2) can be used without increasing the size of the coupling matrix. Finally, we end up with an ordinary time-dependent differential equation that can be solved through integration by any standard method available, for example a Runge-Kutta method [17] as used here.
Fig. 2 Polynomial approximation of the field representation for a 3D element. As a visualization of the polynomials the distribution of their nodal points is sketched for the 1st, 2nd and 4th orders (P).
Curved elements (CEs) are triangles/tetrahedrons that possess (at least) one side/face which is no longer straight but deformed (curved) in order to provide an optimal match with the surface geometry of the nanoobject under consideration (see Fig. 3). In 2D the construction of a curved element is straight-forward: as the triangle takes the form of a pie slice, the surface nodes are just displaced and the volume nodes are mapped uniformly [9]. In 3D the situation becomes more difficult as for a curvature in two directions no perfect mapping to a flat surface exists. Hence, the integration over the curved edge/surface (needed for determining the numerical flux) becomes inexact. To still achieve the same accuracy as with linear elements, more surface nodes are needed and the integration has to be checked for stability. We found that a linear projection is stable and accurate enough for our purpose.
Fig. 3 Schematics of the nodal distribution for linear and curved elements: (a) 2D linear, (b) 2D curved, (c) 3D linear and (d) 3D curved elements. The blue dots represent the (surface) nodes.
For discretizing the volume into basic elements a mesh generator can be used, which easily provides meshes for a large class of problems. For example, we used the mesh generator ‘Net-Gen” [18]. When curved elements are included into these meshes, the object boundary needs to be parameterized using polynomials in order to reflect its curved nature. The order B of these polynomials is visualized in Fig. 4 and is crucial for the accuracy of the calculations. In our investigations we considered 1st-, 2nd- and 4th-order boundary approximations.
Fig. 4 Geometric description of a 3D sphere using a 4th-order icosahedron and a 1st-, 2nd-and 4th-order boundary approximation (B). Edges of the elements are visible for both B=1 (linear elements) and for B=2 (CEs), but not for CEs with a 4th-order approximation (B=4).
3. Linear calculations
As a test system we use a silver nanoparticle of 80 nm diameter, because it represents the typical dimension of a plasmonic system and can be described analytically in terms of multipole expansions and Bessel functions in the linear regime [8]. For the sake of simplicity, we applied the method of multiple multipoles (MMP) [19] to determine the expansion coefficients, which served as a reference solution. The 80-nm silver sphere is treated as a Drude metal (ωp = 1.30 × 1016 s −1, γ = 3.23 × 1013 s −1), embedded in air (cf. Fig. 1) and excited with a differentiated Gaussian pulse (see [14], center wavelength 250 nm) in order to obtain the whole spectral response within a single calculation. For analyzing the influence of CEs in 3D scattering problems, 4 adaptive meshes (M=1...4) (see Fig. 5) having different mesh sizes, 2 polynomial orders (P=2,4) and 3 boundary representation orders (B=1,2,4) were prepared (see Figs. 2 and 4, respectively). In order to investigate the performance of these configurations, we calculated both the maximal near-field enhancement and the scattering cross section.
Fig. 5 Visualization of the meshes used (M=1...4) by means of the grid on the sphere. The four meshes consist of 6675, 7803, 10154 and 33226 elements, respectively, having 176, 385, 780 and 2278 boundary elements.
The results are depicted for linear and curved elements in Fig. 6 (B=1,2, P=2). As expected, the performance/convergence of the linear elements worked out to be very low. For the scattering cross section the coarsest mesh (M1) does not show any coincidence with the MMP reference. With finer meshes the situation improves. However, even for the finest mesh (M4) only two of the three resonances are reproduced.
Fig. 6 Maximum field enhancement and scattering cross section (normalized to the geometrical cross section of the particle) for the different meshes (M=1...4) as well as for linear and curved elements (B=1,2, P=2).
In contrast, with CEs the outcome is much better: Although the M1 mesh shows obvious deviations, already M2 matches the reference spectrum. The difference in performance is even larger for the near-field enhancement: with linear elements there is no convergence to the reference solution at all (with the meshes used), whereas CEs show a good agreement with the reference already for M2 and negligible deviations for the finer meshes. Hence, a M2 mesh with CEs is superior to the best linear mesh.
The reason for this is that with linear elements the smooth surface of the scatterer is replaced by an angular shape, which causes edge plasmons (visible in in the field plots further down) to be excited and leads to an overestimation of the field enhancement in the proximity of the sphere and in an underestimation of the scattering. Therefore, similar to the 2D case [9], much finer meshes are needed to achieve the same accuracy with linear elements as with CEs.
Now we will investigate the importance of the boundary approximation. Again, the scattering cross section and the field enhancement are plotted for curved elements, but this time with the polynomial order P=4 (see Fig. 7). The results are in general quite stable independently of the order of boundary approximation. Only at the lowest resonance (at 222 nm – an octopole mode) and with the coarsest mesh with B=2 the field enhancement deviates visible from the reference solution. For all other cases the curves look almost identical; only a closer look at the relative error shows small deviations and reveals the following dependences: (i) with finer meshes the error decreases, (ii) the near-field error is much larger than the error of the scattering cross section – this is also true for the reference MMP solution – and (iii) B=4 gives slightly better results than B=2. However, as already M2 gives a satisfying accuracy, the boundary approximation seems to play only a minor role here.
Fig. 7 Maximum field enhancement and normalized scattering cross section for the two coarsest meshes, 4th-order polynomial degree (P=4) and 2nd- and 4th-order of boundary approximation. In the lower two plots also the error of the MMP reference is shown, (left) averaged relative error and (right) maximal relative error. Note, for (M=2,B=4) or (M=2, B=4) the geometrical description is rather good and other effects such as reflections from the PML or the detailed discretization of the mesh become dominant. They cause the spectral variations of the error profile.
To finish of our linear-regime investigations, we will take a closer look at the near-field distribution of the sharp resonance at 251 nm. This quadrupole mode is depicted in Fig. 8 for one linear- and two curved-element setups. While for the linear case the finest mesh M4 was used, the curved cases have a higher polynomial order but only M2 was applied. In general, the field away from the particle surface looks similar for all three plots; however, for the linear elements one can clearly see artifacts stemming from the discretization. These artifacts are the edge plasmons and disappear when curved elements are used. In contrast, the two CE plots show no artifacts and look identical even though the boundary approximations are different.
Fig. 8 Linear near-field distribution of the Ag sphere after excitation by a plane wave from the left. The averaged E field at a wavelength of 251 nm is plotted for three different configurations as indicated. All plots are scaled equally and the color range is: white/yellow – high field strength; black – low field strength.
This shows that CEs are not only more robust and converge faster than linear elements with regard to spectral features, but also are less prone to unphysical artifacts in the near-field distribution.
4. SHG calculation
After having verified and characterized our system in the linear optical regime, we will now perform nonlinear calculations for the silver sphere exemplified by second-harmonic generation. In the frequency domain, SHG calculations are often carried out via a first-order perturbation-theoretical ansatz, i.e., the linear field distribution is first calculated by analytical or numerical methods (see e.g. [20, 21]) and then second-order susceptibilities are used to calculate the second-harmonic polarization. The susceptibilities are free parameters which are commonly determined by adapting the numerical far-field results to experimental observations [21].
In contrast, for DG – like for all time-domain methods – SHG effects can be directly incorporated by using a material model that supports higher-harmonic responses. Hence, the back-action of the second-harmonic field on the linear field is automatically included, susceptibilities are replaced by more fundamental parameters and time-dependent parameters such as the shape of the pulse can be investigated. Various material models exist to date with different strengths and weaknesses [22–24]. A detailed comparison is out of scope of this paper, as we are only interested here in the impact of curved elements on nonlinear calculations, which should be similar for all models. We have chosen an extended Lorentz model (ELM) [24] for describing the linear and nonlinear responses, because it is simple to implement, has a relatively low computational burden and is numerically robust. Furthermore, its simplicity will help us to understand the crucial parts of SHG calculations.
The ELM is based on the Drude-Lorentz model as commonly used for describing the linear response of metals in time domain calculations. To derive the ELM, one has to write down the general response of the polarization to the electromagnetic field [24]
with E being the electric field, H the magnetic field, Nd the dipole density, q the charge of an electron, m the effective electron mass, μ the permeability, γ the damping coefficient and ω 0 the dipole eigenfrequency. Now, E and P × H can be approximated with a 2nd-order Taylor expansion around P = 0 (note that the standard Drude-Lorentz model only uses a 1st-order expansion): By inserting these expressions in Eq. (1) we obtain the ELM: with being the plasma frequency. Note that no other free parameters are needed besides the standard Drude-Lorentz parameters in order to obtain a nonlinear response. For our SHG simulations, Eq. (4) is incorporated in DG as an auxiliary differential equation which is updated simultaneously with the electric field computations, i.e., the polarization at each nodal point of the mesh is updated at each time step. Similarly to the linear case, the sphere was illuminated from the left, but now with a Gaussian pulse having its center wavelength at 400 nm and a duration of 10 fs.First, the convergence was investigated for different meshes (M), boundary approximations (B) and polynomial degrees (P) with respect to the maximal near-field enhancement. As no analytical solution exists for this SHG scattering problem, the combination (M=3, B=4, P=4) was used as a reference and the error (deviation from the reference) was determined for the fundamental (400 nm) and the second-harmonic (200 nm) frequency. The results of the CE calculations are plotted in Fig. 9 – the results for linear elements are omitted, as they have too large errors (the SHG scattering is predicted 100% too high).
Fig. 9 Influence of the parameter set (mesh M, boundary approximation B, polynomial degree P) on the accuracy of the maximal near-field enhancement. The errors for the fundamental and second-harmonic frequency are shown for CEs and the different sets.
For the fundamental frequency the maximal near-field enhancement converges quite rapidly both with finer meshes and with increasing orders B and P. In the best case (M=2, B=4, P=4) the error reaches 0.06%, which is comparable to the linear near-field errors (cf. Fig. 7). However, for the SHG results the situation looks quite different: although the error seems to decrease with finer meshes for B=P=2, it explodes when the polynomial degree is increased to P=4. Only after the boundary approximation has been adjusted, does the error approach values around 0.6% (M=2, B=4, P=4). This is similar to the near-field error at the fundamental frequency and is a clear success. But, why do we obtain such huge deviations for (B=2, P=4)?
To shed some light on this problematics, we need to look at the SHG near-field distribution. In Fig. 10 the electric field is plotted for linear elements (B=1), CEs having a 2nd-order boundary approximation (B=2) as well as CEs having a 4th-order boundary approximation (B=4). Although the linear case uses the finest mesh, the field is overestimated by a factor of 100 and the plot is not acceptable – it looks more like the corona of the sun than a SHG pattern of a silver sphere. For B=2 the general field pattern looks better (and is comparable to perturbation-theoretical calculations – see [25]) but still shows unphysical artifacts at the edges of the CEs. Only for B=4 all artifacts disappear. These observations are significantly different than in the linear regime (cf. Fig. 8) where already for B=2 all visible glitches disappeared.
Fig. 10 SHG field of a Ag sphere illuminated with a femtosecond pulse from the left. The parameters M, B, P are the same as in Fig. 8 and all plots are scaled equally.
The origin of this can be understood from Eq. (2) in which the electric field was expanded in a second-order Taylor series. In the series the first term reflects the linear response and, hence, is only proportional to E 0; the second term provides the nonlinear response and is proportional to the gradient of E 0 times the polarization. So, in contrast to the linear case, the SHG calculations are not only dependent on the field in the previous time step but also on the gradient of that previous field. Therefore, an as-good-as-possible geometrical description gets even more important.
Besides by the mesh, the geometrical accuracy is determined by the boundary approximation order B as well as the polynomial order P. Increasing B leads to a better boundary approximation of the CEs and increasing P to a better field approximation within each element. Since DG shows its strengths when higher-order basis functions are used, one is eager to use as large Ps as possible. However, this can result in difficulties as shown in Fig. 11: When the polynomial representation is finer than the boundary approximation, some of the nodal boundary points do not lie on the geometrical boundary. This results in fields and field gradients located at slightly wrong positions. This fact is certainly of minor importance in the linear regime, but for nonlinear calculations, which are much more field sensitive, it results in unphysical field artifacts and large errors. For SHG calculations this means that only parameter sets with P=B lead to reliable results. Note that Niegemann et al. chose B=2 and up to P=6 in [10] for their linear calculations.
Fig. 11 Visualization of the geometrical mismatch for different boundary approximations (B) and polynomial orders (P). Black line: boundary to be described; blue line: shape of the CE; red dots: nodes of the polynomial representation. Note, B < P leads to a wrong boundary description.
5. Numerical effort
Finally, we will discuss the performance of CEs for the SHG calculations. A direct comparison of linear and curved elements would be the easiest way. However, the results of the linear-element SHG calculations are too inaccurate to be used as a benchmark. Thus, we firstly will determine the influence of CEs in the 3D linear regime and then investigate the impact of the extended Lorentz model.
When curved instead of linear elements are used for the same mesh, the numerical cost increases because the computational effort per time step increases (CEs are much more complex) and the number of time steps grows too, as the minimal node distance and, hence, time-step width decreases. On the other hand, much coarser meshes are sufficient to obtain the same accuracy as with linear elements. Hence, we will compare the overall computational effort needed to gain the same (sufficient) accuracy in the spectral calculations of Fig. 6. The normalized computation times of these calculations are as follows:
| mesh | M1 | M2 | M3 | M4 |
|---|---|---|---|---|
| linear elements | 0.25 | 0.57 | 1.04 | 8.75 |
| curved elements | 0.85 | 1 | 1.69 | 10.20 |
The gap between CEs and linear elements is largest for M1, but it decreases with increasing mesh number (with finer meshes). This is caused by the smaller curvature of CEs for finer meshes which leads to a larger possible time-step width. For comparing the performance at similar accuracies, we choose the CE-M2 and the linear-element-M4 cases and obtain a speed increase by a factor of 8.75 for CEs. Even the finer CEs-M3 mesh still leads to a performance gain of 5.15, which proves the advantage of CEs for 3D plasmonic calculations in the linear regime.
Finally, we will determine the impact of the extended Lorentz model on the computation speed. As the additional differential equations needed for the ELM are only slightly more complex than for the Drude-Lorentz model, we expect small changes. Indeed, we found an increase of only around 1%. This means that 3D plasmonic calculations using CEs and including SHG are actually faster than calculations on the same system in the linear regime using linear elements.
6. Conclusions and outlook
In summary, we successfully implemented CEs in 3D-DG and investigated the linear and nonlinear response of a Ag sphere. We calculated far-field scattering cross sections, the maximal near-field enhancements and near-field plots. In the linear regime we found that CEs with a 2nd-order boundary approximation are sufficient. Compared to linear elements they show a 8.75-times faster computation for the same accuracy, a much faster convergence and fewer un-physical artifacts in the near field. For our DG-SHG calculations, CEs are required to achieve any physically meaningful results at all. The material model used led to an increase of the computation time by just 1% and the maximal near-field error was around 0.6%. Furthermore, it was found that the boundary approximations need to have the same order as the polynomials of the basis functions to provide artifact-free results.
In general, the field patterns of our SHG simulations are comparable to those found in perturbation-theoretical calculations [25] and the model can easily be extended e.g. to obtain spectrally resolved polar scattering distributions, to investigate the influence of the chirp of the excitation pulse or to optimize the SHG scattering via pulse-shaping techniques similar to [26]. Furthermore, also other nonlinear effects such as third-harmonic generation (THG), four-wave mixing or the interactions with molecules can be implemented by using appropriate material models.
Finally, CEs can also be used for simulating more complex plasmonic systems in the nonlinear regime. In contrast to the linear case, where they turned out to be advantageous mainly for structures having radii of curvature r >5 nm [9], we expect them to be relevant also for smaller radii in this case, as nonlinear processes are more sensitive to the boundary approximation.
Acknowledgments
This work was supported by the Deutsche Forschungsgemeinschaft (DFG) GR 1288/10-1 as well as the Specific Target Research Project PLAISIR in the European Union (EU) Framework Program 7.
References and links
1. V. Sandoghdar, E. Klotzsch, V. Jacobsen, A. Renn, U. Håkanson, M. Agio, I. Gerhardt, J. Seelig, and G. Wrigge, “Optical detection of very small nonfluorescent nanoparticles,” Chimia 60, 761–764 (2006). [CrossRef]
2. M. T. Wenzel, T. Härtling, P. Olk, S. C. Kehr, S. Grafström, S. Winnerl, M. Helm, and L. M. Eng, “Gold nanoparticle tips for optical field confinement in infrared scattering near-field optical microscopy,” Opt. Express 16, 12302–12312 (2008). [CrossRef] [PubMed]
3. V. Deckert, “Tip-enhanced Raman spectroscopy,” J. Raman Spectrosc. 40, 1336–1337 (2009). [CrossRef]
4. P. Olk, J. Renger, T. Härtling, M. T. Wenzel, and L. M. Eng, “Two particle enhanced nano Raman microscopy and spectroscopy,” Nano Lett. 7, 1736–1740 (2007). [CrossRef] [PubMed]
5. J. Renger, R. Quidant, N. V. Hulst, and L. Novotny, “Surface-enhanced nonlinear four-wave mixing,” Phys. Rev. Lett. 104, 046803 (2010). [CrossRef] [PubMed]
6. T. Hanke, G. Krauss, D. Träutlein, B. Wild, R. Bratschitsch, and A. Leitenstorfer, “Efficient nonlinear light emission of single gold optical antennas driven by few-cycle near-infrared pulses,” Phys. Rev. Lett. 103, 257404 (2009). [CrossRef]
7. M. Ringler, A. Schwemer, M. Wunderlich, A. Nichtl, K. Kürzinger, T. A. Klar, and J. Feldmann, “Shaping emission spectra of fluorescent molecules with single plasmonic nanoresonators,” Phys. Rev. Lett. 100, 203002 (2008). [CrossRef] [PubMed]
8. G. Mie, “Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen,” Ann. Phys. 25, 377–445 (1908). [CrossRef]
9. A. Hille, R. Kullock, S. Grafstrom, and L. M. Eng, “Improving nano-optical simulations through curved elements implemented within the discontinuous Galerkin method,” J. Comput. Theor. Nanosci. 7, 1581–1586 (2010). [CrossRef]
10. J. Niegemann, M. Konig, C. Prohm, R. Diehl, and K. Busch, “Using curved elements in the discontinuous Galerkin time-domain approach,” AIP Conf. Proc. 1291, 76–78 (2010). [CrossRef]
11. J. Hesthaven and T. Warburton, “Nodal high-order methods on unstructured grids I. time-domain solution of Maxwell’s equations,” J. Comput. Phys. 181, 186–221 (2002). [CrossRef]
12. T. Lu, P. Zhang, and W. Cai, “Discontinuous Galerkin methods for dispersive and lossy Maxwell’s equations and PML boundary conditions,” J. Comput. Phys. 200, 549–580 (2004). [CrossRef]
13. K. Stannigel, M. König, J. Niegemann, and K. Busch, “Discontinuous Galerkin time-domain computations of metallic nanostructures,” Opt. Express 17, 14934–14947 (2009). [CrossRef] [PubMed]
14. J. Niegemann, M. König, K. Stannigel, and K. Busch, “Higher-order time-domain methods for the analysis of nano-photonic systems,” Photonics Nanostruct. Fundam. Appl. 7, 2–11 (2009). [CrossRef]
15. J. Niegemann, W. Pernice, and K. Busch, “Simulation of optical resonators using DGTD and FDTD,” J. Opt. A, Pure Appl. Opt. 11, 114015 (2009). [CrossRef]
16. J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications (Springer, 2008). [CrossRef]
17. R. Diehl, K. Busch, and J. Niegemann, “Comparison of low-storage Runge-Kutta schemes for discontinuous-Galerkin time-domain simulations of Maxwell’s equations,” J. Comput. Theor. Nanosci. 7, 1572–1580 (2010). [CrossRef]
18. J. Schöberl, “NETGEN an advancing front 2D/3D-mesh generator based on abstract rules,” Comput. Visual. Sci. 1, 41–52 (1997). [CrossRef]
19. C. Hafner, Post-Modern Electromagnetics: Using Intelligent MaXwell Solvers (John Wiley & Sons, 1999).
20. J. I. Dadap, J. Shan, K. B. Eisenthal, and T. F. Heinz, “Second-harmonic Rayleigh scattering from a sphere of centrosymmetric material,” Phys. Rev. Lett. 83, 4045–4048 (1999). [CrossRef]
21. G. Bachelier, J. Butet, I. Russier-Antoine, C. Jonin, E. Benichou, and P. Brevet, “Origin of optical second-harmonic generation in spherical gold nanoparticles: Local surface and nonlocal bulk contributions,” Phys. Rev. B 82, 235403 (2010). [CrossRef]
22. J. E. Sipe, V. C. Y. So, M. Fukui, and G. I. Stegeman, “Analysis of second-harmonic generation at metal surfaces,” Phys. Rev. B 21, 4389–4402 (1980). [CrossRef]
23. M. Scalora, M. A. Vincenti, D. de Ceglia, V. Roppo, M. Centini, N. Akozbek, and M. J. Bloemer, “Second- and third-harmonic generation in metal-based structures,” Phys. Rev. A 82, 043828 (2010). [CrossRef]
24. M. C. Larciprete, A. Belardini, M. G. Cappeddu, D. de Ceglia, M. Centini, E. Fazio, C. Sibilia, M. J. Bloemer, and M. Scalora, “Second-harmonic generation from metallodielectric multilayer photonic-band-gap structures,” Phys. Rev. A 77, 013809 (2008). [CrossRef]
25. J. Butet, G. Bachelier, I. Russier-Antoine, C. Jonin, E. Benichou, and P. Brevet, “Interference between selected dipoles and octupoles in the optical second-harmonic generation from spherical gold nanoparticles,” Phys. Rev. Lett. 105, 077401 (2010). [CrossRef] [PubMed]
26. M. Aeschlimann, M. Bauer, D. Bayer, T. Brixner, F. J. G. de Abajo, W. Pfeiffer, M. Rohmer, C. Spindler, and F. Steeb, “Adaptive subwavelength control of nano-optical fields,” Nature 446, 301–304 (2007). [CrossRef] [PubMed]
