1. Introduction

The last decades have seen a surge in research studies on generation and applications of terahertz radiation, which typically refers to electromagnetic waves in the spectral range from 0.1 to 3 THz. Among the wide range of applications, spectroscopy [1–3], spin dynamics control [4,5] and linear electron acceleration highly benefit from high power terahertz sources [6–8]. In particular, applications like particle accelerations place steep requirements of a few millijoules of single/multi-cycle terahertz pulse energy to enable bunch manipulation in the relativistic regime. Such performance is largely contingent on increasing the terahertz generation efficiency towards the percent level and beyond.

Different approaches have addressed the generation of few-cycle terahertz pulses, including four wave mixing in ambient air [9], plasma driven effects in air [10], terahertz emission from photoconductive switches [11, 12], tilted-pulse-fronts [13] and echelons [14–16]. In addition, there are techniques available for the generation of multi-cycle terahertz radiation. Quantum cascade lasers are attractive due to their compactness and production simplicity in its use in the spectral range > 1 THz despite their limited tunability [17–19]. Molecular gas lasers can provide high terahertz energies but are less tunable in terms of the generated terahertz frequency [20]. Free electron based terahertz sources are ideal in providing high energy multi-cycle terahertz pulses but are less accessible and are even more challenging for a compact implementation [21,22]. Laser based multi-cycle terahertz generation can leverage on developments in solid-state laser technology to enable compact coherent terahertz sources with high conversion efficiency at low terahertz frequency (<1 THz) [23–25]. This study revolves around the generation of high energy terahertz pulses in the multi-cycle (narrowband) regime.

Laser-driven terahertz generation utilizes the second order nonlinearity of the nonlinear material to perform difference frequency generation. GaAs, CdTe, InP [26–28] HMQ-TMS [29], LiNbO₃ [30], GaP, ZnTe [31,32] and DAST [33,34] are nonlinear materials that have received significant attention. Multi-cycle terahertz generation in periodically poled Lithium Niobate (PPLN) was initially demonstrated by Lee et al. [30]. Recent research efforts focused on terahertz generation in PPLNs with various pump pulse formats such as chirp and delay, pulse-train and optical parametric amplification [23, 35–37]. To our knowledge, the highest experimentally achieved conversion efficiency for multi-cycle terahertz generation in PPLN is 0.13% [38].

Increasing the optical to terahertz generation/conversion efficiency (referred to as efficiency throughout this paper) in laser-based schemes is often hindered by several challenges. The first limitation is introduced by the Manley-Rowe relation [39] obtained after the assumption of single photon conversion. This limitation can be overcome by a cascaded second order process (repeated energy down-conversion of pump photons to terahertz photons) [36,40]. The damage threshold of the nonlinear crystal imposes another challenge. Large pump electric field strength induces optical breakdown and damage to the material, which also limits the achievable terahertz efficiency [41]. Another challenge is introduced by the terahertz absorption of the material. This challenge can be overcome by cryogenic cooling [42,43], since the temperature decrease impedes thermal phonon excitation and, therefore, considerably reduces the terahertz absorption. Additionally, it was found by Lee et al. that cooling of the nonlinear crystal increases the terahertz generation efficiency [44].

Our previous work based on a 1D-propagation model [35], partially addressed these issues and showed the possibility of obtaining high efficiencies of several percent for optimal pump pulse formats. However, in practice, complex pump pulse formats are hard to engineer. Therefore, further innovations are necessary to achieve the desired percent level efficiencies. Here, we aim to alleviate these limitations by recycling the optical pump and separating the generated terahertz radiation in a staged approach, before significant terahertz absorption occurs. In this multi-stage architecture, consecutive stages utilize the same optical pump pulse. The terahertz pulse generated at each PPLN stage is then coupled out, circumventing the excessive terahertz absorption. A similar concept was experimentally tested by Chen et al. [45] with a two-stage system for single-cycle terahertz generation. However, in their system, the terahertz pulse generated in the first stage is also coupled into the second stage together with the pump pulse [45].

In this paper, section 2 describes the coupled wave equations governing terahertz generation and the utilized numerical approach based on a finite difference method in cylindrical coordinates. The effects of dispersion and mitigation mechanisms are discussed. Furthermore, the critical aspect of isolating the co-propagating terahertz and optical beams at each stage is resolved with a suitably designed quartz coupler. Section 3 analyzes the dependence of terahertz generation on pump pulse duration, effective length and beam size. Also, the influence of phase front distortions and spatial intensity fluctuations of the pump pulse on terahertz generation are studied. Detailed information about the terahertz spatial profile and terahertz beam combination is presented.

2. Multi-Stage terahertz generation: concept and physics

The schematic illustration of the multi-stage system is presented in Fig. 1. A set of PPLN crystals are placed in a serial arrangement, through which an optical pump pulse is consecutively recycled. At each stage, the terahertz beam is out-coupled by a designed quartz coupler (see section 2.2.2) to avoid excessive terahertz absorption. The out-coupling and separation of the optical and terahertz beams may also be achieved by a combination of anti-reflection coatings for optical and terahertz beams at the exit of the crystal, and a dichroic mirror. In this article, we focus on the quartz coupler due to its simplicity in manufacturing and the low loss for the terahertz beam.

 

Fig. 1 Schematic illustration of the multi-stage system. Detailed information about the quartz coupler (QC) can be found in Fig. 7. In the sketch, the angular deflection of the optical pump caused by the refraction at the output of each stage is not delineated for illustration purposes. The dark arrows on the PPLN crystals represent the poling of the nonlinear material.

Download Full Size | PDF

2.1. Theory

Terahertz radiation is generated in PPLN through difference frequency generation (DFG) processes. The collinear propagation of terahertz and pump beams suggest considering the cylindrical symmetry in solving the governing equations for terahertz generation. For the quartz coupler, however, this cylindrical symmetry is violated. We solve for the propagation of terahertz beam in quartz using an angular spectrum method [46]. For this purpose, the terahertz beam profile is reconstructed in 3-D Cartesian coordinates from the 2-D cylindrical coordinates. It is assumed that the quartz coupler has no influence on the pump pulse (detail in section 3.3). To precisely simulate the evolution of the terahertz and optical pulses, all dominant effects such as DFG, self-phase modulation (SPM), self-focusing, beam diffraction, dispersion and material absorption of the terahertz radiation, are included. In our studies, stimulated Raman scattering is neglected because of its relatively small influence [47], since the desired terahertz frequency range (< 1 THz) is far away from the phonon resonance frequency of Lithium Niobate (LiNbO₃), occurring around 7.5 THz [48].

2.1.1. Coupled wave equations in cylindrical coordinate

We choose to perform calculations on a slice along the radial axis in cylindrical coordinates as shown in Fig. 2. For high energy terahertz generation, the usage of large crystals is mandatory due to the limited damage threshold of the nonlinear materials. Using the state-of-the-art fabrication technology, PPLN crystals as large as 1 cm × 1 cm transverse dimensions and few centimeters in length are realizable [49]. The propagation of the pump and terahertz beams in bulk PPLN crystals resemble the propagation of linearly polarized beams in a square waveguide. However, the waveguide dimensions are far larger than the wavelength of both pump and terahertz in the crystal. As a result, the effect of the waveguide boundaries become negligible.

 

Fig. 2 Schematic illustration of the simulated geometry: the dark thick arrows represent the polarization direction of both pump and terahertz beams. The polarization direction is aligned with the extraordinary optical axis of PPLN. The origin of the cylindrical coordinate is at the center of the beam. r and z represent the transverse and propagation directions, respectively.

Download Full Size | PDF

Consequently, the entire problem reduces to solving the scalar coupled wave equations with the ansatz E(ω, r, z) = A(ω, r, z)e^−ik(ω)z, considering the slowly varying amplitude approximation. The governing equations of the field amplitude of pump and terahertz beams read as:

In Eqs. (1) and (2), A_op and A_THz stand for the electric field amplitude of the optical pump and the terahertz beam respectively; ℱ denotes the Fourier transform operator; ω and Ω are the optical and terahertz angular frequencies respectively; k represents the wavenumber; α denotes terahertz absorption; n₂ represents the nonlinear refractive index and n represents the refractive index. The effective nonlinear coefficient is written as ${\chi }_{\text{eff}}^{(2)}={\chi }_{333}^{(2)}(z)$ in the last terms on the right hand side of Eqs. (1) and (2). The term e^iΔkz = e^{i(k(Ω)+k(ω)−k(ω+Ω))z} evaluates phase mismatch between the propagation of the terahertz beam and the optical pump. The term e^iΔkz together with the term in the first order Fourier expansion series of ${\chi }_{333}^{(2)}(z)$, leads to the phase-matching condition [35].

