Introduction

Modern photonic devices rely on nonlinear optical effects to carry out functionalities such as parametric tuning of the laser light spectrum, generation of ultra-short pulses, all-optical signal processing, and ultrafast switching1,2. Furthermore, researchers are constantly extending the reach of nonlinear optical phenomena to a growing variety of practical applications, such as nonlinear microscopy3, ultrasensitive-optical shape characterization4, monitoring the processes in chemical synthesis of nanostructures5, and nonlinear plasmonic sensing6.

The inherent weakness of optical nonlinearities however has strongly hindered the chip integration of all-optical devices. Conventionally, to enhance the intrinsic nonlinear response of natural materials, one has to rely on phase-matching or quasi-phase-matching techniques, which require relatively long propagation distances. Distances can be shortened to several microns by employing photonic band-gap crystals, which is in general detrimental for the band width. A different route to increase nonlinear optical response can be achieved through plasmonic effects—collective oscillations of conduction electrons—whose near-field enhancement properties can effectively increase the nonlinear susceptibilities. Plasmons are not limited by diffraction and allow to focus the light in regions where dielectric nonlinearities can be embedded, or such that metallic nonlinearities themselves can be enhanced. The advantage of such approach is that one can effectively increase local nonlinear susceptibilities in nanoscale devices. The drawback, however, is that it is hard to combine building-up mechanisms that allow to increase the overall conversion efficiencies. In order to overcome this limitation then it is necessary to artificially increase nonlinear susceptibilities by engineering optical properties at the quantum level. At THz frequencies, for example, this can be obtained using semiconductor quantum wells which are known to offer the largest nonlinear susceptibilities7. Increasing the carrier density, however, causes a transition from size-quantization to the classical regime of plasmon oscillations, in the visible range of frequencies.

Since the pioneering works on nonlinear optics8,9,10,11, second-harmonic generation (SHG) has been continuously studied both experimentally12,13,14,15,16,17,18,19,20,21,22,23 and theoretically24,25,26,27,28,29,30,31,32,33,34,35,36. In centrosymmetric materials, such as noble metals, second-order nonlinearities in the bulk region are intrinsically forbidden3; nevertheless, the symmetry can be broken at the material surface, giving rise to a surface SHG mechanism37,38. Typically, the SHG surface process occurs within a very thin layer of the order of few angstroms (a few Thomas–Fermi screening lengths) at the surface where the induced charges are confined. At such near-atomic length scales, classical electrodynamics fails to address the microscopic details and consideration of nonlocal and quantum mechanical effects may become crucial for accurately characterizing optical behavior of a metallic system.

In this context, the simplest tool dealing with nonlocal electron pressure is the Thomas–Fermi hydrodynamic theory (TFHT)39,40,41,42,43,44,45. The TFHT is often combined with the hard-wall boundary condition, which overlooks the essential quantum mechanical effects like electron spill-out and quantum tunneling39,40,41,42,46. While some efforts have been made to include in the hydrodynamic description space-dependent electron density profiles, such as linear and quadratic shape functions25,26,29,30, the TFHT remains inadequate to deal with the realistic density profile, as it yields spurious modes originating from the exponentially decaying tail of the density46,47. The significance of an accurate description of the charge density profile in SHG has already been pointed out in the low-frequency regime both using a full quantum mechanical treatment32,33,34,35,36 as well as an orbital-free approach31.

Recently, a variety of methods based on effective local parameters have been proposed to deal with electron spill-out effects. These approaches are generally developed for the linear response regime48,49,50 and their generalization to the nonlinear dynamics is not trivial. The TFHT, on the other hand, is intrinsically nonlinear and it has long been used to describe the nonlinear response of metals43,44,45. In order to overcome its limitation, one should consider ∇n-dependent corrections (n being the electron density) to the free-electron gas kinetic energy, namely adding the von Weizsäcker term to the TF kinetic energy47,51,52. Such approach is usually referred to as quantum hydrodynamic theory (QHT). Over the past few years, QHT has been applied to study linear optical properties of individual nanoparticles as well as of extremely coupled nanostructures (in the tunneling regime) and a very good agreement with the time-dependent density functional theory calculations has been reported47,52,53,54. The QHT can also incorporate classical nonlinear effects, such as Lorentz and convective contributions, in multiscale plasmonic systems with deep sub-wavelength features, where the quantum effects cannot be neglected.

