Multicycle terahertz pulse generation by optical rectification in LiNbO$_3$, LiTaO$_3$, and BBO crystals

We report multicycle, narrowband, terahertz radiation at 14.8 THz produced by phase-matched optical rectification of femtosecond laser pulses in bulk lithium niobate (LiNbO$_3$) crystals. Our experiment and simulation show that the output terahertz energy greatly enhances when the input laser pulse is highly chirped, contrary to a common optical rectification process. We find this abnormal behavior is attributed to a linear electro-optic (or Pockels) effect, in which the laser pulse propagating in LiNbO$_3$ is modulated by the terahertz field it produces, and this in turn drives optical rectification more effectively to produce the terahertz field. This resonant cascading effect can greatly increase terahertz conversion efficiencies when the input laser pulse is properly pre-chirped with additional third order dispersion. We also observe similar multicycle terahertz emission from lithium tantalate (LiTaO$_3$) at 14 THz and barium borate (BBO) at 7 THz, 10.6 THz, and 14.6 THz, all produced by narrowband phase-matched optical rectification.


Introduction
Intense, singlecycle, broadband terahertz (THz) sources are essential for many applications including THz-driven acceleration of electrons and protons [1,2], molecular alignment [3], high harmonic generation [4], and material sciences [5]. In particular, femtosecond laser-based optical rectification (OR) in χ (2) nonlinear materials is considered to be one of the most efficient methods for energyscalable THz generation [6]. OR can be highly effective when the group velocity of the laser pulse is matched to the phase velocity of the THz wave in the nonlinear medium-called phase matching. As an OR-based THz source, lithium niobate (LN) is widely used due to its excellent material properties such as high nonlinearities (d 33 = 168 pm/V at 1THz) [7], high transparency at 0.4∼5 µm [8], and well-developed poling techniques [6]. For efficient phase matching in LN, tilted-pulse-front (TPF) schemes can be used to generate intense singlecycle THz pulses [7,9,10,11].
Multicycle narrowband THz sources are also of great interest owing to many emerging applications including waveguide-based electron acceleration [12], coherent X-ray generation [13], resonant pumping of materials [3], and narrowband spectroscopy [6]. Multicycle narrowband THz radiation is often produced by OR in periodically-poled lithium niobate (PPLN) crystals [14,15,16,17,18]. Cryogenic cooled PPLN crystals are also used to suppress strong THz absorption in LN, lately providing a laser-to-THz conversion efficiency up to 0.1% [17]. Another approach is to drive OR with intensity-modulated laser pulses such that the produced THz waveform can follow the intensity envelope of the modulated laser pulse [19]. Other methods include transient polarization gratings [20], TPF planer waveguides [21], and cascaded second-order processes [22].
Recently, we have observed a new type of multicycle radiation at ∼15 THz emitted from a bulk LN crystal when irradiated by femotsecond laser pulses [23]. High-energy THz radiation up to 0.7 mJ has been also produced from a large diameter (75 mm) LN wafer with 80 TW laser pumping [24]. This type of radiation originates from a narrow phase matching condition naturally satisfied in between two phonon resonance frequencies in LN [23,24]. Previously, similar narrowband radiation around 15 THz was produced by difference frequency generation (DFG) in LN, in which two separate laser pulses with different frequencies are mixed to generate THz radiation at the difference frequency [25]. By contrast, our THz generation method is based upon OR of a single laser pulse. This OR process is expected to produce higher THz energy with reducing laser driver's pulse duration. However, certain LN crystals exhibit enhanced THz radiation when driven by highly chirped laser pulses [24], contrary to our understanding of OR. Moreover, in the previous experiments, the radiation spectrum was poorly characterized with THz bandpass filter sets [24] or incompletely studied [23].
In this paper, we present a comprehensive study of multicycle narrowband THz generation around 15 THz from LN crystals. Experimentally, we measure THz field autocorrelation and spectral power under various laser conditions, especially when the laser driver is chirped with third order dispersion. To explain our experimental observation, we carry out numerical calculations on THz generation and propagation in LN. We also describe experimental measurements of chirp-dependent narrowband THz generation from lithium tantalate and barium borate crystals.

