Analysis of terahertz generation using tilted-pulse-fronts

A 2-D spatio-temporal analysis of terahertz generation by optical rectification of tilted-pulse-fronts is presented. Closed form expressions of terahertz transients and spectra in two spatial dimensions are furnished in the undepleted limit. Importantly, the analysis incorporates spatio-temporal distortions of the optical pump pulse such as angular dispersion, group-velocity dispersion due to angular dispersion, spatial and temporal chirp as well as beam curvature. The importance of the radius of curvature to the tilt-angle and group-velocity dispersion due to angular dispersion to terahertz frequency, conversion efficiency and peak field is revealed.In particular, the deterioration of terahertz frequency, efficiency and field at large pump bandwidths and beam sizes is analytically shown.


I. INTRODUCTION
The generation of terahertz radiation, i.e. electromagnetic radiation in the frequency range spanning 0.1 − 10 THz has experienced a recent surge in interest. The large peak fields that can be generated at these frequencies in combination with long wavelengths offers unique opportunities to manipulate the motion of charged particles. As a result, a number of applications such as compact charged particle acceleration [1][2][3] and streaking [4][5][6], control of emission from nanotips [7] and higherharmonic generation in solids [8] have emerged. In addition, the proximity of these frequencies to lattice vibrations makes them uniquely amenable to probing a number of fundamental scientific phenomena [9].
Among various methods of terahertz generation [10,11], laser-driven nonlinear optical methods have gained ground owing to improvements in solid-state laser technology [12] as well as the intrinsic synchronization they offer-which is valuable to scientific investigations. Furthermore, the long assumed efficiency limitation of this approach due to the large disparity between optical and terahertz photons has been dispelled by experimental demonstrations of cascaded difference frequency generation/optical rectification [13]. Specifically, this class of approaches enables the repeated energy down conversion of optical pump photons to yield energy conversion efficiencies at the η = 1% level [14].
Of various nonlinear optical approaches to generate terahertz radiation, optical rectification of angularly dispersed beams or tilted-pulse-fronts [15] in lithium niobate has become ubiquitous. The approach is accessible by commercially available Titanium Sapphire lasers at 800 nm and has yielded very high conversion efficiencies and pulse energies to date [16].
However, the complex spatio-temporal shaping of the pump pulse introduces many subtleties to the physics governing the approach. The initial treatment of the * koust@mit.edu problem by 1-D spatial models [16][17][18] with effective parameters, while informative, do not account for the intrinsic non-collinearity of the problem.Numerical models considering multiple spatial dimensions [19] as well as the effects of pump depletion [20] furnish accurate quantitative predictions but do not provide the intuition and understanding that arises from analytic approaches. Initial analytic or semi-analytic 2-D models were previously reported [21][22][23][24].
Here, we further develop analytic methods which accurately account for the various spatio-temporal distortions that accompany pulses with tilt-pulse-fronts. These include the effects of group velocity dispersion due to angular dispersion (GVD-AD), spatial and temporal chirp as well as effects of beam curvature. These models in the undepleted limit straddle the middle ground between 1-D models and 2-D numerical models incorporating pump depletion. The analysis sheds light on various spatio-temporal aspects of terahertz generation by tilted-pulse-fronts. Closed form expressions of terahertz spectra in transverse momentum, spatial and temporal domains are furnished.
The analysis reveals the importance of group-velocity dispersion due to angular dispersion. In particular, it reduces the terahertz frequency, conversion efficiency and peak field at large bandwidths and beam sizes. This has ramifications on energy scaling and use of short pump pulses. While this phenomenon has been partially revealed through experiments and numerical studies, here this is shown via analytic approaches and closed-form expressions. The presented expressions are quantitatively useful for room temperature terahertz generation in lithium niobate while they provide qualitative insights at cryogenic temperatures. In particular,this undepleted analysis underscores the importance of considering cascading effects or depletion effects in the low absorption limit.
In Section II, we present the theoretical formulation. In Section III, detailed calculations of terahertz spectra, conversion efficiency and peak field are furnished. We conclude in Section IV. Appendices are provided for selfconsistency. A tilted-pulse-front setup comprising of a diffraction grating and an imaging system. The angularly dispersed pulse produces a tilted-pulse-front (red ellipse) with tilt-angle γ, resulting in terahertz radiation propagating at an angle γ with respect to the direction of pump pulse propagation.