In the present work, we expand the QHT equations beyond the linear approximation to probe second-order nonlinearities in plasmonic structures of arbitrary shapes and sizes. Our method can efficiently handle realistic profiles of the charge density and indeed predicts three contributions in the SHG process that emerge from the spatial dependence of the charge density. We examine the spectral dependence of SHG from Na and Ag thick films and show that the SHG efficiency spectra can exhibit a strong resonance due to the electron spill-out at the metal surface. Our results show that quantum electron spill-out could be used to boost nonlinear efficiencies up to four orders of magnitude. If successfully proven experimentally, our work could enable efficient nonlinear optics at the nanoscale.

Results

In QHT the pressure term takes the more general form of \(\nabla \frac{\delta G[n]}{\delta n}\) where G[n] is the energy functional of the free-electron gas:

$$G[n]={T}_{{\rm{TF}}}[n]+\lambda {T}_{{\rm{W}}}[n,\nabla n]+{E}_{{\rm{xc}}}[n],$$
(1)

where TTF[n] is the TF kinetic energy functional, TW[n, ∇n] is the gradient-dependent correction term (von Weizsäcker term), and Exc[n] is the exchange-correlation energy functional in the local density approximation. The parameter λ, weighting the von Weizsäcker functional, is of extreme importance as it defines the decay of the electron density and it is usually taken as 1/9 ≤ λ ≤ 1 (ref. 47). The nonlinear QHT expansion up to second-order terms introduces three groups of nonlinearities in addition to the classical nonlinear hydrodynamic terms55,56 (see the “Methods” section for more details): (i) a convective term, proportional to ∇n0/n0, arising from the spatial dependence of the electron equilibrium density n0(r); (ii) a quantum pressure term given by the product of the linear induced density n1 and the gradient of the linear quantum pressure, i.e., \({n}_{1}\nabla {\left(\frac{\delta G[n]}{\delta n}\right)}_{1}\); (iii) a second-order quantum pressure term, \(\nabla {\left(\frac{\delta G[n]}{\delta n}\right)}_{2}^{{\rm{NL}}}\), that emerges from the nonlinear dependence of the kinetic functional from the electron density n. The last two terms generalize the nonlinear quantum pressure term in refs. 55 and 56 and introduce a set of contributions that are inversely proportional to some powers of n0. Note that because n0 → 0 near the metal surface, these contributions, as well as the convective contributions, are expected to be important at the surface. Further details about the nonlinear QHT are reported in the “Methods” section and Supplementary Note 1.

We employ the nonlinear QHT to investigate the linear and second-order nonlinear optical properties of Na and Ag films, illuminated by a fundamental field with energy E = ℏω1 impinging at an angle θi with respect to the normal unit vector \(\hat{n}\) at the metal surface, as depicted in Fig. 1. We assume that the electric field is transverse magnetic (TM) polarized, i.e., laying in the incidence plane, as a transverse electric (TE) polarized wave would not induce any accumulation of charge at the metal interface. To make the system more realistic, we place the metal film on a glass substrate, although its impact on SHG is negligible for sufficiently thick films.

Fig. 1: Schematic representation of the numerical setup for SHG.
figure 1

A thick metallic film placed on a glass substrate and illuminated by a transverse magnetic polarized plane wave making an angle θi with the unit normal vector \(\hat{n}\). The fundamental field (FF) and the second-harmonic field (SHF) correspond to the fields of the pump and the generated signal, respectively. A zoom-in at the metal–air interface depicts spatial variation of the static and induced charge densities at the fundamental and second harmonics.

Full size image

The SHG efficiency at ω2 = 2ω1 can be defined as

$${\mathrm{{SHG}}} \; {\mathrm{{efficiency}}}\,=\frac{{P}_{{\rm{SH}}}({\omega }_{2})}{{P}_{{\rm{FF}}}({\omega }_{1})},$$
(2)

where PSH(ω2) is second-harmonic reflected power from the metal interface and PFF(ω1) represents the input power of the exciting (fundamental) field. In the following, we consider the pump intensity of I0 = 80 MW/cm2.

A simple Drude-like metal

