1. Introduction

Terahertz (THz)-field-induced second-harmonic generation (TFISH), a nonlinear optical process driven by THz fields that is inherently sensitive to surfaces and interfaces^[1], has emerged as a powerful tool for probing ultrafast symmetry breaking and lattice dynamics in solids. By enabling direct mapping of interfacial inversion symmetry, carrier redistribution, and ionic polarization on sub-picosecond timescales, TFISH provides critical capabilities for high-speed photonic modulation and phase-sensitive detection. Such capabilities are closely relevant to the study of ultrafast excitation and detection of nonequilibrium states and emergent phenomena in materials systems^[2]. When combined with an infrared (IR) probe, TFISH generates a four-wave mixing signal via third-order nonlinear effects $P_{\text{TFISH}}^{\left(2\right)}\left(2\omega \right)\propto {\chi }^{\left(3\right)}{E}_{\text{TFISH}}{E}_{\text{IR}}^2\left(\omega \right)\propto {\chi }_{\text{eff}}^{\left(2\right)}{E}_{\text{IR}}^2\left(\omega \right)$^[3], this which is equivalent to introducing an additional field-induced second-order polarization, thereby realizing the generation of a THz-controlled effective second-order nonlinear coefficient ${\chi }_{\mathrm{eff}}^{\left(2\right)}$. The TFISH-induced ${\chi }_{\mathrm{eff}}^{\left(2\right)}$ is strongly enhanced and spatially concentrated at polarized surfaces where inversion symmetry is naturally broken and local fields are strongly enhanced^[4]. In wide-bandgap polar crystals where carrier contributions are suppressed and lattice nonlinearities dominate, TFISH thus emerges as a uniquely clean probe of THz–optical interactions and interface-specific nonlinear polarization.

While TFISH has been extensively studied in centrosymmetric semiconductors such as Si and Si/SiO₂ structures^[5,6], where bulk ${\chi }_0^{\left(2\right)}$ vanishes and the second-harmonic generation (SHG) signal arises almost exclusively from interfaces, these platforms are not ideal for exploring the interplay between intrinsic bulk nonlinearity, ionic responses, and THz-driven surface phenomena. Moreover, existing TFISH studies have rarely addressed how surface-to-surface interference depends on crystal thickness and dispersion, nor how this interference can be actively manipulated for sensing and device applications.

Wurtzite ZnO offers an ideal model system: ZnO is a wide bandgap semiconductor (~3.37 eV) with a non-centrosymmetric hexagonal wurtzite structure (space group P6₃mc)^[7,8]. Its second-order susceptibility ${\chi }_0^{\left(2\right)}$ gives rise to a nonlinear polarization, which in turn leads to SHG^[9,10]. The broken symmetry at polar surfaces (e.g., Zn- and O-terminated facets) induces localized nonlinear polarization amplification, rendering SHG a potent method for ultrafast surface/interface dynamics^[11,12]. Quasi-resonant excitation of ZnO near its bandgap with 400 nm pulses has been shown to launch coherent lattice vibrations involving the E₂, A₁(LO), and E₁ modes, which in turn modulate the nonlinear polarization and give rise to temporal oscillations in the transient reflectance signal^[13]. Recent works have shown that not only phonons but also phonon polaritons (PhPs) can be excited under THz pulses^[4,14,15]. Analogous quasiparticle dynamics, such as Fermi polarons in transition metal dichalcogenides, have been extensively studied using ultrafast spectroscopy, revealing valley coherence and decoherence on sub-picosecond timescales^[16]. Due to its wide bandgap structure, it can suppress the contribution of free carriers under IR detection energy while providing a broad THz transmission window^[17,18]. The phase-matching conditions for PhPs in ZnO can be tuned over a broad spectral range by adjusting the THz frequency and IR probe wavelength, enabling flexible access to different PhP branches. Besides, its anisotropic crystal structure enables SHG conversion efficiency to be optimized through crystal orientation design^[19-21]. These characteristics make ZnO an ideal model system for investigating how THz-driven nonlinear polarizations at different crystal surfaces interfere with each other, and how surface-specific TFISH signals emerge on top of a finite bulk ${\chi }_0^{\left(2\right)}$ background.

Here, we experimentally demonstrate that THz fields can actively modulate the coherent second-harmonic (SH) interference between the front and rear surfaces of ZnO crystals, revealing that THz-driven PhPs can govern the characteristics of surface SHG. This mechanism can be extended to manipulate surface-to-surface SH interference via THz fields for practical applications. Using millimeter-scale ZnO crystals under single-cycle THz excitation, we observe energy-resolved SHG interference fringes originating from the two surfaces and demonstrate that they arise from THz-field-induced surface-to-surface SHG interference. By correlating the temporal separation between the front and rear-surface TFISH bursts with the spectral interference pattern, the time-domain dynamics are directly linked to the energy-domain fringe structure. A phase-accumulation model and numerical simulations are developed that quantitatively reproduce the observed fringes and clarify the distinct roles of THz-driven nonlinearities, which determine the fringe visibility and temporal gating, and optical dispersion, which fixes the fringe spacing via the propagation between the two surfaces. These results establish TFISH-based surface-to-surface SHG interferometry as a sensitive route to phase-resolved thickness metrology and in situ monitoring of refractive-index changes, while also enabling ~50 fs IR probes to access THz-driven phase dynamics far beyond the naive time-resolution limit.