On the right hand side of Eqs. (1) and (2), the first terms correspond to the transverse spatial dependence of the propagation, describing beam divergence due to diffraction. In Eq. (1), the second term describes the impact of SPM on the optical beam which is independent of the phase-matching condition. As one can see, we consider a frequency dependent SPM effect which implicitly includes the effect of self-steepening (see Appendix 5.2). Additionally, by considering radially dependent SPM, self-focusing is considered. The third term captures two important effects. One is the DFG term $A_{\text{op}}(\omega +\mathrm{\Omega })A_{\text{THz}}^{*}(\mathrm{\Omega })$ term which describes the frequency down-shift of the pump, and the other is the A_op(ω − Ω)A_THz(Ω) term which describes the back conversion from terahertz to the pump beam (sum-frequency generation) [50]. The second and third terms on the right hand side of Eq. (2) evaluate the effects of terahertz absorption and DFG processes respectively.

Our 2-D calculations in cylindrical coordinates allow for the possibility of including the spatial variation of the beam, diffraction and self-focusing. This can not be addressed by a 1-D calculation.

2.1.2. Numerical method and simulation parameters

Various finite difference methods can be used for solving Eqs. (1) and (2). In this study, a split step Fourier method is used to compute the integral terms corresponding to the nonlinear polarization terms containing n₂ and ${\chi }_{\text{eff}}^{(2)}$. The second order derivative in the r direction is expressed by a conventional 3-point finite difference method as shown by Eq. (15) in Appendix (5.3). Ultimately, the values of A_op and A_THz are updated throughout the propagation distance z using the low-storage Runge-Kutta method [51]. The entire method leads to an explicit update algorithm along z rather than in time t. This is only valid when there is no reflection in the system (i.e only forward propagation along z is involved). The 3-point finite difference method takes advantage of the tridiagonal matrix of the discretization along the r direction, thus reducing computational cost. A combination of split step, explicit finite difference, and low-storage Runge-Kutta methods highly enhance the computational performance. Using other finite difference methods such as the Crank-Nicolson method is also possible. However, this implicit method necessitates the inversion of a finite difference matrix for each frequency and increases computational cost. Additionally, an angular spectrum method based on Hankel transformation [52] can be used but the approach suffers from a lack of the development of the fast Fourier transform (FFT) in Bessel series [53] and the boundary conditions are less straightforward to be implemented. Our calculation is performed in C++ with MPI and OpenMP for code parallelization. Details of the numerical method can be found in Appendix (5.3).

In this work, we base all simulations for terahertz generation centered at 0.3 THz due to its amenability for electron acceleration. In principle, the architecture of the system can be extended to various terahertz frequencies. However, it would require redesigning the PPLN and the coupler. Within the same system, the frequency tunability is limited by the phase-matching bandwidth which varies within a range of few GHz. Terahertz generation with pulse trains is effective in mitigating walk-off and yields high conversion efficiencies [35]. In the remainder of this paper, the simulations are based on pump pulses consisting of two narrow spectral lines (referred to as two-lines throughout this paper) separated by the desired terahertz frequency (Ω₀ = 2π×0.3 THz). This forms a pulse train (beat note) with each sub-pulse separated by a time 2π/Ω₀, i.e. the terahertz electric field period. Several approaches could be used to generate the two-lines spectrum. By angularly spreading a broad-band pulse, the desired spectral lines can be obtained with a spatial filter. Optical parametric amplification with a narrow bandwidth seed is another possibility. The signal and idler together form the two-lines spectrum. The simulation parameters are tabulated in Table 1.

Table 1. Simulation Parameters with Two-lines Input

View Table | View all tables in this article

The refractive index and absorption (α) of terahertz radiation at 75 K and 0.3 THz are obtained from linear interpolation of the experimental data between temperatures 50 K and 100 K from [54]. The material absorption at optical wavelengths, i.e. the pump absorption, is neglected and the refractive index is fitted to a Sellmeier equation [55]. We assume a super-Gaussian spatial distribution (M = 5) and a Gaussian temporal profile for the optical pump. This is also the case for high power lasers. Additionally, a super-Gaussian spatial distribution ensures that the generated terahertz radiation is homogeneous in the transverse r dimension and also minimizes the damage due to peak fluence of the pump compared to a Gaussian spatial profile. The input pump pulse energy, which follows from spatial and temporal integrations of the electric field defined in Table. 1, is given by

2.2. Terahertz efficiency enhancement with dispersion compensation

We suggest compensating the dispersion accumulated in the material of the pump pulse generated at the end of each stage before recycling it to the subsequent PPLN stage in order to enhance the conversion efficiency. In other words, before the pump pulse is recycled to stage N, an opposite dispersion is added to the pump pulse at the end of stage N-1 to cancel the dispersion obtained by propagation through the stage N-1. For the first and second stages, dispersion compensation can be realized by a prism or grating (i.e second and third order dispersion). For many stages, the cascading process broadens the spectrum drastically. Dispersion compensation for a large bandwidth beam can be challenging in practice since higher order dispersion needs to be taken into consideration. However, it still could be realized using double chirped mirrors [59]. Detailed analysis and design of such a system is out of the scope of this paper.

The simulation results comparing the terahertz spectra at the end of each stage with and without dispersion compensation are shown in Figs. 3(a)–3(b). The resulting conversion efficiencies are shown in Figs. 3(c)–3(d).

 

Fig. 3 Comparison of a 4 stage system with input pump pulse parameters σ = 5 mm, τ_FWHM = 150 ps. The effective length of each stage is chosen to be the distance where efficiency saturates. (a,c,e) depict the terahertz spectra generated after each stage, the conversion efficiency, and the terahertz electric field by recycling the pump pulse directly without dispersion compensation (No Comp.). (b,d,f) show the terahertz spectra, the conversion efficiency, and the terahertz electric field with pump pulse dispersion compensation (Comp.) after each stage.

Download Full Size | PDF

One can see that the terahertz conversion efficiency reduces subsequently after each stage without dispersion compensation. Since the terahertz generation broadens the pump spectrum via the cascading process, part of the broadened pump spectrum can not be phase matched for further cascading. In other words, different frequency elements in the pump spectrum pick up different phases due to dispersion in the PPLN via propagation. Consequently, the terahertz radiation generated by DFG from different spectral ranges of the pump pulse carries different phases, leading to partial destructive interference of the total generated terahertz radiation. Thus, the efficiency is degraded due to the dispersion. It can be seen in Figs. 3(b)–3(d) that the efficiency can be significantly enhanced by dispersion compensation. Besides, Figs. 3(e)–3(f) show that dispersion compensation can enhance the terahertz temporal profile.

2.2.1. Spectral dynamics of optical pump and terahertz pulses

Figure 4 depicts the temporal profile of the pump pulse with the short time Fourier transform (STFT) (i.e instantaneous spectrum) of the corresponding time range. Figure 4(a) shows the input pump pulse which consists of two spectral lines. Figure 4(b) shows that the broadening of the pump pulse spectrum causes a reduction of the pulse duration of each sub-pulse and an increase in peak intensity. Besides, the STFT in Fig. 4(b) indicates that the instantaneous spectrum in each sub-pulse forms a ’U’ shape due to the cascading effect, leading to a drastic spectral variation in time. The maximum cascading (spectral down-shift) occurs where the highest peak intensity is present. The STFT in Fig. 4(c) suggests that the dispersion compensation alters the spectral distribution with respect to time. If a transform limited pump pulse, which has a uniform spectrum distributed over time, could be obtained after each stage, the terahertz efficiency could be greatly enhanced. However, the realization of such a process is challenging, since perfect compression of the ’U’ shaped spectrum induced by SPM, the second order nonlinear effect, and dispersion, is not straightforward.

 

Fig. 4 (a) shows the input pump pulse (b) shows the output pump pulse after one 2.5 cm PPLN stage. (c) shows the output pump pulse after dispersion compensation by adding the opposite second order dispersion of Lithium Niobate at the pump center frequency (GDD= 6300 fs²) to (b). The contour plot in the second and third rows are the short time Fourier transforms of the corresponding selected temporal range (indicated by the dotted red box).

