We investigate second harmonic generation (SHG) from hexagonal periodic arrays of triangular nano-holes of aluminum using a self-consistent methodology based on the hydrodynamics-Maxwell–Bloch approach. It is shown that angular polarization patterns of the far-field second harmonic response abide to threefold symmetry constraints on tensors. When a molecular layer is added to the system and its parameters are adjusted to achieve the strong coupling regime between a localized plasmon mode and molecular excitons, Rabi splitting is observed from the occurrence of both single- and two-photon transition peaks within the SHG power spectrum. It is argued that the splitting observed for both transitions results from direct two-photon transitions between lower and upper polaritonic states of the strongly coupled system. This interpretation can be accounted by a tailored three-level quantum model, with results in agreement with the unbiased numerical approach. Our results suggest that the hybrid states formed in strongly coupled systems directly contribute to the nonlinear dynamics. This opens new directions in designing THz sources and nonlinear frequency converters.

The research in plasmonics found various applications in chemistry,1,2 biology,3 and engineering.4,5 All are due to the ability of plasmonic materials to sustain plasmon–polariton resonances with a very small mode volume.6 Even though characteristic Q-factors of such modes are relatively small, large local field enhancement7 makes such systems very attractive for numerous applications, including nonlinear nano-optics8,9 and plasmon–polariton chemistry.10 The latter combines small plasmon mode volumes with molecular excitons.11 When the coupling strength of a local field and an exciton surpasses all the damping rates in a system, the strong coupling regime is reached. In such a regime, new hybrid states called upper and lower polaritons are formed.12 Under strong coupling conditions, these states have a mixed plasmon–exciton character.

While the linear optics of such systems is being actively explored, the nonlinear phenomena such as second harmonic generation (SHG) recently began attracting considerable attention.13–17 What would make nonlinear plasmonics an exciting playground in terms of fundamental as well as applied objectives? From a fundamental viewpoint, interactions of organic and inorganic systems, including molecules of various complexity with metals, challenges our abilities to properly describe such systems using conventional semi-classical methodology. At the application side, advanced sensing and bio-sensing devices relying on the sensitivity of nonlinear effects to perturbations and on the unlimited possibilities for targeted surface sensitization open-up new horizons in diagnosis and therapy.

In this manuscript, we report theoretical investigations of the nonlinear properties of a strongly coupled exciton–plasmon system. It is comprised of a planar array of metallic nanoscale cavities with a contacting molecular layer of a general nature. Whereas SHG has been chosen as the central target at this stage, other phenomena can be predicted along the same theoretical lines, such as other three- and four-wave mixing schemes, among which are third-harmonic generation, the intensity dependent index of refraction, and parametric effects including near IR and THz generation.

Our approach is based on the numerical exploitation of a non-perturbative hydrodynamic-Maxwell model. It can be therefore viewed to some extent as a “numerical experiment,” then calling, very much as in the case of an actual experiment, for simple physical interpretations of numerical findings. Interpretations of the numerical data are based on a semi-classical quantum picture relying on a minimal three-level system that appears sufficient to account for the main features obtained from simulations. The proposed analytical model allows for simple expressions that are prone to straightforward interpretations in terms of photon mediated coupling between the basis states. The purpose of our studies is to theoretically explore in a relatively generic way the potential for nonlinear optics of mixed metal–molecule structures where strong interactions between both species are expected to display new features with the potential to enhance their efficiency based on well identified parameters.

Early theoretical works in nonlinear plasmonics combined either a nonlinear material with otherwise linear plasmonic interface,18 employed ab initio calculations limited due to obvious reasons to small sizes,19 or utilized a semi-classical hydrodynamic model for conduction electrons.20 The latter has recently proved to be highly effective even on a quantitative level.21 It was recently demonstrated15 that by combining the hydrodynamic model for metal with Bloch equations for molecules in the time domain, one can truly explore nonlinear optics of strongly coupled systems. The major question that we attempt to address in this work is as follows: in a strongly coupled system with two well-defined polaritonic states, can we expect to observe the hybrid states directly participating in nonlinear dynamics? And if so, are there direct lower-to-upper polariton transitions? The positive answer to this question would lead to a large body of exciting applications ranging from THz generation to electromagnetic (EM) induced transparency and slow light.22 

In the section titled Theory, we overview the main features of the hydrodynamic-Maxwell–Bloch model. Next, we discuss our main results based on this model in terms of linear and nonlinear spectroscopic as well as polarization properties. Finally, we propose a simple analytical three-level quantum model, confront it to the former “numerical experiment,” and then draw conclusions. Additionally, we provide the  Appendix, which discusses the polarization properties of SHG based on abstract tensor and symmetry considerations as well as a more visual geometric argumentation.

We consider periodic arrays of triangular nanoholes arranged in a hexagonal lattice, as schematically depicted in Fig. 1(a). The optical response of this system is evaluated using Maxwell’s equations in the space–time domain,

Ḃ=×E,ε0Ė=1μ0×BṖpl,
(1)

where the macroscopic polarization Ppl is determined by the hydrodynamic model of conduction electrons for a given plasmonic system [our notation with the subscript “pl” (plasmon) is to distinguish between polarizations in the plasmonic system and the molecular layer as discussed below].23 The equation of motion describing Ppl combines both the continuity equation relating the polarization current and the number density, n, and the equation determining dynamics of the electron velocity field.24 Such an equation can be written in a concise form as25,26

P̈pl+γeṖpl=n0e2me*E+eme*Ṗpl×BEPpl1n0eṖplṖpl+ṖplṖpl,
(2)

where me* is the effective mass of a conduction electron, n0 is the equilibrium number density of electrons, and γe is the phenomenological damping rate. In the simulations, we parameterize the plasma frequency, Ωp=n0e2/ε0me*, the number density, and the damping rate. In this paper, we consider aluminum with the following parameters: Ωp = 10.85 eV, γe = 0.4641 eV, and n0 = 1.81 × 1029 m−3. The geometrical parameters defining the array are shown in the caption of Fig. 1. A few general remarks regarding the hydrodynamic model employed are in order. We first note that this approach is not perturbative and thus, in principle, contains all possible nonlinear phenomena. The terms leading to a nonlinear dependence of the optical response on the field intensity are the Lorentz magnetic term, Ṗpl×B, the quadrupole term, EPpl, and the convective terms [last two terms in Eq. (2)]. All the nonlinear terms were shown to be important in SHG simulations of periodic arrays of nanoholes.15 It was also demonstrated24 that the second order perturbation theory applied to Eqs. (1) and (2) can be used to quantitatively describe the SHG process in plasmonic systems.25 We also note that the hydrodynamic model coupled to Maxwell’s equations results in the SHG signal at flat metal interfaces. This is because the incident k-vector breaks the symmetry leading to SHG. All these features make the hydrodynamic model superior compared to approaches that rely on considering the second-order nonlinear polarization using phenomenological second-order susceptibility tensor.21 

