Second harmonic generation in metasurfaces with multipole resonant coupling

2.1 Multipole lattice resonances

In general, the polarizability of the single particle is determined by the nanoparticle size and shape as well as material and spectral range of consideration, and thus can be tuned to resonant conditions. However, the optical properties of nanoparticle arrays are different, and the spectral position of the lattice resonance is determined not only by polarizability of the individual nanoparticle but also by lattice geometry. Specifically, lattice resonances occur at the wavelength where the real part of the inverse polarizability of the single nanoparticle cancels a lattice sum of the contributions of the multipoles of all other nanoparticles in the lattice. An effective polarizability of nanoparticles in the array can be introduced [13], [27]. Finally, the width of the lattice resonance is determined by the imaginary part of the inverse nanoparticle polarizability [26], [29], [30].

Analytical study of the scattering of one- and two-dimensional arrays of nanoparticles has been well established using the coupled-dipole approximation [31] and more recently coupled-multipole models have been developed [27], [32], [33]. The coupled-dipole model is effective for plasmonic nanoparticle array when the size of the nanoparticle is small compared to the wavelength of light and array period [29], and higher multipoles need to be introduced for larger nanoparticles. Some of the key takeaways from the analytical lattice resonance analysis are that the multipole resonances excited in each nanoparticle and their relative coupling coefficients are sensitive to structure parameters, such as lattice spacing, arrangement, and geometry, polarization of the incident light, the size and shape of nanoparticles, and their material. When the RA wavelength is red-shifted and far away from the single-particle resonance, a lattice resonance appears considerably narrower than conventional single-particle resonances in plasmonic nanoparticles.

In the visible spectral range, isolated plasmonic nanoparticles with radius larger than about 50 nm can support additional higher order multipoles, such as quadrupole resonances, which are relatively narrow and weakly couple to the light field of the plane wave. However, in ordered nanoparticle arrays, the electric quadrupole resonances are greatly enhanced due to the in-phase oscillations and inter-particle interactions, and narrow collective resonances become well-pronounced [26], [30], [34]. These conditions have been theoretically studied by extending the coupled-dipole method to solve coupled dipole-quadrupole equations where the dipole and quadrupole polarizabilities of the individual particles are derived from Mie theory and an infinite periodic array of identical nanoparticles is considered [30].

2.2 Metasurfaces and artificial magnetic response

A traditional assumption in optics has been that the interaction between electrical multipoles dominates material optical properties, whereas magnetic multipoles provide at best a negligible contribution. With the development of metamaterials and the ability to purposefully design their optical response, many types of magnetic multipole response have been found in both plasmonic and all-dielectric nanostructures [35], [36], [37], [38], [39], [40], [41]. Metasurfaces are the two-dimensional counterpart of bulk (three-dimensional) metamaterials and similarly allow for artificially designed, efficient magnetic response from a thin layer of nanostructure. The periodic nanoparticle array we consider in this work enable excitation of both electric and magnetic resonances and can serve as effective metasurface.

The need to additionally include the MD multipole to accurately model the suppression of total reflectance from plasmonic periodic arrays has been demonstrated in Ref. [26]. Importantly, the work has shown that for specific lattice geometries and illumination conditions (details below) the EQ and MD multipoles are coupled and have contributions to far-field reflection and transmission comparable to that of the induced ED. Magnetic multipoles can play a significant role in the scattering characteristics of the effective material and make an important contribution to the nonlinear optical processes [42], [43], [44], [45], [46].

2.3 Symmetry breaking in multipole lattice

