Theory of THz generation by Optical Rectification using Tilted-Pulse-Fronts

A model for THz generation by optical rectification using tilted-pulse-fronts is developed. It simultaneously accounts for (i) the spatio-temporal distortions of the optical pump pulse, (ii) the nonlinear coupled interaction of THz and optical radiation in two spatial dimensions (2-D), (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 (kx) domain as a consequence of the spectral broadening caused by 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. It is seen that optimization of the setup is highly dependent on optical pump conditions. Trade-offs of optimizing the optical-to-THz conversion efficiency on the spatial and spectral properties of THz radiation is discussed to guide the development of such sources.

Of various high field THz generation modalities, optical rectification (OR) of femtosecond laser pulses with tilted-pulse-fronts in lithium niobate has emerged as the most efficient THz generation technique. Among various efforts, it was developed in [9][10][11][12][13] 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 phasematching is achieved. OR using TPF's has resulted in optical-to-THz conversion efficiencies (henceforth referred to as conversion efficiency) in excess of 1% [14][15] and the highest THz pulse energy of 0.4 mJ [16] to date. Therefore, the approach is promising for the development of laboratory scale THz sources with pulse energies greater than 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's, 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 termed spatio-temporal distortions [17] 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 which exceed the Manley-Rowe limit. On the other hand, in the presence of group velocity dispersion due to angular dispersion (GVD-AD) and material dispersion (GVD-MD), this spectral broadening inhibits further THz generation [18][19][20]. 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 distortions imparted by the TPF setup would also have to be considered.
Previously presented models broadly comprise 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) and (ii) 1-D spatial models which account for cascading effects. Of works in category (i), a one-dimensional (1-D) spatial model including the effects of material dispersion and GVD-AD was presented in [21][22]. In [23], a 1-D model considering material dispersion, GVD-AD and self-phase modulation (SPM) was presented. In [24], a 2-D model which took into account material dispersion, GVD-AD and crystal geometry was developed. In category (ii), an effective 1-D model with cascading effects was first presented in [25] and improved in [18].
In this paper, we present the formulation of a 2-D model which simultaneously accounts for spatio-temporal distortions of the optical pump pulse, cascading effects, SPM and stimulated Raman scattering (SRS), material dispersion, THz absorption as well as geometry of the nonlinear crystal. The developed model is applicable to the simulation of a variety of OR systems with different TPF setups, crystal geometries and optical pump pulse formats.
In Section 2, we outline our general approach. In Section 3, we describe the physics of THz generation using OR with TPF's. In particular, a discussion from transverse momentum (k x ) and time domain viewpoints is introduced. It is seen that the generation of THz results in the broadening of the optical pump pulse in both frequency and transversemomentum domains. In the presence of this increased frequency and transversemomentum 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, 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 the experimentally measured 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 an overview for constructing sources customized optimally for various applications. to generate a tilted-pulse-front. The model accounts for the angular dispersion of various spectral components which can generate THz radiation inside the nonlinear crystal by satisfying the appropriate phase-matching condition for optical rectification. From a time-domain viewpoint, the angularly dispersed pulse forms a tiltedpulse-front shown by the red ellipses. THz radiation is generated perpendicular to this tilted-pulse-front (red arrow). (b) 2-D computational space for solving coupled nonlinear wave equations for optical rectification. Nonlinear crystal geometry is accounted for by delineating an appropriate distribution of χeff (2) (x,z). Edges of the distribution are smoothed out to avoid discontinuities. The refractive index is homogeneously distributed throughout the computational space. The optical beam is centered at a distance h from the apex of the crystal which sets the limits to the computational region. The optical field at the beginning of the lattice is calculated analytically using dispersive ray pulse matrices. The THz field profile can be calculated at a distance zd from the crystal after Fresnel reflection is taken into account.