Let us first consider a simple Drude-like metal system, a Na film placed in free space. Optical response of Na is easier to understand as the effect of interband transitions can be neglected at visible frequencies. We numerically analyze the spectral dependence of a 400-nm-thick Na film excited by a TM-polarized plane wave. The self-consistent equilibrium charge density, calculated using Eq. (14) for different values of the spill-out parameter λ (see the “Methods” section), is shown in Fig. 2a. λ =  1/9 is associated with a smaller spill-out and a faster decay of the electron density whereas λ = 1/4 represents a larger spill-out and a slower decay, as can be seen in the inset of Fig. 2a. The equilibrium density in the case of TF approximation is also plotted. In this case, electrons are not permitted to spill from the metal surface into the free space (hard-wall boundary condition). Since the impact of interband transitions in Na can be neglected, we set the core dielectric constant, ε∞ = 1. In Fig. 2b, we plot the linear reflection spectra. As expected, because of the absence of any structuring on the Na film surface, the TFHT does not present any particular feature. On the other hand, the spatial dependence of the charge density in the QHT produces a dip in the reflection spectra that depends on the equilibrium density decay rate at the metal–air interface. As the electron spill-out increases (for larger values of λ), the resonance dip moves to lower energies. This means that, in principle, by controlling the amount of electron spill-out from the metal surface, one could tune the resonance to the frequency of interest. These resonances are due to excitation of multipolar surface plasmons and are a consequence of electron spill-out exclusively. In the TF approximation (as well as in the standard local case not reported here), in fact, no resonant behavior can be observed.

Fig. 2: Linear and nonlinear response of a 400-nm-thick Na film under transverse magnetic field excitation.
figure 2

a Equilibrium electron density calculated self-consistently for different values of spill-out parameter λ. For Thomas–Fermi hydrodynamic theory (TFHT) a constant density inside the metal with no spatial dependence is assumed. The inset shows the decay of the electron density for each λ and the solid gray line indicates the position of the metal interface. b Linear reflectance and c SHG efficiency as a function of exciting (fundamental) energy of the incident field at θi = 75∘ for different types of charge density profiles. We consider the peak intensity of the incident field I0 = 80 MW/cm2 and the electron damping rate ℏγ = 0.066 eV.

Full size image

The existence of this multipolar collective mode was predicted several decades ago for sufficiently diffused electron density profiles57 and was also observed directly in the electron energy loss spectroscopy (EELS) experiments58. This mode is expected to appear between the ordinary (Ritchie’s) surface plasmon and bulk plasmon modes. Experiments performed on a variety of simple metals for the investigations of multipole plasmons accord well with the theoretical predictions59,60,61,62. However, contrarily to the free-electron metals, there has been a controversy on the existence of the multipole mode in noble metals (i.e., Ag), which show complicated behavior due to filled 4d bands. EELS experiments on Ag does not show any clear evidence of the multipole plasmons. Nevertheless, experiments performed on Ag using high-resolution energy loss spectroscopy low-energy electron diffraction revealed a signature between the surface and bulk plasmon modes and it was attributed to the excitation of multipole plasmons63. Later on, Liebsch64 argued that the interpretation of the multipole mode in ref. 63 might not be correct and the multipole mode in Ag, instead of appearing before, should exist after the bulk plasma frequency. The existence of Ag multipole plasmon was also questioned due to the fact that the surface and bulk plasmons lie in a close proximity. In quest of the multipolar plasmon in noble metals, experiments were performed on Ag using a more suitable and sophisticated experimental technique, namely the angle- and energy-resolved photoyield, which confirmed the existing of the multiple plasmon mode both in the Ag bulk single-crystal surfaces and its adlayers65,66. The mode predicted by Liebsch64, however, was not observed in these experiments. On the other hand, in refs. 67,68, the authors performed experiments on thin Ag films grown on Ni and showed that they found a mode in agreement with the Liebsch’s prediction; however, this mode can only be excited under certain experimental conditions (i.e., enhanced surface sensitivity). The reader may refer to ref. 69 and the references therein for a detailed overview of the theoretical and experimental studies of higher-order multipolar surface plasmons.

These spill-out induced resonances, as shown in Fig. 2b, can be exploited to enhance the SHG efficiency by several order of magnitude (see Fig. 2c). The nonlinear QHT predicts, in fact, a very large enhancement of the SHG efficiency when the generated frequency matches the multipolar resonance in the linear spectrum. This kind of peak in the SHG intensity was also predicted in the nonlinear analysis of potassium surface within the TFHT30. Although the authors in this work used unrealistic equilibrium density, nonetheless, they pointed out that the position and strength of this peak strongly depend on the shape of the charge density profile. Our QHT calculations show, as depicted in Fig. 2c, that SHG enhancement can exceed four orders of magnitude than the conversion efficiency away from resonance.