It is known that SHG is inhibited in the bulk of a centrosymmetric media under the electric dipole approximation, while it is allowed at interfaces where the inverse symmetry is broken [47], [48]. Within a classical phenomenological description, the nonlinear polarization Ps(2ω) at the interfaces is given by Ps(2ω)=χ↔s(2): E(ω)E(ω). For a single sphere with a centrosymmetric material, using linearly polarized plane-wave illumination, the most efficient mechanism is a dipole second harmonic emission, which can only come from the nonlocal interactions of ED(ω)+EQ(ω)→ED(2ω) and ED(ω)+MD(ω)→ED(2ω) [49], [50], [51]. In this notation, the two terms on the left of the arrow describe the two exciting modes at fundamental frequency, and the third term refers to the second harmonic emission mode. The ED(ω)+ED(ω)→ED(2ω) excitation process is forbidden in a centrosymmetric object, while the quadrupole moment can be excited through a local interaction ED(ω)+ED(ω)→EQ(2ω). As the particle size increases, the generation of higher order multipoles, e.g., octupoles, was theoretically predicted and experimentally observed for gold nanoparticles with size of about 70 nm [22]. However, a more detailed analysis indicates that at this nanoparticle size, the octupole contribution is not substantial and accounts only for 8% changes in 100 nm nanoparticles.

Periodic arrangement of nanoparticles drastically changes the optical properties of the array, and the mechanism for SHG can be different for spheres in a lattice. Coupling of EQ and MD in the nanoparticle lattice indicates symmetry breaking and the possibility of enhanced SHG under these conditions. This situation is different from the case of a single nanosphere where multipoles are orthogonal, formed independently, and do not couple to each other. With the enhanced EQ and MD resonances in the lattice, high-order mechanisms for SHG are also possible, for example, MD(ω)+EQ(ω)→EQ(2ω), MD(ω)+EQ(ω)→MD(2ω), MD(ω)+MD(ω)→EQ(2ω), EQ(ω)+EQ(ω)→EQ(2ω). Beyond the electric dipole approximation, spatial field variation can also break the centrosymmetry and produce a nonlocal nonlinear polarization through field gradients [9], [52], [53]. In the present work, we employ the hydrodynamic model for the metal material, which considers the nonlinearity from the motion of the electrons and includes the bulk contribution naturally. We show that the presence of strong EQ–MD coupling provides significant advantages for spectrally narrow SHG in plasmonic lattices.

2.4 Numerical model

We employ full-vectorial electromagnetic simulations of the plasmonic metasurfaces using the Finite-Difference Time-Domain (FDTD) method. The FDTD method is extended to include a hydrodynamic nonlinear material model for the metal, which introduces complex nonlinear response in the nanoparticle array and self consistently accounts for all the induced multipoles. This approach extends the optical response of bulk and macroscopic metal objects beyond the standard linear Drude description of the metal properties. For subwavelength metal nanoparticles or metasurfaces composed of arrays of subwavelength nanoparticles, deviations from the standard Drude model become more pronounced due to the increased relative contribution of additional nonlinear terms in the behavior of the conduction-band electrons [54], [55], [56].

The FDTD method provides direct access to near-field properties of the nanoparticles and their arrays such as local field intensity, wave phase, as well as current and charge distributions and also the contributions of individual nonlinear terms. The electromagnetic fields are governed by the Maxwell’s equations and the conduction electrons inside the metal are approximated as a continuous electron gas. In this model, strong external-field excitation of the free-electron gas results in the generation of higher harmonics in the current, J, which acts as a source term for the generation of nonlinear harmonic fields [56].

The resulting coupled fluid-Maxwell system of equations is approximated numerically, and the fluid equations for the current density, J, are solved using a time-split, semi-implicit finite difference algorithm. Computation of electron charge density, (∇·E), or terms of the form (∇·J) at the dielectric–metal interface are inherently ambiguous due to the discontinuity of the normal electric field component at the dielectric–metal interface. A number of sophisticated surface treatments and their roles in nanoparticle SHG have been proposed to alleviate this problem [57], [58], [59]. We employ a transition layer with smoothed ion distribution between metal and the surrounding dielectric media [56]. This regularization provides continuous normal electric field so that the computation of (∇·E) is more physical and the detrimental aspects of the singularity are mitigated and provide physically sound results.

Our numerical FDTD model solves Maxwell’s equations for electromagnetic fields in the nanostructure, where the electrons inside the metal are described by a free electron gas. The motion of the electrons is governed by the continuity equation and Newton’s equation:

where n_e, u_e, q_e, m_e are the electron number density, velocity field, electron charge, and electron mass, respectively. The equation for current density J [56]:

where γ is the phenomenological damping constant, ωp=qe2n0/(ϵ0me) is the plasma frequency, n₀, ρ, ε₀ are positive ion density, charge density, and vacuum permittivity, respectively. Equation (3) summarizes the electron fluid response under the influence of the electromagnetic field. The first two terms correspond to the linear Drude response, the third term represents electric force and magnetic Lorentz force. The last term is the convection term. In results presented here we do not include the pressure term proportional to ∇P, where

is the quantum pressure. As has been shown in the earlier works [60], [61], this approximation is justified for nanoparticles with dimensions of ∼100 nm and larger, which is the case in our work. Equation (3) is solved through a time-splitting, semi-implicit finite difference method as detailed in Liu et al. [56]. In our case, the gold permittivity is fitted by ϵAu=ϵ∞−ωp2/(ω2+iγω), with ϵ∞=2.4023, ω_p = 1.2122 × 10¹⁶ rad/s, and γ = 3.9941 × 10¹⁴ rad/s. To obtain the linear response of the structure in proximity to the fundamental frequency, we illuminate the structure with a 20 fs full width at half maximum (FWHM) pulse. Throughout the work, we extract frequency signal from Fourier transformation of time-dependent signal in our FDTD modeling.

To exclude the contribution of sum and difference frequency generation of other wavelengths, we determine the nonlinear response of the structure by performing a sequence of continuous wave (CW) runs with a field amplitude of E = 10⁷ V/m across the spectral range of interests. A unit cell of the lattice, containing a single metal sphere, is simulated with periodic boundary conditions in the x- and y-directions and uniaxial perfectly matched layer (UPML) boundary condition in the z-direction, which is the propagation direction.

We define the forward and backward emitted second harmonic signal as the integral of the Poynting vector over the field monitor planes located in the far-field zones above and below the lattice plane and sum them to get the total SHG amount. Here, the integral in the forward and backward direction measures the SHG radiation in that half space. In all figures throughout the paper, SHG intensity results are normalized to the corresponding fundamental input field intensity.

3.1 Comparison of single particle and lattice linear responses

We study the surface lattice resonances of a plasmonic metasurface composed of gold nanoparticles arranged in a two-dimensional periodic array of infinite extent [26], [ 62]. Since our primary aim is to study the mechanism of SHG and contribution of nanoparticle multipole resonances in nonlinear processes, we consider an idealized case and assume the nanoparticles are suspended in a homogeneous medium with constant refractive index, which is realizable as a substrate and index-match top coating of the array and common in experimental setups (see e.g., [13], [18]).

The structure analyzed here, Figure 1, is an infinite array of gold nanospheres embedded in a medium with refractive index n extending periodically in the x- and y-directions. Throughout the paper, we consider only gold nanospheres with the radius r = 100 nm in a surrounding medium with n = 1.47.

A plane wave incident on the nanostructure propagates along the z-axis and the light is polarized along the x-axis. We refer to this as parallel polarization, and unless stated otherwise it is the light polarization for the results presented in the following sections. It will excite x-oriented ED and y-oriented MD. The interparticle spacing in the periodic array is initially chosen as p_x = 510 nm along the x-axis and p_y = 250 nm along the y-axis. The RA wavelength is given by λ=2nπ/|q|, where q is the reciprocal lattice vector of the array. We will use the notation 〈l,m〉 to label the different Rayleigh anomalies corresponding to the reciprocal lattice vector q=2π(lxˆ/px+myˆ/py). The choice of x-spacing p_x = 510 nm corresponds to the 〈1,0〉 RA appearing at the wavelength (λ〈1,0〉RA=npx=749.7 nm) while the initial choice of y-spacing p_y = 250 nm avoids coupling to the y-dipole lattice resonance responsible for the 〈0,1〉 RA appearing at the wavelength (λ〈0,1〉RA=npy=367.5 nm). For the case of incident E_x polarization, the EQ and MD multipoles of the nanoparticles are mutually coupled in the nanoparticle lattice. Collective lattice resonances of nanoparticle multipoles result in significant reduction of plasmon radiative damping, and consequently a dramatic narrowing of the plasmon resonance.

