All-optical control of second-harmonic generation in β-BaB_2O_4 via coherent, terahertz-driven acentric lattice displacement

Domain: neuroscience
Published: 2026-04-07
Editorial outcome: minor_revision

Abstract
Dynamical control of the nonlinear optical properties of solids – with light itself – will be essential for future ultrafast photonic technologies. Previously, methods to modulate nonlinear processes including second-harmonic generation (SHG) have relied primarily on non-resonant light-matter interaction or photo-generation of hot electrons in nanoscale materials. However, these approaches are typically constrained by limited interaction lengths and the initial frequency conversion is relatively weak under equilibrium conditions. Here, a ~30% modulation of efficient phase-matched SHG in bulk beta-barium borate (β-BaB O ) is achieved through transient lattice deformation by intense terahertz (TCHz) pulses that are tuned to resonance with an infrared-active phonon mode. The effect originates from modification of the index of refraction ellipsiod and the corresponding nonlinear phase-matching conditions, rather than from direct modulation of the nonlinear susceptibility through THz-mediated 𝜒(3) processes. This mechanism, of resonant selective lattice excitation, points toward novel THz-control schemes to tune the nonlinear optical response in materials.

Introduction
Nonlinear optical frequency conversion – including second harmonic generation (SHG) – greatly extends the spectral range of conventional lasers, underpinning key scientific applications across a wide range of fields from high-resolution microscopy to precision optical metrology to quantum information science1–6. For communications and technological applications, integrated photonic circuits have emerged7–10, in which the nonlinear optical response can be precisely modulated. Control over optical nonlinearities, including SHG, using applied electric fields has now been demonstrated in optical crystals11, polymers12, semiconductors13–15, atomic-layer systems and heterostructures,16,17 and artificial photonic structures18,19. However, modulation of the nonlinear optical response with this approach is fundamentally limited by the speed of conventional electronics to GHz modulation rates, creating a functional barrier to the next generation of optoelectronic devices. In principle, orders-of-magnitude higher speeds can be achieved through all-optical approaches, which would enable ultrafast data transfer and real-time optical signal processing.

Indeed, ultrafast, all-optical SHG modulation has been demonstrated in nanostructures and two-dimensional materials20–23. However, these achievements appear impractical for applications in solid-state devices, as the modulation generally relies on either photoexcitation of hot electrons22,24 with significant collateral heating, or non-resonant light-matter interaction,25–29 which is limited by material damage thresholds. Further, while high contrast has been achieved, the absolute SHG signals themselves remain relatively weak due to limited interaction lengths. On the other hand, while bulk three-dimensional materials can deliver strong SHG signals suitable for applications, high-contrast, ultrafast, all-optical modulation has not yet been reported.

In this work, high-contrast modulation of strong SHG is achieved by resonant excitation of the β-barium borate (β-BaB₂O₄, BBO) crystal lattice using intense terahertz (THz) pulses tuned to an infrared (IR)-active phonon mode. Theoretical considerations — based on group theory and ab-initio calculations30–32 — indicate that both the linear properties and the second-order optical nonlinearity in BBO originates from localized electronic density within the B₃O₆ anionic rings in the crystal unit cell. Therefore, displacement of the anionic ring can be predicted to significantly influence SHG and targeted displacements might lead to strong, controlled modulation of the SHG signal. While this parameter cannot be selectively accessed with conventional stimuli (e.g., strain, static electric fields), selective displacement can be achieved with intense THz pulses, resonantly coupled to IR-active phonons, as illustrated in Fig. 1a. With this approach, THz pulses can be expected to enable ultrafast control over the linear and nonlinear optical properties of BBO, enabling efficient modulation of the SHG process.

The potential for control by THz-driven resonant excitation is investigated here through the SHG of a femtosecond NIR probe laser. By measuring the time-resolved SHG efficiency as a function of the fundamental input wave polarization and delay with respect to the THz excitation, the modulated SHG signal yields the overall out-of-equilibrium nonlinear optical response. And importantly, under appropriate input conditions, high-contrast SHG modulation exceeding 30% can be reached.

Macroscopically, this effect can be explained by an electro-optic, THz-induced rotation of the principal optical axes (perpendicular to the direction of propagation), which alters the SHG phase matching conditions. Density functional perturbation theory (DFPT) calculations validated by the experiment can then be used to obtain an accurate estimate of the frequency dispersion of the electro-optic coefficient in the THz range, as well as the phonon normal-mode amplitude of 0.67 Å(amu)1/2, which displaces the B O anionic unit from its equilibrium position. This quantitative link between the atomic-scale structural dynamics and macroscopic nonlinear optical response may facilitate engineering of materials with specific transient nonlinear optical responses in the future.

