High-harmonic generation from subwavelength silicon films

1 Introduction

The last decade has witnessed a renewed quest for efficient nonlinear frequency conversion implementing micrometric and nanometric size materials and devices that have highlighted the importance of nanostructures and metasurfaces. This pursuit has been marked by a dramatic increase in theoretical and experimental work [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12] that has led to a natural reassessment of the optical properties of semiconductors and metals in optical regimes that were previously deemed impractical or inaccessible [13], [14], [15], [16], [17], [18], [19], [20], [21]. This kind of research has generally been limited to second and third harmonic generation (SHG and THG, respectively) primarily due the applicability of these processes to surface sensing, label-free bio-imaging, and quantum optics [22], [23], [24], [25]. However, nonlinear high harmonic generation (HHG) encompassing visible light, ultraviolet (UV), and even shorter wavelengths, is also of fundamental interest for other key applications such as ultra-fast physics, and X-UV photolithography [10], [26]. Most work on HHG has been carried out in transparent conductive oxides (TCOs) [27], [28], [29], [30] and all dielectric structures [5], [6] for the purpose of avoiding intrinsic losses associated with semiconductors and metals [31]. This has led researchers to not only set metals aside but also to focus almost exclusively on the transparency spectral range of dielectric and semiconductor materials.

Generally, the study of nonlinear optics at the nanoscale requires modifications of the material equations of motion to account for effects that manifest themselves only at the atomic level, given that the average atomic diameter is of order 3 Å–5 Å. Examples are linear and nonlinear nonlocal effects due to pressure and viscosity of the free electron fluid and spill-out effects [32], [33], [34], [35], [36], as well as tunnelling [37], [38], screening and surface and magnetic phenomena [39], [40], [41]. Even in transparent materials, most of the work has focused on bulk nonlinearities, while neglecting crucial surface and magnetic phenomena that in centrosymmetric materials have been traditionally associated with the generation of even-order harmonics [39]. Most semiconductors have broad absorption resonances deep in the UV range, perhaps suggesting that absorption may not be circumvented, with negative dielectric constants and epsilon near zero (ENZ) crossing points in the 100 nm range. However, it has been demonstrated both theoretically and experimentally that SHG and THG can occur despite absorption in the visible and UV ranges and can be transmitted through half millimeter thick wafers of silicon (Si), gallium arsenide (GaAs), and gallium phosphide (GaP) with the pump tuned in the transparency regime, via a phenomenon referred to as phase-locking [13], [14], [15], [16]. Although for nanometer-thick materials the conversion process is inefficient due to the absence and loss of meaning of phase matching, efficiency can be improved dramatically if a resonant mechanism comes into play as demonstrated for semiconductor cavities [17], and at least in theory in GaAs-filled metal gratings [18] and silicon nanowire arrays [19]. This phenomenon plays a critical role in nonlinear optical phenomena for harmonic generation in all semiconductors in the opacity range, and so we now briefly delve into the fundamental reasons that make this possible.

As originally pointed out by Bloembergen and Pershan [42] for SHG in transparent media, Maxwell’s equations in a nonlinear material have two solutions: a homogeneous component that propagates with a wave vector k 2 ω hom = n 2 ω 2 ω / c , and an inhomogeneous component having wave vector k 2 ω i n h = 2 k ω = 2 n ω ω / c . These distinct solutions can be generalized to the jth harmonic. In the plane wave regime and oblique incidence, these waves refract at different angles, interfere, and give rise to Maker fringes [43], with the inhomogeneous portion propagating at the same refraction angle and velocity as the pump beam, a phenomenon more recently clearly observed in bulk lithium niobate [44]. However, if the interaction occurs in such a way that the pump is tuned in the transparency range [ n ω is real], and the harmonic field is tuned in the opaque region [ n 2 ω is complex], then the homogenous component is absorbed, leaving behind the inhomogeneous component, setting up an anomalous situation, i.e., a harmonic signal that propagates with the same dispersion and group velocity of the pump, a phenomenon that has been referred to as phase locking, given the phase relationships between pump and the relevant inhomogeneous harmonic solutions [13], [14], [15], [16], [17], [18], [19], [20], [44].