Download Full Size | PDF

Figure 5 compares the spectral density (spatially dependent spectra in radial dimension r) of the terahertz generated at each stage for two cases with and without dispersion compensation. The terahertz beams generated with dispersion compensation follow a Gaussian-like spectrum, whereas those generated by direct pump pulse recycling have distorted spectral shapes, particularly in the last stage. This is due to the fact that the phase of the pump perturbed by dispersion influences the phase-matching condition of terahertz generation, and thus leads to a broader terahertz spectrum. Figure. 6 shows the comparison of the spectral density of the output optical pump at each stage without and with dispersion compensation. It can be seen that dispersion compensation boosts the cascading process and leads to greater pump spectral broadening. The results depicted in Figs. 3–6 confirm that the terahertz conversion efficiency can be largely enhanced by pump pulse recycling with dispersion compensation.

 

Fig. 5 (a–d) show the terahertz spectral density generated by 4 consecutive stages respectively. The first and second rows show the terahertz spectra generated without and with dispersion compensation, respectively.

Download Full Size | PDF

 

Fig. 6 (a–d) represent pump pulse spectral density after each corresponding stage. The first and second rows show the pump pulse spectra at the output of each stage generated without and with dispersion compensation, respectively.

Download Full Size | PDF

2.2.2. Quartz output coupler for beam separation

Figure. 7(a) shows one possible arrangement of the multi-stage setup accounting for the subtlety of terahertz and optical beam separation. The output pump beam from the PPLN propagates at a certain angle with respect to the input pump pulse direction due to the refraction from the quartz coupler (QC). However, in the sketch, we do not depict this for simplicity of illustration. Incidence at Brewster angle plays a crucial role in this scheme in order to couple out and separate both pump and terahertz beams with high transmission efficiency from the PPLN to air (n = 1). Out-coupling of the terahertz beam from Lithium Niobate to air through a silicon coupler was first demonstrated experimentally and theoretically by Kawase et al. [60]. In this article, crystalline quartz is chosen as the material of the output coupler owing to its low absorption for terahertz radiation at 0.3 THz (0.01/cm) [61, 62], low nonlinearity compared to silicon or Lithium Niobate [63], and large bandgap (12.5 eV) [64] which prevents the ionization effect due to the pump pulse. Additionally, the dispersion of quartz at both pump and terahertz frequencies is negligible for few centimeters of propagation [62,64]. We neglect the small influence of the quartz coupler on the optical pump and fully focus on the out-coupling of the terahertz radiation.

 

Fig. 7 (a) Schematic illustration of the output coupler and the angle definitions. In (b) and (c), comparison of the terahertz transmission from PPLN to air for ’Air’ and ’QC’ cases are shown. Additionally, the transmission of the generated flat-top terahertz beam with different beam sizes σ_t = 1 mm (dashed lines) and σ_t = 5 mm (solid lines) are analyzed. (b) shows terahertz transmission versus the refractive index of the PPLN at the given Brewster angles (BA) at S₁ surface. The BA is calculated with the PPLN refractive index at terahertz frequency marked by the black dash-dotted line. (c) shows terahertz transmission versus the incident angle at surface S₁ with given refractive index where the black dash-dotted lines mark the BA.

Download Full Size | PDF

In Fig. 7(a), the inner angles of the bulk quartz coupler are represented by α₁, α₂, and α₃. This ensures that the terahertz beams experience Brewster incidence angles π − α₁ and γ₂ at S₁ and S₂ surfaces, respectively. β₁ and β₂ are the output angles of the pump beam at the S₁ surface and the pump beam Brewster incidence angle at the S₃ surface, respectively. The refraction loss of the pump is negligible. S₃ can be chosen to be parallel to surface S₁ to further reduce angular deviation. Additionally, L (the length of the top edge of the quartz coupler) is the minimum distance that ensures sufficient spatial separation of pump and terahertz beams. Details of the equations used for calculating the involved angles can be found in Table 2, where σ = 5 mm, n_LN(ω₀) = 2.1, n_Q(ω₀) = 1.4, n_LN(Ω₀) = 5, n_Q(Ω₀) = 2. The subscripts ’LN’ and ’Q’ correspond to PPLN and Quartz respectively.

Table 2. Quartz Coupler Parameters

View Table | View all tables in this article