Experimental setup
The schematic of our experimental setup is shown in Fig. 1 (a). Femtosecond laser pulses from a Ti:sapphire amplifier operating at 1 kHz are loosely focused onto a LN crystal by a lens with a focal length of 1.5 m. The laser (pump) beam size (3.4∼6.8 mm in 1/e 2 diameter) and fluence (3.2∼28.7 mJ/cm 2 ) are varied by translating the LN crystal along the beam propagation direction and/or controlling the laser energy. The pump pulse provides energies up to 2.6 mJ at a central wavelength of 800 nm with a 30 nm full-width half-maximum (FWHM) bandwidth as shown in Fig. 1(b). In our measurements, x-cut congruent LN crystals of 10 mm × 10 mm × 0.5 mm or 1 mm (thickness) are used for THz generation. The LN crystal is oriented such that its extraordinary axis is parallel to the laser polarization for maximal THz generation.
To decouple output THz pulses from the copropagating pump beam, three optical windows coated with 180-nm-thick indium thin oxide (ITO) are placed after the LN crystal. The ITO window allows high optical transmission (>85% each) in the visible and near-infrared regions with strong reflection (∼80% each) at <15 THz [26]. Any pump leakage after the three ITO windows is completely blocked by a 280-µm-thick high-resistivity (>10 kΩ·cm) silicon (Si) window in the downstream beamline.
The resulting THz pulses are characterized by a lab-built Michelson-type Fourier-transform infrared (FTIR) interferometer combined with a pyroelectric detector (PD) (Spectrum Detector Inc., API-A-62-THz). The incoming THz beam is split and recombined with a variable time delay by a 280-µm thick Si wafer in the interferometer and then focused by a 90 • off-axis parabolic (OAP) mirror onto the PD detector. The THz signal from the PD detector is fed into a lock-in amplifier that is phase-locked to an optical chopper modulating the input laser beam at 10 Hz. A delay scan in the interferometer provides a THz field autocorrelation from which the spectral power can be obtained by the Fourier transform.
The pump pulse duration is varied by tuning the distance between the grating pair in the pulse compressor (see Appendix). This effectively changes the group delay dispersion (GDD) of the pump pulse. Figure 1(c) shows the pump pulse duration as a function of the grating distance. With a Gaussian spectral assumption, the Fourier transform-limited pulse duration is calculated to be τ = 0.44λ 2 0 /(c∆λ) ≈ 31 fs, where λ 0 = 800 nm and ∆λ = 30 nm are the central wavelength and bandwidth in FWHM, respectively. With GDD = 0 fs 2 , the pump pulse duration measured by a single-shot second-harmonic autocorrelator is about 67 fs in FWHM [dotted line in Fig. 1(c)]. This is longer than the transform-limited pulse duration of 31 fs. The difference is explained by uncompensated third order dispersion (TOD) and higher-order dispersion (HOD) of the pump pulse. Those result in an asymmetrical plot of the pulse duration as shown in Fig. 1(c). Figure 1(d) shows a typical THz intensity profile captured at the focus by a room-temperature microbolometer focal plane array (FLIR, Tau 2-336) [27,28]. It shows that the focused THz radiation is confined within a spot size of 160 µm in FWHM diameter.   For a series of GDD values, a two-dimensional (2-D) spectral power plot (color scale) is obtained and shown in Fig. 2(c). Figure 2 clearly shows two types of THz radiation emitted from LN-broadband singlecycle emission at 0∼8 THz and narrowband multicycle emission around 15 THz. Interestingly, Fig. 2(c) shows that the narrowband emission is greatly enhanced with a properly stretched laser pulse duration (GDD = −800 fs 2 , 1,500 fs 2 ) whereas the broadband radiation is maximally produced with the shortest pulse duration (GDD ≈ 0 fs 2 ).
The broadband THz emission has been previously observed and explained by non-phase-matched OR in LN through a second-order nonlinear χ (2) process [15,16,29]. This yields singlecycle THz pulses emitted from both the front and rear layers of LN with a thickness of one coherence length l c = λ THz /(2 |n g − n THz |) for each layer. Here n g = 2.3 is the optical group index of LN at 800 nm, n THz is the refractive index of LN at THz frequencies, and λ THz is the THz wavelength. At 1 THz, n THz = 5.1 [30,8] and λ THz = 300 µm gives l c = 53.5 µm. This dual THz pulse generation explains the fast modulations observed in the broadband THz spectrum in Figs. 2(b) and 2(c). They arise due to interference between two temporally separated THz pulses generated from the front and rear surfaces of the LN crystal.

