Nonlinear plasmonic response in atomically thin metal films

2.1 Nonlinear response

We extend the linear QM model outlined above to describe nonlinear processes within the self-consistent field approximation, following a similar procedure as previously applied to nanostructured graphene [12, 56] (see details in Methods and SI). Aside from the linear contribution, the field induced at the surface by a nearby emitter has contributions E ⁽ⁿ⁾(sω) ≡ E ^(n,s) from nonlinear optical processes characterized by their perturbation order n = 1, 2, … and harmonic index s (with |s| ≤ n). Up to third order, these include linear response (n = s = 1), second- (n = s = 2) and third-harmonic (n = s = 3) generation, and an intensity-dependent correction to the response at the fundamental frequency from the optical Kerr effect (n = 3, s = 1). These processes are illustrated schematically in Figure 2(a) for evanescent fields at the surface of an ultrathin crystalline Ag(111) film with thickness d = Na _s , where N denotes the number of stacked atomic planes and a _s = 0.236 nm the Ag(111) interlayer spacing. We consider N = 4 and also plot the atomic layer potential V(z) and the corresponding background charge distribution ρ ⁽⁰⁾, along with the external field and the induced charge densities ρ ^(n,s) for incident Q and ω within the SPP dispersion curve. Following the QM prescription in Methods to calculate the nonlinear induced fields, we find the results presented in Figure 2(b)–(g) as functions of optical in-plane wave vector Q and frequency ω for films consisting of N = 12 (Figure 2(b)–(d)) and N = 30 (Figure 2(e)–(g)) Ag(111) atomic layers. In all cases, prominent features in the generated nonlinear near field follow the plasmon dispersion curve ω _sp(Q), with a more subtle peak emerging in second- and third-harmonic generation at ω _sp(Q/2)/2 and ω _sp(Q/3)/3, respectively (i.e., when the generated frequency matches the plasmon, indicated by red arrows). We note that vertical transitions between occupied and unoccupied bands result in weaker dispersionless features, which may be artificially enhanced by the assumption of parabolic bands in the QM model employed here [54]; signatures of such features with similar characteristics have been reported in experiment [57, 58], warranting further investigation with more rigorous treatments of the electronic band structure in future studies.

To quantify the strength of the plasmon-driven nonlinearity in the near field, we introduce a figure of merit (FoM) for evanescent fields,

obtained from the ratio of the induced potential ϕ ^(n,s)(Q, z ₀) = ∫ d z v ^(s)(Q, z ₀, z)ρ ^(n,s)(Q, z) (computed from the induced nonlinear charge density ρ ^(n,s) and evaluated at the source position z = z ₀) to the external potential ∝ e − n Q ( z − z 0 ) , and recognizing that the Coulomb interaction v ^(s)(Q, z ₀, z) introduces an additional e − s Q z 0 factor at the frequency and wave vector generated by a nonlinear process of order n and harmonic s. In this manner, the surface nonlinearity for various thin films can be directly compared at specific frequencies and wave vector components, independently of the excitation origin z ₀.

In Figure 3 we apply the FoM to quantify second- and third-order nonlinear optical processes in Ag(111) films. Figure 3(a)–(c) shows SHG, THG, and the optical Kerr effect at a fixed in-plane wave vector Q = 0.1 nm⁻¹ for various film thicknesses, as predicted from the QM model (solid curves) and contrasted with CL (dashed curves) predictions based on a purely two-dimensional treatment of ultrathin metal films (see Methods). As expected, we find that the response is dominated by the plasmonic peak, while its echoes in harmonic generation, indicated by the red arrows in Figures 2 and 3, appear with much lower intensity. The CL model predicts the spectral positions of the dominant features and reasonably captures their amplitudes for SHG. These maxima are featured in Figure 3(e)–(g) as functions of the number of atomic layers N at selected in-plane wave vectors Q, and are confirmed to generally increase with decreasing film thickness, particularly for films comprised of less than N = 20 atomic planes.

Figure 2:

Nonlinear response of crystalline silver films. (a) Schematic illustration of the atomic layer potential V(z) (black curve), background electron density ρ ⁽⁰⁾ (orange curve), external potential with fixed frequency ω and in-plane wave vector Q, and induced charge density components ρ ^(n,s) corresponding to processes of order n and harmonic s in a film consisting of N = 4 Ag(111) atomic layers. We take Q = 0.1 nm⁻¹ and ℏω = 1.85 eV to lie on the SPP curve. (b)–(g) Normal component of the induced fields E z ( n , s ) associated with (b) and (e) SHG (n = s = 2), (b) and (f) third-harmonic generation (n = s = 3), and (d) and (g) the optical Kerr effect (n = 3, s = 1) as functions of in-plane wave vector and photon energy for films consisting of N Ag(111) atomic layers with (b)–(d) N = 12 and (e)–(g) N = 30.

Figure 3:

Nonlinear figure of merit. (a)–(c) Calculated figure of merit F ( n , s ) as a function of the input frequency with fixed in-plane wave vector Q = 0.1 nm⁻¹ for (a) SHG, (b) THG, and (c) the Kerr effect. We consider a range of film thicknesses as indicated by the color scale. Solid and dashed curves correspond to QM and CL simulations, whose peak maxima are tracked with solid and open symbols, respectively. Red arrows in (a,b) highlight the SHG and THG replica of the plasmonic peak, equivalent to those of Figure 2. (e)–(g) Variation of the maxima of F ( n , s ) as a function of film thickness for a set of in-plane wave vectors (see colored labels in (e)).

Inspecting Figure 3, we note that SHG undergoes a significant increase with decreasing film thickness, presumably due to the reliance of second-order optical nonlinearities on surface symmetry-breaking in an otherwise centrosymmetric medium, combined with the higher surface-to-volume ratio in thin films. In contrast, third-order processes dominated by the bulk response undergo a more modest increase in thinner films, and even dropping significantly in the Kerr nonlinearity for N ≲ 8. The differing nature of second- and third-order nonlinear processes in ultrathin films corroborates our finding that the CL method, based on a two-dimensional electron gas, cannot satisfactorily predict the amplitudes of THG and the Kerr nonlinearity. Additionally, the dimensionless quantity Qd characterizes the nonlocal response of the system and the spectral position of the plasmon, and consequently, the maxima in the nonlinear response for low in-plane momenta converge to their semi infinite values more slowly with thickness (e.g., comparing the blue and purple curves for Q = 0.1 nm⁻¹ and Q = 1.0 nm⁻¹, respectively). In general, excitation with small Q results in more dramatic changes in optical nonlinearity with respect to film thickness.

The CL model assuming a two-dimensional film is in good agreement with the QM model, but neglects finite-size effects, which become important for films comprised of N ≲ 20 atomic layers, as indicated by the additional features appearing in the spectra of Figure 3(a)–(c) in the QM model, originating in vertical single-electron transitions (see horizontal lines in Figure 2). The spectral positions of these features correlate with discrete levels in the electronic spectrum of ultrathin films, and thus vary dramatically with film thickness, leading to a complex dynamics for the prediction of the peaks, in agreement with experimental SHG measurements [38–41]. These findings underscore the importance of quantum finite-size effects in the thin film nonlinear response, which tend to cause increasing oscillatory behavior in the yield of higher-order processes, particularly those involving the generation of new frequencies, due to the involvement of additional vertical electronic transitions.

2.2 Manifestation of surface states in the nonlinear optical response