Since the refractive index of the PPLN varies with temperature and manufacturing process, an output coupler with a reasonable tolerance is desired. In Figs. 7(b)–7(c), comparison of two out-coupling cases ’Air’ (red curves) and ’QC’ (blue curves) is shown. In the case of ’Air’ where no output coupler is utilized, the PPLN crystal is cut such that the terahertz pulse is incident at Brewster’s angle (tan⁻¹(1/n_LN(Ω₀)) = tan⁻¹(1/5) = 11.3°) from PPLN to air directly. In the case of ’QC’, the terahertz pulse is coupled out from the PPLN to air through the designed quartz coupler shown in Fig. 7(a) with Brewster incidence angle (tan⁻¹(n_Q(Ω₀)/n_LN(Ω₀) = tan⁻¹(2/5) = 21.8°) on S₁ surface.

In both Fig. 7(b) and Fig. 7(c), the transmission window of the ’QC’ case at σ_t = 5 mm (flat interval with terahertz energy transmission T ≈ 1) is due to total internal reflection (TIR) at both the S₁ and S₂ surfaces. This window range reduces for the smaller terahertz beam size σ_t = 1 mm (i.e terahertz beam with larger angular divergence). The terahertz beam components with angular divergence larger than the TIR are filtered out, leading to lower energy transmission. This TIR induced transmission loss is particularly pronounced for the ’Air’ case since the refractive index varies from n_LN(Ω₀) = 5 to n_Air(Ω₀) = 1. This large refractive index difference causes a sharp transmission drop at the incidence angle larger than Brewster angle. This explains why the transmission peak occurs at an angle smaller than Brewster’s angle in the ’Air’ case in Fig. 7(c). With larger beam size σ_t = 5 mm, the transmission peak approaches the Brewster angle due to smaller angular divergence. It can be seen that from PPLN to air, the terahertz energy transmission at σ_t = 1 mm is around 45% and 55% for ’Air’ case and ’QC’ case respectively. The terahertz energy transmission at σ_t = 5 mm is around 63% and 96% for ’Air’ case and ’QC’ case respectively. With the given refractive index at the corresponding Brewster incident angle, a quartz coupler leads to higher transmission. More discussion can be found in section (3.3).

3. Analysis

3.1. Single stage: terahertz efficiency in terms of pump pulse duration

The crucial parameters in designing a complete terahertz source are the pulse duration of the pump and the necessary length of the PPLN crystal. As one can see in Eq. (3), longer pulse durations allow higher input energy, but simultaneously necessitate longer interaction length, which in turn intensifies the effects caused by terahertz absorption. The optimal pump pulse duration is found to be 150 ps with conversion efficiency η = 1.05 % as shown in Fig. 8. In order to define the optimal PPLN length L_eff, we introduce the length parameter L₀ where the efficiency increase is most rapid. This is obtained at d²η(z)/dz²|_z=L0 = 0. We define L_eff as:

 

Fig. 8 Numerical results of the terahertz generation in a single stage: (a) shows the conversion efficiency versus the pump pulse duration. The maximum efficiency η = 1.05 % is obtained at τ_FWHM = 150 ps. The input pump energy is calculated from Eq. (3) with the pump beam size σ = 5 mm. (b) shows the L₀, L_eff versus pump pulse duration. (c) shows the terahertz efficiency versus the propagation distance for τ_FWHM = 150 ps, where the steepest efficiency increase rate (L₀) and effective length (L_eff) are labeled.

Download Full Size | PDF

If we define δ ≈ [n(Ω₀) − n_g(ω₀)] /(cτ_FWHM), Eq. (4) results in the following analytic equations for L_eff and L₀ (see Appendix 5.5 for details):

3.2. Single stage: effects of transverse spatial variations of pump

High power terahertz generation demands high power lasers which usually suffer from phase front distortions and intensity modulations. In the following section, we investigate the influence of the spatial variation of the pump pulse on the terahertz generation process.

3.2.1. Beam size dependence

Here, we aim at exploring the dependence of the terahertz generation efficiency on the pump beam size. For this purpose, we perform simulations for a super-Gaussian (M = 5) pump beam with beam sizes σ = 1 mm, 3 mm and 5 mm, respectively. Figure 9 suggests that the terahertz beam size and spatial profile resemble the corresponding pump properties for large pump beam sizes.

 

Fig. 9 (a) shows the efficiency of terahertz generation versus propagation distance for various beam sizes. (b) shows the spatial profile of the input pump pulse (dashed lines) and output pump pulse (solid lines). (c) shows the output terahertz beam spatial profile. For different beam sizes, the peak fluence is kept constant. The output spatial profiles of both terahertz radiation and the pump are obtained at the end of the interaction length, i.e. z = 25 mm.

Download Full Size | PDF

It can be seen in Fig. 9(a) that the terahertz conversion efficiency is lower when σ = 1 mm compared to σ = 3 mm and 5 mm. The diffraction has a minor influence on the spatial profile of the pump pulse since the Rayleigh range (λ ≈ 1 μm) is about few tens of meters. In a super-Gaussian beam, smaller beam size leads to larger spatial phase gradient at the edge of the spatial profile. Therefore, as observed in Fig. 9(b), the self-focusing effect leads to the peak appearing at the edge of the 1 mm pump beam. The change of the spatial profiles of the terahertz beams in Fig. 9(c) is due to the diffraction of the super-Gaussian beam [65,66]. It can be seen that for a 1 mm terahertz beam, diffraction modifies the super-Gaussian profile to a Gaussian profile after 25 mm generation distance.

3.2.2. Spatial energy and phase front variations

The above results are obtained with perfect super-Gaussian pump spatial profiles and a flat spatial phase front. In order to examine the influence of spatial energy and phase front variations of the pump pulse on terahertz generation with the two-lines pump spectrum, we studied 4 cases: (I) flat phase front with one of the pump spectral lines having a 5% modulation in the electric field strength. (II) flat phase front with both pump spectral lines having a 5% modulation in the electric field strength. (III) flat spatial profile with one of the pump spectral lines having modulation from −π to π in the phase front. (IV) flat spatial profile with both of the pump spectrum lines having modulation from −π to π in phase front.

In the following context, E(ω, r) denotes the electric field of the input pump pulse and E₀(ω, r) is the electric field with a perfect super-Gaussian spatial profile (σ=5 mm) and flat phase front. The interaction length is chosen to be z = 25 mm as in Fig. 9. r_max = 8.5 mm is the maximum distance of the calculated radial dimension.

Case I: one line with cosine electric field amplitude modulation

It can be seen in Fig. 10 that terahertz generation is not degraded by the intensity modulation in one of the spectral lines. Though the transverse intensity modulation can also induce spatial phase variations due to self-focusing, with our pump pulse parameters, the induced spatial phase variation is negligible.

 

Fig. 10 (a) and (b) show the spectral density (first row) and flat phase front (second row) of the input pump spectral lines. (c) shows the spatial profile of the input and output pump pulse. (d) shows the terahertz spatial profile.

Download Full Size | PDF

Case II: both lines with cosine electric field amplitude modulation

Fig. 11 suggests that the terahertz generation is not degraded when both lines are subjected to the same intensity modulation, which is expected from the Case I study.

 

Fig. 11 (a) and (b) show the spectral density (first row) and flat phase front (second row) of the input pump spectral lines. (c) shows the spatial profile of the input and output pump pulse. (d) shows the terahertz spatial profile.

Download Full Size | PDF

Case III: both lines with cosine phase modulation

Fig. 12 suggests that if both of the pump lines are identically modulated in phase, the terahertz generation process is not degraded. The phase modulation induces self-focusing which can be observed in Fig. 12(c).

 

Fig. 12 (a) and (b) show the spectral density (first row) and phase front (second row) of the input pump spectral lines. (c) shows the spatial profile of the input and output pump pulse. (d) shows the terahertz spatial profile.

Download Full Size | PDF

Case IV: one line with cosine phase modulation

Figure 13 indicates that a distortion in the phase front, not unusual for high power laser beams, strongly reduces the conversion efficiency. This is due to the nature of the difference frequency generation process, which transfers the phase difference between the two frequency lines of the optical pump to the terahertz beam at a thousand times longer wavelength. In turn, the modulated terahertz phase front suffers from strong diffraction, leading to destructive interference. Consequently, the terahertz conversion efficiency is greatly reduced and a strongly spatially modulated terahertz beam appears as shown in Fig. 13(d). The peaks in Fig. 13(c) are due to self-focusing, making it different from case III. In case III, since the phase fronts of both spectral lines are identical, the terahertz generated inherits a flat phase front. As a result, the cascaded frequencies newly generated in the optical region maintain the same phase front as the initial two-lines spectrum. However, in case IV since DFG creates a terahertz pulse with the phase difference of the pump spectral lines, the new pump frequency lines generated by the cascading process carry different phases. The later generated cascaded lines inherit multiple times of the original phase difference, forming an enormously curved phase front. Thus, strong self-focusing occurs.

 

Fig. 13 In (a,b) the first row shows the spatial profiles of the pump spectral lines and the second row represents the phase front of the corresponding spectral lines. (c) shows the pump input and output spatial profile. (d) shows the generated terahertz spatial profile.

Download Full Size | PDF

3.3. Terahertz spatial profile after quartz coupler

The quartz coupler breaks the cylindrical symmetry of the PPLN. This asymmetric spatial structure can influence the terahertz generation process inside the PPLN crystal. However, with the 25 mm effective length, the influence of the transverse variation due to the asymmetric structure on terahertz generation is negligible. In order to calculate the terahertz spatial profile after the quartz coupler, we reconstruct the terahertz beam profiles in x–y coordinates from the numerical results in the cylindrical coordinates inside the PPLN before the incidence on the S₁ surface (see Fig. 7(a)). The refraction at the S₁, S₂ surfaces (see Fig. 7(a)) for p polarized light and the beam propagation in the quartz coupler is calculated by the angular spectrum method [46]. Figs. 14(a)–14(b) and Figs. 14(c)–14(d) correspond to terahertz spatial profiles generated with a pump beam size σ = 5 mm and σ = 1 mm, respectively by a two-stage setup.

 

Fig. 14 (a,b) and (c,d) represent terahetz spatial profiles in PPLN in x − y coordinate before S₁ surface (see Fig. 7(a)) generated by a two-stage system with pump beam size σ = 5 mm and σ = 1 mm respectively. (e–h) represent coupled out terahertz spatial profiles of (a–d) after S₂ surface in x′ − y coordinate in the air. The color bar represents fluence in the unite of J/m². Each figure window size is 15 mm × 15 mm. Terahertz out-coupling efficiency from (a–d) to (e–h) is η_a = 96.1%, η_b = 93.7%, η_c = 62.5%, η_d = 60.3% respectively.

Download Full Size | PDF

Note that with a pump beam size σ = 5 mm, the terahertz spatial profile retains its flat top nature, whereas with σ = 1 mm, the terahertz spatial profile reduces to a Gaussian as shown in Fig. 9. The corresponding transmitted spatial profiles after the S₂ surface in x′ − y coordinates in the air (see Fig. 7(a)) are shown in Figs. 14(e)–14(h) with the terahertz out-coupling efficiency η_a = 96.1%, η_b = 93.7%, η_c = 62.5%, η_d = 60.3% respectively. The influence of the S₁ and S₃ surfaces together with the propagation in the quartz coupler has a negligible effect on the pump spatial profile, i.e the spatial profile of the pump stays flat top but with a reduced size in the x′ direction i.e σ′ = σ tan(β₂) cos(β₁)/sin(α₁) ≈ 0.6σ. The beam size of the terahertz radiation in the x′ direction follows the same relation σ′_t = σ_t tan(β₂)/tan(α₁) ≈ 0.2σ_t. On the other hand, the refraction at S₁ and S₂ surfaces greatly changes the energy distribution of the flat top terahertz beams spatially (see Figs. 14(a)–14(b) and Figs. 14(e)–14(f)) while the Gaussian spatial profile remains unchanged (see Figs. 14(c)–14(d) and Figs. 14(g)–14(h)). The flat top function is a superposition of Hermite Gaussian modes in Cartesian coordinates. Thus, the final spatial profile of terahertz is a result of a superposition of all the Hermite Gaussian modes after refraction.

The terahertz beams coupled out by the QC could be combined by a power combiner with a specifically designed input coupler which matches the terahertz beam profile to the mode of the power combiner [67]. This investigation is beyond the scope of this article and will be presented in a separate study. For an ideal symmetric power combiner, the combination efficiency η_combine can be written as the following [68]:

4. Conclusion

We presented a scheme for generating terahertz pulses using a consecutive arrangement of PPLN crystals (multi-stage system). For this purpose, a 2-D numerical scheme assuming cylindrical symmetry is developed. The simulation includes the effects of DFG, SPM, self-focusing, beam diffraction, dispersion and terahertz absorption.

For a single stage, the optimal pump beam pulse duration for a 2-line pump spectrum configuration is 150 ps with effective length 25 mm and conversion efficiency η = 1.05%. With this pump pulse configuration, the terahertz spatial profile resembles the pump when the pump waist σ ≥ 3 mm. Additionally, we found that the intensity profile modulation has negligible influence on the terahertz generation efficiency, whereas a distortion in phase front strongly reduces the conversion efficiency. This suggests that care must be taken in circumventing phase induced efficiency deterioration for high energy terahertz generation.

For a multi-stage system, we found that by recycling the pump pulse with dispersion compensation, the efficiency can be greatly enhanced. Subsequently, a quartz output coupler for separating the terahertz beam and optical pump with high terahertz transmission is designed.

Specifically, for a two stage system, we predict generation of a terahertz pulse of 17.6 mJ energy with total conversion efficiency η_total = 1.6% at 0.3 THz for a 1.1 J input pump energy, where the generation efficiencies are η₁ = 1.0%, η₂ = 0.8%, the out-coupling efficiencies are η_a = 96.1%, η_b = 93.7% and combination efficiency is η_combine = 93.0%.

Appendix

5.1. Phase-Matching

In PPLN, the pump pulse is chosen to be linearly polarized along the extraordinary optical axis of LiNbO₃ in order to maximize the effective second order nonlinear coefficient $d_{33}=\frac{1}{2}{\chi }_{333}^{(2)}=-\frac{1}{4}{n}^4(\omega ){\gamma }_{333}=168\hspace{0.17em}\text{pm}/\text{V}$ [23,35,57,58].

(12)$$\left[\begin{array}{c}{P}_1^{{\mathit{NL}}_{(2)}}\\ P_2^{{\mathit{NL}}_{(2)}}\\ P_3^{{\mathit{NL}}_{(2)}}\end{array}\right]=2{∊}_0\left[\begin{array}{cccccc}0& 0& 0& 0& d_{15}& -d_{22}\\ -d_{22}& d_{22}& 0& d_{15}& 0& 0\\ d_{15}& d_{15}& d_{33}& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\left|E_1\right|}^2\\ {\left|E_2\right|}^2\\ {\left|E_3\right|}^2\\ E_2{E}_3^{*}+E_2^{*}{E}_3\\ E_1{E}_3^{*}+E_1^{*}{E}_3\\ E_1{E}_2^{*}+E_1^{*}{E}_2\end{array}\right]$$