Characteristics of narrowband radiation at 15 THz
The narrowband multicycle radiation in Fig. 2 peaks at 14.8 THz with a 0.94 THz FWHM bandwidth at GDD = 1,500 fs 2 . Its output energy dependence on the pump GDD is also measured and plotted in Fig. 3. Here the THz energy is measured by the pyroelectric detector (PD) along with a THz bandpass filter (BPF) providing a central frequency of 15 THz placed in front. This allows to measure only narrowband THz radiation around 14.8 THz.
As shown in Fig. 3, the output THz energy peaks at two GDD ranges of 900∼1,600 fs 2 (positive chirp) and −1,200∼ −800 fs 2 (negative chirp). More interestingly, the positive GDD range yields more output THz energy with increasing pump energy and/or LN thickness. Furthermore, the output  THz energy is abnormally suppressed at GDD ≈ 0 fs 2 , where the highest nonlinearity is expected due to the shortest pump pulse duration. Previously, similar results were observed and explained by THz screening and absorption by free charge carriers produced by multi-photon laser absorption in LN [11]. In our experiment, the nonlinear absorption coefficient by free carrier absorption (FCA), α f c , is estimated to 31.63 cm −1 at 14.8 THz under laser conditions of 470 GW/cm 2 peak intensity and 800 nm wavelength [31,32]. The value, however, is much smaller than the intrinsic absorption coefficient α = 1,440 cm −1 at 14.8 THz in LN [23,24]. Therefore, three-photon-absorption followed by FCA is not believed to cause the suppressed THz emission at GDD ≈ 0 fs 2 . Instead, the odd GDD-dependence can be explained by a THz-induced cascaded effect on the pump pulse that has nonzero TOD as will be described in Section 5.
Another interesting feature observed with the narrowband THz radiation is its energy scaling. Figure 4(a) shows the measured output THz energy emitted from the 1-mm-thick LN crystal as a function of the pump energy with the beam diameter of D = 3.4 mm, corresponding to Fig. 3(f). In phase-matched OR, the output THz energy is expected to increase quadratically with the pump intensity, i.e., e THz ∝ |I pump | 2 . This is indicated by the red-dotted line in Fig. 4(a). The observed THz energy, however, increases much faster than the theoretical prediction at the pump energy exceeding 1.2 mJ. This unexpected behavior can be also explained by a THz-induced cascaded effect as will be explained in Section 5.
The resulting laser-to-THz conversion efficiency is shown in Fig. 4(b). The maximum THz output energy of ∼92.6 nJ and efficiency of 3.7 × 10 −5 are achieved. This provides a maximum field strength of 0.4 MV/cm at the focus, estimated from the measured energy, pulse duration (∼470 fs), and beam spot size (∼160 µm).

Narrowband THz radiation from LT and β-BBO crystals
Multicycle narrowband THz emission is also observed from lithium tantalate (LT) and beta-barium borate (β-BBO). These two nonlinear materials including LN are commonly used inorganic χ (2) nonlinear crystals and have a trigonal structure with point group 3m. LT and β-BBO are also tested for GDD-dependent THz generation, and the result is shown in Fig. 5. Figures 5(a) shows THz autocorrelation signals obtained from a 0.5-mm-thick LT crystal at two different pump GDD values of 0 and −1,100 fs 2 . The corresponding spectra are shown in Fig. 5(b). From a GDD scan from −4,200 fs 2 and 4,200 fs 2 , a 2-D plot of THz spectrum (color scale) is obtained and displayed in Fig. 5(c). Figures 5(d-f) shows experimental data obtained with a 0.1-mm-thick β-BBO crystal. For both measurements, the pump energy fluence is fixed at 28.7 mJ/cm 2 . Clearly, both crystals exhibit multicycle THz waveforms and consequent narrowband THz emission. In the case of LT, its narrowband emission is centered at 14 THz, consistent with the refractive index [33] and phase matching condition for LT. Similar to the LN crystals tested before, LT shows both singlecycle (broadband) and multicycle (narrowband) radiation depending on the pump chirp condition.
The spectral power 2-D plot (color scale) of β-BBO shown in Fig. 5(f) exhibits narrowband emission at 7 THz, 10.6 THz, and 14.6 THz. Our result is consistent with a previous study reporting narrowband emission at 4.3 THz, 7 THz, and 10.6 THz [34]. Interestingly, 4.3-THz emission is not seen in our experiment possibly due to its weak spectral power. Instead, our experiment reveals new narrowband emission at 14.6 THz. We note that this observation was possible due to our FTIRbased detector's capability of measuring high-frequency THz emission beyond 10 THz. Contrary to commonly used electro-optic sampling (EOS) methods, in our scheme the detection bandwidth is not limited by the laser pulse duration or THz absorption/dispersion in the electro-optic (EO) material. Also, our detector is independent from the source and not affected by any pump chirp. This allows us to characterize the radiation spectrum without being distorted or restricted by the pump GDD.

