Nonlinear plasmonic response in atomically thin metal films

Abstract Nanoscale nonlinear optics is limited by the inherently weak nonlinear response of conventional materials and the small light–matter interaction volumes available in nanostructures. Plasmonic excitations can alleviate these limitations through subwavelength light focusing, boosting optical near fields that drive the nonlinear response, but also suffering from large inelastic losses that are further aggravated by fabrication imperfections. Here, we theoretically explore the enhanced nonlinear response arising from extremely confined plasmon polaritons in few-atom-thick crystalline noble metal films. Our results are based on quantum-mechanical simulations of the nonlinear optical response in atomically thin metal films that incorporate crucial electronic band structure features associated with vertical quantum confinement, electron spill-out, and surface states. We predict an overall enhancement in plasmon-mediated nonlinear optical phenomena with decreasing film thickness, underscoring the importance of surface and electronic structure in the response of ultrathin metal films.

We provide a detailed derivation of the main equations used in the main text, based on a classical electrostatic description of the fields, a self-consistent treatment of conduction electrons in the films, and perturbation theory beyond linear order to compute each of the nonlinear optical response functions. Additional details on the calculation of the effective masses are also provided.
Contents S1. Electronic states and density matrix 1 S2. Perturbation series expansion 2 A. Linear response 2 B. Nonlinear response 2 C. Self-consistent optical response calculations through expansion in an orthonormal basis 3 S3. Electronic structure of crystalline metal films 4 A. Atomic layer potential 4 B. Characterization of effective masses 5 C. First-principles description of effective masses 5 References 6

S1. ELECTRONIC STATES AND DENSITY MATRIX
We consider a metallic film that is translationally invariant in the R = (x, y) plane and finite in the zdirection, so that electrons in the unperturbed film are confined by a potential V (z) entering the Hamiltonian which admits a complete set of separable wave functions where A is a normalization area, k is the in-plane electron momentum, and ϕ j (z) are localized wave functions * Electronic address: cox@mci.sdu.dk † Electronic address: javier.garciadeabajo@nanophotonics.es from which the obtained energy eigenvalues ε ⊥ j contribute to the total electron energy ε j,k = ε ⊥ j + 2 k 2 /2m j after having introduced band-dependent effective masses m j to the parabolic in-plane dispersion.
In the quasistatic approximation, the thin film interacts with an electrostatic potential φ = φ ext +φ ind b +φ ind , comprised of the external perturbation a background polarization potential due to core electron screening, as reported in Ref. 1, and the induced potential arising from the polarization of conduction electrons. We characterize the external potential by a fixed in-plane wave vector component Q and frequency ω, while the induced potential is written in terms of the induced charge ρ ind (r, t) = −2e[ρ(r, r, t) − ρ (0) (r, r)], which is in turn obtained from the diagonal elements of the single-particle density matrix where ρ jj ,k k (t) are time-dependent coefficients, the dynamics of which is governed by the Liouville von Neumann equation The above expression incorporates inelastic electron scattering by relaxing the density matrix at a phenomenological rate γ to its equilibrium value where f j,k are occupation factors following Fermi-Dirac statistics. In practice we adopt the zero-temperature limit, so that f j,k = Θ(E F − ε j,k ), with Θ(x) denoting the Heaviside step function and E F the Fermi energy.

S2. PERTURBATION SERIES EXPANSION
Inserting Eq. (S6) into Eq. (S7), we obtain where i = {j, k } denotes a multiplexed state index that simplifies the notation, and the matrix elements of the potential are The density matrix and total potential are each expanded in perturbation series with respect to the external potential according to where the latter is obtained from the former through Eq. (S5). Inserting Eq. (S11) into Eq. (S9) and equating terms of equal perturbation order, we find which reveals that the equation of motion for the density matrix at perturbation order n is constructed according to (i.e., self-consistently through the induced potential). The latter depends on the density matrix through the induced charge, and also trough the density matrix obtained at lower perturbation orders n < n. We note that the right-hand side of Eq. (S12) disappears at n = 0 order, thus reducing the expression to ρ (0) ii = f i δ ii and satisfying Eq. (S8).