In Eq. (12), subscripts 1 and 2 represent the ordinary axes and 3 represents the extraordinary axis [69]. Axes 1 and 2 are identical at the first order (n_o1 = n_o2), which is known as uniaxial birefringence. However, in the second order, the o axes are not identical. This explains why d_ij matrix elements are not identical on 1 and 2 axes. With the linearly polarized input pump pulse Ê = (0, 0, E₃) along the extraordinary axis, the second order polarization term driving the terahertz field is also linearly polarized along the same direction ($P_3^{{\text{NL}}_{(2)}}=2{∊}_0{d}_{33}{\left|E_3\right|}^2$).

The phase-matching condition shown in Eq. (13) ensures that the terahertz pulses generated at different positions of the PPLN crystal can add up in phase provided:

(13)$$\begin{array}{l}\mathrm{\Omega }\left[1/v_{\text{THz}}(\mathrm{\Omega })-1/v_{\text{g}}(\omega )\right]\frac{\mathrm{\Lambda }}{2}=\pi +2\pi N\\ N=0,1,2,3\dots \end{array}$$

where Λ is the period of PPLN, Ω is the terahertz angular frequency, ω is the pump angular frequency, v_THz is the phase velocity of terahertz pulse and v_g is the group velocity of the pump pulse. Note that for a given PPLN structure, there is more than one terahertz frequency fulfilling the phase-matching condition (Ω_N = Ω₀ × (1 + 2N))

5.2. Frequency dpendent self phase modulation

Frequency dependent SPM includes self-steepening which can be shown using the following derivation:

(14)$$\begin{array}{l}{\omega }_0ℱ\left\{{\left|A_{\text{op}}(t,r,z)\right|}^2{A}_{\text{op}}(t,r,z)-\frac{i}{{\omega }_0}\frac{\partial {\left|A_{\text{op}}(t,r,z)\right|}^2{A}_{\text{op}}(t,r,z)}{\partial t}\right\}\\ ={\omega }_0ℱ\left\{{\left|A_{\text{op}}(t,r,z)\right|}^2{A}_{\text{op}}(t,r,z)\right\}+(\omega -{\omega }_0)ℱ\left\{{\left|A_{\text{op}}(t,r,z)\right|}^2{A}_{\text{op}}(t,r,z)\right\}\\ =\omega ℱ\left\{{\left|A_{\text{op}}(t,r,z)\right|}^2{A}_{\text{p}}(t,r,z)\right\}.\end{array}$$

where the first and second terms correspond to the SPM effect at ω₀ and self-steepening, respectively.

5.3. Finite difference method

The finite difference methodology for solving Eq. (1) is presented. Equation (2) can be solved by following the same procedure. By defining $r_j^{+}=\frac{1}{2}(r_{j+1}+r_j)/\left(\mathrm{\Delta }{r}^2{r}_j\right)$, $r_j^{-}=\frac{1}{2}(r_{j-1}+r_j)/\left(\mathrm{\Delta }{r}^2{r}_j\right)$, where Δr is the step size in the r dimension and r_j represents a mesh point at a specific radial position, one can obtain the discretized form of the second order derivative as shown in Eq. (15).

(15)$$\begin{array}{ll}\frac{1}{r_j}\frac{\partial }{\partial r}\left(r_j\frac{\partial A_{\text{op}}({\omega }_m,r_j,z_k)}{\partial r}\right)\hfill & =r_j^{-}{A}_{\text{op}}({\omega }_m,r_{j-1},z_k)-\left(r_j^{-}+r_j^{+}\right)A_{\text{op}}({\omega }_m,r_j,z_k)\hfill \\ \hfill & +r_j^{+}{A}_{\text{op}}({\omega }_m,r_{j+1},z_k)\hfill \end{array}$$

Equation (15) can be written in a matrix multiplication form as shown on the right hand side of the last column in Fig. 15 where the r_i related matrix is denoted as a tridiagonal matrix M_R as shown in Eq. (16).

Equation (1) can be written as Eq. (17) where P^NL(ω_m, r_j, z_k) represents all the second and third nonlinear polarization terms.

(17)$$\begin{array}{ll}\frac{\partial A_{\text{op}}({\omega }_m,r_j,z_k)}{\partial z}\hfill & =P^{\text{NL}}({\omega }_m,r_j,z_k)+\frac{i}{2k({\omega }_m)}[-r_j^{-}{A}_{\text{op}}({\omega }_m,r_{j-1},z_k)\hfill \\ \hfill & +\left(r_j^{-}+r_j^{+}\right)A_{\text{op}}({\omega }_m,r_j,z_k)-r_j^{+}{A}_{\text{op}}({\omega }_m,r_{j+1},z_k)]\hfill \end{array}$$

Equation (17) can be written in matrix form as Fig. 15 where the 2-D matrix of angular frequency (ω) and radius (r) on the left hand side is updated by a low-storage Runge-Kutta method [51]. The electric field elements at the same position r_j with different frequencies are mixed by the nonlinear polarization P^NL. The electric field elements of the same frequency ω_m at different positions in r are mixed by M_R (e.g diffraction). The totally number of the mesh points in r dimension is N_r + 2. The extra 2 points A_op(ω, r₀, z), A_op(ω, r_Nr+1, z), marked with dashed boxes, are “ghost points”. The “ghost points” together with the two columns marked by dashed lines in M_R in Eq. 16 are used to construct boundary conditions.

 

