## Abstract

A model for terahertz (THz) generation by optical rectification using tilted-pulse-fronts is developed. It simultaneously accounts for in two spatial dimensions (2-D) (i) the spatio-temporal variations of the optical pump pulse imparted by the tilted-pulse-front setup, (ii) the nonlinear coupled interaction of THz and optical radiation, (iii) self-phase modulation and (iv) stimulated Raman scattering. The model is validated by quantitative agreement with experiments and analytic calculations. We show that the optical pump beam is significantly broadened in the transverse-momentum (*k _{x}*) domain as a consequence of its spectral broadening due to THz generation. In the presence of this large frequency and transverse-momentum (or angular) spread, group velocity dispersion causes a spatio-temporal break-up of the optical pump pulse which inhibits further THz generation. The implications of these effects on energy scaling and optimization of optical-to-THz conversion efficiency are discussed. This suggests the use of optical pump pulses with elliptical beam profiles for large optical pump energies. Furthermore, it is seen that optimization of the setup is highly dependent on optical pump conditions. Trade-offs in optimizing the optical-to-THz conversion efficiency on the spatial and spectral properties of THz radiation are discussed to guide the development of such sources.

© 2015 Optical Society of America

## 1. Introduction

Terahertz (THz) sources are characterized by wavelengths roughly hundred times larger than optical and ten times smaller than radio frequency sources. In particular, short pulses of THz light with high peak field strengths are attractive for a number of unprecedented applications. They are uniquely amenable to probing and controlling material properties in a variety of systems such as superconductors, ferro/anti-ferromagnets etc [1–4]. In addition, they are attractive for compact charged particle-acceleration [5–9] and as drivers for other sources of radiation such as in high-harmonic generation [10,11].

Of various high field THz generation modalities, difference frequency generation between spectral components of a femtosecond optical pulse in nonlinear crystals or optical rectification (OR) has emerged as the most efficient approach [12,13] and resulted in the highest THz pulse energies [13] to date. Within this category, OR using tilted-pulse-fronts in lithium niobate is of particular interest due to its compatibility with optical pumping by widely available 800 nm and 1 μm sources. This approach was developed in [14–18] among other works as a means to achieve phase-matching in materials with large disparities between THz and optical refractive indices.

In this approach, an optical pump pulse is angularly dispersed to produce an intensity front which is tilted with respect to its propagation direction. THz radiation propagating perpendicular to this tilted intensity front, or tilted-pulse-front (TPF), is then generated. Since the optical and THz radiation travel different distances in the same time, the difference between optical and THz refractive indices is compensated and phase-matching is achieved. OR using TPFs in lithium niobate has resulted in optical-to-THz conversion efficiencies (henceforth referred to as conversion efficiency) in excess of 1% [12,19] and THz pulse energies of 0.4 mJ [20].Therefore, the approach is promising for the development of laboratory scale THz sources with THz pulse energies much greater than the mJ level. Comprehensive theoretical models to aid understanding and quantitatively predict the performance of such systems are therefore of interest. The requisites of a physically accurate model and the current state of theory are described below.

As a consequence of the angular dispersion of the optical pump pulse in OR using TPF, various frequency components of the optical pump pulse spectrum are spatially separated. This is tantamount to having different spectral bandwidths, pulse durations and average frequency at each spatial location. These effects are referred to as spatio-temporal variations [21] and affect the properties of the generated THz radiation. Secondly, since the generated THz propagates perpendicular to the TPF, the optical pump and THz radiation propagate non-collinearly. Most importantly, as THz radiation is generated, it is accompanied by a dramatic cascaded frequency down-shift and spectral broadening of the optical pump pulse spectrum (cascading effects). On one hand, cascading is responsible for conversion efficiencies that exceed the Manley-Rowe limit. However, in the presence of group velocity dispersion due to angular dispersion (GVD-AD) and material dispersion (GVD-MD), the increased optical bandwidth inhibits further THz generation [22–24]. A comprehensive theoretical model should therefore be able to account for all of the above effects. This would require a simultaneous solution of optical and THz electric fields (henceforth referred to as field) in at least two spatial dimensions (2-D). In addition, spatio-temporal variations imparted by the specific TPF setup would also have to be considered.

Previously presented models are broadly comprised of (i) 1-D and 2-D spatial models without the inclusion of cascading effects (i.e., nonlinear coupling between THz and optical radiation is not considered) [25–29] and (ii) 1-D spatial models which account for cascading effects [22,30]. Models in category (i) overestimate the possible conversion efficiencies by not considering the spectral re-shaping of the optical pump pulse [22]. While models in category (ii) can phenomenologically address this shortcoming, they are still too simplistic to accurately describe experiments (e.g. implications of pump beam size, spatial frequency variations of THz radiation, THz beam propagation properties cannot be accounted for).

In this paper, we present the formulation of a 2-D model which simultaneously considers the spatio-temporal variations of the optical pump pulse, cascading effects, SPM and stimulated Raman scattering (SRS), angular and material dispersion, THz absorption as well as geometry of the nonlinear crystal. It can therefore account for the effects of finite beam size, spatial walk-off, spatial frequency variations and beam propagation which is not possible with our previous 1-D formulation [22]. The model is applicable to the simulation of a variety of OR systems with different TPF setups, crystal geometries and optical pump pulse formats. Despite its complexity, the model is formulated so that it can be solved efficiently to enable parametric studies. This is achieved by the use of appropriate simplifying assumptions, co-ordinate transformations and Fourier decomposition as outlined in Sections 2 and 4.

In Section 2, we provide an overview of THz generation using TPF and the general approach of our model. In Section 3, we provide an illustrative simulation result using the developed model. In particular, a discussion from transverse momentum (*k _{x}*) and time domain viewpoints is introduced for the first time. It is seen that the generation of THz results in the broadening of the optical pump pulse in both frequency and transverse-momentum domains. In the presence of this increased frequency and transverse-momentum spread, GVD-AD and GVD-MD cause a spatio-temporal break-up of the optical pump pulse which inhibits further THz generation. These descriptions serve to motivate the theoretical formulation of the complete TPF THz generation system, which is presented in Section 4. In Section 5, we validate the model by comparisons to experiments and analytic calculations. The impact of imaging errors on conversion efficiency is shown quantitatively. It is seen that small perturbations to the optimal imaging configuration can result in sizeable degradation of conversion efficiency. Insights into broadening of the optical spectrum are provided. In Section 6, we discuss the meaning of effective propagation length in 2-D. This is then used to discuss scaling to large optical pump energies and optimization of conversion efficiency. It is seen that the optimization of conversion-efficiency is highly dependent on the optical pump parameters. Finally, we highlight the trade-offs incurred while optimizing the conversion efficiency on spatial and spectral properties of THz radiation. We conclude in Section 7. This paper thus provides a basis for constructing sources customized optimally for various applications.

## 2. THz generation using tilted-pulse-fronts: overview of physics and modelling approach

The overall schematic of our approach is depicted in Fig. 1. An optical pump pulse with the input electric field described by the complex variable${E}_{op}^{in}(\omega ,{x}_{0},{z}_{0})$at angular frequency *ω*, propagating in the *z _{0}* direction is incident on a TPF setup. In Fig. 1(a), a commonly employed setup incorporating a diffraction grating and single lens is depicted. However, the model is applicable to a variety of TPF setups (e.g. employing telescope and diffraction grating, contact grating etc.). In Fig. 1(a), electric fields at two angular frequencies

*ω*and

*ω + Ω*are depicted for convenience although there are many more frequency components.