Results
To investigate the SHG modulation by selective lattice excitation on ultrafast timescales, intense THz pulses are used to drive the material while its linear and nonlinear response are probed with an overlapping femtosecond near-infrared (NIR) pulse (Fig. 1a-b and Secs.1 and 2 of the Supplementary Information for details). Prior to probing the dynamics, the BBO crystal is characterized in the chosen experimental geometry in which the THz pulse is polarized along the ordinary axis (o-axis) of the BBO crystal. Steady-state optical spectroscopy (see Supplementary Fig. 3 and Sec. 3 of the Supplementary Information for details) yields the imaginary part of the dielectric function 𝜀 (𝜔) and the corresponding IR-active phonon modes (Fig 1d). The dominant IR-active 𝐸-symmetry phonon mode is found at a frequency of 𝜔 = 4.32 THz, consistent with other recent independent characterizations by THz time-domain spectroscopy33. This mode is of particular interest, as density functional perturbation theory (DFPT) calculations indicate that it leads to significant co-planar motion of anionic [B₃O₆]³⁻ rings while the Ba²⁺ cations will exhibit in nearly motion in the opposite direction with comparatively smaller displacements, due to their large atomic mass31 (Fig. 1a).

Electro-optic sampling (EOS) of the driving THz pulse is used to verify spectral overlap with the 𝐸 phonon mode required for targeted excitation. In-situ EOS measurements of the THz temporal waveform 𝐸 (𝑡) yield a peak electric field strength of 𝐸 = 8.5 MV/cm, and Fourier analysis reveals a spectrum heavily weighted around 4.3 THz, with additional weaker frequency components ranging from 0−10 THz (Fig. 1c-d). The spectral content ensures efficient coupling of the THz pulse to the desired 𝐸 phonon mode. While there are additional IR-active phonon modes spanned by the broadband THz pulse, they have lower absorption cross-sections33 and are not expected to contribute significantly to the driven lattice dynamics (Fig. 1e).

Following initial characterization, time-resolved experiments are performed by focusing the THz pulse into a 300 µm thick BBO crystal cut for type-I phase-matched SHG of a ~50 femtosecond NIR laser pulse at the fundamental wavelength of 800 nm. In the experiment, 20 nJ NIR pulses are focused to a spot size of 35 µm at the BBO crystal with corresponding intensity of ~ 8×1010 W/cm2, and at equilibrium a maximum SHG energy conversion efficiency of ~ 4 % is observed.

By varying the delay between the THz pump pulse and polarization of the femtosecond NIR probe pulse, the influence of the phonon dynamics on the SHG response is investigated in the time domain as function of polarization. A key observable is the time-dependent, normalized relative strength of the SHG intensity, Δ𝐼 = 𝛥𝐼 /𝐼 . Here, 𝛥𝐼 = 𝐼 − 𝐼 , 𝐼 , and the equilibrium SHG intensity, 𝐼 . By convention, time-zero occurs when the peak of the femtosecond NIR pulse coincides with the maximum of the THz electric field. Under optimized input polarization (detuned from peak conversion efficiency) and delay the THz pulse induces a modulation in the relative SHG strength exceeding 30% (Fig. 1f), whereas no significant THz-induced modulation is observed in the NIR probe intensity (Supplementary Fig. 2).

The time-dependent SHG exhibits a damped oscillatory behavior persisting for over 4 ps, with Fourier analysis revealing a distinct peak slighlty below the phonon frequency, 𝜔 , coincident with the positive maximum of real part of the dielectric fuction, ε , see Fig. 1g. Tuning the spectrum of the THz driving field with spectral filters to above and below the phonon frequency 𝜔 . When the spectral components of the driving THz field are above 𝜔 , the SHG modulation is suppressed. Similarly, when the THz cut-off frequency is below 𝜔 , modulation of the SHG is also strongly suppressed and can be attributed to weak excitation of lower frequency phonon modes (Fig. 2a).

The phase dependence of the SHG modulation is further investigated by returning to the nominal THz pulse (tuned to resonance) and inverting its polarity. A corresponding inversion in the relative SHG intensity is observed, indicating that the dynamics are resonantly driven and phase-locked to the THz field (Fig. 2b). Additionally, scaling the THz field strength results in a linearly dependent effect, suggesting that the SHG modulation is directly proportional to the phonon amplitude, which is itself linearly dependent on the THz driving field strength34–36 (Fig. 2c).

At equilibrium, the expression for the type-I SHG in the undepleted pump approximation in BBO of length 𝐿 by an electromagnetic wave of intensity 𝐼, frequency 𝜔, and polarisation angle 𝛼 is given by (see Supplementary Information Sec. 4):

𝐼 (𝛼) = 𝑁𝐿 ( 2𝜔2𝑑2 𝐿2 sin( Δ𝑘 𝐿) / 𝑛 𝑛2𝑐3𝜖 1 ) cos4(𝛼)𝐼2(𝜔), (1)

