Introduction

Phonon polaritons (PhPs)—hybrid quasiparticles arising from strong photon-phonon coupling1—have emerged as a transformative platform for probing nonlinear optics and lattice dynamics in polar crystals2,3,4. Predominantly residing at terahertz (THz) frequencies, these modes enable unprecedented THz field manipulation, offering critical advantages for high-speed photonic modulation5,6. Particularly in the low-THz regime ( <5 THz), PhPs exhibit enhanced spatial delocalization and temporal coherence, amplifying ionic polarization to boost nonlinearity by 2-3 orders of magnitude7,8. Their subwavelength light confinement further drives innovations in ultra-compact photonic circuits and quantum devices2.

While spontaneous PhPs exist intrinsically in polar crystals1,8, coherent excitation requires external field coupling. Current methods face fundamental limitations: impulse stimulated Raman scattering (ISRS) necessitates non-centrosymmetric crystals with concurrent infrared (IR)/Raman activity9,10, difference-frequency generation imposes stringent phase-matching constraints11,12. Recent breakthroughs demonstrate direct THz-pulse excitation through angle-polarization engineering13,14,15, bypassing these restrictions.

The interplay between phonon modes and ionic nonlinear polarizability governs PhP dynamics. THz-driven PhPs uniquely probe spatiotemporal symmetry breaking in nonequilibrium systems—a regime where ionic contributions dominate over electronic nonlinearities7,8,16. Wide-bandgap polar crystals like zinc oxide (ZnO) provide an ideal testbed by: (1) minimizing carrier-induced complications; (2) enhancing phonon-phonon anharmonicity; (3) enabling clean isolation of phonon-polariton (PhP) interactions. Though prior studies mapped ZnO’s electron-phonon coupling17,18,19, PhP generation mechanisms remain underexplored. Furthermore, the optical nonlinearities mediated by PhPs, such as second harmonic generation (SHG), exhibit fundamental distinctions from transverse optical (TO) phonon-driven effects in ferroelectric systems20,21.

In this work, THz pulses applied to a ZnO crystal coherently generate PhPs that interact with optical waves, inducing THz-frequency modulation of optical SHG. Owing to multiple reflections of THz-driven PhPs within the 1-mm-thick crystal, the SHG modulation persists for >90 ps—a duration directly governed by polariton propagation dynamics. Narrowband optical probing conclusively confirms THz-wave excitation of PhPs, as the limited optical bandwidth precludes direct photoexcitation of these quasiparticles. By employing broadband optical pulses interacting with co- and counter-propagating PhPs, we demonstrate phase-matched interactions between polaritons and optical fields. This nonlinear coupling generates broadband SHG signals exhibiting spectral oscillations synchronized with the characteristic frequency of phase-matched PhPs, thereby establishing a direct probe of polaritonic wavevector-selective dynamics.

Results

The experimental configuration in Fig. 1 enables spectral-temporal tracking of PhP dynamics through THz-IR synergistic spectroscopy. Collinearly focused THz (0.1–6 THz spectrum, 100 kV/cm peak field) and tunable IR pulses (50 fs duration, 1400–2100 nm broadband) illuminate the <0001 > -oriented 1-mm-thick ZnO crystal, with time-resolved SHG spectroscopy executed via precision delay scanning between pulses. THz generation employs an optical parametric amplifier (OPA)-driven DSTMS crystal (Fig. 1b)22, while IR pulses from the same OPA system (Fig. 1a) drive the SHG process (see Methods). Electromagnetic excitation induces non-zero dipole oscillations, enabling coherent phonon detection via SHG modulations13. By analyzing the SHG signals from 1450 nm-centered infrared probe pulses, we observe sustained oscillatory dynamics persisting over tens of picosecond timescales (Fig. 1c). Zooming into a representative temporal window reveals THz-rate nonlinear optical modulation at 3.6 THz frequency with 18 dB extinction ratio—performance metrics surpassing conventional lithium niobate (LN) electro-optic modulators in both bandwidth and switching contrast23. This THz-mediated modulation mechanism demonstrates unprecedented photonic control capabilities through engineered phonon-photon interactions.

