1. INTRODUCTION

Single- to few-cycle (broadband) high energy terahertz pulses have many promising applications such as spectroscopy [1], strong field terahertz physics [2,3], particle acceleration [4], electron spin manipulation [2], and phonon resonance studies [5]. All of these applications require well-characterized terahertz fields.

There are many possible ways to generate terahertz radiation. Free electron lasers and synchrotron radiation have a high degree of tunability and are capable of delivering high peak-power coherent terahertz pulses [6]. Gyrotrons, based on the principle of electron cyclotron radiation, are able to generate watt-to-megawatt-level terahertz continuous wave radiation [7] at low terahertz frequencies (0.3–1.3 THz) [8,9]. These devices, however, have limited accessibility to the larger scientific community and can be difficult to synchronize to laser sources with high (femtosecond) precision.

Alternatively, single- to few-cycle terahertz generation, based on table-top optical laser systems, brings the advantages of high accessibility and intrinsic synchronization, but suffers from limited optical-to-terahertz conversion efficiency. Although various schemes exist for optical-to-terahertz conversion, difference frequency generation using the tilted pulse-front (TPF) method has proven one of the most useful. Using the “tilted-pulse-front” technique, first proposed and demonstrated by Hebling et al. [10], pulse energies in the millijoule range can be reached [11]. The ease of this setup, the high pulse energies, and the controllability of the terahertz properties have made this last approach ubiquitous for high-field applications.

However, the noncollinear geometry of the phase matching and the spatial asymmetry of the interaction, in combination with the cascading effect, result in terahertz beams with nonuniform spatial distribution. This nonuniformity due to the $ {{\rm LiNbO}_3} $ (LN) crystal geometry was explained by Tóth et al. based on a 1D model neglecting nonlinear effects [12]. This effect was also studied by Bakunov et al. [13,14] via a ${\rm 2D} + 1$ numerical model. However, this model does not include the back conversion of the terahertz to the optical pump (OP), i.e., the cascading effect. Consequently, the effective length and the conversion efficiency are overestimated. Later on, the interaction between the optical pump and the terahertz pulse was included into the ${\rm 2D} + 1$ model along with a one lens imaging system by Ravi et al. [15]. The one-lens imaging system has larger imaging errors and induces more terahertz divergence [16], compared to the telescope imaging system. A hybrid-type terahertz generation setup omitting the imaging system was proposed, reducing the imaging errors even further [12]. In order to fully investigate the nonuniform spatial distribution of the generated terahertz pulse, a robust ${\rm 3D} + 1$ numerical tool, including the cascading effect and the telescope imaging system, is necessary. Additionally, we systematically compare ${\rm 1D} + 1$, ${\rm 2D} + 1$, and ${\rm 3D} + 1$ numerical models regarding their practical relevance for the first time to our knowledge.

The numerical tool is based on the fast Fourier transform beam propagation method (FFT-BPM) [17] and split-step Fourier method. The combination of these two methods reduces computational cost compared to the finite difference time-domain (FDTD) method, which is very accurate but requires a massive computational effort. In Section 2, we show analytically that the higher-order dispersion introduced by the grating can cause spatial inhomogeneity of the generated terahertz pulses. In Section 3, the numerical model including the nonlinear interactions that lead to further spatial inhomogeneity is introduced. In addition, the differences of the ${\rm 1D} + 1$, ${\rm 2D} + 1$, and ${\rm 3D} + 1$ models are presented. Finally, in Section 4, the dependence of optical pump beam sizes on the generated terahertz fields is discussed.

2. THEORETICAL MODEL

The tilted-pulse-front setup modeled and simulated is shown in Fig. 1. Often used nonlinear materials are LN, CdTe, GaAs, GaSe, GaP, and ZnTe. Here, we focus on LN due to its large second-order nonlinear coefficient, high damage threshold, and easy accessibility. The results, however, can be extended to other materials. In this paper, we focus on the impact of the OP beam size on the spatio-temporal properties of the generated terahertz beam.

 

Fig. 1. Illustration of the simulated tilted-pulse-front setup. The optical pump pulse is noted by OP, and the $ {{\rm LiNbO}_3} $ crystal is represented by LN. The OP propagates along the $ {z^\prime} $ direction. The $ x - y - z $ coordinates denote the pulse-front-tilt frame (terahertz frame) inside the LN crystal. The $ y $ and $ {y^\prime} $ axes are equivalent. The optical system is chosen such that the image of the grating is parallel to the pulse-front-tilt inside the LN [18].