where 𝑛 (𝜔) is the ordinary refractive index at the fundamental (𝜔) and 𝑛 (2𝜔) is the extraordinary refractive index at second harmonic (2𝜔), 𝑑 is the nonlinear optical coefficient, and Δ𝑘 represents the momentum mismatch between the fundamental and second harmonic waves (Δ𝑘 = 0 in perfect phase matching conditions). The polarisation angle 𝛼 defines the orientation of the fundamental wave with respect to the ordinary axis, as illustrated in Fig. 3a. According to Eq. (1), and consistent with equilibrium experimental observations, the SHG intensity 𝐼 (𝛼) exhibits a cosine-to-the-fourth power dependence on input polarisation, resulting in a two-fold symmetric angular distribution with maximum intensity at 𝛼 = 0 and 𝛼 = 180° (Fig. 3b).

Out-of-equilibrium, the SHG can in principal be modulated by two distinct effects. The dynamics can result from a time-dependent variation in the second-order nonlinear coefficient, which is proportional to the driving THz field and can be considered effectively as a third-order 𝜒(3) effect. Notably, under certain circumstances, this 𝜒(3 )-driven effect has been used to describe THz-Field Induced Second Harmonic Generation (TFISH) 37–40. Alternatively, the dynamics can result from phonon-mediated electro-optic effects in the linear index of refraction (a 𝜒(2) effect involving THz and phonon frequencies) – and subsequent modification of the phase-matching conditions (Δ𝑘), affecting the SHG efficiency, without inducing a time-dependent 𝜒(2) or 𝜒(3) effect, through the second-order nonlinear coefficient involving the fundamental field frequency (Secs. 5 and 6 of the Supplementary Information). This mechanism has been described as a cascade of 𝜒(2) processes26,41. Under strong-field THz excitation, both 𝜒(2) and 𝜒(3) effects are always present, however, it can be expected that only one of the effects dominates the nonlinear response depending on the specific equilibrium phase-matching conditions.

To identify the origin of the observed dynamics, the NIR probe pulse polarization is scanned at a time delay of 𝑡 = 2.4 ps - corresponding to a local maximum of the SHG dynamics. Notability, the angular dependence of the measured SHG modulation, ∆𝐼 (𝛼) = 𝐼 (𝛼)−I𝐼 (𝛼), now exhibits a four-fold symmetry, with the two lobes peaking at 𝛼 ≈ 30° and 𝛼 ≈ 330° having a positive modulation (∆𝐼 (𝛼) > 0), corresponding to an increase in SHG intensity (Fig 3) and vice-versa for the lobes observed at 𝛼 ≈ 150° and 𝛼 ≈ 210°. Further, as the temporal delay between the THz and NIR pulses is scanned over one oscillation cycle, the sign of ∆𝐼 (𝛼) is inverted in all lobes, passing through a common zero-crossing where no modulation is observed (∆𝐼 (𝛼) = 0), as shown in Fig. 3d.

Importantly, this angular and temporal dependence of the SHG signal in the type-I phase matched geometry is incompatible with a TFISH interpretation involving an effective 𝜒(3)-effect resulting from a THz-induced modulation of the second-order 𝜒(2) susceptibility. Evidence of this inconsistency is obtained from comparative measurements of SHG modulation in BBO where the phase-matching conditions are not met (see Supplementary Fig. 5). Here, due to the limited coherence length, the the 𝜒(3) effects to be isolated, as electro-optically induced changes in the index of refraction ellipsoid are negligible in the overall SHG response. In the phase-mismatched conditions, the 𝜒(3)-effects dominate, resulting in a polarization dependence in the SHG modulation with features fundamentally different from those observed in the phase matched case. Additionally, the equilibrium signal strength and magnitude of the modulated intensity is orders of magnitude weaker (see Secs. 7 and 8 of the Supplementary Information for full explanation). Therefore, in the current experiment, the strength of the effect and the 4-fold symmetry leads to the conclusion that the dynamic effect is dominated by the phonon-mediated alteration of the type I phase-matching conditions – a cascaded 2nd-order effect.

Resonant excitation of the 𝐸 phonon mode by the THz pulse displaces the (B O )3− anionic rings far from equilibrium positions (Fig 1a), thereby inducing an electro-optic modulation of the refractive index along the crystal’s principal axes (𝑥′,𝑦′,𝑧′; see Fig. 1a). Under equilibrium conditions, BBO is a negative uniaxial crystal (𝑛 = 𝑛 ≠ 𝑛 ), but the phonon excitation breaks its axial symmetry, leading to biaxial behaviour where all three refractive indices become distinct (𝑛 ≠ 𝑛 ≠ 𝑛 ) – see Sec. 5 of the Supplementary Information for details. The dynamical refractive indices along the 𝑥 and 𝑦 axes in the crystal-frame become: 𝑛𝑇𝐻𝑧 = 𝑛 +𝛿𝑛 and 𝑛𝑇𝐻𝑧 = 𝑛 −𝛿𝑛, where the magnitude of the phonon-induced refractive index modulation is:

𝛿𝑛 = 𝑛3 𝑟 𝑄/2 (2)

Here, 𝑄 is the phonon amplitude and 𝑟 is the phonon-optic coupling coefficient, relating the modulation of the refractive index to the phonon amplitude. To interpret the index of refraction ellipsoid in the context of phase matching, which governs the efficiency of the SHG 𝜒(2)(2𝜔,𝜔,𝜔)-process, a projection normal to the direction of propagation in the laboratory frame is taken from which the indices along the ordinary and extraordinary axes, known as the principal axes, are found.

In the experiment, the BBO crystal was cut for optimal type-I phase matching at equilibrium (SHG at 800 nm), with its crystal axes oriented at θ = 29.3° and φ = 90° (following the conventions in Ref. 42). The cut angle defines the direction of propagation in the crystal system and ensures that the refractive index of the ordinary fundamental wave (𝑛 ) matches that of the extraordinary second harmonic wave (𝑛 ). In the dynamical case, the phonon-induced modulation of the index of refraction ellipsoid (𝛿𝑛 in the crystal coordinate system) leads to a rotation of the ordinary and extraordinary axes in the laboratory frame by an angle 𝜾 = 𝜅(𝜃) 𝑟 𝑄, as illustrated in Fig. 4a. Here, 𝜅(𝜃) is a geometrical parameter that depends on the crystal orientation 𝜃.

While the resonant phonon excitation rotates the ordinary and extraordinary axes, impacting the phase-matching conditions in the SHG process, the magnitudes of the ordinary and extraordinary refractive indices – apart from a third-order correction – remain nearly unchanged. (Derivations of the modulation in the refractive index 𝛿𝑛, the effective rotation angle 𝜾 of the principal axes, and all relevant details are given in Secs. 5 and 6 of Supplementary Information).

To quantitatively model the dynamics, it is also essential to consider the spatio-temporal propagation of the phonon inside the BBO crystal. The phonon amplitude decays exponentially along the direction of propagation. This spatial dependence of the phonon amplitude can be described by 𝑄(𝑧) = 𝑄 𝑒 −𝑧/𝛿𝑝ℎ, where 𝑄 is the initial phonon amplitude at the crystal surface, and 𝛿 is its penetration depth. By accounting for the phonon-induced rotation of the ordinary and extraordinary optical axes, as well as the penetration depth of the phonon, the variation in SHG intensity as a function of the input polarization angle 𝛼 (at fixed THz-NIR delay) is the given by:

Δ𝐼 (𝛼) = 4 𝐼 (0) 𝜾 Λ(𝛿 ,𝐿) / 𝐿 cos3𝛼sin𝛼, (3)

where 𝐼 (0) is the equilibrium SH intensity at 𝛼 = 0, and 𝜾 = 𝜅(𝜃) 𝑟 𝑄 is the maximum rotation angle of the optical axes induced by the phonon amplitude 𝑄 at the fixed delay. The 4-fold dependence on the input polarization, ~cos3𝛼sin𝛼, is consistent with the alternating sign and angular symmetry observed experimentally, and the strength of the effect depends on the BBO thickness, 𝐿 and Λ(𝛿 ,𝐿), which includes the phonon penetration depth.

A key parameter is the phonon penetration depth 𝛿 , which can be estimated independently using finite-difference time-domain (FDTD) methods to calculate the spatio-temporal dynamics of the phonon amplitude. With FTIR measurements of the linear optical properties of the BBO crystal, EOS measurements of the THz driving electric field, and the phonon equation of motion, the normalized amplitude 𝑄(𝑧,𝑡) of the phonon at 𝜔 is simulated using FDTD and shown as a function of the crystal depth and time in Fig. 4b. Integrating the phonon amplitude in time and fitting the depth-dependent profile with an exponential decay (Fig. 4c) yields a phonon penetration depth of 𝛿 =30.2 µm.

Using 𝛿 , angularly dependent SHG measurements (at fixed time delay) for THz pulses of varying intensity are then fit with Eq. 3, keeping 𝐼 (0) as a constant parameter. These fits yield 𝜾 as a function of the THz field amplitude (Fig. 4d-e). The rotation of the ordinary and extraordinary axes, 𝜾 , is found to vary linearly with applied THz field amplitude, consistent with the linear dependence of the IR-active phonon amplitude 𝑄 on the THz field strength. The value of 𝜾 achieved for the maximum THz field strength in the experiment is 𝜾 = 0.216 𝑟𝑎𝑑.

