Modal Analysis of the Ultrafast Dynamics of Optical Nanoresonators

Mar 24, 2017 - The FDTD calculations are performed with the FDTD Solutions software from Lumerical. Details about the PEEM modeling calculations can b...
6 downloads 15 Views 2MB Size
Download PDF
Article pubs.acs.org/journal/apchd5

Modal Analysis of the Ultrafast Dynamics of Optical Nanoresonators Rémi Faggiani,† Arthur Losquin,‡ Jianji Yang,† Erik Mårsell,‡ Anders Mikkelsen,‡ and Philippe Lalanne*,† †

Laboratoire Photonique, Numérique et Nanosciences (LP2N), UMR 5298, CNRS-IOGS-Univ Bordeaux, Institut d’Optique d’Aquitaine, 33400 Talence, France ‡ Department of Physics, Lund University, P.O. Box 118, SE-221 00 Lund, Sweden S Supporting Information *

ABSTRACT: We propose a semianalytical formalism based on a time-domain resonant-mode-expansion theory to analyze the ultrafast temporal dynamics of optical nanoresonators. We compare the theoretical predictions with numerical data obtained with the FDTD method, which is commonly used to analyze experiments in the field. The comparison reveals that the present formalism (i) provides deeper physical insight onto the temporal response and (ii) is much more computationally efficient. Since its numerical implementation is easy, the formalism, albeit approximate, can be advantageously used to both analyze and design ultrafast nano-optics experiments. KEYWORDS: optical resonators, ultrafast plasmonics, quasi-normal modes, nanoplasmonics, electron microscopy, spatiotemporal control

O

ptical nanoresonators illuminated by femtosecond laser pulses confine optical energy in subwavelength volumes and short periods of time. In particular, nanoresonators made up of metallic (or semiconductor) nanoparticles undergo ultrafast dynamics at time scales down to subfemtosecond due to their broad spectral bandwidth. As a result, the excitation of plasmonic nanoresonators by femtosecond lasers allows to control light−matter interaction in space and time with nanometric and subfemtosecond precisions.1,2 This capability has repercussions in fields as diverse as material science, cell biology and quantum optics.3 In addition to fundamental interests, it has enabled strong nonlinear optical phenomena to build intense high energy light sources,4 ultrashort electron sources5 or high contrast nanoscale probes.6 To make the most of this enormous potential, one needs to monitor the temporal response of the plasmonic nanoresonator under illumination by spatially varying electric field with complex temporal evolution, see Figure 1. The prospect of accessing experimentally this temporal response has stimulated the development of various, increasingly sophisticated spatially resolved techniques based on optical microscopy,6−9 electron microscopy,10−12 or scanning near-field optical microscopy.13−15 Such development is continuously fed by new theoretical proposals and discussions.16,17 Ultimately, a measurement would provide a quantitative estimation of the time evolution of the induced plasmonic field at nanometer and attosecond scales. Achieving optimum spatiotemporal control of optical energy at these scales using pulse shaping and learning algorithms2 would then be possible. © 2017 American Chemical Society

Figure 1. Light scattering by a plasmonic nanoresonator driven by an incident optical pulse Einc(r,t) defined in the absence of nanoresonator (scattered-field formulation). The scattered field is denoted by Esca(r,t) and the total field is simply E(r,t) = Einc(r,t) + Esca(r,t).

However, despite this strong experimental effort, theoretical modeling remains limited. The measurements reported so far are compared to simple damped harmonic oscillator analytical models,7−10,14,15,17 which provide intuitive pictures, but are only valid for simple samples and experimental conditions. Alternatively, the measurements are compared to classical numerical simulations (performed either in the spectral domain18 or in the time domain,19 mostly using FiniteDifference Time-Domain (FDTD) solvers), which can model arbitrary experiments but do not always provide a clear insight into the physics driving the temporal response. Indeed, in classical Maxwell’s equation solvers, the interaction between light and nanoresonators is described via continuum (scattering) theory, which does not access directly the core physical concepts attached to the response: the excitation and decay of a Received: December 13, 2016 Published: March 24, 2017 897

DOI: 10.1021/acsphotonics.6b00992 ACS Photonics 2017, 4, 897−904

ACS Photonics

Article