The features emerging at low energies ℏω < 0.5 eV in Figures 1(b) and 2(b)–(g), which increase slowly with in-plane wave vector, originate in Shockley surface states (SSs) characteristic of the (111) metallic crystallographic orientation [60]. Unlike the weaker dispersionless features appearing in the nonlinear response due to vertical electronic transitions, the SSs are nearly independent of film thickness, each behaving as a two-dimensional electron gas that hybridizes with the bulk three-dimensional electron gas to produce intrinsic acoustic plasmons [59]. In Figure 4, we explore the low-energy peaks associated with SSs in the nonlinear optical response of thick films comprised of N = 55 Ag(111) atomic layers, which we regard as a well-converged limit representing a semi-infinite film. Specifically, we present the normal second-harmonic induced field in Figure 4(a), which reveals a broad, low-energy acoustic plasmon feature that rapidly disappears with increasing in-plane wave vector Q, presumably due to the onset of Landau damping. The acoustic plasmon is accompanied by two narrower features that persist at larger Q, which we attribute to single- and two-photon transitions involving the SSs. For comparison, we plot the linear acoustic dispersion (SS-Lin, black solid line) predicted in a previous work [59] and covering the energy range under consideration. In Figure 4(b), we compare the FoM for the SHG nonlinear response at selected values of Q, for which the SSs are clearly resolved at large in-plane wave vectors. Qualitatively similar behavior is observed for third-order nonlinear optical processes, as shown in Figure 4(c) and (d), where we present the nonlinear near field and FoM for THG, respectively. We note that a third parabolic feature emerges in the dispersion diagram of THG, corresponding to three-photon single-electron excitations involving the SSs. The commented SHG and THG response features are parametrized in Table 1 for THG and THG.

Figure 4:

Shockley states in the nonlinear optical response. (a) Calculated dispersion diagram of SHG in the low Q and ω corner, revealing a nonlinear enhancement that involves Shockley states. The solid black line (SS-Lin) indicates the position of the SS acoustic plasmon demonstrated in the linear response and obtained from ref [59], whereas in SHG two resonances emanate at a different position than SS-Lin and are indicated by the dashed blue and red parabolic curves. (b) Low-frequency nonlinear FoM spectra for selected Q, indicated by the color-coded curves, in which surface-state features are clearly resolved. Panels (c) and (d) show the THG field and corresponding FoM in the same parameter space as panels (a) and (b), respectively. A third parabolic feature is indicated by the dashed yellow curve in (c). All curves are parametrized in Table 1.

4.1 Near-field excitation by a point dipole

With near-field experiments in mind, we consider the external field produced by an electric dipole p = − p z ̂ oriented along the z-direction (e.g., representing an SNOM tip) a distance z ₀ > 0 above the upper surface of a thin film occupying the R = (x, y) plane and extending up to a distance z = d. Assuming a monochromatic time dependence e^−iωt with frequency ω and working in the electrostatic limit, the external potential ϕ ext ( r , ω ) = ϵ d − 1 p ∂ z | r − z 0 z ̂ | − 1 is decomposed in wave vector Q components using the Weyl identity according to

where ϵ _d is the dielectric function of the homogeneous medium that surrounds the film (we set ϵ _d = 1 for vacuum in the calculations here presented). The external potential excites the metal thin film, inducing fields that act back on the dipole (position R = 0, z = z ₀), where the self induced field is given by

Here, we have separated the induced potential ϕ ^(1,1)(Q, ω) in Q-dependent components, which we compute following either classical (CL) or quantum-mechanical (QM) frameworks, as outlined below.

4.2 Electronic states of atomically-thin metal films

Following the prescription of ref [51], we consider a metal film of finite thickness d in the z direction, so that an electron propagating in the R = (x, y) plane with 2D momentum ℏ k _∥ is characterized by the single-particle wave function Ψ j , k ∥ ( r ) = A − 1 / 2 e i k ∥ ⋅ R φ j ( z ) , with A denoting the quantization area and φ _j (z) the out-of-plane wave function component obtained by solving the 1D Schrödinger equation H 0 φ j ( z ) = ℏ ε j ⊥ φ j ( z ) for the Hamiltonian H 0 = ( − ℏ 2 / 2 m e ) ∂ z 2 + V ( z ) . The associated energy eigenvalues ℏ ε j ⊥ of the 1D problem are supplemented in the full electron dispersion ℏ ε j , k ∥ = ℏ ε j ⊥ + ℏ 2 k ∥ 2 / 2 m j by a parabolic term describing free-electron motion with effective mass m _j in directions parallel to the film surface. The choice of 1D potential V(z) determines the nature of the metal under consideration. We elaborate below and in the SI on the procedure used to determine the effective masses m _j .