FIG. 1.

Linear spectra of the periodic hexagonal array of triangular nanoholes in the Al metal film. Panel (a) shows the schematic setup of the array. Holes have a shape of equilateral triangles with a side of 230 nm. Red dashed lines indicate periodic boundaries. Panels (b) and (c) show linear transmission and reflection as functions of the incident frequency, respectively. Black lines show results for the 250 nm thick metal. Red lines are for 300 nm thickness, and blue lines are for the 350 nm thick film. The period of the array is 395 nm.

FIG. 1.

Linear spectra of the periodic hexagonal array of triangular nanoholes in the Al metal film. Panel (a) shows the schematic setup of the array. Holes have a shape of equilateral triangles with a side of 230 nm. Red dashed lines indicate periodic boundaries. Panels (b) and (c) show linear transmission and reflection as functions of the incident frequency, respectively. Black lines show results for the 250 nm thick metal. Red lines are for 300 nm thickness, and blue lines are for the 350 nm thick film. The period of the array is 395 nm.

Close modal

In addition to the metal, we also consider two-level molecules forming a thin layer on the input side of the array. The molecular layer is simulated using rate equations12 that determine the dynamics of the number density of molecules in the excited state, N2, and the corresponding macroscopic polarization of the molecular layer, Pmol,

Ṅ2+γ21N2=1ω12EṖmol,P̈mol+ΓṖmol+ω122Pmol=23ω12μ122N02N2E,
(3)

where γ21 is the population relaxation rate, Γ = γ21 + 2γd, γd is the pure dephasing rate, ω12 and μ12 are the transition frequency and the transition dipole, respectively, and N0 is the total number density of molecules. The initial conditions correspond to all molecules in the ground state. The parameters used in this work are γ21 = 4.14 meV, γd = 7.03 meV, μ12 = 10 D, and N0 = 2 × 1026 m−3, and the thickness of the layer is 20 nm. In addition, in Eq. (3), we assumed that the total number of molecules is conserved, N0 = N1 + N2.

Equations (1)(3) are discretized in space and time following the finite-difference time-domain method.15 When evaluating both linear and nonlinear responses, the system is excited by a planewave at normal incidence. We follow the conventional total field/scattered field approach27 to generate plane waves. Periodic boundary conditions are applied along the X and Y axes [see Fig. 1(a)], and convolution perfectly matched layers28 are set on top and bottom of the grid to absorb outgoing waves. The EM energy flux is detected on the input (reflection) and output (transmission) sides of the array in the far-field zone as the corresponding Poynting vector integrated over a given detection plane. In order to speed the simulations up, we employ a three-dimensional domain decomposition technique and message-passing interface subroutines. Our home-build codes are parallelized for 12 × 12 × 8 = 1152 processors with a spatial resolution of 1.5 nm and a time step of 2.5 attoseconds. Linear spectra shown in Figs. 1 and 3 are usually converged after propagating equations for 300 fs. The average execution times vary between 15 and 30 min. The nonlinear simulations are performed using a 500 fs long pump. To eliminate spurious nonlinear contributions that always occur due to spatial discretization, we apply the Blackman–Harris time window to all detected fields.19 Such a procedure results in clear power spectra with well-resolved harmonics. The execution times of our nonlinear codes are between 2 and 3.5 h. The nonlinear response is calculated the following way. The system is driven by a 500 fs strong laser pulse at a fixed pump frequency and the peak amplitude of 5 × 107 V/m. The latter is chosen such that the intensity of the second harmonic scales quadratically with the pump intensity. The signal is calculated as an outgoing EM energy flux integrated over detection planes on the input (reflection) and output (transmission) sides of the system. At each grid point along the detection planes, we compute corresponding EM field amplitudes as functions of time and space. Once the time propagation is complete, we perform Fast Fourier Transform (FFT) of all detected EM field components. The spatial integration of the corresponding Poynting vector is performed after FFT. The obtained signal is obtained as a function of the outgoing frequency. The SHG signal is extracted from the total signal by integrating over outgoing frequency in the vicinity of the second harmonic of the pump.

Figures 1(b) and 1(c) present linear transmission and reflection spectra of the nanohole array in vacuum for different thicknesses of the metal film. We note that both transmission and reflection exhibit multiple resonances with the broad mode clearly demonstrating the splitting, which is more pronounced in the reflection. The splitting is gradually disappearing with increasing thickness indicating that the physics of this is due to the coupling between surface modes on the input and output sides of the array, forming symmetric and antisymmetric states with respect to the center of the film (such modes were first discussed in Ref. 29 and then later in Ref. 30, and experimental implications of these modes and their properties were examined in Ref. 31). The broad mode with the splitting centered at 2.45 eV corresponds to a localized plasmon mode with an EM field primarily localized near the vertices of holes. The other modes observed near 3 eV are Bragg-plasmon modes. In this paper, we consider the localized plasmon mode, its coupling to molecules, and its influence on SHG. We note that the localized plasmon enhances SHG efficiency notably higher compared to Bragg plasmons, as demonstrated in Refs. 15 and 17.

Figure 2 shows angular diagrams of SHG response detected in the transmitted signal. We perform two types of detection: (1) calculating the SHG signal polarized along the pump polarization and perpendicular to it [Figs. 2(a) and 2(b)] and (2) calculating the SHG signal polarized along the X and Y axes in the fixed (laboratory) system of coordinates [see Fig. 1(a) where the X and Y axes are defined].

FIG. 2.

Angular properties of the second harmonic generation enhanced by a localized plasmon mode. Angular diagrams correspond the SHG signal detected on the output side of the array as a function of the pump polarization. Panels (a) and (b) correspond to the detector that rotates together with the pump. SHG polarized parallel to the pump is shown in panel (a). Panel (b) shows SHG polarized perpendicular to the pump. Panels (c) and (d) correspond to the stationary detector in the laboratory system of coordinates. Horizontally polarized SHG is shown in panel (c), and the vertically polarized SHG is shown in panel (d). Simulations are performed for the array with the period of 395 nm. The metal thickness is 350 nm. The pump frequency is 2.45 eV.

FIG. 2.