ωm is the complex frequency of the mth QNM and Δε(r,ω) the permittivity difference between the resonator and background medium.25 In principle, to be as accurate as possible, the expansion requires to consider many QNMs.21,22 However, considering a few dominant QNMs often achieves very good accuracy,24 the net advantage being that interpretation, understanding and design become simple and very fast. Referring to Figure 1, we now consider that the driving field is an optical pulse, Einc(r,t), that is, a wave packet that can be Fourier transformed Einc(r,t) = ∫ +∞ −∞Einc(r,ω) exp(iωt)dω, with Einc(r,ω) the frequency spectrum of the pulse. Driven by the incident pulse, the resonator scatters a time-dependent electromagnetic field, Esca(r,t). Every infinitesimal frequency component of the driving field Einc(r,ω)dω gives rise to an infinitesimal harmonic scattered field dEsca(r,ω). The latter is conveniently expanded as a sum of QNM modes dEsca(r,ω) ≈ ∑mdβm(ω)Em(r), with

few resonance modes. Furthermore, these calculations are extremely time-consuming and cannot be effectively used to design spatiotemporal control experiments.2 To the best of our knowledge, semianalytical tools simultaneously general, intuitive, accurate, and efficient do not presently exist. Hereafter, by relying in recent advances on the modeling of electromagnetic resonances with quasi-normal-mode (QNM) expansion techniques,20−25 we develop a semianalytical modal formalism to model the dynamics of nanoresonators. QNMs are eigensolutions of Maxwell’s equations, that is, solutions without external source. As such, they solely reflect the intrinsic properties of the nanoresonator, and have a very intuitive physical meaning: QNM’s are the natural resonant modes of the system (e.g., localized surface plasmon modes for metallic particles or Mie leaky modes for dielectric objects). Thus, by expending the scattered field scattered in the QNM basis, the scattering process is naturally mapped into an intuitive formalism, and this greatly facilitates design and interpretation. Recently, key theoretical works have nicely established the completeness of QNM expansions for simple open systems, for example, 1D and 3D spherical resonators.21 Interestingly, numerical results22 have also shown that the expansion provides remarkably accurate predictions for the scattering cross-section of 2D dielectric (nondispersive) resonators and TE polarization, provided that the “continuum” of radiation QNMs is included in the expansion.23 The present work is not concerned by such theoretical issues, but rather by the derivation of approximate expansions. If the scattering process can be reduced to the excitation of a few dominant resonant modes, as is often the case, simple expansions may be used to interpret experimental results. Intuition and understanding are then made available, in sharp contrast with more widespread classical scattering theories that do not intrinsically rely on the natural modes of the resonator. Thus, the general and insightful character of the modal expansion is expected to advantageously help engineering new spatiotemporal control schemes,1,2,7,26 design original geometries and develop new applications in the field of ultrafast nano-optics. In the following, we will emphasize the gain brought by the modal formalism to interpret experimental measurements. We will then quantitatively test its accuracy by comparing the results obtained on a complex plasmonic structure with those obtained using widespread numerical methods. The formalism shows similarities with earlier works on quasistatic theories based on modal decompositions.1,27 However, in sharp contrast, its scope of validity is not restricted to tiny deep-subwavelength resonator sizes as it relies on fullwave computations of QNMs.

dβm(ω) =

− ω dω ω − ωm

∫ ∫ ∫ Δε(r , ω)Einc(r , ω)·Em(r )dr 3 (1)

Assuming that the scattered field can be inverse Fourier transformed, we obtain the scattered field in the time domain by summing up all the frequency components, Esca(r,t) = +∞ ∫ −∞ dEsca(r,ω) exp(iωt). The latter can be conveniently expressed as Esca(r , t ) ≈

∑ βm(t )Em(r )

(2)

m

with +∞

βm(t ) = +∞

=

∫−∞

dβm(ω) exp(iωt )

ω exp(iωt ) ωm − ω 3 Em(r )dr dω

∫−∞

∫ ∫ ∫ Δε(r , ω)Einc(r , ω)· (3)

The total field induced by the optical pulse is thus E(r,t) = Esca(r,t) + Einc(r,t), with Esca(r,t) given by eqs 2 and 3. Because of the linearity in the time-domain, the scattered field remains a simple sum of independent contributions from every single QNM retained in the expansion. These contributions are weighted by the time-dependent excitation coefficients βm(t). The latter depend on the resonant nature of the interaction through the 1/(ωm − ω) terms and the coupling to the driving field through the inner product Einc(r,ω)·Em(r). Equations 2 and 3 provide a very simple and intuitive description of the temporal response of nanoresonators with decoupled time and space dependences, as we will see in the experimental modeling section. Let us note that the formalism is valid for any analytical driving field, for example, a field scattered by another resonator, a dipole emission in the near-field of the resonator or a far-field focused laser pulse, possibly polarization-shaped.2 An experimental spectrum could also be fed into the calculation. Hereafter we will consider a plane-wave driving field propagating toward the positive z-axis, defined in the spectral domain by Einc(r,ω) = Wg(ω − ω0) exp(−ik(ω)z), with ω0 the central frequency of the pulse, k(ω) the wavevector modulus in the background medium, g the frequency spectrum of the pulse,