Fig. 15 Illustration of the matrix representation of Eq. (17). The numbers of the mesh grids at r and ω dimensions are represented by N_r and N_ω respectively. Each color of the corresponding column vector represents the coupled wave equation at a specific frequency ω_m.

Download Full Size | PDF

5.4. Boundary conditions

The boundary conditions are enforced to the solution by the elements in the dashed boxes in Fig. 15, and M_R matrix in Eq. (16). For the mesh points in r dimension, the r = 0 singularity is avoided by defining r₀ = 3Δr, and thus $r_1^{-}=\frac{1}{2}(r_0+r_1)/\left(\mathrm{\Delta }{r}^2{r}_1\right)$. The “boundary condition” near r₁ is confined by the symmetry at the origin of the cylindrical coordinate itself (e.g ${\frac{\partial A_{\text{op}}(\omega ,r,z)}{\partial r}|}_{r=r_1}=0$). Thus,

(18)$$A_{\text{op}}(\omega ,r_0,z)=A_{\text{op}}(\omega ,r_1,z)$$

For the boundary condition at r_Nr, three different conditions can be applied. Note that as shown in Fig. 2, in the calculated area the polarization of the field is perpendicular to the boundary. Of course with the same numerical method, any polarization at the boundary can be considered. However, the boundary condition has to be changed accordingly.

1. Dirichlet boundary.

   This condition represents the interface from PPLN crystal (n_op) to air (n). Due to Gauss’s law in maxwell equation ∬D⃗ · dS⃗ = 0. One can get

   (19)$$A_{\text{op}}(\omega ,r_{N_r+1},z)=A_{\text{op}}(\omega ,r_{N_r},z)\frac{n_{\text{op}}^2{r}_{N_r}}{n^2{r}_{N_r+1}}$$
   If the field polarization is parallel to the contact surface of the dielectric (PPLN) and the metal cladding outside the PPLN, one can also apply a reflective boundary condition.
   (20)$$A_{\text{op}}(\omega ,r_{N_r+1},z)=0$$
   This is equivalent to the situation where electric field is reflected by a mirror.

2. Transparent boundary [70].

   This boundary represents the situation where the boundary has no effect on the field as the field propagates through the calculation boundary. This kind of boundary is also equivalent to the situation when Eq. (1) and (2) are solved by angular spectrum method. The condition can be written as [70].

   (21)$$A_{\text{op}}(\omega ,r_{N_r+1},z_{k+1})=A_{\text{op}}(\omega ,r_{N_r},z_{k+1})\frac{A_{\text{op}}(\omega ,r_{N_r+1},z_k)r_{N_r}}{A_{\text{op}}(\omega ,r_{N_r},z_k)r_{N_r+1}}$$

5.5. Pulse duration in terms of effective length

As section (3.1) suggests, the effective length is related to the pump pulse duration. In order to get more insight into this problem, we begin with a 1-D undepleted model. Consequently, Eq. (2) can be written as Eq. (22).

(22)$$\begin{array}{ll}{A}_{\text{THz}}(\mathrm{\Omega },z)=\hfill & \frac{-i\mathrm{\Omega }{\chi }_{333}^{(2)}}{\pi n(\mathrm{\Omega })c}{\int }_0^{\infty }{A}_{\text{op}}(\omega +\mathrm{\Omega },r,z)A_{\text{op}}^{*}(\omega ,r,z)d\omega \left[\frac{e^{i\delta z}-e^{\frac{-\alpha (\mathrm{\Omega })z}{2}}}{\frac{\alpha (\mathrm{\Omega })}{2}+i\delta }\right]\hfill \\ =\hfill & P^{\text{NL}}(\mathrm{\Omega })\left[\frac{e^{i\delta z}-e^{\frac{-\alpha (\mathrm{\Omega })z}{2}}}{\frac{\alpha (\mathrm{\Omega })}{2}+i\delta }\right]\hfill \end{array}$$

where,

(23)$$\delta =\left[k(\mathrm{\Omega })+k(\omega )-k(\omega +\mathrm{\Omega })-\frac{2\pi }{\mathrm{\Lambda }}\right]=\left[\frac{\mathrm{\Omega }}{c}\left(n(\mathrm{\Omega })-n_{\text{g}}(\omega )\right)-\frac{2\pi }{\mathrm{\Lambda }}\right]\approx \frac{\left(n({\mathrm{\Omega }}_0)-n_{\text{g}}({\omega }_0)\right)}{c{\tau }_{\text{FWHM}}}.$$

In order to study the distance L₀ and L_eff, we calculate the second order derivative of |A_THz|² which is proportional to the efficiency η. By calculating d²|A_THz(Ω, z)|²/d²z = 0, one can get