Fig. 1: Collinear experimental geometry for THz-phonon polariton interaction studies.
figure 1

THz pulses (green beam, 0.1–6 THz spectrum) and tunable infrared (IR) pulses (red beam, 100 fs bandwidth) are incident at normal incidence on the <0001> surface of a 1-mm-thick ZnO crystal. Time-delay-scanned interferometric detection of second harmonic generation (SHG) signals using a spectrometer-coupled CCD. Insets: Spectral characteristics of excitation sources a broadband/narrowband IR profiles in red shadow/colored curves; b THz spectral amplitude in green shaded curve. c Time-delay-dependent SH intensity at 725 nm central wavelength demonstrates 3.6-THz modulation bandwidth with 18-dB extinction ratio.

Full size image

To investigate phonon excitation dynamics and energy-scale evolution, we employed narrowband IR pulses (centered at 1400, 1450, and 1490 nm) and THz to drive SHG, as illustrated by the spectral curves in Fig. 1a. Figure 2a demonstrates the comparative SHG spectra with (solid lines) and without (colored shadows) THz pulse assistance, revealing significant enhancement of SHG intensity under THz irradiation. Notably, a series of discrete peaks emerge in the THz-assisted spectra, exhibiting precise correspondence with the LO, E2L, and E1(TO) phonon modes. Temporal evolution of the phonon modes was monitored through time-resolved SHG measurements (Fig. 2b). The synchronized excitation of IR and THz pulses at the ZnO crystal surface (t = 0 ps) initiates a THz-assisted Hyper-Raman process that coherently activates multiple phonon modes. This dual-beam synergistic approach enables comprehensive Raman spectroscopic signatures that are distinctly different from single-beam excitation phenomena24, establishing a methodology for spectral-temporal tracking of phonon mode evolution dynamics with unprecedented temporal resolution.

Fig. 2: THz-phonon coupling dynamics probed by wavelength-resolved SHG spectroscopy under a narrowband probe light set-up.
figure 2

a THz-induced SHG spectral modification: Solid curves show SHG intensities under THz excitation (0 ps delay) versus reference spectra without THz (color-shaded regions). Dashed vertical lines mark TO/LO phonon mode energies. b Time-delay-dependent SHG spectral evolution, revealing oscillatory features persisting beyond 60 ps. c Wavelength-selective SHG oscillation traces extracted from b, demonstrating pump-energy-dependent modulation depths. d Fourier transform analysis of oscillations SH in c, identifying dominant frequency components at 3.3, 3.6 and 3.9 THz, with shaded regions indicating ±0.15 THz spectral resolution.

Full size image

The LO phonon modes exhibit a characteristic lifetime of 2.8 ps, consistent with the lifetimes reported in previous studies25,26. Subsequent temporal evolution shows exclusive persistence of the E2 mode beyond 4 ps, confirming its long-lived nature19. Figure 2c and d illustrates the sustained intensity oscillations in SHG signals and their corresponding frequency spectrum, respectively. These oscillations offer quantitative insights into the coherent phonon dynamics within the THz frequency range. The frequency centers are located at 3.3, 3.6, and 3.9 THz, respectively. This modulation frequency demonstrates that short-wave infrared corresponds to high frequencies, while long-wave infrared corresponds to low frequencies, and this modulation behavior persists over an extended period. Although the E2 mode also exhibits an exceptionally long lifetime, this modulation resembles a wavelength-dependent matching relationship, significantly distinguishing the operational principles and nature of the E2 mode. Given its characteristics, it is necessary to validate this matching relationship using a broadband probe light with a wide spectral range.