Thanks to phase locking, a harmonic field tuned in the opaque regime can resonate along with the pump, as if absorption were not present for that component, provided the pump is tuned in the transparency range. This is the fundamental mechanism that explains the apparent suppression of absorption observed, for example in references [8], [9], and most recently highlighted and reported in the opacity range of a chalcogenide glass grating at 354 nm [45], well inside the band gap, and in silicon etalons [46], with THG conversion efficiencies of order 10⁻⁷ at 266 nm, with a pump intensity of 1.5 GW/cm². The remarkable facts here are that one is exploiting the intrinsic properties of silicon without local field amplification [see Supplementary Information], and that silicon is metallic [Re(ε) < 0] between 100 nm and 300 nm. This efficiency should be contrasted with purely experimental observations of THG from three-dimensional (3-D), resonant silicon metasurfaces characterized by guided mode resonances and bound states in the continuum [47], [48]; where similar peak power densities yield THG conversion efficiencies of order 10⁻⁴. In the present case, using the same theory illustrated in references [13], [19], [46], we analyze experimental results that yield high harmonic generation in silicon using a time-domain, hydrodynamic-Maxwell pulse propagation method that considers linear and nonlinear bulk material dispersions, nonlinearities triggered by the mere traversal of surfaces, and by the intrinsically nonlinear magnetic portion of the Lorentz force. As we will see below, the confluence of the above-listed physical phenomena enables an accurate theoretical description of the relevant qualitative and quantitative aspects of nonlinear frequency conversion in the perturbative regime. We note that it is possible to extend the method to the non-perturbative regime, which we define as a harmonic oscillator far from its ground state, requiring adopting additional polarization and inversion components generated by a system of Bloch equations that describe a two- or multi-level atom, as shown in references [49], [50], where it was applied to stimulated Raman scattering near plasmonic structures.

The paper is organized as follows: in Section 2 we will first introduce equations of motion (Eqs. (1) and (2)) that describe the two absorption (Lorentzian) resonances that characterize silicon in the UV range. In Eq. (3) we then introduce the dynamics of conduction electrons, whose density is proportional to the incident field intensity, and which are mostly responsible for even harmonic generation at large peak power densities via surface (Coulomb, convective terms, and Lorentz magnetic force) and volume (Lorentz magnetic force) nonlinear sources. Unlike other approaches where high harmonic generation is also tackled theoretically [51], [52], our approach includes surface and magnetic nonlinearities, also described in some detail for both free and bound charges, which are essential for the generation of even harmonic fields in centrosymmetric materials. Finally in Section 2 we also describe our experimental measurements, compare with the results of our simulations, and offer a closing summary of our present effort.

2 Theory and results

The suspended sub-wavelength silicon films used in our experiments are pictured and schematically depicted in Figure 1(a) and (b). The resonance near 530 nm gives the sample a characteristic green hue. They were purchased from Norcada (Alberta, Canada) and consist of ∼200 nm-thick <100> silicon film with doping level of order 10¹⁶/cm³, etched out of a silicon carrier wafer, and held by a frame of proprietary composition. This density corresponds to a plasma frequency in the THz range and does not affect the dynamics at low intensities. We will return to the subject later under peak intensity conditions. In Figure 1(c) we report both measured and fitted transmittance at normal incidence. The etalon displays a series of Fabry–Perot resonances whose amplitude expectedly decreases rapidly with decreasing wavelength. Our measurements are well represented by Palik’s data [31], which we show in Figure 1(d) along with a fit obtained using two detuned Lorentzian functions that we will use in our simulations. The sample is completely opaque below 400 nm, erroneously suggesting that we should not expect any meaningful nonlinear response in this range. The dielectric function displays rapid variations below 500 nm. Undoped silicon is metallic below 300 nm (shaded region), and its dispersion is marked by two adjacent resonances in the UV range [13], [31], [46], which according to Miller’s rule [53], [54] give rise to resonant nonlinearities in the opaque range. The system may then be described by two polarization components having two closely spaced resonance frequencies and separate linear and nonlinear spring constants, as follows:

Figure 1:

Schematic and linear properties of the silicon membrane. (a) Picture of a 200 nm-thick silicon film. The surface area of the suspended portion of the film (the area with greenish hue) is approximately 5 mm × 5 mm. (b) A pulse is incident from the top. We monitor only transmitted fields. (c) Linear transmittance measured at normal incidence (red markers) and fitted (blue curve) using Palik’s data in reference [31]. Some wavelengths of interest are marked in the figure. (d) Dielectric function of silicon as reported in reference [31] (markers) and as fitted by two Lorentzian functions (solid curves.)

Here, m b j * is the effective bound electron mass for the jth polarization component; ω ̃ 0 , b j = ω 0 , b j λ r / c , ω ̃ p j = ω p j λ r / c , γ ̃ b j = γ b j λ r / c are the scaled resonance and plasma frequencies, respectively, and damping coefficient; λ _r = 1 μm is a reference wavelength that scales ξ = z/λ _r , ζ = x/λ _r , ς = y/λ _r where x and y are the transverse coordinates, and z is the longitudinal coordinate; τ = ct/λ _r is the scaled time, and c is the speed of light in vacuum. Spatial and temporal derivatives on the fields are performed with respect to these new scaled spatial and temporal coordinates. The left-hand side of Eq. (1) contains internal linear and nonlinear restoring forces, while on the right-hand side we have forces resulting from externally applied fields.

Since we consider up to 7th harmonic generation, the bulk nonlinear response from bound electrons must reflect at least a 7th order nonlinearity, and may be written as:

Material response is assumed to be isotropic but could easily be modified to account for anisotropies [55] by isolating the spatial coordinates and by assigning the coefficients different values in different directions. In reference [55], which to our knowledge is the first report of reflected THG from silicon in the UV range, the authors investigated the question of THG using a nanosecond pump tuned to 1,064 nm and found an effective anisotropy in the third order nonlinear response. In the low intensity regime (which we arbitrarily set at 100 GW/cm² or less) integration of Eq. (1) together with Maxwell’s equations account for linear and nonlinear dispersions, and allow us to describe the process by considering the consequences of terms like e λ r m 0 , b j * c 2 P b j ⋅ ∇ E , which represents surface nonlinearities, and the magnetic Lorentz contribution e λ r m 0 , b j * c 2 P ̇ b j × H , which contains both surface and volume nonlinear bound currents. Both terms are also expanded up to their 7th harmonic contributions (not shown here), and partially account for even harmonic generation in centrosymmetric materials that act as insulators.

We now outline a novel approach to solving Eq. (1). The parameters β ̃ b j ≈ ω 0 , b j 2 λ r 2 / L 2 n 0 , b j 2 e 2 c 2 , δ ̃ b j = β ̃ b j / L 2 n 0 , b j 2 e 2 , and θ ̃ b j = β ̃ b j / L 4 n 0 , b j 4 e 4 are higher order scaled coefficients in a perturbative expansion that are derived directly from a nonlinear classical oscillator model, and are usually assumed to be constant, n _0,bj is the bound electron density, and L is usually interpreted as the lattice constant. From the point of view of a classical spring and classical macroscopic electrodynamics, the lattice constant L represents the spring’s maximum allowed extension, which is generally taken to be constant across the entire spectral range. From an atomic point of view, L may be interpreted as either an atomic/orbital diameter or inter-particle distance, which for solids can vary from a fraction of 1 Å to several Ås, a disparity that will be reflected in the particle density, and that will substantially affect the magnitudes of β ̃ , δ ̃ , and θ ̃ . However, since harmonic generation occurs over a wide wavelength range, all the parameters outlined above are unlikely to remain constant, including electronic effective masses, due to the different energy levels that are found deeper and deeper inside the atom, with complex and primarily unknown orbital radii, population density, and band curvatures. Therefore, orbitals and energy levels that are excited at the pump wavelength are different compared to states excited at each harmonic wavelength, will have different diameters and effective particle densities depending on the number of electrons in each orbital, so that at least the parameters L and n ₀ at the pump wavelength may be different at each of the harmonics. From a classical, macroscopic point of view, retrieval of β ̃ , δ ̃ and θ ̃ should be done in the same manner that the dielectric constant is retrieved, i.e. via ellipsometric techniques for a given thickness and wavelength. Then, it is reasonable to expect that ɛ, β ̃ , δ ̃ , and θ ̃ will display dispersive behavior not only as a function of wavelength, but also as a function of geometrical dispersion (resonances) [56]. Therefore, while it is generally acceptable to assume that β ̃ is constant under most circumstances [46], [53], [54], we may well have β ̃ = β ̃ 0 f β ̃ λ . Here, β ̃ 0 = ω 0 , b j 2 λ r 2 / L 0 2 n 0 , b j 2 e 2 c 2 is defined using the known lattice constant and particle density for the material, and f β ̃ λ is a parameter that generalizes the third order nonlinear coefficient in Eq. (1) that reflects the changing nature of the material across the spectrum. By the same token, we will assume that δ ̃ = δ ̃ 0 f δ ̃ λ and θ ̃ = θ ̃ 0 f θ ̃ λ , with δ ̃ 0 = β ̃ 0 / L 0 2 n 0 , b j 2 e 2 and θ ̃ 0 = β ̃ 0 / L 0 4 n 0 , b j 4 e 4 . Finally, given the proximity of the material resonances in Figure 1(d), for simplicity we will assume that the coefficients β ̃ b j , δ ̃ b j , and θ ̃ b j have similar amplitudes at both resonance wavelengths.