An optical pulse with a TPF is typically angularly dispersed [21]. Therefore, various frequency components of the emergent optical field described by ${E}_{op}^{out}(\omega ,{x}_{0},{z}_{0})$ propagate at different angles. This is depicted in Fig. 1(a) as spectral components at *ω* and *ω + Ω* emerge with wave vectors (henceforth referred to as momentum) $\overrightarrow{k}(\omega )$and $\overrightarrow{k}(\omega +\Omega )$ respectively inside the nonlinear crystal which is in the shape of a right-angled prism. In the time domain, such an angularly dispersed pulse has an intensity profile tilted with respect to its propagation direction as displayed by the red ellipses in Fig. 1(a). This gives rise to the terminology of ‘tilted-pulse-front ’. These red ellipses make an angle *π/2-γ* with respect to the propagation direction of the optical pulse, where *γ is* termed the pulse-front-tilt angle. Since OR is intra-pulse difference frequency generation (DFG), in the phase-matched condition, the momentum of the generated THz (at angular frequency *Ω*) is given by $\overrightarrow{k}(\Omega )=\overrightarrow{k}(\omega +\Omega )-\overrightarrow{k}(\omega )$ as depicted by the red arrow in Fig. 1(a). The generated THz then emerges at an angle *γ* with respect to the direction of propagation of the optical pump pulse and exits approximately normal to the output facet of the nonlinear crystal that is specifically cut to enable this. Since the various frequency components are angularly separated, there is a spatial variation in the average frequency (spatial-chirp) as well as amplitude across the beam profile. As a consequence, the pulse durations and bandwidths at various points in space are different, which affects the spatial and spectral properties of the generated THz pulse. In our model, we consider these effects by applying an analytic formulation of dispersive ray pulse matrices from [31] in Section 4.1. By invoking this analytic approach, large propagation distances of arbitrary TPF setups are accounted for by straightforward matrix multiplications which efficiently models the various spatio-temporal variations associated with the optical pulse. It then supplies the optical pump field${E}_{op}^{out}(\omega ,{x}_{0},{z}_{0}=0)$at the input face of the nonlinear crystal as shown in Fig. 1(b). By appropriately merging this analytic description of the TPF setup with the numerical solution of the coupled nonlinear wave equations, the model directly maps experimental conditions to the properties of the generated THz radiation. For example, in Section 5 we are able to quantitatively predict the trade-off in efficiency when the system is not at the optimal imaging condition.

The incident optical field excites a polarization in the nonlinear material to drive the generation of THz radiation. The generated THz in turn influences the propagation of the optical field and vice versa. Thus, the evolution of the optical and THz fields is described by a solution of a system of 2-D coupled nonlinear wave equations in the (*z-x*) co-ordinate system depicted in Fig. 1(b). We effectively reduce this system of 2-D wave equations to 1-D by employing Fourier decomposition. This enables a massive parallelization of the solution and improves the efficiency of the solution. In our approach, the (z*-x*) co-ordinate system is rotated with respect to (*z _{0}-x_{0}*) by an angle

*α*, which is the apex angle of the nonlinear crystal. The angle

*α*is approximately equal to the pulse-front-tilt angle

*γ*. The rotated co-ordinate system then has two key advantages. Firstly, in this set of axes, the THz radiation has small transverse-momentum components, i.e.

*k*, which relaxes the constraints on spatial resolution

_{x}~0*Δx*and consequently alleviates computational cost. This is particularly useful because the dimensions of the nonlinear crystals are on the order of ~cm, which results in large computational domains. Secondly, it makes it convenient to include the transmission of THz radiation at the crystal boundary.

In Fig. 1(b), we delineate how we consider the geometry of the crystal. An extended Cartesian space in the (*z-x*) co-ordinate system, uniformly filled with material of refractive index$n(\omega )$is considered. This is an auxiliary configuration for the numerical computation which mimics the situation in the actual experiment where the nonlinear material is surrounded by air. Only regions of the computational space physically occupied by the crystal would have a non-zero value of second order susceptibility${\chi}_{eff}^{(2)}(x,z)$ as shown in Fig. 1(b). If the length of the input crystal face is *L* and the optical field is incident at a distance *h* from the apex, the computational region extends from $\left(-(L-h)\mathrm{sin}\alpha ,-(L-h)\mathrm{cos}\alpha \right)$to$\left(h\mathrm{sin}\alpha ,h\mathrm{cos}\alpha \right)$ as shown in Fig. 1(b). The initial optical field profile along the line$z=-(L-h)\mathrm{sin}\alpha $can be calculated analytically by back-propagating the optical field calculated at the input crystal face at *z _{0}* = 0. In this model, we assume that the THz field is zero at the beginning of the computational space and consider only a single passage of the optical and THz beams through the crystal. In the limit of relatively thick crystals where the reflected THz energy is absorbed (absorption length at 300K is~2 mm) or with the use of THz anti-reflection coatings, this approximation is well justified.

## 3. Spatio-temporal break-up of the optical pump pulse due to THz generation

In order to illustrate the essential physics governing THz generation via OR using TPF, we first provide as an example, a solution of our model using the following parameters. A transform limited optical pump pulse, centred at a wavelength ${\lambda}_{0}$ = 1030 nm with pulsewidth of 0.5 ps full width at half maximum, peak intensity of 40 GW/cm^{2} and e* ^{−2}* beam radius ${w}_{in}$ = 2.5 mm is incident at

*h*= 1.5 mm from the apex of the nonlinear crystal with crystal temperature

*T*= 300 K. Figure 2(a), depicts the normalized optical fluence $\underset{0}{\overset{\infty}{\int}}{\left|{E}_{op}(\omega ,x,z)\right|}^{2}d\omega$ (where ${E}_{op}(\omega ,x,z)$is the optical field at angular frequency ω) and THz fluence $\underset{0}{\overset{\infty}{\int}}{\left|{E}_{THz}(\Omega ,x,z)\right|}^{2}d\Omega$ (where ${E}_{THz}(\Omega ,x,z)$is the THz field at angular frequency Ω). The crystal orientation and geometric parameters are as defined in the

*z-x*co-ordinate system described in Fig. 1(b). The cross sectional area of the various beams are defined by the two transverse directions x and

*y*(out of the plane of the paper). In our model, we solve for the frequency-domain optical field ${E}_{op}(\omega ,x,z)$and THz field${E}_{THz}(\Omega ,x,z)$. Fluences which are the aggregate of the magnitude squared of these fields over all frequencies are readily related to conversion efficiency and are hence depicted in Fig. 2(a). In Fig. 2(a),

**t**he region within the heavy green lines has non-zero${\chi}_{eff}^{(2)}(x,z)$and represents regions of the nonlinear crystal. The optical pump pulse propagates with momentum ${k}_{op}$(with an angular spread due to angular dispersion) as indicated by the black arrow in Fig. 2(a).The corresponding optical fluence is indicated by the cyan/light-blue colormap in Fig. 2(a). In lithium niobate, the group refractive index at a wavelength of ~1μm is ~2.25 and the THz refractive index is ~5. To achieve phase-matching, this would then require the optical and THz beams to propagate at an angle of ${\mathrm{cos}}^{-1}(2.25/5)$ = 63° to each other. The generated THz pulse in the simulation then propagates in the

*z*direction at an angle

*γ*~63° to the optical pump pulse and emerges perpendicular to the output face. The momentum of the THz radiation,${k}_{THz}$is delineated by the red arrow in Fig. 2(a). The corresponding THz fluence is represented by the red colormap.

As the optical pump pulse propagates, it generates THz photons at $\Omega $and simultaneously suffers a frequency downshift by the same amount due to the intra-pulse DFG process. With successive generation of THz photons it repeatedly experiences a ‘cascaded’ frequency down-shift which leads to large spectral broadening of the optical pump pulse spectrum as seen in Fig. 2(b). Some upshifting of the optical frequencies is observed, but the THz waves that are generated by DFG are phased relative to the generating optical waves so as to produce further THz generation and optical downshifting. Thus downshifting of the optical spectrum is favored over upshifting. Even if the total depletion of the optical pump energy is only 1%, the drastic spectral reshaping renders undepleted pump approximations inaccurate.