When the polarization of the incident wave is perpendicular to the direction of a periodicity of interest, i.e., y-polarized, lattice resonances have been shown to be the result of ED coupling between nanoparticles and have been extensively studied numerically and analytically using dipole approximation models mentioned earlier. Conventionally, lattice resonances are always narrow when compared to single-particle resonances of lowest order (e.g., ED). At the same time, lattice resonances that arise from dipole coupling in the lattice are relatively broad in comparison to lattice resonances due to the excitation of higher order multipoles. Under y-polarized illumination conditions with array parameters we choose in this work, a relatively broad dipole–dipole lattice resonance is excited and the ultra-narrow lattice resonances features do not appear in the reflection, transmission, or absorption spectra.

As a baseline reference, the simulated scattering and absorption cross-sections of an isolated nanoparticle are shown in Figure 2A. The isolated sphere exhibits an EQ resonance at the wavelength 550 nm and an ED resonance at the wavelength 870 nm, and the MD has negligible contribution to the total scattering and absorption cross-sections.

Figure 2:

(Color online) (A) Scattering and absorption cross-sections for a single gold sphere with r = 100 nm numerically calculated with FDTD method. Reflection, transmission, and absorption spectra for the nanoparticle array with (B) x-polarized and (C) y-polarized incident electric fields. The vertical dotted line denotes the wavelength of the 〈1,0〉 Rayleigh anomaly. Note that the spectral ranges are different and adjusted to enlarge regions containing narrow resonant features. In (B), a minimum in the reflection nearly coincides with a maximum in the transmission and both are close to the maximum of absorption, where the EQ lattice resonance is excited. In the panel C, a comparatively broad ED resonance is excited around the wavelength 950 nm with a small feature present near the Rayleigh anomaly wavelength at the wavelength 750 nm.

In contrast to the isolated sphere, the resonant behavior of an infinite rectangular lattice of gold nanospheres excited for two orthogonal polarization modes has different spectral features originating from collective effects and is shown for p_x = 510 nm and p_y = 250 nm in Figure 2B and C. The linear reflection, transmission, and absorption spectra when the electric field of the incident wave is polarized parallel to the x-direction (x-polarized) are shown in Figure 2B. Under x-polarized illumination, an EQ is the dominant multipole mode excited in each nanoparticle in the lattice. A minimum in the reflectance nearly coincides with a maximum in the transmittance, and both are close to the maximum of absorbance, where the EQ lattice resonance is excited. Figure 2C shows the spectrum when the electric field of the incident wave is y-polarized, no narrow features are present, and instead a comparatively broad ED resonance is excited around the wavelength 950 nm. Note that even under strictly y-polarization excitation, weak features associated with the related x-spacing (λ〈1,0〉RA=npx) are still evident at the wavelength about 750 nm [30].

3.2 SHG at the lattice resonance

Now let us consider the nonlinear response of the nanoparticle array. Figure 3 shows effective polarizabilities of multipoles at fundamental wavelength and the generated second harmonic. The generated second harmonic exhibits a narrow resonance peak, approximately 10 nm in width, which is commensurate to the spectral width of the peak in fundamental harmonic. This narrow resonance peak corresponds to the excitation of EQ- and MD-multipole resonance around the wavelength 750 nm. The off-resonance SHG is mainly due to the contribution of the broad ED lattice resonance at fundamental wavelengths. The asymmetric line shape of total the SHG suggests important contributions from the MD. The total second harmonic peak is slightly blue-shifted from the expected wavelength of 375 nm due to an asymmetry between the forward and backward scattered light, as shown in Figure S7.

Figure 3:

(Color online) Effective polarizabilities of (A) EQ, MD and (B) ED in the periodic nanoparticle array with lattice periods p_x = 510 nm and p_y = 250 nm. Multipoles are shown in panels A and B separately for clarity. (C) Relative strength of the SHG in the periodic array. The SHG with x-polarized incident electric field exhibits a narrow resonance peak around the wavelength 375 nm corresponding to the excitation of EQ- and MD-multipole resonance around the wavelength 750 nm. The vertical dotted line denotes the wavelength of the 〈1,0〉 and 〈2,0〉 Rayleigh anomaly. The dashed line in panel C denotes the ED contribution at fundamental wavelength to the SHG.

