Modes interplay and dynamics in the second harmonic generation of plasmonic nanostructures

The full wave surface integral equation computation of the second harmonic generation (SHG) dynamics for metal spheres and nanorods presented as multimedia files is performed to reveal the dynamics of the modes supported by the nanostructure. We demonstrate that the interplay between different modes controls the nonlinear response and that the size-induced redshift of the eigenmodes can be manipulated by adjusting the nanostructure geometry, so that the SHG signal can be boosted at specified frequencies. We show that the SHG radiation is not necessarily quadrupolar in spherical nanoparticles, as it is often assumed. Finally, we introduce an efficient way to reduce the SHG calculation time. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
After the invention of the first pulsed Ruby laser [1,2], it became possible to effectively compress monochromatic laser radiation into tightly confined beams with extremely intense electric fields. This led to the first experimental observation of the second harmonic generation (SHG) from a quartz crystal [3]. Numerous interesting theoretical and experimental ideas have then been developed on this topic [4,5]. Particularly interesting effects appear in nonlinear light conversion with plasmonic nanostructures because of their ability to localize light beyond the diffraction limit and create hot spots, which boost nonlinear conversion [6]. These hot spots are associated with localized surface plasmon resonances (LSPR) supported by the nanostructures [7][8][9][10]. Such resonances appear first at the fundamental frequency and can then be excited by the nonlinear fields at the doubled frequency [11][12][13][14][15][16][17][18][19][20].
The study of such resonances can be performed in two different regimes. Generally, in order to achieve a large enough instantaneous nonlinear response, the illumination is required to be pulsed. Depending on the duration of the pulse with respect to the lifetime of the modes supported by the system (typically several tens of femtoseconds), the response will be different. In the case of long pulses (say, with a duration above 100 fs) the corresponding spectrum is narrow (the spectral width for a 100 fs pulse centered at λ=800 nm is about 3 nm) and can be considered as a quasi-monochromatic laser excitation, which allows retrieving the response at a selected frequency (see theoretical [11,[21][22][23][24][25][26] and experimental [19,[27][28][29][30][31][32][33] publications on this topic). The far-field SH response in such experiments can be decomposed in terms of multipolar moments, which oscillate harmonically at the pulse's central frequency. On the other hand, the other regime involves pulses with broad spectra and correspondingly short durations, comparable or smaller than the lifetime of plasmonic modes. Such approach was used to investigate linear and nonlinear response of plasmonic systems [33][34][35][36][37][38][39]. In the case of a short duration pulse, in addition to the interference of the modes typical for both regimes [33], one can also have the situation in which one mode persists while the other one has already decayed [38]. Such experiments provide new degrees of freedom, for example to dynamically manipulate SH energy scattering directivity [33]. In this article, we study the second harmonic (SH) response of a single nanoparticle to a plane-wave Gaussian pulse, with emphasis on the second regime, such that we can investigate the modes interplay that is often observed experimentally.
Analysis of the response to a spectrally broad pump can take advantage from the theoretical and experimental data of a quasi-monochromatic response of the nanoparticle because the latter one provides important information about the supported modes. Before discussing the results of this paper, we must first detail the mechanisms of SHG. The linear scattering of the system under harmonic illumination (the distribution of the electromagnetic field at the surface of the scatterer) represents the response at the fundamental excitation frequency ω, which can then be used to compute the scattering at the SH 2ω; in the following, we will call this approach SH scattering. It is suitable for long pulses (above 100 fs), for which it is sufficient to consider only SH scattering at twice the central frequency if the pulse is spectrally narrow. However, as the pulse duration reduces and the spectrum width increases this quasi-monochromatic approximation becomes inaccurate and a more complex computation is required to capture the subtleties of SHG. The nonlinear process responsible for SHG is of the second order and can occur with any combination of two frequencies present in the pulse [4]. In that process, which is called sum frequency generation (SFG), two photons ω 1 and ω 2 are combined to produce a nonlinear photon at Ω = ω 1 + ω 2 . In the following, we will use SH scattering spectrum to support the discussion of the SFG far-field scattering spectrum and dynamics.
Theoretical [11,[21][22][23][24][25][26] and experimental [19,[27][28][29][30][31][32] works have been performed with narrow spectrum quasi-monochromatic laser excitations to study the mechanisms of SHG from spherical nanoparticles of various sizes and materials. Theoretically, dipolar and quadrupolar modes provide comparable contributions to the SH far-field distribution in such quasi-monochromatic experiments, even for the smallest aluminum nanoparticles of 5 nm diameter [11,22]. Experiments with 32 nm diameter silver spheres confirm distinguishable dipolar and quadrupolar SH radiation [28]. However, some experimental works find the dipolar response to be dominating from small silver spheres up to 30-40 nm [30], which was explained by the deviation from a perfect spherical shape [26]. However, the scaling dependence of the SH signal does not unambiguously prove that dipolar pattern appears only due to these deviations [26]. When the diameter of the sphere increases, the retardation effect enhances the quadrupolar contribution [11]. Although, in the theoretical full wave analysis of the scattering from a 20 nm diameter gold sphere, a dipolar response was found to be dominating [23]. Obviously, the interplay between different optical modes in the nonlinear response of plasmonic structures is quite intricate. This is the focus of this manuscript, where we show that SH quadrupolar and dipolar responses can effectively be obtained from a 40 nm diameter silver sphere and outline the necessary conditions to get them. The analysis of larger 120 nm length nanorod is also preformed, showing enhanced quadrupolar and additional small octupolar components, which are expected when the retardation effects become more noticeable as the object size increases. The same effect was observed for a 100 nm gold sphere [32].
In a previous study, we unveiled some mode interplay dynamics for the SFG response for a nanorod [38], the aim of this paper is to provide a comprehensive study of the influence of the size and shape of the structures and take advantage of the animations possible in this Journal. The results are organized in the following manner. In the Methods section, we describe the numerical frequency-and time-domain methods to calculate SFG and SH scattering from metallic nanoparticles of various shapes and sizes. We also demonstrate a new approach to reduce the number of computational steps for SFG dynamics calculations. In the Results section, we present the SFG scattering from silver nanostructures, namely 20 nm and 40 nm diameter spheres as well as a 40×40×120 nm 3 nanorod, and study them using the multipole decomposition. Enhancement of the quadrupolar radiation pattern contribution is observed for the abovementioned nanoparticles as the size of the system increases. For the sphere, the resonance redshift caused by the size increment is observed. We show the ability to control the resonance positions by considering a nanorod with the aspect ratio adjusted in such a way that the positions of spectral peaks appear close to the one observed for the 20 nm sphere [40]. We also propose the empirically derived mechanisms for the SH modes formation in nanorods, which is very similar to that of a sphere.

Methods
In this paper, we study the SFG time-domain response from plasmonic nanoparticles driven by Gaussian pulses with fixed duration and different central frequencies. We calculate linear and SH scattering spectra as well as the SFG response. All calculations in the linear and nonlinear regimes use the surface integral equation (SIE) approach based on the T-PMCHWT formulation [41,42]. Such a formulation has been proven to be very accurate for scattering calculations from isolated structures [43]. We also use the eigenmode analysis to confirm the origin of the scattering peaks found in the linear and nonlinear regimes [44]. This allows finding the eigenfrequencies of the modes along with the corresponding surface charge distribution. The analysis of the far-field radiation pattern for linear, SH and SFG signals, is performed by first calculating electric field values on an imaginary 20 µm diameter sphere and then by obtaining the Mie coefficients by projecting the scattered fields on the vector spherical harmonic (VSH) basis [45].
To calculate the second order time-domain nonlinear response, one needs to take into account the SFG process. For materials, this process combines two photons at the fundamental frequencies ω 1 and ω 2 to create the nonlocal-bulk (depending on the gradient of electric field, appearing inside the volume of the nanoparticle) and the local-surface (appearing just above the surface of the nanoparticle) polarization sources at the sum frequency Ω = ω 1 + ω 2 [46]. Hereafter we assume that the nonlocal-bulk SFG polarization is much weaker than the local-surface SFG. This assumption is supported by two beams sum frequency measurements on plasmonic metals [47]. Note that for the SHG, the phase matching condition is not required because the nanoparticle is much smaller than the excitation wavelength [48]. A surface tensor susceptibility χ (2) relates the nonlinear polarization sources to the magnitude of the fundamental field at the surface of the nanoparticle. It was shown that the χ (2) n,nn component of the tensor is dominant for metallic nanoparticles, the subscript n referring to the normal of the surface component [49]. Thus, we calculate the sum frequency nonlinear polarization sources just above the surface with the scalar equation: Here, E n (r, ω) is the normal component of the electric field just below the surface at the position r. We also consider the frequency range far away from the electronic resonances of the material. Additionally, we assume that the second-order process is parametric. It leads to the real valued and frequency independent value of χ (2) n,nn [4]. Without loss of generality we set χ (2) n,nn = 1. Equation (1) is an auto-convolution of the fundamental electric field spectrum, which gives nonzero contributions at SH frequencies Ω. Indeed, in time domain, the second order nonlinearity leads to a squaring of the time-dependent field values. Due to the convolution theorem, this self-multiplication in time domain translates into an auto-convolution in frequency domain. The SFG field resulting from a fundamental Gaussian pulse can be obtained by first finding all nonlinear polarization sources using Eq. (1) and then finding the SFG fields with the SIE calculations [20]. The time dynamics of the system is then obtained with the inverse Fourier transform of the frequency signal. We refer the interested reader to [38] for additional details. In [38] it was mentioned that for a set of N different equally spaced frequencies utilized to model the Gaussian excitation pulse, (N 2 + N)/2 nonlinear computations were required to account for all possible combinations of the SFG processes. Here, we propose to use the SFG's matrix linear properties to reduce the computer processing time of the SFG dynamics. Indeed, the discretization of the surface of the nanoparticle utilized to find the SFG field at different positions of the surface mesh implies that the scattered field at frequency Ω can be expressed as a vector x (Ω) . The surface polarization excitation can also be represented as a vector y (Ω) , which express the SFG polarization sources P (2) n (r, Ω), Eq. (1). These two vectors are related via the SFG matrix C (Ω) [42]: The brute force way to compute the SFG dynamics is to consider all possible combinations of (ω k , ω l ) present in the pulse and find the nonlinear polarization P (2) n (r, ω k + ω l ) for all of them. Then for each pair of frequencies (ω k , ω l ), leading to the frequency Ω = ω k + ω l , the matrix C (Ω) should be found and inverted to find x (Ω) . Such a process requires, indeed, (N 2 + N)/2 different matrices C (Ω) to be filled and inverted because the SFG processes (ω k , ω l ) and (ω l , ω k ) are identical [38]. However, the values of the elements of the SFG matrix depend only on the frequency Ω and are independent on ω k and ω l ; only the excitation vector depends directly on the specific fundamental frequencies. It stems from the fact that the matrix C (Ω) is effectively the discretization of an integral operator relating the field and currents on the surface of the scatterer through the Green's function, at one frequency. Thus, for arbitrary positive frequencies ω 1 , ω 2 , ω 3 , ω 4 such that ω 1 + ω 2 = ω 3 + ω 4 , the following equality holds: In terms of linear algebra, C (Ω) can be viewed as a linear operator for a particular frequency Ω, which has the following property: for any two vectors x (Ω) 1 and x (Ω) 2 , Consequently, the time-consuming calculation of the SFG matrix for all frequency pairs Ω = ω n + ω m is not required. Instead, all excitation vectors y (Ω) i are calculated first. Those, which have the same sum frequency Ω are summed up to form y (Ω) . Then, the SFG matrix is calculated for this remaining set of frequencies, which requires only 2N + 1 SIE matrix calculations.
Silver was chosen for our simulations due to its low losses in the optical range in comparison to gold, leading thus to longer mode lifetime. For the SH and SFG computations we use experimental permittivity data of silver from [50]. For the eigenmode analysis, on the other hand, such experimental data is not suitable, since the analytical continuation of permittivity function in the complex plane is required; this is done with a Drude model that fits the experimental values from [50] in the range 1.0 < ω < 3.0 eV, using the following model parameters: plasma frequency ω p = 9.3 eV, damping constant γ = 0.03 eV and permittivity ε ∞ = 4.3 at ω → ∞. We considered water as a background with refractive index n = 1.33.
We also provide open access data containing results of our SIE calculations and the codes for vector spherical harmonic analysis [51].

20 nm sphere
In this section, we study the SFG produced by a 20 nm diameter silver sphere, illuminated by a Gaussian pulse. Before studying the dynamics under pulsed illumination, we first look at the linear and SH scattering under monochromatic plane wave illumination. As mentioned in the introduction, only frequency doubling is considered at this stage for the nonlinear calculation. In the linear scattering cross section spectrum, Fig. 1(a), only one peak is resolved. Please note that in the following the plots are normalized to the maximum values. In order to understand the origin of the peak and find the position of the resonances, the multipolar decomposition was performed; this will also be done for the SH and SFG spectra, as well as for the time dependent nonlinear signals. Two peaks are observed at 3.2 eV, respectively 3.38 eV, in the multipolar decomposition, Fig. 1(a), and found to be caused by the excitation of the longitudinal dipolar mode (LD), respectively the transversal quadrupolar mode (TQ). Although the amplitude of the TQ resonance is negligible in the linear regime, this resonance is needed to explain some features of the SH signal. The SH scattering spectrum, Fig. 1(b), has peaks at 3.2 eV and at 3.36 eV that are associated to the resonance of the transversal dipolar (TD) and longitudinal quadrupolar (LQ) modes at the SH frequency. Following the standard excitation-radiation notation [11], these modes are associated with the E1 + E2→E1 and E1 + E1→E2 channels. Here, the coupling between a linear dipolar (E1) and a quadrupolar (E2) modes gives rise to a dipolar (E1) mode at SH (LD + TQ→TD) and the coupling of two linear dipolar (E1) modes leads to the formation of the quadrupolar (E2) mode at SH (LD + LD→LQ). Do note that the first channel involves the TQ mode (E2), whereas, the second channel involves the LQ mode (E2). In the case of a sphere, these two modes are degenerate [52,53]; this will not be the case for the nanorod studied later. These two channels, along with the eigenmode's surface charge distribution are schematized in the inset of Fig. 1(c). It is interesting to note that in the linear regime, only the modes which possess odd charge distribution (with respect to the propagation direction), can be excited. However, as is seen in Fig. 1(c), the modes with even charge distribution appear at SH. This happens because the second harmonic process essentially squares the electric field, such that the polarization sources are distributed symmetrically with respect to the propagation direction. Thus, only modes with even charge distribution can be excited at SH stage.
With the domination of the dipolar response in the linear regime, one could expect the channel E1 + E1→E2 to be the dominant one over the entire SH spectrum. However, it appears to be the case only in a short energy interval, from 3.35 eV to 3.4 eV, Fig. 1(b), the second channel that involves only once the dipolar mode being dominant in the rest of the spectrum. To understand this, a few things need to be elaborated. A peak in the SHG spectrum can be due to a resonance either in the linear regime (at half of the energy of the scattering peak), or at the second harmonic. In the present case, the linear frequency range leading to the 3.0-3.6 eV SH emission falls below all the resonances of the sphere. Thus, the two peaks in the SHG spectrum are indeed due to resonances at the second harmonic, Fig. 1(b), and the multipole nature of those peaks is given by the corresponding eigenmodes. At first, it thus seems strange that the E2 SH, which comes from twice the strong linear E1, is only dominant around the LQ resonance. For small particles, by its very nature, the quadrupolar modes radiates less efficiently than the dipolar one (as the second term in multipolar decomposition), which results in a leading dipolar term at the SH. However, as the frequency reaches the LQ mode, the magnitude of quadrupolar oscillations becomes sufficient to radiate effectively in the far field.
Let us now consider the dynamics of the SH scattered light under a Gaussian pulse illumination. We use a pulse consisting of plane waves with the following electric field dependence on time and space: Here ∆t is the pulse width, ω 0 is the central frequency, k is the wave vector with |k| = ω 0 /c and c is the speed of light in water. For the central frequency of the pulse, we choose here half the resonance frequency of the LQ mode: ω 0 =3.36/2 eV, with the goal to effectively enhance the quadrupolar mode in the SH emission. The width of the pulse is ∆t = 22 fs, corresponds to a spectral width of 0.071 eV. In Fig. 1(c) we present the resulting SFG scattering spectra.
Comparing the SH and SFG scattering spectra, Figs. 1(b) and 1(c), we observe that for both panels the quadrupolar mode is leading around 2ω 0 and stronger in the SFG than in the SH scattering (relatively to the dipolar component); note the different scales in these two panels. For clarity we added vertical dashed lines at ω 0 and 2ω 0 in the figures. However, quite surprisingly, in the time dynamics of the SFG response, shown in Fig. 1(d), the quadrupolar mode is rather weak and the signal is dominated by the dipolar mode. This is caused by the Gaussian pulse's tail, which excites the strong dipolar peak at 3.2 eV. Since this dipolar peak is two orders of magnitude larger than the quadrupolar one, the small overlap still leads to a dipolar-dominated SFG signal, Fig. 1(d). To get an estimate of the ratio between quadrupolar and dipolar intensity patterns, we introduce the ratio between the integrals over the SFG spectrum F(ω) for two modes A and B as For the quadrupolar and dipolar modes of the 20 nm sphere, Fig. 1(c), this ratio gives R F (E 2 , E 1 ) = 11.5%. According to Parseval's theorem, we know that the abovementioned integrals are proportional to the integrated intensity of the temporal signal. Thus, introducing the ratio between the intensity components I A(B) (t) in time domain for modes A and B R I (A, B) = ∫ +∞ −∞ I A (t)dt/∫ +∞ −∞ I B (t)dt and computing it for the quadrupolar I E2 (t) and dipolar I E1 (t) intensities from Fig. 1(d), we obtain R I (E 2 , E 1 ) = 11.5%, which is, indeed, the same value as that obtained for R F (E 2 , E 1 ).
We must also detail interesting features of the dynamics observed in Fig. 1(d). One can notice that the relative phase between E1 and E2 components changes with time. This causes constructive and destructive interferences in the far field intensity. In Visualization 1 we present the dynamics of the far-field radiation pattern at the oscillations maxima, which clearly illustrates the dipolar nature of the far-field radiation. This result is in good agreement with experiments [30,32]. In addition, we see that only a small part comes from the quadrupolar mode, appearing as a slight backward/forward asymmetry of the radiation pattern. From the visualization, it is clearly seen that the radiation pattern changes its shape due to the above-mentioned destructive/constructive interference of strong dipolar and relatively weak quadrupolar modes. This effect can be explained by the fact that the frequencies of the two modes are different, Fig. 1(c), resulting in a beating effect in the total intensity. This effect is indeed only observable for a spectrally broad excitation.

40 nm sphere
Let us now consider a larger silver sphere with a diameter of 40 nm. The linear and SH scattering spectra are presented in Figs. 2(a) and 2(b), respectively. As expected, the position of the resonances slightly redshifts compared to the case of the 20 nm sphere, Figs. 1(a) and 1(b). This effect is due to the increased distance needed to separate the charges on the nanoparticle surface that forms the eigenmodes; a larger distance leading to smaller electrostatic restoring forces and hence interaction energy. At the same time, we observe in the SH scattering spectrum, Fig. 2(b), that the relative ratio between the peak values of the quadrupolar LQ resonance (3.34 eV) and the dipolar TD resonance (3.1 eV) increased from 0.012 (for the 20 nm sphere) to 0.023 for the 40 nm one. This happens because the larger the sphere, the better the longitudinal quadrupolar mode radiates thanks to increased retardation. Thus, the channel E1 + E1→E2, with the E2 corresponding to the LQ at the SH stage, benefits from stronger retardation, which is supported by experimental observations [26,27]. The dynamics of the larger sphere is obtained with a Gaussian pulse centered at half the LQ resonance (ω 0 =3.34/2 eV). The width of the pulse remains the same, i.e. 22 fs. The SFG spectrum and the far-field intensity dynamics are presented in Figs. 2(c) and 2(d). A behavior similar to that of the 20 nm sphere is observed, with the dipolar mode dominating the dynamics in spite of the comparable spectral strength of the quadrupolar mode at 2ω 0 . Again, coherence between time and spectral domains for the ratio of the dipolar and quadrupolar strengths is verified: R F (E 2 , E 1 ) = R I (E 2 , E 1 ) = 18.8%. We note that compared to the small sphere, here the quadrupolar contribution is noticeably enhanced with respect to the dipolar one. In Visualization 2 we present the far-field intensity dynamics. At 49 fs, the quadrupolar radiation pattern is clearly visible. From Fig. 2(d) we again notice the variation of relative phase between E1 and E2 modes due to the E1 and E2 components being maximal at different, well-separated, frequencies. From Visualization 2, the directivity flip of the SH radiation pattern after approximately 63 fs is also the manifestation of this effect.
It is interesting to note that a pure dipolar SH radiation pattern can also be obtained when this sphere is exited with the Gaussian pulse at half the dipolar resonance (ω 0 =3.1/2 eV), with the same width as previously. The corresponding SFG spectrum is shown in Fig. 3(a) and clearly exhibits a pure dipolar component. The intensity dynamics is given in Fig. 3(b), where we observe only a dipolar component. This is attributed to the E1 + E2→E1 channel; the quadrupolar component being negligible. This behavior is also consistent with experimental results obtained with 50 nm gold spheres [32]. The observation of a single maximum in the dynamics of the SFG in Fig. 3(d) where only the dipolar resonance is excited at the SH, leads to the conclusion that the local minima of intensity observed around 50 fs for the 20 nm and 40 nm (excited at half the quadrupolar resonance) is due to interferences between the two types of radiation [33,38].

Nanorod
In the previous sections, we have been able to increase the SH quadrupolar contribution to the far-field intensity by increasing the sphere diameter; i.e. increasing retardation. However, in doing so, we have also detuned the resonances positions. A better control over the individual modes positions can be achieved by considering a nanorod instead of a sphere. Let us study the dynamic response of a 40×40×120 nm 3 silver nanorod with hemispherical ends, excited with a wave propagating along its longest axis, as sketched in Fig. 5(d). In order to classify the different peaks in linear and SH scattering we perform an eigenmode analysis of the nanorod [44], which gives the eigenfrequencies and the corresponding eigencharge distributions. When modes are excited at their resonance frequency, they scatter light effectively, leading to peaks in the scattering spectrum. The surface charge distributions, attributed to the different eigenmodes of interest here, are presented in Fig. 4. With the excitation configuration chosen here, only the modes with odd charge distribution with respect to the propagation axis are excited in linear scattering. On the other hand, the modes with even charge distribution can only be obtained at the SH stage (a tilted illumination with respect to the long axis of the nanorod would eliminate this limitation). Fig. 4. Surface eigencharge distribution for the modes, excited in linear, panels (a-c), and in SH regimes, panels (d-f). The channel of (4,1) mode formation (g). The first pair of integers in short captions represent the mode definition, the second row stands for the real part of the eigenvalue of the mode. The corresponding rows represent the far-field multipole contributions obtained by vector spherical harmonic analysis. The notation of multipoles is adapted from [11]. (g) The (4,1) mode formation channel at the SH regime. Up to now we used the notation for the modes from [11], indicating E1 and E2 as dipolar, respectively quadrupolar modes. For the nanorod, we have higher order modes, which are more conveniently classified by the number of surface charge extrema (islands). Henceforth, each mode is classified by a pair of two numbers (m,n), which represent the number of charge islands distributed along the long and the short axes of the nanorod. Thus, the dipolar modes with the dipole moment aligned perpendicular or along the long axis of the nanorod are termed (1,2) or (2,1), respectively. For example, the E1 + E2→E1 channel that we studied for the spheres would correspond to (1,2)+(2,2)→(2,1) with that notation. In this case the mode (2,1) (the longitudinal dipole) exists, but its resonance energy 1.74 eV lies out of the range of interest. The next mode with the same symmetry that can be excited is (4,1) with the channel (1,2)+(4,2)→(4,1), see Figs. 4(d) and 4(g).
Before proceeding to the classification of different peaks in scattering cross section spectra by associating them to eigenmodes, we would like to point out that the multipoles coefficients obtained by VSH analysis do not allow to unambiguously make conclusions about the nature of the corresponding surface charge distribution, which is responsible for a given far-field intensity distribution. For example, while the two charged islands of the E1 mode, Fig. 4(a), radiate like a dipole and the four charged islands of E2 mode (2,2), Fig. 4(b), radiate like a quadrupole, but quite surprisingly, the (3,2) mode, Fig. 4(c), radiates like a dipole. Indeed, higher order modes do not exhibit radiation patterns that one could expect from a radio antenna analogy. This effect has already been observed in plasmonics [54], where it was found that modes, higher than quadrupolar, usually give lower-order VSH components (electric or magnetic dipoles and quadrupoles). For instance, for an even value m, the radiation pattern of the (m,1) modes (for example, (4,1) or (6,1), Figs. 4(d) and 4(f)) is mostly dipolar. On the other hand, for odd values m > 1, the radiation pattern of the (m,1) modes (for example, (5,1), Fig. 4(e)) is mostly quadrupolar.
Having explained the radiation features of the modes, let us proceed to the classification of the scattering peaks. In the linear scattering spectrum, Fig. 5(a), the peak at 3.2 eV is associated to the (2,2) mode, Fig. 4(b). The dip in the dipolar trace E1, is associated to a Fano resonance due to the interference of the bright (1,2) and dark (2,2) modes, Figs. 4(a) and 4(c).
Four peaks are observed in the SH scattering spectrum, Fig. 5(b); looking at the eigenfrequencies and the multipolar nature of the modes in Fig. 4, we can associate each peak to the resonant excitation of one of them. The dipolar peak at 3.17 eV is attributed to the (4,1) mode, Fig. 4(d). The dipolar peak at 3.35 eV is attributed to the (6,1) mode, Fig. 4(f). The octupolar peaks are also observed, matching the positions of dipolar peaks due to the abovementioned radiation properties of high-order modes. The quadrupolar peak at 3.3 eV is attributed to the (5,1) mode, Fig. 4(e). The quadrupolar peak at 3.53 eV is attributed to the combination of higher order modes. The slight discrepancy (about 0.1 eV) between the peak values indicated in Fig. 5 and the eigenvalues obtained from the eigenmode analysis (Fig. 4) originates from the Drude model used for eigenmode analysis: this model is optimized to fit the experimental data used for scattering calculations in the range 1.0 eV < ω < 3.0 eV. For higher energies, the real parts of Drude model and Johnson and Christy data deviate as the energy increases with an underestimation of the real part of the permittivity for the Drude model. This underestimation of the permittivity is coherent with an overestimation of the eigenfrequency.
In the SH scattering spectrum, we clearly see a dominant quadrupolar component between 3.23 eV and 3.57 eV. When we excite the nanorod with a Gaussian pulse centered at 3.38/2 eV with a width 0.071 eV / 22 fs, i.e. at half the center of the quadrupolar-dominated scattering window, we expect to have a significantly enhanced quadrupolar far-field component. The SFG scattering spectrum is presented in Fig. 5(c). The ratio between quadrupolar and dipolar components gives R F (E 2 , E 1 ) = 176.8%. Also an octupolar component is present in the SH signal because of the (4,1) mode (3.23 eV) and the (6,1) mode (3.48 eV), Figs. 4(d) and 4(f). Actually, the emergence of such octupolar radiation was observed experimentally in quasi-monochromatic experiments with 100 nm diameter gold spheres [32]. The ratio between integrated intensity components of octupolar and dipolar patterns is R F (E 3 , E 1 ) = 6.1%. The intensity dynamics of the SFG signal and its multipolar components are presented in Fig. 5(d), where a dominating quadrupolar radiation pattern is indeed clearly seen. The far-field radiation pattern is presented in Visualization 3. After 52 fs a clear quadrupolar pattern is observed. Despite the fact that the SFG spectrum, Fig. 5(c), is distributed along a relatively large frequency range, it is interesting to note in Fig. 5(d) that the relative phase between E1 and E2 components is approximately zero. This stems from the fact that, in contrast with the previous cases in Fig. 1(c) and Fig. 2(c), no sharp and separated peak is present in the spectra, Fig. 5(c). Thus the amplitudes of the E1 and E2 frequency components are quite close, which ensures a constant phase difference between them. One can notice from Visualization 3 that the asymmetric pattern does not experience the directivity flip that was observed in Visualizations 1 and 2.
Finally, let us discuss the dynamics of the SFG sources P (2) n (r, t) appearing above the surface of the nanorod due to the same pulsed illumination. As is seen from Fig. 5(a), the excitation at 1.69 eV of the nanorod at the linear stage is far from the resonances appearing close to 3.2 eV. However, even though not resonant, the modes in that region influence the response of the system excited at low frequency. Indeed, we have found that the dynamics of the SFG sources can be very closely approximated by a combination of the first three modes excited in the linear regime, i.e. the modes (1,2), (2,2) and (3,2), Figs. 4(a)-4(c). This is obtained by summing the complex-valued distributions of the electric field just below the surface of each eigenmodes and then squaring the resulting fields. The time dependence is realized by varying the phase of the fields. The separate eigenmodes, their sum and their sum squared are shown in Visualization 4(a)-4(e). The colors represent the magnitude and direction of normal electric field distribution associated with the eigenmodes of the sphere, blue for electric field pointing from background to metal and red for the opposite direction. For the Visualization 4(e), the electric field is squared leading to the polarization pointed always away from the nanostructure. The constructive and destructive interference of the fields associated with the eigenmodes explains the asymmetrical polarization dynamics of the nanorod. Please note that the dynamics of polarization sources does not reflect the self-consistent SFG plasmon dynamics. Instead, it represents the dynamics of the SFG driving sources.

Conclusion
We have numerically studied the linear, SH and SFG scattering for different silver nanostructures, namely 20 nm and 40 nm diameter spheres as well as a 40×40×120 nm 3 nanorod.
For small 20 and 40 nm spheres, we showed that the dipolar response is expected at the SH due to the poorly radiating nature of the quadrupolar mode. However, a quadrupolar-dominated response is possible when long pulses (above approximately 100 fs) close to the resonance of the quadrupolar mode drive a silver sphere larger than 40 nm. On the other hand, for short pulses (22 fs), the associated broad spectrum leads to the excitation of the neighboring strongly radiating dipolar mode. Thus, the SH far-field signal is mostly dipolar for short pulse excitations. The appropriate choice of the pulse central frequency enables obtaining purely even dipolar responses, as was shown in Sec. 3.2.
A nanorod was used to demonstrate how the increased retardation effect enables using the radiation properties of higher order modes that radiate like a quadrupole at the SH, thus leading to a strong quadrupolar signal at approximately the same energy as that of the 20 nm sphere. The eigenmode analysis of the nanorod demonstrated the preponderance of the low order multipolar components in the scattering peaks.
Simulations of the nanorod response to 22 fs pulse showed major quadrupolar together with dipolar and octupolar responses. Furthermore, the asymmetric dynamics of the SFG polarization sources could be explained by the interference of three nanorod eigenmodes.
Due to different lifetimes, a situation of pure extinction of one mode while another one still radiates can take place in the signal dynamics, as was observed for the 20 nm sphere and the nanorod.
Finally, we also proposed a new computation procedure based on the linear properties of the sum frequency matrix to significantly reduce the SFG calculation time.