2. Experimental Setup

A pump-probe experiment is implemented, whereby the pump THz beam and the probe IR beam are incident collinearly into ZnO crystals. A Ti:sapphire laser with a central wavelength of 800 nm, a repetition rate of 1 kHz, and a pulse duration of 35 fs, is employed to drive an optical parametric amplification (OPA) that produces IR pulses with a center wavelength of 1,450 nm, as shown in Figure 1c. The THz radiation, shown in Figure 1b, is generated through optical rectification within an organic crystal, 4-N,N-dimethylamino-4’-N’-methyl-stilbazolium 2,4,6-trimethylbenzenesulfonate (DSTMS), which is pumped by the IR beam from the OPA, as shown in Figure 1a. The residual IR light that pumps the DSTMS was filtered out by a 2 mm-thick high-density polyethylene (HDPE) window, while the generated THz radiation could almost completely transmit through the HDPE to excite the sample. The THz electric field is characterized via electro-optical sampling with a 50 μm-thick GaP crystal, reaching a peak field of ~100 kV/cm. The TFISH signal is acquired using a spectrometer. Time-resolved TFISH measurements are performed by varying the temporal delay between the THz and IR pulses. The experiment is conducted using [0001]-cut ZnO crystals with thicknesses of 0.3 mm, 0.4 mm, 0.5 mm, 0.6 mm, and 1 mm. The pump and probe beams were incident perpendicularly onto the sample, with their temporal overlap adjusted via a motorized delay stage. The SHG signal generated by the probe was recorded using a fiber-coupled spectrometer, enabling time-resolved measurements of the pump-induced dynamics across the sample.

Display Full Size

Download

Figure 1. Experimental optical setup. (a) Schematic experimental setup for time-resolved THz-induced SHG interference fringe measurements in ZnO. A perforated off-axis parabolic mirror was used to combine the THz pump and IR probe collinearly at the ZnO sample position. The IR probe passed through the central hole, while the THz beam was reflected and focused by the mirror onto the sample; (b) Corresponding frequency spectrum of the THz pulse, obtained by FFT of the time-domain signal; (c) The IR light with a central wavelength of 1,450 nm used in the experiment. THz: terahertz; SHG: second-harmonic generation; IR: infrared; FFT: fast Fourier transform; OPA: optical parametric amplification; DSTMS: 4-N,N-dimethylamino-4’-N’-methyl-stilbazolium 2,4,6-trimethylbenzenesulfonate; HDPE: high-density polyethylene.

3. Results and Discussion

When the pump THz pulse is temporally synchronized with the probe IR beam in the ZnO sample, the TFISH signal is observed. The THz-induced second-order nonlinear polarization is given by:

(1)$P_{\mathrm{TFISH}}^{\left(2\right)}(2\omega ,t)={\epsilon }_0{\chi }^{\left(3\right)}{E}_{\mathrm{THz}}\left(t\right)E_{\mathrm{IR}}(\omega ,t)E_{\mathrm{IR}}(\omega ,t)$

where ${\chi }^{\left(3\right)}$ denotes the additional nonlinear susceptibility arising from symmetry breaking induced by THz polarization, E_IR(ω,t) is the fundamental frequency electric field, and E_THz(t) is the THz electric field. As zinc oxide is inherently a non-centrosymmetric crystal, it possesses the intrinsic capability to generate SHG, and its own second-order nonlinear coefficient is denoted by ${\chi }_0^{\left(2\right)}$. Therefore, the experimentally observed surface SHG arises from two channels: (i) the intrinsic second-order susceptibility ${\chi }_0^{\left(2\right)}$, and (ii) an effective second-order response ${\chi }_{\mathrm{eff}}^{\left(2\right)}$, induced by the THz field, arising from the third-order susceptibility ${\chi }^{\left(3\right)}$ coupled to the THz electric field.

Figure 2a,b,c,d,e display the time-domain SHG signals from ZnO samples with five different thicknesses (0.3, 0.4, 0.5, 0.6, and 1.0 mm) under THz electric field excitation, clearly demonstrating the pronounced thickness dependence of THz-induced SHG modulation. For all thicknesses, two distinct groups of SHG signal peaks can be identified along the time axis (white dashed lines). The first group, centered at τ = 0 ps, originates from the ZnO front surface (air–ZnO interface), while the second group arises from the rear surface. The corresponding time delays are τ = 0.95, 1.25, 1.65, 1.9, and 3.3 ps for sample thicknesses of d = 0.3, 0.4, 0.5, 0.6, and 1.0 mm, respectively. As indicated by the two white dashed lines, the front and rear-surface signals originate from the TFISH processes occurring at the respective interfaces, i.e., the front-surface TFISH and the rear-surface TFISH, respectively. This observation establishes the temporal separation between the two surfaces. Such a delay further provides the basis for the phase accumulation that later manifests as interference fringes in both the time and energy domains. By comparing with the THz waveform (Figure 2f), it is further confirmed that the front-surface signal at τ = 0 ps and the rear-surface signal peaks coincide with the extrema of the THz electric field, indicating that both contributions are directly driven by the instantaneous THz field.