Angular properties of the second harmonic generation enhanced by a localized plasmon mode. Angular diagrams correspond the SHG signal detected on the output side of the array as a function of the pump polarization. Panels (a) and (b) correspond to the detector that rotates together with the pump. SHG polarized parallel to the pump is shown in panel (a). Panel (b) shows SHG polarized perpendicular to the pump. Panels (c) and (d) correspond to the stationary detector in the laboratory system of coordinates. Horizontally polarized SHG is shown in panel (c), and the vertically polarized SHG is shown in panel (d). Simulations are performed for the array with the period of 395 nm. The metal thickness is 350 nm. The pump frequency is 2.45 eV.

Close modal

Figure 2 displays two strikingly different SHG polarization patterns. Panels (a) and (b) exhibit six lobes in a sixfold symmetry configuration, whereas panels (c) and (d) exhibit four lobes in a fourfold symmetry configuration. These angular symmetry features are discussed in great detail in the  Appendix based on the tensor symmetry considerations along with a pictorial discussion of the non-intuitive fourfold emission from a threefold symmetric nonlinear response.

It is informative to examine the SHG signal as a function of the pump frequency. In general, any dispersive material with a second order nonlinearity exhibits the second-order susceptibility, which can be written as a product of the first-order susceptibility at the fundamental and second harmonic frequencies following Miller’s rule,32 

χ22ω;ω,ωχ12ωχ1ω2.
(4)

We thus should expect to observe peaks in the SHG signal corresponding to one- and two-photon emission. Figure 3 presents the conversion efficiency in transmission and reflection as a function of the pump frequency. The conversion efficiency defined as the ratio of the intensity of the second harmonic signal, I2ω, and the intensity at the fundamental frequency, Iω.

FIG. 3.

Conversion efficiency defined as a ratio of the intensity of the second harmonic and the intensity at the fundamental frequency as a function of the pump frequency. Panel (a) shows data for the transmission. Panel (b) shows the reflection. Both panels present data on a double-Y scale. SHG efficiency is shown as black (in a logarithmic scale) lines, and the linear transmission/reflection are shown as red lines. The red dashed lines correspond to transmission/reflection plotted as a function of the half-frequency. Material parameters are the same as in Fig. 2.

FIG. 3.

Conversion efficiency defined as a ratio of the intensity of the second harmonic and the intensity at the fundamental frequency as a function of the pump frequency. Panel (a) shows data for the transmission. Panel (b) shows the reflection. Both panels present data on a double-Y scale. SHG efficiency is shown as black (in a logarithmic scale) lines, and the linear transmission/reflection are shown as red lines. The red dashed lines correspond to transmission/reflection plotted as a function of the half-frequency. Material parameters are the same as in Fig. 2.

Close modal

The results presented in Fig. 3 are compatible with Miller’s rule for plasmonic systems. There are clear signatures of both one- and two-photon emission enhanced by corresponding plasmon modes of the array. As expected, the largest enhancement of the conversion efficiency is due to the localized plasmon at 2.45 eV, which is higher than that due to the Bragg plasmon at 3.1 eV by more than an order of magnitude as was reported elsewhere.15 The two-photon emission peaks also have an obvious correspondence to the plasmon modes of the array when compared with the red dashed line. Interestingly, Wood’s anomaly, i.e., the peak associated with the period of the array,33 influences the conversion efficiency as seen near 2 eV as a knee shape of the conversion efficiency. Another observation worth emphasizing is the splitting of the plasmon peak discussed in Figs. 1(b) and 1(c), which is mostly pronounced in the reflection near 2.4 eV. One can note that the two-photon emission peak at 1.2 eV also evidences a clear signature of the splitting. Interestingly, the splitting due to the hybridization of the surface modes on the input and output sides of the array has a different structure in the one-photon and two-photon peaks. Note that the low frequency mode has a higher conversion efficiency at the one-photon peak compared to the higher frequency mode. The two-photon peak has this reversed. This again reinforces the hypothesis that the local symmetry plays an important role in the SHG process.

We now add a thin molecular layer placing it on top of the array and examine its influence on the SHG. Figure 4 explores the linear optical response of such a system. The transmission with molecules being resonant with the localized plasmon mode exhibits a clear Rabi splitting indicative of strong coupling between molecular excitons and the corresponding plasmon mode. Figure 4(b) shows frequencies of the upper and lower polaritons as functions of the periodicity of the array exhibiting expected avoided crossing. The Rabi splitting observed for the array with a period of 395 nm is measured to be 179 meV. Next, we explore the SHG process in such a system under strong coupling conditions.

FIG. 4.

The inset shows a schematic energy diagram describing the strong coupling between an emitter and an optical cavity with a formation of two polaritonic bands. The main panel shows the linear optical response of a 20 nm thin molecular layer on top of the periodic array of triangular nano-holes. Panel (a) shows transmission as a function of the incident frequency. The black line shows transmission for the array with no molecules; the red line shows data for molecules with ω12 = 2.45 eV; the blue line is for ω12 = 2.48 eV; and the green line is for ω12 = 2.50 eV. Panel (b) shows frequencies of the upper (blue) and lower (red) polaritons as functions of the period of the array. The black line shows the frequency of the localized plasmon mode, and the black dashed line indicates the transition frequency of the molecules (fixed at 2.45 eV).

FIG. 4.

The inset shows a schematic energy diagram describing the strong coupling between an emitter and an optical cavity with a formation of two polaritonic bands. The main panel shows the linear optical response of a 20 nm thin molecular layer on top of the periodic array of triangular nano-holes. Panel (a) shows transmission as a function of the incident frequency. The black line shows transmission for the array with no molecules; the red line shows data for molecules with ω12 = 2.45 eV; the blue line is for ω12 = 2.48 eV; and the green line is for ω12 = 2.50 eV. Panel (b) shows frequencies of the upper (blue) and lower (red) polaritons as functions of the period of the array. The black line shows the frequency of the localized plasmon mode, and the black dashed line indicates the transition frequency of the molecules (fixed at 2.45 eV).

Close modal

The main result of this manuscript is shown in Fig. 5. Here, we plot the SHG efficiency for the molecule/metal system under strong coupling conditions and compare it directly with the results obtained without molecules. Unexpectedly, the two-photon peak exhibits a clear signature of the Rabi splitting. Previously, it was argued that the Rabi splitting observed in the second harmonic signal for strongly coupled systems can be explained by expanding the two coupled oscillator model.14,16With molecules treated as simple two-level emitters with no diagonal elements of the transition dipole, the only source of the second harmonic is the metal. The model based on a strongly driven oscillator with a second order nonlinearity (plasmons) coupled to a linear oscillator (molecular excitons) results in the modified second-order susceptibility (5).14 Only the second term in (5) carries the signature of the strong coupling,