When broadband IR pulses are applied (red shadows in the inset of Fig. 1a), multiple phonon modes are coherently excited through the THz-assisted Hyper-Raman process upon synchronization of THz and IR pulses at the ZnO crystal surface, consistent with the narrowband experimental observations. The THz wave assistance facilitates coherent excitation of LO, TO, and E2 phonons on the crystal surface, generating pronounced anharmonic phonon responses across a broad frequency range (Fig. 3a). To improve the detection of low-frequency phonon modes and to determine whether PhPs are generated, a p-polarized slit aligned with the in-plane E1(TO) vibrational mode - sharing the same plane as the E2L vibration - is positioned before the spectrometer, with the collected data presented in Fig. 3b. This configuration selectively detects the E2L mode and PhPs while excluding the out-of-plane A1(TO) mode27. Temporal analysis reveals the E2L phonon (blue arrows) excitation within 1 ps, followed by the emergence of a distinct mode (orange arrows) between 3-4 THz at approximately 3 ps, persisting for tens of picoseconds. This mapping methodology deciphers fundamental spectral-temporal correlations through geometrically distinct signatures: vertical features (e.g., at 1 and 3 ps in Fig. 3b) encode intrinsic phonon mode dynamics, while the oblique trajectory emerging beyond 4 ps and persisting over 40 ps (Fig. 3b) directly visualizes propagating PhP wavepackets. Crucially, these dual geometric markers establish a unified analytic framework for real-time tracking of both local phonon excitations and nonlocal polariton-mediated energy transfer processes at the nanoscale, bridging atomic-scale vibrations to macroscopic photonic responses through their spatiotemporal interplay14.

Fig. 3: Time-resolved FFT spectra of SHG and schematic illustration of phonon and phonon-polariton (PhP) generation and detection mechanisms.
figure 3

a FFT spectra showing distinct phonon modes (indicated by arrows) excited at the temporal coincidence of THz and IR pulses on the crystal surface. b Polarization-resolved FFT spectra obtained using a p-polarization slit positioned before the spectrometer, with colored arrows indicating the spectral positions of E2L (blue) and PhP (orange) modes. c Schematic diagram illustrating the experimental configuration for PhP generation and detection.

Full size image

Figure 3c illustrates the generation and detection mechanism of PhPs, where THz waves (green arrows) initially excite PhPs (brown waveform) that subsequently couple with the IR field. The propagating PhPs then interact with the incident IR pulse (red waveform/arrows), producing a SHG signal at \({\omega }_{{{\mbox{SHG}}}{{{\_}}}{{\mbox{PhP}}}}\) via ionic displacement contributions. Concurrently, the inherent polarization modulation induced by PhP generating an additional SHG signal at \({\omega }_{{{\mbox{SHG}}}{{{\_}}}{{{{\boldsymbol{\chi }}}}}^{\left(2\right)}}\). Interference between these two phase-correlated SHG components—one arising from optical nonlinearity and the other from polariton-mediated ionic oscillations—drives the characteristic intensity oscillations in the detected SHG signal. In the inset of Fig. 3c, we demonstrate how the SHG oscillation period in the time-___domain experimental data reflects the PhP frequency. To qualitatively describe the TFISH effect and PhP-IR light interaction, we derive the effective second-order nonlinear susceptibility \({{{{\boldsymbol{\chi }}}}}_{{eff}-{THz}}^{\left(2\right)}\) and \({{{{\boldsymbol{\chi }}}}}_{{eff}-{PhP}}^{\left(2\right)}\):