Theoretical background
The narrowband THz radiation observed in Figs. 2 and 5 is fundamentally characterized by phasematched (n g = n THz ) OR. Here n g and n THz are the group and refractive indices at the pump and THz frequencies, respectively. For example, Fig. 6(a) shows the refractive index n THz of congruent LN as a function of frequency [8,35,36,37] (see Appendix). The optical group index n g = 2.3 at 800 nm is also plotted in Fig. 3(a) with a gray dotted line. It shows that the phase-matching condition, n g = n THz , is satisfied at 14.8 THz in between two strong transverse-optical (TO) phonon resonance frequencies in LN (7.4 THz and 18.8 THz). These resonance frequencies are clearly shown by the absorption coefficient α plotted in Fig. 3(a) with a red solid line. Note that there are two additional phase-matched frequencies occurring at 8.3 THz and 19.3 THz. However, little or no emission is expected at those frequencies because of their strong absorption in LN. Figure 6(b) shows the effective nonlinear coefficient d eff in the extraordinary direction of congruent LN [35,36]. It shows d eff reaches its local minimum value of 10 pm/V at 11 THz while peaking to 1,424 pm/V and 543 pm/V at the two phonon resonance frequencies, 7.4 and 18.8 THz, respectively. At frequencies <7.4 THz, d eff asymptotically approaches 168 pm/V, which is consistent with the reported value at 1 THz in Ref. [7]. At 14.8 THz, d eff = 82 pm/V, which is still sufficiently large to generate strong THz radiation [7,32].
We emphasize this type of narrow phase matching can occur in many nonlinear crystals including LT and BBO. In general, the refractive index n THz changes so large in between two phonon resonance, and often there exist one or multiple frequencies at which n THz becomes equal to the optical group index of refraction n g . Absorption is also expected to be relatively low in between phonon resonance frequencies. In addition, the phase-matched frequency can be tuned by varying the pump laser wavelength although it provides a narrow tuning range. For example, optical pumping at 0.4∼1.9 µm in LN can yield phase-matched emission at 14.4∼15.7 THz.
At the phase-matched frequency of 14.8 THz in LN, the absorption coefficient reaches its minimum value of α = 1,440 cm −1 as shown in Fig. 6(a). This value, however, is still large enough to attenuate the emission significantly. In phase-matched OR, the effective length for maximal THz generation is generally given by L eff = (α/2 − α L ) −1 ln [α/(2α L )] ≈ 160 µm, where α L = 0.0078 cm −1 is the laser absorption coefficient near 800 nm [24,38]. This means that only a layer of 160 µm thickness from the front surface of the LN crystal can maximally generate 14.8 THz radiation when the incident pump pulse is unchirped. In this case, a thin LN crystal is best suited for efficient THz generation as demonstrated in our previous experiment [24]. For thicker ( L eff ) crystals as in the current experiment, negatively-chirped pump pulses are generally preferable as they can be compressed with propagation and effectively produce THz radiation from near the rear surface of the LN crystal.
Interestingly, our experiment shows that both positive and negative chirps yield enhanced THz emission as shown in Figs. 2 and 3. So it is the stretched pulse duration, not chirp, that matters in our narrowband THz generation. In general, too long pulses are not good as they do not provide enough bandwidths to generate 14.8 THz radiation by OR. For example, at GDD = 1,600 fs 2 that yields enhanced 14.8 THz radiation, the corresponding pump pulse duration is estimated to ∼130 fs in FWHM from Fig. 1(c). This is about twice longer than the period of 14.8 THz radiation, 68 fs. This implies that the pump pulse must have certain intensity modulations (or pulse splitting) within the pulse envelope to provide a sufficient bandwidth to generate 14.8 THz. Such modulations can be initially made by applying a proper combination of GDD and TOD onto the pump pulse.
Pump intensity modulations can also arise and be amplified from a nonlinear process. For instance, a pump pulse propagating through LN can be distorted by the THz field it produces via a linear electro-optic (or Pockels) effect [39,40]. This is a second-order χ (2) nonlinear process and can lead to spectral shifts, broadening, and modulations of the pump pulse [10,39,38]. Thus, in order to explain the peculiar GDD dependence of the narrowband radiation, one needs to include the Pockels effect, as well as laser-THz dispersion and absorption in conducting numerical simulations.