Fitting experimental data of SHG from Ag films

Simple Drude-like metals, such as Na, are convenient from a theoretical point of view since they present an ideal optical behavior. In practice, however, they are very hard to work with because of their high reactivity. It would be interesting to analyze the possibility to observe a similar phenomenon in commonly used metals in plasmonics, such as noble metals. Noble metals’ optical response is complicated by the contribution of the interband transitions, which can act on many levels. They modify the linear susceptibility, contribute to the nonlinear response70, and affect both the equilibrium charge density and the induced charge distribution at the metal surface51,64,71. In the following, we neglect the nonlinear contribution from interband transitions and account for all the other effects by considering a polarizable medium of ε∞ > 1 permeating the whole metal volume.

In order to evaluate the relevance of the nonlinear QHT model, we consider a smooth Ag film of 400 nm thickness. We account for core electrons’ response with a local permittivity ε∞ = 5.9 and compare the numerical results with the experimental data already available in the literature21. The value of core permittivity ε∞ = 5.9 is taken from ref. 72, which is obtained by fitting the classical Drude model with the experimentally measured data. The SHG efficiency as a function of angle of incidence for different values of λ, computed within the full self-consistent nonlinear QHT, is compared against the experimental data, as shown in Fig. 3a. It can be seen that for smaller λ (lower electron spill-out) the efficiency is largely underestimated. However, as the value of λ increases (eventually spill-out increases), the SHG efficiency steps up and for λ = 1/2.5 it coincides well with the peak of the data, although the angular dependence is not thoroughly reproduced. The experimental data show a smooth angular dependence while the spectra predicted by the nonlinear QHT exhibit a minimum in the SHG efficiency in the range of 40∘–50∘. Although different from the measured data, this trend, and in particular the presence of a minimum, are in agreement with the time-dependent density functional theory calculations33. In order to fit the experimental data, we introduce a phenomenological parameter α weighting the nonlinear part of the von Weizsäcker energy functional given in Eq. (1) such that

$$\lambda {T}_{{\rm{W}}}[n,\nabla n]=\lambda \left[{T}_{{\rm{W}}}^{{\rm{L}}}+\alpha {T}_{{\rm{W}}}^{{\rm{NL}}}\right],$$
(3)

where the superscripts L and NL indicate linear and nonlinear contributions, respectively. By properly choosing an optimal value of α, the experimental data can be nicely approximated qualitatively as well as quantitatively for all angles of incidence. Figure 3b reports a comparison of the SHG efficiency between the experimental data and the nonlinear QHT when λ = 1/4 and α = 0.24(=1/4.15). It is important to remark that since the full self-consistent QHT does not result in a good agreement with the SHG measurements due to nature of the approximate functionals involved. The introduction of the parameter α, however, tunes the nonlinear contributions due to the von Weizsäcker terms and allows to fit the measured data. As a result, there exists more than one combination of λ and α that reproduces experimental data. The experimental results can also be approximated for other values of λ provided that a suitable value of the parameter α is chosen (see Supplementary Note 2). We would also like to point out that the choice of the parameter λ and α is not exact and it is a consequence of the non-universality of the kinetic energy functionals, the Thomas–Fermi-von Weizsäcker functionals, used to describe the free-electron internal energy. Ciracì and Della Sala47 have shown that the spill-out parameter can be fixed by matching the DFT equilibrium density asymptotic decay. The price to pay, however, is giving up self-consistency of the QHT model.

Fig. 3: SHG efficiency as a function of angle of incidence from a 400-nm-thick Ag film.
figure 3

The film is illuminated at E = 1.1653 eV with a peak intensity I0 = 80 MW/cm2 and considering the electron damping rate ℏγ = 0.030 eV: a for different values of the spill-out parameter λ and assuming α = 1, the parameter weighting the nonlinear part of the von Weizsäcker energy functional; b when λ = 1/4 and α = 0.24. The results are compared against the experimental data taken from ref. 21.

Full size image

In ref. 21, the authors performed experiments on Ag film for three different pumping-detection polarization conditions: (i) TM-incident and TM-detected, (ii) TE-incident and TM-detected, and (iii) 45∘ TE/TM-incident and TE-detected. The results shown in Fig. 3 correspond to the first case, i.e., when the incident and detected signals are TM-polarized. Note that in this condition, most of the surface contributions are activated. The SHG efficiency spectra for the other polarizations are presented in the Supplementary Note 3.