Display Full Size

Download

Figure 2. Time-resolved SHG signals from ZnO samples with varying thickness. (a-e) SHG signals for ZnO thickness of (a) 0.3 mm, (b) 0.4 mm, (c) 0.5 mm, (d) 0.6 mm, and (e) 1 mm. The time scales marked as I, II, and III correspond to the situations where the THz pulse has not yet arrived, is propagating inside the sample, and has transmitted through the rear surface, respectively; (f) Corresponding THz waveform in time domain (the peak field intensity of THz pulse is ~100 kV/cm); (g-i) schematically illustrate the interference of the two SHG beams when the THz pulse arrives at the front surface, inside, and after the rear surface of the sample, respectively; (j) In the region near the rear surface, the originally horizontal fringes become tilted for the first time under the influence of the PhP. SHG: second-harmonic generation; THz: terahertz; PhP: phonon polariton; SH: second-harmonic; IR: infrared; TFISH: THz-field-induced second-harmonic generation.

The temporal separation between the front and rear-surface SHG signals is governed by the group velocity mismatch between the THz pump pulse and the IR probe pulse propagating through the ZnO crystal. This delay τ_th can be quantitatively expressed as:

(2)${\tau }_{\mathrm{th}}=(n_{\mathrm{THz}}-n_{\mathrm{IR}})\frac{d}{c}$

where n_THz (near 3.5 THz, corresponding to the phase-matched frequency for the SHG of the 1,450 nm probe at its central wavelength)^[17] and n_IR ≈ 1.93 (at 1,450 nm) are the refractive indices of ZnO for the pump THz pulse and probe IR pulse, d denotes the thickness of ZnO samples, and c is the speed of light. Experimental measurements of τ show good agreement with theoretical predictions, with discrepancies remaining relatively low across all sample thicknesses. Using Eq. (2) together with the refractive indices n_THz and n_IR of ZnO, the expected delay times for the THz–IR overlap at the rear surface are calculated to be 0.96, 1.28, 1.6, 1.92, and 3.2 ps for crystal thicknesses of 0.3, 0.4, 0.5, 0.6, and 1.0 mm, respectively, in good agreement with the experimental peak positions shown in Figure 2a,b,c,d,e. The use of a nearly single-cycle THz pulse ensures that the TFISH bursts from the front and rear surfaces are well separated in time, which enables a direct correlation between their temporal separation and the spectral interference fringes discussed below.

Noting that in the transmission geometry used here, the detected SHG signal has passed sequentially through both the front and rear surfaces, it is therefore crucial to clarify the microscopic origins of the SHG contributions throughout the entire spectral range. Microscopically, three distinct contributions to the observed SHG signal can be identified: (i) the intrinsic ${\chi }_0^{\left(2\right)}$ driven SHG of the IR probe, (ii) the TFISH component launched at the front surface and propagating through the crystal, and (iii) the TFISH component generated locally at the rear surface by the THz-driven PhP field.

Based on the temporal evolution of the SHG signals, three distinct stages can be identified, as schematically illustrated in Figure 2g,h,i. These stages represent different interference scenarios occurring sequentially in time. The interference between the front-surface’s own SHG signal and the TFISH-generated SHG (Figure 2g), the interference between the SHG generated by the forward-propagating PhP acting synchronously with the probe and the intrinsic SHG (Stage II, Figure 2h), and the interference between the SHG generated by the PhP at the rear surface (including its reflected backward-propagating component) acting synchronously with the probe and the intrinsic SHG (Stage III, Figure 2i).

In particular, the regime corresponding to Stage III (Figure 2i) is directly evidenced by the experimental data shown in Figure 2j, where a clear splitting along the SHG energy axis is observed. Compared with the signal in Stage II, the SHG intensity in the delay region after the rear-surface overlap is noticeably stronger. This indicates that, although both Stage II and Stage III involve phase-matching processes associated with forward and backward-propagating PhPs, respectively, the backward propagating phase-matching condition after the rear surface is more favorable in the present experiment. This point will be discussed in detail below. More importantly, Figure 2j exhibits a distinct delay-dependent tilted-fringe pattern, which does not appear in the absence of THz excitation. This THz-induced tilting is a key signature that the interference is no longer governed solely by the static front–rear optical phase accumulation, but is additionally modulated by the PhP-related backward-matching contribution that becomes dominant after the rear-surface delay.

Having established the temporal origin and interference mechanism of the SHG signals in Figure 2, we now turn to a detailed analysis of the spectral interference features observed along the SHG energy axis. In particular, the splitting and modulation patterns provide direct insight into the underlying phase accumulation and PhP-assisted matching processes.