Download Full Size | PDF

First, we show analytically that, regardless of the noncollinear phase matching and the prism geometry, the second-order dispersion generated by the grating can already cause inhomogeneity of the terahertz pulses [see Eq. (4)]. We start from the grating relation $ \sin ({\theta _1}) + \sin ({\theta _2}) = 2\pi c/\omega d $, where $ {\theta _1} $, $ {\theta _2} $ are the incidence and output angles with respect to the normal to the grating surface, respectively. The optical frequency is represented by $ \omega $ while $ d $ denotes the grating groove period. Assuming that the input OP is collimated ($ d{\theta _1}/d\omega = 0 $), one can get Eq. (1):

One can obtain the angular dispersion $ \Delta {\theta _2} $ with respect to the OP propagation direction, by performing a Taylor expansion to the second order. By inserting Eq. (1) and setting the first-order angular dispersion $ {F_1} = - 2\pi c/(\omega _0^2d\cos ({\theta _{02}})) $ where $ {\omega _0} $ is the center frequency and $ {\theta _{02}} $ is the grating output angle at the center frequency, we obtain

The second-order dispersion is denoted by $ {F_2} $. After the propagation through the telescope system, the angular dispersion is magnified by a factor of $ - {f_1}/{f_2} $. Thus, the angular dispersion becomes $ - \Delta {\theta _2}{f_1}/{f_2} $. Note that the transverse k-vector $ {k_{{x^\prime}0}}(\omega ) = - \Delta {\theta _2}{f_1}{\omega _0}/(c{f_2}) $ remains the same before and after entering the LN crystal due to the Fresnel law. As a result, by assuming _SSS, where $ n $ is the refractive index, the OP electric field inside the LN crystal is given in Eq. (3) in the $ {x^\prime} - {y^\prime} - {z^\prime} $ coordinates. In Eq. (3), the pulse duration is $ \tau = {\tau _0}/\sqrt {2\log 2} $, $ {A_0} $ is the electric field amplitude, and $ \sigma _x^\prime,\sigma _y^\prime$ are beam waists (1/e) of the OP at the incidence surface of the LN crystal in $ {x^\prime} $ and $ {y^\prime} $ dimensions, respectively. Equation (3) is as follows:

The imaging system was chosen such that the image of the grating is parallel to the pulse-front-tilt plane inside the nonlinear crystal [18]. The expression of the OP electric field given in Eq. (3) is only valid at the imaging plane of the telescope. It is very important to notice that the higher-order angular dispersion ($ {F_2} $) induced by the grating leads to temporal broadening of the OP [19,20], as it adds the nonlinear phase [see Eqs. (2) and (3)]. If the second-order angular dispersion is neglected [$ {F_2} = 0 $ in Eq. (2)], the OP pulse duration reduces to the transform limited case at each spatial point along the $ {x^\prime} $ dimension.

Using Eq. (3), the second-order polarization that is responsible for the terahertz generation can be expressed as in Eq. (4), where $ \gamma $ is the angle between the terahertz pulse propagation direction and the OP propagation direction, $ {\chi ^{(2)}} $ is the second-order nonlinear susceptibility, and $ \Omega $ is the terahertz angular frequency. The group refractive index at the center frequency of the optical pump is denoted by $ {n_g} $. Equation (4) is as follows:

The last two exponential phase terms in Eq. (4) represent the phase-matching condition:

The term in the third exponential ($ {( {{F_2}/{F_1}} )^2}16{x^{\prime 2}}n_g^2$ $\tan {(\gamma )^2}/{\tau ^4}{c^2} $) is due to the second-order angular dispersion that has been discussed in more detail in Ref. [21]. It can be seen that the second-order angular dispersion leads to a spatial dependence of the generated terahertz bandwidth along the $ {x^\prime} $ dimension. At the center of the pump pulse ($ {x^\prime} = 0 $), the generated terahertz pulse possesses its largest bandwidth. However, toward the sides of the OP beam, the bandwidth of the terahertz pulse reduces. In other words, because of the second-order angular dispersion ($ {(\omega - {\omega _0})^2} $ related term), the OP experiences a temporal chirp and, thus, the pulse duration varies with respect to $ {x^\prime} $. The effective instantaneous bandwidth of the OP reduces toward the sides of the beam, leading to a narrower terahertz spectrum (multicycle pulses).