The polarization-resolved measurements in conjunction with the analytical analysis yield the rotation of the optical axes due to the resonant THz-phonon excitation. This, in turn, permits quantification of the associated modulation of the refractive indices 𝛿𝑛 in the crystal coordinate system, but does not immediately allow for independent retrieval of the electro-optic coefficient 𝑟𝐸 or a determination of the phonon displacement 𝑄. At the microscopic level, the effect of the THz-induced coherent lattice vibrations on the local electronic environment within the anionic units remains hidden, as does the electro-optic coefficient 𝑟𝐸.

To access this information density functional theory (DFPT) can be used in the frozen-phonon approximation, to calculate the electronic band structure32 - and, by extension, the dielectric function corresponding to a lattice displacement 𝑄 along the phonon coordinates. The resulting dielectric function can then be used to extract absolute values for the refractive indices, 𝑛𝑇𝐻𝑧 and 𝑛𝑇𝐻𝑧 in the (𝑥′,𝑦′,𝑧′) coordinate system of the crystal as a function of the phonon displacement 𝑄 (Fig. 5a). Starting from equilibrium (𝑄 = 0), both indices are found to be 𝑛 = 𝑛 = 1.705, in good agreement with experimental measurements and previous ab-initio calculations32. With a finite phonon amplitude 𝑄₀, the refractive indices respond linearly as expected: 𝑛𝑇𝐻𝑧 increases, while 𝑛𝑇𝐻𝑧 decreases with increasing 𝑄₀ (see Fig. 5a).

Using the DFPT-calculated dependence of 𝛿𝑛 on 𝑄₀, the phonon-optic coupling coefficient 𝑟₂₂ can be estimated from the relation, 𝑟 = 2𝛿𝑛(𝑄 )/𝑛3 𝑄 , yielding a value of 𝑟 = 0.0049 Å−1(amu)−1/2. This coefficient is especially relevant, as it can be used to calculate the maximum phonon amplitude at maximum THz field strength using the experimentally determined rotation of the optical axes 𝜾 at the crystal surface, 𝑄 = 𝜾 /(𝜅 𝑟 )= 0.67 Å(amu)1/2(Fig.5a). This value is consistent with expectations based on experiments targeting different effects in various materials that are nevertheless similar in that they also require THz excitation of large amplitude phonons. In these independent experiments, performed at X-ray free-electron lasers, it was possible to make a direct structural determination of the phonon amplitude35,43.

Finally, the electro-optic coefficient 𝑟𝐸(𝜔) can be deduced from phonon-optic coupling coefficient 𝑟₂₂ using the conventional linear relation between the THz field and 𝑄 in the frequency domain, which leads to the following expression:

𝑟𝐸(𝜔) = 𝑟𝐸(0) (𝜔2 −𝜔2) / ( (𝜔2 −𝜔2) −(Γ𝜔)2 ) , (4)

where 𝜔 and Γ are the phonon linewidth and the 𝑟𝐸(0) is the DC electro-optic coefficient (see Sec. 9 of the Supplementary Information for details). Approaching the phonon resonance, the 𝑟𝐸(𝜔) peak value is 55 𝑝𝑚/𝑉 at 𝜔 = √𝜔2 −Γω = 4.22 THz. Notably, the spectrum of the time-dependent SHG modulation dynamics closely matches the EO coefficint, with shared maxima just below 𝜔 (dashed line Fig. 5b). This frequency dependence provides further evidence that the SHG modulation arises from the electro-optic effect. Indeed, in phonon-mediated EO processes, the response is governed not by the imaginary part of the dielectric function, 𝜀 (𝜔), which peaks at the phonon frequency 𝜔 , but by the real part, 𝜀 (𝜔), whose maximum occurs at 𝜔 .

The peak value of the vibrational electro-optic coefficient can also be estimated experimentally by combining the analytical expression for the SHG modulation with the phonon equation of motion (see Sec. 9 of the Supplementary Information). Using this approach and the broadband THz waveform applied in the experiment, a peak value of 𝑟𝐸 of ~140 𝑝𝑚/𝑉 is obtained in relatively good agreement with that obtained by DFPT (see Supplementary Fig. 7). To exclude the involvement of other phonon modes, additional experiments were performed using a narrowband THz pulse, which yields a similar value of ~170 pm/V (see Supplementary Fig. 8). The discrepancy of less than a factor of 3 between theory and experiment, which exists in both the broadband and narrowband case could therefore stem from the calibration of the absolute THz field amplitude, phonon parameters, and on the optical constants involved in the analysis—including their large frequency dispersion near the phonon resonance and dependence on the geometrical parameters in the experiment.