As shown in Figure 3a, well-defined interference fringes emerge along the SHG energy axis, and their spacing exhibits a clear dependence on the crystal thickness. Figure 3b further presents the delay-dependent SHG map of the 0.6 mm-thick sample, showing the spectral response without THz excitation (stage I) and with THz excitation (stage II). In the pre-zero-delay region, where the THz pump is absent, the interference pattern is governed only by the intrinsic propagation-induced phase accumulation. This is analyzed more explicitly in Figure 3c, which shows that the peak-to-peak energy spacing of the intrinsic SHG remains nearly unchanged as the delay varies. This indicates that, within the spectral window relevant to our measurements, the observed ΔE is primarily determined by the optical phase accumulation between the front- and rear-surface SHG contributions, rather than by delay-dependent absorption of the IR probe or the SHG signal. Meanwhile, the signal shown in Figure 3a corresponds to the SHG response at the rear-surface delay. Therefore, it can be regarded as the starting point at which the backward phase-matching contribution becomes involved after the arrival of the THz field, and the energy splitting ΔE at this onset remains consistent with that in the absence of THz excitation. The behavior at delays beyond the rear-surface position will be discussed in the following section. Due to the difference in refractive indices between the fundamental IR field and the generated SHG field, a phase mismatch is accumulated along the propagation path, leading to an intrinsic phase difference between these two contributions. The characteristic energy spacing ΔE between adjacent maxima therefore follows the fundamental relation:

Display Full Size

Download

Figure 3. Interference fringe characteristics of the intrinsic SHG signal of varying thicknesses. (a) Time-integrated SHG spectra of samples with thicknesses of 0.3, 0.4, 0.5, 0.6, and 1.0 mm exhibit clear interference fringes along the photon-energy axis, and these fringes are observed at the delay of rear surface. The measured energy spacings between adjacent peaks near 1.71 eV (the centra photon energy of SHG) are: ΔE₁ = 0.030 eV, ΔE₂ = 0.0250 eV, ΔE₃ = 0.0210 eV, ΔE₄ = 0.0154 eV, ΔE₅ = 0.0094 eV respectively; (b) The SHG spectrum of the 0.6 mm-thick crystal without (stage I) and with (stage II) the THz excitation; (c) Analysis of the intrinsic SHG spectral data for the region before zero delay in panel (b). Without the THz pump, the peak energy spacing of the intrinsic SHG remains nearly unchanged as the delay increases. With the THz pump, the PhP-induced SHG observed at the rear surface of the sample (which marks the onset of backward-propagating PhP) exhibits a shift in the peak energy, while the energy spacing remains unchanged. SHG: second-harmonic generation; THz: terahertz; PhP: phonon polariton; SH: second-harmonic.

(3)$\Delta E\left(E\right)\approx \frac{hc}{d\left(|n_{2\omega }-n_{\omega }|\right)+E\frac{\partial \left(|n_{2\omega }-n_{\omega }|\right)}{\partial E})}$

where E = 2ħω (which expresses the photon energy of SHG), d denotes the ZnO thickness, n_ω and n_2ω are the refractive indices of ZnO at the fundamental frequency ω and the SH frequency 2ω respectively. The refractive index difference |n_2ω - n_ω| accounts for material dispersion. This inverse-thickness dependence (ΔE ∝ 1/d) is quantitatively confirmed by our measurements: as the sample thickness increases, the accumulated phase difference between the front and rear surfaces grows, leading to progressively denser interference fringes in energy.

The relation for ΔE can be understood by considering the phase difference Δφ(E) between SHG₁ and SHG₂ as a function of photon energy. For a given energy, the phase difference can be written as:

(4)$\Delta {\phi }_{\mathrm{SHG}}={\phi }_{{\mathrm{SHG}}_1}-{\phi }_{{\mathrm{SHG}}_2}=2({\phi }_f+\omega d\frac{n_{2\omega }}{c})-2({\phi }_f+\omega d\frac{n_{\omega }}{c})$

where φ_f is the phase of the fundamental-frequency IR light, while ${\phi }_{{\mathrm{SHG}}_1}$ and ${\phi }_{{\mathrm{SHG}}_2}$ are the phases of the two interfering SH waves, respectively. The factor of 2 in their phases reflects the quadratic dependence of the SHG polarization on the fundamental light field $P^{\left(2\right)}\left(2\omega \right)\propto E_{\mathrm{IR}}^2\left(\omega \right)$, so that the SHG field inherits twice the optical phase of the driving field. The phase difference between the two interfering SHG fields can be written as follows.

(5)$\Delta {\phi }_{SHG}\left(E\right)=\left(\frac{2d}{hc}\right)\left(|n_{2\omega }-n_{\omega }|\right)E$

Therefore, according to the condition for neighboring maxima, Δφ(E + ΔE)- Δφ(E) = 2π, together with the Sellmeier relation, Eq. (3) is obtained. According to the Sellmeier relation $n^2=2.81418+\frac{0.87968{\lambda }^2}{{\lambda }^2-{0.3042}^2}-0.00711{\lambda }^2$^[22], |n_2ω - n_ω| is not a strict constant, but varies weakly with wavelength and therefore with the SHG photon energy. Since the probe pulse used in our experiment has a specific spectral bandwidth, as shown in Figure 1c, the generated SHG correspondingly covers a certain photon-energy range. As a result, the interference oscillation period is not strictly identical across all spectral components, but varies slightly over the SHG spectrum. Within the energy window relevant to our measurements, however, this variation remains relatively small. Therefore, for a direct comparison with the experiment, we selected several representative peak-to-peak spacings around the central SHG photon energy of ~1.71 eV. As shown in Figure 3a, ΔE₁, ΔE₂, ΔE₃, ΔE₄ and ΔE₅ denote the SHG energy splittings for ZnO crystals with thicknesses of 0.3, 0.4, 0.5, 0.6, and 1.0 mm, respectively. Substituting the thickness parameters into Eq. (3) yields energy splitting intervals near the center energy of 1.71 eV of 0.0321 eV (0.3 mm), 0.0241 eV (0.4 mm), 0.0193 eV (0.5 mm), 0.0161 eV (0.6 mm), and 0.0096 eV (1.0 mm). The calculated values are in excellent agreement with the experimental values measured in Figure 3.