3. COMPARISON OF THE ${\rm 1D} + 1$, ${\rm 2D} + 1$, AND ${\rm 3D} + 1$ SIMULATIONS

Owing to the trapezoid geometry of the LN crystal, the interaction length varies with respect to $ {x^\prime} $. Since the generated terahertz pulses act back onto the OP, cascading occurs. This causes a spatially dependent OP spectrum, which further enhances spatial inhomogeneities of the generated terahertz pulses. Furthermore, at the desired terahertz frequency range ($ \lt {4}\;{\rm THz}$), the material absorption ($ \alpha $) increases with respect to frequency. This favors lower terahertz frequencies toward the base of the LN crystal due to a longer interaction length. The aforementioned effects contribute to the spatial inhomogeneity of the terahertz pulses, in addition to the higher-order angular dispersion caused by the grating (Section 2). The combined effect can only be investigated by a robust numerical model.

Our numerical model solves the coupled wave equations with slowly varying amplitude approximation in the terahertz coordinates (${x} - {y} - {z}$). By setting the electric field of the OP $ {\cal F}[{E_{{\rm op}}}(t,x,y,z)] = {E_{{\rm op}}}(\omega ,x,y,z) = E(\omega ,x,y,z){e^{ - i[{k_{z0}}(\omega )z + {k_{x0}}(\omega )x]}} $ and the electric field of the terahertz $ {\cal F}[{E_{{\rm THz}}}(t,x,y,z)] = {E_{{\rm THz}}}(\Omega ,x,y,z) = E(\Omega ,x,y,z){e^{ - i{k_0}(\Omega )z}} $, respectively, one can get Eqs. (5) and (6). The operator $ {\cal F} $ represents the Fourier transform. Equations (5) and (6) are as follows:

In Eqs. (5) and (6), $ \Delta {k_x} = {k_{x0}}(\omega ) - {k_{x0}}(\omega + \Omega ) $, $ \Delta {k_z} = {k_{z0}}$ $ (\omega )- {k_{z0}}(\omega + \Omega ) + {k_0}(\Omega ) $. The $ {\chi ^{(2)}} $-related terms are responsible for the second-order nonlinear effects, i.e., the terahertz generation and back conversion processes. In Eq. (6), the third-order nonlinear effects, including self-phase-modulation, self-steepening, and the stimulated Raman effect, are represented by the term $ {n_2}(\tau ) = {\cal F}[{n_2}(\omega - {\omega _0})] $ [22], where $ {n_2} $ is the nonlinear refractive index. The phonon resonances at terahertz frequencies [23] are implemented by considering the stimulated Raman effect at the optical frequency region together with the frequency-dependent refractive index in the terahertz frequency region. In the simulation, the refractive index and terahertz absorption are frequency dependent. The parameters used are listed in Table 1. The peak fluence of the OP at the input LN crystal surface is chosen to be right beneath the estimated damage threshold $ 70.7\;{\rm mJ}/{{\rm cm}^2} $ based on our previous studies [24].

Table 1. Simulation Parameters

View Table

 

Fig. 2. Comparison of the results obtained from ${\rm 1D} + 1$, ${\rm 2D} + 1$, and ${\rm 3D} + 1$ simulations. (a), (b), and (c) are the output OP spectra, the output terahertz spectra, and the efficiencies, respectively. The solid curves correspond to the results generated by OP with $ \sigma _x^\prime = 0.44\;{\rm mm} $. The dashed curves in (c) represent the results generated by OP with $ \sigma _x^\prime = 1.32\;{\rm mm} $. (a) and (b) are plotted at the location of maximum efficiency [marked by the vertical gray lines in (c)].

Download Full Size | PDF