χ21ωSPP2+iγSPP2ω2ω2×ωexc2+iγexcωω2g4ωSPP2+iγSPPωω2ωexc2+iγexcωω22,
(5)

where g is the coupling strength between the plasmon (SPP) and the corresponding molecular exciton (exc). Equation (5) evidently predicts that the Rabi splitting does not appear in the two-photon peak since only plasmon enhancement is observed.

FIG. 5.

Effective second-order susceptibility χ2 as a function of the pump frequency. Panel (a) shows χ2 for metal only (black line) and metal with molecules (red line). Panel (b) shows χ2 (blue circles) and the fitting using expression (6) (red line). Panels (c) and (d) show χ2 (red line) and linear absorption due to the molecular layer only (black line). Panel (c) is a zoomed-in view of the frequency region near one-photon emission, and panel (d) corresponds to the two-photon peak (linear absorption here plotted as a function of the half-frequency). The molecular transition frequency is 2.45 eV, and the periodicity is 395 nm.

FIG. 5.

Effective second-order susceptibility χ2 as a function of the pump frequency. Panel (a) shows χ2 for metal only (black line) and metal with molecules (red line). Panel (b) shows χ2 (blue circles) and the fitting using expression (6) (red line). Panels (c) and (d) show χ2 (red line) and linear absorption due to the molecular layer only (black line). Panel (c) is a zoomed-in view of the frequency region near one-photon emission, and panel (d) corresponds to the two-photon peak (linear absorption here plotted as a function of the half-frequency). The molecular transition frequency is 2.45 eV, and the periodicity is 395 nm.

Close modal

In contrast with the coupled oscillator model, our results shown in Fig. 5(a) indicate that the molecules also contribute to the two-photon emission. One may argue that the splitting seen near 1.2 eV may correspond to the absorption of two 1.2 eV photons by the molecules. However, the splitting observed at 1.2 eV equals exactly the half of the Rabi splitting, while the absorption linewidth for molecules is significantly smaller. This is also illustrated in Figs. 5(c) and 5(d), where we plot the linear absorption of the molecular layer right next to the corresponding one- and two-photon peaks. One can see that the linear absorption is significantly narrower and notably red-shifted compared to the observed resonances. Another interpretation of the observed splitting at 1.2 eV is thus needed to account for the 2ω splitting.

We propose a three-level quantum model in the semi-classical approximation to describe a two band plasmon–polariton spectrum close to the avoided crossing point. 0 designates the ground state, 1 denotes the lower polaritonic branch with energy ℏω1, and 2 denotes the upper one with energy ℏω2. We limit ourselves to the component χyyy2 for threefold symmetry as the only other component χyxx2 can be obtained via χyyy2+χyxx2=0. We also assume that the transitions are only weakly populating the excited states, making the population dominantly in the ground state. In the density matrix formulation,32 the only remaining diagonal term is ρ000=1, thus neglecting excited state populations. In the following, we shall omit the Cartesian indices, which will all be implicitly taken as yyy.

We can now proceed to adopt the general expressions given [page 174, Eq. (3.6.16) in Ref. 32] in the case of a non-dipolar three-level system. For an N-level system, we have

χ(2)(2ω;ω,ω)=N22n,m(a+b),a=μ0nμnmμm0(ωn2ωiγn)(ωmωiγm),b=μ0nμnmμm0(ωnm+2ω+iγn)(ωmωiγm),
(6)

where each term in a and b corresponds to a three-photon transition from the ground state and back, with any pair of intermediate excited states chosen among the N states.

Reduction to a three-level system limits the summation to n = 1, 2, likewise for m. Moreover, threefold symmetry and the absence of dipolar properties in the molecular layer lead to the cancellation of all participating dipoles in the summation, resulting for the case of a three-level system in μ00μ11μ22 = 0. The a and b terms take then the following simplified expressions:

a=μ01μ12μ20(ω12ωiγ1)(ω2ωiγ2)+μ02μ21μ10(ω22ωiγ2)(ω1ωiγ1),b=μ01μ12μ20(ω12+2ω+iγ12)(ω2ωiγ2)+μ02μ21μ10(ω21+2ω+iγ21)(ω1ωiγ1),
(7)

where ω12 = ω1ω2 < 0 and ω21 = ω2ω1 > 0.

The two denominators of the term clearly exhibit two- and one-photon resonances, respectively, at ω1 and ω2 in agreement with the a priori non-perturbative fully vectorial hydrodynamic-Maxwell–Bloch model. This behavior corresponds to one- and two-photon resonant second harmonic generation processes for a pair of two-level systems, namely, 0,1 and 0,2. Using simplified Eq. (7), we performed fitting of the second-order susceptibility to numerical results. The result is shown in Fig. 5(b). Although we note that the fitting is obviously not perfect, the model adequately describes the splitting observed at one- and two-photon peaks.

The four possible corresponding schemes are sketched in Fig. 6 under (a1,2) and (b1,2). Following conventions in Fig. 6, ωa1 and ωb1 are in one-photon resonance with levels 1 and 2, whereas 2ωa2 and 2ωb2 are in two-photon resonance with the same levels. In contrast with expression a, the two denominators in b exhibit only one-photon resonances, respectively, at both ω1 and ω2.

FIG. 6.

Transition schemes for SHG and DFG processes within a three-level quantum system as discussed in the text. Schemes (a1) and (b1) correspond to one photon resonances in SHG whereas (a2) and (b2) depict two-photon resonances. Scheme (c) describes a difference frequency mixing process, which sustains parametric emission or difference frequency generation depending on the experimental configuration (namely, only one incoming beam at frequency ω02 for parametric processes and two beams, respectively, at ω01 and ω02 for the difference mixing process, which can also be considered as a stimulated parametric processes). Scheme (c) is seen to provide the benefit of a doubly resonant process.

FIG. 6.

Transition schemes for SHG and DFG processes within a three-level quantum system as discussed in the text. Schemes (a1) and (b1) correspond to one photon resonances in SHG whereas (a2) and (b2) depict two-photon resonances. Scheme (c) describes a difference frequency mixing process, which sustains parametric emission or difference frequency generation depending on the experimental configuration (namely, only one incoming beam at frequency ω02 for parametric processes and two beams, respectively, at ω01 and ω02 for the difference mixing process, which can also be considered as a stimulated parametric processes). Scheme (c) is seen to provide the benefit of a doubly resonant process.

Close modal

