Retrieving Linear and Nonlinear Optical Dispersions of Matter: Combined Experiment-Numerical Ellipsometry in Silicon, Gold and Indium Tin Oxide

Introduction

The study of nonlinear optical effects is crucial in many applications in photonics covering many different research fields, including quantum electronics, quantum optics, telecommunications, surface science, metrology, microscopy and many others. In order to properly account for nonlinear dispersion and its effects, including for example, self-phase modulation, nonlinear pump absorption, and harmonic conversion efficiencies, one must determine as accurately as possible the values of the different susceptibilities of the medium, which are ultimately related to material parameters like electron density and effective mass, lattice constant, and damping rates. The inclusion of nonlinear terms in the equations of motion thus leads to a total polarization density that accounts for linear and nonlinear responses. The nonlinear contributions of the total polarization at the frequency corresponding to the incident fundamental field give rise to a modified index of refraction, or dielectric constant of the medium, which for third order processes, under certain conditions may be written in terms of the intensity simply as: $n=n_0+n_{2}I$. The determination of $n_2$ and it dispersive properties are crucial to properly describe the nonlinear interactions and to obtain the nonlinear phase shift ${\varphi }_{NL}$ induced by the field intensity. When nonlinear interactions occur in a bulk dielectric, for instance, the value of $n_2$ can be directly related to the third-order susceptibility.

The issue of accurately determining the values of the nonlinear coefficients has been considered for many years. The different techniques involve the use of nonlinear interferometry, degenerate four-wave mixing, ellipse rotation, beam distortion and deflection measurements. The classical method to determine the value (magnitude and sign) of $n_2$ is the z-scan technique, developed by Sheik-Bahae et al. (1989), Sheik-Bahae et al. (1990). By moving the sample and measuring transmittance along the longitudinal direction on the focal plane of a focused Gaussian beam, this technique allows one to infer the nonlinear coefficient both in amplitude and sign.

Different z-scan theories have been proposed in the literature. The first theories considered thin samples and a Kerr nonlinearity to derive an analytical formula for weak nonlinearities (${\varphi }_{NL}<0.2\pi$), and a numerical estimation for larger nonlinear phase shifts considering beam propagation based on the nonlinear paraxial wave equation under the parabolic approximation. Other theories have been proposed to increase the range of applicability to longer samples based on the nonlinear paraxial wave equation, completed by the Huygens-Fresnel propagation (Pálfalvi et al., 2009).

Extensions of the z-scan technique have also been implemented in order to unravel the different contributions to the nonlinear response of the material. For instance, in time-resolved z-scan methods, the introduction of a temporal delay into the two-color z-scan device allows one to separate nonlinear contributions having different temporal responses (Wang et al., 1994). More recently, methods based on beam deflection have been implemented using pump-probe configurations, whereby the pump generates an index gradient experienced by the probe, which in turn is deflected (Ferdinandus et al., 2013; Ferdinandus et al., 2017). However, a detailed theory for either z-scan or beam deflection technique should take into account the origin and nature of all relevant nonlinearities present in the problem, in order to extract accurate information from experimental measurements, and thus obtain the desired, accurate values of the nonlinear coefficients not only for the pump beam, which either z-scan or beam deflection are mostly concerned with, but also for the simultaneous determination of nonlinear dispersion of the generated harmonic fields.

In more complex situations, where the material consists of metal layers, or perhaps semiconductor or conductive oxide layers, new linear and nonlinear sources become relevant, including nonlocal effects, magnetic dipole and electric quadrupole (surface) nonlinearities, convection, hot electrons, pump depletion, and phase-locking. As a result, the behavior of a harmonic component cannot be extrapolated based on the behavior of a pump or probe field grounded on mere transmission or deflection assessments. At the same time, one should be able to identify and distinguish between competing second and/or third order nonlinearities separately triggered by free and bound electrons. Therefore, different terms introduced in the equations of motion contribute to the nonlinear susceptibility and add to beam dynamics, terms that typically are either not distinguished or accounted for in the theories behind either z-scan or beam deflection methods. Consequently, the results may be either inaccurate or misleading in situations where these contributions may be significant. In what follows we describe a numerical technique that utilizes the constitutive relations to extract complex nonlinear dispersions and related coefficients of various orders. The approach is based on a microscopic, hydrodynamic representation of the material equations of motion coupled to the macroscopic Maxwell’s equations, to perform what amounts to numerical ellipsometry.

In the classical realm, electrons move according to Newtonian principles. Free electrons are described by a modified Drude equation of motion, one that contains linear and nonlinear, nonlocal, magnetic, convective, and quadrupole-like surface contributions (Scalora et al., 2010). In addition to similar external electric and magnetic forces, multiple bound electron species may also contribute to the macroscopic dielectric constant, and are subject to linear and nonlinear restoring forces (Scalora et al., 2012; Vincenti et al., 2011). For example, at near-IR, visible and UV wavelengths, in metals d-shell (bound) electrons may dominate the dielectric response, and their influence should be introduced in the dynamics, also because they contribute to the total balance of momentum and energy transfer. Each electron species may experience its own linear and nonlinear response, have different resonant frequencies, effective electron masses, and damping constants. The resulting nonlinear polarizations are then added together to form the total polarization, which in turn is inserted into Maxwell’s equations.