SHG resonances in Ag film with a dielectric coating

We have seen in the previous section that the Na films exhibit a strong enhancement in the SHG spectra due to the presence of resonances induced by the electron density decays at the metal surface. Contrarily to the Na case, however, we could not see any resonance structures neither in the linear nor in the SHG spectra for Ag films. The lack of such resonances is due to the combination of two factors. On the one hand, a larger ion density in Ag (smaller Wigner-Seitz radius) compared to Na leads to a greater work function that keeps the electron tighter at the metal surface, producing a much smaller spill-out. As we have seen by varying λ, a faster decay of the electron density pushes the resonance to higher energies. On the other hand, the presence of a polarizable background in Ag (due to interband transitions) shifts the effective plasma frequency to lower energies. The result is that the Ag film becomes transparent before any multipolar resonance can be excited, completely suppressing the behavior observed for Na.

One possible way to circumvent this issue is to increase the electron spill-out from the Ag surface by reducing the work function at the metal interface. The work function predicted for the QHT can be expressed as46

$$W=\lambda \frac{\hslash {k}_{\mathrm{D}}^{2}}{2{m}_{e}},$$
(4)

where kD describes the decay of the charge density proportional to \(\sim {e}^{-{k}_{\mathrm{D}}x}\). It can be seen from the above equation that when kD is smaller, i.e., the decay is slower, which implies that the spill-out is higher and therefore the work function decreases. The potential barrier at the metal surfaces can be lowered by reducing the difference in the dielectric constant between the metal background and the surround region. This can be obtained by modifying the dielectric environment surrounding the metal. In order to do so, we coat the film with a sub-wavelength thin layer of a dielectric material with a dielectric constant εr, as depicted in Fig. 4a.

Fig. 4: Linear and nonlinear response of a Ag film coated with a sub-wavelength thin layer of a dielectric material.
figure 4

a Simulation setup for SHG from the Ag film placed on a glass substrate. b Linear reflectance and SHG efficiency from a 400-nm-thick Ag film coated with a 10-nm-thick dielectric layer with dielectric constant εr = 5.9 (solid lines). The film is excited with a transverse magnetic polarized plane wave impinging at an angle θi = 75∘ by considering the spill-out parameter λ = 1/4 and the weight factor α = 0.24. The results for the Ag film without any coating are also shown for comparison (broken curves). The vertical dotted gray line represents the bulk plasma frequency. c Self-consistent equilibrium density n0 and d induced charge density n2 at E = 1.28 eV for coated and uncoated films, plotted along the normal direction to the metal interface. The vertical solid gray line marks the metal interface.

Full size image

The impact of the dielectric coating material on the equilibrium charge density is shown in the Supplementary Fig. 3 (see Supplementary Note 4). The linear reflectance and SHG efficiency spectra for a 400 nm Ag film, coated with a 10-nm-thick layer of a dielectric material with εr = ε∞ and placed on a glass substrate with refractive index n = 1.45, are plotted in Fig. 4b. Note that the figure also shows, for the sake of comparison, the spectra for the bare Ag film without any coating (as in Fig. 1). We observe that in the case of dielectric coating of Ag film, the linear response shows a dip in the reflectance which will result into a very strong resonance in the SHG efficiency, however, the Ag film without a coating does not show any resonance structure. This result is, in fact, very interesting as the spill-out induced resonances, which are suppressed in the commonly used setups in plasmonics, could be excited through a dielectric coating and, therefore, can be exploited for enhancing SHG.

The reason why the dielectric-coated Ag surface supports such resonant feature can be better understood by looking at the charge densities of the coated and uncoated Ag films. Figure 4c presents a comparison between the equilibrium charge density of a clean and dielectric-coated Ag film. It can be observed that that due to presence of the dielectric coating, the electron spill-out is more pronounced as compared to the uncoated slab. The dielectric material actually provides a higher electron density of states channel, allowing for the electrons to move further away from the metal surface, and as a consequence of this higher spill-out, the resonances can be observed in the linear and nonlinear spectra. Figure 4c presents the electron density both for the coated and uncoated films at λ = 1/4. Thus, the dielectric coating significantly assists the electron to spill-out of the metal surface and consequently producing a strong enhancement in the SHG efficiency. The induced charge density at the SHG, plotted at E  = 1.28 eV, for the coated and uncoated case is given in Fig. 4d. In both cases, the centroid of the charge density lies outside of the positive background edge. Although the spatial variation is quite similar, the induced charge density for the coated film is comparatively larger in width and magnitude, and is shifted farther away from the metal surface.

