Pith. sign in

REVIEW 3 major objections 5 minor 61 references

A THz pump–Raman probe technique maps the momentum-resolved frequency and damping of phonon-polaritons in LiNbO3 and traces a step-like drop in the intrinsic phonon damping to anharmonic decay.

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

T0 review · deepseek-v4-flash

2026-08-01 20:08 UTC pith:D2SF5ZCO

load-bearing objection A genuinely new THz pump–Raman probe scheme with a solid dispersion measurement, plus an interesting but under-supported damping step that needs error bars before it becomes a result. the 3 major comments →

arxiv 2607.16718 v1 pith:D2SF5ZCO submitted 2026-07-18 cond-mat.mtrl-sci physics.optics

Measuring momentum-resolved dissipation of phonon-polaritons in LiNbO₃ with terahertz driving

classification cond-mat.mtrl-sci physics.optics
keywords phonon-polaritonsTHz pump–Raman probepolariton dispersionLiNbO3phonon dampinganharmonic couplingphase-matchingthree-wave mixing
verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

This paper establishes a new spectroscopic route, THz pump–Raman probe (TP-RP), for mapping the full complex dispersion of phonon-polaritons—both frequency and damping as functions of momentum—in non-centrosymmetric polar crystals, demonstrated on the lowest E-symmetry mode of lithium niobate. The key idea is that a broadband THz pump launches polaritons that propagate forward or backward through a thick crystal, while a wavelength-tunable near-infrared Raman probe enforces a phase-matching condition that selects specific momenta; the forward and backward signals are separated in time, giving cleaner linewidths than collinear all-optical schemes. Combining the spectra with a many-body theory that includes propagation effects and uses the bare-phonon damping as the only free parameter, the authors find that no constant damping fits the data. Instead, the intrinsic phonon damping falls from about 3.15 THz to about 0.96 THz in a smooth step centered near 2.8 THz, which they interpret as a crossover between decay into a low-frequency continuum—possibly acoustic phonons—and suppression of that channel. If correct, the result shows that momentum-resolved nonlinear THz experiments can access intrinsic phonon dissipation in a regime invisible to conventional linear optics.

Core claim

On its own terms, the paper claims that TP-RP can reconstruct both the real and imaginary parts of the E(TO1) phonon-polariton dispersion in LiNbO3 with better accuracy than earlier all-optical ISRS methods, because the incident THz pump and its rear-interface echo produce temporally separated forward- and backward-propagating polariton signals. From the phase-matched peak frequencies and linewidths, with momenta assigned by Eqs. (1) and (2), the real part of the dispersion matches the standard Lorentz-oscillator prediction, while the imaginary part deviates from any constant-damping calculation. Fitting the full theoretical spectra—Eqs. (5)–(7) with the bare-phonon damping γ as the sole fre

What carries the argument

The load-bearing elements are the phase-matching equations for three-wave mixing, k+ = −(ne/c)ωpr + (no/c)(ωpr + Ω+) and k− = (ne/c)ωpr − (no/c)(ωpr − Ω−), which map each measured peak frequency to a wavevector; the many-body second-order current whose interaction kernel K(2)(ω′,Ω) ∝ Z*R/(Ω2 + 2iγΩ − ωTO2) contains the bare-phonon propagator and the single free parameter γ; and a generalized Maxwell–Fresnel treatment that propagates pump and probe through the sample, including transmission, reflection, and Fabry–Pérot effects. The temporal separation of Signal 1 (forward polaritons) and Signal 2 (backward polaritons) is what makes the linewidth extraction cleaner, since the two branches do n

Load-bearing premise

The central damping result rests on the assumption that the measured FFT peak linewidths are described by the model of Eqs. (5)–(7) with the bare-phonon damping γ as the only free parameter; if any unmodeled, wavelength-dependent broadening (finite probe bandwidth, FFT window choice, sample inhomogeneity, or residual Signal-1 leakage into Signal-2 windows) contributes, the extracted step in γ(Ω) would not be an intrinsic property of the phonon.

What would settle it

Compare TP-RP spectra taken with different FFT time windows and sample thicknesses: if the step in γ(Ω) is intrinsic, the extracted γ values must be independent of window choice and thickness; the current data already hint at fragility, since Signal-2 peak frequencies for 1250 and 1300 nm shift by roughly 0.3–0.5 THz between full and reduced windows. A direct test would measure the E(TO1) polariton damping in the 1–2 THz range using narrowband THz transmission or time-domain spectroscopy; if no step near 2.8 THz appears in the intrinsic damping, the step is an artifact of the TP-RP linewidth m