The b term does not exhibit resonance in the optical regime since ω21ω, 2ω. However, a modified b term could be of interest in the context of a possible difference-frequency generation (DFG) process ωω,Ω=ωω, inverse to the corresponding ω,Ωω+Ω sum-frequency generation (SFG) process, whereby a coherent beam at frequency Ω in the THz regime could be generated in resonance with the ω21 transition. The corresponding susceptibility can be evaluated in two steps: the first one by way of generalizing the SHG susceptibility as detailed below onto SFGωp,ωqωp+ωq and then replacing ωp by −ωp corresponding to DFG.

For SFG, we have

χ(2)(ωpωq;ωp,ωq)=N22(a1+a2+b1+b2),a1=μ01μ12μ20(ω2ωpωqiγ2)(ω1ωpiγ1)+μ02μ21μ10(ω1ωpωqiγ1)(ω2ωpiγ2),a2=μ01μ12μ20(ω2ωpωqiγ2)(ω1ωqiγ1)+μ02μ21μ10(ω1ωpωqiγ1)(ω2ωqiγ2),b1=μ01μ12μ20(ω21+ωp+ωq+iγ21)(ω1ωpiγ1)+μ02μ21μ10(ω21+ωp+ωq+iγ21)(ω2ωpiγ2),b2=μ01μ12μ20(ω21+ωp+ωq+iγ21)(ω1ωqiγ1)+μ02μ21μ10(ω21+ωp+ωq+iγ21)(ω2ωqiγ2).
(8)

For DFG, we obtain

χ(2)(ωpωq;ωp,ωq)=N22(α1+α2+β1+β2),α1=μ01μ12μ20(ω2+ωpωqiγ2)(ω1+ωpiγ1)+μ02μ21μ10(ω1+ωpωqiγ1)(ω2+ωpiγ2),α2=μ01μ12μ20(ω2+ωpωqiγ2)(ω1ωqiγ1)+μ02μ21μ10(ω1+ωpωqiγ1)(ω2ωqiγ2),β1=μ01μ12μ20(ω21ωp+ωq+iγ21)(ω1+ωpiγ1)+μ02μ21μ10(ω21ωp+ωq+iγ21)(ω2+ωpiγ2),β2=μ01μ12μ20(ω21ωp+ωq+iγ21)(ω1ωqiγ1)+μ02μ21μ10(ω21ωp+ωq+iγ21)(ω2ωqiγ2).
(9)

We further chose the incoming frequencies close to each other such that the frequency difference ωqωp be much smaller than optical frequencies, namely, ωqωpω1, ω2 and of the order of the frequency shift between the higher and lower polaritonic branches, i.e., ωqωpω21, which can span the THz domain.

Under such conditions, the α1 expression cannot support any resonances at optical frequencies as it would otherwise contradict the above assumptions. The α2 expression can be resonant, respectively, at ωq = ω1 for its first term or at ωq = ω2 for the second one.

More interesting contributions are arising from the β1 and β2 expressions, i.e., resonances at ωqωp = ω21. In addition, ωq = ω2 in the second term of β2 can lead to a second resonance. Such a doubly resonant scheme is displayed under (c) in Fig. 6. Whereas such a scheme would be probably of limited interest in a bulk medium with the benefit of a double resonance offset by the cost of corresponding absorptions for each resonance, the situation appears much more favorable for surface driven nonlinearities at the nanoscale where propagative absorption is negligible at the scale of the structure.

We note that the doubly resonant process described above could be of major interest for the efficient generation as well as detection of THz radiation, the structure studied here playing the role of a nanoscale emitter–receptor antenna in the THz range.

At a more fundamental level, such a process would also be of major spectroscopic interest as it could permit to resonantly access the otherwise elusive μ12 transition dipole that couples levels 1 and 2. In the generic octupolar molecular system described in Ref. 34, levels 1 and 2 then referred to as 1̄ are degenerate for essential symmetry reasons resulting from threefold symmetry. Our present plasmonic system is mimicking such a system with the advantage of providing a handle for tuning a strong field induced split between 1̄ and 1.

Utilizing the non-perturbative hydrodynamic-Maxwell–Bloch approach, we examined the properties of the SHG process in a system comprised of two-level emitters strongly coupled to a localized plasmon mode supported by a periodic array of nano-holes in a thin metal film. It is shown that the effective second-order susceptibility exhibits the Rabi splitting at both one- and two-photon peaks. We proposed a simple analytical three-level quantum model, which allows setting on a firm basis the following main findings:

  1. The Rabi splitting observed for both one- and two-photon transitions provides the evidence for strong field interaction between plasmonic and molecular parts of the system. Whereas the occurrence of such splitting and its sensitivity to the lattice period have been observed and reported for linear properties, it is to our knowledge a first prediction in the context of nonlinear properties. Such features are recognized in both approaches at comparable spectral positioning and intensity levels.

  2. A set of four resonances corresponding to one- and two-photon driven transitions from the ground level to the two plasmon–polariton branches is evidenced in the simulated SHG spectrum for both models.

  3. Polarization properties from the simulation properties show well-defined threefold and fourfold patterns depending on the experimental consideration. These are accounted for in simple algebraic and geometric terms in the  Appendix.

Our theoretical studies are surely not to be seen by any means as replacement for actual experiments currently under way but can be considered of a high predictive value in terms of interesting and possibly unexpected phenomena well worthy of a subsequent experimental endeavor. Our results are to be seen as a preliminary step that paves the way onto the design, nanofabrication, and optical testing of original configurations of considerable fundamental and applied interests.

This work was sponsored by the Air Force Office of Scientific Research under Grant No. FA9550-19-1-0009. The authors also acknowledge computational support from the Department of Defense High Performance Computing Modernization Program. A.S. acknowledges the ISF (Grant No. 2525/17).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

We proceed to clarify here to account for the patterns of SHG angular polarization plots (APPs) from nonlinear entities with threefold symmetry.35 Such patterns may exhibit fourfold and sixfold, depending on the experimental configuration used, namely, the respective orientations of the incoming ω polarizer vs the 2ω outgoing analyzer, together with the sample rotation scheme, either rotating in a fixed polarizer–analyzer configuration or conversely with a fixed sample in a rotating polarizer–analyzer configuration. As expected from adequate overall symmetry considerations that encompass both the intrinsic symmetry of the microscopic nonlinear entity at stake together with that of the measurement configuration at the macroscopic laboratory scale, the former microscopic symmetry can be lowered or modified by its convolution with the macroscopic level. Such a response has been already observed in many instances from optically poled polymers36 to molecular crystals37 or semiconductor nano-crystals with tetrahedral symmetry38 onto metallic holes or metallic nano-particles with triangular shapes.14,39,40