Theoretical model
In order to simulate THz generation by OR, we first present the electric field of a laser (pump) pulse by using a Gaussian envelope shape given by where ω 0 is the center frequency, ∆ω = 4 ln 2/∆t is the FWHM spectral bandwidth, and ∆t is the pulse duration in FWHM. The amplitude E 0 is determined by the pump fluence F as To account for chirped pump pulses, the spectral phase of the pump pulse is expanded in a Taylor series about ω 0 as where the fourth and higher-order terms are neglected for the sake of simplicity. The input pump pulse is obtained in the frequency domain as where F T denotes the Fourier transform. Then we solve one-dimensional (1-D) coupled forward Maxwell equations (FME) [41,29,24] self-consistently for both THz and optical pump pulses in the frequency domain as where E T and E P are the electric fields of THz and optical pump, respectively, which propagate in the coordinate ξ = z − ct/n g moving at the group velocity of the pump. The first terms on the right hand side of both equations (α terms) correspond to absorption of both fields. The second terms correspond to material dispersion D (ω T , ω) = ω (ω T , ω) [n (ω T , ω) − n g ] /c. The third term in Eq. (5) represents the second-order nonlinear polarization due to OR, a source term for THz radiation. The third term in Eq. (6) describes the Pockels effect on the pump pulse induced by the produced THz field. The last term in Eq. (6) corresponds to self-phase modulation (SPM) of the pump pulse via the Kerr effect, where the third-order nonlinear susceptibility χ (3) is derived from the nonlinear refractive index n 2 = 3χ (3) / 4c 0 n 2 0 . Here n 2 = 10 −6 cm 2 /GW is used in our calculation [42]. In the simulation, we ignore THz-induced pump modulation via the χ (3) nonlinear process. This is because the χ (3) -based pump phase shift, ∆ϕ (3) , is much smaller than ∆ϕ (2) induced by the Pockels effect, i.e., ∆ϕ (3) /∆ϕ (2) 1 [39]. The numerical integrals in Eqs. (5) and (6) are solved by a 4th-order Runge-Kutta method with spatial resolution of 500 nm in order to achieve required numerical convergences. Note that the model used here considers only 1-D space along the propagating direction z. This is justified for a large pump beam size, where any transverse beam effects such as self-focusing and diffraction can be ignored.
In the simulation, the input pump pulse is assumed to be Gaussian with a 30-nm FWHM spectral bandwidth at 800 nm to be consistent with our experimental condition. Experimentally, nonzero third order dispersion (TOD) arises from two sources; one is from the compressor's tuning for GDD control (see Appendix). The other one comes from the amplifier itself but not properly compensated by the compressor. This residual TOD is implicitly shown in Fig. 1(c) by the asymmetric slope. Here the total TOD can be expressed as TOD = TOD g + TOD i , where the subscripts g and i denote "grating" and "initial (residual)". In the simulation, we used TOD g ≈ −2(fs) · GDD and TOD i = 3,800 fs 3 for our compressor system (see Appendix).

