We simulated nonlinear light–matter interactions using a theoretical model described in detail elsewhere.17,19,29,35,39–4117. M. Scalora, J. Trull, C. Cojocaru, M. A. Vincenti, L. Carletti, D. de Ceglia, N. Akozbek, and C. De Angelis, “Resonant, broadband, and highly efficient optical frequency conversion in semiconductor nanowire gratings at visible and UV wavelengths,” J. Opt. Soc. Am. B 36(8), 2346–2351 (2019). https://doi.org/10.1364/josab.36.00234619. T. S. Luk, D. de Ceglia, S. Liu, G. A. Keeler, R. P. Prasankumar, M. A. Vincenti, M. Scalora, M. B. Sinclair, and S. Campione, “Enhanced third harmonic generation from the epsilon-near-zero modes of ultrathin films,” Appl. Phys. Lett. 106, 151103 (2015). https://doi.org/10.1063/1.491745729. M. Scalora, M. A. Vincenti, D. de Ceglia, V. Roppo, M. Centini, N. Akozbek, and M. J. Bloemer, “Second- and third-harmonic generation in metal-based structures,” Phys. Rev. A 82, 043828 (2010). https://doi.org/10.1103/physreva.82.04382835. L. Rodríguez-Suné, J. Trull, C. Cojocaru, N. Akozbek, D. de Ceglia, M. A. Vincenti, and M. Scalora, “Harmonic generation from gold nanolayers: Bound and hot electron contributions to nonlinear dispersion,” Opt. Express 29, 8581–8591 (2021). https://doi.org/10.1364/oe.41545939. M. Scalora, M. A. Vincenti, D. de Ceglia, C. M. Cojocaru, M. Grande, and J. W. Haus, “Nonlinear duffing oscillator model for third harmonic generation,” J. Opt. Soc. Am. B 32, 2129–2138 (2015). https://doi.org/10.1364/josab.32.00212940. L. Rodríguez-Suné, M. Scalora, A. S. Johnson, C. Cojocaru, N. Akozbek, Z. J. Coppens, D. Perez-Salinas, S. Wall, and J. Trull, “Study of second and third harmonic generation from an indium tin oxide nanolayer: Influence of nonlocal effects and hot electrons,” APL Photonics 5, 010801 (2020). https://doi.org/10.1063/1.512962741. M. Scalora, J. Trull, D. de Ceglia, M. A. Vincenti, N. Akozbek, Z. Coppens, L. Rodríguez-Suné, and C. Cojocaru, “Electrodynamics of conductive oxides: Intensity-dependent anisotropy, reconstruction of the effective dielectric constant, and harmonic generation,” Phys. Rev. A 101, 053828 (2020). https://doi.org/10.1103/physreva.101.053828 The classical description of free and bound electron current densities, defined as Jf=Ṗf and Jbj=Ṗbj, respectively, leads to the following scaled equations of motion for conduction and core electrons:P̈f+γ̃fṖf=e2λ02n0,fm0*c2E+Λ̃(E⋅E)E−eλ0m0*c2(∇⋅Pf)E+eλ0m0*c2Ṗf×H+35EFm0*c2∇(∇⋅Pf)+12∇2Pf−1n0,feλ0[(∇⋅Ṗf)Ṗf+(Ṗf⋅∇)Ṗf],(1)P̈bj+γ̃bjṖbj+ω̃0,bj2Pbj−β̃j(Pbj⋅Pbj)Pbj+Θ̃j(Pbj⋅Pbj)2Pbj−Δ̃j(Pbj⋅Pbj)3Pbj=n0,be2λ02mbj*c2E+eλ0mbj*c2(Pbj⋅∇)E+eλ0mbj*c2Ṗbj×H.(2)

The subscript b stands for bound, and the counter j represents the jth bound oscillator species. The linear dielectric response of gold is modeled using one Drude component and three separate Lorentzian functions to represent the dielectric response more accurately in the 200–300 nm range.