On the other hand, the theoretical treatment could also be slightly flawed, as extrapolating the theoretical electro-optic coefficient to the DC limit gives 𝑟𝐸(0) = 6.6 𝑝𝑚/𝑉, which is larger than the commonly reported experimental DC value of approximately 2.5 pm/V for BBO 45.

Regardless of the small discrepancy that remains, relative agreement between the experimental and theoretically determined electro-optic coefficient lends credence to the DFPT calculations, by which additional parameters can be obtained in the vicinity of the phonon resonance. In particular, the interaction between the phonon and the local electronic environment can be explored computationally. Here, DFPT-calculated differential charge density maps Δρ(𝑥′,𝑦′) confirm that the phonon mode predominantly drives cooperative displacements of the anionic [B₃O₆]³⁻ rings within the lattice plane, while the Ba ²⁺ cations remain largely stationary (Fig 5c). This selective atomic motion produces a marked redistribution of charge density around the planar [B₃O₆]³⁻ units. The electronic density around the oxygen atoms is largely influenced by the phonon dynamics, suggesting that the local electronic environment of these atoms plays a central role in modulating the refractive indices 𝑛𝑇𝐻𝑧 and 𝑛𝑇𝐻𝑧.

Methods
Intense THz pulses were generated through optical rectification in a DSTMS organic crystal using the signal beam of an optical parametric amplifier, driven by 50 fs, 800 nm pulses at a 100 Hz repetition rate. The SHG dynamics in β-BaB₂O₄ were investigated using a time-delayed replica of the 800 nm pulses. The experimental setup is shown in Supplementary Fig. 1; further details are provided in Sec. 1 of the Supplementary Information. All experiments were performed at room temperature. The sample used was a commercially available β-BaB₂O₄ single crystal (5 mm × 5 mm × 0.3 mm), oriented at 𝜃 = 29.3° and 𝜙 = 90°. Crystal orientation and crystallographic quality were verified by Laue X-ray diffraction (Supplementary Fig. 4).

The THz-driven phonon propagation along the crystal depth has been simulated in time domain by computing Maxwell’s equations based on FTDT method48. The material optical response along the o-axis at THz frequencies was modeled including the dielectric function of IR-active phonons, whose parameters have been experimentally determined using FTIR spectroscopy (see Sec. 3 of the Supplementary Information). The dominant feature is the phonon at 𝜔 , but we include additional 13 phonons for a faithful description of the optical response (phonon parameters are listed in Supplementary Table 1). For each 𝑖-th phonon mode, the equation of motion in the time domain is given by35

𝑄̈ +2Γ𝑄̇ +𝜔2 𝑄 = 𝑍 𝐸 (𝑡) (5)

Where Γ is the damping rate, 𝜔 is the phonon frequency and 𝑍 = 𝜔 √𝜖 (𝜀 −𝜀 ) is the phonon effective charge, being 𝜀 and 𝜀 the dielectric constant at the low and high-frequency limit respectively and 𝜖 the vacuum permittivity.

Lattice equations and the electric field component from Maxwell’s equations in media are related by displacement field as:

𝐷 = 𝜖 𝜀 𝐸 +∑𝜔 √𝜖 (𝜀 −𝜀 ) 𝑄 (6)

By FTDT calculation we numerically solve the above equation at points of the grid through leapfrog time-step scheme. In the simulation, the THz electric field 𝐸 (𝑡) = 𝐸 cos(𝜔 𝑡)𝑒 − 2𝜎2 is a Gaussian field having a width 𝜎 = 100 fs and a center frequency 𝜔 = 5 THz. The plot of the spatiotemporal propagation of the phonon at 𝜔 is show in Fig. 4b.

To quantify the linear and nonlinear optical responses due to phonon displacement, first-principles calculations were perfomed in the framework of density-functional perturbation theory (DFPT) implemented using PHonon package49 in Quantum ESPRESSO50,51, with PBE SG15 Optimized Norm-Conserving Vanderbilt (ONCV) pseudopotentials52,53.

DFPT calculations were conducted within the structural unit cell [42 atoms, R3c (No. 161) space group] of the BBO crystal at room temperature has been taken from Ref. 54. For the Brillouin zone sampling, 4x4x4 k-point grid and a plane-wave energy cutoff of 80 Ry have been used and the structural parameters have been relaxed until the force acting on each atom was less than 0.0025 eV/Å and the pressure was below 0.5 kbar. Though phonon energies and corresponding displacement vectors at the Brillouin zone center, the macroscopic dielectric tensor, and Born effective charges were computed.

The refractive index ellipsoid were computed as the square-root of the macroscopic dielectric tensor, calculated with DFPT in the long-wavelength limit (𝑞 = 0). This is valid in non-absorbing materials, in which the refractive index dispersion is small, as in the case of BBO in the optical optical region.