GDD-dependent THz spectrum
We compare our simulation results with the experimental ones shown in Fig. 2. Figure 7(a) shows a simulated THz spectral power plot (color scale) obtained from a 1-mm-thick LN at laser fluence of 13.6 mJ/cm 2 . Clearly, it reproduces multicycle narrowband THz radiation around 15 THz. Also, it is most efficiently produced at GDD = 1,600 fs 2 and −800 fs 2 . At GDD ≈ 0 fs 2 , it is greatly suppressed while the broadband radiation at 0∼8 THz is maximally enhanced. This is in good agreement with our experimental results.
For comparison, we repeated the simulation with a chirped pump with TOD i = 0 fs 3 to better understand the role of TOD in THz generation. The result is shown in Fig. 7(b). Interestingly, it shows the narrowband emission singly peaks at GDD = −800 fs 2 , and the maximal THz energy increases by 18.9% compared to Fig. 7(a). This does not agree with our experiment but is consistent with a general OR process, in which the radiated THz energy increases with decreasing pump pulse duration at fixed laser energy. Here a negatively chirped pump pulse is more favorable for THz generation because it can be compressed as it propagates through LN that possesses positive material dispersion. In addition, we repeated the simulation without including the Pockels effect but keeping the original TOD-the result is not shown in Fig. 7. In this case, maximal THz radiation occurs at GDD = −400 fs 2 , but the generated THz energy decreases to 62% compared to the first simulation result shown in Fig. 7(a). All these suggest that the strange GDD dependence of 15 THz observed in our experiment can be attributed to a combined action of both the Pockels and TOD effects.

Evolution of THz and laser fields with propagation
For a detailed understanding of multicycle THz generation, we simulated the evolution of THz electric fields and laser intensity envelopes as they propagate through LN. Here all simulation parameters are the same as in Fig. 7(a). First, Figure 8(a) shows the cumulative THz energy (color scale) plotted as a function of the initial GDD (vertical axis) and the propagation distance z (horizontal axis). Here the pump TOD varies as TOD = −2(fs) · GDD + 3,800 fs 3 . Two line-outs at GDD = −800 fs 2 and 1,600 fs 2 are plotted in Fig. 8(b). In the case of GDD = −800 fs 2 , the generated THz energy slowly increases and then decreases beyond z = 0.6 mm. However, with GDD = 1,600 fs 2 , the energy noticeably increases after z = 0.3 mm. Figure 8(c) shows the THz electric fields (black lines) calculated at z = 0, 0.3, 0.6, and 0.9 mm with an initial GDD value of 1,600 fs 2 . Also co-plotted (red lines) are the temporal derivatives of the pump intensity profiles, −dI (t, z) /dt. For reference, the input pump intensity profile I(t) is shown (blue line) at z = 0 mm. Clearly, it exhibits damped intensity modulations on its tail due to nonzero GDD and TOD. Initially, this type of intensity modulations is not best suited for multicycle THz generation because its oscillation frequency is chirped and not fully matched to 14.8 THz. However, with a propagation, the pump envelope becomes synchronously modulated by the co-propagating 14.8 THz field via the Pockels effect. This results in a series of pump pulses (pulse splitting) separated by the THz period, which in turn drives OR resonantly to generate multicycle 14.8 THz radiation. This is evidently shown by the time derivative of the pump intensity envelope in Fig. 8(c) (red lines). Its Fourier spectral power at various z is plotted in Fig. 8(d). At z = 0.9 mm, it peaks at ∼15 THz. It also produces two side bands. The left one is responsible for singlecycle broadband THz radiation at <10 THz, whereas the right (weak) one is believed to be the source of 20∼23 THz shown in Fig. 7(a) although it was not observed in our experiment.
Also, the input I 0 (ω) and transmitted I t (ω) pump spectra are computed and plotted in Fig.  8(e) along with experimentally measured ones. The transmitted one corresponds to a 1-mm-thick LN crystal pumped at 22 mJ/cm 2 with GDD = 1,500 fs 2 . As shown in Fig. 8(e), both simulated and measured pump spectra show spectral modulations and small frequency shifts. For negative GDD values, blue-shifted spectra are observed (not shown here). Figure 9 shows the simulated output energy of multicycle THz radiation as a function of the pump GDD for 0.5-mm and 1-mm thick LN crystals. It shows that more output THz energy is produced at negative pump GDD values with a thinner (0.5 mm) LN. With increasing laser fluence and LN thickness (1.0 mm), however, the peak moves to the positive GDD side, consistent with our measurement in Fig. 3. Also, our simulation provides a laser-to-THz conversion efficiency of 2.1 × 10 −5 from a 1-mm-thick LN pumped at 13.6 mJ/cm 2 . This is in reasonable agreement with our experimental value ∼10 −5 obtained under similar laser fluence conditions in Fig. 4(b).