As the optical pump pulse spectrum ${\left|{E}_{op}(\omega ,x,z)\right|}^{2}$ is modified between locations (i)-(iii), the subsequent THz spectrum ${\left|{E}_{THz}(\Omega ,x,z)\right|}^{2}$ is also modified as shown in Fig. 2(c) and vice-versa. As a result, there is significant spatial variation in both optical and THz spectra. The generated THz spectra are broadband, extending from 0 to 1 THz, consistent with earlier experiments and theory.

For an angularly dispersed pulse, each spectral component of the optical pulse at *ω* has a well-defined transverse-momentum *k _{x}* (smaller

*ω’*s have a more negative

*k*value as shown in Fig. 1(a)). Therefore, spectral broadening of the optical pulse also directly leads to re-distribution of optical pulse energy among various transverse momentum values

_{x}*k*as seen in the plot of the quantity ${{\displaystyle {\int}_{0}^{\infty}\left|{E}_{op}(\omega ,{k}_{x},z)\right|}}^{2}d\omega $ in Fig. 2(d). This quantity represents the total energy in a small transverse-momentum interval

_{x}*Δk*or energy density in momentum space. The spread in transverse momentum in Fig. 2(d) is on the order of 10

_{x}^{4}m

^{−1}, which is still much smaller than the optical wave number (~10

^{6}m

^{−1}) which means the beam is still relatively paraxial.

In the presence of this large spectral broadening one would expect a temporal break-up of the pulse due to group velocity dispersion. In addition, due to broadening in transverse momentum, there is also a spatial break-up of pulses. In combination, a rapid spatio-temporal break-up of the optical pulse occurs. Figure 3 shows snapshots of the electric field intensity of the optical pulse${\left|{E}_{op}(t,x,z)\right|}^{2}$at three different time instants *t* in the *z-x* spatial co-ordinate system. Thus an initially clean TPF in Fig. 3(a) suffers a spatio-temporal break-up upon propagation in Fig. 3(c) as it propagates over very short distances on the order of ~2 mm. Due to this spatio-temporal break-up, different parts of the optical pulse arrive at different times and the generated THz no longer builds up coherently. Neglecting the spectral re-shaping of the optical pump spectrum would result in significantly longer interaction lengths as in [28,29], leading to much larger conversion efficiencies than experimentally achieved [22]. Therefore, exclusion of these cascading effects would not result in the drastic pulse break-up observed here. Furthermore, the spatio-temporal break-up of the optical pump pulse has implications on the optimization of conversion efficiency as well as spatial and spectral properties of the THz radiation which are discussed in Section 6.

## 4. Theoretical formulation for complete THz generation system

#### 4.1 Analytic description of optical pump field propagation through tilted-pulse-front setup

As described in Section 2 and depicted in Fig. 1, a TPF setup imparts a number of spatio-temporal variations to the optical pump pulse which influences the properties of the generated THz radiation.

In this section, we show how to account for these effects for an arbitrary TPF setup by employing dispersive ray pulse matrices [31]. Although developed for passive optical elements, our application of this approach to OR using TPF results in a powerful model, closely connected to experiments. An explicit expression for the electric field of the optical pump pulse at the entrance facet of the crystal is obtained which allows the calculation to be performed rapidly. Note that alternate ray-pulse matrix approaches such as [32] are also applicable. Since the beam size of the optical pump used in OR is much larger than the optical wavelength, paraxial approximations of ray-pulse matrix schemes are valid for the optical pump. Each spectral component ${E}_{op}^{in}(\omega ,{x}_{o},{z}_{o})$ of the optical pump pulse at angular frequency *ω*, with input beam centred at a transverse position ${x}_{in}(\omega )$and propagation direction ${x}_{in}^{\text{'}}(\omega )$ emerges with transverse position ${x}_{out}(\omega )$ and propagation direction${x}_{out}^{\text{'}}(\omega )$ from the optical setup, just before entering the nonlinear crystal.