(24)$$\begin{array}{l}{e}^{-\frac{\alpha z}{2}}=\left[\frac{1}{2}\text{cos}(\delta z)+2\frac{\delta }{\alpha }\text{sin}(\delta z)-2\frac{{\delta }^2}{{\alpha }^2}\text{cos}(\delta z)\right]\\ \{\begin{array}{l}\text{for}\hspace{0.17em}\text{short}\hspace{0.17em}\text{pump}\hspace{0.17em}\text{pulse},\hspace{0.17em}\frac{\delta }{\alpha }\gg 1,2\frac{\delta }{\alpha }\text{sin}(\delta z)=2\frac{{\delta }^2}{{\alpha }^2}\text{cos}(\delta z)\hfill \\ \Rightarrow L_0={\text{tan}}^{-1}\left(\frac{\delta }{\alpha }\right)/\delta ,\hfill \\ \text{for}\hspace{0.17em}\text{long}\hspace{0.17em}\text{pump}\hspace{0.17em}\text{pulse},\hspace{0.17em}\frac{\delta }{\alpha }\ll 1,\hspace{0.17em}\hspace{0.17em}{e}^{-\frac{\alpha z}{2}}=\frac{1}{2}\text{cos}(\delta z)\approx \frac{1}{2}\hfill \\ \Rightarrow L_0=2\text{ln}(2)/\alpha ,\hspace{0.17em}\hspace{0.17em}{L}_{\text{eff}}=\frac{2}{\alpha }\text{ln}\left(\frac{2}{1-\sqrt{1-e^{-1}}}\right)\hfill \end{array}\end{array}$$

Funding

European Union’s Seventh Framework Programme (FP7/2007-2013) through the Synergy Grant AXSIS (609920); The Hamburg Center for Ultrafast Imaging – Structure, Dynamics and Control of Matter at the Atomic Scale (CUI, DFG-EXC1074).

Acknowledgments

We thank the DESY HPC team for providing access and computation time on the Maxwell cluster. Lu Wang would like to thank IMPRS (International Max Planck Research School for Ultrafast Imaging & Structural Dynamics). The authors additionally would like to thank Dr. Nicholas Matlis and Dr. Steve Aplin for the supportive and helpful discussions.

References

1. A. Markelz, A. Roitberg, and E. J. Heilweil, “Pulsed terahertz spectroscopy of DNA, bovine serum albumin and collagen between 0.1 and 2.0 THz,” Chem. Phys. Lett 320, 42–48 (2000). [CrossRef]  

2. M. Nagel, P. Haring Bolivar, M. Brucherseifer, H. Kurz, A. Bosserhoff, and R. Büttner, “Integrated THz technology for label-free genetic diagnostics,” Appl. Phys. Lett 80, 154–156 (2002). [CrossRef]  

3. K. Young, Q. Wen, S. Hanany, H. Imada, J. Koch, T. Matsumura, O. Suttmann, and V. Schütz, “Broadband millimeter-wave anti-reflection coatings on silicon using pyramidal sub-wavelength structures,” J. Appl. Phys 121, 213103 (2017). [CrossRef]  

4. T. Kampfrath, A. Sell, G. Klatt, A. Pashkin, S. Mährlein, T. Dekorsy, M. Wolf, M. Fiebig, A. Leitenstorfer, and R. Huber, “Coherent terahertz control of antiferromagnetic spin waves,” Nat. Photonics 5, 31–34 (2011). [CrossRef]  

5. T. Kubacka, J. Johnson, M. Hoffmann, C. Vicario, S. de Jong, P. Beaud, S. Grübel, S.-W. Huang, L. Huber, L. Patthey, Y.-D. Chuang, J. J. Turner, G. L. Dakovski, W.-S. Lee, M. P. Minitti, W. Schlotter, R. G. Moore, C. P. Hauri, S. M. Koohpayeh, V. Scagnoli, G. Ingold, S. L. Johnson, and U. Staub, “Large-amplitude spin dynamics driven by a THz pulse in resonance with an electromagnon,” Science 343, 1333–1336 (2014). [CrossRef]   [PubMed]  

6. E. A. Nanni, W. R. Huang, K.-H. Hong, K. Ravi, A. Fallahi, G. Moriena, R. D. Miller, and F. X. Kärtner, “Terahertz-driven linear electron acceleration,” Nat. Commun 6, 8486 (2015). [CrossRef]   [PubMed]  

7. L. J. Wong, A. Fallahi, and F. X. Kärtner, “Compact electron acceleration and bunch compression in THz waveguides,” Opt. Express 21, 9792–9806 (2013). [CrossRef]   [PubMed]  

8. D. Zhang, A. Fallahi, M. Hemmer, X. Wu, M. Fakhari, Y. Hua, H. Cankaya, A.-L. Calendron, L. E. Zapata, N. H. Matlis, and F. X. Kaertner, “Segmented terahertz electron accelerator and manipulator (STEAM),” Nat. photonics pp. 336–342 (2018). [CrossRef]   [PubMed]  

9. X. Xie, J. Dai, and X.-C. Zhang, “Coherent control of THz wave generation in ambient air,” Phys. Rev. Lett 96, 075005 (2006). [CrossRef]   [PubMed]  

10. M. Kress, T. Löffler, S. Eden, M. Thomson, and H. G. Roskos, “Terahertz-pulse generation by photoionization of air with laser pulses composed of both fundamental and second-harmonic waves,” Opt. Lett. 29, 1120–1122 (2004). [CrossRef]   [PubMed]  

11. P. U. Jepsen, R. H. Jacobsen, and S. Keiding, “Generation and detection of terahertz pulses from biased semiconductor antennas,” J. Opt. Soc. Am. B 13, 2424–2436 (1996). [CrossRef]  

12. A. Dreyhaupt, S. Winnerl, T. Dekorsy, and M. Helm, “High-intensity terahertz radiation from a microstructured large-area photoconductor,” Appl. Phys. Lett 86, 121114 (2005). [CrossRef]  

13. J. Hebling, G. Almasi, I. Z. Kozma, and J. Kuhl, “Velocity matching by pulse front tilting for large-area THz-pulse generation,” Opt. Express 10, 1161–1166 (2002). [CrossRef]   [PubMed]  

14. B. K. Ofori-Okai, P. Sivarajah, W. R. Huang, and K. A. Nelson, “THz generation using a reflective stair-step echelon,” Opt. Express 24, 5057–5068 (2016). [CrossRef]   [PubMed]  

15. K. Ravi, B. K. Ofori-Okai, P. Sivarajah, W. R. Huang, F. X. Kärtner, and K. A. Nelson, “Broadband terahertz generation with a stair-step echelon,” in Conference on Lasers and Electro-Optics, (Optical Society of America, 2017), p. STu3J.1.

16. L. Pálfalvi, G. Tóth, L. Tokodi, Z. Márton, J. A. Fülöp, G. Almási, and J. Hebling, “Numerical investigation of a scalable setup for efficient terahertz generation using a segmented tilted-pulse-front excitation,” Opt. Express 25, 29560–29573 (2017). [CrossRef]   [PubMed]  

17. K. Vijayraghavan, Y. Jiang, M. Jang, A. Jiang, K. Choutagunta, A. Vizbaras, F. Demmerle, G. Boehm, M. C. Amann, and M. A. Belkin, “Broadly tunable terahertz generation in mid-infrared quantum cascade lasers,” Nat. Commun 4, 2021 (2013). [CrossRef]   [PubMed]  

18. B. S. Williams, “Terahertz quantum-cascade lasers,” Nat. Photonics 1, 517–525 (2007). [CrossRef]  

19. M. A. Belkin, F. Capasso, A. Belyanin, D. L. Sivco, A. Y. Cho, D. C. Oakley, C. J. Vineis, and G. W. Turner, “Terahertz quantum-cascade-laser source based on intracavity difference-frequency generation,” Nat. Photonics 1, 288–292 (2007). [CrossRef]  

20. Y.-S. Lee, Principles of terahertz science and technology (Springer Science & Business Media, 2009).

21. C. P. Hauri, “Intense laser-based THz sources for XFEL experiments,” in Conference on Lasers and Electro-Optics, (Optical Society of America, 2017), p. JTu3L.3.

22. Z. Wu, A. S. Fisher, M. C. Hoffmann, S. Bonetti, D. Higley, and H. Durr, “THz light source at slac facet user facility,” in 39th International Conference on Infrared, Millimeter, and Terahertz waves (IEEE, 2014), pp. 1–2.

23. K. L. Vodopyanov, “Optical generation of narrow-band terahertz packets in periodically-inverted electro-optic crystals: conversion efficiency and optimal laser pulse format,” Opt. Express 14, 2263–2276 (2006). [CrossRef]   [PubMed]  

24. Z. Chen, X. Zhou, C. A. Werley, and K. A. Nelson, “Generation of high power tunable multicycle teraherz pulses,” Appl. Phys. Lett 99, 071102 (2011). [CrossRef]  

25. M. Tonouchi, “Cutting-edge terahertz technology,” Nat. Photonics 1, 97–105 (2007). [CrossRef]  

26. K. Vodopyanov, M. Fejer, X. Yu, J. Harris, Y.-S. Lee, W. C. Hurlbut, V. Kozlov, D. Bliss, and C. Lynch, “Terahertz-wave generation in quasi-phase-matched GaAs,” Appl. Phys. Lett 89, 141119 (2006). [CrossRef]  

27. A. Rice, Y. Jin, X. Ma, X.-C. Zhang, D. Bliss, J. Larkin, and M. Alexander, “Terahertz optical rectification from< 110> zinc-blende crystals,” Appl. Phys. Lett 64, 1324–1326 (1994). [CrossRef]  

28. K. L. Vodopyanov, “Optical THz-wave generation with periodically-inverted GaAs,” Laser & Photonics Rev. 2, 11–25 (2008). [CrossRef]  

29. J. Lu, H. Y. Hwang, X. Li, S.-H. Lee, O.-P. Kwon, and K. A. Nelson, “Tunable multi-cycle THz generation in organic crystal hmq-tms,” Opt. Express 23, 22723–22729 (2015). [CrossRef]   [PubMed]  

30. Y.-S. Lee, T. Meade, V. Perlin, H. Winful, T. Norris, and A. Galvanauskas, “Generation of narrow-band terahertz radiation via optical rectification of femtosecond pulses in periodically poled lithium niobate,” Appl. Phys. Lett 76, 2505–2507 (2000). [CrossRef]  

31. D. Dietze, K. Unterrainer, and J. Darmo, “Dynamically phase-matched terahertz generation,” Opt. Lett. 37, 1047–1049 (2012). [CrossRef]   [PubMed]  

32. P. Nugraha, G. Krizsán, G. Polónyi, M. Mechler, J. Hebling, G. Tóth, and J. Fülöp, “Efficient semiconductor multicycle terahertz pulse source,” J. Phys. B: At. Mol. Opt. Phys. 51, 094007 (2018). [CrossRef]  

33. K. Kawase, T. Hatanaka, H. Takahashi, K. Nakamura, T. Taniuchi, and H. Ito, “Tunable terahertz-wave generation from dast crystal by dual signal-wave parametric oscillation of periodically poled lithium niobate,” Opt. Lett. 25, 1714–1716 (2000). [CrossRef]  

34. A. Schneider, M. Neis, M. Stillhart, B. Ruiz, R. U. Khan, and P. Günter, “Generation of terahertz pulses through optical rectification in organic dast crystals: theory and experiment,” J. Opt. Soc. Am. B 23, 1822–1835 (2006). [CrossRef]  

35. K. Ravi, D. N. Schimpf, and F. X. Kärtner, “Pulse sequences for efficient multi-cycle terahertz generation in periodically poled lithium niobate,” Opt. Express 24, 25582–25607 (2016). [CrossRef]   [PubMed]  

36. G. Cirmi, M. Hemmer, K. Ravi, F. Reichert, L. E. Zapata, A.-L. Calendron, H. Çankaya, F. Ahr, O. D. Mücke, N. H. Matlis, and F. X. Kärtner, “Cascaded second-order processes for the efficient generation of narrowband terahertz radiation,” J. Phys. B 50, 044002 (2017). [CrossRef]  

37. G. Tóth, J. A. Fülöp, and J. Hebling, “Periodically intensity-modulated pulses by optical parametric amplification for multicycle tunable terahertz pulse generation,” Opt. Express 25, 28258–28272 (2017). [CrossRef]  

38. F. Ahr, S. W. Jolly, N. H. Matlis, S. Carbajo, T. Kroh, K. Ravi, D. N. Schimpf, J. Schulte, H. Ishizuki, T. Taira, A. R. Maier, and F. X. Kärtner, “Narrowband terahertz generation with chirped-and-delayed laser pulses in periodically poled lithium niobate,” Opt. Lett. 42, 2118–2121 (2017). [CrossRef]   [PubMed]  

39. J. Manley and H. Rowe, “Some general properties of nonlinear elements-part I. general energy relations,” Proc. IRE 44, 904–913 (1956). [CrossRef]  

40. M. Cronin-Golomb, “Cascaded nonlinear difference-frequency generation of enhanced terahertz wave production,” Opt. Lett. 29, 2046–2048 (2004). [CrossRef]   [PubMed]  

41. M. C. Hoffmann, K.-L. Yeh, J. Hebling, and K. A. Nelson, “Efficient terahertz generation by optical rectification at 1035 nm,” Opt. Express 15, 11706–11713 (2007). [CrossRef]   [PubMed]  

42. J.-i. Shikata, M. Sato, T. Taniuchi, H. Ito, and K. Kawase, “Enhancement of terahertz-wave output from LiNbO₃ optical parametric oscillators by cryogenic cooling,” Opt. Lett. 24, 202–204 (1999). [CrossRef]  

43. S.-W. Huang, E. Granados, W. R. Huang, K.-H. Hong, L. E. Zapata, and F. X. Kärtner, “High conversion efficiency, high energy terahertz pulses by optical rectification in cryogenically cooled lithium niobate,” Opt. Lett. 38, 796–798 (2013). [CrossRef]   [PubMed]  

44. Y.-S. Lee, T. Meade, M. DeCamp, T. Norris, and A. Galvanauskas, “Temperature dependence of narrow-band terahertz generation from periodically poled lithium niobate,” Appl. Phys. Lett 77, 1244–1246 (2000). [CrossRef]  

45. C.-W. Chen, Y.-S. Lin, J. Y. Huang, C.-S. Chang, C.-L. Pan, L. Yan, and C.-K. Lee, “Generation and spectral manipulation of coherent terahertz radiation with two-stage optical rectification,” Opt. Express 16, 14294–14303 (2008). [CrossRef]   [PubMed]  

46. J. W. Goodman, Introduction to Fourier optics (Roberts and Company Publishers, 2005).

47. K. Ravi, W. R. Huang, S. Carbajo, X. Wu, and F. Kärtner, “Limitations to THz generation by optical rectification using tilted pulse fronts,” Opt. Express 22, 20239–20251 (2014). [CrossRef]   [PubMed]  

48. U. Schwarz and M. Maier, “Asymmetric raman lines caused by an anharmonic lattice potential in lithium niobate,” Phys. Rev. B 55, 11041–11044 (1997). [CrossRef]  

49. H. Ishizuki and T. Taira, “Half-joule output optical-parametric oscillation by using 10-mm-thick periodically poled mg-doped congruent LiNbO₃,” Opt. Express 20, 20002–20010 (2012). [CrossRef]   [PubMed]  

50. K. Ravi, W. R. Huang, S. Carbajo, E. A. Nanni, D. N. Schimpf, E. P. Ippen, and F. X. Kärtner, “Theory of terahertz generation by optical rectification using tilted-pulse-fronts,” Opt. Express 23, 5253–5276 (2015). [CrossRef]   [PubMed]  

51. J. Williamson, “Low-storage Runge-Kutta schemes,” J. Comput. Phys 35, 48–56 (1980). [CrossRef]  

52. A. Siegman, “Quasi fast Hankel transform,” Opt. Lett. 1, 13–15 (1977). [CrossRef]   [PubMed]  

53. M. Cree and P. Bones, “Algorithms to numerically evaluate the Hankel transform,” Comput. & Math. with Appl. 26, 1–12 (1993). [CrossRef]  

54. X. Wu, C. Zhou, W. R. Huang, F. Ahr, and F. X. Kärtner, “Temperature dependent refractive index and absorption coefficient of congruent lithium niobate crystals in the terahertz range,” Opt. Express 23, 29729–29737 (2015). [CrossRef]   [PubMed]  

55. D. H. Jundt, “Temperature-dependent sellmeier equation for the index of refraction, n e, in congruent lithium niobate,” Opt. Lett. 22, 1553–1555 (1997). [CrossRef]  

56. R. DeSalvo, A. A. Said, D. J. Hagan, E. W. Van Stryland, and M. Sheik-Bahae, “Infrared to ultraviolet measurements of two-photon absorption and n₂ in wide bandgap solids,” IEEE J. Quantum Electron 32, 1324–1333 (1996). [CrossRef]  

57. R. Weis and T. Gaylord, “Lithium niobate: summary of physical properties and crystal structure,” Appl. Phys. A: Mater. Sci. & Process. 37, 191–203 (1985). [CrossRef]  

58. J. Hebling, K.-L. Yeh, M. C. Hoffmann, B. Bartal, and K. A. Nelson, “Generation of high-power terahertz pulses by tilted-pulse-front excitation and their application possibilities,” J. Opt. Soc. Am. B 25, B6–B19 (2008). [CrossRef]  

59. F. Kärtner, N. Matuschek, T. Schibli, U. Keller, H. Haus, C. Heine, R. Morf, V. Scheuer, M. Tilsch, and T. Tschudi, “Design and fabrication of double-chirped mirrors,” Opt. Lett. 22, 831–833 (1997). [CrossRef]   [PubMed]  

60. K. Kawase, M. Sato, K. Nakamura, T. Taniuchi, and H. Ito, “Unidirectional radiation of widely tunable THz wave using a prism coupler under noncollinear phase matching condition,” Appl. Phys. Lett 71, 753–755 (1997). [CrossRef]  

61. M. Naftaly and R. E. Miles, “Terahertz time-domain spectroscopy for material characterization,” Proc. IEEE 95, 1658–1665 (2007). [CrossRef]  

62. D. Grischkowsky, S. Keiding, M. Van Exter, and C. Fattinger, “Far-infrared time-domain spectroscopy with terahertz beams of dielectrics and semiconductors,” J. Opt. Soc. Am. B 7, 2006–2015 (1990). [CrossRef]  

63. I. Shoji, T. Kondo, A. Kitamoto, M. Shirane, and R. Ito, “Absolute scale of second-order nonlinear-optical coefficients,” J. Opt. Soc. Am. B 14, 2268–2294 (1997). [CrossRef]  

64. G. Ghosh, “Dispersion-equation coefficients for the refractive index and birefringence of calcite and quartz crystals,” Opt. Commun 163, 95–102 (1999). [CrossRef]  

65. A. Parent, M. Morin, and P. Lavigne, “Propagation of super-gaussian field distributions,” Opt. Quantum Electron. 24, S1071–S1079 (1992). [CrossRef]  

66. A. Chafiq, Z. Hricha, and A. Belafhal, “Flat-topped mathieu-gauss beam and its transformation by paraxial optical systems,” Opt. Commun 278, 142–146 (2007). [CrossRef]  

67. N. I. Landy and W. J. Padilla, “Guiding light with conformal transformations,” Opt. Express 17, 14872–14879 (2009). [CrossRef]   [PubMed]  

68. M. Gupta, “Power combining efficiency and its optimisation,” IEE Proc. H (Microwaves, Antennas Propagation) 139, 233–238 (1992). [CrossRef]  

69. R. W. Boyd, Nonlinear optics (Elsevier, 2003).

70. G. R. Hadley, “Transparent boundary condition for the beam propagation method,” IEEE J. Quantum Electron 28, 363–370 (1992). [CrossRef]