Conclusion
In conclusion, we have demonstrated efficient multicycle narrowband THz generation at 14.8 THz from bulk LN crystals by using chirped optical pump pulses. The generation mechanism is explained by phase-matched OR naturally occurring in between two phonon resonance frequencies in LN. In our experiment, we have observed enhanced multicycle THz emission when the pump pulse is highly chirped. This anomalous behavior is also observed in our numerical simulations and explained by resonant intensity modulations of the pump pulse by self-produced THz fields through the Pockels effect. The modulated pump pulse can in turn produce multicycle THz radiation efficiently with propagation. This cascaded effect becomes highly efficient when the pump pulse is pre-modulated with proper second and third order dispersion. We also report the first demonstration of multicycle THz pulse generation at 14 THz from LT and 14.6 THz from β-BBO crystals. This new type of narrowband phase matching scheme is universal and can be applied to many nonlinear materials with potential to provide robust, efficient, and tabletop multicycle THz sources.

Dispersion control in a dual grating compressor
In this experiment, the laser pulse duration is controlled by adjusting the separation between a pair of diffraction gratings in the pulse compressor. Right after the compressor, the second-order and third-order spectral phases of the laser pulse are given by [43] where L is the perpendicular distance between the two parallel gratings, θ in is the incident laser angle, and d is the grating groove spacing. In our compressor system, we use 1,500 grooves/mm gratings at θ in = 58.4 • around a central wavelength of λ = 800 nm. The second-order φ (2) and third-order φ (3) are generally referred to GDD and TOD, respectively. Figure 10 plots how GDD and TOD vary as a function of the grating separation L in our pulse compressor. For simplicity, Eq. (8) can be rewritten as TOD = A · GDD, where A is a constant that depends on both the laser wavelength and incident angle, i.e., A = A(λ, θ in ). In our case, it is estimated to A ≈ −2 fs, and the shortest pulse duration is achieved with GDD = −1.229 ps 2 . Thus, a "negative" or "positive" chirp is defined when GDD becomes less or greater than −1.229 ps 2 , respectively. 7.2 Calculation of (ω) and d eff (ω) of lithium niobate The complex dielectric function (ω) and nonlinear coefficient d eff (ω) of congruent LN are given by [37,36] ( where S j , ω j , and Γ j are the oscillator strength, the resonance frequency, and the width of the jth transverse-phonon mode; ∞ = 4.6 is the frequency-independent bound electronic dielectric function; d e and d Qj are the electronic and ionic nonlinear coefficients, respectively. The real (n) and imaginary (k) parts of the square root of the complex dielectric function, (ω) = n + ik, are related to the refractive index n (ω) = n and absorption coefficient α (ω) = 2kω/c. All parameter values necessary to calculate Eqs. (9) and (10) are listed in table 1.

TOD effects on multicycle THz pulse generation
To investigate the TOD effects, we repeated the simulation with various pump TOD values. In our experiment and simulation, the pump TOD varies with GDD as TOD = TOD g + TOD i = −2(fs) · GDD + TOD i . Figure 11 shows 2-D plots of narrowband THz energy (color scale) obtained with TOD i of −4000, −2000, 0, 2000, and 4000 fs 3 . The pump fluence is set to 13.6 mJ/cm 2 . For comparison, Fig. 8(a) is obtained with TOD i = 3,800 fs 3 . As shown in Fig. 11, the output THz energy strongly depends on both GDD and TOD. With TOD i = 0 fs 3 , efficient THz generation occurs near GDD = 0 fs 2 where the pump pulse duration remains relatively short. In this case, the energy rapidly increases but also quickly drops due to large THz absorption. Due to normal material dispersion in LN (368 fs 2 /mm at 800 nm), more negative GDD is necessary to keep the pump pulse duration short with increasing propagation distance (or LN thickness). This is why the color plot in Fig. 11(c) has a single stripe with a negative slope. With large TOD i values (either positive or negative), optimal THz generation occurs with relatively large GDD values (positive or negative) depending on the TOD i sign and the propagation distance z. In this regime, the pump intensity envelop can be pre-modulated by a proper combination of GDD and TOD, and the modulation can be resonantly amplified through the cascaded Pockels effect.