Watch this falsifier — get emailed when new claim-graph text bears on it.

If this is right

  • TP-RP yields both the real and imaginary parts of the polariton dispersion in one experiment, with forward and backward signals separated in time so their linewidths are not blended as in collinear ISRS.
  • The extracted intrinsic phonon damping γ(Ω) shows a smooth step from roughly 3.15 THz to 0.96 THz centered near 2.8 THz, implying the phonon has a decay channel at low frequencies that switches off above the step.
  • The step location matches zone-boundary acoustic-phonon frequencies in LiNbO3 (2.5–3 THz), making anharmonic phonon–phonon coupling the paper's proposed microscopic mechanism.
  • Because the measurement region lies far from the TO resonance, TP-RP can access damping behavior that linear-response spectroscopies cannot see.
  • The technique can be extended to other non-centrosymmetric polar crystals and to higher E(TO2)/E(TO3) branches by shaping the THz pump spectrum.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • If the step in γ(Ω) is intrinsic, a temperature-dependent TP-RP study should be revealing: the step's position and height would trace the thermal occupation of the acoustic-phonon bath, whereas an artifact would be temperature-independent.
  • The same phase-matched formalism could be turned around: rather than fixing the probe wavelength to extract a single (Ω, k) point, a chirped broadband probe could reconstruct an arc of the dispersion in one shot, effectively imaging the complex dielectric function.
  • The anomaly around k ≈ 6000 cm−1 suggests a second damping channel near 3.2–3.8 THz; extending TP-RP to probe wavelengths between 800 and 1200 nm would test whether γ(Ω) rises again.
  • Should the method transfer to materials with weaker damping, the temporal separation of forward and backward signals could enable direct measurements of polariton group velocity and propagation losses rather than relying on fits.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit.

Referee Report

3 major / 5 minor

Summary. The paper introduces THz pump–Raman probe (TP-RP) as a method to map the momentum-resolved dispersion and damping of phonon-polaritons in non-centrosymmetric crystals, demonstrated for the E(TO1) mode in LiNbO3. A broadband THz pump resonantly drives forward- and backward-propagating polaritons, and a tunable NIR probe provides phase-matched detection via three-wave mixing. The real part of the dispersion is reconstructed from peak frequencies and phase-matching equations, and agrees with a Lorentz-oscillator model using literature constants. The imaginary part is obtained from FFT peak widths and a many-body propagation model with the bare phonon damping γ as a free parameter per spectrum. The extracted γ(Ω) shows a step-like decrease from ~3.15 THz to ~0.96 THz near 2.8 THz, attributed to coupling to acoustic phonons. The paper claims TP-RP enables accurate extraction of both frequency and damping with reduced uncertainty compared to collinear ISRS.

Significance. If the central claims hold, the technique is a genuine advance: it provides direct THz excitation and phase-matched Raman detection with temporal separation of forward and backward signals, and the theoretical framework, Eqs. (5)–(7) with the Maxwell–Fresnel propagation treatment of Appendix G, is substantial and reproduces many experimental features. The real-dispersion reconstruction is checked against fixed literature constants and is not circular; the reported frequency points are consistent with the Lorentz-oscillator dispersion. However, the damping result is load-bearing for the paper's main novel conclusion, and the manuscript does not currently secure it. The fitted γ values carry no reported uncertainties, the manual choice of FFT windows produces large shifts in the low-frequency Signal-2 peaks, and the final 'reproduction' of polariton damping using the fitted γ(Ω) is not an independent validation. The stress-test concern therefore lands: unmodeled broadening or windowing artifacts could masquerade as the step in γ(Ω). The paper would be publishable after a major revision that quantifies these uncertainties and tests the robustness of the step.