$${{{{\boldsymbol{\chi }}}}}_{{eff}-{THz}}^{\left(2\right)}={{{{\boldsymbol{\chi }}}}}_{0}^{\left(2\right)}+{{{{\boldsymbol{\chi }}}}}_{0}^{\left(3\right)}\cdot {{{{\boldsymbol{E}}}}}_{{THz}}$$
(1)
$${{{{\boldsymbol{\chi }}}}}_{{eff}-{PhP}}^{\left(2\right)}= {{{{\boldsymbol{\chi }}}}}_{0}^{\left(2\right)}+\frac{1}{{\epsilon }_{0}}{{{{\boldsymbol{\chi }}}}}_{0}^{\left(3\right)}\cdot {{{{\boldsymbol{P}}}}}_{{PhP}}+\frac{\partial {{{{\boldsymbol{\chi }}}}}_{0}^{\left(2\right)}}{\partial {{{\boldsymbol{Q}}}}}\cdot {{{\boldsymbol{Q}}}}\\= {{{{\boldsymbol{\chi }}}}}_{0}^{\left(2\right)}+\left(\frac{1}{{\epsilon }_{0}}{{{{\boldsymbol{\chi }}}}}_{0}^{\left(2\right)}{\omega }_{0}\sqrt{{\epsilon }_{0}\left({\varepsilon }_{0}-{\varepsilon }_{\infty }\right)}+\frac{\partial {{{{\boldsymbol{\chi }}}}}_{0}^{\left(2\right)}}{\partial {{{\boldsymbol{Q}}}}}\right)\cdot {{{\boldsymbol{Q}}}}$$
(2)

Equation (1) governs the surface-localized THz electric field-induced TFISH effect, whereas Eq. (2) describes the bulk-mediated SHG process driven by PhPs. In these formulations: \({{{{\boldsymbol{\chi }}}}}_{0}^{\left(2\right)}\) is the intrinsic second-order susceptibility of ZnO crystal, \({{{{\boldsymbol{\chi }}}}}_{0}^{\left(3\right)}\cdot {{{{\boldsymbol{E}}}}}_{{THz}}\) is TFISH effect term, \({{{{\boldsymbol{\chi }}}}}_{0}^{\left(3\right)}\cdot {{{{\boldsymbol{P}}}}}_{{PhP}}/{\epsilon }_{0}\) describes the participation of PhPs in four-wave mixing process, and \(\left(\partial {{{{\boldsymbol{\chi }}}}}_{0}^{\left(2\right)}/\partial {{{\boldsymbol{Q}}}}\right)\cdot {{{\boldsymbol{Q}}}}\) is the coherent Hyper-Raman scattering term. \({{{\boldsymbol{Q}}}}\) represents the amplitude of the ionic displacements, \({{{{\boldsymbol{P}}}}}_{{PhP}}={\omega }_{0}\sqrt{{\epsilon }_{0}\left({\varepsilon }_{0}-{\varepsilon }_{\infty }\right)}{{{\boldsymbol{Q}}}}\)1,2. \({\omega }_{0}\) represents the circular frequency of E1(TO) phonon mode in ZnO, \({\epsilon }_{0}\) is the vacuum permittivity, and \({\varepsilon }_{0}/{\varepsilon }_{\infty }\) donates the static/high-frequency dielectric constant.

Notably, while the Hyper-Raman term originates from spontaneous PhP behavior, the incident THz field synchronizes initially incoherent ionic oscillations16, enabling SHG modulation through a four-wave-mixing-like mechanism. Both processes strictly obey energy and momentum conservation laws.

We establish four distinct phase matching configurations based on the relative propagation directions of PhPs and IR pulses. Using scalar representations for the IR electric field \({E}_{{{\mbox{IR}}}}\) and the PhP field \({Q}_{{{\mbox{PhP}}}}\) to clarify these conditions, considering both forward (\(+z\)) and backward (\(-z\)) propagation scenarios. For the forward-propagating infrared electric field \({E}_{{{\mbox{IR}}}}={E}_{0}exp \left(-i{\omega }_{{{\mbox{IR}}}}t+i{k}_{{{\mbox{IR}}}}z\right)\), the coupling between the backward-propagating PhP field QPhP- and its complex conjugate \({Q}_{{{\rm{PhP}}-}}^{ * }\) results in:

$$\begin{array}{c}{P}_{{{\mbox{bwd}}}}^{{{\mbox{as}}}}\propto {E}_{{{\mbox{IR}}}}{E}_{{{\mbox{IR}}}}{Q}_{{{\mbox{PhP}}}-}\propto exp \left[-i\left(2{\omega }_{{{\mbox{IR}}}}+{\omega }_{{{\mbox{PhP}}}}\right)t+i\left(2{k}_{{{\mbox{IR}}}}-{k}_{{{\mbox{PhP}}}}\right)z\right]\\ {P}_{{{\mbox{bwd}}}}^{s}\propto {E}_{{{\mbox{IR}}}}{E}_{{{\mbox{IR}}}}{Q}_{{{\mbox{PhP}}}-}^{ * }\propto exp \left[-i\left(2{\omega }_{{{\mbox{IR}}}}-{\omega }_{{{\mbox{PhP}}}}\right)t+i\left(2{k}_{{{\mbox{IR}}}}+{k}_{{{\mbox{PhP}}}}\right)z\right]\end{array}$$
(3)

While the conservation of momentum and energy is considered as:

$$\begin{array}{c}{\omega }_{{{\mbox{SHG}}}}=2{\omega }_{{{\mbox{IR}}}}\pm {\omega }_{{{\mbox{PhP}}}}\\ {k}_{{{\mbox{SHG}}}}=2{k}_{{{\mbox{IR}}}}\,\mp\, {k}_{{{\mbox{PhP}}}}\end{array}$$
(4)

For the forward-propagating PhP, the conservation of momentum and energy is considered as:

$$\begin{array}{c}{\omega }_{{{\mbox{SHG}}}}=2{\omega }_{{{\mbox{IR}}}}\pm {\omega }_{{{\mbox{PhP}}}}\\ {k}_{{{\mbox{SHG}}}}=2{k}_{{{\mbox{IR}}}}\pm {k}_{{{\mbox{PhP}}}}\end{array}$$
(5)

The oscillatory signal originates from the interference between the inherent second-order susceptibility \({{{{\boldsymbol{\chi }}}}}_{0}^{\left(2\right)}\) induced SHG and the PhP-induced SHG. Here, we consider two frequency components of fundamental IR spectrum participant in stokes-like backward-propagating PhP scenario, \({\omega }_{0}\) with phase \({\varphi }_{0}\) generate \({\omega }_{{{\mbox{SHG}}}{{{\_}}}{{\mbox{PhP}}}}=2{\omega }_{0}-{\omega }_{{{\mbox{PhP}}}}\) through \({{{{\boldsymbol{\chi }}}}}_{{{\mbox{PhP}}}}^{\left(2\right)}\) and \({\omega }_{0}-{\omega }_{{{\mbox{PhP}}}}/2\) with phase \({\varphi }_{0}+\varDelta \varphi\) generate \({\omega }_{{{\mbox{SHG}}}{{{\rm{\_}}}}0}=2{\omega }_{0}-{\omega }_{{{\mbox{PhP}}}}\) through \({{{{\boldsymbol{\chi }}}}}_{0}^{\left(2\right)}\). The phases of these two SHG fields can be expressed as:

$$\begin{array}{c}{\varphi }_{{{\mbox{SHG}}}\_{0}}=2\left({\varphi }_{0}+\varDelta \varphi \right)=2{\varphi }_{0}+2\varDelta \varphi \\ {\varphi }_{{{\mbox{SHG}}}\_{{\mbox{PhP}}}}=2{\varphi }_{0}+{\varphi }_{{{\mbox{PhP}}}}\left(\tau \right)=2{\varphi }_{0}+{\omega }_{{{\mbox{PhP}}}}\tau+{\varphi }_{{{\mbox{PhP}}}}\end{array}$$
(6)