As we will see below, observation of the 7th harmonic at 210 nm using our detection system will require peak power densities of order 3 TW/cm². At these power densities and 10 Hz repetition rate our incident beam delivers an energy density of approximately 270 mJ/cm². Nevertheless, we expect single-pulse effects given the temporal displacement between consecutive pulses, each containing approximately 27 mJ/cm². The threshold for plasma formation in 3 μm-thick Si layers has been reported to be of order 150 mJ/cm² when it is delivered within ∼200 fs using optimized pulses and yielding remarkably high free electron densities of order 10²²–10²³/cm³ [57]. In our case, the sample survives undamaged if we keep peak power densities below 4 TW/cm². Our calculations suggest free carrier densities of order 10¹⁹–10²⁰/cm³. Plasma formation requires an additional, transient polarization component P _f to include a description of the dynamics of free carriers. This component triggers additional effects like pump shielding and an intensity dependent plasma frequency (increasing carrier density) that tends to blueshift as a function of incident peak power density. In turn this lowers the effective dielectric constant dynamically, affecting the structure and position of the geometrical resonances, including modifying the coefficients β ̃ , δ ̃ , and θ ̃ , which continue to be dispersive and have to be modified accordingly. Put another way, the increased free charge density comes at the expense of a decreasing bound charge density, causing the coefficients β ̃ , δ ̃ , and θ ̃ to change dynamically as a redistribution of charge occurs. Given these factors, the essential equation of motion that supplements Eq. (1) and describes the dynamics of free carriers may be written as follows [13], [40], [41], [46]:

Therefore, we account for surface and magnetic phenomena in both free and bound charges. Plasma frequency ω ̃ p , f and damping coefficients γ ̃ f are scaled similarly to their bound electron counterparts. The presence of the coefficient κ _{f ,} which may depend on fluence and may also be dispersive, allows the scaled plasma frequency ω ̃ p , f to change as a function of intensity (blueshift) and time, with contributions to all the harmonics, as revealed by an expansion of the fields up to 7th harmonic. The term proportional to the Fermi energy, E F = ℏ 2 2 m f * 3 π 2 n 0 f 2 / 3 , contains the effects of linear pressure and viscosity, m f * and n _0f are the free electron mass and background density, respectively. The terms in Eq. (3) thus reflect a continuity equation that yields contributions to an effective plasma frequency that depends on the local charge density ∇ ⋅ P _f and a third order nonlinearity as follows [58]

Nonlocal terms ∇ ∇ ⋅ P f + 1 2 ∇ 2 P f then make the dielectric function a highly dynamic variable.