major comments (3)
  1. [Appendix H / Fig. 4c] The extracted γ values are presented without error bars, and the tanh parameters (γ0, γ∞, Ω0, Δ) are fitted to those values without propagating any uncertainty. Since γ is the only free parameter in Eqs. (5)–(7), every unmodeled broadening mechanism (probe bandwidth, FFT windowing, sample inhomogeneity, residual Signal-1 leakage) is absorbed into γ. This is not a cosmetic issue: Tables E5 and E6 show that for Signal 2 at 1250 nm the peak frequency shifts from 1.36±0.37 THz (full window) to 1.91±0.05 THz (reduced window), and the FWHM standard deviations at 1300–1350 nm are as large as 0.45 THz. These low-frequency points anchor the γ0 plateau, so a wavelength-dependent window artifact would mimic exactly the reported step. The authors should report per-wavelength γ values with uncertainties and show that the step survives systematic variation of the FFT windows and probe-bandwidth modeli
  2. [Fig. 4d / Discussion] The claim that Eq. (3) with the fitted γ(Ω) 'reproduces' the measured polariton damping is not an independent check. The γ(Ω) values were obtained by fitting Eqs. (5)–(6) to the same experimental spectra whose FWHM define the polariton damping points in Fig. 3d. The two-layer fit (γ per spectrum, then tanh to those γ values, then Eq. (3)) is self-consistent by construction, not validated. A stronger test would be to hold out some probe wavelengths, fit the remaining spectra, and predict the omitted damping points, or to compare the fitted line shapes at all wavelengths with a quantitative goodness-of-fit metric.
  3. [Eqs. (5)–(6) and Appendix G] The theoretical line-shape model is the basis for extracting γ, but several of its assumptions are not tested against the low-frequency Signal-2 data that carry the step. The probe is modeled as a Gaussian pulse with assumed 100–120 fs duration and the manuscript states a ±10 nm bandwidth without giving a measured spectrum; the model also does not reproduce the satellite features near 3 THz in Signal 2 (main text, gray arrow). Because these omissions directly affect the fitted linewidth, the sensitivity of γ to the assumed probe duration, bandwidth, and the treatment of the unmodeled satellites should be quantified. At present the reader cannot tell whether the step is intrinsic or an artifact of the incomplete forward model.
minor comments (5)
  1. [Fig. 4c] The label 'γ∞ΤΡΤΥΙΓ' in Fig. 4c appears to contain corrupted characters; it should read 'γ∞'.
  2. [Author list] The corresponding-author email 'elsabreu@pyhs.ethz.ch' contains a typo ('pyhs' should be 'phys').
  3. [Methods / main text] The statement that the probe has '±10 nm' bandwidth is not supported by the Methods, which give only the OPA pulse duration (120 fs). If the bandwidth was measured, provide the value and uncertainty; otherwise temper the claim.
  4. [Fig. 3d] The legend of Fig. 3d uses 'Eq. 2' and 'Eq. 1' to label red dashed lines that correspond to phase-matching conditions, but panel d shows damping data and theory from Eq. (3). Please clarify the legend to avoid confusion.
  5. [Main text, Discussion] The sentence referring to the gray arrow in 'Fig. 1c' should cite Fig. 2c, where Signal-2 satellite features are shown.

Circularity Check

1 steps flagged

Damping 'prediction' is a two-stage fit: γ is fitted per spectrum and then re-inserted into Eq. (3) to reproduce the same measured polariton damping.

specific steps
  1. fitted input called prediction [Main text, Theory section (Fig. 4c,d; Eqs. (3)–(7)) with Appendix H]
    "The extracted values for γ are shown in Fig. 4c, as a function of the frequency of the main peak of the spectra, revealing a non-constant behavior. ... The curve γ(Ω) is then used in Eq. (4) to re-evaluate the polariton dispersion via Eq. (3). While the real part of the dispersion is unaffected by the frequency-dependent damping rate, the imaginary part displays an anomalous behavior, as shown in Fig. 4d."

    Appendix H says 'The phonon damping rate γ is left as a free parameter in the calculations and changes between different spectra': each γ point is obtained by matching the measured FFT spectrum to Eqs. (5)–(6), with γ in both K^(2) (Eq. 7) and n_THz (Eq. 4). A tanh is fitted to those extracted γ values; re-inserting γ(Ω) into Eq. (4) and solving Eq. (3) gives the Im Ω(k) curve in Fig. 4d, compared with the same experimental FWHM dots used to get γ. The agreement is a consistency check of the fit, not independent confirmation; unmodeled linewidth is absorbed into γ by construction.

full rationale