A. General approach
In the most general case, the generation of terahertz radiation by tilted-pulse-fronts follows from a solution of the coupled nonlinear wave equations for terahertz and optical waves [20]. However, since the approach outlined in this paper is analytic, we restrict ourselves to an undepleted solution. A typical tilted-pulse-front setup is depicted in Fig.1. It consists of a diffraction grating off which the input pump pulse is scattered and then imaged into a crystal. As shown in Fig.1, various wave vectors k(ω) of optical components at corresponding angular frequencies ω are angularly dispersed in the x − z plane, producing an intensity front which is tilted with respect to the propagation direction by an angle γ.
The angular dispersion imparted to the pump pulse produces spatio-temporal coupling effects only in one plane, hence making a two-dimensional (x, z) spatial model sufficient to capture the essential physics of the system. Various x − z slices in the y direction, are replicas weighted by the intensity profile of the pump pulse, which can be accounted for by appropriate scaling factors.
The scalar wave equation governing terahertz spectral components at angular frequencies Ω ≥ 0 delineated by the phasor E T Hz (Ω, x, z) is given by Eq.1a (For additional details on transformation between phasor quantities and real fields, the reader is referred to the Appendix). The first term on the right hand side (RHS) of Eq.1a corresponds to terahertz absorption while the second term on the RHS of Eq.1a delineates the nonlinear polarization term P T Hz driving the generation of terahertz radiation at angular frequency Ω. Essentially, P T Hz is an ensemble of difference frequency generation processes proportional to the second order nonlinear susceptibility χ (2) between various spectral components of the optical pump E op (ω, x, z) and is given by Eq.1b.
Performing a Fourier decomposition on Eq.1a, we obtain the ordinary differential equation with respect to z for E T Hz (ω, k x , z) in Eq.2. Here, k x corresponds to the transverse momentum in the x direction.
To solve the above, we set E T Hz (Ω, k x , z) = A T Hz (Ω, k x , z)e −jkz(Ω,kx)z , where x is the z-component of the terahertz wave vector k(Ω).In using the above ansatz , a 1/cosγ pre-factor is obtained for the absorption term due to the pre-factor k(Ω)/k z (Ω, k x ) ≈ 1/cosγ. Here, γ is the angle at which the generated terahertz wave propagates with respect to the pump beam and is also equal to the tilt-angle. Furthermore, the initial condition A(Ω, k x , z = 0) = 0 is assumed to delineate the absence of any terahertz field at the entrance boundary of the crystal. The general solution for E T Hz (Ω, k x , z) is subsequently obtained as follows : In Eq.3b, ∆k(Ω, k x ) represents phase-mismatch between the pump and terahertz fields. Notice that setting Eq.3b to zero,yields the well-known relationship for the tilt angle γ = cos −1 (n g /n T Hz ), where n g and n T Hz are the optical group and terahertz phase refractive indices respectively.

B. Nonlinear polarization
Having laid out the general framework Eq.(1a), we proceed to obtain solutions for the generated terahertz field for an optical pump pulse with spatio-temporal distortions. We proceed by defining the optical pump spectrum E op (ω, x, z) in Eq.4. Due to the large number of variables, we provide a glossary in Table.III in the appendix for quick reference.
delineates the spectral components of the optical field at angular frequencies ∆ω = ω − ω 0 (where ω 0 is the center frequency of the pump) with various spatio-temporal distortions. The first exponential term describes a spectrum with transform limited (TL) pulse duration τ 0 . The second exponential term on the RHS of Eq.4, represents a spatially chirped transverse beam with e −2 radius w 0 . Here ζ(z) = dx c (z)/dω and delineates that different spectral components are centered at different spatial locations x c (z). The z dependence of ζ, x c arises from angular dispersion in the pump beam, which causes their spectral positions to change. The third and fourth exponential terms represent temporal chirp with group delay dispersion (GDD) φ 0 that may be imparted to the beam and phases accrued due to the finite radius of curvature R 0 (z) of the beam respectively. The fifth exponential term represents the z-directed momentum of the optical field. As shown, this is represented by a polyno-mial expression about the wave number k 0 = ω 0 n(ω 0 )/c at the central angular frequency ω 0 , accounting for the group velocity v g and material dispersion via the parameter k m . Here n(ω 0 ) represents the refractive index of material at the central angular frequency ω 0 . The penultimate and final terms are the most important terms for a beam forming a tilted-pulse-front (TPF). These correspond to the angular dispersion term β 0 = dθ/dω and the group-velocity dispersion due to angular dispersion (GVD-AD) term denoted by k T . The latter represents the fact that the angles of various spectral components are not distributed equally for equal increments in frequency ∆ω.
Using, the above expression for the optical spectrum E op (ω, x, z) in the nonlinear polarization term defined by Eq.1b, one obtains the following expression for P T Hz (Ω, x, z) in Eq.5a.
Equation 5a exhibits a number of spatio-temporal coupling effects. The first exponential term in Eq.5a shows how the terahertz pulse duration is given by τ / √ 2 fork T , k m , φ 0 = 0. The second exponential term depicts an increase in terahertz pulse duration due to the temporal chirp φ 0 as well as the GVD-AD (k T ) and GVD-MD terms (k m ). Specifically, while GVD-AD causes an increase in the pulse duration in the transverse (x) direction, the duration in the propagation direction z increases due to material dispersion.
A further illustration of spatio-temporal coupling is evident by the change in effective pulse duration from τ 0 to τ in Eq.5b. Here, the spatial-chirp ζ is seen to increase the effective duration. In a spatially chirped beam, all spectral components do not overlap in space. therefore, only a certain fraction of the original bandwidth is overlapped to produce terahertz radiation and this manifests via an increased effective pulse duration.
The third exponential term simply shows that the beam size of the generated terahertz radiation would correspond to a beam radius w 0 / √ 2, as with any second order nonlinear process.
The fourth exponential term shows that the polarization contains a radius of curvature which is different from that of the incident pulse. The schematic in Fig.2, depicts why this is the case. In the presence of a finite k T , the phase-matched terahertz directions are not all the same and thus produce a phase-front with some curvature. Since, the extent of this spread in terahertz directions shall be proportional to the incident pump bandwidth, the situation is more adverse for shorter pump pulse durations. This is yet another example of a spatiotemporal coupling effect that arises in TPFs that can only be accounted for via an appropriate modelling of the optical pump field. Therefore, the above equations already indicate the great importance of the k T term in determining the properties of the generated terahertz radiation. While prior work has illustrated the importance of GVD-AD numerically [25], here we are able to provide an insight into its effects analytically.