Here we consider \({\omega }_{{{\rm{SHG}}}\_{0}}\) and \({\omega }_{{{\rm{SHG}}}\_{{\rm{PhP}}}}\) are \({{{{\boldsymbol{\chi }}}}}_{0}^{(2)}\)-induced and \({{{{\boldsymbol{\chi }}}}}_{{{\rm{PhP}}}}^{(2)}\)-induced SHG signals respectively, interference between the two SHG components requires strict frequency matching \({\omega }_{{{\rm{SHG}}}\_{0}}={\omega }_{{{\rm{SHG}}}\_{{\rm{PhP}}}}\). The phase difference between these two fields leads to an interference term that varies with time delay \(\tau\):

$${I}_{{{\mbox{oscillation}}}}\left(\tau \right) \propto {\left|{E}_{{{\mbox{SHG}}}\_{0}}+{E}_{{{\mbox{SHG}}}\_{{\rm{PhP}}}}\right|}^{2}-{\left|{E}_{{{\mbox{SHG}}}{{{\_}}}0}\right|}^{2}-{\left|{E}_{{{\mbox{SHG}}}\_{{\rm{PhP}}}}\right|}^{2} \\ \propto cos \left({\varphi }_{{{\mbox{SHG}}}{{{\_}}}{{\mbox{PhP}}}}-{\varphi }_{{{\mbox{SHG}}}\_{0}}\right)=cos \left({\omega }_{{{\mbox{PhP}}}}\tau+{\varphi }_{0}^{{\prime} }\right)$$
(7)

\({\omega }_{{{\mbox{PhP}}}}\tau\) quantifies the temporal modulation of the interference signal induced by PhPs, corresponding to the oscillation frequency observed in the SHG signal, while φPhP governs the initial phase of the PhP mode that determines the interference phase offset. The cosinusoidal dependence explicitly demonstrates periodic modulation of the interference term, directly manifesting the coherent coupling between the two SHG components and generating time-dependent oscillations in the SHG intensity as a function of pump-probe delay. Systematic expansion of the phase terms enables mechanistic interpretation of distinct physical phenomena, including time-delay-induced phase accumulation and PhP-mediated nonlinear contributions (Fig. 3c). Ultimately, the measured interference pattern encodes both the mutual coherence of the participating SHG waves and the phase-resolved characteristics of their interaction dynamics, providing critical insights into the energy transfer pathways and phase relationships within the system. We have therefore developed a methodology for sustaining modulated SHG waveform outputs, where temporal persistence and spectral characteristics obey the following governing relationship. The SH intensity obeys as follows,

$${I}_{{{\mbox{SHG}}}}\left(N,\tau \right)\sim {I}_{0}+{I}_{1}{R}^{{\,N}}exp \left(-\frac{\tau -N{t}_{0}}{{\tau }_{0}}\right)+2\sqrt{{I}_{0}{I}_{1}{R}^{N}}exp \left(-\frac{\tau -N{t}_{0}}{2{\tau }_{0}}\right)cos \left({\omega }_{{{\mbox{PhP}}}}\tau \right)$$
(8)

while \(N\) is the reflection times of PhP, \(R\) is the reflectivity of PhP, \({\tau }_{0}\) is the lifetime of E1(TO) phonon, \({t}_{0}=d/{\upsilon }_{{{\mbox{PhP}}}}\) is the time of PhP propragating from front to rear surface, is the time delay, \({I}_{0}\) is the intensity of SHG filed induced by \({{{{\boldsymbol{\chi }}}}}_{0}^{\left(2\right)}\), \({I}_{1}{R}^{N}exp \left[-\left(\tau -N{t}_{0}/{\tau }_{0}\right)\right]\) is the SHG intensity induced by PhP. We can caculate the durtion of SHG signal through the relation

$$\tau \sim N{t}_{0}+2{\tau }_{0}{ln}\left(2\right)-{\tau }_{0}{ln}\left(\frac{{I}_{0}}{{I}_{1}{R}^{N}}\right)\sim N{t}_{0}+2{\tau }_{0}{ln}\left(2\right)$$
(9)