In Figure 2 we show the measured and simulated spectra generated by a ∼3 TW/cm², 85 fs pump pulse tuned to 1,475 nm incident normal to the surface to temporarily suppress even harmonics. The measured conversion efficiencies are 4.2 × 10⁻⁵, 6.2 × 10⁻⁷, and 4.57 × 10⁻⁹ for 3rd, 5th, and 7th harmonic, with peak power densities of order 126 MW/cm², 1.8 MW/cm² and 14 kW/cm², respectively, despite being tuned in the absorption and metallic ranges. Agreement between the measured and simulated data is remarkable. It is worth mentioning that the experimentally recorded side lobes contain less than 0.1 % of their respective harmonic energy and can be attributed to small variations in temporal and spatial profiles of the experimental source when compared to the perfectly formed Gaussian pulses used in simulations. Phase-locking is the fundamental mechanism that makes this kind of interaction possible. Our simulations suggests that results like those depicted in Figure 2 can be obtained by tuning the pump field at each of the resonances shown in Figure 1(c), pushing the highest harmonic we can presently simulate below 100 nm. See Supplementary Information for a scan of 7th harmonic generation efficiency for a range of wavelength that includes the resonances shown in Figure 1(c).

Figure 2:

Measured and simulated spectra of transmitted 3rd, 5th, and 7th harmonic conversion efficiencies for a single pulse 85 fs in duration, having carrier wavelength tuned to 1,475 nm, peak power density of 3 TW/cm², and normally incident on the silicon film. Relevant parameters: β ̃ ω = 5 × 1 0 − 10 , β ̃ 3 ω = 2 × 1 0 − 10 , β ̃ 5 ω = 5 × 1 0 − 11 , β ̃ 7 ω = 5 × 1 0 − 11 , δ ̃ ω = 1.25 × 1 0 − 19 , δ ̃ 3 ω = 5 × 1 0 − 20 , δ ̃ 5 ω = 1.25 × 1 0 − 20 , δ ̃ 7 ω = 1.25 × 1 0 − 20 , ϑ ̃ ω = 6.25 × 1 0 − 29 , ϑ ̃ 3 ω = 2.5 × 1 0 − 29 , ϑ ̃ 5 ω = 6.25 × 1 0 − 30 , ϑ ̃ 7 ω = 6.25 × 1 0 − 30 , k f ω = 1.6 × 1 0 − 9 , k f 3 ω = 6.4 × 1 0 − 10 , k f 5 ω = 1.6 × 1 0 − 10 , k f 7 ω = 1.6 × 1 0 − 10 .

We then pump the sample at 1,550 nm at normal incidence and perform an intensity scan to detect third and fifth harmonic signals at 516 nm and 310 nm, which is our detection limit for the low-intensity setup. In Figure 3(a) and (b) we show the result of both scans, along with the result of our simulations. Conversion efficiency is defined by dividing the total energy of a given harmonic by the total incident pump energy. Once the dielectric constant is established and fitted using the linearized version of Eq. (1) [13], β ̃ and δ ̃ are the only remaining parameters to be determined. The reconstruction of nonlinear dispersion when δ ̃ is not zero is complicated because the δ ̃ terms tend to quench the β ̃ contributions at both the pump and its harmonic wavelengths, given their opposite signs. This effect is discussed in reference [59]. The detection system in Figure S1(a) can be used to detect signals down to approximately 310 nm. Therefore, a 7th harmonic was also observed when the pump pulse was tuned to a carrier wavelength of 2,150 nm. We show the measured 7th harmonic efficiency versus incident peak power density in Figure 3(c). We note that the transmission function in Figure 1(c) expresses the behavior of a pump field tuned to that wavelength and says nothing about the behavior of a harmonic signal generated and tuned to the same wavelength. This conclusion is fully supported by the fact that the generated 5th harmonic signal is transmitted by the etalon, notwithstanding the fact that the linear transmission function suggests that wavelength should be fully suppressed. We emphasize that transmission of the 210 nm and 310 nm signals occurs because they are the inhomogeneous components of the generated harmonics, which propagate under phase locking conditions, i.e., with the dispersive properties of the pump field tuned in the transparency range of the material.