A. Linear response
To first order, the right-hand side of Eq. (S12) only has a contribution ∝ φ (1) ρ (0) , while the external potential in Eq. (S4) depends on in-plane position and time according to e is(Q·R−ωt) , where s = ±1 is the harmonic index. Then, the integral over R in Eq. (S10) for φ ext effectively yields a term for each s proportional to the Kronecker delta δ k ,k −sQ , directly enforcing parallel momentum conservation, and containing a time dependence e −isωt . We thus conclude that the density matrix at linear order in the external potential inherits the same momentum and frequency dependence, allowing us to write Inserting the above expression for the total potential into Eq. (S10), we isolate the coefficients Here, we have introduced the sQ in-plane Fourier component of the Coulomb interaction v (s) (z, z ) and the linear induced charge density at harmonic s We now insert Eq. (S13) into Eq. (S12) to isolate the coefficients from which the linear induced charge density is obtained.

B. Nonlinear response
The method employed in the previous section to extract the linear response of the thin film apply generally to nonlinear processes of order n > 1 and harmonic index s, with |s| ≤ n. In particular, following Eq. (S13), we write which upon insertion into Eq. (S10) yields the matrix elements Here, are the amplitudes of the potential and induced charge, respectively, which depend on the film confinement coordinate z.
Following analogous steps as those used to obtain the first-order density matrix, Eq. (S12) becomes where the matrix elements of the potential are We now follow a computationally efficient method to calculate ρ (n,s) by first isolating the self-consistent term in the summation of Eq. (S21) according to where we have used Eq. (S8) to write ρ The latter collects the dependence on perturbation orders n < n, which drive the nonlinear response. C. Self-consistent optical response calculations through expansion in an orthonormal basis In practical calculations, we expand quantities depending on z in a complete basis set of orthonormalized sine functions where l is a positive integer and L is the length of the simulation domain. In terms of the above basis functions, we write the induced charge as where the expansion coefficients are Here, we have transformed the sum over k into an integral according to the prescription A −1 k → dk /(2π) 2 . In a similar fashion, we write the potential as where and We identify the χ (0,s) as the noninteracting susceptibility in the random-phase approximation (RPA) evaluated at harmonic index s.
In practice, the nonlinear potential at order n and harmonic s is computed self-consistently by first constructing η (n,s) from lower-order terms by means of Eqs. (S24) and (S32). Following a standard procedure in the context of RPA calculations, we obtain the induced potential by expressing Eq. (S31) in matrix form and multiplying by v (s) to isolate where the dot denotes matrix multiplication and we have introduced the Coulomb interaction in the sine basis v (s) as well as the quantity which acts as a source potential.

S3. ELECTRONIC STRUCTURE OF CRYSTALLINE METAL FILMS
The theoretical formalism presented thus far to describe the self-consistent-field optical response of thin films is independent of the confinement potential V (z) appearing in Eq. (S1). In what follows, we elaborate on the methods employed to treat electrons in crystalline metal films incorporating details of atomic corrugation in the out-of-plane direction.

A. Atomic layer potential
To describe a silver crystalline noble metal film of thickness d, we employ the so-called atomic layer potential introduced in Ref. [1], where the potentials V surf are those reported in Ref. [2], which capture the electronic structure of noble metal surfaces in a specific crystallographic orientation by treating them as a semi-infinite stack of atomic planes with interatomic layer spacing a s = 0.2361 terminating at z = 0: Here, the coefficients A 1 = 4.30 eV and A 10 = −9.640 eV are chosen to reproduce the width and position of the noted energy gap, respectively; the space between z = z i and z = 0 represents the transition from the solid bulk to free-space, where electron spill out takes place; and the parameters A 2 = 3.8442 and η = 48.470 nm −1 determine the positions of the Fermi level and the surface states relative to the vacuum level. Five of the remaining six parameters are determined by imposing the continuity of the potential and its first derivative, so that A 20 = A 2 − A 1 − A 10 , A 3 = −A 20 + A 2 cos (ηz t ), α = ηA 2 sin (ηz t )/A 3 , λ = 2α, and z i = α −1 log −λ/4A 3 +z t , while the intermediate point z t = −5π/4η is fixed with respect to the surface atomic layer in such a way that z i corresponds to the image plane position, which is important for describing image states [2]. The atomic layer potential of Eq. (S37) then describes a film with outermost atomic planes located at z = a s /2 and z = d − a s /2, and is continuous at z = d/2 by construction. Obviously, d must be taken to be an integer multiple of the atomic-layer spacing a s .