while \({I}_{0}\sim {I}_{1}{R}^{N}\). For the 1 mm thick ZnO crystal, \({t}_{0}\sim 10\) ps, \({\tau }_{0}\sim 3\) ps, \(R\sim 0.27\), \({I}_{1}/{I}_{0}\sim 500\), so we can get \(N\sim 5\), at this time \(\tau \sim 50-60\) ps. Our experimental realization achieves sustained modulated waveforms with exceptional ~90 ps coherence duration, outperforming fundamental-frequency modulation schemes (e.g., GHz-range electro-optic modulators23,28) by two orders of magnitude in temporal span. This phase-engineered approach uniquely enhances modulation depth at second-harmonic frequencies through tailored nonlinear optical interactions. Scalability analysis indicates substantial bandwidth expansion potential via crystal orientation optimization and THz field enhancement, effectively transforming ZnO—a prototypical wide-bandgap semiconductor—into a high-performance platform for THz waveform synthesis and control. This crystalline reengineering paradigm establishes transformative design principles for ultra-broadband photonic modulators, overcoming intrinsic limitations of conventional material systems through symmetry-engineered nonlinear optics.

Through combined experimental and theoretical analysis, we elucidate the physical origin of long-lived PhPs in ZnO, which arise when intersecting dispersion relations between optical fields and material polar excitations achieve simultaneous energy-momentum matching. This strong coupling induces light-matter hybridization, forming PhPs as coherent superposition states rather than independent optical or vibrational modes. To validate this mechanism, we conduct full-wave simulations incorporating both momentum and energy conservation criteria. Figure 4a presents time-resolved SHG data at t = 50 ps under THz-pump/broadband 1450 nm probe configuration with 3 ps temporal window, where the oscillatory features are exclusively governed by PhP dynamics. The vertical IR incidence on the ZnO (0001) surface inherently satisfies the scattering phase-matching condition, while the dielectric function analysis focuses on the ordinary optical axis \({\varepsilon }_{\perp }\) for normal incidence. Through E-Q coupling formalism applied to Hyper-Raman and four-wave mixing processes, the refractive index and the propagation vector \(k\left(\omega \right)=n\left(\omega \right)\cdot \omega /c\) are substituted into the four energy and momentum conservation conditions of the PhPs dispersion relationship \(\omega \left(k\right)\), the momentum of SHG can be derived as follows,

$${k}_{{{\rm{SHG}}}}=n\left(2{\omega }_{{{\mbox{IR}}}}-{\omega }_{{{\mbox{PhP}}}}\right)\cdot \left(2{\omega }_{{{\mbox{IR}}}}-{\omega }_{{{\mbox{PhP}}}}\right)$$
(10)
Fig. 4: Phase-matching characteristics of PhPs in ZnO.
figure 4

a Dispersion analysis showing theoretical PhP branches (red and purple curves) superimposed with experimental data (colored spheres). The blue and green demarcation lines indicate momentum-energy matching boundaries under Stokes scattering conditions b Numerically simulated PhP branches from momentum-energy conservation theory, with the experimentally accessible spectral range highlighted by the gray-shaded region. c Experimental verification of PhP-optical wave phase-matching using broadband probes, demonstrating spectral alignment between observed SHG signals (markers) and theoretical predictions (solid curves).

Full size image