C. Tilt angle
A critical aspect of Eq.5a is delinated in the two final exponential phase terms, which represent a line of constant phase z + tanγx = 0. Here, γ is the pulse-front-tilt angle and is given by Eq.6.
Each term in Eq.6 represents a different source of pulse-front-tilt. The first term in Eq.6 is tilt that is obtained from angular dispersion and is the one most widely used in the context of terahertz generation [26]. The second term in Eq.6 describes pulse-front-tilt due to simultaneous spatial and temporal chirp. It indicates that if different colors are located at different transverse locations and each color arrives at a different time, then one obtains a pulse-front-tilt. This term has been described by prior work [27], but has not been examined in the context of terahertz generation.
The third term in Eq.6 is pulse-front-tilt introduced due to a finite radius of curvature and spatial chirp. While, prior work has suggested the relevance of the radius of curvature to pulse-front-tilt [28], here we provide an explicit expression. In the results and discussion section and Fig.8, we specifically highlight the relevance of this term in conventional tilted-pulse-front setups.

D. Polarization in the transverse momentum domain
To solve Eq.3a, we first need to obtain an expression for P T Hz (Ω, k x , z), which is obtained by performing a Fourier transform of Eq.5a. This yields the following expression and associated effective parameters.
In Eq.7a, the first exponential term is identical to the first exponential term in Eq.5a and has already been discussed. The second term in Eq.7a is due to the group delay dispersion term φ 0 +k m z. From the third exponential term in Eq.7a, we see that the polarization is distributed about k x = k(Ω)sinγ. This indicates that the polarization term drives the radiation of terahertz waves at an angle γ with respect to the pump direction, as expect.
Various effective parameters in Eq.7a are delineated in Eqs.7b to 7f. Barring the net group delay dispersion term φ in Eq.7b, every term illustrates the importance of the role played by spatial-chirp and GVD-AD (k T ) terms in modifying the effective radius of curvature R (Eq.7c), beam radius w (Eq.7d) and beam position x 0 (Eq.7f). Furthermore, the modifications are more adverse for short durations τ or equivalently, large bandwidths. This point will be revisited repeatedly throughout the remainder of this paper.
E. Terahertz spectra 1. Closed form expressions for ET Hz (Ω, kx, z) While Eq.3a is valid in general , it requires numerical evaluation. However, it may be reduced to a closed-form expression in frequency, transverse momentum and longitudinal space, i.e.(Ω, k x , z) as shown in Eq.8a(See Appendix for derivation).
Upon obtaining Eq.8a, one may then obtain the spatial profile by taking an inverse spatial Fourier transform as shown in Eq.8b. Subsequently, the real terahertz field E T Hz (t, x, z) may be obtained by an inverse temporal Fourier transform as delinated in Eq.8c. Here Θ(Ω) represents the Heaviside function. It delineates the fact that E T Hz (Ω, x, z) calculated for Ω ≥ 0 has to be reflected about Ω = 0 in order to obtain the real field. By noting that E T Hz (−Ω, x, z) = E * T Hz (Ω, x, z), it is evident that E T Hz (t, x, z) shall be a real valued function. Equation 8d corresponds to optical-to-terahertz conversion efficiency η obtained by integrating the spectral intensity over space and angular frequency. By Parseval's theorem, this could also be calculated in transverse momentum k x and/or temporal t domains (See Appendix).
The validity of Eq.8a maybe verified by noticing that, in the absence of loss, it reduces to the all familiar form proportional to the sinc function presented in Eq.9.
Equation 8a accounts for effects of loss, dispersion in both optical and terahertz frequency ranges, spatiotemporal distortions such as GVD-AD and spatial chirp as well as spatial walk-off between terahertz and optical pulses. Importantly, since no assumption is made on limiting the range of k z , this expression does not contain paraxial approximations and is thus generally valid for all pump beam sizes and propagation distances (i.e near and far-field). Furthermore, Eq.8a is equally valid for conditions of low and high absorption. It can be thus used quite generally to obtain meaningful predictions of terahertz efficiency, field strength, frequency and spatiotemporal profiles in the undepleted limit. In order to obtain physical intuition for the spatio-temporal properties of the generated terahertz radiation,closed-form expressions of the terahertz spectrum E T Hz (Ω, x, z) is desirable.
However, due to the presence of functions of k x in the denominator of Eq.8a, a general closed-form expression for E T Hz (ω, x, z) appears intractable. However, one may be obtained if the following conditions are satisfied.
If the transverse-momentum distribution P T Hz (Ω, k x , z) in Eq.7a is highly localized about k x = k(Ω)sinγ or when the phase-mismatch ∆k = 0, then the k x dependencies in the denominator of Eq.8a maybe eliminated. Formally, this constraint maybe delineated by setting ∆k α at k x = k(Ω)sinγ ± 2 √ 2/w, or when the transverse momentum distribution in Eq.7a reduces to e −1 of its value. Furthermore, accounting for the Ω dependence of w as given by Eq.7d , we obtain the following constraints : Equations 10a and 10b suggest that while on one hand beam sizes need to be large enough , they must also satisfy the condition that 2k T w 0 /τ 2 be small enough. In ensuing sections we will establish that very important ex-perimental situations fulfill these conditions, thus making the forthcoming expressions relevant to practical situations.
Since Eq.10a requires transverse momentum distributions to be highly localized (possess large beam sizes), k z (Ω) maybe approximated paraxially as where k x is a slight transverse momentum variation about k x = k(Ω)sinγ. Invoking the above , one obtains the following expression for E T Hz (Ω, x, z) in Eq.11a upon taking the inverse Fourier transform of Eq.8a in the paraxial limit.
Examining Eq.11a, the complex exponential factors outside the square brackets clearly delineate a terahertz pulse propagating at an angle γ relative to the pump pulse. Further, note that the spectral intensity is inversely proportional to the absorption coefficient. The first term inside the square brackets represents the source term driven by the nonlinear polarization term in Eq.5a. The effective duration of this nonlinear polarization term is given by τ 1 as defined in Eq.11b. As is evident, τ 1 has a transverse variation, which delineates spatial chirp of the generated terahertz pulse proportional to the GVD-AD term k T . The extent of this transverse variation is accentuated at larger bandwidths or shorter τ . Furthermore, the effective terahertz pulse duration is enlarged by the effective GDD φ (Eq.7b). The x 0 term defined in Eq.7f does not feature in the above equations due to the assumption of small |k T |w 0 /τ 2 , which causes x 0 ≈ 0 according to Eq.7f. From Eq.11a, a finite radius of curvature R(z) is also imparted to the terahertz spectrum which is non-zero even for an infinite pump radius of curvature due to spatio-temporal coupling effects as delineated in Eq.7c.
The second term within square-brackets in Eq.11a on the subsequent line represents the propagating terahertz wave. The total field is thus a difference of the nonlinear polarization term and the propagating terahertz wave. Here, notice that since the terahertz pulse propagates at an angle γ with respect to the pump beam, the spatial profile is distributed about the line x − ztanγ = 0. The difference in the spatial distributions is what creates the walk-off effect. The effective duration of the propagating component in Eq.11a is given by τ 2 as defined in Eq.11c. Naturally, the distribution of τ 2 spatially is different from τ 1 but in addition, it increases due to group velocity dispersion in the terahertz region β T .

