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 ϕ N L 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 ( ϕ N L < 0.2 π ), 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 χ ω ( 3 ) 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 χ 3 ω ( 3 ) 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.

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: P b j ¨ + γ ˜ b j P ˙ b j + ω ˜ 0 , b j 2 P b j − β ˜ ( P b j ⋅ P b j ) P b j = π ω ˜ p j 2 E + e λ 0 m b j ∗ c 2 ( P b j ⋅ ∇ ) E + e λ 0 m b j ∗ c 2 P ˙ b j × H (1)

P b j 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, λ 0 = 1 μ m , so that z ˜ = z / λ 0 , x ˜ = x / λ 0 , y ˜ = y / λ 0 (1 longitudinal and 2 transverse coordinates, respectively, as shown in Figure 2) and time ,   τ = c t / λ 0 . Fields and polarizations are assumed to be invariant along the transverse 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 P N L = − β ˜ ( P b j ⋅ P b j ) P b j . The parameter β ˜ ≈ ( ω 0 , b 1 2 + ω 0 , b 2 2 ) λ 0 2 / ( 2 L 2 n 0 b 2 e 2 c 2 ) is a scaled, unique real coefficient derived from a nonlinear oscillator model. ω ˜ 0 , b 1 and ω ˜ 0 , b 2 are the two resonance frequencies, n 0 b ∼ 10 22 c m − 3 is the bound electron density, m b j ∗ 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, γ ˜ b j is a phenomenological damping coefficient, and π ω ˜ p j 2 = n 0 , b e 2 λ 0 2 m b j * 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 χ ω ( 3 ) and χ 3 ω ( 3 ) , which are in fact proportional to β ˜ (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 0 b . For example, the wave function of valence (bound) electrons in silicon peaks near 1.32 Å (https://www.infoplease.com/semimetallics/silicon), suggesting that L s i l i c o n ∼ 2.6 × 10 − 8 c m . However, the highest d-orbital wave function for bound electrons in gold peaks near 0.5 Å (Mann, 1968; Kaupp, 2007), with L g o l d ∼ 10 − 8 c m . One should then be mindful that the magnitude of β ˜ 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 e λ 0 m b j * c 2 ( P b j ⋅ ∇ ) E is a surface nonlinearity, and the magnetic Lorentz contribution e λ 0 m b j * c 2 P ˙ b j × 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 P N L = − β ˜ ( P b j ⋅ P b j ) P b j .

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., χ ω ( 3 ) and χ 3 ω ( 3 ) ) 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 〈 ε ( λ ) L i n e a r 〉 = 1 + 4 π 〈 P ( λ ) L i n e a r 〉 〈 E ( λ ) L i n e a r 〉 , at any given carrier wavelength λ, where P ( λ ) L i n e a r and E ( λ ) L i n e a r 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 〈 ε ( λ ) N o n l i n e a r 〉 = 1 + 4 π 〈 P ( λ ) N o n l i n e a r 〉 〈 E ( λ ) N o n l i n e a r 〉 . 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 χ ω ( 3 ) , in the undepleted pump approximation we may write: P ω N L = χ ω L E ω N L + χ ω ( 3 ) | E ω N L | 2 E ω N L (2)

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 χ ω L = P ω L E ω L , where P ω L and E ω 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: χ ω ( 3 ) = P ω N L − P ω L E ω L E ω N L | E ω N L | 2 E ω N L = ε ω N L − ε ω L 4 π | E ω N L | 2 (3)

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 χ 3 ω ( 3 ) , the procedure mirrors that used to recover χ ω ( 3 ) . We first write the expression for the third order polarization as: P 3 ω N L = χ 3 ω L E 3 ω N L + χ 3 ω ( 3 ) ( E ω N L ) 3 (4)

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 χ 3 ω L = P 3 ω L E 3 ω L , where P 3 ω L and E 3 ω L are the linear, spatially-averaged, complex polarization and field at the third harmonic wavelength, respectively. In contrast, E 3 ω N L is the field generated at the third harmonic wavelength when pumping at the fundamental frequency, which is necessarily nonlinear, while E ω N L is now the nonlinear pump field. Therefore, we may write: χ 3 ω ( 3 ) = P 3 ω N L − P 3 ω L E 3 ω L E 3 ω N L ( E ω N L ) 3 (5)

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 χ ω ( 3 ) and χ 3 ω ( 3 ) 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: χ ω ( 3 ) = 3 β ˜ 4 π 2 ( ω ˜ p 1 2 4 π ) 3 ( ω ˜ 01 2 − ω ˜ 2 − i γ ˜ b 1 ω ˜ ) 3 ( ω ˜ 01 2 − ω ˜ 2 − i γ ˜ b 1 ω ˜ ) + 3 β ˜ 4 π 2 ( ω ˜ p 1 2 4 π ) 3 ( ω ˜ 02 2 − ω ˜ 2 − i γ ˜ b 2 ω ˜ ) 3 ( ω ˜ 02 2 − ω ˜ 2 − i γ ˜ b 2 ω ˜ ) (6) χ 3 ω ( 3 ) = β ˜ 4 π 2 ( ω ˜ p 1 2 4 π ) 3 ( ω ˜ 01 2 − ω ˜ 2 − i γ ˜ b 1 ω ˜ ) 3 ( ω ˜ 01 2 − 9 ω ˜ 2 − 3 i γ ˜ b 1 ω ˜ ) + β ˜ 4 π 2 ( ω ˜ p 1 2 4 π ) 3 ( ω ˜ 02 2 − ω ˜ 2 − i γ ˜ b 2 ω ˜ ) 3 ( ω ˜ 02 2 − 9 ω ˜ 2 − 3 i γ ˜ b 2 ω ˜ ) (7) where ω ˜ = ω ω r is a scaled frequency, ω r = 2 π c λ r , and ω ˜ p 2 = n e 2 m b ∗ ( λ r c ) 2 is the scaled plasma frequency. In MKS units χ ω , 3 ω , M K S ( 3 ) = 4 π ( 3 × 10 4 ) 2 χ ω , 3 ω ( 3 ) .

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⁻¹⁰.

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).