Simulations utilizing refractive index29 show excellent agreement between experimental data (colored circles) and phase-matching trajectories (inclined lines) in Fig. 4a. Intersections of blue/green theoretical curves with the PhP dispersion mark strong coupling regions, while the current experimental configuration prohibits upper branch access due to insufficient THz spectral coverage for Stokes resonance (requiring \({\omega }_{{{\mbox{SHG}}}}=2{\omega }_{{{\mbox{IR}}}}-{\omega }_{{{\mbox{PhP}}}}\)). Consequently, PhP measurements are exclusively recorded along the lower dispersion branch. Figure 4b numerically reconstructs the SHG spectral distribution through momentum-energy matched PhP combinations (PA, PB and PA + PB), while Fig. 4c cross-validates these predictions using dual experimental configurations (THz-1450 nm/THz-1750 nm IR), demonstrating quantitative consistency between simulated and observed PhP signatures. The experimental PB-related spectral feature exhibits minor deviations from theoretical predictions, primarily attributed to its sensitivity in the high-wavevector (high-k) regime where minor dielectric variations (particularly in the \({\varepsilon }_{0}/{\varepsilon }_{\infty }\) ratio) amplify spectral discrepancies, and further complicated by unresolved coupling between E1(TO) and A1(TO) vibrational modes. This investigation establishes a fundamental framework for polariton manipulation through THz excitation based on the demonstrated phase-matching principles in polar crystals. The methodology proves universally applicable beyond ZnO systems, with analogous polar crystals exhibiting comparable phase-matching trajectories and excitation protocols. The conceptual advancement extends significantly beyond conventional second-harmonic analysis—our spectral mapping reveals that higher-order nonlinear phenomena (including n-th harmonic generation30 and multi-wave mixing processes) maintain viable phase-matching conditions under appropriate multi-chromatic excitation schemes. This design principle opens transformative possibilities for developing ultra-broadband photonic devices through crystallographic engineering, particularly enabling advanced architectures for high-speed optical signal processing and quantum coherent control applications.

In summary, our findings establish a complete framework for PhP formation and nonlinear interaction dynamics in ZnO through time-resolved THz-SHG spectroscopy. Real-time tracking of narrowband IR frequency shifts reveals THz-phonon coupling sustaining SHG signals beyond 90 ps, while variation of SHG oscillation frequencies with broadband probes directly maps phase-matching conditions. These results not only decode PhP formation mechanisms but also demonstrate their exceptional coherence properties for enhanced nonlinear optical processes. The achieved fundamental understanding enables targeted phonon engineering through dispersion relation manipulation and opens avenues for designing THz-band nonlinear devices exploiting PhP-mediated light control. This work bridges the gap between material spectroscopy and practical photonic applications, providing both microscopic insights into polariton dynamics and macroscopic guidelines for developing next-generation THz photonic devices.

Methods

Experimental setup

The experimental configuration employs a Ti:sapphire amplifier generating 30 fs laser pulses at 800 nm central wavelength with 7 mJ pulse energy and 1 kHz repetition rate. This drives an optical parametric amplifier (OPA), producing tunable infrared (IR) pulses (1200–2100 nm spectral range, ~50 fs duration, ~1.3 mJ pulse energy). Approximately 90% of the OPA signal output pumps a 0.4-mm-thick DSTMS organic crystal, generating single-cycle terahertz (THz) pulses spanning 0.1–6 THz bandwidth with 2 μJ pulse energy22. The remaining 10% of signal pulses and the idler pulses are temporally synchronized with THz pulses to perform time-resolved second harmonic generation (SHG) spectroscopy measurements (broadband IR detection setup). The second harmonic of IR pulses is subsequently detected and recorded by a spectrometer with 0.2 nm spectral resolution. To implement the narrowband IR detection setup, we further insert a 12 nm bandwidth bandpass filter into the IR path, using the filtered output as the probe beam for narrowband SHG measurements.

Material parameters

The ZnO crystal used in our experiments has a hexagonal unit cell with lattice parameters: a = 3.2500 Å and c = 5.2060 Å. The key phonon modes and corresponding frequency in hexagonal ZnO: E2L(99 cm−1), E2H (438 cm−1), A1(TO) (375 cm−1), E1(TO) (410 cm−1), A1(LO) (575 cm−1), E1(LO) (590 cm−1).