At zero displacement along the phonon coordinates (𝑄 = 0), the DFPT-calculated dielectric tensor in the crystal frame is:

ξ = ( ξ0 0 0 / 0 ξ0 0 / 0 0 ξ0 ) (7)

which corresponds to a refractive 𝑛 = √ξ = 1.705 and 𝑛 = √ξ = 1.705. Previous DFPT studies have reported a refractive index of 𝑛 = 𝑛 =1.692 at 800 nm, see Ref.32, which is within 0.5% of our DFPT-calculated value. This close agreement validates the use of the long-wavelength approximation in our approach. Additionally, the small discrepancy might by due to the fact our DFPT method, Ref.49, includes local-field effects and goes beyond the independent-particle approximation employed in Ref. 32.

Experimentally, the refractive index has been measured to be 𝑛 = 1.66, which is within 3% of our calculated results.

The dielectric tensor as a function of the phonon amplitude, was computed with the same approach but assuming structures displaced along the relevant IR-active phonon mode displacement vectors.

The dielectric tensor as a function of the phonon amplitude has the following form:

ξ(𝑄) = ( ξ ξ 0 / ξ ξ 0 / 0 0 ξ ) → 𝑅(𝜋) ( ξ̃ 0 0 / 0 ξ̃ 0 / 0 0 ξ ) (8)

which diagonalized with a 𝑥’−𝑦’ plane rotation around the 𝑧’ axis, as in the analytical model (see Sec. 6 of the Supplementary Information). For the dielectric tensor in Eq. 7, the refractive indeces as a fuction of phonon amplitude were calculated.