The dependence of the position and propagation direction on frequency accounts for spatial-chirp and angular dispersion. The propagation direction here is not the physical direction but one normalized by the refractive index of the medium [31] which makes it convenient to apply to systems with interfaces of mismatched refractive index. The relationship between${x}_{in}(\omega )$, ${x}_{in}^{\text{'}}(\omega )$ and${x}_{out}(\omega )$,${x}_{out}^{\text{'}}(\omega )$ is given by Eq. (1). Here, $\overline{\overline{M(\omega )}}$ is the overall ray pulse matrix obtained by the product of ray pulse matrices of individual optical components $\overline{\overline{{M}_{i}(\omega )}}$ in reverse order of incidence. The beam positions and propagation directions after the *i ^{th}* optical element are ${x}_{ou{t}_{i}}(\omega )$ and ${x}_{ou{t}_{i}}^{\text{'}}(\omega )$ and are connected to the input parameters as follows.

The ray pulse matrix of the *i ^{th}* optical component $\overline{\overline{{M}_{i}(\omega )}}$is described by a 3x3 matrix for a single transverse spatial dimension given by

Equation (2) shows that the upper 2x2 matrix is the standard ABCD matrix for Gaussian beams. However, in order to account for dispersion, there are two additional terms *E _{i}* and

*F*which correspond to the terms ${\partial {x}_{out}/\partial \omega |}_{i}(\omega -{\omega}_{0})$ and ${\partial {x}_{out}\text{'}/\partial \omega |}_{i}(\omega -{\omega}_{0})$. These terms refer to the shift in output beam position and output beam propagation direction in response to a shift in frequency. Here, we calculate

_{i}*F*upto the fourth order in frequency and accounts for GVD-AD and higher order terms. Note that the last row of $\overline{\overline{{M}_{i}(\omega )}}$ is [0 0 1] as the source frequency does not change. The matrices for common optical elements are presented in Table 1.

_{i}With the knowledge of the input and output beam positions and propagation directions, a Huygen’s integral can be used [31] to calculate the electric field of the emergent optical pump pulse at the crystal entrance facet after it has passed through the TPF setup as shown below.

In Eq. 3(a), ${w}_{in}(\omega ),{w}_{out}(\omega ,{z}_{0})$ are the *e ^{−2}* beam radii of the spectral component at

*ω*at the input and output of the TPF setup respectively. The pre-factor$\sqrt{{w}_{in}(\omega )/{w}_{out}(\omega ,{z}_{0})}$in Eq. 3(a) represents the change in the optical field intensity due to change in beam size. ${E}_{0}(\omega )$represents the initial spectral amplitude of the electric field at angular frequency

*ω*and may include phase terms describing chirp. The first exponential term in Eq. 3(a), represents the transverse variation of the beam along

*x*, where

_{0}*q*is the usual

_{out}*q*parameter associated with Gaussian beams. Notice that each spectral component is centered at a different beam position${x}_{out}(\omega ,{z}_{0})$, which accounts for spatial-chirp. The second exponential term ${e}^{-jk(\omega ){z}_{0}}$is the phase associated with the roughly

*z*propagating paraxial field. The wave number $k(\omega )=\omega n(\omega )/c$includes the material dispersion of the medium at optical frequencies. The third exponential term${e}^{-jk(\omega ){x}_{out}^{\text{'}}(\omega ){x}_{0}}$ represents angular dispersion which directly influences phase matching. Since the optical pump beam is treated paraxially, ${x}_{out}^{\text{'}}(\omega )\ll 1$. Since, the value of ${x}_{out}^{\text{'}}(\omega )$ is mapped to the configuration of the TPF setup (grating angle, imaging distances etc.), Eq. 3(a) can quantify the impact of imaging errors on the properties of THz radiation (spectrum, spatial profile, conversion efficiency). The various phase terms in Eq. 3(a), $\varphi (\omega )={\varphi}_{1}(\omega )+{\varphi}_{2}(\omega )+{\varphi}_{3}(\omega )$ are given by Eqs. 3(b)-3(d). The ${\varphi}_{1}(\omega )$term represents the phase correction due to a Gaussian beam. The remaining terms ${\varphi}_{2}(\omega )$and ${\varphi}_{3}(\omega )$account for additional phases introduced by the TPF setup. Thus using these dispersive ray pulse matrices, one obtains an explicit expression for the electric ${E}_{op}^{out}(\omega ,{x}_{0},{z}_{0})$ field for every spectral component of the optical pump pulse spectrum. This calculated ${E}_{op}^{out}(\omega ,{x}_{0},{z}_{0})$ value is used as an initial condition to solve the nonlinear coupled system of wave equations. This expression accounts for spatial variations, material dispersion, angular dispersion including GVD-AD, spatial frequency-variations and spatial variations in pulsewidth. It allows for analysing the properties of the produced THz radiation as a function of characteristic parameters of the TPF setup.

_{0}#### 4.2 Nonlinear polarization due to optical rectification

In Section 4.1, we obtained the electric field of the optical pump pulse inside the crystal in the co-ordinate system (*z _{0}-x_{0}*) as shown in Fig. 1. In this section, we calculate the nonlinear polarization terms which drive the optical and THz fields. The nonlinear polarization will be calculated in the (

*z-x*) co-ordinate system introduced in Fig. 1(b). The transformation between (z

*) and (*

_{0}-x_{0}*z-x*) co-ordinate systems is easily obtained via Eq. (4).

We use Eq. (4) in Eq. (3), to obtain the electric field of the optical pump pulse in the transformed co-ordinates i.e. ${E}_{op}^{out}(\omega ,x,z)$. The value of ${E}_{op}^{out}(\omega ,x,z)$ at $z=-(L-h)\mathrm{sin}\alpha $ serves as the initial condition for the evolution of wave equations described in Section 4.3.

Equation (5) describes the nonlinear polarization term ${P}_{THz}(\Omega ,x,z)$due to OR which drives the THz electric field ${E}_{THz}(\Omega ,x,z)$at the angular frequency *Ω*.

In Eq. (5), ${E}_{op}(\omega +\Omega ,x,z)$ corresponds to the electric field of the spectral component of the optical pump pulse spectrum at angular frequency *ω* and spatial location (*z,x*). The nonlinear polarization term at each spatial location (*z,x*) in Eq. (5) can be seen to be an aggregate of all possible DFG processes between the spectral component ${E}_{op}(\omega +\Omega ,x,z)$and${E}_{op}(\omega ,x,z)$. In Eq. (5), ${\chi}_{eff}^{(2)}(x,z)$ is the effective second order nonlinear susceptibility for OR at each spatial location and${\epsilon}_{0}$ is the free space permittivity. The spatial dependence of the effective non-linear susceptibility is used to account for the geometry of the nonlinear crystal as was shown in Fig. 1(b).

If one substitutes the expression for the electric field from Eq. 3(a) in Eq. (5), various factors which affect phase-matching effects become evident. For instance, because of the angular dispersion term in Eq. 3(a), a term proportional to $\mathrm{exp}\left\{-j\left[k(\omega +\Omega ){x}_{out}^{\text{'}}(\omega +\Omega )-k(\omega ){x}_{out}^{\text{'}}(\omega )\right]x{}_{0}\right\}$ appears in Eq. (5) which enables us to quantify the effect of phase-mismatch due to imaging errors. This term also accounts for the effects of GVD-AD. Similarly, a term proportional to $\mathrm{exp}\left\{-j\left[\frac{\omega +\Omega}{c}\frac{{\left({x}_{0}-{x}_{out}(\omega +\Omega )\right)}^{2}}{2{q}_{out}(\omega +\Omega )}-\frac{\omega}{c}\frac{{\left({x}_{0}-{x}_{out}(\omega )\right)}^{2}}{2{q}_{out}(\omega )}\right]\right\}$ describes a spatial variation in the magnitude of ${P}_{THz}(\Omega ,x,z)$leading to spatial variation of the generated THz frequency. This term also describes the effects of the finite radius of curvature of the optical pump pulse wavefront which affects phase-matching. In addition, there is a term of the form $\mathrm{exp}\left\{-j\left[\varphi (\omega +\Omega )-\phi (\omega )\right]\right\}$which introduces phase mismatch due to various phase accumulations through the optical setup.

Similarly, each spectral component ${E}_{op}(\omega ,x,z)$ of the optical pump pulse spectrum is driven by a nonlinear polarization term ${P}_{op}(\omega ,x,z)$described in Eq. (6).

The first term in Eq. (6) is the analogue term on the right hand side of Eq. (5). It signifies that an optical photon at angular frequency *ω* is created by an aggregate of DFG processes between optical photons at angular frequency *ω + Ω* and THz photons at angular frequency *Ω*. It represents the red-shift of the optical spectrum depicted in Fig. 2(b). The second term corresponds to an aggregate of sum frequency generation (SFG) processes between optical photons at angular frequency *ω-Ω* and THz photons at angular frequency *Ω*. This term partially contributes to the blue-shift of the optical spectrum that was seen in Fig. 2(b). The third term in Eq. (6) represents the SPM term. Here, ${E}_{op}(t,x,z)$is the time-domain electric field of the optical pump pulse and **F*** _{t}* represents the Fourier transform between time and frequency domains. The intensity dependent refractive index coefficient is given by${n}_{2}(x,z)$. Since the SPM term in Eq. (6) contains details of the spatial distribution of the optical field, it also accounts for self-focusing effects which are not accounted for in the 1-D case [22]. The final term models Stimulated Raman Scattering. This term is related to the SPM term but includes the effects of a Raman gain lineshape given by${h}_{R}(\omega \text{'})$.

#### 4.3 Solving the 2-D coupled non-linear wave equations using Fourier decomposition

In this section, we present our approach for solving the coupled system of nonlinear wave equations. The nonlinear polarization terms defined in Eqs. (5) and (6) will drive the corresponding THz and optical fields. The nonlinear scalar wave equation for the evolution of the THz field ${E}_{THz}(\Omega ,x,z)$ is presented in Eq. (7). A single scalar wave equation, treating one vector component of the THz beam would suffice for lithium niobate since the *d*_{33} element of the second order nonlinear tensor is much larger than other *d _{ij}* and THz generation scales as${d}_{ij}^{2}$.

In Eq. (7), $k(\Omega )=\Omega n(\Omega )/c$is the wave number at the THz angular frequency *Ω* and *n(Ω)* is the corresponding refractive index. Similar to Eq. (7), one can also write the corresponding wave equation for the optical fields at various angular frequencies *ω* in Eq. 8(a).

In the (*z-x*) co-ordinate system, the optical radiation is propagating at a large angle of ~63° relative to the THz radiation. Furthermore, the optical wave number is a hundred times larger than the THz wave number. Therefore, the spread in transverse momentum required to account for the oblique propagation of the optical field would be very large. Direct solutions of Eqs. (7) and 8(a) will therefore be cumbersome because of the fine spatial resolution that will be required to account for this large spread in transverse momentum. Also, solving Eqs. (7) and 8(a) would require knowledge of the initial value of the first derivative of the electric fields in addition to the initial electric field values. In order to provide an efficient solution to the problem, we define a solution of the form${E}_{op}(\omega ,x,z)={A}_{op}(\omega ,x,z){e}^{j{k}_{x0}(\omega ).x}{e}^{-j{k}_{z0}(\omega ).z}$. Here, ${k}_{x0}(\omega )$and ${k}_{z0}(\omega )$ are the momentum components in the *x* and *z* directions and are functions of *ω* due to angular dispersion. This expression signifies that the optical beam has a narrow spread in transverse-momentum in comparison to${k}_{x0}(\omega )$. Equation 8(a) then reduces to Eq. 8(b) for the evolution of the optical field. It only requires the initial value of the electric field and not the initial value of the first derivative of the electric field. The second derivative in *z* is dropped since$\left|{k}_{z0}{A}_{op}^{\text{\'}}(\omega ,z)\right|>>\left|{A}_{op}"(\omega ,z)\right|$. The second derivative in *x* is retained to account for the radius of curvature of the Gaussian beam wavefront.

Equation (7) drives the THz field which in turn affects the optical field via the nonlinear polarization term ${P}_{op}(\omega ,x,z)$defined in Eq. (6). However, the optical field directly influences the THz field via the nonlinear polarization term ${P}_{THz}(\Omega ,x,z)$ defined in Eq. (5). Thus, Eqs. (7) and 8(a) form a coupled system of wave equations for THz and optical fields. An elegant solution to this system can be obtained via spatial Fourier decomposition. This effectively breaks up Eqs. (7) and 8(b) into a system of coupled 1-D first order differential equations which can be solved highly efficiently as will be shown below.

Applying Fourier transforms to both the left and right hand sides of Eq. (7), we obtain Eq. 9(a), which is a 1-D differential equation in *z* for${E}_{THz}(\Omega ,{k}_{x},z)$. The evolution of ${E}_{THz}(\Omega ,{k}_{x},z)$ at each *Ω* and *k _{x}* can be updated in parallel.

Equation 9(a) can be further simplified to Eq. 9(b) by assuming a solution of the form$E(\Omega ,{k}_{x},z)=A(\Omega ,{k}_{x},z){e}^{-j{k}_{z}(\Omega ,{k}_{x})z}$ which is a slowly varying envelope approximation, where${k}_{z}(\Omega ,{k}_{x})=\sqrt{{k}^{2}(\Omega )-{k}_{x}^{2}}$. This is not to be confused with a paraxial approximation for the THz field since each component has a different${k}_{z}(\Omega ,{k}_{x})$, which allows for deviations of $\pm \pi /2$ from the *z* direction. However, back-propagating THz components are ignored as they would not be well phase matched. The THz absorption coefficient at angular frequency *Ω* is given by *α(Ω)*.

Similar to Eq. 9(b), the corresponding 1-D differential equation in *z* for the optical field in is

In Eq. (10), the first term on the right hand side represents the oblique propagation of the optical field. The second term corresponds to the quadratic phase associated with a Gaussian beam and the final term is the spatial Fourier transform of the nonlinear polarization in Eq. (6). Equations 9(b) and (10) form a system of coupled 1-D differential equations for each value of *k _{x}*. Equations (5), (6) and 9(b) and (10) are thus solved progressively up to$z=h\mathrm{sin}\alpha $.

#### 4.4 Transmission and propagation of THz radiation through exit facet of crystal

To calculate the transmitted THz field, we use standard Fresnel reflection coefficients as a function of the transverse-momentum *k _{x}*. For THz field polarization perpendicular to the plane of Fig. 1, the Fresnel reflection coefficients are presented in Eq. (11).

*z*away from the exit surface is calculated without paraxial approximations:

_{d}*k*and

_{x}*ω*and Ω. Numerical integration was performed using a 4th order Runge-Kutta method. The evaluation of Fourier transforms was accelerated by GPU parallelization.

**5.** Validation of the model: comparison to experiments and analytic calculations

In this section, the developed model is validated against analytic theory and experiments [12]. Simulations assume the setup shown in Fig. 1(a). The ABCDEF matrices for the various optical components are presented in Table.1. The *F* term for the diffraction grating contains dispersive terms up to the 4th order in frequency and accounts for GVD-AD and higher order dispersion terms. The full-width at half maximum (FWHM) pulse duration was assumed to be 0.5 ps with a fluence of 20 mJ/cm^{2} and a peak intensity of 40 GW/cm^{2}. The central wavelength of the pulse is 1030 nm. This closely resembles the output of the commercial diode-pumped Yb:KYW chirped pulse amplification system (s-Pulse, Amplitude Systemès) used in [12] with pulse energies upto 2 mJ at 1 kHz repetition rate and 2.6 nm spectral bandwidth. The nonlinear material is assumed to be congruent lithium niobate doped with 5% magnesium oxide. The effective second order susceptibility was assumed to be ${\chi}_{eff}^{(2)}(x,z)$ = 360 pm/V [34] and the intensity dependent refractive index coefficient is *n _{2}* = 10

^{−15}cm

^{2}/W [35]. The crystal temperature was

*T*= 300 K. The optical beam with input

*e*radius of ${w}_{in}=2.5$mm is incident at

^{−2}*h*= 2.2 mm from the apex of the crystal. The apex angle of the prism is assumed to be

*α*= 58°. The focal length of the lens was 23 cm. The refractive indices and Raman gain lineshape for congruent lithium niobate are taken from [36] and [37] respectively. We calculate the conversion efficiency as a function of the grating incidence angle (

*θ*) and the grating to lens (

_{i}*s*) and lens to crystal (

_{1}*s*) distances. The values of

_{2}*θ*and

_{i}, s_{1}*s*obtained for maximum conversion efficiency are compared to analytic calculations [25].

_{2}The simulation results are plotted in Fig. 4(a). A maximum conversion efficiency of 0.8% is obtained at *s _{1}* = 60.89 cm,

*s*= 36.84 cm and

_{2}*θ*46.5°. It can be seen that these values of

_{i}=*s*satisfy the imaging condition for the setup, i.e. ${s}_{1}^{-1}+{s}_{2}^{-1}={f}^{-1}$. The optimal imaging condition is determined by the magnification required to produce the optimal pulse-front-tilt angle inside the crystal. Analytic calculations were developed in [25] to supply optimal imaging conditions and are presented in Fig. 4(b). The blue and red curves in Fig. 4(b) correspond to the values of

_{1}, s_{2}*s*and

_{1}*s*for various pulse-front-tilt angle values

_{2}*γ*and are plotted along

*y*-axis on the left. The grating incidence angle

*θ*as a function of

_{i}*γ*is given by the black curve plotted along the

*y*-axis on the right. We know that for pumping at 1030 nm, the optimum pulse-front-tilt angle $\gamma ={\mathrm{cos}}^{-1}({n}_{g}({\lambda}_{0})n{(\Omega )}^{-1})$~63°. For this value of

*γ*, it is seen that the corresponding imaging conditions closely match the simulation results, i.e.

*s*= 60.06 cm,

_{1}*s*= 37.1 cm and

_{2}*θ*46°. Since the simulation is

_{i}=*ab-initio*, i.e it does not make any prior assumptions about the optimal imaging conditions, this agreement is a strong validation of the developed model. As

*s*or

_{1}*s*are varied, the direction ${x}_{out}^{\text{'}}(\omega )$ at which each optical component at

_{2}*ω*emerges from the setup changes which affects phase-matching via Eq. (5). Thus, the presented formalism can map the performance of the system directly to experimental conditions.

In Fig. 4(a), we see how small deviations in imaging conditions can lead to sizeable degradation of conversion efficiency. For example, *Δs _{2}*~1 mm, leads to drop in conversion efficiency by about 40%. In Fig. 4(c), we show experimental scans of conversion efficiency for similar parameters. The white spaces in the figure indicate regions where data was not collected. For similar displacements

*Δs*and

_{1}*Δs*from the optimum values, the conversion efficiency reduction agrees well with the calculations in Fig. 4(a). The slight difference in the tilt of the ellipse in Fig. 4(c) can be attributed to a different grating incidence angle from that presented in Fig. 4(a).

_{2}Further verification of the model is provided by comparisons to experiments [12]. For the optical pump conditions described for Fig. 4(a), a conversion efficiency of 0.8% (*w _{in}* = 3.5 mm

*h*= 2.2 mm) which is in reasonable agreement with the experimentally reported value of 1.15%. The larger number may be partially owed to uncertainties in THz absorption coefficients below 0.9 THz [36]. When the conversion efficiency is 0.8%, the absorption coefficient below 0.9 THz was ~10 cm

^{−1}. When this was adjusted to 5 cm

^{−1}, the resulting conversion efficiency was 0.9%, which is closer to the experimental result.

In Fig. 5(a), the experimentally obtained THz spectrum is compared to theoretical calculations. The calculation presented is the spatially averaged spectrum $\u3008{\left|{E}_{THz}(\Omega )\right|}^{2}\u3009$$={\displaystyle \underset{-\infty}{\overset{\infty}{\int}}{\left|{E}_{THz}(\Omega ,x,h\mathrm{sin}\alpha )\right|}^{2}dx}$ at the output facet of the crystal. The theoretical and experimental spectra are in agreement and peak at ~0.45 THz, in line with expectations from prior models [22], [25]. In Fig. 5(b), the experimentally reported optical spectrum is compared to theoretical calculations. The black dotted line represents spatial averaging of the optical spectrum over a single transverse dimension, i.e the 2-D average${\u3008{\left|{E}_{op}(\omega )\right|}^{2}\u3009}_{2-D}={\displaystyle \underset{-\infty}{\overset{\infty}{\int}}{\left|{E}_{op}(\omega ,x,h\mathrm{sin}\alpha )\right|}^{2}dx}$. The theoretically obtained optical spectrum shows significant broadening and red-shift similar to the experiments. Some blue-shift is also observed due to the much less dominant SFG process between the optical pump and THz radiation. The theoretically obtained spectrum is larger in extent in relation to the above estimate as well as the experimental spectrum. We simulated several 2-D slices, each with different optical pump intensity, to mimic a 3-D simulation. Based on this quasi 3-D simulation, the optical spectrum was averaged over both transverse spatial dimensions (*x* and *y*) to yield the 3-D average ${\u3008{\left|{E}_{op}(\omega )\right|}^{2}\u3009}_{3-D}={\displaystyle \underset{-\infty}{\overset{\infty}{\int}}{\displaystyle \underset{-\infty}{\overset{\infty}{\int}}{\left|{E}_{op}(\omega ,x,y,h\mathrm{sin}\alpha )\right|}^{2}dx}}dy$. This is shown in the solid black curve and is in better agreement with experiments compared to ${\u3008{\left|{E}_{op}(\omega )\right|}^{2}\u3009}_{2-D}$ . The reason for this is that for lower optical intensities, the extent of spectral broadening will be less as the nonlinear polarization term in Eqs. (5)-(6) would be smaller and therefore the spatial average shows less red-shift and spectral broadening. Another reason for the disparity between experiments and theory could be attributed to uncertainties in the measurement of the optical spectrum. Extremities of the frequency spectrum may not have been collected due to their larger divergence. The internal conversion efficiency${\eta}_{\mathrm{int}}~{(1-R)}^{-1}[{\alpha}_{THz}L{(1-{e}^{-{\alpha}_{THz}L})}^{-1}]{\eta}_{ext}$, where *R* is the power reflection and the ${\alpha}_{THz}L{(1-{e}^{-{\alpha}_{THz}L})}^{-1}$ factor accounts for THz absorption. For *L* = 2mm, *α _{THz}* = 5-10cm

^{-1,}R~0.5, this amounts to about 4-5%. The red-shift in wavelength $\Delta \lambda ={\lambda}_{0}{\eta}_{\mathrm{int}}$ is therefore ~40-50 nm, which translates to a wavelength spread upto 1070-1080 nm, which is in reasonable agreement with the simulated values.

## 6. Optimizing conversion efficiency

#### 6.1 Discussion of effective length in two dimensions

It is useful to understand what the effective propagation length *L _{eff}* in a 2-D geometry is. In general, absorption and dispersion determine the optimal value of

*L*. A longer

_{eff}*L*leads to more absorption and dispersion. A shorter effective length would mean less of both but would also translate to lesser THz generation. Therefore, there must exist an optimum value where the amount of THz generation is sufficiently large but absorption and dispersion are small.

_{eff}In a 2-D non-collinear geometry, two parameters influence the extent of absorption and dispersion. These include the beam radius ${w}_{in}$(or${w}_{out}$) and beam position *h* of the optical pump beam (with respect to the apex of the crystal, see Fig. 1 (b) for definitions). Therefore, the effective length parameter maps to both *h* and${w}_{in}$, i.e.${L}_{eff}\to g(h,w)$where *g* is some function. Since THz is generated only in the region where there is optical fluence, (i) a larger value of *h* would mean that the THz would propagate over a longer absorptive region (i.e. regions where there is no optical fluence, e.g. see Figs. 6(a) vs. 6(b)). Mathematically, the regions without optical fluence correspond to regions where the nonlinear polarization term ${P}_{THz}(\Omega ,{k}_{x},z)=0$ in Eq. 9(b). (ii) For larger *h,* the optical beam would propagate longer distances which would mean it would suffer greater spectral broadening due to cascading effects and SPM. This would in turn make dispersive effects more acute. Thus a larger *h* increases both dispersion and absorption. At the same time, a small *h* also leads to lesser amount of THz being generated. (iii) A smaller value of beam radius ${w}_{in}$ would lead to a larger region without optical fluence and consequently a larger amount of absorption. (iv) A larger value of ${w}_{in}$ would mean that parts of the optical pump beam propagate a longer distance (since different sections of the beam propagate different distances), which would lead to a greater amount of cascading and dispersive effects. Therefore, it is reasonable to expect an optimal value of ${w}_{in}$ and *h* for a given set of optical pump conditions. Increasing the intensity or initial bandwidth of the optical pulse will lead to more rapid spectral broadening (w.r.t length), which would require a readjustment of *h* and${w}_{in}$.

#### 6.2 Implications on energy scaling

These effects are illustrated in Figs. 6(a)-(c). The fluence, bandwidth, pulsewidth, material parameters, setup conditions, crystal geometry and temperature are the same as those used in Fig. (4). The optical fluence ($\underset{0}{\overset{\infty}{\int}}|{E}_{op}(\omega ,x,z){|}^{2}d\omega$) and THz fluence ($\underset{0}{\overset{\infty}{\int}}|{E}_{THz}(\Omega ,x,z){|}^{2}d\Omega$) in Fig. 6 are plotted in the *z-x* co-ordinate system defined in Fig. 1 (b). In Fig. 6(a), a beam with ${w}_{in}$ = 2.5 mm, is incident at *h* = 1.5 mm from the crystal apex. Figure 6(a) uses the same parameters as that used in Fig. 2 (a) but is normalized to a different scale so that it can be appropriately compared to Figs. 6(b), 6(c). It can be seen how the THz and optical beams have good overlap which reduces absorptive effects and results in a relatively high conversion efficiency of 0.7%. In Fig. 6(b), the same beam is displaced further down the crystal. One sees that increased THz absorption, as delineated in Fig. 6(b), causes the conversion efficiency to drop to 0.3%. In Fig. 6 (c), a larger beam size with ${w}_{in}$ = 10 mm is used at *h* = 5 mm. We see that only a small portion of the optical pump beam cross-section produces THz radiation, resulting in a conversion efficiency of 0.5%. This is because, after initial THz generation, subsequent parts of the beam are spectrally broadened due to cascading effects to an extent that prevents further THz generation in the presence of GVD-AD and GVD-MD. This has an important implication for the scaling of these systems to large pump energies. In Fig. 6(d), we plot the maximum conversion efficiencies for various values of ${w}_{in}$while keeping the peak optical pump intensity constant. The top *x*-axis depicts the corresponding optimal values of beam position *h*. Note that the size of the beam at the input crystal face is${w}_{out}~0.6{w}_{in}$. In Fig. 6(d), for${w}_{in}$> 3.5 mm, there is a drastic drop in the maximum achievable conversion efficiency. If cascading effects were not considered, this drop in conversion efficiency would not be observed. For instance, in the undepleted model of [29], saturation rather than a drop in conversion efficiency for large beam sizes is predicted. There is an initial increase in the maximum achievable conversion efficiency for ${w}_{in}$< 3.5mm, because of reduced absorption due to the increase in beam size (See section 6.1 for explanation). In Fig. 6(d), it is seen that optimal values of *h* increase with larger ${w}_{in}$ and that they are present relatively close to the apex of the crystal. The degradation of conversion efficiency can be circumvented by using an elliptical pump beam with its major axis perpendicular to the plane of tilting (i.e out of the plane of the paper).

#### 6.3 Effects of pump intensity

In Fig. 7(a), we present experimental results for the conversion efficiency as a function of fluence for two different cases. In each case, the fluence was varied by attenuating the pulse energy of the pump source. In the first case, shown by the red curve in Fig. 7(a), the fluence was first reduced to a minimum. The conversion efficiency was then locally optimized at this lowest fluence by only adjusting the beam position *h*. Subsequently, the fluence was progressively increased by reducing the attenuation and the corresponding conversion efficiency was recorded. In the second case, shown by the black curve in Fig. 7(a), the fluence was first maximized. Conversion efficiency was optimized at this maximum fluence by re-adjusting the beam position *h*. Subsequently, the fluence was progressively decreased and the corresponding conversion efficiency was recorded. In each case, an imperfectly compressed pulse with duration of 1.39 ps and transform-limited pulsewidth of 0.5ps was used. The central wavelength of the pulse was 1030 nm. The material parameters and crystal geometry were the same as that described in Fig. 4 and the crystal temperature was 300K. The input beam radius${w}_{in}=2.5mm$. The curves show differing saturation behaviour or hysteresis with a cross-over at a fluence of 22 mJ/cm^{2}. In Fig. 7(b), we simulate these experiments using the developed model. The simulated results show qualitative and quantitative agreement with the experiments. For lower fluences, a larger value of *h* is required to optimize the conversion efficiency. This is because, the smaller peak intensity in this case leads to a slower rate of spectral broadening of the optical pump, thereby causing dispersive effects to be ‘delayed’ with respect to length in their appearance. This leads to a longer effective length or larger *h* for optimum efficiency. Note that the optimal values of *h* are larger compared to those depicted in Fig. 6(d) in Section 6.2 because of the smaller peak intensity of the stretched pulse. Thus we see that the conditions for the optimal conversion efficiency are highly sensitive to optical pump conditions. This could be one reason why, various experiments report very different saturation curves.

#### 6.4 Trade-offs of optimizing conversion efficiency

Finally, we highlight some trade-offs of optimizing only conversion efficiency in Figs. 8(a) and 8(b). The material parameters, TPF setup conditions and crystal temperature are the same as in Fig. 4. Here, the THz spectrum as a function of transverse spatial coordinate (*x*) are shown for two different fluences. In Fig. 8(a), the optical fluence is 10 *mJ/cm ^{2}* which results in a conversion efficiency of 0.5%. The THz spectrum is virtually identical across the beam cross-section (

*x*co-ordinate

*)*as seen in Fig. 8(a). In Fig. 8(b), a conversion efficiency of 0.9% is achieved with a fluence of 35

*mJ/cm*. However, the THz spectrum is spatially chirped across the beam-cross section and has an effectively smaller spot-size. As the optical pump beam generates THz, it suffers spectral broadening due to cascading effects. Phase mismatch is accentuated in the presence of this larger spectral bandwidth by GVD-AD and GVD-MD, which causes the coherent growth of THz to cease at the more negative values of

^{2}*x*(See Fig. 2(a): optical beam propagates towards negative

*x*values). In the absence of coherent growth, THz absorption dominates and the THz spectrum red-shifts (since absorption is less for lower THz frequencies). Thus, while conversion efficiency is increased, a spatial chirp is introduced in the THz beam profile along with a reduction in the THz beam spot-size which may not be suitable for certain applications.

*7.* Conclusion and outlook

In conclusion, a new approach to modelling THz generation via optical rectification (OR) using tilted-pulse-fronts (TPF’s) was presented and discussed. The approach was formulated to consider (i) spatio-temporal variations of the optical pump pulse, (ii) coupled non-linear interaction of the THz and optical fields in 2-D as well as (iii) self-phase-modulation and (iv) stimulated Raman scattering. The formulation was done in a way to circumvent challenging numerical issues. It was validated by comparisons to experiments and analytic calculations, with good quantitative agreement in both cases. We described the physics of OR using TPF’s. In particular, we discussed the problem from transverse-momentum and time domains. It was seen that the large spectral broadening which accompanies THz generation also leads to broadening in transverse-momentum since each frequency component in an angularly dispersed beam has a well-defined transverse-momentum value. Thus, the optical pump has an increased frequency and transverse-momentum spread. Group velocity dispersion due to angular and material dispersion then causes a spatio-temporal break-up of the optical pump pulse which inhibits further coherent build-up of THz radiation. It is seen that THz conversion efficiency reduces for very large beam sizes. This suggests the use of optical pump beams that are elliptically shaped for high energy pumping. Guidelines to optimize the setup are provided. Imaging errors were shown to be critical and careful alignment is required to optimize efficiency. It is seen that optimal setup conditions are different for different pump conditions. Finally, we show how optimizing the conversion efficiency could lead to other trade-offs such as a deterioration of the spatial THz beam profile. This work lays the foundation for optimizing such sources for various applications.

## Acknowledgment

This work was supported by AFOSR under grant A9550-12-1-0499, the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 609920, the Center for Free-Electron Laser Science at DESY and the excellence cluster “The Hamburg Centre for Ultrafast Imaging- Structure, Dynamics and Control of Matter at the Atomic Scale” of theDeutsche Forschungsgemeinschaft. We acknowledge support from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 609920. W. R. Huang acknowledges support by a NDSEG graduate fellowship.

## References and Links

**1. **K. A. Leitenstorfer, K. A. Nelson, K. Reimann, and K. Tanaka, “Focus on nonlinear terahertz studies,” New J. Phys. **16**(4), 045016 (2014). [CrossRef]

**2. **T. Kampfrath, K. Tanaka, and K. A. Nelson, “Resonant and nonresonant control over matter and light by intense terahertz transients,” Nat. Photonics **7**(9), 680–690 (2013). [CrossRef]

**3. **H. Y. Hwang, S. Fleischer, N. C. Brandt, B. G. Perkins Jr, M. Liu, K. Fan, A. Sternbach, X. Zhang, R. D. Averitt, and K. A. Nelson, “A review of non-linear terahertz spectroscopy with ultrashort tabletop-laser pulses,” J. Mod. Opt. **62**, 1–33 (2014). [CrossRef]

**4. **S. Fleischer, Y. Zhou, R. W. Field, and K. A. Nelson, “Molecular orientation and alignment by intense single-cycle THz pulses,” Phys. Rev. Lett. **107**(16), 163603 (2011). [CrossRef] [PubMed]

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

**6. **L. Pálfalvi, J. A. Fulop, G. Toth, and J. Hebling, “Evanescent-wave proton postaccelerator driven by intense THz pulse,” Phys. Rev. Spec. Top. **17**(3), 031301 (2014). [CrossRef]

**7. **E. A. Nanni, W. S. Graves, K.-H. Hong, W. R. Huang, K. Ravi, L. J. Wong, G. Moriena, A. Fallahi, D. Miller, and F. X. Kärtner, “Linear electron acceleration in THz waveguides,” Int. Part. Accel. Conf., Dresden (2014).

**8. **W. R. Huang, E. A. Nanni, K. Ravi, K.-H. Hong, L.-J. Wong, P. D. Keathley, A. Fallahi, L. Zapata and F.X.Kärtner, “A terahertz-driven electron gun,” arXiv 1409.8668 (2014).

**9. **L. Wimmer, G. Herink, D. R. Solli, S. V. Yalunin, K. E. Echternkamp, and C. Ropers, “Terahertz control of nanotip photoemission,” Nat. Phys. **10**(6), 432–436 (2014). [CrossRef]

**10. **E. Balogh, K. Kovacs, P. Dombi, J. A. Fulop, G. Farkas, J. Hebling, V. Tosa, and K. Varju, “Single attosecond pulse from terahertz-assisted high-order harmonic generation,” Phys. Rev. A **84**(2), 023806 (2011). [CrossRef]

**11. **O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek, C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W. Koch, and R. Huber, “Sub-cycle control of terahertz high-harmonic generation by dynamical Bloch oscillations,” Nat. Photonics **8**(2), 119–123 (2014). [CrossRef]

**12. **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**(5), 796–798 (2013). [CrossRef] [PubMed]

**13. **C. Vicario, A. V. Ovchinnikov, S. I. Ashitkov, M. B. Agranat, V. E. Fortov, and C. P. Hauri, “Generation of 0.9-mJ THz pulses in DSTMS pumped by a Cr:Mg₂SiO₄ laser,” Opt. Lett. **39**(23), 6632–6635 (2014). [CrossRef] [PubMed]

**14. **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**(21), 1161–1166 (2002). [CrossRef] [PubMed]

**15. **K. L. Yeh, M. C. Hoffmann, J. Hebling, and K. A. Nelson, “Generation of 10 μJ ultrashort terahertz pulses by optical rectification,” Appl. Phys. Lett. **90**(17), 171121 (2007). [CrossRef]

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

**17. **H. Hirori, A. Doi, F. Blanchard, and K. Tanaka, “Single-cycle terahertz pulses with amplitudes exceeding 1 MV/cm generated by optical rectification in LiNbO_{3},” Appl. Phys. Lett. **98**(9), 091106 (2011). [CrossRef]

**18. **J. A. Fülöp, L. Pálfalvi, S. Klingebiel, G. Almási, F. Krausz, S. Karsch, and J. Hebling, “Generation of sub-mJ terahertz pulses by optical rectification,” Opt. Lett. **37**(4), 557–559 (2012). [CrossRef] [PubMed]

**19. **W. R. Huang, S. W. Huang, E. Granados, K. Ravi, K. H. Hong, L. Zapata, and F. X. Kärtner, “Highly efficient terahertz pulse generation by optical rectification in stoichiometric and cryo-cooled congruent lithium niobate,” J. Mod. Opt. **62**, 1–8 (2014). [CrossRef]

**20. **J. A. Fülöp, Z. Ollmann, C. Lombosi, C. Skrobol, S. Klingebiel, L. Pálfalvi, F. Krausz, S. Karsch, and J. Hebling, “Efficient generation of THz pulses with 0.4 mJ energy,” Opt. Express **22**(17), 20155–20163 (2014). [CrossRef] [PubMed]

**21. **S. Akturk, X. Gu, E. Zeek, and R. Trebino, “Pulse-front tilt caused by spatial and temporal chirp,” Opt. Express **12**(19), 4399–4410 (2004). [CrossRef] [PubMed]

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

**23. **X. Wu, S. Carbajo, K. Ravi, F. Ahr, G. Cirmi, Y. Zhou, O. Mucke, and F. Kärtner, “Terahertz Generation in Lithium Niobate Driven by Ti:Sapphire Laser Pulses and its Limitations,” Opt. Lett. **39**(18), 5403–5406 (2014). [CrossRef]

**24. **F. Blanchard, X. Ropagnol, H. Hafez, H. Razavipour, M. Bolduc, R. Morandotti, T. Ozaki, and D. G. Cooke, “Effect of extreme pump pulse reshaping on intense terahertz emission in lithium niobate at multimilliJoule pump energies,” Opt. Lett. **39**(15), 4333–4336 (2014). [CrossRef] [PubMed]

**25. **J. A. Fülöp, L. Pálfalvi, G. Almási, and J. Hebling, “Design of high-energy terahertz sources based on optical rectification,” Opt. Express **18**(12), 12311–12327 (2010). [CrossRef] [PubMed]

**26. **J. A. Fülöp, L. Pálfalvi, M. C. Hoffmann, and J. Hebling, “Towards generation of mJ-level ultrashort THz pulses by optical rectification,” Opt. Express **19**(16), 15090–15097 (2011). [CrossRef] [PubMed]

**27. **S. Bodrov, A. Murzanev, Y. A. Sergeev, Y. A. Malkov, and A. N. Stepanov, “Terahertz generation by tilted-front laser pulses in weakly and strongly nonlinear regimes,” Appl. Phys. Lett. **103**(25), 251103 (2013). [CrossRef]

**28. **M. I. Bakunov, S. B. Bodrov, and E. A. Mashkovich, “Terahertz generation with titled-front laser pulses: dynamic theory for low absorbing crystals,” J. Opt. Soc. Am. B **28**(7), 1724–1734 (2011). [CrossRef]

**29. **M. I. Bakunov and S. B. Bodrov, “Terahertz generation with tilted-front laser pulses in a contact-grating scheme,” J. Opt. Soc. Am. B **31**(11), 2549–2557 (2014). [CrossRef]

**30. **M. Jewariya, M. Nagai, and K. Tanaka, “Enhancement of terahertz wave generation by cascaded chi_2 processes,” J. Opt. Soc. Am. B **26**(9), A101–A106 (2009). [CrossRef]

**31. **O. E. Martinez, “Matrix Formalism for Pulse Compressors,” IEEE J. Quantum Electron. **24**(12), 2530–2536 (1988). [CrossRef]

**32. **A. G. Kostenbauder, “Ray-pulse matrices: a rational treatment for dispersive optical systems,” IEEE J. Quantum Electron. **26**(6), 1148–1157 (1990). [CrossRef]

**33. **M. I. Bakunov, A. V. Maslov, and S. B. Bodrov, “Fresnel Formulas for the forced electromagnetic pulses and their application for optical-to-terahertz conversion in nonlinear crystals,” Phys. Rev. Lett. **99**(20), 203904 (2007). [CrossRef] [PubMed]

**34. **T. Fujiwara, M. Takahashi, M. Ohama, A. J. Ikushima, Y. Furukawa, and K. Kitamura, “Comparison of electro-optic effect between stoichiometric and congruent LiNbO_{3},” Electron. Lett. **35**(6), 499–501 (1999). [CrossRef]

**35. **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/sub 2/ in wide bandgap solids,” IEEE J. Quantum Electron. **32**, 1324–1333 (1996). [CrossRef]

**36. **L. Pálfalvi, J. Hebling, J. Kuhl, A. Peter, and K. Polgar, “Temperature dependence of the absorption and refraction of Mg-doped congruent and stoichiometric LiNbO_{3} in the THz range,” J. Appl. Phys. **97**(12), 123505 (2005). [CrossRef]

**37. **C. R. Phillips, C. Langrock, J. S. Pelc, M. M. Fejer, I. Hartl, and M. E. Fermann, “Supercontinuum generation in quasi-phasematched waveguides,” Opt. Express **19**(20), 18754–18773 (2011). [CrossRef] [PubMed]