Spatial terahertz profiles
Using Eqs.11a-11c, we proceed to understand the spatial distribution of the generated terahertz radiation in this section. In Fig.3(a), we plot the terahertz fluence (or spectral intensity integrated over all terahertz angular frequencies) for a lithium niobate crystal with w 0 = 2mm and τ 0 = 500fs at T=300 K in the x − z plane. The parameters from Table.II are used for calculations.
The white-dotted line in Fig.3(a) delineates the crystal boundary and the optical pump beam is centered about x = 0. The fluence distribution depicts an asymmetric beam profile, which may be understood by examining Eq.11a and Figs.3(b)-(d).As previously discussed, in Eq.11a, the terahertz spectrum E T Hz (Ω, x, z) is the difference between two transverse spatial profiles. The first term inside the square brackets of Eq.11a is centered The distribution of ET Hz (Ω) in x is the difference of two profiles centered at x = 0 and x = ztanγ evident from Eq.15a. At z=0, the two terms cancel out. (c). For larger z = 2cosγ/α, the displaced and attenuated second term produces an asymmetric total transverse profile of larger strength. This is the case for z < 1mm in (a).(d) For large enough z, the second term in Eq.11a is attenuated sufficiently and displaced faraway from x = 0 to result in the ET Hz (Ω, x, z reaching a maximum value. This is the case for z > 2mm in (a). about x = 0 and maybe viewed as source (blue curve in Fig.3(b)-(d)) or near-field term( [22]. The second term inside the square brackets of Eq.11a, is centered about x = ztanγ and corresponds to the radiated or propagating term (red curve in Fig.3(b)-(d)), further evident by the fact that it suffers attenuation as it propagates. The total transverse profile is the difference of these two profiles.
As illustrated in Fig.3(b), at z = 0, the two transverse profiles are identical and cancel each other out, resulting in the expected value of E T Hz (Ω, x, z = 0) = 0 (blackdotted curve in the insets of Fig.3). As z > 0, the propagating term is displaced to positive x > 0, which results in imperfect cancellation to produce an asymmetric total profile as seen by the black-dotted curve in Fig.3 (b) and also in the fluence map for z < 1mm. A growth in intensity from 0 is seen due to the attenuation of the propagation term in Eq.11a, which then reaches a maximum beyond the absorption length ∼ 2cosγ/α. While for small z, the asymmetry in the transverse profile is a result of the physics of the nonlinear process described in Figs.3(b)-(d), for large z, only parts of the beam for x < 0 lie within the crystal.

Efficiency variation with length
Insights into the optimal interaction length can be gleaned by evaluating the conversion efficiency η from Eq.8d with the terahertz spectral profile E T Hz (Ω, x, z) obtained from Eq.11a as follows : The first z dependent term inside the square brackets of Eq.12a, represents saturation of conversion efficiency due to absorption. The second term represents limitations due to walk-off. For beam sizes w 0 α −1 , the second term dominates and the system would be walkoff limited. The optimal interaction length would thus be z opt ∼ w 0 /tanγ. For large beam sizes w 0 α −1 , the first term would dominate and the system would be primarily absorption limited with z opt ∼ 2cosγ/α. In practice, cascading effects reduce the interaction length further, which are not considered in this analytic study. For lithium niobate, at room temperature, conversion efficiencies are lower and the situation is closer to the undepleted case.
While Eq.12b does not have an explicit form in general, it may be evaluated for constant α and |k T |w 0 /τ 2 1. For large enough z, the conversion efficiency η maybe approximated by the following : Note, the reduction of conversion efficiency for large k T (GVD-AD) and beam radii. In the absence of GVD-AD, the conversion efficiency should improve with pump bandwidth or shorter τ . However, due to the effects of GVD-AD, an optimal pulse duration τ opt ≈ 1.62(k T w 0 ) 1/2 exists.

Terahertz frequency
Maximizing Eq.12b with respect to terahertz angular frequency Ω, one may obtain the following expression for the central terahertz frequency Ω max for the case when k T 2w 0 /τ 2 1 as shown below : FIG. 4. Comparison of conversion efficiencies η calculated from expressions for spectra in Eq.11a and Eq.8a for a lithium niobate crystal. (a) At T=300 K, ∆η ≤ 20% if constraints in Eq.10a (Black-dotted) and Eq.10b (Blue-dotted) are satisfied. (b) At T=100 K, relative error is large predominantly due to violation of constraint in Eq.10a due to reduced absorption at cryogenic temperatures.
Once again, the reduction of average frequency for larger bandwidths and beam radii is evident for systems with non-zero GVD-AD. This trend is consistent with experimental observations for tilted-pulse-front experiments in lithium niobate.

Validity of analytic expressions of ET Hz (Ω, x, z)
While Eq.8a is generally valid and does not incur significantly more computational cost compared to evaluating Eq.11a, the latter is appealing due to being explicit and more intuitive. We have already laid out the conditions when Eq.11a can be employed in Eqs.10a-10b. Here, we plot the relative error in conversion efficiency ∆η between Eqs.11a and 8a for terahertz generation in lithium niobate in Fig.4. Parameters for the calculations are provided in Table.2.
From Fig.4(a), it can be seen that for T=300 K or in the large absorption limit, the closed-form expression in Eq.11a is within 20% of the general expression in Eq.8a for regions approximately bounded by the constraints in Eq.10a and 10b. They are particularly accurate for ≈ cm beam sizes, making them attractive for analyzing high energy terahertz generation setups. For instance, at τ 0 ≈ 0.5ps ,f T Hz ≈ 0.5 THz, we obtain the threshold beam radius w 0 ≈ 7 mm from Eq.10a, which is consistent with Fig.4(a). However, Eq.11a is less efficacious for very small beam radii or short pulse durations.
We set α = α(2/τ 0 ) in Eq.10a, since the central terahertz frequency is approximately Ω max ≈ 2/τ 0 (Eq.14). Therefore, the threshold beam size is larger for longer pump pulses due to the reduced terahertz frequency and hence smaller absorption.
Naturally, the discrepancy between Eqs.11a and 8a is much larger for T=100 K due to a much smaller absorption coefficient of lithium niobate at cryogenice temperatures. For 0.5 THz, the absorption coefficient at T=100 K ≈ 130/ cm which is > 7 times smaller compared to that at T=300K (≈ 960/cm). This translates to the critical beam size for low relative error at T=100 K being about 5.5 cm, which is out of bounds in the parameter space depicted in Fig.4(b).

F. Temporal profiles
Here, we present closed-form expressions for terahertz transients when the constraints supplied by Eqs.10a-10b are satisfied. This may be obtained by taking an inverse Fourier transform in the temporal domain of the expression provided by Eq.11a. To obtain, a closed-form expression, the dispersive properties of the absorption coefficient are neglected and α is set to α(2/τ 0 ), since the central terahertz frequency is roughly given by 2/τ 0 (Eq.14).The results are presented in Eqs.15a-15c below.
Firstly, note that the intensity along the plane defined by t = t − xsinγv −1 T Hz − zcosγv −1 T Hz is constant due to an obliquely propagating terahertz pulse. Secondly, the finite radius of curvature slightly modifies this plane to a curved surface as is evident in Eq.15b. Secondly, as already delineated in Eqs.11b-11c, the pulse duration of the terahertz transients varies along transverse spatial dimension x. The degree of spatial variation increases with pump bandwidth or shorter τ .
We employ Eq.15a to plot the evolution of the terahertz electric field in lithium niobate at T = 300K for for two different durations τ = 500fs in Fig.5 and τ = 50fs in Fig.6. As can be seen in Fig.5, the field grows as it propagates at an angle γ = 63 • with respect to the pump, while evolving into a single-cycle terahertz field. Across the tilt-plane, the duration of the pulse does not vary noticeably, with relatively uniform properties.
However, for τ 0 = 50fs, in Fig.6, the field is seen to grow as was the case for Fig.5 but the degree of asymmetry in pulse duration along the tilt-plane is greater due to the effects of GVD-AD as expected from Eqs.11b and 11c.
By inspecting Eq.15a, it is clear the peak field E T Hz,max occurs at the maxima of the function t e −2t 2 /τ 2 which is at t = τ /2. Thus the peak field is given by the following expression in Eq.16 : n T Hz α(2/τ 0 )cτ 3 0 (16)

Validity of closed-form expression for transients
In Fig.7, we test the proximity of the bounds provided by the peak field in Eq.16 to those obtained from the general expressions in Eqs.8a-8c for lithium niobate . Similar to Fig.4, the agreement is within 20% at T=300 K for regions bounded approximately by the constraints in Eqs.10a -10b. However, just as in the case of Fig.4(b), the relative error is large for T=100 K in the parameter space of w 0 − τ 0 depicted. Hence, the closed form expression of transients are also more useful at room temperature for high energy tilted-pulse-front setups (i.e large beam radii). In Table.I, we summarize the various equations and their domains of applicability.

III. RESULTS AND DISCUSSION
In this section, we analyze important features of terahertz generation with tilted-pulse-fronts in lithium niobate. We first illustrate the importance of spatial-chirp and radius of curvature from Eq.6 to the tilt angle γ. Subsequently, we employ Eqs.8a-8c to evaluate terahertz spectra E T Hz (Ω, x, z), conversion efficiency η as well as peak electric fields at T=300 K and T= 100 K for various values of beam radius w 0 and pump pulse duration τ 0 . We compare cases with and without GVD-AD k T to illustrate its effects on terahertz frequency, conversion efficiency and peak electric field. We find that for larger beam radii and shorter pump pulse durations, there is a decrease in all these quantities due to the detrimental effects of GVD-AD. Overall the undepleted models provide good qualitative predictions and understanding of terahertz generation using tilted-pulse-fronts. However, while quantitative predictions for T=300 K are reasonable , those at T=100 K are overestimated. Although absorption reduces significantly at cryogenic tempera-tures, cascading effects induce limitations which necessitates full depleted calculations for these conditions.

A. Tilt angle due to spatial-chirp and radius of curvature
The importance of the contribution of the radius of curvature and spatial-chirp to the tilt-angle is particularly evident in considering the effect of the lens-tocrystal distance s 2 (see Fig.1) on conversion efficiency. Experimentally, the displacement of the crystal from the optimal imaging distance s 2 results in a dramatic loss in conversion efficiency [20]. In Fig.8, we plot the tiltangle of the pump beam just inside the crystal for various values of s 2 using Eq.6. A tilted-pulse-front setup with grating of p = 1500 lines/mm, grating incidence angle θ i = 56.16 • , lens focal length f , magnification M = 0.4775, s 2 = (1+M )f and s 1 = s 2 /M was assumed. Further, the pulse properties corresponded to a duration τ 0 = 50 fs and beam radius w 0 = 1 mm. The values of ζ, β 0 , R 0 (z) were calculated using dispersive ray-pulse matrices, following Martinez [30].
As expected,the angular dispersion of the beam does not change after the focus (which lies before the crystal) as is evident by the flat blue line in Fig.8. Secondly, the contribution to the tilt-angle from the spatio-temporal chirp is found to be infinitesimally small (on the order of 10 −6 even for group delay dispersion φ 0 on the order of 5×10 5 fs 2 (for τ 0 = 50fs). However, the contribution from the final term of Eq.6 comprising of spatial-chirp and the radius of curvature is responsible for a significant change in tilt-angle as can be seen in Fig.8. This variation arises from changes in both spatial chirp ∆ζ = n(ω 0 )β 0 ∆s 2 as well as changes in the radius of curvature with displacement of the crystal ∆s 2 . As Fig.8 suggests, the use of longer focal lengths (green curve) reduces the sensitivity of the setup to displacements ∆s 2 , making it easier to obtain optimal performance.
Inside the crystal, this term has a much smaller impact since the variation of spatial chirp with distance z is smaller due to reduced angular dispersion inside the material (by Snell's law), i.e. ∆ζ = β 0 z. Furthermore, the variation of the radius of curvature is also reduced. Thus, the impact on interaction length within the crystal due to varying tilt-angle is much lesser (black,dashed line in 8).If a near collimated beam (with R = ∞ can be produced, then the impact of this effect will of course be negligible. Thus, the use of telescopic setups for imaging may indeed prove to be advantageous as prior experiments [31] seem to suggest.

B. Spectra and frequency
Using Eq.8a and 8b and parameters from Table.II, we first plot the spatial distribution of terahertz spectra in Fig.9. All calculations do not consider the effects of the FIG. 5. Spatio-temporal evolution of terahertz transients for a τ0 = 500 fs, w0 = 2mm beam radius in lithium niobate at T=300 K. The propagation of the terahertz transient at an angle γ ≈ 63 • is evident in panels (a)-(d). A continuous growth in intensity is seen. Due to the large duration of the pump pulse, spatial inhomogeneities are relatively low.
FIG. 6. Spatio-temporal evolution of terahertz transients for a τ = 50fs, w0 = 2mm beam radius in lithium niobate at T=300 K. Due to the short duration of the pump pulse, spatial inhomogeneities are obviously present in (d).
FIG. 7. Relative error of peak field calculations between Eqs.15a and 8c. Consistent with Fig.4, the relative error is ≤ 20% at T=300 K upon satisfying the constraints in Eqs.10a (black-dotted) and 10b (blue-dotted). (b) At T=100 K, the relative error is large, analogous to Fig.4(b) due to violation of Eq.10a resulting from small absorption at cryogenic temperatures.
crystal boundary on the spectral properties. In Fig.9(a), the spatial distribution of the terahertz spectrum for a beam with radius w 0 = 2mm and pump bandwidth τ 0 = 500fs is depicted at the point of maximum efficiency. It can be seen that the spatial distribution of the spectrum is rather uniform in relation to the case of τ 0 = 50fs in Fig.9(b), where the spectrum and center frequency reduce for larger values of x. The greater degree of inhomogeneity for τ = 50fs, is consistent with the spatio-temporal snapshots of terahertz transients contrasted in Figs.5 and 6.In Fig.9(c), the average spectrum corresponding to Fig.9 is depicted. In addition, the average spectrum for a larger beam radius of w 0 = 1cm is also shown. Clearly, the average frequency drops for larger beam sizes as is evident from larger terahertz durations anticipated for larger beam radii in Eqs.11b-11c. In Fig.9(d), average spectra for τ 0 = 50fs are shown. Here, the beam radius of 1 cm produces a greater reduction in average frequency due to the greater impact of GVD-AD for shorter pump durations.
As is clearly evident, the strongest influence on central frequency Ω max is the bandwidth and beam radius w 0 . While a shorter duration τ 0 produces higher terahertz frequencies of Ω max ≈ 2/τ 0 for k T = 0, the effect of GVD-AD causes this value to reduce upon continuous increase of either bandwidth of beam radius w 0 .

Equation Description
Validity Eqs.8a-8c Closed form for ET Hz (ω, kx, z) High(Quantitative) and low absorption (Qualitative). Eqs.11a,12a Closed-form for ET Hz (Ω, x, z) High (Quantitative) absorption for Eqs.10a-10b Eqs.15a, 16 Closed-form for ET Hz (t, x, z) High (Quantitative) absorption for Eqs.10a-10b.  FIG. 8. Variation of pulse-front-tilt angle γ just inside a lithium niobate crystal with lens-to-crystal distance variations ∆s2. Angular dispersion(blue) does not change with s2 as expected. Variations of pulse-front-tilt due to spatial-chirp ζ(z) and radius of curvature R(z) (third term, Eq.6) however produces significant variations in tilt-angle. The variations are more severe for short focal length of f = 7.5cm(red) compared to longer focal length of f = 15cm(green). Variation of tiltangle due to spatial-chirp and radius of curvature inside the crystal is negligible due to reduction in growth of ∆ζ(z) due to reduced angular dispersion (reduces by a factor nω0 due to Snell's law) as well as slower change of radius of curvature.
In Fig.10, the central terahertz frequencies for T=300,100K with and without the effects of GVD-AD are depicted using calculations employing Eq.8a.
For T=300 K and k T = 0, the maximum frequency occurs at the shortest durations and show no change upon varying beam radius w 0 . However, for finite values of GVD-AD, the maximum frequency no longer occurs at the shortest duration but at a slightly longer pump duration as is evident in Fig.10)(b). Furthermore, the dramatic drop in peak frequency with beam radius is ev- FIG. 9. Terahertz spectra obtained from Eqs.8a-8b in lithium niobate at T=300 K. (a) Spatial distribution of the terahertz spectrum for a pump pulse duration of τ0 = 500fs and w0 = 2 mm shows relatively homogeneous properties consistent with transients in Fig.5.(b) Spatial distribution of the terahertz spectrum for a pump pulse duration of τ0 = 50 fs and w0 = 2 mm shows significant asymmetry analogous to Fig.6.(c) Average spectra |ET Hz (Ω, x, z)| 2 dx for varying beam sizes shows the reduction in frequency bandwidth for larger beam radii due to the effects of group-velocity dispersion due to angular dispersion (GVD-AD).(d) Average spectra for τ0 = 50fs shows even greater reduction of frequency due to enhancement of GVD-AD effects at short durations.(See Eq.14).
ident. In Figs.10(c)-(d), the situation for T=100 K is shown. Firstly, in the absence of GVD-AD, the average frequency shows appreciable increase in contrast to the T=300 K case due to reduction in absorption. However, for the case of finite GVD-AD, the reduction in absorption by cryogenic cooling does not appear to influence the peak frequency much, indicating that the frequency is mainly GVD-AD limited.

C. Conversion Efficiency
Using Eq.8a, the conversion efficiency η is then evaluated for the four cases outlined above, i.e. for T=300,100 K with and without the inclusion of the effects of GVD-AD for various values of τ 0 , w 0 . The results are plotted in Fig.11 below.
For T=300 K and k T = 0 (or no GVD-AD), conversion efficiency saturates after a certain minimum value of w 0 , which is consistent with the calculations presented in [22,23]. An optimal duration exists since for shorter durations, the peak frequency is larger which corresponds to larger absorption coefficients. For the case when GVD-AD effects are included in Fig.11(b), firstly conversion efficiency does not saturate with w 0 but instead shows a decline beyond a certain value of w 0 . Furthermore, the optimal pump durations shift to larger values. In Fig.11(c), the case for T=100 K and k T = 0 is shown. Compared to Fig.11(a), the beam radius at which efficiency saturation has increased due to a reduction in absorption, which yields larger interaction lengths. A reduction in the optimal pulse duration accompanying an increase in frequency due to reduced absorption is also evident. In Fig.11(d), the reduction of absorption has K, effects of GVD-AD increase optimal duration as alluded to by Eq.13a. A reduction rather than saturation at large w0 is seen (c) For kT = 0, T = 100K, saturation occurs at much larger w0 due to increased interaction lengths by virtue of reduced absorption.(d) For |kT | > 0, peak efficiency occurs at larger τ compared to (c), consistent with (b) increased the interaction length and moved the optimum to larger beam sizes but due to marginal change in frequency, the optimal pulse durations have not changed much. The latter trend is consistent with experiments but the former is not, since at cryogenic temperatures, terahertz generation is limited by cascading effects more adversely compared to absorption.
Similar to the previous section, the influence of GVD-AD on η is evident from the closed-form expression for η presented for the case of k 2 T w 2 0 /τ 4 1 in Eq.13a. It is therefore evident that in the absence of absorption, conversion efficiency should increase with decreasing τ . However, for finite k T , the second term in the square brackets reduces the efficiency for larger beam size and reducing τ , counteracting the τ −3 factor outside the brackets.
If one accounts for a scaling factor of 1/ √ 2 due to the third spatial dimension y and an additional factor of 1/2 due to Fresnel losses, then the maximum conversion conversion efficiency in Fig.11(b) for T=300 K is ≈ 0.88% at τ ≈ 400fs and w 0 = 5 mm. If the expression for optimal pulse duration from Eq.13b is used, then one obtains the value to be τ opt = 413 fs, which is very close to the simulated results. This is also consistent with experimental trends [13]. However, the conversion efficiencies from this undepleted model are a bit higher even at room temperatures since cascading effects are not accounted for. These cascading effects in combination with GVD-AD deteriorate phase-matching and reduce the conversion efficiency [25]. At T=100 K, the conversion efficiency is signifi- cantly overestimated. Here, the impact of cascading shall limit performance even further. In summary, undepleted calculations provide good qualitative understanding and fair quantitative predictions at room temperature while overestimating conversion efficiency at cryogenic temperatures.

D. Peak field
These values are proportional to the ratio of the product of η and Ω max from Figs.11 and 10. For T=300K and k T = 0, ω max is largest for shortest pump durations while η is optimized for slightly durations. Therefore, the peakfield strengths show a flatter profile as a function of τ in Fig.12(a). In Fig.12(b)-(d), a similar effect is observed. Consistent with Fig.11(d), the values of peak-field for T=100 K in Fig.12(d) are inflated due to operating in the undepleted limit.

IV. CONCLUSION
An undepleted analysis of terahertz generation by optical rectification of tilted-pulse-fronts was performed. The analysis accounted for the effects of various spatiotemporal distortions imparted to the optical pump beam by the tilted-pulse-front setup such as angular dispersion, spatial and temporal chirp, group velocity dispersion due to material and angular dispersion as well as the radius of curvature of the pump beam. Closed form expres-sions to evaluate the properties of the generated terahertz radiation were provided in the transverse momentum domain, which are generally valid for all bandwidths and beam sizes. In addition, closed form expressions for terahertz spectra and transients for small k 2 T w 2 0 /τ 4 were provided. These closed-form expressions in spatial and temporal domains provide clear qualitative understanding of the physics at work and accurate quantitative predictions at room temperature for high energy terahertz generation systems. The effect of radius of curvature on tilt angle was formally derived and presented along with illustrative calculations relevant to widely employed experimental setups. The overarching point shown via the analysis was the detrimental effect of GVD-AD on beam profiles, terahertz frequency, efficiency and scaling with beam-size. While, this was numerically shown in previous work, here we show this conclusively using analytic formulae in multiple spatial dimensions. The undepleted analysis here is able to provide useful insights into terahertz generation using tilted-pulse-fronts. In the large absorption limit, such as the case for terahertz generation at room temperature in lithium niobate, the model provides reasonable quantitative predictions. In the low absorption limit, the effects of cascading become more important, which would require full numerical simulations.

A. Fourier Transform relations
We begin by decomposing the real scalar electric in space (r) and time (t) as E(r, t) = [f (r, t)e −jω ref t + f * (r, t)e jω ref t ]/2. The corresponding Fourier transform of the real field is given by E(r, ω) = [f (r, ω − ω ref ) + f * (r, ω + ω ref )]/2. The reference angular frequency ω ref is used here only to delineate the fact that the spectrum of the real field comprises of positive and negative frequency components, symmetric and conjugate about ω = 0. It may be dropped henceforth. Furthermore, the dependency on r shall be assumed to be implicit and will also be dropped. From the above it is clear that f (ω) is a Fourier transform of f (t). The conventions in Eq.17 are assumed for Fourier transforms between temporal and spectral domains.
Similarly, for spatial Fourier transforms in the transverse dimension x are provided by Eq.18. Input angular dispersion rads-s β Effective angular dispersion rads-s kT Group velocity dispersion due to angular dispersion s 2 m −1 x0 Central beam position m The energy in each space/time is equal to that in transverse momentum/angular frequency according to Parseval's theorem in Eq.19.
In evaluating Eq.1a, it would suffice to book-keep for terahertz angular frequencies Ω > 0 for E T Hz (Ω, x, z), which corresponds to the f (ω −ω ref ) term above. Therefore, in evaluating the real electric field E T Hz (t, x, z), we adopt Eq.8a. In practice, since E T Hz (−Ω, x, z) = E * T Hz (Ω, x, z) readily, one may evaluate E T Hz (t, x, z) = F −1 t E T Hz (Ω, x, z)/2 (for positive and negative values of Ω).
The energy per unit length of real electromagnetic field vectors E(electric) and H (magnetic) passing through an area of cross-section with unit vector dxdyẑ and refractive index n is given by For integrals of the following form : z 0 P T Hz (Ω, k x , z )e j∆kz + α 2 z dz = P T Hz (Ω, k x , z)e j∆kz+ α 2 z α/2 + j∆k − P T Hz (Ω, k x , z)e j∆kz+ α 2 z (α/2 + j∆k) 2 + P T Hz "(Ω, k x , z)e j∆kz+ α 2 z (α/2 + j∆k) 3 ..
Naturally, for relatively large α and P T Hz (Ω, k x , z)/P T Hz (Ω, k x , z), only the first term in the expansion shall be significant. The variation along space of the Polarization term due to material is typically small while absorption coefficient values even at cryogenic temperatures for lithium niobate 126/m are quite large to result in the validity of this approximation.