By comparing different models, one can see that the ${\rm 1D} + 1$ calculation neglects diffraction effects and, more importantly, the spatial walk-off between the terahertz and OP beams (the operator $ 2i{k_{x0}}\frac{\partial }{{\partial x}} $). The neglected terms are marked in Eqs. (5) and (6) by “neglected by ${\rm 1D} + 1$.” In Fig. 2, the solid curves correspond to the results generated by OP with $ \sigma _x^\prime = 0.44\;{\rm mm} $, and the dashed curves correspond to the results generated by OP with $ \sigma _x^\prime = 1.32\;{\rm mm} $. Figure 2 suggests that, by neglecting any spatial walk-off effect, the ${\rm 1D} + 1$ model overestimates the OP spectral broadening leading to a more pronounced stimulated Raman effect. Consequently, higher terahertz frequencies can be generated within the OP bandwidth via difference-frequency generation. This explains why the ${\rm 1D} + 1$ model predicts higher terahertz frequency content as compared to the ${\rm 3D} + 1$ model shown in Fig. 2(b). The ${\rm 2D} + 1$ model, though neglecting diffraction in the dimension $ y $ (labeled by “neglected by ${\rm 2D} + 1$”), captures the broadening of the OP spectrum better than the ${\rm 1D} + 1$ case. Accordingly, the terahertz spectrum is predicted in very good agreement with the ${\rm 3D} + 1$ model.

Figure 2(c) shows the computed conversion efficiency along the terahertz propagation direction ($ z $). Two dashed curves correspond to the results generated by OP with $ \sigma _x^\prime = 1.32\;{\rm mm} $. The effective length is strongly dependent on the OP beam size. Furthermore, the calculated efficiency agrees well with the existing experiment [25]. Additionally, the ${\rm 1D} + 1$ model drastically underestimates the optimal interaction length, which also leads to a significant overestimation of the conversion efficiency. The interaction length is captured better by the ${\rm 2D} + 1$ model, while the conversion efficiency is still overestimated by about 25% compared to the value obtained by a ${\rm 3D} + 1$ calculation, since the reduction of the OP fluence along the $ y $ dimension is not accounted for in this model and only properly captured in the ${\rm 3D} + 1$ case.

It can be concluded that in order to capture the key characteristics of the OP and terahertz spectra, at least a ${\rm 2D} + 1$ model should be used while for a proper prediction of the conversion efficiency both the ${\rm 1D} + 1$ and ${\rm 2D} + 1$ models overestimate the efficiency significantly. Here, a ${\rm 3D} + 1$ model is recommended.

4. SPATIAL DEPENDENCE OF THE TERAHERTZ ELECTRIC FIELD

Without loss of generality, the nonlinear interaction between the OP and the LN crystal is numerically implemented in the $ x - y - z $ coordinate frame where $ x = 0 $ represents the apex location of the LN crystal. Note that the OP beam size in the $ x - y - z $ coordinates $ {\sigma _x} = {\sigma _{{x^\prime}}}/\cos (\gamma ) $ is due to the projection onto the plane of the tilted pulse front. The simulations suggest that within an OP beam size range $ {\sigma _y} = [0.5,4.5] \;{\rm mm} $ (not shown), diffraction has a negligible effect on the terahertz generation process, and the terahertz beam size scales as $ {\sigma _y}/\sqrt 2 $. This agrees well with the analytic result in Eq. (4). In the following simulations, $ {\sigma _y} $ is chosen to be 3.5 mm.

In Fig. 3, the maximum terahertz generation efficiency is plotted against the OP beam size $ {\sigma _{{x^\prime}}} $ for two different OP peak fluences. For this calculation, the ${\rm 2D} + 1$ model is chosen due to the high computational cost of the ${\rm 3D} + 1$ model.

 

Fig. 3. With the input pump fluence $ 70.7\;{\rm mJ}/{{\rm cm}^2} $ (blue dots) and $ 35.3\;{\rm mJ}/{{\rm cm}^2} $ (orange dots), the maximum terahertz generation efficiencies versus the OP beam size $ {\sigma _{{x^\prime}}} $, calculated by the 2D model, are presented. The black circles indicate three beam sizes chosen as examples in the following ${\rm 3D} + 1$ calculations.

Download Full Size | PDF

 

Fig. 4. Spatial dependence of the generated terahertz beams along $ x $ and $ y $ dimensions. (a–c) and (d–f) represent the OP and the terahertz beam profiles at the output surface of the LN crystal, respectively. (g–i) and (j–l) represent the OP and terahertz fluence, respectively, at a given position $ y = 0 $ (black curve) and $ y = {\sigma _y}/\sqrt 2 = 2.47 $ (red curve). The OP beam sizes at the input LN-crystal surface are $ \sigma _x^\prime = 0.44\;{\rm mm} $, 0.88 mm, and 1.32 mm in the $ {x^\prime} - {y^\prime} - {z^\prime} $ frame, respectively. The center position of the OP beam is marked by the dashed line. The OP beam size in the $ y $ dimension is chosen to be $ {\sigma _y} = 3.5\;{\rm mm} $. The apex of the LN crystal is located at $ x = 0 $.