B. Characterization of effective masses
In order to determine the band-dependent effective masses m j for a thin metal film comprised of a specific number of atomic planes, we examine the energy dispersion bands E j (k ) extracted from ab initio calculations (see below). The optical response of the metal is mainly determined by the behavior of electrons at the Fermi level E F , and consequently, we compute the effective mass at the Fermi wave vector k = k F by means of a parabolic fit according to where the second derivative is calculated numerically from the generated bands. Figure S1a shows the electronic binding energies ε ⊥ j for Ag(111) films with respect to the vacuum level as a function of the number of (111) atomic planes, as obtained from the diagonalization of the Hamiltonian in Eq. (S3) when using the atomic layer potential [1,2]. Here, electronic states are represented by color-coded dots, with bulk states following the thickness-dependent color scale in (b) and surface states (SS) indicated in light blue and red. The binding energies of the latter exhibit a different evolution with thickness that we attribute to the concentration of their wave functions in the vicinity of the metal-dielectric interface. In films consisting of just a few atomic layers, the surface states undergo strong hybridization and energy splitting that is substantially reduced as the film thickness increases to 20 atomic planes, with the energy increasingly approaching the semi-infinite surface value of 0.026 meV below the Fermi energy, indicated by the solid black line.
The effective masses computed for each of the considered thicknesses in Figure S1b, are observed to converge towards a particular shape as a function of their binding energies determined in Figure S1a. In Figure S1b, the effective masses of films with thickness ranging from 5 to 40 Ag(111) atomic planes are presented, showing convergence toward the upper end.
Within the valence band, states closer to the Fermi level exhibit a linear variation of the effective mass, in agreement with previous experiments [3]. However, the variation in effective mass for states with energies −6 eV deviates from the linear trend, motivating the use of density functional theory in the QM calculations with fitted masses. For an arbitrary number of atomic planes, we find a polynomial fit of the effective mass as a function of the number of atomic planes N to exhibit extreme sensitivity to the precision of the fitting coefficients, whereas a more robust result is obtained by using the functional form with coefficients A 1 = 0.457, A 2 = −1.442, B 1 = 0.831 eV −1 , B 2 = 0.319 eV −1 , and C = −0.113, and with binding energies ε ⊥ j computed by diagonalizing the Hamiltonian in Eq. (S3). This expression is plotted as a black solid curve in Figure S1b.
For surface sates, the variation of effective mass as a function of the number of atomic planes is shown in Fig-ure S1c, revealing a decrease in effective mass with increasing number of layers that converges to m/m e = 0.39 for Ag(111) films thicker than 20 layers, in agreement with reported experimental values [4]. In this work, we fit the effective mass variation as a function of the number of atomic planes N as m ss m e = a j e −bj N + c, with j referring the j = 1 or j = 2 surface state (SS), for which we have coefficients a 1 = 0.4996, a 2 = 1.391, b 1 = 0.7675, b 2 = 0.2689, and c = 0.39. We note that the color code of panel Figure S1c coincides with the SS representation in Figure S1a.

C. First-principles description of effective masses
We perform density functional theory (DFT) calculations using the Vienna ab initio simulation package (VASP) [5][6][7] using the generalized gradient approximation (GGA) within the Perdew-Burke-Ernzerhof (PBE) formalism [8] for the electron exchange and correlation potential. The interatomic distances are all fixed at the optimized values for the bulk crystal, since relaxation of the silver-silver interatomic distance does not significantly change the electronic band structure, even in films comprised of only a few (N ∼ 5) atomic planes. Such effects are more important in the outermost surface layers where the interfaces of the film are located, and it is precisely the atomic layer potential model introduced above that incorporates in a phenomenological way the description of the surface fitted to reproduce experimental band parameters.
The cutoff energy of the plane-wave basis set is set to 600 eV. The convergence criterion for the energy is taken as 10 −5 eV and the total force in the unit cell is reduced to a value of less than 10 −4 eV/Å in the optimization of the bulk interlayer distance. A Γ-centered k-point mesh of size 18×18×1 is used for the Brillouin zone sampling.
A vacuum space of 10Å is appended in order to prevent interactions between neighboring layers in the vertical supercell.