Figure 3:

Transmitted third (a) and fifth (b) harmonic generation efficiency versus incident peak power density. An intensity scan is performed with a ∼100 fs pulses tuned to 1,550 nm with repetition rate 2 Hz. The 3rd (516 nm) and 5th (310 nm) harmonic data are fitted almost perfectly by the simulations. Relevant parameters: β ̃ 3 ω = β ̃ 0 = 8.5 × 1 0 − 9 , δ ̃ 3 ω = δ ̃ 0 = 4.25 × 1 0 − 18 , ϑ ̃ 3 ω = ϑ ̃ 0 = 2.125 × 1 0 − 27 , β ̃ 5 ω = β ̃ 0 , δ ̃ 5 ω = 9 δ ̃ 0 , ϑ ̃ 5 ω = 5 ϑ ̃ 0 , β ̃ 7 ω = 4 β ̃ 0 , δ ̃ 7 ω = δ ̃ 0 , ϑ ̃ 7 ω = ϑ ̃ 0 . (c) Measured 7th harmonic conversion efficiency versus peak power density for the setup up in Figure S1(a). The fit follows the 5th power of the peak power density. The field is tuned off resonance to 2,150 nm, is incident at 30°, with 20 Hz repetition rate.

We now focus on the geometrical resonance centered near 1,400 nm and scan the carrier wavelength of the pump field across it. The result is plotted in Figure 4, where we show measurements and simulations of the THG conversion efficiency as a function of wavelength, as well as the dispersive character of β ̃ across the resonance – Figure 4(b). Figure 4(a) contains two simulated curves, one obtained with a constant β ̃ at all wavelengths, and one calculated using the dispersive coefficient β ̃ 3 ω λ shown in Figure 4(b). While a constant β ̃ yields good qualitative agreement across the resonance, it is only when using our full dispersive model that qualitative and quantitative agreement can be achieved. This is clearly shown in Figure 4(a) by a direct comparison of experiment and the two models with constant and dispersive β ̃ 3 ω , respectively. At resonance, the dispersive coefficient β ̃ 3 ω is nearly 1.5 times larger compared to its unitary value off resonance, shown in Figure 4(b). The experimental results clearly point towards a dispersive variation of the model and thus reveal that at least from a classical, macroscopic electrodynamics point of view, theory can be easily reconciled with measurements if we assume that β ̃ is dispersive across the resonance.

Figure 4:

Transmitted THG conversion efficiency for (a) constant β (dashed blue curve) and with a β _3ω that varies as a function of wavelength (solid blue curve) as in (b). The red markers represent the measured data. (b) Actual functional dependence of the coefficient β _3ω used at the third harmonic wavelength. Relevant parameters: β ̃ 3 ω = β ̃ 0 = 8.5 × 1 0 − 9 ; δ ̃ 3 ω = δ ̃ 0 = 4.25 × 1 0 − 18 ; ϑ ̃ 3 ω = ϑ ̃ 0 = 2.125 × 1 0 − 27 , β ̃ 5 ω = β ̃ 0 ; δ ̃ 5 ω = 9 δ ̃ 0 ; ϑ ̃ 5 ω = 5 ϑ ̃ 0 , β ̃ 7 ω = 4 β ̃ 0 ; δ ̃ 7 ω = δ ̃ 0 ; ϑ ̃ 7 ω = ϑ ̃ 0 .

In Figure 5(a) we show 4th harmonic generation conversion efficiency as a function of incident peak power density for fixed angle. In Figure 5(b) we fix the incident peak power density and perform an angular scan to reveal the angular dependence of the 4th harmonic. Once again, the agreement between experiment and theory is remarkable and, to the best of our knowledge, unprecedented. The only free parameters are the effective free electron masses at the different harmonic wavelengths, and the bulk coefficients in Eq. (2). SHG results for similar samples were previously reported [46], while the 6th harmonic was not detectable for the intensities we used in this setup (up to 100 GW/cm²).

Figure 5:

Transmitted (a) 4th harmonic generation conversion efficiency versus incident peak power density with angle of incidence fixed at 30°, carrier wavelength of 1,550 nm and repetition rate of 20 Hz. The red solid markers denote the data, while the simulations are reported as a solid blue curve. (b) Measured and simulated angular dependence of the 4th harmonic signal as a function of incident angle for fixed incident peak power density of 75 GW/cm². Data are shown with red markers and simulations are shown as blue solid curves. Relevant parameters: m b * 2 ω = 0.24 m e ; m b * 4 ω = 0.002 m e ; β ̃ = 1.5 × 1 0 − 9 ; δ ̃ = 2.25 × 1 0 − 18 ; ϑ ̃ = 2.25 × 1 0 − 29 .

In Figure 6(a) we plot transmitted THG conversion efficiency versus incident peak power density for the 85 fs pulse tuned to 1,475 nm. Both measured data and simulation display a maximum near 2 TW/cm², followed by a decrease in efficiency. We refer to Figure 6(b) and (c) to explain this behavior. In Figure 6(b) we simulate and plot the effective dielectric constant as a function of time as the pulse traverses the thin film, for a peak power density of approximately 1 TW/cm². The simulation reveals a maximum change of order R e δ ε ∼ − 4 when the peak of the pulse reaches the center of the layer. The extraction of this parameter and the method is described in detail in reference [13] and requires knowledge of the constitutive relations. Suffice it to say here that in essence we perform linear and nonlinear numerical ellipsometry on the sample, by calculating the macroscopic dielectric constant as ε t = 1 + 4 π P ( t ) E ( t ) , where the brackets stand for spatial average inside the layer at each instant of time. This definition coincides with the result of an experimental ellipsometric measurement, as demonstrated in reference [13]. Then, the complex change of the dielectric constant may be calculated as the difference between the high and the low intensity dielectric functions as follows: δ ε = ε N L − ε L = 4 π P N L E N L − P L E L , where the subscript L and NL stand for linear and nonlinear, respectively. A dynamic reduction of the dielectric constant from approximately 13 to 9 at the pump wavelength causes an intensity dependent blueshift of the resonance, shown in Figure 6(c), effectively detuning the pump out of resonance and reducing the conversion efficiency.

Figure 6:

Measured and simulated (a) transmitted THG efficiency at high intensities. Pump is tuned at 1,475 nm with a pulse duration of 85 fs; (b) simulated extraction of the effective dielectric constant from the constitutive relations following the procedure developed in reference [13], for peak power densities of order 1 TW/cm². (c) Linear (blue curve) and nonlinear (red curve) transmission functions showing the blueshift of the resonance for a 200 nm thick etalon. The relevant parameters are: β ̃ ω = 5 × 1 0 − 10 ; β ̃ 3 ω = 2 × 1 0 − 10 ; β ̃ 5 ω = 5 × 1 0 − 11 ; β ̃ 7 ω = 5 × 1 0 − 11 , δ ̃ ω = 1.25 × 1 0 − 19 ; δ ̃ 3 ω = 5 × 1 0 − 20 ; δ ̃ 5 ω = 1.25 × 1 0 − 20 ; δ ̃ 7 ω = 1.25 × 1 0 − 20 ; ϑ ̃ ω = 6.25 × 1 0 − 29 ; ϑ ̃ 3 ω = 2.5 × 1 0 − 29 ; ϑ ̃ 5 ω = 6.25 × 1 0 − 30 ; ϑ ̃ 7 ω = 6.25 × 1 0 − 30 ; k f ω = 1.6 × 1 0 − 9 ; k f 3 ω = 6.4 × 1 0 − 10 ; k f 5 ω = 1.6 × 1 0 − 10 ; k f 7 ω = 1.6 × 1 0 − 10 .

Figure 6(b) reveals an exceptionally powerful spatio-temporal dynamics, i.e., time varying capability that may rival and even surpass that of transparent conductive oxides by exploiting resonant nonlinearities and plasma formation. While conductive oxides may seem to hold the advantage in transparency, the transmission resonances and field localization that can be created in a Fabry–Perot cavity or in nanowire arrays [19] may also lead silicon and other semiconductors to serve as alternative building blocks of time varying metamaterials.