THz Output Field
The overall schematic of our approach is depicted in Fig. 1. An optical pump pulse with input electric field described by the complex variable 0 0 ( , , ) in op E x z ω at angular frequency ω, propagating in the z 0 direction is incident on a TPF setup. In Fig. 1, 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. 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 infinitely many frequency components.
An optical pulse with a TPF is typically angularly dispersed [17]. Therefore, various frequency components of the emergent optical field described by 0 0 ( , , ) out op E x z ω 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) ( ), ( ) k k ω ω+ Ω   respectively. In the time domain, such an angularly dispersed pulse has an intensity profile tilted with respect to its propagation direction as shown 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 ( ) ( ) ( ) k k k ω ω Ω = + Ω −    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.
Since the various frequency components of the optical pump pulse are angularly separated, there is a spatial variation in the average frequency (spatial-chirp) as well as amplitude across the optical 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 [26] in Section 4.1. This approach efficiently models the various spatio-temporal distortions associated with the optical pulse for an arbitrary TPF setup and supplies 0 0 ( , , 0) out op E x z ω = at the input face of the nonlinear crystal as shown in Fig. 1(b).
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). 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 coordinate system then has two key advantages. Firstly, in this set of axes, the THz radiation has small transverse-momentum components, i.e. k x~0 , which relaxes the constraints on spatial resolution Δx and consequently alleviates computational cost. 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 ω is considered. Only regions of the computational space physically occupied by the crystal would have a non-zero value of second order susceptibility (2) ( , ) eff 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 ( ) α as shown in Fig. 1(b). The initial optical field profile along the line 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.

Optical Fluence
THz Fluence Fig. 2(a). Spatial distribution of the optical and THz fluences: The THz field propagates in the z direction with the optical field at an angle γ~63° with respect to it. THz is only generated over a small portion of the optical field due to spectral broadening of the optical pulse by cascading which results in disruption of phase-matching due to group velocity dispersion (by angular and material dispersion) for subsequent portions of the beam. (b) Optical spectrum is broadened between (i)-(iii) due to cascading effects. (c) THz spectra at locations (i)-(iii) show significant spatial-chirp due to spatial variations of the optical electric field in (c). (d) Since each frequency component has a certain value of transverse-momentum in an angularly dispersed beam, spectral broadening also necessarily results in broadening in transverse-momentum kx. As the optical spectrum broadens, there is a broadening in transverse-momentum between z = -0.3 mm and 1.4 mm.
In Fig. 2, we provide a sample solution using the developed model. This will bring out the essential physics and put the formulation in Section 4 into context. Figure 2(a), depicts the optical and THz fluences as a function of space. The region within the green lines has non-zero (2) The optical pump pulse propagates obliquely in the nonlinear crystal as indicated by the cyan/light-blue colormap in Fig. 2(a). The THz pulse however propagates in the z direction, at an angle γ~63° to the optical pump pulse and emerges perpendicular to the output face as shown by the THz fluence in the red colormap.
As the optical pump pulse propagates, it generates THz photons at Ω and simultaneously suffers a frequency downshift by the same amount. With successive generation of THz photons it repeatedly experiences a frequency down-shift or 'cascaded' frequency down-shift which leads to large spectral broadening of the optical pump pulse spectrum as seen in Fig. 2(b).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 is modified between locations (i)-(iii), the subsequent THz spectrum 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 x 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 k x as seen in Fig. 2(d). The spread in transverse momentum in Fig. 2(d) is on the order of 10 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 breakup 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 as shown in Fig. 3. Thus an initially clean TPF in Fig. 2(a) suffers a spatio-temporal break-up upon propagation in Fig. 2(c) as it propagates over very short distances on the order of ~ 2 mm. Due to this spatiotemporal break-up, different parts of the optical pulse arrive at different times and the generated THz no longer builds up coherently.