■

THEORETICAL METHOD We start by considering recent results obtained for QNMexpansions in the spectral domain for harmonic fields. In the scattering formulation, the total field is decomposed as a sum of the driving field Einc(r,ω) at frequency ω and the scattered field Esca(r,ω) by the nanoresonator. Esca can be expanded in a finite set of normalized QNMs with spatial electric distributions Em(r) labeled by the integer m Esca(r,ω) ≈ ∑mβm(ω)Em(r). Provided that a few resonance modes are dominantly excited, the expansion admits accurate closed-form expressions even for strong radiation leakage, absorption, and material dispersion. The coupling coefficients βm(ω) are simply obtained as overlap integrals between the normalized modes and the driving field βm(ω) = −ω∫ ∫ ∫ Δε(r,ω)Einc(r,ω). Em(r)dr3/(ω − ωm), where 898

DOI: 10.1021/acsphotonics.6b00992 ACS Photonics 2017, 4, 897−904

ACS Photonics

Article

and W the 3 × 1 electric-field vector of the incident plane wave. Substituting into eq 3, we obtain ∞

βm(t ) =

∫−∞

ω exp(iωt )g (ω − ω0) ωm − ω

∫ ∫ ∫ Δε(r , ω) exp(−ik(ω)z)W ·Em(r )dr 3dω (4)

The βm(t)’s are obtained by first computing the volume integral that depends on the frequency and then the frequency integral that is conveniently performed using Fast Fourier Transform algorithms. The present formalism is very flexible. Indeed, once the few dominant QNMs are computed, any change in the driving field (nature, polarization, propagation direction, central frequency, pulse duration, etc.) is straightforwardly taken into account without rerunning the Maxwell’s equations solver because the formalism decouples the intrinsic properties of the resonator from the driving field. The latter could thus be efficiently adjusted with a feedback loop for designing coherent control experiments for example.2 This markedly contrasts with classical numerical methods such as the FDTD method that requires to redo the entire computation for every instance of the driving field. To compute and normalize the QNMs, we use a freeware25 that relies on a Finite Element Method (FEM) solver. A few minutes, typically 2−3 min, depending on the mesh, are required to calculate every QNM with high accuracy on a low performance workstation. Details about the calculations can be found in the Methods section.

Figure 2. Difficulty encountered when interpreting PEEM measurements of silver nanorice with FDTD calculations. (a) Left: Experimental photoemission intensity induced by two delayed (τ = 5.5 fs) laser pulses at the two ends of a silver nanorice, measured by interferometric time-resolved PEEM in ref 19. Each curve is normalized to its maximum. The black arrows evidence a duplication of the intensity peaks, and suggest a markedly different and complex time-evolutions of the fields. Right: Scanning electron microscope image of the nanorice. The colored arrows indicate the rice ends that correspond to the colored curves. The slight asymmetry with respect to τ = 0 is an experimental artifact. Details about the experiment can be found in ref 19 and in Figure 3a. (b) Normalized time-evolution of the fields at the two ends of the nanoparticle, obtained from FDTD computations. (c) Normalized nonlinear autocorrelation function 6 ∫ +∞ −∞|E(r,t) + E(r,t + τ)| dτ of the fields shown in (b). Consistent with the experimental data in (a), the autocorrelation functions markedly differ and exhibit complex oscillations for delays ≥5 fs (see the black arrows).

■

PEEM EXPERIMENTAL DATA To evidence the benefit brought by QNM expansions to interpret measurements in ultrafast spectroscopy, we start by reminding some difficulties encountered when interpreting experiments with available modeling tools. For the sake of illustration, we use time-resolved PhotoEmission Electron Microscopy (PEEM) experiments as a prototypical example for probing the temporal dynamics of plasmonic fields. Specifically, we consider interferometric time-resolved PEEM measurements of silver rice shaped nanoparticles (Figure 2a) recently performed by some of the authors of the present work. In the experiment reported in ref 19, the nanoparticles are deposited on a glass substrate covered with ITO and are illuminated by 5.5 fs laser pulses. The driving pulses are splitted into two identical delayed replicas, and PEEM images of the particles are recorded by varying the delay. An important observation of the experiments was the marked and systematic difference of the PEEM intensity oscillations at the two ends of the particles (indicated by the blue and red arrows on the Scanning Electron Microscope image of Figure 2a) for delays ≥5 fs. For some particles, we even measured complicated oscillations (see black arrows in Figure 2a), suggesting a complex time evolution of the plasmonic fields despite the high degrees of simplicity and symmetry of the geometry.19 Because the intensity oscillations as a function of delay τ in PEEM experiments only indirectly reflect the time evolution of the plasmonic fields, it is difficult to interpret the measurements with simple analytical models. Thus, consistently with the literature, we compared the experimental oscillations with those obtained from a third-order autocorrelation-function model:

+∞ ∫ −∞ |E(r,t) + E(r,t + τ)|6dτ, where E(r,t) is a classical plasmonic field.12,17 In the initial interpretation in ref 19, E(r,t) was exclusively computed with the FDTD method. As shown in Figure 2c, the third-order autocorrelation oscillations qualitatively reproduce the PEEM measurements, confirming the existence of a complex temporal evolution of the plasmonic fields. In ref 19, we tentatively attributed this complex temporal evolution to the interplay between (i) the broadband character of the driving field, (ii) its oblique incidence, and (iii) the particle, themselves, which operate well beyond the quasistatic regime and support several eigenmodes. This interplay lead us to interpret the marked difference of the PEEM intensity oscillations as an intrictate delocalized induction of surface wave packets of different natures, which interfere with each other as well as with the driving pulse. However, we could not capture the details of the interference from the field calculated by FDTD, and our interpretation of the experiment remained limited. In contrast, by disentangling the intrinsic properties of the object from the excitation, the present QNM expansion method will allow an in-depth interpretation of the temporal evolution of the field probed by the PEEM experiment, as we will see in the next section.

■

INTERPRETATION OF THE PEEM EXPERIMENTAL DATA To compute the QNMs, we model the nanorice as an ellipsoid of length 600 nm and radius 50 nm positioned above a substrate of permittivity 2.1 (see Figure 3a). The permittivity of 899

DOI: 10.1021/acsphotonics.6b00992 ACS Photonics 2017, 4, 897−904

ACS Photonics

Article

and scattering) at a decay rate given by Im(ωm). Referring to eq 2, since the coupling coefficients βm do not depend on position r, the spatial dependence of the response of every individual mode is solely determined by the spatial distribution of the mode. Accordingly, the field scattered by a mode shows the same normalized envelop and oscillations at all spatial locations, and only amplitudes or phases may differ. In the present case, the fields scattered by mode 1 at both ends are identical by symmetry, and the fields scattered by mode 2 are π-phase shifted (see Figure 3b). The total field (Figure 3f) is given by the sum of the three contributions. It is highly different at the two ends. Notably, the signal envelops strongly differ, both in amplitude and shape. In addition, the field oscillations are distorded for times ranging from 10 to 25 fs. Overall, the QNM expansion corroborates the observations in ref 19 and reproduces the FDTD calculations of Figure 2b (the slight differences between Figures 2b and 3f are mostly due to slightly different calculation parameters). However, compared to the FDTD calculations, the decomposition into different contributions greatly eases the interpretation and brings new insights into the underlying phenomena. For instance, the difference in the response at both ends is understood as resulting from the spatial dependence of the interferences between the different contributions, namely the difference in the arrival time of the driving pulse and the antisymmetric character of mode 2, whose contributions are out-of-phase at point A and B (Figure 1c,e). It also clearly appears that the distortions in the field oscillations result from the interferences between the two QNMs that oscillate at their respective frequencies Re(ωm) for times greater than 15 fs. The interferences induce a beating at the frequency Re(ω2 − ω1). However, since the beating frequency is only one-third of the central frequency of the driving pulse, the usual beating pattern with marked nodes and antinodes is not observed and is replaced by a distortion appearance. By identifying the individual driving contributions, the present method thus provides unique insights on the temporal response. Future directions for quantitative matching between experiment and theory may include a better pulse characterization28 to implement accurate excitation pulses in the calculation, the use of refined semiclassical models for photoemission29 or photoluminescence,30 and the inclusion of ultrashort changes in the metal dielectric permittivity due to nonthermal electron− hole pairs and electron−electron scattering for instance.31−33 As we show in the Supporting Information, where we calculate autocorrelation functions, the QNM-expansion method also reproduces the PEEM experimental results. Our calculations show that the complex PEEM oscillations as a function of delay (Figure 2a) are due to the simultaneous excitation of two plasmon modes by the ultrashort laser pulses. Indeed, the complex oscillations as a function of delay reflect complex field dynamics at the particle boundaries that happen only when two modes interfere. We note that quasi-static theories cannot be used to quantitatively interpret the experiments, since they predict that mode 1, which has a zero dipole moment, cannot be excited27 (see Supporting Information for further details). The QNM-expansion method is already very relevant for nanoresonators as simple as nanorices and also accounts for the presence of the substrate, in contrast with analytical field expressions of simple nanoparticles. Its relevancy should further grow with complexity. Additionally, the method is not restricted to PEEM experiments and can be used to model a