In the presence of the electrostatic potential ϕ, the electron dynamics is governed by the Liouville–von Neumann equation

supplemented by a term that accounts for relaxation of the density matrix ρ at a phenomenological rate γ to its t → −∞ equilibrium state ρ ( 0 ) ( r , r ′ ) = ∑ i f i Ψ i ( r ) Ψ i * ( r ′ ) , where i ≡ {j, k _∥} denotes the multiplexed state index. Following Fermi–Dirac statistics at zero temperature, the equilibrium state is determined by the Fermi energy E _F through filling factors f i = Θ E F − ℏ ε i , where Θ(x) denotes the step function. In practice, the Fermi energy is computed by populating states according to

where d is the film thickness (e.g., d = Na _s for a film of N atomic layers with interlayer separation a _s ), n _eff is the effective electron density, and M is the highest partially occupied band, determined by the condition ℏ ε M ⊥ < E F < ℏ ε M + 1 ⊥ .

To compute the linear and nonlinear optical response of a metal film, we consider its interaction with an external potential ϕ ^ext that describes excitation by an evanescent field originating at a distance z ₀ above the film. Translational symmetry in R then facilitates the decomposition of ϕ ^ext in Fourier components of the in-plane optical wave vector Q. The optical response of conduction electrons in the film is characterized by the total electrostatic potential ϕ = ϕ ext + ϕ b ind + v ⋅ ρ ind , where ϕ b ind accounts for the potential produced by a background induced polarization due to core electron screening, as discussed in ref [51], while the rightmost term in the potential results from the self-consistency of the induced charge,

mediated by the Coulomb interaction v(Q, z, z′). We use an analytical expression for v including the contribution of the polarization background [51], whereby the latter is calculated from the response of slab described by a local permittivity that is in turn obtained from the measured bulk dielectric function after eliminating the contribution of conduction electrons in the Drude model. The above expression is a sum over the diagonal elements of the real-space density matrix ρ(r, r, t) = ∑ _ii ρ _ii (t)Ψ _i (r)Ψ _i (r), expanded in the state basis {i, i′} with time-dependent matrix elements ρ _ii (t), and including a factor of 2 to account for spin degeneracy.

4.3 Perturbative solution of the density matrix

Expanding the density matrix as a perturbation series in ϕ ^ext, which is assumed to have a harmonic time dependence with frequency ω, we isolate the matrix elements of the nth-order contribution as

where, invoking in-plane momentum conservation, we have introduced the anzatz ρ j j ′ , k ∥ k ∥ ′ ( n , s ) → ρ j j ′ , k ∥ ( n , s ) δ k ∥ ′ , k ∥ − s Q , and have expanded in harmonics s up to |s| ≤ n. Given the self consistency of ϕ due to its dependence on ρ, we similarly write the potential as

Inserting the above expressions into Eq. (5) and equating terms of the same order n and harmonic s, we find a general expression for the matrix elements

where the first term describes the contribution from the self-consistent potential and

acts as a source term constructed from the response at lower perturbation orders. The induced charge density characterizing the optical response at order n and harmonic s is constructed according to

In practice, we compute the matrix elements ρ j j ′ k ∥ ( n , s ) by projecting in sinusoidal basis functions s _l (z) ≡ (2/L)^−1/2 sin(πlz/L), indexed by l = 1, 2, … and spanning a simulation domain in z of width L > d = Na _s extending beyond the actual film thickness. The induced charge density can be expressed in this basis as ∑ l ρ l ( n , s ) s l ( z ) e i s ( Q ⋅ R − ω t ) , with coefficients

where we have defined s l , j j ′ ≡ ∫ d z s l ( z ) φ j ( z ) φ j ′ * ( z ) . We obtain ρ l ( n , s ) by expressing Eq. (6) in terms of matrices indexed by l and l′ as

where

and

is the noninteracting RPA susceptibility evaluated at frequency sω and optical wave vector s Q, while