It is interesting to note that a similar behavior can be observed for different amounts of the electron spill-out and it is not the product of a specific choice of parameters. In Fig. 5a and b, we report linear and SHG spectra, respectively, for the dielectric-coated Ag slab for different values of λ and compare results against the TFHT (with no spill-out). For each λ, we observe a large resonance in the SHG spectrum which is related to the dip in the corresponding linear reflectance. Neglecting the electron spill-out from the metal surface, as in the TFHT, prevents such features from appearing neither in SHG nor in linear spectra.

The enhancement in the SHG efficiency is therefore strongly depends on the surface conditions. About three decades ago, monochromatic-light experiments were performed on dielectric (C60−  and SiOx− ) coated Ag island films to study the SHG and the sensitivity of the generated harmonic to the surface properties was emphasized. An order of magnitude enhancement in the SH response was observed in the coated Ag island films as compared to the uncoated ones73.

The SHG spectra of Fig. 5b display a noticeable minimum just before the peak resonance for each λ, which is not seen in the situation shown in Fig. 4b. To examine this feature, we consider λ = 1/4 and compute the SHG efficiency as a function of the angle of incidence and the fundamental energy for different values of α. We find that when α = 1, there exists a minimum at almost all angles of incidence (Fig. 5c) whereas for α = 0.24, there is no such minimum before the resonance; nevertheless, in this scenario a minimum also appears but after the resonance (Fig. 5d). In the latter case, the minimum is not present for all angles of incidence and since the spectrum shown in Fig. 4b is plotted at θi = 75∘, therefore, no such minimum is detected there. This minimum in the SHG efficiency can be associated to the interplay between the different nonlinear contributions that involve the first- and second-order terms of the von Weizsäcker energy functional, respectively. While varying λ has a simultaneous effect on both terms, varying α modifies their respective weights, modifying the condition for which the two contributions cancel out.

Fig. 5: Optical reflection from the dielectric-coated Ag film considering different amounts of electron spill-out from the metal surface.
figure 5

a Linear reflectance and b SHG efficiency as a function of exciting energy (in eV) of the incident field for different values of the spill-out parameter λ and considering the weight factor α = 1. The film is coated with a 10-nm-thick dielectric layer with a dielectric constant εr = 5.9 and excited with a transverse magnetic polarized light incident at an angle θi = 75∘. The results are compared against the Thomas–Fermi hydrodynamic theory (TFHT) spectra represented by the dash-dotted black curves. The vertical gray lines indicate the bulk plasma frequency ωp at the fundamental and generated harmonic. SHG efficiency of the dielectric-coated Ag slab, for λ = 1/4 considering c α = 1 and d α = 0.24, plotted as a function of exciting energy and angle of incidence.

Full size image

So far, we have fixed the refractive index of the coating material such that it corresponds to ε∞ in order to reproduce a similar condition to Na where ε∞ = 1. It would be interesting at this point, studying the influence of the dielectric constant of the coating material on the SHG efficiency. In Fig. 6, we report numerical results for linear response and SHG of the Ag film considering λ = 1/4, α = 0.24, and various values of the dielectric constants εr of the coating material. The spectra for silver–air interface (εr = 1, i.e., without any coating) are also shown for completeness. Contrarily to the Ag–air interface, when εr ≠ 1, a dip in the linear reflectance can be seen for each εr, although it becomes less intense for the smaller values of εr and completely disappears for εr = 1, as shown in Fig. 6a. This structure in the linear spectrum yields a large efficiency of the generated signal, as depicted in Fig. 6b. It is interesting to note that as the value of the dielectric constant of the coating layer increases, the resonance shifts towards the lower energies. This result, in fact, could be very important in many optical applications as the resonance of the SHG process can be tuned by coating the Ag film with a suitable dielectric material. A minimum in the SHG spectra for a lower value of the dielectric constants (with the exception of εr = 1) is also visible in the SHG spectra, which appears after the resonance in the vicinity of bulk plasma frequency and corresponds to the similar situation as was discussed in Fig. 5b. However, this minimum vanishes for the higher values of the dielectric constant.

Fig. 6: Impact of the dielectric constant of the coating material on the reflection properties of Ag film.
figure 6