We shall further explain here the counter-intuitive difference between the threefold symmetry shape of the nonlinear entities and the fourfold symmetry of the resulting SHG angular polarization plot in relation with the chosen polarizer–sample–analyzer (PSA) experimental configuration. A different setup configuration in which the sample is rotating leads to a more intuitive threefold symmetric APP polar plot that directly reflects the symmetry of the sample.

We shall consider in the following the two configurations depicted in Fig. 7 [note that Fig. 7(b) has been demonstrated elsewhere41,42] and provide a tensor-based justification for their APPs as well as a more physically intuitive rationale. Our ultra-thin sample can be considered to a good approximation as quasi two-dimensional thereby canceling all tensor coefficients with a z out-of-plane index, thereby limiting the number of non-vanishing tensor coefficients to that with coefficients limited to combinations of the sole x and y indices. Under this approximation, only eight non-zero tensor coefficients survive, out of which four are independent coefficients of χijk2 as a result of index permutation symmetry

χyyy2,χxxx2,χxxy2=χyxx2=χxyx2,χxyy2=χyxy2=χyyx2.
(A1)

In the case of isosceles triangular symmetry, an additional mirror plane symmetry perpendicular to the symmetry axis of the triangle leads to the cancellation of additional tensor elements that are anti-symmetric with respect to the mirror plane, e.g., featuring an odd number of x or y coefficients, namely,

χxxx2=χxyy2=0.
(A2)

This further reduces the number of non-zero coefficients from eight to four coefficients, out of which only two are independent,

χyyy2,χyxx2=χxxy2=χxyx2.
(A3)

A further increase in symmetry elements from that of isosceles triangles to equilateral triangles, that is, from D2h to D3h imposes an additional relation that results from the very definition of octupolar symmetry, namely, the cancellation of all dipolar tensors attached to the entity of interest. Such symmetry induced cancellation applies to proper physical vectors such as their permanent dipole moment or adequate combinations of higher order tensor coefficients behaving as vectors under point-wise rotations. In the case of a three rank-order tensor such as χijk2 of interest here, the resulting χi2 irreducible vector component of the full χijk2 tensor is given by χi2=χijj2 (adopting Einstein’s convention for repeated index summation). This leads here to χijj2=0 = 0 for all cartesian indices I = x, y, z. Implementation of this condition in the case of a planar threefold object leads to

χyxx2=χxyx2=χxxy2=χyyy2=χ2.
(A4)

This leaves only four mutually dependent non-zero coefficients, which all relate to a single independent χ2 coefficient, the ensuing relations between the transverse incoming field, E, and the induced polarization components being

Px=P2ωex=χxyx2+χxxy2ExEy=2χ2ExEy,Py=P2ωey=χyyy2Ey2+χyxx2Ex2=χ2Ex2Ey2.
(A5)

The projection operation along the unit vectors ex and ey is the mathematical expression of our experimental configuration (a) in Fig. 7, whereby the output harmonic polarization analyzers are set along those same ex and ey orientations. This is an important remark that will further account for the fourfold symmetry of the measured APPs that convolves the symmetry of the nonlinear entity proper, that is, threefold, together with that for the chosen instrumental PSA configuration.

FIG. 7.

(a) represents the configuration with a rotating fundamental beam polarizer and fixed harmonic analyzers, respectively, along the X and Y axes; (b) has been demonstrated elsewhere,41,42 comprising co-rotating parallel linear polarizer and analyzer. In both configurations, the sample is fixed with its symmetry elements as per the central triangles, regardless of scale considerations.

FIG. 7.

(a) represents the configuration with a rotating fundamental beam polarizer and fixed harmonic analyzers, respectively, along the X and Y axes; (b) has been demonstrated elsewhere,41,42 comprising co-rotating parallel linear polarizer and analyzer. In both configurations, the sample is fixed with its symmetry elements as per the central triangles, regardless of scale considerations.

Close modal

Taking a linearly polarized input beam with a transverse fundamental field making an angle α with the X axis, with components Ex = E0 cos α, where E0 stands for the field modulus, the second order harmonic polarization is then given by

Px2ωα=χ2E02sin2α,Py2ωα=χ2E02cos2α.
(A6)

The far-field SHG intensity is, to a multiplicative constant, the above squared expressions of the polarization components

Ix2ωαχ22Iω2sin22α,Iy2ωαχ22Iω2cos22α.
(A7)

These harmonic intensities reported here in Fig. 7 correspond to experimental configurations with a rotating input linear polarizer on a rotated mount allowing it to cover a full circular rotation in the 0 ≤ α ≤ 2π interval and two fixed linear harmonic analyzers along the X and Y axes, with the sample itself being set on a fixed mount. These expressions are in agreement with the experimentally observed as well as calculated fourfold symmetry pattern in Figs. 2(c) and 2(d) as from their respective sin22α and cos22α angular dependences. It also reflects the observed π/4 angular shift inbetween the two x- and y-analyzed APPs, namely,

Ix2ωα=Iy2ωπ4α
(A8)

is a simple consequence. Figure 8 shows specific configurations that illustrate the ensuing fourfold pattern of the APP for this configuration.

FIG. 8.

Two specific configurations in (a) (from Fig. 7) with parallel and crossed polarizer and analyzer, with the latter set along Y in both cases. These can be viewed as equivalent due to the identity of the mediating χYXX2=χYYY2=χ2 tensor coefficients in both cases (to a minus sign that disappears after squaring).

FIG. 8.

Two specific configurations in (a) (from Fig. 7) with parallel and crossed polarizer and analyzer, with the latter set along Y in both cases. These can be viewed as equivalent due to the identity of the mediating χYXX2=χYYY2=χ2 tensor coefficients in both cases (to a minus sign that disappears after squaring).

Close modal

One can further illustrate the π/4 shift when the analyzer is set along the X axis, with respect to the former configuration, as shown in Fig. 9.

FIG. 9.

Two specific configurations in (a) (from Fig. 7) with parallel and crossed polarizer and analyzer, with the latter set along Y in both cases. These can be viewed as equivalent due to the identity of the mediating χyxx2=χyyy2=χ2 tensor coefficients in both cases (to a minus sign that disappears after squaring).

FIG. 9.

Two specific configurations in (a) (from Fig. 7) with parallel and crossed polarizer and analyzer, with the latter set along Y in both cases. These can be viewed as equivalent due to the identity of the mediating χyxx2=χyyy2=χ2 tensor coefficients in both cases (to a minus sign that disappears after squaring).