Download Full Size | PDF

Owing to the nature of the noncollinear phase-matching condition, the terahertz generation process requires different sections of the beam along the $ x $ dimension to add up coherently in the emission direction $ z $. In contrast with the OP beam size in the $ y $ dimension, a small beam size along the $ x $ dimension cannot produce high generation efficiencies due to the walk-off between the OP and the terahertz beam. On the other hand, if the beam size is too large, the terahertz radiation generated by the side of the OP at the farther side from the LN-crystal apex suffers from more absorption compared to the part closer to the apex. Thus, the generation efficiency shows a maximum as a function of OP beam size in the $ {x^\prime} $ direction. Additionally, since lower pump fluence leads to longer interaction length, the optimal pump beam size increases (see the orange dots in Fig. 3).

As in experiments different OP beam sizes may be required in order to optimize the use of the available pump energy and limited crystal aperture; we select three pump sizes ($ {\sigma _{{x^\prime}}} = 0.44\;{\rm mm}$, 0.88 mm, and 1.32 mm, marked by black circles in Fig. 3) for studying the spatiotemporal properties of the terahertz field using the ${\rm 3D} + 1$ model.

In order to compare the terahertz fields generated by different spatial positions of the OP beam, the center of the OP (highest peak fluence $ y = 0 $) and the side region of the OP (low fluence, $ y = {\sigma _y}/\sqrt 2 = 2.47\;{\rm mm} $) within the optical beam are chosen. Figures 4(a)–4(f) show the OP and the terahertz beam profiles at the output LN-crystal surface with different input OP beam sizes. For the beam sizes $ \sigma _x^\prime = 0.44 \;{\rm mm} $, 0.88 mm, and 1.32 mm, the conversion efficiencies are 0.57%, 0.54%, and 0.46%, respectively. Figures 4(g)–4(i) present the fluence of the OP at $ y = 0 $ and $ y = {\sigma _y}/\sqrt 2 $ with the center of the OP marked by the dashed lines. The fluence distribution indicates an energy shift towards the base of the LN crystal. This shift is due to the noncollinear nature of the phase-matching mechanism, which causes the new optical frequencies of the OP, generated via the cascading effect, to propagate toward the base of the crystal. Surprisingly, the size of the OP does not strongly influence the size of the generated terahertz beams. Toward the base of the LN, the interaction length increases, and the terahertz absorption and the nonlinear effect become more pronounced. The terahertz fields generated close to the base is absorbed before they reach the output LN surface. Thus, for large OP beam sizes, the OP farther from the LN apex is wasted. Figures 4(j)–4(l) indicate that at a given OP beam size, a higher pump fluence leads to a smaller terahertz beam size compared to the case using a lower pump fluence. This finding agrees with the experimental results of Lombosi et al. [27]. The terahertz beams are found to be symmetric along the $ y $ dimension.

 

Fig. 5. Spatial dependence of the generated terahertz spectra and temporal profiles along the $ x $ dimension. (a–c) and (d–f) are the terahertz electric field and the corresponding terahertz spectra with respect to $ x $ at $ y = 0 $. (g–i) and (j–l) show the terahertz electric field and the corresponding terahertz spectra with respect to $ x $ at $ y = {\sigma _y}/\sqrt 2 $.

Download Full Size | PDF

Figures 5(a)–5(c) and 5(g)–5(i) show the spatially dependent electric fields at $ y = {\sigma _y}/\sqrt 2 $ and $ y = 0 $, respectively, with the corresponding terahertz spectra shown in Figs. 5(d)–5(f) and 5(j)–5(l). It can be seen that the few-cycle terahertz electric fields are only generated at the vicinity of the apex of the LN crystal. The inhomogeneity along the $ x $ dimension can lead to large terahertz diffraction [27,28]. This effect is seen to increase with larger pump beam sizes, as a larger fraction of the terahertz is generated farther from the apex and the terahertz electric fields deviate from a single-cycle waveform.