a Linear reflectance and b SHG efficiency as a function of exciting energy plotted at an incident angle θi = 75∘ for various values of the dielectric constant εr considered for the coating of Ag film. The result for the clean Ag film without any coating (εr = 1) is also reported. The linear and nonlinear spectra are arbitrarily shifted in the vertical direction for the sake of a better visualization.

Full size image

Discussion

We have presented a theoretical microscopic model to probe second-order nonlinearities at metal surfaces. Our approach is based on the QHT, which can efficiently handle the realistic profiles of the ground-state electron density. The nonlinear QHT predicts a very strong resonance in the SHG spectra that can boost the conversion efficiency by several orders of magnitude. This resonance shows a strong dependence on the spill-out of the electron density from the metal surface and is related to the excitation of multipolar surface plasmons in the linear response. Interestingly, these spill-out induced resonances can be tuned by controlling decay of the electron density from the metal surface with the aid of a dielectric coating.

Our results suggest that an accurate description of the spatial variation of charge density at a metal surface is crucial in characterizing an accurate nonlinear optical response. Although we have limited our study to films, the spill-out induced enhancement of SHG efficiency could be combined with plasmon-based field enhancement techniques to further increase nonlinear conversion rates. Moreover, because the nonlinear process is confined to the metal surface, such processes could be embedded in a hybrid metal-dielectric waveguide configuration that allows the build up of the nonlinear signals. We believe that the analysis presented in this article will be a source of motivation for further analytical and experimental research.

Methods

Nonlinear QHT

Materials with a very high density of free charges can be treated as charged fluids (in a positively charged background) and their optical properties can be characterized under the hydrodynamic description. The full nonlinear hydrodynamic equation of motion for an electronic system under the influence of electromagnetic fields E and B can be expressed as47

$${m}_{e}\left(\frac{\partial }{\partial t}+{\bf{v}}\cdot \nabla +\gamma \right){\bf{v}}=-e\left({\bf{E}}+{\bf{v}}\times {\bf{B}}\right)-\nabla \frac{\delta G[n]}{\delta n},$$
(5)

where me is the electron mass, v is the electron velocity field, γ is the damping rate, and e is the absolute value of the electron charge; n(r) represents the electron density and G[n] = TTF[n] + λTW[n, ∇n] + Exc[n] is the energy functional incorporating the TF kinetic energy, the gradient-dependent correction term and the exchange-correlation energy functional respectively. In specifying the energy functionals in Eq. (1), we do not include the effect of nonlocal broadening (see ref. 52 for more details) to avoid further numerical complexity. Considering the electron current J = −nev, Eq. (5) can be expressed as

$$\frac{\partial {\bf{J}}}{\partial t}+\gamma {\bf{J}}=\frac{{e}^{2}n}{{m}_{e}}{\bf{E}} \, -\frac{e}{{m}_{e}}{\bf{J}}\times {\bf{B}}+\frac{en}{{m}_{e}}\nabla \frac{\delta G[n]}{\delta n}\\ \,+\frac{1}{e}\left(\frac{{\bf{J}}}{n}\nabla \cdot {\bf{J}}-{\bf{J}}\cdot \nabla \frac{{\bf{J}}}{n}\right).$$
(6)

A similar equation can also be directly derived from the single-particle Kohn–Sham equation52.

We proceed by writing the fields as the sum of few harmonics:

$$n({\bf{r}},t)={n}_{0}({\bf{r}})+{n}_{1}({\bf{r}}){e}^{-i{\omega }_{1}t}+{n}_{2}({\bf{r}}){e}^{-i{\omega }_{2}t}+{\rm{c.c.}},$$
(7)
$${\bf{E}}({\bf{r}},t)={{\bf{E}}}_{0}({\bf{r}})+{{\bf{E}}}_{1}({\bf{r}}){e}^{-i{\omega }_{1}t}+{{\bf{E}}}_{2}({\bf{r}}){e}^{-i{\omega }_{2}t}+{\rm{c.c.}},$$
(8)
$${\bf{J}}({\bf{r}},t)={{\bf{J}}}_{1}({\bf{r}}){e}^{-i{\omega }_{1}t}+{{\bf{J}}}_{2}({\bf{r}}){e}^{-i{\omega }_{2}t}+{\rm{c.c.}},$$
(9)
$${\bf{B}}({\bf{r}},t)={{\bf{B}}}_{1}({\bf{r}}){e}^{-i{\omega }_{1}t}+{{\bf{B}}}_{2}({\bf{r}}){e}^{-i{\omega }_{2}t}+{\rm{c.c.}},$$
(10)

