Enhancing second-harmonic generation with electron spill-out at metallic surfaces

Second-order nonlinear optical processes do not manifest in the bulk of centrosymmetric materials, but may occur in the angstroms-thick layer at surfaces. At such length-scales, quantum mechanical effects come into play which could be crucial for an accurate description of plasmonic systems. In this article, we develop a theoretical model based on the quantum hydrodynamic description to study free-electron nonlinear dynamics in plasmonic systems. Our model predicts strong resonances induced by the spill-out of electron density at the metal surface. We show that these resonances can boost second-harmonic generation efficiency up to four orders of magnitude and can be arbitrarily tuned by controlling the electron spill-out at the metal surface with the aid of thin dielectric layers. These results offer a possibility to artificially increase nonlinear susceptibilities by engineering optical properties at the quantum level.


INTRODUCTION
M odern 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 switching 1,2 . Furthermore, researchers are constantly extending the reach of nonlinear optical phenomena to a growing variety of practical applications, such as nonlinear microscopy 3 , ultrasensitive-optical shape characterization 4 , monitoring the processes in chemical synthesis of nanostructures 5 and nonlinear plasmonic sensing 6 .
The inherent weakness of optical nonlinearities however has strongly hindered the chip integration of alloptical 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 * muhammad.khalid@iit.it † cristian.ciraci@iit.it 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 susceptibilities 7 . Increasing the carrier density, however, causes a transition from size-quantization to the classical regime of plasmon oscillations, in the visible range of frequencies.
In centrosymmetric materials, such as noble metals, secondorder nonlinearities in the bulk region are intrinsically forbidden 3 , nevertheless, the symmetry can be broken at the material surface, giving rise to a surface SHG mechanism 37,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 spillout and quantum tunneling [39][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 functions 25,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 density 46,47 . The significance of an accurate description of the charge density profile in SHG has already been pointed out in the lowfrequency regime both using a full quantum mechanical treatment 32-36 as well as an orbital-free approach 31 .
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 regime [48][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 metals [43][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 energy 47,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 reported 47,52-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 4 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 ∇ δG [n] δn where G[n] is the energy functional of the freeelectron gas: where the spatial dependence of the electron equilibrium density n 0 (r); ii) a quantum pressure term given by the product of the linear induced density n 1 and the gradient of the linear quantum pressure, i.e., n 1 ∇ δG[n] iii) a second-order quantum pressure term, ∇ δG[n] δn NL 2 , 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 n 0 . Note that because n 0 → 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 =hω 1 impinging at an angle θ i with respect to the normal unit vectorn 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.
The SHG efficiency at ω 2 = 2ω 1 can be defined as: where P SH (ω 2 ) is second-harmonic reflected power from the metal interface and P FF (ω 1 ) represents the input power of the exciting (fundamental) field. In the following, we consider the pump intensity of I 0 = 80 MW/cm 2 .

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 selfconsistent equilibrium charge density, calculated using equation (14) for different values of the spill-out parameter λ (see the Methods section), is shown in Fig. 2a. λ = 1/9 is associated to 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. The existence of this multipolar collective mode was predicted several decades ago for sufficiently diffused electron density profiles 57 and was also observed directly in the electron energy loss spectroscopy (EELS) experiments 58 . 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 predictions [59][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 plasmons 63 . Later on, Liebsch 64 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 singlecrystal surfaces and its adlayers 65,66 . The mode predicted by Liebsch 64 , however, was not observed in these experiments. On the other hand, in Ref. 67 and 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 TFHT 30 . 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 4 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 response 70 and affect both the equilibrium charge density and the induced charge distribution at the metal surface 51,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 literature 21 . 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 shows 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 calculations 33 . 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 equation (1) such that: 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 Sala 47 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. In Ref. 21, the authors performed experiments on Ag film for three different pumping-detection polarization conditions: (i) TM-incident and TM-detected, (ii) TEincident 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 as 46 : where k D describes the decay of the charge density proportional to ∼ e −k D x . It can be seen from the above equation that when k D 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. 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 magni- tude, and is shifted farther away from the metal surface.
It is interesting to note that a similar behavior can be observed for different amount of the electron spill-out and it is not the product of a specific choice of parameters. In Figs. 5a and 5b, 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 (C 60 − and SiO x −) 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 ones 73 .
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 firstorder 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.
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.

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 effi-ciency 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.

Nonlinear quantum hydrodynamic theory
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 as 47 : where m e 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] = T TF [n] + λT W [n, ∇n] + E xc [n] is the energy functional incorporating the TF kinetic energy, the gradient-dependent correction term and the exchangecorrelation energy functional respectively. In specifying the energy functionals in equation (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, equation (5) can be expressed as: A similar equation can also be directly derived from the single-particle Kohn-Sham equation 52 . We proceed by writing the fields as the sum of few harmonics: n(r, t) = n 0 (r) + n 1 (r)e −iω1t + n 2 (r)e −iω2t + c.c., (7) where ω 2 = 2ω 1 and E 0 and n 0 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 J =Ṗ, in view of equations (7)-(10), equation (6) results into following set of equations: where ω p (r) = e 2 n 0 (r)/(m e ε 0 ) is the space-dependent plasma frequency and the nonlinear source term is:  55,56 , i.e., Coulomb, Lorentz and convective terms (note that there is one extra convective term, ∝ ∇n 0 /n 0 , 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 approximation 74 , which assumes that the electrons in a metal are confined by a constant positive background charge n + = ( 4 3 πr 3 s ) −1 , where r s indicates the Wigner-Seitz radius (r s = 4 for Na and r s = 3 for Ag). The space-dependent equilibrium charge density can be calculated by combining the zero-th order QHT equation with the Gauss's law, which leads to the following nonlinear static equation 52 : where ε 0 symbolizes permittivity of free space. ε ∞ is the core dielectric constant presenting a local contribution to the metal permittivity. We solve equations (11)- (14) coupled to two frequencydependent wave equations in a finite-element-method based numerical simulator comsol Multiphysics 75 , which offers a quite flexible implementation of these expressions. The numerical implementation of these equations is given in the Supplementary Note 5.

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

CODE AVAILABILITY
The codes that support the findings of this study are available from the corresponding author upon reasonable request.
Focusing only on the second-order terms and neglecting the higher harmonics (ω > ω 2 ), we are left with: Expanding n 0 + n j e −iωj t −2 in the above equation and doing some algebra results in the following nonlinear contribution: Similarly, the nonlinear terms for exchange-correlation potential can be obtained by using Supplementary Equation (S8) into Supplementary Equation (S4), resulting in the following expression: where, The nonlinear quantum hydrodynamic theory (QHT) can also approximate the experimental data 3 for Ag film for a given λ, representing the electron spill-out from the metal surface, provided that the value of the parameter α, weighting the nonlinear part of the von-Weizsäcker energy functional (see equation (3) in the main text), is properly chosen. As we discussed in the main text, since the QHT is based on the approximate energy functionals, the full self-consistent approach does not reproduce well the experimental data. This is because the QHT shows a strong dependence on the equilibrium charge density which is affected by the choice of λ and this in turn affects the SHG efficiency. The parameter α, on the other hand, changes the relative weight of the nonlinear contributions emerging from the von Weizsäcker terms. Therefore, there is more than one combination of λ and α that enables us to fit the experimental data. Supplementary Figure 1  E/TM-incident and TE-detected. Case (i) has already been reported in the main text (Fig. 3b). Here we attempt to fit the experimental data for the other two cases within the framework of QHT. A comparison of the SHG efficiency between the measured and simulated data is shown in Supplementary Figure 2. The choice of the parameters is the same as used in Fig. 3b of the main text. In case (iii), we used an input intensity I 0 = 160 MW/cm 2 . It is interesting to note that even under these polarization conditions, the QHT approximates the experimental data quite well.  Figure 3 shows the impact of the dielectric constant of the coating material on the equilibrium charge density. It can be observed how the electron density is more loosely localized (slower decay) at the metal surface when increasing values of the dielectric constant of the coating material are considered. This corresponds to lowering the work function W (see Eq. (4) in the main text).

Supplementary Note 5: Weak formulation
In finite-element-method based approaches, the differential equations are usually expressed in their weak form to relax the continuity requirements for the test functions. We numerically implement equations (11)-(14) of the main text and two frequency-dependent wave equations in finite-element-method based solver comsol Multiphysics 4 . In the following, we derive the weak formulation of these equations. Rewriting equations (11) and (12)  − ω 2 j + iγω j P j = ε 0 ω 2 p E j + S j , where, j = 1, 2 (S14) where S 1 = 0 and S 2 = S NL , given by equation (13) of the main text, as: S NL = e 2 n 1 m e E 1 + iω 1 µ 0 e m e P 1 × H 1 − ω 2 1 en 0 P 1 ∇ · P 1 + P 1 · ∇P 1 − P 1 · ∇n 0 n 0 P 1 + + . (S15) The weak formulation of Supplementary Equation (S14) can be obtained by multiplying it with a test functionP j and integrating over the simulation domain.
e m e ∇ δG[n] δn j + 1 n 0 ω 2 j + iγω j P j + ε 0 ω 2 p E j + 1 n 0 S j ·P j dV = 0, (S16) − e m e δG[n] δn j ∇ ·P j + 1 n 0 ω 2 j + iγω j P j + ε 0 ω 2 p E j ·P j + 1 n 0 S j ·P j dV = 0. (S17) Since S 1 = 0, the integral of the nonlinear source (last) term in the Supplementary Equation (S17) can be expressed as: (S18) The last two terms appearing in the above equation emerge from the spatial dependence of the charge density. Thus, its weak form can be finally written as: (S19) The advantage of the weak formulation is that there is no need to compute the gradients of the energy functionals, and the derivatives are distributed to the test functions. The energy functionals also contain second-order derivatives, and we introduced working variables, such that, F j = ∇n j where n j = 1 e ∇ · P j which implies that ∇ · F j = ∇ 2 n j . Thus, the complete set of nonlinear QHT equations, along with the two wave equations and a static QHT equation for computing equilibrium density, in their weak formulation reads, − e m e δG[n] δn j ∇ ·P j + 1 n 0 ω 2 j + iγω j P j + ε 0 ω 2 p E j ·P j + 1 n 0 S j ·P j dV = 0, (S21) · ∇ξ + e 2 ε 0 ξ 2 − n + ξ dV = 0, where j = 1, 2 and S 1 = 0. In Supplementary Equation (S23), the transformed variable ξ = √ n 0 was introduced as in this way the solution converges more quickly.