Propagation of Optical Pump Fields through the Optical Setup
As described in Section 2 and depicted in Fig. 1, a TPF setup imparts a number of spatio-temporal distortions 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 [26]. Although developed for passive optical elements, our application of this approach to OR using TPF's results in a powerful model, closely connected to experiments. An explicit expression for the electric field of the optical pump pulse is obtained which allows calculations to be performed rapidly. Note that alternate ray-pulse matrix approaches such as [27] 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 The ray pulse matrix of the i th optical component ( ) i M ω is described by a 3x3 matrix for a single transverse spatial dimension as shown in Eq. (2).
( ) Equation (2) shows that the upper 2x2 matrix is nothing but the standard ABCD matrix for Gaussian beams. However, in order to account for dispersion, there are two additional terms E i and F i which correspond to the partial derivatives 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 F i upto the fourth order in frequency and accounts for GVD-AD and higher order terms. Note that the last row of ( ) i M ω is [0 0 1] as the source frequency does not change. With the knowledge of the input and output beam positions and propagation directions, a Huygen's integral can be used [26] to calculate the electric field of the emergent optical pump pulse after it has passed through the TPF setup as shown in Eq.(3). x z ω 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 a one-to-one correspondence between the TPF setup configuration and the properties of THz radiation to be established.

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 0 -x 0 ) and (z-x) co-ordinate systems is easily obtained via Eq.  x z ω . In Eq. (5), (2) ( , ) eff x z χ is the effective second order nonlinear susceptibility for OR at each spatial location and 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) 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 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, ( , , ) op E 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 2 ( , ) n 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. 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 ( ') R h ω .

Solving the 2-D non-linear wave equation 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) In Eq. (7), ( ) ( ) / k n 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).
Thus using Eqs. (1)-(12), we can model an arbitrary TPF setup for THz generation by OR. Equations 9 and 10 are effectively a 1-D system of coupled equations and can be solved in parallel for various k x and ω and Ω. Numerical integration was performed using a 4 th order Runge-Kutta method. The evaluation of Fourier transforms was accelerated by GPU parallelization.

Optical Element ABCDEF Matrix
Diffraction grating cos cos In this section, the developed model is validated against analytic theory and experiments.
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 upto the 4 th 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 effective second order susceptibility was assumed to be (2) ( , ) eff x z χ = 360 pm/V [29] and the intensity dependent refractive index coefficient was n 2 =10 -15 cm 2 /W [30]. The optical beam with input e -2 radius of 2.5 in w = mm and is incident at h = 1.5 mm from the apex of the crystal. The focal length of the lens was 23 cm. The refractive indices and Raman gain lineshape are taken from [31] and [32] respectively. We calculate the conversion efficiency as a function of the grating incidence angle (θ i ), grating to lens (s 1 ) and lens to crystal (s 2 ) distances. The values of θ i , s 1 and s 2 obtained for maximum conversion efficiency are compared to analytic calculations [21]. Fig.4(a): Simulation of conversion efficiency as a function of imaging conditions. The surface plot shows the conversion efficiency versus displacements from optimal imaging distances Δs1 and Δs2 . s1, s2 are the lens-tograting and lens-to-crystal distances respectively. The inset shows conversion efficiency versus incidence angle to diffraction grating. As s1,s2 are varied, there is variation of the pulse-front-tilt angle which leads to a change in conversion efficiency. Careful optimization of the experimental setup is required to identify the optimal conversion efficiency point. (b) Theoretical calculations of optimal imaging conditions for various pulse-fronttilt angles based on analytic theory from [21]. For a pulse-front-tilt angle of 63°, the imaging conditions are in close agreement with the simulation results, validating the accuracy of the presented model. (c) Experimental scans of conversion efficiency vs displacements Δs1 and Δs2 agree well with the simulations in Fig. 4(a).
The simulation results are plotted in Fig. 4(a). A maximum conversion efficiency is obtained at s 1 = 60.89 cm, s 2 =36.84 cm and θ i =46.5°. It can be seen that these values of s 1 , s 2 satisfy the imaging condition for the setup, i.e. 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 [21] 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 s 1 and s 2 for various pulse-front-tilt angle values γ and are plotted along y-axis on the left. The grating incidence angle θ i as a function of γ 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 x ω at which each optical component at ω 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 1 and Δs 2 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).
Further verification of the model is provided by comparisons to experiments [14]. 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 [31]. 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. Fig. 5 (a) The experimental and theoretically calculated THz spectra are presented. The theoretical calculation is spatially averaged over x and is centred at 0.45 THz, in close agreement with experiments [14] (b) The experimental output optical spectrum (red) is presented along with calculations. The theoretical calculation averaged over a single transverse spatial dimension x (black-dotted) are broadened significantly more than the experiments. However, if the spatial averaging is performed over both transverse spatial dimensions x and y by simulating numerous 2-D slices, the output spectrum matches experiments more closely. The disparity may also be partially explained by the possibility of incomplete collection of extreme optical frequency components with large divergence.
In Fig. 5(a), the experimentally obtained THz spectrum is compared to theoretical calculations. The calculation presented is the spatially averaged spectrum 2 ( ) = Ω ∫ 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 [18], [21]. In Fig. 5(b), the experimentally reported optical spectrum is compared to theoretical calculations.
. This is shown in the solid black curve and is in better agreement with experiments. 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.

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 eff . A longer L eff 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. In a 2-D non-collinear geometry, two parameters influence the extent of absorption and dispersion. These include the beam radius in w (or out w ) 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 in w , i.e. 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 in w 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 in w .