At this point, a note about parameter space is in order. The total nonlinear bulk polarization that is inserted into Maxwell’s equation is the sum of the combined solutions of Eqs. (1) and (3), while the instantaneous nonlinear potential may be thought of as a combination of free and bound polarizations as follows:

When the fields and polarizations are expanded up to the 7th harmonic frequency, each term in Eq. (5) contributes to each of the harmonics, including even harmonics. For instance, some of the terms that oscillate at the fundamental wavelength are proportional to κ f E ω 2 E ω , β ̃ j P ω 2 P ω , δ ̃ b j P ω 4 P ω , and θ ̃ b j P ω 6 P ω . Similarly, just a few of the terms that oscillate at the second harmonic are κ f E ω 2 E 2 ω , β ̃ j P ω 2 P 2 ω , δ ̃ b j P ω 4 P 2 ω * , and ϑ ̃ b j P ω P ω 4 P 2 ω * . At the third harmonic wavelength some of the terms are proportional to κ f E ω 3 , β ̃ j P ω 3 , δ ̃ b j P ω 2 P ω 3 , θ ̃ b j P ω 4 P ω 3 . The leading cubic terms can balance each other, while 5th and 7th powers of the fields tend to either quench or reinforce terms of lower order, depending on the magnitude of the coefficients, as well as material dispersion, since the polarizations intrinsically contain the susceptibilities, which may be either positive or negative depending on tuning. Therefore, the set of parameters that we use in our simulations can be modified somewhat and still yield similar results.

3 Conclusions

In summary, we have reported odd and even nonlinear frequency conversion up to the 7th harmonic from a subwavelength silicon film. Several of the harmonics are generated in the opacity range, where the material displays metallic behavior with a negative dielectric constant. Those harmonics are transmitted thanks to a phase locking mechanism that allows the inhomogeneous component to effectively defeat nonlinear absorption as it resonates and propagates through thick layers, if the pump is tuned in the transparency range. Our measurements and simulations show that harmonic generation deep in the UV range is possible, can be efficient under resonant conditions [19], 46], and Supplementary Information], and the sample can withstand peak power densities well above 1 TW/cm². It is remarkable that the reported conversion efficiencies are quite high deep in the UV range, where the material is metallic and absorptive. It is notable that the efficiencies observed for silicon membranes, although impressive and somehow unexpected, are still at least one order of magnitude lower (across all harmonics) than those reported in transparent conductive oxides films pumped near the ENZ crossover wavelength [60]. It should also be pointed out, that here we achieve these results by exploiting a resonant structure and more advanced fabrication processes [61]. Our simulations and measurements also demonstrate that silicon exhibits extraordinary time-varying capabilities that could potentially rival or even exceed those of transparent conductive oxides by leveraging resonant nonlinearities and plasma formation. While conductive oxides might have an edge in terms of transparency, any cavity environment can yield high transparency and large local fields, thus suggesting that silicon and other semiconductors could serve as viable alternative building blocks for time-varying metamaterials. Finally, we note that while our theoretical and experimental results show a remarkable degree of agreement, we are unable to find an apt comparison with other published work on the subject. For instance, high harmonic generation in solids and silicon is typically studied using time-dependent density functional theory, or TDDFT [62], [63], [64], or in environments where plasmonic antennas are placed on top of a silicon layer to enhance performance [65]. In references [62], [63], [64] odd harmonics are reported and featured, but no reference is made to even harmonic generation. Even harmonic generation is briefly mentioned in reference [65] but only to point out that it is not observed. Some speculation is offered to describe circumstances where it might appear, for example if inversion symmetry is broken for sufficiently large fields across electron-hole pairs. However, in all cases even harmonic generation is never quantified or described in terms of surfaces, intrinsic magnetism triggered by the Lorentz force, or convection. Therefore at least in this regard a direct comparison is not possible. In addition, while some TDDFT methods assume the interaction occurs in a single cell where the field is spatially uniform [62], other TDDFT approaches neglect propagation effects [64]. In our approach we do neither: fields are not uniform in cavity environments, since they are characterized by field localization. Therefore, a direct comparison with odd harmonic generation is currently unavailable.