Close modal

For the sake of illustrating the contribution of the PSA configuration onto the symmetry pattern of the APP in that configuration, let us now investigate an alternative one where the linearly polarized input polarizer and output analyzer are parallel and co-rotating in the 0 ≤ α ≤ 2π interval, as shown in Fig. 10, whereby the α=0,2π3,4π3 angular positions are being singled-out. It comes clearly out from inspection of Fig. 10 that those three positions (as well as any other three related to this one by an arbitrary in-plane rotation) are equivalent, thus suggesting that the APP must exhibit a threefold symmetry pattern or rather a sixfold one when squaring the outgoing field into the corresponding considering intensity. Such a configuration has been used, for example, in Refs. 41 and 42.

FIG. 10.

Co-rotating parallel polarizer and analyzer at three symmetric angular positions in configuration (b) from Fig. 7.

FIG. 10.

Co-rotating parallel polarizer and analyzer at three symmetric angular positions in configuration (b) from Fig. 7.

Close modal

The above projections along the ex and ey unit vectors become then irrelevant, and we turn now to the more appropriate projection onto the rotating unit vector u=excosα+eysinα,

Puα=P2ωuα=χ2E02sin3α,Iu2ωαχ2Iω2sin23α.
(A9)

This PSA configuration leads to the more “natural” sixfold [due to the squaring from the threefold Puα to a then sixfold Iu2ωα angular dependence].

Likewise, one can straightforwardly show that the alternative co-rotating perpendicular polarizer–analyzer configuration [Fig. 2(b)] leads to a sixfold symmetry pattern rotated by 30° with respect to the former parallel perpendicular polarizer–analyzer configuration [Fig. 2(a)].