The real-part dispersion reconstruction is self-contained: FFT peak frequencies are converted to momenta by the phase-matching equations (1)–(2) and compared with the wave-equation dispersion (3) using fixed literature values (ε∞=22.47, ωTO/2π=4.44 THz, ωLO/2π=5.94 THz from Ref. [15]); no circularity there. The circularity concerns the central damping claim. Figure 4c is obtained by fitting γ as the only free parameter of Eqs. (5)–(6) to each measured spectrum, so the resulting γ(Ω) is an inverse solution, not an independent observable. Fitting a tanh to these points and then using that γ(Ω) in Eq. (4)/(3) to draw Fig. 4d, which is overlaid on the same experimental FWHM data, means the agreement is enforced by the fitting procedure rather than being a prediction. The paper does acknowledge γ is a free fitting parameter, but its language ('revealing', 'uncover') treats the fit output as a discovery. This is partial circularity: the real dispersion and the qualitative anomalous damping relative to constant-γ curves retain independent content, but the specific γ(Ω) step and its Fig. 4d confirmation reduce by construction. Unmodeled broadening (FFT-window choice, probe bandwidth) is a correctness risk, not a circularity argument; self-citation of Ref. [32] is not flagged because it is an externally published framework and no uniqueness theorem is invoked.

Axiom & Free-Parameter Ledger

5 free parameters · 5 axioms · 0 invented entities

The real-dispersion reconstruction uses fixed literature inputs (n_o, n_e, ε∞, ω_TO, ω_LO) and is not circular. The damping story, however, is built from two layers of fitting: γ is fitted per experimental spectrum, and the step-like γ(Ω) is a four-parameter tanh fit to those fitted values. No independent observable (e.g., an anharmonic calculation, sub-1 THz measurements, or 800-1200 nm probe data) is used to fix the step. The model assumptions are standard perturbative three-wave-mixing and single-Lorentz-oscillator propagation; they are reasonable but not separately validated in this paper.