To quantify the deviation from a single-cycle waveform, the root-mean-square pulse duration $ \Delta t $ is chosen to evaluate the electric field distribution [see Eq. (6)]. In Eq. (6), $ \delta t = |t(x,y) - t{(x,y)_p}| $, $ t{(x,y)_p} $ is the time coordinate of the peak of the electric field at position $ (x,y) $ and $ I $ represents the intensity. The reference $ \Delta t({x_p},{y_p}) $ is chosen at the position ($ {x_p},{y_p} $) where the terahertz peak fluence is located. Equation (7) is as follows:

The resulting map of $ \Delta t(x,y) $ can be used to determine the portion of the terahertz beam, where the electric field deviates significantly from a single-cycle format (see Fig. 6). Figure 6 shows the terahertz beam generated by an OP with $ {\sigma _x} = 1.32\;{\rm mm} $, where the part with $ \Delta t(x,y) \gt 2\Delta t({x_p},{y_p}) $ is shaded gray to indicate the non-single-cycle content.

 

Fig. 6. Example shown is for $ \sigma _x^\prime = 1.32 \;{\rm mm} $, where the non-single-cycle region, $ \Delta t(x,y) \gt 2\Delta t({x_p},{y_p}) $, is indicated by the gray region. The terahertz beam under the gray region contains up to 25% of the total terahertz energy.

Download Full Size | PDF

The gray region in Fig. 6 shows that high-quality few-cycle pulses are only generated close to the apex of the crystal ($ x = 0 $). The terahertz beam profiles are symmetric along the $ y $ direction. We find that even though higher pump fluence ($ y = 0 $) leads to a higher terahertz energy the portion of non-single-cycle content in the beam increases. It can be seen from Figs. 5 and 6 that the side of the OP beam (for example, $ y = {\sigma _y}/\sqrt 2 \;{\rm mm} $) possesses lower fluence, which leads to less OP spectral broadening in the terahertz generation process. Such a less broadened OP spectrum leads to a relatively homogeneous terahertz distribution along the $ x $ dimension. For $ {\sigma _{{x^\prime}}} = 0.44\;{\rm mm} $, 0.88 mm, and 1.32 mm the $ \Delta t(x,y) \gt 2\Delta t({x_p},{y_p}) $ region (gray region) takes up 4%, 20%, and 25% of the total terahertz energy, respectively. It can be seen that, with the increase of the OP beam size, the single-cycle region of the generated terahertz pulses reduces significantly. Additionally, due to the geometry of the nonlinear crystal, it is inevitable that the generated terahertz possesses a spatial inhomogeneity. In order to resolve this problem, a setup that combines a conventional tilted-pulse-front setup and a transmission stair-step echelon is promising to generate spatially homogeneous terahertz pulses [29].

5. CONCLUSION

We have studied the spatial dependence of the terahertz electric field properties generated in a tilted-pulse-front setup. By comparing the ${\rm 1D} + 1$, ${\rm 2D} + 1$, and ${\rm 3D} + 1$ numerical models we found that the ${\rm 1D} + 1$ calculation cannot capture certain key features of the terahertz generation process. Additionally, the ${\rm 3D} + 1$ calculation shows that within the OP beam size range $ {\sigma _y} = [0.5,4.5] \;{\rm mm} $, diffraction in the dimension perpendicular to the pulse-front-tilt plane does not have much influence on the terahertz generation process. For those cases, the 2D calculation (${x} - {z}$ coordinates) is a good approximation.

Perpendicular to the pulse front tilt plane, in the $ y $ direction, the terahertz beam waist relates to the OP waist with $ {\sigma _y}/\sqrt 2 $, which is as expected for the second-order nonlinear process. However, in the $ x $ dimension, the OP beam size does not have a large impact on the generated terahertz beam size. The terahertz pulses are generated close to the apex of the crystal and the single-cycle region is located close to the vicinity of the crystal apex. Thereby, large OP beam sizes lead to a reduced percentage of single-cycle content of the generated terahertz pulses. Attention must be paid to the fact that the terahertz generated farther from the apex of the crystal, which can possess a significant fraction of the total generated terahertz energy, suffers from poor electric field quality and temporal chirp.

In order to maximize the single-cycle content, the OP beam size along the $ x $ dimension should be kept reasonably small, while the size along the $ y $ dimension can be used to enlarge the pump beam if necessary. Another option is to reduce the OP input fluence. These findings are of particular relevance to carrier-envelope-phase sensitive terahertz applications and strong filed terahertz physics.