Next, we will discuss the influence of the THz field on the interference fringes when it arrives. When the THz field is introduced, an additional third-order nonlinear contribution ${\chi }^{\left(3\right)}$ is incorporated into the interference framework that is otherwise driven solely by the IR field. Within this picture, the thickness dependence of the spectral fringe spacing on the energy scale can be described by an effective nonlinear polarization, which includes both the intrinsic second-order susceptibility ${\chi }_0^{\left(2\right)}$ and an effective second-order contribution mediated by the third-order susceptibility ${\chi }^{\left(3\right)}$ in the presence of the PhP field E_PhP.

(6)${\chi }^{\left(2\right)}(t,z)={\chi }_0^{\left(2\right)}+{\chi }^{\left(3\right)}{E}_{\mathrm{PhP}}(t,z)$

In analyzing the TFISH response, we therefore focus on the nonlinear mixing pathways that contribute to the SHG signal. Besides the bare 2ω_IR component from the static ${\chi }_0^{\left(2\right)}$, the presence of the THz-driven PhP field opens additional Stokes and anti-Stokes channels described by ${\chi }^{\left(3\right)}$, in which one phonon-polariton is absorbed or emitted and photons are generated: ω_s = 2ω - Ω, and ω_as = 2ω + Ω, where Ω is the PhP frequency, ω is a spectral component of the IR optical frequency, and ω_s and ω_as are the components of the four-wave mixing SHG. These inelastic channels must be included because they govern both the spectral position and the phase-matching condition of the observed SHG fringes. The IR probe and PhP are represented as follows.

(7)$\begin{array}{cc}& E_{{\mathrm{IR}}_{\_}\omega }(t,z)={\epsilon }_{{\mathrm{IR}}_{\_}\omega }\left(z\right)e^{i(k_{\omega }z-\omega t)}\\ & E_{{\mathrm{PhP}}_{\_}\Omega }(t,z)=\mathrm{Re}\{{\epsilon }_{{\mathrm{PhP}}_{\_}\Omega }{e}^{i(-k_{\mathrm{PhP}}z-\Omega t)}+{\epsilon }_{\mathrm{PhP}\_\Omega }^{\ast }{e}^{-i(-k_{\mathrm{PhP}}z-\Omega t)}\}\end{array}$

Second-order polarization source is expressed as:

(8)$P^{\left(2\right)}(t,z)={\epsilon }_0{\chi }^{\left(2\right)}(t,z)E_{{\mathrm{IR}}_{-}\omega }^2(t,z)={\epsilon }_0{\chi }_0^{\left(2\right)}{E}_{{\mathrm{IR}}_{-}\omega }^2(t,z)+{\epsilon }_0{\chi }^{\left(3\right)}{E}_{\mathrm{PhP}}(t,z)E_{{\mathrm{IR}}_{-}\omega }^2(t,z)$

where $E_{\mathrm{IR}\_\omega }^2(t,z)={\epsilon }_{\mathrm{IR}\_\omega }^2\left(z\right)e^{i(2k_{\omega }z-2\omega t)}$. Multiplying by the PhP field yields two cases, the Stokes and anti-Stokes, as follows: $e^{-i(2\omega +\Omega )t}\Rightarrow {\omega }_{\mathrm{as}}=2\omega +\Omega$ (anti-Stokes) and spatial phase factor is 2k_ω - k_PhP; $e^{-i(2\omega -\Omega )t}\Rightarrow {\omega }_{\mathrm{as}}=2\omega -\Omega$ (Stokes) and spatial phase factor is 2k_ω + k_PhP.

We consider phase matching, which involves the counter-propagating PhP and the forward-propagating probe light. For clarity, we explicitly write the Stokes contribution; the corresponding anti-Stokes channel can be treated in an analogous manner and is omitted here and in what follows. The nonlinear polarization associated with ${\chi }^{\left(3\right)}$ is given by the Eq. (9).

(9)$P_s^{\left(2\right)}(t,z)\propto {\epsilon }_0{\chi }^{\left(3\right)}\left(\frac{1}{2}{\epsilon }_{\mathrm{PhP}\_\Omega }^{\ast }\right){\epsilon }_{\omega }^2{e}^{i[(2k_{\omega }+k_{\mathrm{PhP}})z-(2\omega -\Omega \left)t\right]}$

The output field is obtained by integrating the nonlinear polarization along the propagation direction, with the relative time delay between the fields explicitly included. This integration yields:

(10)$E_s^{\text{out}}({\omega }_s,\tau ,t,d)\propto \underset{0}{\overset{d}{\int }}{P}_{3,s}^{\left(2\right)}(t,z,\tau )e^{i{k}_s(d-z)}dz=A_s({\omega }_s,\tau ,t,d)\exp \left[\frac{i\Delta k_{s}d}{2}\right]\frac{\sin \left(\frac{\Delta k_{s}d}{2}\right)}{\Delta k_s}$