For the electronic density calculations in Fig. 4b, a plane-wave kinetic energy cutoff of 320 Ry has been used. To evaluate the nonlinear optical coefficients of BBO (see Supplementary Fig. 6) finite-field approach was used55. In this method, the electric field is explicitly included in the DFPT calculations, and the nonlinear coefficients are obtained numerically by computing the second derivative of the polarization with respect to the applied electric field.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Present manuscript reports the observation of THz field modulated SHG effect in the BBO crystal that has potential applications for fast optical control in photonic technologies. The modulated response was convincible attributed to the THz electric filed induced birefringence, i.e. electro-optic effect, which in turn modulates the SHG phase matching conditions through the phase mismatchΔk. Since the local field of THz waves could be very high and close to resonance with certain vibration mode that make major contributions to the EO coefficients, the enhancement of the EO effect could be substantial. These findings are novel and have strong tendency towards applications hence deserve to be published in NC. Before being accepted, there are following concerns or suggestions to the authors:\n1. Estimation of the THz field, a number of 8.5 MV/cm was mentioned in the text. Presumably that’s for the whole spectral range, how much it will be around 4.3 THz (e.g., ±0.1 THz) it could be potentially important to estimate the EO coefficient close to phonon resonance.\n2. In fig. 1d and f and throughout the text, the authors seemed to believe that the spectrum of modulated SHG signal has a peak at exact the phonon frequency of 4.3 THz. However, at the exact resonance the penetration depth of the THz field could be very limited (perhaps ~1 um since the absorption coefficient is normally on the order of 10^4 cm^{-1} at resonance). And above the resonance, contribution from the phonon mode will go negative to the EO effect. Actually, it looks like that the peak position of the SHG spectrum (Fig. d and f) is slightly below that of the peak in  ε_2, please confirm that.\n3. Related to 2, the  ε_2 curve and its peak positions were obtained by fitting the FTIR absorption spectrum of thin BBO crystal (Fig. S3). However, from 4 to 8 THz the transmittance is practically zero. How much reliable will be the phonon frequencies spectral weights determined this way in that range (e.g. 4.3 THz-144.83 cm^{-1}, and 188.78, 245.15 cm^{-1}).\n4. In the text, line 268, the refractive indices by DFT calculation were given as n_x = n_y = 1.705, yet in the in line 374, n_x = 1.76 and n_y = 1.76 were listed. At what wavelength of those values are referred to, it is about 1.66 at 800 nm from measured values.\n5. The vibrational electro-optic coefficient r_{22} was estimated to be 0.0014 Å^{-1}(amu)^{-1/2}, which is not in the standard unit. It must be converted to the SI unit (m/V) for comparison to the electronic contributions or the contributions from all the phonons at off-resonance frequencies.", "reviewer_2": "This paper reports on a TFISH experiment, in which the TeraHertz modulation of the second order optical nonlinearity of the material is achieved through excitation of an infrared active phonon, which can potentially provide greater response than nonlinear excitation.\nThe work is potentially interesting and although it is not so new (there has been work in phononics probed with SHG, for example in the study of ferroelectrics or cuprates). Nevertheless, the nonlinear optical implications have not been discussed in detail.\na suitably revised version of this manuscript would in principle be appropriate for Nature Communications.\nYet, there are some major flaws in the present discussions that the authors must address. I must say that the authors do not strike me as having a particularly deep understanding of the microscopic physics behind their experiment. I suggest that the authors think about this work more deeply and consult with experts in phononics, Nonlinear optics and relate fields.\nI list some of the most important flaws I have found when studying this manuscript.\n1) The material is not centrosymmetric, hence the mode cannot be of Eg symmetry. In a non centrosymmetric medium one cannot separate modes in even and odd, since they are all infrared and Raman active. Parenthetically, if the mode were of Eg symmetry it would not be Infrared active because it would be even. I think that the authors should take a look at their nomenclature and correct accordingly.\n2) The phonon mode being driven is being displaced by the electric field. Hence, the relevant quantity is a Polarization P, which is proportional to Chi and hence to epsilon1 (not epsilon2). The figure compares the response to epsilon2, which is a measure of coherent dissipation, and not of the induced phonon dipole. They should plot the loss function of the material (1/Im(epsilon1)) and compare the observed oscillation frequencies to existing phonons. Please correct accordingly.\n3) More importantly. What are the authors measuring in their experiment? There are at least two contributions to the second harmonic signal when modulated by an IR phonon QIR. One contribution, which is a homodyned response, will be proportional to the square of the SH field and therefore one expects oscillations at twice the fundamental frequency. There will also be a heterodyne contribution, which will presumably rely on the interference between static SHG form the BBO crystal and the induced lattice oscillations. Have the authors studied the effect in detail and are they sure which made they are looking at? Oscillations at 4.3 THz could be a homodyned SHG from a 2.15 THz phonon or heterodyne oscillations from a 4.3 THz phonon. The best way to solve this issue would be to add a local oscillator to make sure that the heterodyned component dominates (se for example Mankowsky et al. PRL 118, 197601 (2017))\n4) Results of figure 4. I was confused in reading about this part of the paper. If the mode is an E mode, it means that it is doubly degenerate and I imagine that there will be two sets of mode displacements (lets call them E1 and E2) which will respond to the THz drive. If I look at figure 1 I imagine that a second set of vibrations will be found at some angle. Can the authors clarify this aspect? Maybe this is already contained in their analysis and hence this simply needs to be made clearer in the text.\n5) The authors should maybe mention that this effect is in fact a third order (Chi(3)) effect, in which two photons form the mid infrared pulse are mixed with one photon from the THz modulator. Indeed, the modulated SHG will be upshifted and downshifted by one phonon, which may be important for the phase matching conditions. In this sense, the propagation may have to include polariton equations, since on and off resonant excitation of a phonon will lead to a modulating field that propagates at vastly different group and phase velocities.\nIn summary, this is a paper with some potential, but in its current form, the paper is unpublishable. I am happy to help cleaning this paper up, but the manuscript requires a lot of work."}, "author_response": "We thank the reviewer for time spent evaluating our manuscript, for constructive feedback, and for their provisional recommendation for publication in Nature Communications. We are pleased that the reviewer finds our observations of all-optical control of phase-matched SHG in BBO to be both novel and relevant for future photonics applications.\nIn response to both reports, we have revised the manuscript substantially, which we believe has improved its clarity and overall quality. Specifically in response to the first report, we (i) estimated the magnitude of the vibrational electro-optic coefficient in SI units experimentally using broadband and narrowband THz pumping, which required additional measurements, and computationally via density-functional perturbation theory (DFPT); (ii) compared phonon frequencies obtained from FTIR with literature values from THz-TDS and Raman spectroscopy; (iii) provided additional details on the refractive-index calculations performed with DFPT; and (iv) compared the frequency-dependent features of the SHG modulation dynamics with the frequency dependence of the vibrational electro-optic coefficient.\nWe believe that these revisions and additional information, included in direct response to the referee’s suggestions, substantially clarify and strengthen our findings. Below, we provide a point-by-point response to the reviewer’s comments (in italics) and list the corresponding changes made in the manuscript.\n[... full rebuttal text provided in source ...]"}, {"round": 2, "reviewer_comments": {"reviewer_1": "After carefully reading there revised manuscript, replies to comments and SI, it is found that all my previous concerns were considered and corresponding revisions were made satisfactorily. Therefore I have no objection to accept present manuscript.\nps. there is a missing reference R5 in the SI, please update it.", "reviewer_2": "The authors have done a good job and convincingly rebutted the vast majority of my concerns. There is no point in litigating a few things that I am still not 100% clear about. These are marginal in the context of the claims made here and would add unnecessary delays.\nI am happy for this paper to be published as is and I would like to congratulate the authors for the nice results and to thank them for their comprehensive clarifications."}, "author_response": ""}]

Decision letter