The depiction of nonlinear phenomena of metals, for example, is usually limited and focused only on the second order response of free electrons, with effective surface and bulk second order nonlinear coefficients (Xiang Wang et al., 2009; Krause et al., 2004) that are chosen and typically adjusted to fit experimental results (Coutaz et al., 1987) without distinguishing between surface (Coulomb and convection) and magnetic (Lorentz) contributions. When third order effects are introduced, the ${\chi }_{\omega }^{\left(3\right)}$ usually lacks information about dispersion (Suresh et al., 2021), is not always obtained under conditions that are directly pertinent to the problem at hand, and ${\chi }_{3\omega }^{\left(3\right)}$ remains undetermined. For instance, nanosecond pulses (Yang et al., 2004) are not representative of a fast electronic response (Rodríguez-Suné et al., 2021). Similar considerations apply to insulators, semiconductors, and conductive oxides.

Material characterization usually begins with the deposition of relatively thin layers of a given material in order to ascertain its linear, frequency-dependent complex dielectric response (Palik, 1985). The dielectric function so obtained amounts to a value averaged over the thickness of that layer, and cannot account for spatial inhomogeneities due to density fluctuations or surface roughness that may occur during the fabrication process, whose effects are already folded into the measurement and for practical purposes assumed to be negligible. Once the average value of the dielectric constant is obtained, it is used to extract effective particle density, electron mass, plasma frequency, and damping coefficients. The material’s temporal response to external excitation may then be recast in the context of the hydrodynamic model to simulate arbitrary shapes of that material. Thin layers of ITO or Cadmium Oxide (CdO) may display a robust nonlocal response (De Ceglia et al., 2018), so that the measured dielectric function may depend on both frequency and wave-vector, which translates into dependence on spatial derivatives of the fields and boundary conditions. In essence, under the action of radiation pressure, the free electron gas undergoes longitudinal oscillations that trigger standing waves and a resonant spectral response that depends on material thickness in one dimension (De Ceglia et al., 2018), and more complex geometrical and topological considerations in higher dimensions. As another example, the method may be used to predict that in a flat ITO layer, at oblique incidence, nonlocal effects trigger an anisotropic, linear dielectric response in the spectral region where the magnitude of the real part of the dielectric “constant” approaches zero (Scalora et al., 2020). This near-singular behavior is mitigated by a necessary degree of absorption dictated by causality, but that can nevertheless induce a seemingly large nonlinear response, whose quantitative and qualitative aspects (i.e., dispersion) have not yet been clearly understood or quantified (Alam et al., 2016), or appropriately compared to other materials. The same may be said of the optical response of most semiconductors in their respective opacity ranges (Burns and Bloembergen, 1971). Therefore, it is natural to regress and pose questions about the nature and magnitude of nonlinear dispersion in any material, and to identify the intrinsic properties before any hypothetical, metasurface-induced enhancement may be quantified. The method allows us to easily extract bulk nonlinear dispersions from the equations of motion, like the third-order nonlinear susceptibilities, based solely on a comparison with experimental THG conversion efficiencies. Extraction of the bulk, second order response is straightforward and is done in similar fashion. However, for the moment we will refrain from extracting effective second order surface and volume nonlinearities that are triggered by quadrupole-like sources and the magnetic Lorentz contributions, since they cannot be mapped easily to analytical solutions of the nonlinear oscillator model.

Model: Retrieval of Complex Linear and Nonlinear Dispersions From the Equations of Motion

Our first example consists of an investigation of the dielectric constant of undoped crystalline silicon (Palik, 1985) grown in the <100> direction. It is plotted in Figure 1A (markers) and displays two resonances in the UV range. Also in Figure 1A, the two resonances are fitted using two separate Lorentzian functions (solid curves), appropriately detuned and with different plasma frequencies, densities, effective masses and damping coefficients. Therefore, in this case the dynamical model includes two separate bound electron species having different effective parameters. For further details, we refer the reader to reference Rodríguez-Suné et al. (2019), where SHG and THG are discussed for a gallium arsenide wafer, and reference (Scalora et al., 2019), where SHG and THG are studied theoretically for the specific case of a silicon-based, nanowire grating. For completeness, here we limit ourselves to reproducing the basic equations of motion in Gaussian units, and to briefly discuss their content:

FIGURE 1. (A) Real (circle markers) and imaginary (square markers) parts of the complex dielectric constant of silicon as reported in reference Palik (1985). The solid curves are Lorentzian fits to the data as follows: $\mathit{\epsilon }\left(\tilde{\mathit{\omega }}\right)=1-\frac{{\tilde{\omega }}_{p1}^2}{\left({\tilde{\omega }}^2-{\tilde{\omega }}_{01}^2+\mathit{i}{\tilde{\mathit{\gamma }}}_{\mathit{b}1}\tilde{\omega }\right)}-\frac{{\tilde{\omega }}_{p2}^2}{\left({\tilde{\omega }}^2-{\tilde{\omega }}_{02}^2+\mathit{i}{\tilde{\mathit{\gamma }}}_{\mathit{b}2}\tilde{\omega }\right)}$, with $\left({\tilde{\mathit{\omega }}}_{\mathit{p}1},{\tilde{\mathit{\omega }}}_{01},{\tilde{\mathit{\gamma }}}_{\mathit{b}1}\right)=\left(\mathrm{3,2.75},0.1\right)$ and $\left({\tilde{\mathit{\omega }}}_{\mathit{p}2},{\tilde{\mathit{\omega }}}_{02},{\tilde{\mathit{\gamma }}}_{\mathit{b}2}\right)=\left(\mathrm{11,3.3},0.75\right)$. The scaled frequency $\tilde{\mathit{\omega }}$ is in units of 1/microns. (B) Once we have analytical functions that approximate well the measured dispersion (solid curves), we insert the parameters in Eq. 1 to predict the dispersion (markers) calculated as $<\mathit{\epsilon }{\left(\mathit{\lambda }\right)}_{\mathit{Linear}}>=1+4\mathit{\pi }\frac{<\mathit{P}{\left(\mathit{\lambda }\right)}_{\mathit{Linear}}>}{<\mathit{E}{\left(\mathit{\lambda }\right)}_{\mathit{Linear}}>}$ using incident pulses a few tens of femtoseconds in duration.