where ω2 = 2ω1 and E0 and n0 represent the steady-state electric field and ground-state electron density. The subscripts 1 and 2 indicate the field quantities at the fundamental and the second harmonic, respectively. Assuming the undepleted pump approximation, that is, the fundamental field is not affected by the conversion process, and recalling that \({\bf{J}}=\dot{{\bf{P}}}\), in view of Eqs. (7)–(10), Eq. (6) results into following set of equations:

$$-\frac{e{n}_{0}}{{m}_{e}}\nabla {\left(\frac{\delta G[n]}{\delta n}\right)}_{1}-\left({\omega }_{1}^{2}+i\gamma {\omega }_{1}\right){{\bf{P}}}_{1}={\varepsilon }_{0}{\omega }_{p}^{2}{{\bf{E}}}_{1},\,\,\qquad$$
(11)
$$-\frac{e{n}_{0}}{{m}_{e}}\nabla {\left(\frac{\delta G[n]}{\delta n}\right)}_{2}-\left({\omega }_{2}^{2}+i\gamma {\omega }_{2}\right){{\bf{P}}}_{2}={\varepsilon }_{0}{\omega }_{p}^{2}{{\bf{E}}}_{2}+{{\boldsymbol{S}}}_{{\rm{NL}}},$$
(12)

where \({\omega }_{\mathrm{p}}({\bf{r}})=\sqrt{{e}^{2}{n}_{0}({\bf{r}})/({m}_{e}{\varepsilon }_{0})}\) is the space-dependent plasma frequency and the nonlinear source term is

$${{\bf{S}}}_{{\rm{NL}}}= \, \frac{{e}^{2}{n}_{1}}{{m}_{e}}{{\bf{E}}}_{1}+i{\omega }_{1}\frac{{\mu }_{0}e}{{m}_{e}}{{\bf{P}}}_{1}\times {{\bf{H}}}_{1}\\ \, -\frac{{\omega }_{1}^{2}}{e{n}_{0}}\left({{\bf{P}}}_{1}\nabla \cdot {{\bf{P}}}_{1}+{{\bf{P}}}_{1}\cdot \nabla {{\bf{P}}}_{1}-\frac{{{\bf{P}}}_{1}\cdot \nabla {n}_{0}}{{n}_{0}}{{\bf{P}}}_{1}\right)\\ \, +\frac{e{n}_{1}}{{m}_{e}}\nabla {\left(\frac{\delta G[n]}{\delta n}\right)}_{1}+\frac{e{n}_{0}}{{m}_{e}}\nabla {\left(\frac{\delta G[n]}{\delta n}\right)}_{2}^{{\rm{NL}}}.$$
(13)

Mathematical expressions for the linear potentials \({\left(\frac{\delta G[n]}{\delta n}\right)}_{1,2}\) are reported in detail in ref. 47 whereas expressions for the nonlinear potentials are derived in the Supplementary Note 1. In addition to the classical nonlinear hydrodynamic terms55,56, i.e., Coulomb, Lorentz, and convective terms (note that there is one extra convective term, ∝ ∇n0/n0, arising from the spatial dependence of the electron density), two quantum terms emerge as discussed in the “Results” section.

Self-consistent equilibrium charge density

In implementing our nonlinear QHT, we follow the jellium approximation74, which assumes that the electrons in a metal are confined by a constant positive background charge \({n}^{+}={(\frac{4}{3}\pi {r}_{s}^{3})}^{-1}\), where rs indicates the Wigner-Seitz radius (rs = 4 for Na and rs = 3 for Ag). The space-dependent equilibrium charge density can be calculated by combining the 0th order QHT equation with the Gauss’s law, which leads to the following nonlinear static equation52:

$$\nabla {\varepsilon }_{\infty }({\bf{r}})\cdot \nabla {\left(\frac{\delta G}{\delta n}\right)}_{n = {n}_{0}}+\frac{{e}^{2}}{{\varepsilon }_{0}}\left({n}_{0}-{n}^{+}\right)=0,$$
(14)

where ε0 symbolizes permittivity of free space. ε∞ is the core dielectric constant presenting a local contribution to the metal permittivity.

We solve Eqs. (11)–(14) coupled to two frequency-dependent wave equations in a finite-element-method-based numerical simulator COMSOL Multiphysics75, which offers a quite flexible implementation of these expressions. The numerical implementation of these equations is given in Supplementary Note 5.