where Δk_s = k(ω_s) = n(ω_s)ω_s/c is the wave vector of Stokes-SH, and Δk_s represents the phase mismatch parameter.

(11)$\Delta k_s=2k_{\omega }\left(\omega \right)-k_s(2\omega -\Omega )+k_{\mathrm{PhP}}\left(\Omega \right)$

Here k_PhP is the corresponding PhP wavevector. The Stokes-sideband SHG complex spectral amplitude A_s(ω_s,τ,t,d) at the exit surface is given by:

(12)$A_s({\omega }_s,\tau ,t,d)=e^{-i\Omega \tau }{e}^{i{k}_{s}d}{e}^{-i{\omega }_{s}t}{\epsilon }_0{\chi }_0^{\left(3\right)}{\epsilon }_{\omega }^2{\epsilon }_{\mathrm{PhP}\_\Omega }^{\ast }$

where A_s(ω_s,τ,t,d) is obtained by integrating the third-order nonlinear polarization along the propagation direction under the slowly varying envelope approximation^[23]. Eq. (10) takes into account both phase matching and energy matching, i.e., the conservation of momentum and energy in the nonlinear interaction, from which the participating PhP frequency is derived. The contribution from the intrinsic second-order susceptibility ${\chi }_0^{\left(2\right)}$ to the SHG field can be written as:

(13)$E_0^{\text{out}}(t,d)\propto \underset{0}{\overset{d}{\int }}{P}_0^{\left(2\right)}(t,z)e^{i{k}_s(d-z)}dz=A_{s\_0}\left({\omega }_s\right)\exp \left[\frac{i\Delta k_{s\_0}d}{2}\right]\frac{\sin \left(\frac{i\Delta k_{s\_0}d}{2}\right)}{\Delta k_{s\_0}}$

where $A_{s_{-}0}\left({\omega }_s\right)=2e^{i{k}_{s}d}{e}^{-i{\omega }_{s}t}{\epsilon }_0{\chi }_0^{\left(2\right)}{\epsilon }_{\omega }^2$, Δk_{s_0} = 2k_ω - k_s. The total output field is thus given as follows.

(14)$\begin{array}{cc}& E_{\text{total}}^{\text{out}}=E_s^{\text{out}}({\omega }_s,\tau ,t,d)+E_0^{\text{out}}(t,d)=A_s({\omega }_s,\tau ,t,d)\exp \left[\frac{i\Delta k_{s}d}{2}\right]\frac{\sin \left(\frac{i\Delta k_{s}d}{2}\right)}{\Delta k_s}\\ & +A_{s\_0}\left({\omega }_s\right)\exp \left[\frac{i\Delta k_{s\_0}d}{2}\right]\frac{\sin \left(\frac{i\Delta k_{s\_0}d}{2}\right)}{\Delta k_{s\_0}}\end{array}$

The superscript “out” denotes SHG fields evaluated at the output side of the crystal (after propagation through the ZnO slab). In this expression, $E_0^{\text{out}}$ and $E_s^{\text{out}}$ represent the contributions from the intrinsic ${\chi }_0^{\left(2\right)}$ response and the TFISH contribution mediated by THz or ${\chi }^{\left(3\right)}{\left|{\epsilon }_{\omega }\right|}^2\left|{\epsilon }_{\mathrm{PhP}\_\Omega }\right|$ is the total detected SHG field. The corresponding SHG intensity is as follows.

(15)$\begin{array}{cc}& I({\omega }_s,\tau )={\left|E_0^{\text{out}}\right|}^2+{\left|E_s^{\text{out}}\right|}^2+2\mathrm{Re}\left[E_0^{\text{out}}{E}_s^{\text{out*}}\right]={\left(d{\epsilon }_0{\chi }_0^{\left(2\right)}{\left|{\epsilon }_{\omega }\right|}^2\right)}^2\sin c^2\left(\frac{\Delta k_{s\_0}d}{2}\right)+\\ & \frac{1}{4}{\left(d{\epsilon }_0{\chi }^{\left(3\right)}{\left|{\epsilon }_{\omega }\right|}^2\left|{\epsilon }_{\mathrm{PhP}\_\Omega }\right|\right)}^2\sin c^2\left(\frac{\Delta k_{s}d}{2}\right)\\ & +d^2{\epsilon }_0^2{\chi }_0^{\left(2\right)}{\chi }^{\left(3\right)}{\left|{\epsilon }_{\omega }\right|}^2\left|{\epsilon }_{\mathrm{PhP}\_\Omega }\right|\sin c\left(\frac{\Delta k_{s\_0}d}{2}\right)\sin c\left(\frac{\Delta k_{s}d}{2}\right)\cos (\Omega \tau +\frac{\Delta k_{\mathrm{PhP}}d}{2}+{\phi }_{\mathrm{PhP}\_\Omega })\end{array}$