Figure 3. Modal analysis of the nanorice temporal response probed by the PEEM experiment. (a) The particle is modeled as a silver ellipsoid deposited on a glass substrate and illuminated by a p-polarized 5.5 fs plane-wave Fourier limited Gaussian pulse of central frequency ω0 = 2.43 × 1015 rad/s−1 propagating with a 25° incidence degree with respect to the substrate. (b) Ez-spatial variations of the normalized dominant QNM modes. (c−f) Temporal response of the Ez fieldcomponent computed at points A (red) and B (blue). (c) Driving pulse, which includes the reflection of the incident plane wave on the substrate in the absence of the particle. (d) Field scattered by mode 1. (e) Field scattered by mode 2. (f) Total field.

silver is approximated by a Drude model with plasma frequency of 1.25 × 1016 rad/s and relaxation rate of 1.6 × 1014 s−1. In order to match the experimental conditions, the driving field is considered to be a 5.5 fs Fourier-limited plane-wave Gaussian pulse with a central frequency ω0 = 2.43 × 1015 rad/s. The pulse is assumed to be p-polarized and obliquely incident with an angle of 25° with respect to the substrate (Figure 3a). In the spectral range spanned by the pulse, the particle supports two dominant QNMs whose spatial distributions in the plane of incidence are shown in Figure 3b. Their respective complex resonance frequencies are ω1 = 2.08 × 1015(1 + i/(2Q1)) rad/s and ω2 = 2.95 × 1015(1 + i/(2Q2)) rad/s, Q1 = 5.4 and Q2 = 8.1 being the quality factors of the resonances. The fundamental dipolar mode is not considered in the expansion, as its resonance frequency does not match the spectral range covered by the excitation pulse. The QNMs are independent from the driving field and are, thus, either symmetric (mode 1) or antisymmetric (mode 2) with respect to the long x-axis of the ellipsoid. Additionally, due to the presence of the substrate, a small vertical asymmetry is present. In Figure 3c−e, we show the three contributions (excitation field, field scattered by mode 1 and field scattered by mode 2) to the total field at the two ends of the particle, denoted by points A and B in Figure 3a. We solely consider the zcomponent of the field as it is the dominant component driving the nonlinear photoemission process recorded in the experiment. Because of the oblique incidence, the driving field exhibits a phase delay corresponding to the propagation time between points A and B, see Figure 3c. The fields scattered by modes 1 and 2, see Figure 3d,e, resemble the usual temporal response of any driven dissipative resonators (or damped harmonic oscillators) with markedly different oscillation frequencies given by Re(ωm) and relaxations (by absorption 900

DOI: 10.1021/acsphotonics.6b00992 ACS Photonics 2017, 4, 897−904

ACS Photonics

Article

23.1, and Q3 = 8.7. Modes 1 and 2 preferably couple to xpolarized driving fields and mode 3 to y-polarized fields. Figure 5a shows the FDTD and QNM-expansion predictions for the x-component of the total electric field, Ex(rA,t)

variety of experiments measuring various nonlinear signals (harmonic generation, photoluminescence or photoemission).

■

QUANTITATIVE TEST OF THE MODAL FORMALISM The previous analysis was meant to demonstrate the potential of the QNM-expansion method to provide intuitive interpretations of experimental measurements. Hereafter we aim at showing the versatility of the method and at testing its capabilities. We thus study a more complex geometry and directly compare the QNM-expansion predictions with numerical data obtained with the FDTD method. We choose a Dolmen nanostructure composed of three gold nanorods, see Figure 4a. The Dolmen provides a richer physics

Figure 5. Temporal and spectral analysis of the Dolmen. (a) Comparison of the temporal responses computed at point A = (−5,−10,0) nm (the origin is shown by the axes) with the FDTD method (blue) and the QNM-expansion formalism (red), for a 12.7 fs Gaussian incident pulse propagating along the ẑ direction and polarized along the x̂ + 2ŷ direction (see the rightmost inset). Bottom-left/center/right insets: Enlarged views of the temporal response. The black curves represent the temporal response obtained by Fourier transforming the FEM “exact” data computed shown in (b). (b) Spectrum of the field scattered by the Dolmen at point A under illumination by a plane wave with unit amplitude and same polarization as in (a). The dots “are” exact FEM data. The colored curves show the contributions of modes 1, 1 + 2, and 1 + 2 + 3; see the main text for more details.

Figure 4. Dolmen QNMs. (a) The Dolmen is composed of three gold rods. The geometrical parameters are hs = 100 nm, ws = g = G = 30 nm, hb = 50 nm, and wb = 128 nm. The rod thickness is 20 nm. (b−d) Field amplitudes of the normalized QNM modes. The small yellow arrows represent the local direction of the electric field. The white arrows represent the equivalent electric dipoles of every rod. For the simulations, the gold permittivity is approximated by a Drude model with a plasma frequency of 1.365 × 1016 rad s−1 and a relaxation rate of 3.2 × 1013 s−1 and the Dolmen is assumed to be in an air background refractive index of 1.

computed at point A = (−5,−10,0) nm. Both methods evidence a clear signature of mode beating, with a beating period much greater than the nanorice one due to the smaller frequency difference between the modes. However, a more attentive inspection reveals two main discrepancies. First, the QNM expansion predicts oscillation amplitudes that are larger than those obtained with the FDTD method, especially for the first beating. Second, the responses are in quadrature for large t’s and slightly shifted for small ones, as respectively shown in the right and left lower insets. To investigate the origin for these deviations, we further perform frequency-domain computations using the same FEM Maxwell’s equations solver as the one used for the QNM computation (see the Methods section). The computed data are shown with dots in Figure 5b. They will be referred to as “exact” hereafter and are obtained for the field Esca,x(rA,ω) scattered by the Dolmen at point A under illumination by an

than the nanorice, with a Fano response resulting from the interference of out-of-phase modes34,35 that can be exploited for sensing and implementing in plasmonic metamaterials.36 In the following computations, the Dolmen is assumed to be illuminated by a 12.7 fs plane-wave Gaussian pulse (central frequency ω0 = 2.9 × 1015 rad/s) propagating in the z-direction and linearly polarized along the x̂ + 2ŷ direction, x̂ and ŷ being the unitary vectors along the x- and y-axis. At visible frequencies, the Dolmen supports three dominant QNMs, shown in Figure 4b−d. For the sake of understanding, the inplane equivalent electrical dipoles of the nanorods, determined from the field lines, are depicted with white arrows in the figure. The QNM complex resonance frequencies are ω1 = 2.70 × 1015(1 + i/(2Q1)) rad/s, ω2 = 2.97 × 1015(1 + i/(2Q2)) rad/s, and ω3 = 3.11 × 1015(1 + i/(2Q3)) rad/s with Q1 = 19.1, Q2 = 901

DOI: 10.1021/acsphotonics.6b00992 ACS Photonics 2017, 4, 897−904

ACS Photonics harmonic plane wave E inc(r , ω) =

Article

(

x̂ 5

,

2y ̂ 5

) (

ω

, 0 exp −i c z

appears much more stabilized than the FDTD ones. Ultimately, the high degree of accuracy together with the high level of simplicity constitute a direct evidence of the force of the present formalism.

)

with the same incident angle and polarization vector as the ones used for the time-domain computations. We first compare these data with those obtained with the QNM expansion, Esca,x(rA,ω) ≈ ∑m=1−3βm(ω)Em(r), with the βm’s given by eq 1 and shown with the red-solid curve. We use the same parameters for both calculations, so that the discrepancies between the two sets of data directly assess inaccuracy issues inherent to the small number of QNMs retained in the QNM-expansion method. Although it only qualitatively predicts the magnitude of the spectral response, the QNM expansion well predicts the complex shape of the response. It is remarkable that the shape can be unambiguously attributed to the interference between modes 1 and 2 (blue solid curve) and that the frequencies of the extrema of the complex spectral response are all accurately predicted with an expansion involving only 3 QNMs (red curve). This evidence that the physics of the electromagnetic interaction between the three nanorods is well captured by the QNM expansion and that the discrepancy between the “exact” (dots) and QNM expansion (red-solid) data can be attributed to higher-order QNMs that are neglected in the expansion and are likely to exhibit rather smooth spectral responses in the spectral range of interest, away from their own resonance frequencies. It is thus tempting to attribute the deviation between the FDTD and QNM methods in Figure 5a to the QNM-expansion inaccuracies revealed in Figure 5b. But it is not that simple. To evidence it, we finally compute an “exact” FEM field Ex(rA,t) by weighing the “exact” data Esca,x(rA,ω) by the pulse spectrum and then Fourier transform the weighted product. A comparison between the time-domain “exact” field obtained from FEM computations, shown with black dots in the lower insets of Figure 5a, and the FDTD- and QNM-method predictions allows to define three regimes. For small t’s < 30 fs (leftmost inset), the “exact” FEM and FDTD results quantitatively agree, evidencing that the initial temporal response is only qualitatively predicted with the QNM formalism. It is reasonable to assume that the pulse excites higher-order QNMs with large complex frequencies, which are neglected in the expansion. For intermediate t’s (central inset), a transition regime is observed, in which the “exact” FEM data start to depart from the FDTD results and become nearly superimposed with the QNM-expansion results. The transition regime occurs for t ≈ 40 fs, which approximately corresponds to the beginning of the relaxation of the three dominant QNMs. It is thus a characteristic time beyond which the impact of higher-order QNMs becomes less significant. Finally, for longer t’s (rightmost inset), the exact FEM and QNM-expansion results perfectly match, whereas the FDTD results are phase-shifted by ≈π/2. Only the fundamental QNMs contribute to the response. Although the FDTD method is clearly a reference method that is widespread for time-domain analysis, its performance deteriorates when implementing metal dispersion with auxiliary differential equations.37 The second-order accuracy of the Yee grid is no longer guaranteed at metal−dielectric interfaces even for straight interfaces and square corners, and the FDTD predictions for metallic structures suffer from numerical errors that are known to far exceed those of FEMs.38 As shown in the Supporting Information, where we study the accuracy of the QNM expansion and FDTD methods for different mesh discretizations, the convergence of the QNM-expansion data

■

CONCLUSION We have presented a semianalytical time-domain QNMexpansion formalism to predict the temporal response of optical nanoresonators. The method relies on an expansion of the scattered field onto a small set of resonance modes, which may be computed for most nanoresonators with effective QNM solvers.25 The method is intrinsically simple and approximate since only a few QNMs are retained in the expansion. Comparison with the FDTD method has evidenced two major assets. First, by providing a direct access to the few control knobs, the QNM excitation coefficients that are driving the temporal response, the QNM expansion eases the interpretation of experiments and provides deeper insights. Second, the QNM expansion provides highly accurate predictions when the temporal response is driven by a few dominant QNMs, as is often the case, and third, the high accuracy is achieved for very small CPU times of the order of a few minutes with a standard workstation. This implies that the QNM expansion may be helpful for inverse design39 to engineer, for instance, nanoscale pulse shapers. In the future, it would be interesting to investigate the limitations imposed by the use of a low number of QNMs. Other interesting perspectives concern the extension of the formalism toward molecular plasmonics,40 that is, hybrid systems composed of molecule ensembles that collectively interact with plasmonic resonances.

■

METHODS The normalized QNMs are computed using a freeware41 implementing a method described in ref 25. They are found using an iterative procedure to fit to a Padé approximated polelike response function. Let us note that the pole search in the complex plane requires analytical functions to describe the materials permittivity. The iterative procedure is carried out utilizing the COMSOL Multiphysics Finite Element Method (FEM) solver driven by a MATLAB code available in ref 41. The QNM spatial distributions are obtained from the field scattered by the resonator at a frequency close to the resonance frequency. The coupling coefficients βm(t) are calculated using eq 4. The spatial integral is performed directly within COMSOL Multiphysics. To simplify, the frequency dependence of Δε(r,ω) and exp(−ik(ω)z) in the spectral window defined by the resonance term 1/(ω − ωm) is neglected, so that the spatial integral is performed only once at a single frequency equal to Re(ωm). We have numerically checked that accounting for the dispersion by computing the spatial integral for every frequency provides nearly identical results, see details in the Supporting Information. Finally, the integral over ω in βm(t) is performed using the Fast Fourier Transform algorithm. A MATLAB code is provided in ref 41 to compute the βm(t) for the Dolmen geometry. The Dolmen QNM are calculated in a spherical simulation space of radius 600 nm, a PML layer of thickness 300 nm and second-order finite elements. An extremely refined mesh (∼ λ/ 300) is used close to the Dolmen surface to finely sample the rapid variations of the field. In total, the simulation space 902

DOI: 10.1021/acsphotonics.6b00992 ACS Photonics 2017, 4, 897−904

ACS Photonics

Article

(4) Kim, S.; Jin, J.; Kim, Y. J.; Park, I. Y.; Kim, Y.; Kim, S. W. High harmonic generation by resonant plasmon field enhancement. Nature 2008, 453, 757. (5) Vogelsang, J.; Robin, J.; Nagy, B. J.; Dombi, P.; Rosenkranz, D.; Schiek, M.; Groß, P.; Lienau, C. Ultrafast Electron Emission from a Sharp Metal Nanotaper Driven by Adiabatic Nanofocusing of Surface Plasmons. Nano Lett. 2015, 15, 4685. (6) Berweger, S.; Atkin, J. M.; Olmon, R. L.; Raschke, M. B. Adiabatic Tip-Plasmon Focusing for Nano-Raman Spectroscopy. J. Phys. Chem. Lett. 2010, 1, 3427. (7) Hanke, T.; Krauss, G.; Trautlein, D.; Wild, B.; Bratschitsch, R.; Leitenstorfer, A. Efficient Nonlinear Light Emission of Single Gold Optical Antennas Driven by Few-Cycle Near-Infrared Pulses. Phys. Rev. Lett. 2009, 103, 257404. (8) Brinks, D.; Castro-Lopez, M.; Hildner, R.; van Hulst, N. F. Plasmonic antennas as design elements for coherent ultrafast nanophotonics. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 18386. (9) Accanto, N.; Piatkowski, L.; Renger, J.; van Hulst, N. F. Capturing the optical phase response of nanoantennas by coherent second-harmonic microscopy. Nano Lett. 2014, 14, 4078. (10) Kubo, A.; Onda, K.; Petek, H.; Sun, Z.; Jung, Y. S.; Koo Kim, H. Femtosecond Imaging of Surface Plasmon Dynamics in a Nanostructured Silver Film. Nano Lett. 2005, 5, 1123. (11) Barwick, B.; Flannigan, D. J.; Zewail, A. H. Photon-induced near-field electron microscopy. Nature 2009, 462, 902. (12) Aeschlimann, M.; Bauer, M.; Bayer, D.; Brixner, T.; Cunovic, S.; Dimler, F.; Fischer, A.; Pfeiffer, W.; Rohmer, M.; Schneider, C.; Steeb, F.; Strü ber, C.; Voroninec, D. V. Spatiotemporal control of nanooptical excitations. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 5329. (13) Onishi, S.; Matsuishi, K.; Oi, J.; Harada, T.; Kusaba, M.; Hirosawa, K.; Kannari, F. Spatiotemporal control of femtosecond plasmon using plasmon response functions measured by near-field scanning optical microscopy (NSOM). Opt. Express 2013, 21, 26631. (14) Nishiyama, Y.; Imura, K.; Okamoto, H. Observation of Plasmon Wave Packet Motions via Femtosecond Time-Resolved Near-Field Imaging Techniques. Nano Lett. 2015, 15, 7657. (15) Kravtsov, V.; Ulbricht, R.; Atkin, J. M.; Raschke, M. B. Plasmonic nanofocused four-wave mixing for femtosecond near-field imaging. Nat. Nanotechnol. 2016, 11, 459. (16) Stockman, M. I.; Kling, M. F.; Kleineberg, U.; Krausz, F. Attosecond nanoplasmonic-field microscope. Nat. Photonics 2007, 1, 539. (17) Aeschlimann, M.; Brixner, T.; Fischer, A.; Hensen, M.; Huber, B.; Kilbane, D.; Kramer, C.; Pfeiffer, W.; Piecuch, M.; Thielen, P. Determination of local optical response functions of nanostructures with increasing complexity by using single and coupled Lorentzian oscillator models. Appl. Phys. B: Lasers Opt. 2016, 122, 199. (18) Sun, Q.; Ueno, K.; Yu, H.; Kubo, A.; Matsuo, Y.; Misawa, H. Direct imaging of the near field and dynamics of surface plasmon resonance on gold nanostructures using photoemission electron microscopy. Light: Sci. Appl. 2013, 2, e118. (19) Marsell, E.; Losquin, A.; Svärd, R.; Miranda, M.; Guo, C.; Harth, A.; Lorek, E.; Mauritsson, J.; Arnold, C. L.; Xu, H.; L’Huillier, A.; Mikkelsen, A. Nanoscale Imaging of Local Few-Femtosecond NearField Dynamics within a Single Plasmonic Nanoantenna. Nano Lett. 2015, 15, 6601. (20) More, R. M.; Gerjuoy, E. Property of resonance wave functions. Phys. Rev. A: At., Mol., Opt. Phys. 1973, 7, 1288. (21) Leung, P. T.; Liu, S. Y.; Young, K. Completeness and orthogonality of quasinormal modes in leaky optical cavities. Phys. Rev. A: At., Mol., Opt. Phys. 1994, 49, 3057. (22) Vial, B.; Zolla, F.; Nicolet, A.; Commandré, M. Quasimodal expansion of electromagnetic fields in open two-dimensional structures. Phys. Rev. A: At., Mol., Opt. Phys. 2014, 89, 023829. (23) Doost, M. B.; Langbein, W.; Muljarov, E. A. Resonant state expansion applied to two-dimensional open optical systems. Phys. Rev. A: At., Mol., Opt. Phys. 2013, 87, 043827.

consists of 329000 tetrahedral elements. To compute the QNM, we use two mirror-symmetry planes. Only 11 min are required to compute and normalize every individual QNM with the very fine mesh, using the iterative pole search, with an Intel Xeon(R) E5−2637 CPU, 3.50 GHz, 64 Go RAM PC. The CPU time to compute the βm(t) coefficients is negligible (