free parameters (5)
  • bare E(TO1) phonon damping γ = frequency-dependent; individual per-wavelength values in Fig. 4c, ranging roughly 0.96-3.15 THz
    Free parameter in K(2)(ω',Ω) and in n_THz(Ω); fitted separately to each TP-RP spectrum in Appendix H. This is the quantity whose frequency dependence is the paper's main physical claim.
  • γ0 (low-frequency damping parameter) = 3.15 THz
    Parameter of the phenomenological tanh curve γ(Ω) fitted to the per-wavelength γ values; no independent microscopic calculation.
  • γ∞ (high-frequency damping parameter) = 0.96 THz
    Parameter of the tanh curve fitted to extracted γ values.
  • Ω0/2π (step threshold) = 2.81 THz
    Crossover frequency of the tanh curve; fit to extracted γ values, and the basis for the acoustic-phonon decay interpretation.
  • Δ/2π (step width) = 0.36 THz
    Width parameter of the tanh curve; fitted, not derived.
axioms (5)
  • domain assumption The measured pump-probe signal is proportional to the first-order nonlinear field E_NL generated by a second-order current, with Fabry-Pérot factors set to 1 and oscillating exponential terms cut.
    Invoked in Appendix G and in Eqs. (5)-(6). The signal is assumed dominated by three-wave mixing and single-pass propagation; unmodeled multiple reflections or higher-order nonlinearities could affect the line shapes from which γ is extracted.
  • domain assumption The THz dielectric response is described by a single damped Lorentz oscillator Eq. (4) with the same γ entering the nonlinear kernel.
    Explicitly stated: Eq. (4) 'is valid only for a single Lorentz oscillator' and describes only the lowest phonon-polariton branch. Other modes are neglected in the fitting of γ.
  • domain assumption A narrowband probe justifies the monochromatic phase-matching equations Eqs. (1)-(2).
    The paper states Eqs. (1)-(2) 'are reliable when a narrowband probe field is used' (±10 nm). The finite probe bandwidth is included only in the full model, and its residual effect enters linewidth extraction.
  • domain assumption The measured transmitted THz field adequately represents the backward-propagating pump A^r_p after reflection at the rear interface.
    In the Theory Model section, A^r_p is modeled using the measured transmitted THz field combined with transmission/reflection coefficients; the internal reflection amplitude is not directly measured.
  • standard math Standard diagrammatic perturbation theory and Maxwell-Fresnel propagation give the correct nonlinear response.
    The framework generalizes Ref. [32] from cubic to uniaxial crystals; accepted theory, but the specific application to TP-RP line shapes is not independently bench-marked beyond this experiment.

pith-pipeline@v1.3.0-alltime-deepseek · 4479 in / 4514 out tokens · 157849 ms · 2026-08-01T20:08:35.289312+00:00 · methodology

0 comments
read the original abstract

Mapping the dispersion of polaritons, hybrid quasiparticles arising from light-matter coupling, can provide key insights into the material dielectric response, coupling strength, and energy transfer pathways with other excitations. In this work, we present THz pump-Raman probe (TP-RP) as a versatile method for mapping the polariton dispersion in polar non-centrosymmetric materials, demonstrated here for the case of phonon-polaritons in LiNbO$_3$. By resonantly driving polaritonic modes with a broadband THz pump and probing them with a tunable NIR Raman pulse, TP-RP allows for the extraction of the momentum-dependence of both their frequency and damping rate with high accuracy. The spectral features observed in the pump-probe signal, including the polaritonic response as well as pulse artifacts, are reproduced within a many-body theoretical approach. Applying the technique to study the E(TO$_1$) phonon of LiNbO$_3$ enables the combined analysis of theory and experiments to uncover a nontrivial frequency dependence of the phonon intrinsic damping rate, revealing possible anharmonic couplings to other modes.

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

61 extracted references · 14 canonical work pages

  1. [1]

    Low, T. et al. Polaritons in layered two- dimensional materials. Nature Materials16, 182–194 (2017). URL https://doi.org/10.1038/ nmat4792

  2. [2]

    N., Fogler, M

    Basov, D. N., Fogler, M. M. & de Abajo, F. J. G. Polaritons in van der Waals materials. Science 354, aag1992 (2016). URL https://www.science. org/doi/abs/10.1126/science.aag1992

  3. [3]

    N., Asenjo-Garcia, A., Schuck, P

    Basov, D. N., Asenjo-Garcia, A., Schuck, P. J., Zhu, X. & Rubio, A. Polariton panorama. Nanophotonics10, 549–577 (2021). URL https: //doi.org/10.1515/nanoph-2020-0449

  4. [4]

    Basov, D. et al. Polaritonic quantum matter. Nanophotonics14, 3723–3760 (2025). URL https: //doi.org/10.1515/nanoph-2025-0001

  5. [5]

    Kuzmenko, A. B. Kramers–Kronig con- strained variational analysis of optical spectra. Review of Scientific Instruments76, 083108 (2005). URL https://doi.org/10.1063/1.1979470

  6. [6]

    Caldwell, J. D. et al. Low-loss, infrared and tera- hertz nanophotonics using surface phonon polari- tons. Nanophotonics4, 44–68 (2015). URL https: //doi.org/10.1515/nanoph-2014-0003

  7. [8]

    Feurer, T. et al. Terahertz Polaritonics. Annual Review of Materials Research37, 317–350 (2007). URL https://doi.org/10.1146/annurev. matsci.37.052506.084327

  8. [9]

    & Nelson, K

    Kampfrath, T., Tanaka, K. & Nelson, K. A. Reso- nant and nonresonant control over matter and light by intense terahertz transients. Nature Photonics 7, 680–690 (2013). URL https://doi.org/10.1038/ nphoton.2013.184

  9. [10]

    Zeng, Z. et al. Photo-induced chirality in a nonchiral crystal. Science387, 431–436 (2025). URL https://www.science.org/doi/abs/10. 1126/science.adr4713

  10. [11]

    Baron, A. Q. R. High-Resolution Inelastic X-Ray Scattering I: Context, Spectrometers, Samples, and Superconductors (2020). URL https://doi.org/10. 1007/978-3-030-23201-6 41

  11. [12]

    Shirane, G., Shapiro, S. M. & Tranquada, J. M. Neutron scattering with a triple-axis spectrometer: Basic techniques (2002). URL https://doi.org/10. 1017/CBO9780511534881

  12. [13]

    & Cook, J

    Ollivier, J., Plazanet, M., Schober, H. & Cook, J. First results with the upgraded IN5 disk chopper cold time-of-flight spectrom- eter. Physica B: Condensed Matter350, 173– 177 (2004). URL https://doi.org/10.1016/j.physb. 2004.04.022

  13. [14]

    Lory, P.-F. et al. Direct measurement of individ- ual phonon lifetimes in the clathrate compound Ba7. 81Ge40. 67Au5. 33. Nature communications 8, 491 (2017). URL https://doi.org/10.1038/ s41467-017-00584-7

  14. [16]

    Luo, T. et al. Time-of-flight detection of tera- hertz phonon-polariton. Nature Communications 15, 2276 (2024). URL https://doi.org/10.1038/ s41467-024-46515-1

  15. [17]

    Henry, C. H. & Hopfield, J. J. Raman Scatter- ing by Polaritons. Phys. Rev. Lett.15, 964–966 (1965). URL https://link.aps.org/doi/10.1103/ PhysRevLett.15.964

  16. [18]

    Planken, P. C. M., Noordam, L. D., Kennis, J. T. M. & Lagendijk, A. Femtosecond time-resolved study of the generation and propagation of phonon 10 polaritons in LiNbO 3. Phys. Rev. B45, 7106– 7114 (1992). URL https://link.aps.org/doi/10. 1103/PhysRevB.45.7106

  17. [19]

    F., Stoyanov, N

    Crimmins, T. F., Stoyanov, N. S. & Nelson, K. A. Heterodyned impulsive stimulated Raman scatter- ing of phonon–polaritons in LiTaO 3 and LiNbO 3. The Journal of Chemical Physics117, 2882–2896 (2002). URL https://doi.org/10.1063/1.1491948

  18. [20]

    Wahlstrand, J. K. & Merlin, R. Cherenkov radi- ation emitted by ultrafast laser pulses and the generation of coherent polaritons. Phys. Rev. B 68, 054301 (2003). URL https://link.aps.org/doi/ 10.1103/PhysRevB.68.054301

  19. [21]

    Leitenstorfer, A. et al. The 2023 ter- ahertz science and technology roadmap. Journal of Physics D: Applied Physics56, 223001 (2023). URL https://doi.org/10.1088/1361-6463/ acbe4c

  20. [22]

    Caruso, F. et al. The 2025 roadmap to ultra- fast dynamics: frontiers of theoretical and compu- tational modeling. Journal of Physics: Materials 9, 012501 (2025). URL https://doi.org/10.1088/ 2515-7639/ae1165

  21. [23]

    S., Hall, J

    Dastrup, B. S., Hall, J. R. & Johnson, J. A. Experimental determination of the interatomic potential in LiNbO 3 via ultrafast lattice control. Applied Physics Letters110, 162901 (2017). URL https://doi.org/10.1063/1.4980112

  22. [24]

    Knighton, B. E. et al. Terahertz waveform consid- erations for nonlinearly driving lattice vibrations. Journal of Applied Physics125, 144101 (2019). URL https://doi.org/10.1063/1.5052638

  23. [25]

    & Blake, G

    Lin, H.-W., Mead, G. & Blake, G. A. Mapping LiNbO3 Phonon-Polariton Nonlinearities with 2D THz-THz-Raman Spectroscopy. Phys. Rev. Lett. 129, 207401 (2022). URL https://link.aps.org/ doi/10.1103/PhysRevLett.129.207401

  24. [26]

    Schwarz, U. T. & Maier, M. Frequency dependence of phonon-polariton damping in lithium niobate. Phys. Rev. B53, 5074–5077 (1996). URL https: //link.aps.org/doi/10.1103/PhysRevB.53.5074

  25. [27]

    Weis, R. S. & Gaylord, T. K. Lithium Niobate: Summary of Physical Properties and Crystal Structure. Appl. Phys. A37, 191–203 (1985). URL https://link.springer.com/article/10. 1007/BF00614817

  26. [28]

    & Mansingh, A

    Dhar, A. & Mansingh, A. Optical proper- ties of reduced lithium niobate single crystals. Journal of Applied Physics68, 5804–5809 (1990). URL https://doi.org/10.1063/1.346951

  27. [30]

    Zhang, C. et al. Bandwidth tunable THz wave generation in large-area periodically poled lithium niobate. Opt. Express20, 8784–8790 (2012). URL https://opg.optica.org/oe/abstract. cfm?URI=oe-20-8-8784

  28. [31]

    Lu, Y. et al. Giant enhancement of THz- frequency optical nonlinearity by phonon polari- ton in ionic crystals. Nature Communications 12, 3183 (2021). URL https://doi.org/10.1038/ s41467-021-23526-w

  29. [32]

    P., Benfatto, L

    Sellati, N., Fiore, J., Villani, S. P., Benfatto, L. & Udina, M. Theory of terahertz pump optical probe spectroscopy of phonon polaritons in non- centrosymmetric systems. npj Quantum Materials 10, 46 (2025). URL https://doi.org/10.1038/ s41535-025-00761-8

  30. [33]

    Wiederrecht, G. P. et al. Explanation of anomalous polariton dynamics in LiTaO 3. Phys. Rev. B51, 916–931 (1995). URL https://link.aps.org/doi/10. 1103/PhysRevB.51.916

  31. [34]

    Parlinski, K., Li, Z. Q. & Kawazoe, Y. Ab initio calculations of phonons in LiNbO 3. Phys. Rev. B 61, 272–278 (2000). URL https://link.aps.org/doi/ 10.1103/PhysRevB.61.272

  32. [35]

    Rader, C. et al. A New Standard in High-Field Terahertz Generation: the Organic Nonlinear Opti- cal Crystal PNPA. ACS Photonics9, 3720–3726 (2022). URL https://pubs.acs.org/doi/10.1021/ acsphotonics.2c01336

  33. [36]

    Sanna, S. et al. Raman scattering efficiency in LiTaO3 and LiNbO 3 crystals. Phys. Rev. B91, 224302 (2015). URL https://link.aps.org/doi/10. 1103/PhysRevB.91.224302

  34. [38]

    & Benfatto, L

    Udina, M., Cea, T. & Benfatto, L. Theory of coherent-oscillations generation in terahertz pump-probe spectroscopy: From phonons to elec- tronic collective modes. Phys. Rev. B100, 165131 (2019). URL https://link.aps.org/doi/10.1103/ PhysRevB.100.165131

  35. [39]

    P., Wiederrecht, G

    Dougherty, T. P., Wiederrecht, G. P. & Nelson, K. A. Impulsive stimulated Raman scat- tering experiments in the polariton regime. J. Opt. Soc. Am. B9, 2179–2189 (1992). URL 11 https://opg.optica.org/josab/abstract.cfm?URI= josab-9-12-2179

  36. [40]

    & Juraschek, D

    Yaniv, O. & Juraschek, D. M. Phonon polari- ton Hall effect (2025). URL https://arxiv.org/abs/ 2509.16100. arXiv:2509.16100

  37. [41]

    L., Knighton, B

    Johnson, C. L., Knighton, B. E. & John- son, J. A. Distinguishing Nonlinear Tera- hertz Excitation Pathways with Two-Dimensional Spectroscopy. Phys. Rev. Lett.122, 073901 (2019). URL https://link.aps.org/doi/10.1103/ PhysRevLett.122.073901

  38. [42]

    Terahertz polariton dis- persion in uniaxial optical crystals

    Kojima, S. Terahertz polariton dis- persion in uniaxial optical crystals. Progress In Electromagnetics Research Letters 77, 109–115 (2018). URL https://doi.org/10. 2528/PIERL18050505

  39. [43]

    Basini, M. et al. Terahertz ionic kerr effect: Two-phonon contribution to the nonlinear opti- cal response in insulators. Phys. Rev. B109, 024309 (2024). URL https://link.aps.org/doi/10. 1103/PhysRevB.109.024309

  40. [44]

    & Benfatto, L

    Sellati, N., Fiore, J. & Benfatto, L. Light-induced faraday effect from dynamical breakdown of klein- man symmetry (2026). URL https://arxiv.org/ abs/2605.27127. arXiv:2605.27127

  41. [45]

    F., Wang, F., Liu, Y

    Huber, L., Maehrlein, S. F., Wang, F., Liu, Y. & Zhu, X.-Y. The ultrafast Kerr effect in anisotropic and dispersive media. The Journal of Chemical Physics154, 094202 (2021). URL https://doi.org/10.1063/5.0037142

  42. [46]

    Fiore, J. et al. Investigating Josephson plas- mons in layered cuprates via nonlinear terahertz spectroscopy. Phys. Rev. B110, L060504 (2024). URL https://link.aps.org/doi/10.1103/PhysRevB. 110.L060504

  43. [47]

    & Benfatto, L

    Fiore, J., Sellati, N., Udina, M. & Benfatto, L. Two-dimensional terahertz spectroscopy in electronic systems: A many-body diagrammatic approach. Phys. Rev. B113, 174524 (2026). URL https://link.aps.org/doi/10.1103/8ffz-rtyf

  44. [48]

    Polyanskiy, M. N. Refractiveindex.info database of optical constants. Scientific Data11, 94 (2024). URL https://doi.org/10.1038/s41597-023-02898-2

  45. [49]

    & G¨ unter, P

    Stillhart, M., Schneider, A. & G¨ unter, P. Optical properties of 4-N,N-dimethylamino-4′-N′-methyl- stilbazolium 2,4,6-trimethylbenzenesulfonate crys- tals at terahertz frequencies. J. Opt. Soc. Am. B 25, 1914–1919 (2008). URL https://opg.optica. org/josab/abstract.cfm?URI=josab-25-11-1914

  46. [50]

    P., Ruchert, C., Vicario, C

    Hauri, C. P., Ruchert, C., Vicario, C. & Ardana, F. Strong-field single-cycle THz pulses generated in an organic crystal. Applied Physics Letters99, 161116 (2011). URL https://doi.org/10.1063/1. 3655331

  47. [51]

    Hine, G. A. & Doleans, M. Intrinsic spatial chirp of subcycle terahertz pulsed beams. Phys. Rev. A 104, 032229 (2021). URL https://link.aps.org/ doi/10.1103/PhysRevA.104.032229

  48. [52]

    & Blake, G

    Lin, H.-W., Hsieh, P.-H., Mead, G. & Blake, G. A. Characterization of the nonlinear thz focus for 2d thz spectroscopy. J. Opt. Soc. Am. B41, 466– 470 (2024). URL https://opg.optica.org/josab/ abstract.cfm?URI=josab-41-2-466

  49. [53]

    Jewell, M. et al. Spatiotemporal characterization of 2d thz systems using a compact electro-optic camera (2025). URL https://ieeexplore.ieee.org/ abstract/document/11319688

  50. [54]

    Planken, P. C. M., Nienhuys, H.-K., Bakker, H. J. & Wenckebach, T. Measurement and calculation of the orientation dependence of terahertz pulse detection in ZnTe. J. Opt. Soc. Am. B18, 313– 317 (2001). URL https://opg.optica.org/josab/ abstract.cfm?URI=josab-18-3-313

  51. [55]

    Aspnes, D. E. & Studna, A. A. Dielectric func- tions and optical parameters of Si, Ge, GaP, GaAs, GaSb, InP, InAs, and InSb from 1.5 to 6.0 eV. Phys. Rev. B27, 985–1009 (1983). URL https: //link.aps.org/doi/10.1103/PhysRevB.27.985

  52. [56]

    & Chirakadze, A

    Berozashvili, Y., Machavariani, S., Natsvlishvili, A. & Chirakadze, A. Dispersion of the linear electro- optic coefficients and the non-linear susceptibil- ity in GaP. Journal of Physics D: Applied Physics 22, 682 (1989). URL https://doi.org/10.1088/ 0022-3727/22/5/017

  53. [57]

    & Zhang, X.-C

    Wu, Q. & Zhang, X.-C. 7 terahertz broadband gap electro-optic sensor. Applied Physics Letters70, 1784–1786 (1997). URL https://doi.org/10.1063/ 1.118691

  54. [58]

    Nelson, D. F. & Mikulyak, R. M. Refractive indices of congruently melting lithium niobate. Journal of Applied Physics45, 3688–3689 (1974). URL https://doi.org/10.1063/1.1663839

  55. [59]

    Barker, A. S. & Loudon, R. Dielectric Properties and Optical Phonons in LiNbO 3. Phys. Rev.158, 433–445 (1967). URL https://link.aps.org/doi/10. 1103/PhysRev.158.433. Acknowledgements We thank Janine Zemp-D¨ ossegger for her assistance with preliminary experiments, and Dominik Juraschek and Michael M. Fechner for their contributions to early- stage discus...

  56. [60]

    (ω,x) (0<x< d), (G15) where J(2)

  57. [61]

    (ω,x) = Z dω1dω2A[0](ω1,x)K (2)(ω1, ω2)A[0](ω2,x)δ(ω−ω 1 −ω 2).(G16) The solution of Eq. (G15) in 0<x< dcan be written as A[1](ω,x) = A u(ω,x) + B(ω)e in(ω)ωx/c + C(ω)e−in(ω)ωx/c ,(G17) where B(ω) and C(ω) are the coefficients of the forward- and backward-propagating fields respectively, determined by boundary conditions, while A u(ω,x) is a unique soluti...

  58. [62]

    (ω, k),(G18) with J(2)

  59. [63]

    (G16) by insertion of the solution Eq

    (ω, k) obtained from Eq. (G16) by insertion of the solution Eq. (G13) as J(2)

  60. [64]

    (ω, k) = Z d 0 dx J(2)

  61. [65]

    The integral in Eq

    (ω,x)e −ikx = Z dω1dω2δ(ω−ω 1 −ω 2)K(2)(ω1, ω2) ±1X α1 ±1X α2 Aα1 (ω1)Aα2 (ω2) 1−e −i(k−α1k1−α2k2)d i(k−α 1k1 −α 2k2) ,(G19) where, to keep a compact notation, we definedk i =n(ω i)ωi/c, A+1(ω) = At(ω) and A−1(ω) = Ar(ω). The integral in Eq. (G18) can be solved with the residue theorem, and the unique solution reads explicitly Au(ω,x) = 4π c Z dω1dω2δ(ω−ω...