This procedure yields an explicit analytical expression for the total SHG field intensity as a function of time delay τ and sample thickness d. At this point, we introduce a key simplifying assumption: focusing on the PhP contribution, we assume perfect phase matching for the Stokes channel, i.e., Δk_sd = 0. Using the data obtained from the 1 mm-thick sample, we performed momentum and energy matching among the probe light, SHG, and the PhP according to: ω_SHG = 2ω_IR ± ω_PhP; k_SHG = 2k_IR ± k_PhP and ω_SHG = 2ω_IR ± ω_PhP; k_SHG = 2k_IR ∓ k_PhP. The spectral peaks obtained by fast Fourier transform (FFT) of the experimentally measured time-domain waveforms were found to lie close to the calculated energy–momentum matching curve of the stimulated PhPs. To quantify the deviation, we evaluated the wave-vector offset of all experimental data points from the calculated energy–momentum matching relation and obtained an average mismatch of approximately Δk ≈ 2 cm^-1. This corresponds to a coherence length of about 1.57 cm, which is much larger than the sample thickness of 1 mm. We therefore conclude that the observed frequency-domain peaks are well described by the energy–momentum matching theory. The detailed derivation of the phase-matching conditions (including energy and momentum conservation, forward/backward propagation, and Stokes/anti-Stokes configurations) has been included in the Supplementary materials.

Under this condition, the total SHG intensity at the signal frequency ω_s reduces to $I({\omega }_s,\tau )=I_1+I_2+2\sqrt{I_1{I}_2}\cos (\Omega \tau +\frac{k_{\mathrm{PhP}}d}{2}+{\phi }_{{\mathrm{PhP}}_{-}\Omega })$, which has the standard form of an interference pattern between two coherent SHG waves. Here φ_{PhP_Ω} collects the residual phase offset introduced by the PhP-driven TFISH channel.

When only the influence of the IR light itself is considered, taking thicknesses of 0.3 and 1 mm as examples, the simulated results are shown as Figure 4a,e. Their corresponding interference oscillation curves are shown in Figure 4c,g. The simulation results obtained by incorporating the time-delay-dependent THz-induced SHG contribution into the model are shown in Figure 4b,f, and the corresponding interference curves are presented in Figure 4d,h. To match the photon energy scale of the experiments and calculations, the spacing between adjacent peak energies around 1.71 eV was compared. For a given thickness, the ΔE value at the rear surface is identical with and without the THz pump, and agrees well with the experimental and calculated results (Figure 2i).

Display Full Size

Download

Figure 4. SHG interference in ZnO with and without THz excitation. (a-d) SHG interference in ZnO with and without THz excitation; Calculated E–τ maps and the corresponding energy-dependent SHG intensity curves for a ZnO crystal with thickness d = 0.3 mm, without THz excitation (a, c) and with THz excitation (b, d); (e-h) Calculated E–τ maps and the corresponding energy-dependent SHG intensity curves for a ZnO crystal with thickness d = 1.0 mm, without THz excitation (e, g) and with THz excitation (f, h). In (a, b) and (e, f), the vertical axis denotes the SH photon energy; in (c, d) and (g, h), the vertical axis denotes the total SHG intensity; (i) Fringe spacing ΔE as a function of thickness, as demonstrated by both experiment and theory. SHG: second-harmonic generation; THz: terahertz; SH: second-harmonic.

To account for the tilted interference fringes in the E–τ maps of Figure 4, we generalize the phase difference between the front- and rear-surface SHG contributions to include both the energy-dependent optical phase accumulation and the PhP-induced temporal modulation. The total phase can then be written as

(16)$\Delta \phi (E,\tau )={\Omega }_{\mathrm{PhP}}\tau +\frac{2d}{hc}(n_{2\omega }-n_{\omega })E+{\phi }_0$

where Ω_PhP is the phonon-polariton frequency, and φ₀ is a constant offset. Constructive interference occurs when Δφ(E,τ) = 2πm(m ∈ Z) which yields a family of straight lines in the E–τ plane.

(17)$E_m\left(\tau \right)\approx -\frac{hc{\Omega }_{\mathrm{PhP}}}{2d\left(|n_{2\omega }-n_{\omega }|\right)}\tau +\mathrm{const}$

Thus each bright fringe follows a linear trajectory whose slope is set jointly by Ω_PhP. The observed tilt fringes in the calculating results in Figure 4 confirm that the spectral interference encodes both the PhP frequency and the optical phase accumulation between the two ZnO surfaces governed by the same dispersion term that determines the energy spacing ΔE.

Figure 4a,b show the 0.3 mm sample without and with THz excitation, respectively, while Figure 4e,f show the corresponding data for the 1.0 mm sample. When THz excitation is applied to the sample, it induces an additional effective second-order nonlinear susceptibility, leading to temporal SHG interference and fringe tilting. The simulation results reproduce the experimentally observed fringe patterns, supporting the above interpretation of the interference as a coherent surface-to-surface SHG process. Figure 4d,h further show that, although THz excitation substantially reshapes the overall SHG spectral envelope, the spacing of the interference fringes along the photon-energy axis remains essentially unchanged. This behavior is quantified in Figure 4i, where the simulated values extracted from Figure 4a,b,e,f and the experimental values from Figure 3a both follow the ΔE ∝ 1/d scaling predicted by Eq. (3), and simulations with and without THz excitation collapse onto the same curve. This confirms that the fringe spacing on the photon energy scale is independent of the THz waveform, whereas the THz field primarily controls the fringe contrast and the temporal gating of the SHG signal via TFISH.