1.
T. W.
Ebbesen
, “
Hybrid light–matter states in a molecular and material science perspective
,”
Acc. Chem. Res.
49
(
11
),
2403
2412
(
2016
).
2.
J.
Feist
,
J.
Galego
, and
F. J.
Garcia-Vidal
, “
Polaritonic chemistry with organic molecules
,”
ACS Photonics
5
(
1
),
205
216
(
2018
).
3.
N.
Jiang
,
X.
Zhuo
, and
J.
Wang
, “
Active plasmonics: Principles, structures, and applications
,”
Chem. Rev.
118
(
6
),
3054
3099
(
2018
).
4.
M. I.
Stockman
, “
Nanoplasmonics: Past, present, and glimpse into future
,”
Opt. Express
19
(
22
),
22029
22106
(
2011
).
5.
T. W.
Ebbesen
,
C.
Genet
, and
S. I.
Bozhevolnyi
, “
Surface-plasmon circuitry
,”
Phys. Today
61
(
5
),
44
50
(
2008
).
6.
J. A.
Schuller
,
E. S.
Barnard
,
W.
Cai
,
Y. C.
Jun
,
J. S.
White
, and
M. L.
Brongersma
, “
Plasmonics for extreme light concentration and manipulation
,”
Nat. Mater.
9
(
3
),
193
204
(
2010
).
7.
W. L.
Barnes
and
W. A.
Murray
, “
Plasmonic materials
,”
Adv. Mater.
19
(
22
),
3771
3782
(
2007
).
8.
M.
Kauranen
and
A. V.
Zayats
, “
Nonlinear plasmonics
,”
Nat. Photonics
6
(
11
),
737
748
(
2012
).
9.
N. C.
Panoiu
,
W. E. I.
Sha
,
D. Y.
Lei
, and
G.-C.
Li
, “
Nonlinear optics in plasmonic nanostructures
,”
J. Opt.
20
(
8
),
083001
(
2018
).
10.
R. F.
Ribeiro
,
L. A.
Martínez-Martínez
,
M.
Du
,
J.
Campos-Gonzalez-Angulo
, and
J.
Yuen-Zhou
, “
Polariton chemistry: Controlling molecular dynamics with optical cavities
,”
Chem. Sci.
9
(
30
),
6325
6339
(
2018
).
11.
P.
Törmä
and
W. L.
Barnes
, “
Strong coupling between surface plasmon polaritons and emitters: A review
,”
Rep. Prog. Phys.
78
(
1
),
013901
(
2015
).
12.
M.
Sukharev
and
A.
Nitzan
, “
Optics of exciton-plasmon nanomaterials
,”
J. Phys.: Condens. Matter
29
(
44
),
443003
(
2017
).
13.
T.
Chervy
,
J.
Xu
,
Y.
Duan
,
C.
Wang
,
L.
Mager
,
M.
Frerejean
,
J. A. W.
Münninghoff
,
P.
Tinnemans
,
J. A.
Hutchison
,
C.
Genet
,
A. E.
Rowan
,
T.
Rasing
, and
T. W.
Ebbesen
, “
High-efficiency second-harmonic generation from hybrid light-matter states
,”
Nano Lett.
16
(
12
),
7352
7356
(
2016
).
14.
E.
Drobnyh
,
R.
Pachter
, and
M.
Sukharev
, “
Harmonic generation by metal nanostructures optically coupled to two-dimensional transition-metal dichalcogenide
,”
J. Phys. Chem. C
123
(
11
),
6898
6904
(
2019
).
15.
E.
Drobnyh
and
M.
Sukharev
, “
Plasmon enhanced second harmonic generation by periodic arrays of triangular nanoholes coupled to quantum emitters
,”
J. Chem. Phys.
152
(
9
),
094706
(
2020
).
16.
C.
Li
,
X.
Lu
,
A.
Srivastava
,
S. D.
Storm
,
R.
Gelfand
,
M.
Pelton
,
M.
Sukharev
, and
H.
Harutyunyan
, “
Second harmonic generation from a single plasmonic nanorod strongly coupled to a WSe2 monolayer
,”
Nano Lett.
21
(
4
),
1599
1605
(
2021
).
17.
M.
Sukharev
,
O.
Roslyak
, and
A.
Piryatinski
, “
Second-harmonic generation in nonlinear plasmonic lattices enhanced by quantum emitter gain medium
,”
J. Chem. Phys.
154
(
8
),
084703
(
2021
).
18.
G. A.
Wurtz
,
R.
Pollard
, and
A. V.
Zayats
, “
Optical bistability in nonlinear surface-plasmon polaritonic crystals
,”
Phys. Rev. Lett.
97
(
5
),
057402
(
2006
).
19.
G.
Aguirregabiria
,
D. C.
Marinica
,
R.
Esteban
,
A. K.
Kazansky
,
J.
Aizpurua
, and
A. G.
Borisov
, “
Electric field-induced high order nonlinearity in plasmonic nanoparticles retrieved with time-dependent density functional theory
,”
ACS Photonics
4
(
3
),
613
620
(
2017
).
20.
C.
Ciracì
,
J. B.
Pendry
, and
D. R.
Smith
, “
Hydrodynamic model for plasmonics: A macroscopic approach to a microscopic problem
,”
ChemPhysChem
14
(
6
),
1109
1116
(
2013
).
21.
H.
Maekawa
,
E.
Drobnyh
,
C. A.
Lancaster
,
N.
Large
,
G. C.
Schatz
,
J. S.
Shumaker-Parry
,
M.
Sukharev
, and
N.-H.
Ge
, “
Wavelength and polarization dependence of second-harmonic responses from gold nanocrescent arrays
,”
J. Phys. Chem. C
124
(
37
),
20424
20435
(
2020
).
22.
M.
Fleischhauer
,
A.
Imamoglu
, and
J. P.
Marangos
, “
Electromagnetically induced transparency: Optics in coherent media
,”
Rev. Mod. Phys.
77
(
2
),
633
673
(
2005
).
23.
M.
Corvi
and
W. L.
Schaich
, “
Hydrodynamic-model calculation of second-harmonic generation at a metal surface
,”
Phys. Rev. B
33
(
6
),
3688
3695
(
1986
).
24.
Y.
Zeng
,
W.
Hoyer
,
J.
Liu
,
S. W.
Koch
, and
J. V.
Moloney
, “
Classical theory for second-harmonic generation from metallic nanoparticles
,”
Phys. Rev. B
79
(
23
),
235109
(
2009
).
25.
M.
Scalora
,
M. A.
Vincenti
,
D.
de Ceglia
,
V.
Roppo
,
M.
Centini
,
N.
Akozbek
, and
M. J.
Bloemer
, “
Second- and third-harmonic generation in metal-based structures
,”
Phys. Rev. A
82
(
4
),
043828
(
2010
).
26.
C.
Ciracì
,
E.
Poutrina
,
M.
Scalora
, and
D. R.
Smith
, “
Second-harmonic generation in metallic nanoparticles: Clarification of the role of the surface
,”
Phys. Rev. B
86
(
11
),
115451
(
2012
).
27.
A.
Taflove
and
S. C.
Hagness
,
Computational Electrodynamics: The Finite-Difference Time-Domain Method
, 3rd ed. (
Artech House
,
Boston
,
2005
).
28.
J. A.
Roden
and
S. D.
Gedney
, “
Convolution PML (CPML): An efficient FDTD implementation of the CFS–PML for arbitrary media
,”
Microwave Opt. Technol. Lett.
27
(
5
),
334
339
(
2000
).
29.
E. N.
Economou
, “
Surface plasmons in thin films
,”
Phys. Rev.
182
(
2
),
539
554
(
1969
).
30.
D.
Sarid
, “
Long-range surface-plasma waves on very thin metal films
,”
Phys. Rev. Lett.
47
(
26
),
1927
1930
(
1981
).
31.
W.
Mu
,
D. B.
Buchholz
,
M.
Sukharev
,
J. I.
Jang
,
R. P.
Chang
, and
J. B.
Ketterson
, “
One-dimensional long-range plasmonic-photonic structures
,”
Opt. Lett.
35
(
4
),
550
552
(
2010
).
32.
R. W.
Boyd
,
Nonlinear Optics
, 3rd ed. (
Academic Press
,
New York
,
2009
), p.
640
.
33.
A.
Hessel
and
A. A.
Oliner
, “
A new theory of wood’s anomalies on optical gratings
,”
Appl. Opt.
4
(
10
),
1275
1297
(
1965
).
34.
M.
Joffre
,
D.
Yaron
,
R. J.
Silbey
, and
J.
Zyss
, “
Second order optical nonlinearity in octupolar aromatic systems
,”
J. Chem. Phys.
97
(
8
),
5607
5615
(
1992
).
35.
J.
Zyss
, “
Molecular engineering implications of rotational invariance in quadratic nonlinear optics: From dipolar to octupolar molecules and materials
,”
J. Chem. Phys.
98
(
9
),
6583
6599
(
1993
).
36.
R.
Piron
,
S.
Brasselet
,
D.
Josse
,
J.
Zyss
,
G.
Viscardi
, and
C.
Barolo
, “
Matching molecular and optical multipoles in photoisomerizable nonlinear systems
,”
J. Opt. Soc. Am. B
22
(
6
),
1276
1282
(
2005
).
37.
V.
Le Floc’h
,
S.
Brasselet
,
J.
Zyss
,
B. R.
Cho
,
S. H.
Lee
,
S. J.
Jeon
,
M.
Cho
,
K. S.
Min
, and
M. P.
Suh
, “
High efficiency and quadratic nonlinear optical properties of a fully optimized 2D octupolar crystal characterized by nonlinear microscopy
,”
Adv. Mater.
17
(
2
),
196
200
(
2005
).
38.
M.
Zielinski
,
D.
Oron
,
D.
Chauvat
, and
J.
Zyss
, “
Second-harmonic generation from a single core/shell quantum dot
,”
Small
5
(
24
),
2835
2840
(
2009
).
39.
A.
Salomon
,
M.
Zielinski
,
R.
Kolkowski
,
J.
Zyss
, and
Y.
Prior
, “
Size and shape resonances in second harmonic generation from silver nanocavities
,”
J. Phys. Chem. C
117
(
43
),
22377
22382
(
2013
).
40.
A.
Salomon
,
Y.
Prior
,
M.
Fedoruk
,
J.
Feldmann
,
R.
Kolkowski
, and
J.
Zyss
, “
Plasmonic coupling between metallic nanocavities
,”
J. Opt.
16
(
11
),
114012
(
2014
).
41.
Y.
Li
,
Y.
Rao
,
K. F.
Mak
,
Y.
You
,
S.
Wang
,
C. R.
Dean
, and
T. F.
Heinz
, “
Probing symmetry properties of few-layer MoS2 and h-BN by optical second-harmonic generation
,”
Nano Lett.
13
(
7
),
3329
3333
(
2013
).
42.
Y. D.
Glinka
,
S.
Babakiray
,
T. A.
Johnson
,
M. B.
Holcomb
, and
D.
Lederman
, “
Resonance-type thickness dependence of optical second-harmonic generation in thin films of the topological insulator Bi2Se3
,”
Phys. Rev. B
91
(
19
),
195307
(
2015
).