In Figures 4 and 5, we present both the fundamental and second harmonic near-field intensities around the peak in the array absorption spectra at λ = 750 nm and its second harmonic. Figure 4A and B shows the electric and magnetic field distributions in the xz- and the yz-planes cut through the nanoparticle center at the fundamental frequency. The corresponding generated electric fields E at the SHG wavelength in the xz- and yz-planes are shown in Figure 5. The generated electric fields we show are the differences between the simulations with and without nonlinear material response taken at the same time moment. We also obtained similar results using sequentially direct and inverse fast Fourier transforms; however the field enhancement in the obtained vector diagrams is weaker because of the use of finite number of frequency harmonics and loss of some information with harmonics cut. Additionally, we present vector plot in Figure S4 in Supplementary Information. At the fundamental harmonic, the enhanced near-field resembles quadrupole profile with high intensity lobes that exhibit long attenuation lengths (asymmetrically) above and below the lattice plane at the antinodes of the lattice standing wave [29].

Figure 4:

(Color online) (A) Electric |E| and (B) magnetic |H| field amplitudes at the fundamental wavelength. The colormap shows averaged field intensity, and white arrows show the direction of the electric field in a particular time moment. The field amplitude is normalized to the fundamental incident field.

Figure 5:

(Color online) The generated second harmonic electric field component (A) E_x and (B) E_z in the xz-plane, and components (C) E_y and (D) E_z in the yz-plane cut through the nanoparticle center. The field amplitude is normalized to the fundamental incident field.

The second harmonic near-field is elongated in the direction transverse to the lattice mode standing wave and its dipolar profile is more evident when viewed in the yz-plane, see Figure 5D. Inferring a multipole nature of the second harmonic field from a field distribution alone is problematic. In the following, we use common notations for vector harmonics N_emn and M_omn from Ref. [63]. Figure 6 shows the scattering cross-section of six different multipoles for a single spherical particle providing each multipole is excited with a plane wave at the second-harmonic wavelength. At the wavelength λ = 375 nm, the ED (N_e11), EQ (N_e12), and electric octupole (EO, N_e13) have non-negligible scattering cross-sections, and EQ is the dominant one. We see very small values for MD (M_o11), electric sedecapole (ES, N_e14) and nearly zero MQ (M_o12) multipoles. However, it is important to note that SHG is driven by the excitation of nanoparticle resonances at the fundamental wavelength, where the near field distribution is different from the plane wave, and second harmonic multipole composition may not coincide with the results of a simple plane wave excitation. Having the multipoles N_e1n and M_o1n coming from the x-polarized plane wave, one can only excite the multipoles with m = 0 or 2 in the second harmonic (see discussion in Ref. [64]). The first possible excited multipoles are: N_e01 as dipolar mode (which is also ED N_e11 but in another direction), N_e02, N_e22, and M_o22 as quadrupolar modes, and N_e03, N_e23, and M_o23 as octupolar modes. In this case, electric dipole mode is the same in magnitude but different in direction (x in fundamental frequency, and z in second harmonic). However, the higher multipoles, of quadrupolar and octupolar kind, are different for m = 1 and m = 0 or 2. Scattering cross-section is defined by Mie coefficients through the plane-wave excitation and related to only m = 1 and not necessary represent excitation for m = 0 or 2. Furthermore, the excited second-harmonic multipoles are arranged in the periodic lattice, and the collective effects may contribute to the effective polarizabilities and scattering cross-sections of nanoparticles in the lattice. However, the analytical theory of normal-angle plane-wave excitation of nanoparticle lattice is not applicable in this case as some of the SHG multipoles are excited out-of-plane. The field pattern in Figure 5 resembles EQ profile in the E_y- and EO in the E_x-field components from the cut-through plots. The E_z component is dominant in the SHG electric field and it is of dipole kind. We also note that the linewidth features of the generated second harmonic are more influenced by the linear multipole decomposition at the fundamental as presented in Figure 3. This indicates that the multipole composition at the fundamental wavelength dominates the character of the second harmonic emission from the lattice.

Figure 6:

(Color online) Scattering cross section σ_sca calculated with Mie theory. EO denotes electric octupole, and ES denotes electric sedecapole.

We note that the field profile at the second harmonic also correlates more strongly with the magnetic field at the fundamental, Figure 4B, than the electric field at the fundamental, Figure 4A. We note that the dominant nonlinear contribution comes from the Lorentz force term, J × B, which can be considered a bulk contribution due to the presence of the magnetic field inside the metal sphere. We provide more details on the SHG contribution of the various nonlinear terms in Equation (3) in the Supplementary material.

In terms of the second harmonic conversion efficiency, it is difficult to provide a direct comparison to other geometries (e.g., nanoparticles of more complex shape, like the commonly used U-shape) due to the problem of finding suitable nanostructure parameters which have comparable resonance wavelength and exhibit similar multipole modes and EQ–MD coupling. The spectral range under consideration defines not only the ratio between the wavelength and characteristic nanostructure dimensions, but also material response such as permittivity value and surface parameters included in the hydrodynamic model and numerical scheme. As the optical properties of nanoparticles are defined by their complex-value polarizability tensors, matching these polarizabilities for several multipoles at the same time would require an independent work of numerical optimization and a broad parameter search, which is beyond the scope of our current work.

As for the comparison of plasmonic and all-dielectric nanostructures, the former in general provide better light localization on subwavelength nanoparticle and result in resonances of a higher quality factor because of the low radiative losses. In turn, the latter has lower non-radiative losses, support magnetic resonances with simple nanoparticle shape, and spans over different material platforms, including well-established silicon and III–V nanofabrication processes. Changing the spacing between nanoparticles in a periodic lattice provides the possibility to control nanostructure resonances over a broad range and tune their linewidth on demand. Thus, both plasmonic and all-dielectric nanoparticle array can be used to achieve comparable polarizability values and resonance quality [65], and in this sense, comparably enhance the nonlinear response in the nanostructure.

3.3 Multipole effect on SHG

While the FDTD with the nonlinear hydrodynamic model can provide near field profiles and contributions of the various nonlinear terms, in the following section we additionally employ an analytical approach in order to disentangle the individual multipole contributions to SHG. Using the analytic dipole–quadrupole approximation model we determine the linear ED, EQ, and MD polarizabilities and the linear reflection, transmission and absorption linear spectra. It has been shown previously that the spacing in the perpendicular direction p_y affects the ED lattice resonance while the spacing in the parallel direction p_x is responsible for the EQ and MD resonances [30], [66]. Thus, the multipole contributions to SHG can be understood by examining the variation of the respective linear multipoles polarizabilities under lattice spacing variation.

To study the ED contributions to the SHG we fix the interparticle spacing in the x-direction at p_x = 510 nm, which keeps the x-direction RA 〈1,0〉 at the wavelength about 750 nm, and we scan the y-direction spacing over a range from 250 to 550 nm. Figure 7A–C shows the ED, EQ, and MD polarizabilities under different lattice spacing. The dashed lines show the 〈1,0〉 and 〈0,1〉 RA wavelengths. For all values of p_y, the EQ and MD resonances coincide and stay locked at the RA wavelength λRA=npx=749.7 nm. On the other hand, for values of p_y above 460 nm, the ED peak begins to broaden and shift to longer wavelengths and transitions to an ED lattice resonance aligning with the 〈0,1〉 RA wavelength. Figure 7D–F illustrates the effect of varying p_y on the linear spectrum. Two separate regimes can be clearly identified: for p_y spacing smaller than 460 nm, the spectra are dominated by the 〈1,0〉 RA, and for larger p_y spacing, the spectra are dominated by the 〈0,1〉 RA. As the p_y spacing increases beyond 460 nm, the lattice resonance mode changes from an EQ- and MD-multipole dominated resonance to the hybrid mode and then to an ED-dominated resonance, which changes its spectral position together with 〈0,1〉 RA. Finally, for p_y > 510 nm, the wavelength of interest 750 nm lies below the diffraction limit and light scattering is not zero-order directional anymore. In this regime, the intensity of the electromagnetic fields on the particle surface drops significantly, and very little field localization is observed on the nanoparticle. We provide more details in Supplementary Information.