In the presence of the THz field, the TFISH and PhP-mediated contributions acquire a delay-dependent phase and amplitude, so that the total phase difference between the front and rear-surface SHG components becomes a function of both photon energy and pump–probe delay Δφ(E,τ). The constructive-interference condition then defines tilted lines in the (E,τ) plane, giving rise to the observed fringe tilting. To verify the contribution of PhPs to the SHG spectra, we performed FFT analysis on the time-domain data of the 1 mm-thick ZnO sample to obtain the tilted interference fringes in the frequency–SH photon-energy map. These were then plotted together with the calculated PhP dispersion curve (white dashed line), obtained based on the method described in the work^[14], as shown in Figure 5. The nearly perfect overlap between the two confirms the crucial role of PhPs.

Display Full Size

Download

Figure 5. Comparison between the calculated dispersion curve and the experimental data. The spectral peaks obtained by FFT of the experimentally measured time-domain waveforms perfectly match the low-frequency branch of the stimulated PhPs. The dashed line represents the calculated dispersion curve; the color map represents the experimentally measured data. FFT: fast Fourier transform; PhPs: phonon polaritons.

According to the interference phase-accumulation model, the main factors affecting the interference are: the initial THz frequency Ω_THz that excites the PhP, the total accumulated phase difference d of the PhP, and the total accumulated phase difference k_s-0d of the intrinsic SHG in the sample. Therefore, even though the THz field can enhance ${\chi }_0^{\left(2\right)}$ and ${\chi }^{\left(3\right)}$ through resonant excitation to some extent^[24], it does not change the temporal oscillation period or the slope of the interference fringes in the frequency domain; it mainly affects their intensity. This finding further verifies the feasibility of modulating the interference pattern from the surface SHG signal by tuning the thickness of the ZnO sample, and provides practical guidance for the design and optimization of nonlinear optical devices. The potential for engineering phase-sensitive THz devices, including ultrafast optical modulators, and high-sensitivity detectors, is demonstrated by the ability to control SHG interference through crystal thickness, IR wavelength, and phase delay.

4. Conclusion

In summary, we have demonstrated that THz fields can actively control coherent SHG at the front and rear surfaces of ZnO crystals, producing well-resolved spectral fringes whose properties are separately governed by THz-driven nonlinearities and optical dispersion. By combining time-resolved TFISH measurements on samples with different thicknesses and spectrally resolved SHG detection, we directly correlated the temporal separation of the front and rear-surface bursts with the fringe pattern in the SHG spectra. The energy spacing ΔE of these fringes follows the scaling ΔE ∝ 1/d and is quantitatively captured by a simple phase-accumulation model that uses the refractive-index difference between the fundamental and SH frequencies, obtained from a Sellmeier-type dispersion relation for ZnO.

The simulations further show that within our model, the TFISH-induced contribution mainly controls the relative amplitudes and temporal separation of the SHG bursts from the two surfaces, thus governing the fringe contrast and the temporal window over which interference is observed. In contrast, the fringe spacing along the energy axis is fixed by the optical dispersion and the propagation distance between the two interfaces. This separation of roles allows the THz-driven phase modulation to be encoded into a static spectral interference pattern: a probe pulse with a duration of only tens of femtoseconds can thus access phase dynamics at frequencies well above the THz carrier without requiring sub-cycle time resolution or high-order harmonic detection. This technology for detecting high-frequency information and ultrafast dynamics is expected to be applicable to a broader class of materials and interfacial systems^[17,25].

The thickness dependence of the fringes highlights the sensitivity of this surface-to-surface SHG interferometry to refractive-index variations under realistic conditions, suggesting a route toward THz–optical hybrid devices for thickness-sensitive sensing, interface-specific diagnostics, and ultrafast surface modulation in polar wide-bandgap materials.

Supplementary materials

The supplementary material for this article is available at: Supplementary materials.

Authors contribution

Tian Y, Song L: Conceptualization, supervision, writing-review & editing.

Zhang X: Investigation, formal analysis, writing-original draft, writing-review & editing.

Ding Y: Investigation, writing-original draft, writing-review & editing.

Sang J: Investigation, writing-review & editing.

Fang Y: Formal analysis, writing-original draft, writing-review & editing.

Hao J: Formal analysis, writing-review & editing.

Conflicts of interest

Ye Tian is an Associate Editor of Light Manipulation and Applications. The other authors declare no conflicts of interest.

Ethical approval

Not applicable.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

Data supporting the findings of this study are available from the Supplementary materials and the corresponding authors upon reasonable request.

Funding

This work was supported by National Key Research and Development Program of China (Grant No. 2022YFA1604400), National Natural Science Foundation of China (Grant Nos 12325409, 12388102, 12304374, U23A6002, 12404400 and 12304477), Shanghai Natural Science Foundation (Grant No. 23ZR1471600), Shanghai Pilot Program for Basic Research, Chinese Academy of Sciences, Shanghai Branch (This program does not have a grant number), CAS Project for Young Scientists in Basic Research (Grant Nos. YSBR-059 and YSBR-127), China Postdoctoral Science Foundation(Grant No. 2025M773338), and Longcheng Talent Plan (Grant No.CQ20230077).

Copyright

© The Author(s) 2026.