${\mathit{P}}_{bj}$ is the bound polarization, j = 1,2 represent two separate atomic species. Eq. 1 describe the behavior of electrons that are not allowed to leave atomic sites. The spatial coordinates and time have been scaled with respect to a convenient reference wavelength, ${\lambda }_0=1\mu m$, so that $\tilde{z}=z/{\lambda }_0$, $\tilde{x}=x/{\lambda }_0$, $\tilde{y}=y/{\lambda }_0$ (1 longitudinal and 2 transverse coordinates, respectively, as shown in Figure 2) and time$, \tau =ct/{\lambda }_0$. Fields and polarizations are assumed to be invariant along the transverse $\tilde{x}$ coordinate. For undoped, centrosymmetric silicon, free charges play no role (densities of order 10¹⁴ cm⁻³ or less) and the second order bulk nonlinearity is absent. We assume that each molecular species exhibits a third order nonlinearity expressed by ${\mathit{P}}_{NL}=-\tilde{\beta }\left({\mathit{P}}_{bj}\cdot {\mathit{P}}_{bj}\right){\mathit{P}}_{bj}$. The parameter $\tilde{\beta }\approx \left({\omega }_{0,b1}^2+{\omega }_{0,b2}^2\right){\lambda }_0^2/\left(2L^2{n}_{0b}^2{e}^2{c}^2\right)$ is a scaled, unique real coefficient derived from a nonlinear oscillator model. ${\tilde{\omega }}_{0,b1}$ and ${\tilde{\omega }}_{0,b2}$ are the two resonance frequencies, $n_{0b}\sim {10}^{22}c{m}^{-3}$ is the bound electron density, $m_{bj}^{\text{∗}}$ is the bound electron’s effective mass, c is the speed of light in vacuum, L is the length of a taut, classical spring that corresponds to the lattice constant, ${\tilde{\gamma }}_{bj}$ is a phenomenological damping coefficient, and $\pi {\tilde{\omega }}_{pj}^2=\frac{n_{0,b}{e}^2{\lambda }_0^2}{m_{bj}^{\text{*}}{c}^2}$ is the scaled plasma frequency. The two resonances are located near 300 and 360 nm, respectively. This parameter completely determines nonlinear dispersion of bound electrons, including self-phase modulation, nonlinear absorption, THG conversion efficiencies, etc… It replaces the usual ${\chi }_{\omega }^{\left(3\right)}$ and ${\chi }_{3\omega }^{\left(3\right)}$, which are in fact proportional to $\tilde{\beta }$ (see below) and are often introduced as dispersionless parameters for qualitative and quantitative comparisons with experimental observations. For solids, L can vary from a fraction of 1 Å to a few Å, a disparity that is reflected on the particle density $n_{0b}$. For example, the wave function of valence (bound) electrons in silicon peaks near 1.32 Å (https://www.infoplease.com/semimetallics/silicon), suggesting that $L_{silicon}\sim 2.6\times {10}^{-8}cm$. However, the highest d-orbital wave function for bound electrons in gold peaks near 0.5 Å (Mann, 1968; Kaupp, 2007), with $L_{gold}\sim {10}^{-8}cm$. One should then be mindful that the magnitude of $\tilde{\beta }$ can change considerably from material to material, and can range between 10⁻⁶ and 10⁻⁸, depending on atomic orbital radii, densities, and resonance frequencies. In Eq. 1 the term $\frac{e{\lambda }_0}{m_{bj}^{\text{*}}{c}^2}\left({\mathit{P}}_{bj}\cdot \nabla \right)\mathit{E}$ is a surface nonlinearity, and the magnetic Lorentz contribution $\frac{e{\lambda }_0}{m_{bj}^{\text{*}}{c}^2}{\dot{\mathit{P}}}_{bj}\times \mathit{H}$ contains both surface and volume nonlinear bound currents (Rodríguez-Suné et al., 2019). We expand these terms and account for pump depletion and down-conversion to occur. These terms are responsible for the generated second harmonic signal in a centrosymmetric system, as well as a TH signal much weaker compared to the TH generated by the term ${\mathit{P}}_{NL}=-\tilde{\beta }\left({\mathit{P}}_{bj}\cdot {\mathit{P}}_{bj}\right){\mathit{P}}_{bj}$.

FIGURE 2. Set-up developed to measure SH and TH signals.

Since the model is described in considerable details elsewhere, here we outline a new technique that we refer to as numerical ellipsometry, that allows us to first verify the measured linear material dispersion, and then predict complex nonlinear dispersion functions (i.e., ${\chi }_{\omega }^{\left(3\right)}$ and ${\chi }_{3\omega }^{\left(3\right)}$) using our dynamical model, which includes Eq. 1, Maxwell’s equations, the macroscopic constitutive relations, and experimentally obtained harmonic conversion efficiencies. In simplified terminology, the method to retrieve nonlinear dispersion consists of taking the following steps: first, we perform a calculation in the linear regime (low power densities) using a pulse only a few tens of femtoseconds in duration, incident normal to the surface, and extract spatially averaged, complex polarizations and fields inside a 20 nm-thick layer of material when the peak of the pulse reaches the layer. We thus estimate a spatially averaged, complex dielectric constant as $\langle \epsilon {\left(\lambda \right)}_{Linear}\rangle =1+4\pi \frac{\langle P{\left(\lambda \right)}_{Linear}\rangle }{\langle E{\left(\lambda \right)}_{Linear}\rangle }$, at any given carrier wavelength λ, where $P{\left(\lambda \right)}_{Linear}$ and $E{\left(\lambda \right)}_{Linear}$ are complex polarization and field amplitudes, respectively. The brackets indicate spatial averages. This straightforward procedure allows us to reconstruct, and verify, the linear, monochromatic material dispersion data quite well even using very short pulses, because sharp spectral features are absent, and absorption resonances are quite broad. We report the results of this calculation in Figure 1B as empty circle and square markers, which are now seen to map quite well and are nearly indistinguishable from the solid, analytical Lorentzian curves that have replaced the complex dielectric data. Other semiconductors display similar behavior in the 100 nm range, where the magnitude of the complex dielectric constant approaches zero.

The second important step consists of performing the same calculations in the nonlinear regime, to once again estimate spatially-averaged polarization and field inside the same layer of material, such that $\langle \epsilon {\left(\lambda \right)}_{Nonlinear}\rangle =1+4\pi \frac{\langle P{\left(\lambda \right)}_{Nonlinear}\rangle }{\langle E{\left(\lambda \right)}_{Nonlinear}\rangle }$. The full nonlinear dispersion of the material is generally obtained in the wavelength range of interest, by taking the difference between linear and nonlinear, spatially-averaged polarizations, at both pump and harmonic wavelengths. For example, if we wish to retrieve ${\chi }_{\omega }^{\left(3\right)}$, in the undepleted pump approximation we may write:

For simplicity, spatial averages are now implied and bracket symbols dropped. Additional considerations and terms are needed on the right side of Eq. 2 if either the pump is allowed to deplete, or if surface, magnetic, and/or higher order nonlinearities come into play. In Eq. 2, a linear calculation is required to obtain ${\chi }_{\omega }^L=\frac{P_{\omega }^L}{E_{\omega }^L}$, where $P_{\omega }^L$ and $E_{\omega }^L$ are the linear, spatially-averaged polarization and field at the fundamental frequency. It is important to distinguish between the fields that are calculated in the linear regime, denoted by L, and the fields calculated in the nonlinear regime, denoted by NL. Eq. 2 may be inverted to yield:

Once again, we emphasize that all fields and polarizations are evaluated and spatially-averaged inside the layer. As such, the numerical procedure that we employ is identical to the ellipsometric extraction of the linear dielectric constant, since it too produces average values of the dielectric constant.

In order to retrieve ${\chi }_{3\omega }^{\left(3\right)}$, the procedure mirrors that used to recover ${\chi }_{\omega }^{\left(3\right)}$. We first write the expression for the third order polarization as:

As was the case for Eq. 2, here too modifications are required if additional nonlinearities become relevant. As in Eq. 3, a linear calculation is required to obtain ${\chi }_{3\omega }^L=\frac{P_{3\omega }^L}{E_{3\omega }^L}$, where $P_{3\omega }^L$ and $E_{3\omega }^L$ are the linear, spatially-averaged, complex polarization and field at the third harmonic wavelength, respectively. In contrast, $E_{3\omega }^{NL}$ is the field generated at the third harmonic wavelength when pumping at the fundamental frequency, which is necessarily nonlinear, while $E_{\omega }^{NL}$ is now the nonlinear pump field. Therefore, we may write:

As we will see below for ITO, at oblique incidence steps must be taken in order to account for the vector nature of the fields. Similarly, if the system displays dielectric anisotropies due to varying effective masses, densities or spring constants in different spatial directions, then Eqs 1–5 should be modified accordingly.

In the case of silicon and other similarly undoped semiconductors that display multiple absorption resonances (for example, GaAs, GaP, Ge, etc…) it is possible to derive analytical expressions for ${\chi }_{\omega }^{\left(3\right)}$ and ${\chi }_{3\omega }^{\left(3\right)}$ from Eq. 1 using the same perturbative approach based on Miller’s rule (Miller, 1964), illustrated in references Boyd (2003), Scalora et al. (2015) for a single oscillator, provided the pump remains undepleted, and no other nonlinearities or effects enter the picture. Under those conditions, the scaled expressions for third order nonlinear susceptibilities for a two-resonance system may be derived and written as follows:

where $\tilde{\omega }=\frac{\omega }{{\omega }_r}$ is a scaled frequency, ${\omega }_r=\frac{2\pi c}{{\lambda }_r}$, and ${\tilde{\omega }}_p^2=\frac{n{e}^2}{m_b^{\ast }}{\left(\frac{{\lambda }_r}{c}\right)}^2$ is the scaled plasma frequency. In MKS units ${\chi }_{\omega ,3\omega ,MKS}^{\left(3\right)}=\frac{4\pi }{{\left(3\times {10}^4\right)}^2}{\chi }_{\omega ,3\omega }^{\left(3\right)}$.

Experimental Set-up

Our experimental set-up is shown schematically in Figure 2. It measures second and third harmonic generated signals as a function of the pump’s angle of incidence and polarization, in both transmission and reflection configurations. For the silicon sample that we investigated, the source is a pulsed fiber laser (FYLA PS50) that emits a train of 13 ps pulses at 1064 nm, with a CW average output power of 2 W and 1 MHz repetition rate, delivering 2 μJ/pulse. This source can be modulated to deliver a train of N pulses at a frequency repetition rate of 1 kHz, thus enabling the detection of weaker SH or TH signals by integrating the response of N pulses on the photomultiplier. A half-wave plate controls the polarization of the fundamental field (FF), enabling illumination of the sample with either TM- or TE-polarized light. Spectral filters are used to remove SH and/or TH signals arising from different optical components placed before the sample. The beam is focused on the sample plane to obtain fundamental beam peak power densities on the order of a few GW/cm². After the FF traverses the crystalline silicon sample, generating SH and TH fields, it is attenuated using a filter to avoid any potentially significant SH and/or TH generation from the surfaces of the optical elements placed after the sample’s position. Another lens collimates the beam and a polarizer analyzes the polarization (TM or TE) of the generated SH and TH fields. A dispersive prism is used in combination with a blocking edge to separate and obscure the remaining FF radiation from the SH and TH paths, so that the only radiation arriving at the detector is either the SH or TH signal arising from the silicon wafer. The detector consists of a photomultiplier tube (Hamamatsu H10722 PMT series) with a mounted narrow-band spectral filter having a 20 nm band pass transmission around either the SH or the TH wavelength. The entire detection system is mounted on a rotary platform that allows measurements in transmission and in reflection. The sample is positioned on a motorized rotation stage so that the SH or TH signals can be measured as a function of the angle of incidence. We placed a BBO crystal at the sample position generating a SH signal that could be measured with a calibrated photodiode in order to estimate conversion efficiencies, as the ratio between the generated SH or TH energy (transmitted or reflected), and the total initial fundamental pulse energy. By measuring the generated NL signal at the sample plane and at the PM position (after traversing the total path across the detection system) we can estimate the losses in the experimental system for each polarization state. By replacing the photodiode with the photomultiplier, after attenuating the signal using neutral density filters, we obtain the relation between the measured signal at the PM and the energy of the generated NL signal at the sample plane. This experimental setup allows the recording of harmonic conversion efficiencies as low as 10⁻¹⁰.

Results

Silicon

Using our set-up we record transmitted and reflected THG conversion efficiencies as functions of the angle of incidence, as the fundamental pump pulse tuned to 1064 nm traverses a silicon wafer 500 microns thick. A focusing lens having focal length f = 100 mm was used to obtain fundamental beam intensities of approximately 6 GW/cm² on the silicon wafer. A TM-polarized TH signal was detected for TM-polarized incident light, shown in Figure 3 plotted in blue, while the TE-TE case is depicted in red. We compared the experimental results (empty circles and squares) with numerical simulations (dashed curves) carried out via integration of Eq. 1 together with Maxwell’s equations. The dashed curves are the corresponding predicted signals. In the absence of slower nonlinearities or other spurious nonlinear effects, the results of the simulations depicted in Figure 3 emerge for pulses only a few tens of femtoseconds in duration, and remain insensitive to pulse duration. Estimates of measured conversion efficiencies for surface and magnetically induced SHG are of order 10⁻¹², and thus inconsequential to the retrieval of third order nonlinear dispersion curves. If the sample is thicker than a few tens of microns, reflected harmonic generation becomes independent of sample thickness, because the pump is absorbed faster than the round trip time necessary to trigger meaningful cavity effects. Figure 3 suggests reasonably good agreement between predictions and measured efficiencies, which are of order 10⁻⁷, obtained by inserting $\tilde{\beta }\approx 3.6\times {10}^{-7}$ in Eq. 1. A transmitted TH signal was also detected at 354 nm for both polarization configurations, with efficiencies also of order 10⁻⁷. The amplitude of the transmitted signal depends on sample thickness because the pump is absorbed as it decays inside the sample, a fact confirmed by our simulations. That notwithstanding, the shape of the transmitted signal as a function of incident angle is not altered regardless of thickness. Our numerical results also confirm this. The TH signal is tuned deep into the absorption range, but it is nevertheless transmitted. The reason that a TH signal tuned in the absorption range emerges from a silicon sample half a millimeter in thickness is due to the generation of an inhomogeneous component at the harmonic wavelength, which propagates with the pump’s dispersion (Burns and Bloembergen, 1971). Since the sample is semi-transparent to the pump, the sample will also be semi-transparent at the harmonic wavelength. This so-called phase-locked component has been previously observed and discussed in details in references Rodríguez-Suné et al. (2019), Scalora et al. (2019), Scalora et al. (2015), Roppo et al. (2007), Centini et al. (2008), Fazio et al. (2009), Roppo et al. (2009), Roppo et al. (2011a), Roppo et al. (2011b).

FIGURE 3. Reflected TM-polarized TH for a TM-polarized FF (blue), TE-polarized TH for a TE-polarized FF (red). Experimental results are plotted with empty circles and squares; numerical simulations are depicted with dashed lines.

In Figure 4 we compare ${\chi }_{\omega ,MKS}^{\left(3\right)}$ calculated using Eqs 3, 6, in the visible and near-IR ranges. The inset shows contrast extended well into the UV range. Both figure and inset show remarkably good agreement between analytical results derived using the nonlinear oscillator model, and the simulations performed using ultrashort pulses. We note that while in the near-IR range the magnitude of ${\chi }_{\omega ,MKS}^{\left(3\right)}$ is of order 10⁻¹⁷ (m/V)² and decreasing at longer wavelengths, near resonance at UV wavelengths its amplitude is catapulted upward by nearly three orders of magnitudes. The comparison between Eqs 5, 7 for ${\chi }_{3\omega ,MKS}^{\left(3\right)}$ is shown in Figure 5, with maximum values occurring near resonance. The agreement between analytical and retrieved values is noteworthy. The nonlinear dispersion exhibited by silicon in Figures 4, 5 is consistent with the THG efficiencies reported above. A check may be performed by inserting the retrieved complex nonlinear dispersion into an independent plane wave model such as the one provided by the frequency-domain, finite-element method (COMSOL Multiphysics). The conversion efficiencies are reproduced just as reported in Figure 3.

FIGURE 4. Real and imaginary parts of the complex ${\mathit{\chi }}_{\mathit{\omega },\mathit{MKS}}^{\left(3\right)}$ calculated analytically viaEq. 6 (solid blue and red curves), and calculated viaEq. 3 (empty blue circles and red square markers).

FIGURE 5. Real and imaginary parts of the complex ${\mathit{\chi }}_{3\mathit{\omega },\mathit{MKS}}^{\left(3\right)}$ calculated analytically viaEq. 7 (solid blue and red curves), and calculated viaEq. 5 (empty blue circles and red square markers).

Gold

Using the same approach outlined above, we are now able to compare analytical results for gold with our simulations, which in turn are based on previously reported experimental results of THG in nanometer-thick gold layers, probed with both picosecond and femtosecond pulses (Rodríguez-Suné et al., 2021). The situation is different for gold compared to undoped silicon, because now a free electron gas and nonlocal effects may contribute to the dielectric constant. This difference is fully delineated in reference Rodríguez-Suné et al. (2021), where the model is extended to include one free and two bound electron species to the dielectric constant. Accordingly, Eq. 6 is modified as follows:

with a similar alteration to Eq. 7. The coefficient $\widetilde{\text{Λ}}$ relates to the nonlinear third order contributions of free (hot) electrons, whose Fermi surface is modified as a result of increased free electron density as a function of applied intensity (Rodríguez-Suné et al., 2021).

In Figure 6A we show the linear gold dielectric constant data found in reference Palik (1985), fitted using one Drude and two Lorentzian functions, along with the data retrieved using our model. The calculations are carried out at normal incidence on a 20 nm-thick gold layer, which was investigated in reference Rodríguez-Suné et al. (2021) in the picosecond and femtosecond regimes. Once again, the agreement between the simulations and the observed data is notable. We remark that for wavelengths longer than 1 μm some disagreements begin to emerge in the Im(ε), that are likely due to slower computational convergence as a result of the large values of the dielectric constant. The imaginary parts of the nonlinear dispersions depicted in Figures 6B,C also display a similar slight divergence from the analytical results. However, the overall agreement between analytical and numerical results is once again quite remarkable. Nonlocal effects play a minor role for these gold layer thicknesses. Once again, we note that the retrieved nonlinear dispersion of gold in Figure 6 is consistent with THG efficiencies reported in reference Rodríguez-Suné et al. (2021).

FIGURE 6. (A) Real and imaginary parts of the complex dielectric constant of gold found in reference Palik (1985) (Solid curves), and as retrieved using our model (markers). (B) Analytical (solid curves) and retrieved (markers) ${\mathit{\chi }}_{\mathit{\omega },\mathit{MKS}}^{\left(3\right)}$; (C) Analytical (solid curves) and retrieved (markers) ${\mathit{\chi }}_{3\mathit{\omega },\mathit{MKS}}^{\left(3\right)}$. In this case, $\tilde{\beta }$ = 1.2 × 10⁻⁶, with $\widetilde{\text{Λ}}$ = 3.7 × 10⁻⁷,$\mathit{\epsilon }\left(\tilde{\mathit{\omega }}\right)=1-\frac{{\tilde{\mathit{\omega }}}_{\mathit{p},\mathit{free}}^2}{\left({\tilde{\mathit{\omega }}}^2+\mathit{i}{\tilde{\mathit{\gamma }}}_{\mathit{free}}\tilde{\mathit{\omega }}\right)}-\frac{{\tilde{\mathit{\omega }}}_{\mathit{p}1}^2}{\left({\tilde{\mathit{\omega }}}^2-{\tilde{\mathit{\omega }}}_{01}^2+\mathit{i}{\tilde{\mathit{\gamma }}}_{\mathit{b}1}\tilde{\mathit{\omega }}\right)}-\frac{{\tilde{\mathit{\omega }}}_{\mathit{p}2}^2}{\left({\tilde{\mathit{\omega }}}^2-{\tilde{\mathit{\omega }}}_{02}^2+\mathit{i}{\tilde{\mathit{\gamma }}}_{\mathit{b}2}\tilde{\mathit{\omega }}\right)}$, with $\left({\tilde{\mathit{\omega }}}_{\mathit{pfree}},{\tilde{\mathit{\gamma }}}_{\mathit{free}}\right)=\left(6.45,0.05\right)$, $\left({\tilde{\mathit{\omega }}}_{\mathit{p}1},{\tilde{\mathit{\omega }}}_{01},{\tilde{\mathit{\gamma }}}_{\mathit{b}1}\right)=\left(\mathrm{3,2.75},0.1\right)$ and $\left({\tilde{\mathit{\omega }}}_{\mathit{p}2},{\tilde{\mathit{\omega }}}_{02},{\tilde{\mathit{\gamma }}}_{\mathit{b}2}\right)=\left(\mathrm{11,3.3},0.75\right)$.

Indium Tin Oxide

The measured local dielectric constant of ITO is displayed in Figure 7A for a 20 nm-thick layer, along with a fit that includes one Drude (free electrons) and one Lorentz oscillator (bound electrons). The absorption resonance becomes discernible in the data near 300 nm. The presence of the Lorentz resonance ascribes an intrinsic, nonlinear third order response described by Eqs 6, 7, usually neglected in typical theoretical treatments, but that supplements and competes with the hot electron nonlinearity (Scalora et al., 2020; Rodríguez-Suné et al., 2020), which in turn accounts for a dynamic (time dependent) redshift of the free electron plasma frequency. The relevance of nonlocal effects in the electrodynamics of conducting oxides becomes conspicuous in Figures 7B,C, where we plot the retrieved transverse and longitudinal dielectric constants using our approach, for an incident angle of 65° (Scalora et al., 2020; Rodríguez-Suné et al., 2020). Additional absorption resonances can be seen to form at wavelengths that correspond to the standing wave conditions of the free electron gas component oscillating in the longitudinal direction, which in turn affect nonlinear dispersion. No analytical solutions equivalent to Eq. 8 are known in the presence of nonlocal effects. At oblique incidence, an expansion of the term $\tilde{\beta }\left({\mathit{P}}_{bj}\cdot {\mathit{P}}_{bj}\right){\mathit{P}}_{bj}$ must be carried out into its vector and frequency components. Assuming the simultaneous presence of the fundamental and its generated harmonics, the basic requirement is to write the total polarization as:

FIGURE 7. (A) Real and imaginary parts of the complex, local dielectric constant of a 100 nm-thick layer of ITO (markers). The solid curves are a Drude-Lorentz fit of the data, which suggest the presence of an absorption resonance in the 200–300 nm range. Retrieved transverse (B) and longitudinal (C) dielectric constants obtained at 65° angle of incidence. Nonlocal effects trigger anisotropic behavior, as outlined in Palik (1985). Nonlocal resonances are visible in both (B) and (C), and are more pronounced in the longitudinal directions. For ITO, $\tilde{\beta }$ = 10^-7, with $\widetilde{\text{Λ}}$ = −1 × 10⁻⁸,$\mathit{\epsilon }\left(\tilde{\mathit{\omega }}\right)=1-\frac{{\tilde{\mathit{\omega }}}_{\mathit{p},\mathit{free}}^2}{\left({\tilde{\mathit{\omega }}}^2+\mathit{i}{\tilde{\mathit{\gamma }}}_{\mathit{free}}\tilde{\mathit{\omega }}\right)}-\frac{{\tilde{\mathit{\omega }}}_{\mathit{p}1}^2}{\left({\tilde{\mathit{\omega }}}^2-{\tilde{\mathit{\omega }}}_{01}^2+\mathit{i}{\tilde{\mathit{\gamma }}}_{\mathit{b}1}\tilde{\mathit{\omega }}\right)}$ with $\left({\tilde{\mathit{\omega }}}_{\mathit{pfree}},{\tilde{\mathit{\gamma }}}_{\mathit{free}}\right)=\left(1.5785,0.1105\right)$ and $\left({\tilde{\mathit{\omega }}}_{\mathit{p}1},{\tilde{\mathit{\omega }}}_{01},{\tilde{\mathit{\gamma }}}_{\mathit{b}1}\right)=\left(11.5,\mathrm{7,0.01}\right).$

We have assumed that fields and polarizations do not vary along the transverse x coordinate. Then, following straightforward manipulation of the vector product, development of a cubic equation, and use of the constitutive relation $\mathbf{P}=\chi \mathbf{E}$, Eqs 3, 5 take the following form:

and

At normal incidence Eqs 10, 11 yield Eqs 3, 5. The curves that correspond to Eq. 10 are plotted in Figure 8, while Eq. 11 are displayed in Figure 9. The strong oscillatory behavior that characterizes especially the longitudinal nonlinear response in Figure 8A is triggered by the nonlocal resonances seen in the linear response of Figure 7. The THG data used to retrieve the nonlinear dispersion depicted in Figures 8, 9 was obtained using femtosecond pulses, and may be found in reference Rodríguez-Suné et al. (2020). As we tune the pump away from the ENZ condition, nonlocal effects are attenuated, and the nonlinear response becomes more Lorentz-like. Our previous report on harmonic generation from an ITO nanolayer using femtosecond pulses suggests that at the ENZ condition THG conversion efficiencies are of similar order of magnitude compared to THG efficiencies presently recorded for the silicon wafer. A comparison between the nonlinear dispersions that we have derived for the various materials suggests that silicon may have the largest intrinsic nonlinear third order response in the visible and UV ranges.

FIGURE 8. Real and imaginary parts of the (A) longitudinal and (B) transverse third order susceptibility experienced by the pump as a function of incident wavelength. The oscillatory behavior reflects the oscillations that characterize the linear dielectric response shown in Figure 7. Just as in Figure 7, the oscillations occur around the local nonlinear response obtained viaEq. 8 (not shown for clarity).

FIGURE 9. Real and imaginary parts of the longitudinal (solid) and transverse (dashed) third order susceptibility experienced by the TH signal, as a function of incident wavelength. The wavelength range is question is far from the ENZ condition, and thus insusceptible to nonlocal effects.

Summary

We have presented a combined experimental and numerical method that can be used to predict complex nonlinear dispersion curves in almost any material, based exclusively on the experimental determination of harmonic generation conversion efficiencies. We used a hydrodynamic approach that faithfully duplicates linear dispersion, and predicted the wavelength dependent nonlinear response. The method is particularly useful when analytical solutions are not available, as is the case for conductive oxides like ITO, which displays nonlocal effects that trigger an effective anisotropy, and extendable to second order bulk nonlinearities. Our experimental results for THG in silicon suggest that it is possible to exploit its large nonlinear response in the visible and UV ranges, thus opening up new prospects for silicon photonics (Scalora et al., 2019).

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

All authors contributed to the conception and design of this study. MS, DC, and MV performed the numerical simulations. LR-S, JT, NA, and CC carried out the experiments. All authors contributed equally to this work.

Funding

RDECOM Grant W911NF-18-1-0126 from the International Technology Center-Atlantic, and from the Spanish Agencia Estatal de Investigación (project no. PID2019-105089GB-I00/AEI/10.13039/501100011033). Army Research Laboratory Cooperative Agreement Numbers W911NF1920279 and W911NF2020078 issued by US ARMY ACC-APG-RTP.

Conflict of Interest

Author NA was employed by the company Bluehalo Inc.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Z. Coppens, K. Hallman, and D. Dement for help assessing and measuring the optical properties of the various samples, and for helpful discussions. We thank M. Centini, M. Sanghadasa, and N. Litchinitser for critical reading of the manuscript.

References