Implications on energy scaling
These effects are illustrated in Figs. 6(a)-(c). The fluence, bandwidth, pulsewidth and material parameters are the same as that used for Figs. (4). In Fig. 6(a), a beam with in w = 2.5 mm, is incident at h = 1.5 mm from the crystal apex. 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 in w = 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 on the scaling of these systems to large pump energies. In Fig. 6(d), we plot the maximum conversion efficiencies for various values of in w while keeping the peak optical pump intensity constant. The top x-axis depicts the corresponding optimal values of h. Note that the size of the beam at the input crystal face is~0.6 out in w w . In Fig. 6(d), for in w > 3.5 mm, there is a drastic drop in the maximum achievable conversion efficiency. There is an initial increase in the maximum achievable conversion efficiency for in w < 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 in w 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).

Effects of pump intensity
In Fig. 7(a), we experimentally determine the conversion efficiency as a function of fluence for two different cases. In these experiments, a pulse with duration of 0.5ps was stretched to 1.39 ps. In the black curve, the setup is optimized to yield maximum conversion efficiency at the highest peak intensity by adjusting h. The fluence is then progressively decreased. In the red curve, the conversion efficiency is optimized for the lowest fluence and then progressively increased. The curves show a 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' 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. Fig. 7(a): Experimentally obtained conversion efficiency saturation curves optimized for different pump intensities. The black curve is optimized for the maximum intensity while the red curve is optimized for the smallest intensity. The optimal experimental conditions are different for different intensities as seen in the hysteresis of the curve (b) Theoretical calculations of conversion efficiency saturation curves for experimental parameters in (a). Good quantitative agreement between experiments and theory is seen. When the intensity is lower, the optimal efficiency occurs at a larger value of h. This is because cascading effects occur at a slower rate and enable a longer effective interaction length. The optimal values of h are larger than that in Fig. 6(d) due to the use of a stretched pulse in the experiments.  Fig. 8 (a). When the fluence is 10 mJ/cm 2 , the conversion efficiency is 0.5%. The THz spectrum as a function of transverse co-ordinate x, is relatively uniform with all points having a broadband THz spectrum centred at ~0.45 THz. (b) As the fluence is increased to 35 mJ/cm 2 , the conversion efficiency increases to 0.9% but the THz beam now contains a large spatial chirp and has an effectively reduced spot size. As the optical beam propagates to more negative values of x, it has been significantly broadened spectrally. Along with dispersive effects, this inhibits further coherent growth of THz radiation. Absorptive effects then dominate, leaving only the lower frequency THz components with smaller absorption intact.

Trade-offs of optimizing conversion efficiency
Finally, we highlight some trade-offs of optimizing only conversion efficiency in Figs. 8(a) and 8(b). Here, the THz spectrum as a function of transverse spatial co-ordinate (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 2 . However, the THz spectrum is spatially chirped across the beam-cross section and has a 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 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.

Conclusion and Future 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 distortions of the optical pump pulse, (ii) coupled non-linear interaction of the THz and optical fields in 2-D as well as (iii) selfphase-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 transversemomentum since each frequency component in an angular dispersed beam has a welldefined transverse-momentum value. Thus, the optical pump as 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 provides an overview for optimizing such sources for various applications.