Figure 7:

(Color online) Effective polarizabilities of (A) ED, (B) EQ, and (C) MD nanoparticle multipoles in the periodic nanoparticle array with fixed lattice periods p_x = 510 nm and variant p_y. The calculations are based on the analytical dipole-quadrupole approximation model. (D) Reflection, (E) transmission and (F) absorption linear spectra in periodic array with fixed lattice periods p_x = 510 nm and variant p_y. The colorbar is the same for (D)–(F) panels. The dashed lines show the 〈1,0〉 and 〈0,1〉 Rayleigh anomaly wavelengths.

The total SHG for different p_y spacing versus wavelength is shown in Figure 8. As p_y increases, initially from 250 nm (see Figure 3) to p_y = 355 nm, the SHG corresponding to the narrow EQ- and MD-multipole resonance near the wavelength 375 nm initially grows in magnitude and begins to exhibit a slight elbow on the long wavelength side due to increasing impact of the ED. As the p_y spacing further increases to 405 nm, the peak has split into a doublet. This trend continues at larger spacing until eventually a distinct and spectrally broad ED resonance emerges on the long wavelength side while the sharp EQ- and MD-multipole spectral peak has drops significantly in amplitude. This trend is mirrored in the absorption in the linear spectrum shown in Figure 7, in which the 〈1,0〉 RA wavelength no longer has the strongest absorption in the linear spectrum. Beyond 460 nm the peak SHG tracks the 〈2,0〉 RA and is proportional to the magnitude of the ED polarizability. As the ED transitions to a pure ED lattice resonance, without the presence of EQ- and MD-multipole resonance, the corresponding SHG becomes spectrally broad and shifts to longer wavelength commensurate with the ED at the fundamental.

Figure 8:

(Color online) Intensity of the second harmonic signal generated for different p_y spacing in colormap plots. The dashed lines show the 〈2,0〉 and 〈0,2〉 Rayleigh anomaly wavelengths.

To study the EQ- and MD-multipole contributions to the SHG strength, we vary the x-spacing in the nanoparticle array. With the lattice period p_y fixed at 250 nm and p_x varied from 400 to 600 nm, the ED, EQ, and MD resonances closely track the 〈1,0〉 RA as shown in top row in Figure 9. As p_x spacing increases, ED resonances increase slightly while EQ- and MD-multipole resonance decreases in magnitude. The ED and EQ resonance peaks for a single particle occur at wavelengths around 900 and 550 nm, respectively (see Figure 2). As we increase p_x the EQ-induced lattice resonance moves to longer wavelength further away from EQ resonance peak of a single particle, consequently the EQ resonance weakens. Simultaneously the ED resonance strengthens as it approaches the single particle at 900 nm and also starts to interact with RA at a longer wavelength. Figure 9D–F shows the linear spectra of the lattice for different p_x spacing. Along the 〈1,0〉 RA the EQ- and MD-multipole resonance weakens and the absorption around the RA wavelength is reduced.

Figure 9:

(Color online) Effective polarizabilities of (A) ED, (B) EQ, and (C) MD nanoparticle multipoles in the periodic array with fixed lattice periods p_y = 250 nm and variant p_x. (D) Reflection, (E) transmission, and (F) absorption linear spectra in periodic array with fixed lattice periods p_y = 250 nm and variant p_x. The dashed lines show the 〈1,0〉 Rayleigh anomaly wavelength for each lattice spacing.

The total SHG for different p_x spacing versus wavelength is shown in Figure 10. As p_x increases, the SHG resonance decreases in intensity and tends towards the broader ED resonance. The reduction in SHG peak intensity is consistent with the weakening of the EQ- and MD-multipole resonance (see Figure 3C). The blue-shifting of the total SHG peak is due to the sudden drop-off of backward SHG scattering at the 〈2,0〉 RA wavelength, as shown in Figure S11, since the generated multipoles at 〈2,0〉 RA wavelength interfere destructively in the backward direction.

Figure 10:

(Color online) Intensity of the second harmonic signal generated for different p_x spacing. The dashed lines show the 〈2,0〉 Rayleigh anomaly wavelengths.