Viscoelastic optical nonlocality of low-loss epsilon-near-zero nanofilms

Optical nonlocalities are elusive and hardly observable in traditional plasmonic materials like noble and alkali metals. Here we report experimental observation of viscoelastic nonlocalities in the infrared optical response of epsilon-near-zero nanofilms made of low-loss doped cadmium-oxide. The nonlocality is detectable thanks to the low damping rate of conduction electrons and the virtual absence of interband transitions at infrared wavelengths. We describe the motion of conduction electrons using a hydrodynamic model for a viscoelastic fluid, and find excellent agreement with experimental results. The electrons’ elasticity blue-shifts the infrared plasmonic resonance associated with the main epsilon-near-zero mode, and triggers the onset of higher-order resonances due to the excitation of electron-pressure modes above the bulk plasma frequency. We also provide evidence of the existence of nonlocal damping, i.e., viscosity, in the motion of optically-excited conduction electrons using a combination of spectroscopic ellipsometry data and predictions based on the viscoelastic hydrodynamic model.

where , , is the distribution function, the microscopic electron velocity, m the electron mass, the Lorentz force acting on the electrons, and is an operator that accounts for collisions (with ions and between electrons). The next step is the conversion of the Boltzmann equation into a set of fluid-dynamic equations that describe the time evolution of macroscopic quantities, i.e., the velocity moments of the distribution function. These moments are measurable, macroscopic variables such as temperature, pressure, mean velocity, current density, and stress tensor. We consider the first three velocity moments. The first is the electron density, which is the zero-order velocity moment obtained by integrating the distribution function in the velocity space, , , , . The second is the first-order velocity moment divided by the electron density, i.e., the average (macroscopic) velocity, , . The (macroscopic) current density is simply , , , , , , and it is approximated to the first order , , in the main text. For higher-order moments, it is useful to transform the distribution function in the Lagrange comoving frame, i.e., , , , , . In this frame, the macroscopic velocity is zero and the second-order velocity moment, i.e., the pressure (stress) tensor, is written as , ⊗ , , .
Next, the equations of motion for these macroscopic quantities (velocity moments) are retrieved by transforming the Boltzmann equation in an infinite chain of coupled equations for the moments of the distribution functions. The chain is truncated to the second order as follows.
The first equation (zero-order moment equation) is obtained by integrating the Boltzmann equation over the velocity space and assuming collisions with conservation of charge, mass and electron density (no recombination), hence 0. This leads to the usual continuity equation, (3) in the main text is the linearization of this equation around the equilibrium, i.e., it is obtained by expanding the density as . The first-order moment equation is obtained by multiplying the Boltzmann equation by and integrating in the velocity space. This yields the transport equation: In the relaxation time approximation, the first moment of the collision term assumes the simple expression , in which represents scattering of electrons with ions and impurities. Since is symmetric, the isotropic part of the pressure tensor (i.e., the scalar pressure P) is separated from the viscous term as follows: , where is the identity matrix, is the trace of and is the traceless part of the stress tensor. Eq. (1) in the main text is the linearization of equation (S2) around the equilibrium, i.e., expanding density and pressure as and . Transport effects due to the magnetic Lorentz force and the convective term ⋅ play a central role in harmonic generation problems, but they can be neglected in our treatment in which only linear, nonlocal phenomena are investigated. Eq. (S2) generalizes classical hydrodynamic theories by introducing viscosity, therefore it describes electrons as a viscous fluid. It is important to remember that other classical approaches to describe the motion of conduction electrons can be derived from the Boltzmann equation by adopting different approximations. For example, the local Drude model treats conduction electrons in plasmonic materials as an incompressible gas, in which both elasticity and viscosity are neglected, i.e., in eq. (S2). The Bloch hydrodynamic theory, another classical and widely adopted model, considers only the isotropic part of the stress tensor, i.e., hence and 0.
Finally, the second-order moment equation (conservation of energy) is retrieved by multiplying the Boltzmann equation, eq. (S1), by ⊗ and integrating over the velocity space: where the second-order collision term is ⊗ . , and the right-hand side of eq. (S3) becomes ⊗ .
Since the Fermi pressure , the collision moment is simply , i.e., it only acts on the traceless part of the stress tensor.
The trace of eq. (S3) reads as while the traceless, remaining part of eq. (S3) is Considering for the kinetic pressure an expression equal to that of a degenerate Fermi gas at T = 0 K, 2 2 2/3 5/3 * (3 ) 5 P n m    , the term ⋅ ⋅ vanishes in eq. (S4). This means that the linearized version of eq. (S4) is equal to eq. 2(a) of the main text.
Terms that depend on the product of by in eq. (S5) vanish in our perturbative, linear approximation, therefore eq. (S5) reduces to eq. 2(b) of the main text.

Supplementary Note 2: Additional observations of higher-order modes
The trend of higher-order resonances has been consistently observed in films with thickness smaller than 25 nm. In Fig. S1 we show the ellipsometric quantity Ψ tan | |/| | for three additional samples having different thicknesses of 16.2 nm, 21.3 nm, and 38.8 nm. The theoretical fit obtained with the nonlocal viscoelastic model is reported together with the measured data. The trend is in agreement with Fig. 2 of the main manuscript: the larger-thickness film of 38.8 nm shows little nonlocal behavior (blue curve), hence it displays no higher-order resonances, while the two thinner films show one or two additional higher-order resonances (green and red curves). Fig. S1. On the left, measured and calculated (via nonlocal theory) spectra of the ellipsometric parameter Ψ (reported in degrees) for three samples. On the right, zoom of the measured spectrum for the d = 16.2 nm to highlight the resonances at λ 3 and λ 5 . The properties of the three samples are summarized in table S1.
We stress that these samples were not prepared for ATR, therefore the spectra were taken from the air side and the modes appear as shallower resonances (see also Fig. 1 in the manuscript). The enlarged view of the thinnest-film spectrum clearly shows the two additional resonances (right panel of Fig. S1). As for the samples presented in the manuscript, film thicknesses were measured with X-ray reflectivity (XRR), while carrier concentrations and mobilities are determined with Hall-effect measurements. The film properties are summarized in Table S1.