is the source potential expressed in the sinusoidal basis, with the external part contributing only to drive the linear response. Further details on the perturbative approach outlined above are provided in the SI.

4.4 Determination of the electron effective mass

In the QM model, we characterize in-plane electron motion in the quantized state j by the parabolic dispersion ℏ ε j , k ∥ = ℏ ε j ⊥ + ℏ 2 k ∥ 2 / 2 m j , adopting band-dependent electron effective masses m _j . These effective masses are fitted to ab-initio band structure calculations using the Vienna ab initio simulation package (VASP) [63–65]. More specifically, the effective mass m _j of state j is obtained by fitting ℏ k ∥ 2 / 2 m j to the dispersion of that state at the Fermi level E _F (i.e., in the neighborhood of electronic transitions dominating the low-energy plasmonic response), since plasmons are primarily comprised of virtual electron–hole pair transitions among one-electron states with small separation compared with the plasmon energy [36]. It should be noted that, while the first principles electronic bands are well captured by fitting parabolas within a few eV of the Fermi level, the bands deviate significantly from parabolic profiles at energies beyond this range, where they affect negligibly the optical response for the range of optical energy and momentum under consideration. Further details on the implementation of the numerical simulation can be found in the SI.

4.5 Classical model of a two-dimensional film

In a classical quasistatic framework, we approximate the ultrathin film as a two-dimensional layer, for which the nth-order optical response generating a harmonic s with in-plane wave vector component Q and frequency ω is characterized by an induced charge density ρ ^(n,s)(z) = ρ ^(n,s) δ(z) (i.e., in the z = 0 plane) that generates the potential ϕ ^(n,s)(z) = ∫dz′v ^(s)(z, z′)ρ ^(n,s)(z′) = (2π/sQ)e^−sQ|z| ρ ^(n,s). Invoking the continuity equation ∇ ⋅J ^(n,s) + ∂ _t ρ ^(n,s) = 0, we can then write the potential as

where J ^(n,s) is the total induced surface current.

To linear order, the total potential resulting from an evanescent source with amplitude ϕ ₀ located at z ₀ > 0 is

where r _p(Q, ω) is the reflection coefficient of the film, obtained in the Fabry–Pérot formalism as

with r dm = − r md = ( ϵ m − ϵ d ) ( ϵ m + ϵ d ) − 1 and t dm t md = 4 ϵ d ϵ m ( ϵ m + ϵ d ) − 2 . In the calculations here presented, we consider self standing films (i.e., ϵ _d = 1) and take the dielectric function of the metal from experimental values [55].

For nonlinear optical processes, we construct the current J = ∂ _t P from the polarization P ( n , s ) = χ ( 1 , s ) E ( n , s ) + P NL ( n , s ) , where the first term accounts for the linear response of the film (via the linear susceptibility χ ^(1,s)) to the nonlinear induced field and the second term is the source of the nonlinearity, proportional to the nonlinear susceptibility χ ^(n,s), and more explicitly, given by

for SHG, THG, and Kerr nonlinearity. The associated potential is then given from Eq. (7) as

where η CL ( n , s ) plays a similar role to that in Eq. (6) in the QM model, and is obtained for SHG, THG, and the Kerr nonlinearity as

with χ ̃ ( n , s ) = χ ( n , s ) d incorporating the film thickness d. For SHG, the use of a two-dimensional model is meaningful because only the surface contributes due to the centrosymmetry of the bulk metal; here, we use a local surface second-order susceptibility of the form

where e > 0 and m _e are the electron charge and mass, respectively, and a(ω) is a standard nonlinear parameter, for which we adopt the dispersionless mean value a = 7.75 reported in ref [66].

To describe third-order processes, we also use a 2D susceptibility and incorporate the reported bulk value [2] χ ⁽³⁾ = 2.9 × 10⁻¹⁹ V²/m². Note that this model for third order processes is dispersionless, neglects nonlocality, and further cannot account for (weak) cascaded nonlinear processes such as ∝ ϕ ²² ϕ ¹¹, which are included to the QM model. However, it should be reasonable when the response is dominated by plasmons at the fundamental frequency.