Equation (2) thus represents three separate equations, each describing an oscillator species. The local dielectric constant may be written asε(ω)=1−ω̃pf2ω2+iγ̃fω−ω̃p12ω2−ω̃012+iγ̃01ω−ω̃p22ω2−ω̃022+iγ̃02ω−ω̃p32ω2−ω̃032+iγ̃03ω,(3)where (ω̃pf,ω̃p1,ω̃p2,ω̃p3)=(7.1,3.4,4.79,6.35), (γ̃f,γ̃01,γ̃02,γ̃03)=(0.05,1.25,1.45,1.25), and (ω̃0f,ω̃01,ω̃02,ω̃03)=(0,2.45,3.45,4.75). Equation (3) is obtained by Fourier transforming Eqs. (1) and (2) after dropping all nonlinear and nonlocal terms. The frequency ω=1λ, where λ is expressed in micrometers. Equations (1) and (2) amount to an expanded hydrodynamic model that includes a description of surface [(∇ · Pf)E], magnetic (Ṗf×H), convective [(∇⋅Ṗf)Ṗf+(Ṗf⋅∇)Ṗf], nonlocal [∇(∇⋅Pf)+12∇2Pf], and hot electron [Λ̃(E⋅E)E] nonlinearities in the free electron polarization component.Furthermore, Eq. (2) describes the surface [(Pbj · ∇)E], magnetic (Ṗbj×H), and bulk [−β̃j(Pbj⋅Pbj)Pbj+Θ̃j(Pbj⋅Pbj)2Pbj−Δ̃j(Pbj⋅Pbj)3Pbj] nonlinearities triggered by bound electrons. In the bulk response, we neglect even order nonlinearities because the system is centrosymmetric. We choose to expand these bound electron bulk terms up to the seventh order (i.e., Pbj7) because the local fields can be quite large near edges and points. For instance, Fig. 2(b) suggests a local electric field intensity enhancement by a factor of nearly 400. As a result, just a few GW/cm2 incident onto the sample can turn into several TW/cm2, thus necessitating going beyond a mere third order response. In fact, the higher order terms shown also contribute to odd harmonics. For instance, a few of the many non-instantaneous terms contributing to the full, complex χω(3) response are proportional to β̃j|Pω|2Pω, Θ̃j|Pω|4Pω, and Δ̃j|Pω|6Pω, where β̃j=ω0,j2λ02L2n0b2e2c2, Θ̃j=β̃jn0b2e2L2, and Δ̃j=β̃jn0b4e4L4. Similarly, some of the terms contributing to the full, complex χ3ω(3) are proportional to β̃jPω3, Θ̃j|Pω|2Pω3, and Δ̃j|Pω|4Pω3. It should be evident that the approach outlined earlier may be considered more accurate than the usual, instantaneous responses of types χω(3)|Eω|2Eω and χ3ω(3)Eω3, which require prior knowledge of nonlinear dispersion. On the other hand, β̃j, Θ̃j, and Δ̃j are real, known, and depend on the linear oscillator parameters like resonance frequency ω0,j and lattice constant L. These considerations about free and bound electron third and higher order nonlinearities are important because they can cause a dynamic shift of the plasmonic resonance, either toward the red or the blue, depending on tuning relative to the bound electron resonances.3535. L. Rodríguez-Suné, J. Trull, C. Cojocaru, N. Akozbek, D. de Ceglia, M. A. Vincenti, and M. Scalora, “Harmonic generation from gold nanolayers: Bound and hot electron contributions to nonlinear dispersion,” Opt. Express 29, 8581–8591 (2021). https://doi.org/10.1364/oe.415459 In Eqs. (1) and (2), n0,f and n0,b are the free and bound electron densities, respectively, in the absence of an applied field; the subscript j in Eq. (2) indicates the presence of multiple bound electron species, three in our case, which are for simplicity assumed to have the same particle density n0,b; e is the electronic charge; Ef is the Fermi energy; λ0 = 1 µm is a conveniently chosen scaling wavelength; c is the speed of light in vacuum; m0* and mbj* are the free and bound electron effective masses; and ξ=zλ0,η=yλ0,ζ=xλ0,τ=ctλ0 are the scaled space and time coordinates that we used to calculate spatial and temporal derivatives. Finally, plasma and resonance frequencies and damping coefficients are scaled as follows: ω̃p2=4πne2m*λ02c2,ω̃0=ω0λ0c,γ̃=γλ0c. Here, we emphasize that in the wavelength range below 1 µm, it is essential to model the metal as a collection of free and bound charges as the incident wave begins to probe atomic orbitals that reside below the free electron cloud. Once the dielectric function data4242. E. D. Palik, Handbook of Optical Constants of Solids (Academic Press, London; New York, 1985). are fitted using the dielectric function in Eq. (3), the ratios n0,fm0* and n0,bmbj* are set and suffice to allow the linear portion of the dynamics to evolve in both Eqs. (1) and (2). However, the nonlinear coefficients require direct knowledge of both n0,f and m0* in Eq. (1) and of mbj* in Eq. (2). For this purpose, we choose the nominal values m0*=mbj*=me and the known gold free electron density n0,f ∼ 5.8 × 1022 cm−3. These coefficients thus determine the relative magnitudes of three competing free-electron sources of SHG: the Coulomb (surface) term, the magnetic Lorentz (both surface and volume) term, convective (surface) contributions, and the Fermi energy: EF=ℏ22m0*(3π2n0,f)23, which lands at ∼5 eV. Once these coefficients are fixed, we determine the lattice constant using average orbital diameters. For example, in the case of gold,3535. L. Rodríguez-Suné, J. Trull, C. Cojocaru, N. Akozbek, D. de Ceglia, M. A. Vincenti, and M. Scalora, “Harmonic generation from gold nanolayers: Bound and hot electron contributions to nonlinear dispersion,” Opt. Express 29, 8581–8591 (2021). https://doi.org/10.1364/oe.415459 the diameter of the p-orbital (which determines the amplitude of the oscillation of bound electrons and thus the magnitude of L) is ∼0.5 Å. Since each gold atom contributes approximately one free electron to the electron density, it is reasonable to assume that the bound electron density may be approximated by the free electron density. This is then sufficient to determine β̃j,Θ̃j,Λ̃j. Finally, Λ̃j is a scaled perturbative parameter that measures changes in the free electron density as a function of incident light intensity. It is generally dispersive and depends on field fluence and absorption. For simplicity, we parameterize it and assume it to be nearly constant in the range of wavelengths of interest. For additional, specific details about the model and method of solution, we refer the reader to Refs. 1717. M. Scalora, J. Trull, C. Cojocaru, M. A. Vincenti, L. Carletti, D. de Ceglia, N. Akozbek, and C. De Angelis, “Resonant, broadband, and highly efficient optical frequency conversion in semiconductor nanowire gratings at visible and UV wavelengths,” J. Opt. Soc. Am. B 36(8), 2346–2351 (2019). https://doi.org/10.1364/josab.36.002346, 1919. T. S. Luk, D. de Ceglia, S. Liu, G. A. Keeler, R. P. Prasankumar, M. A. Vincenti, M. Scalora, M. B. Sinclair, and S. Campione, “Enhanced third harmonic generation from the epsilon-near-zero modes of ultrathin films,” Appl. Phys. Lett. 106, 151103 (2015). https://doi.org/10.1063/1.4917457, 2929. M. Scalora, M. A. Vincenti, D. de Ceglia, V. Roppo, M. Centini, N. Akozbek, and M. J. Bloemer, “Second- and third-harmonic generation in metal-based structures,” Phys. Rev. A 82, 043828 (2010). https://doi.org/10.1103/physreva.82.043828, 3535. L. Rodríguez-Suné, J. Trull, C. Cojocaru, N. Akozbek, D. de Ceglia, M. A. Vincenti, and M. Scalora, “Harmonic generation from gold nanolayers: Bound and hot electron contributions to nonlinear dispersion,” Opt. Express 29, 8581–8591 (2021). https://doi.org/10.1364/oe.415459, and 39–4139. M. Scalora, M. A. Vincenti, D. de Ceglia, C. M. Cojocaru, M. Grande, and J. W. Haus, “Nonlinear duffing oscillator model for third harmonic generation,” J. Opt. Soc. Am. B 32, 2129–2138 (2015). https://doi.org/10.1364/josab.32.00212940. L. Rodríguez-Suné, M. Scalora, A. S. Johnson, C. Cojocaru, N. Akozbek, Z. J. Coppens, D. Perez-Salinas, S. Wall, and J. Trull, “Study of second and third harmonic generation from an indium tin oxide nanolayer: Influence of nonlocal effects and hot electrons,” APL Photonics 5, 010801 (2020). https://doi.org/10.1063/1.512962741. M. Scalora, J. Trull, D. de Ceglia, M. A. Vincenti, N. Akozbek, Z. Coppens, L. Rodríguez-Suné, and C. Cojocaru, “Electrodynamics of conductive oxides: Intensity-dependent anisotropy, reconstruction of the effective dielectric constant, and harmonic generation,” Phys. Rev. A 101, 053828 (2020). https://doi.org/10.1103/physreva.101.053828. The integrations are carried out using a split-step method coded in Fortran that advances the fields and material equations separately in space and time. The field equations are solved in free space using a fast Fourier transform algorithm, while the material equations are solved using a predictor-corrector algorithm.3939. M. Scalora, M. A. Vincenti, D. de Ceglia, C. M. Cojocaru, M. Grande, and J. W. Haus, “Nonlinear duffing oscillator model for third harmonic generation,” J. Opt. Soc. Am. B 32, 2129–2138 (2015). https://doi.org/10.1364/josab.32.002129 We use a spatial step δz = δy = δx = 1.25 nm and a correspondingly causal (δz = cδt) temporal integration step δt = 4 × 10−18 s. However, using δz = δy = δx = 2.5 nm and a correspondingly larger time step yields similar results. The increased spatial step size introduces a degree of realism to the problem that mimics fabrication imperfections such as rounded corners and roughness. The split-step pulse propagation method is notorious for its accuracy in calculating spatial derivatives accurately to all orders, unlike finite difference time domain (FDTD) methods, because it utilizes a Fourier expansion of the fields that uses hundreds of thousands of plane waves. We note that, unlike unconditionally stable spectral methods, ordinary finite difference time domain (FDTD) methods are subject to numerical stability criteria that, in a multidimensional grid, necessarily disrupt the causality condition, leaving open the possibility of introducing phase errors and numerical dispersion,4343. T. Kashiwa, H. Kudo, Y. Sendo, T. Ohtani, and Y. Kanai, “The phase velocity error and stability condition of the three-dimensional nonstandard FDTD method,” IEEE Trans. Magn. 38, 661–664 (2002). https://doi.org/10.1109/20.996172 especially in highly nonlinear environments. To our knowledge, these issues and concerns remain unaddressed.