Terahertz generation by beamlet superposition

We analytically show how a superposition of beamlets produces terahertz radiation with greater spatial homogeneity and efficiency compared to tilted-pulse-fronts generated by diffraction gratings. The advantages are particularly notable for large pump bandwiths and beam sizes, alluding to better performance in the presence of cascading effects and higher energy pumping. A theory of terahertz generation using a superposition of beamlets is developed. It is shown how such an arrangement produces a distortion free tilted-pulse-front. Closed form expressions for terahertz spectra and transients in three spatial dimensions are derived. Conditions for obtaining performance parity and bounds for optimal parameters are furnished.


Introduction
The large peak fields and long wavelengths of terahertz radiation make them attractive for a number of fundamental scientific investigations [1,2] and technological applications. On the technological front, compact particle acceleration [3,4,5,6,7,8], electron microscopy and beam diagnostics [9,10] have been predicted to benefit significantly from advances in terahertz generation. While a number of modalities for terahertz generation exist [11,12], nonlinear down conversion of near infra-red optical radiation using second order nonlinearities in electro-optic materials has achieved great success, producing percent-level conversion efficiencies [13] and milli-joule level pulse energies [14]. In particular, the long assumed efficiency limitation due to the large disparity in optical and terahertz photon energies is overcome by cascaded difference frequency generation, whereby a pump photon can be repeatedly down-converted to yield multiple terahertz photons as long as phase-matching conditions are satisfied [15].
Compatibility with off-the-shelf Titanium Sapphire and Ytterbium laser technology has made terahertz generation using tilted-pulse-fronts [16,17,18,19] in lithium niobate ubiquitous. However, the use of diffraction gratings makes imaging of tilted-pulse-fronts with large pump bandwidths, a challenging task. Furthermore, the large group-velocity dispersion due to angular dispersion in such systems inhibits terahertz generation once significant cascading effects set in [20,21]. In short, tilted-pulse-fronts produced by diffraction gratings are severely bandwidth limited.
An approach to circumvent these limitations using a superposition of discrete beamlets was demonstrated in [22] and variants thereof were studied in [23,24,25]. In this approach, terahertz radiation is generated at the Cherenkov angle by each beamlet due to the spread in transverse momentum created by the small beamlet radius. The superposition of radiation from various beamlets, then combine to form a terahertz 'plane wave'. Due to uniformly distributed transverse momentum vectors, this approach may have the potential to circumvent bandwidth limitations due to cascading effects [26].
Here, we present an analytic formulation in three spatial dimensions in the undepleted limit that establishes the advantages offered by beamlet superposition in relation to tilted-pulse-fronts generated by  Figure 1: Schematic of terahertz generation using a super position of beamlets of duration τ , separated transversely by a distance ∆x and longitudinally by v g ∆t. Each beamlet emits radiation at the Cherenkov angle γ, which sets the condition for coherent super position of radiation to be v g |∆t| = |∆x|tanγ diffraction gratings. In the presented analysis, effects such as absorption, diffraction of optical and terahertz beams, beam curvature effects and interaction within and between beamlets are accounted for .In particular, closed form expressions are derived which predict terahertz spectra without spatial inhomogeneities, larger terahertz frequencies and higher conversion efficiencies. In particular, the relative performance of beamlet superposition is superior at larger pump bandwidths and beam sizes. The superior performance at larger pump bandwidths also indicates greater insensitivity of beamlet superposition to cascading effects. The general physics of the approach is elucidated, conditions for obtaining performance better than conventional tilted-pulse-fronts and the limitations of the approach are discussed. Guidelines to optimize various parameters are furnished.

Tilted-pulse-fronts by beamlet superposition
In Fig.1, a schematic depicting the relevant system is presented. In general, a series of beamlets with duration τ , separated in increments of transverse distance ∆x and delayed in increments of time |∆t| = |∆x|tanγ/v g propagate with group velocity v g in the z direction in an electro-optic material with a second order nonlinear coefficient χ (2) .
The total electric field of all beamlets corresponding to this schematic is presented in Eq. 1a in the frequency domain. Here x m is the position of the m th beamlet, while ∆zv −1 g is the corresponding temporal delay. Based on Fig.1, beamlets at larger transverse positions are located at smaller longitudinal positions. Therefore, we have ∆z = −∆xtanγ.
Material dispersion is accounted for by the k m ∆ω 2 term in Eq.1a, where ∆ω = ω − ω 0 represents displacement from the central pump angular frequency ω 0 . Each beamlet has an e −2 radius defined by σ, where σ is significantly smaller (< 100 µm) than the typical beam size (> mm) used for terahertz generation with conventional tilted-pulse-fronts. Due to the small size of beamlets, they have a Rayleigh length comparable to the typical terahertz absorption length ∼ cm, which makes consideration of longitudinal variations of beamlet size and radius of curvature R(z) necessary as delineated in Eqs.1b-1c. In Eqs.(1c-1d),z 0 represents the location of the beam waist , while z R corresponds to the Rayleigh length k(ω 0 )n(ω 0 )σ 2 /2. The pump Figure 2: (a)E(ω, k x ) of a single beamlet with τ = 50 fs and σ = 50 µm shows a Gaussian spread in transverse momentum and frequency.(b)N = 20 beamlets separated by 100 µm appear as a series of oblique lines in ω −k x space. The slope of the lines are n g tanγc −1 , where γ is the pulse-front-tilt angle. The spacing between lines is inversely proportional to the spacing between successive beamlets ∆x in real space. (c) When the spacing between lines is reduced to σ/2 = 25 µm, the spacing between oblique lines increases to an extent that only a single line falls within the Gaussian spread of k x . This represents a case of a distortion free tilted-pulse-front.
spectrum is assumed to be Gaussian, centered about ω 0 as depicted in Eq.1d.
To visualize the possibilities that result from the above pulse format, let us consider its behavior in the ω − k x plane (where k x corresponds to the x directed component of transverse momentum) by examining the spatial Fourier transform of Eq. 1a. Without loss of generality, we set E m (ω) = A(ω). For illustrative purposes, we also ignore the radius of curvature in Eq.2a,2b and variations along the y dimension, since the key physics involves only the x − z plane.
Equation 2a simplifies to Eq.2b. The ratio of sinusoidal functions in Eq. 2b is maximized along oblique lines of the form k x = lπ/∆x + ∆ωn g tanγc −1 in the ω − k x plane, where l is an integer. Naturally, the magnitude of lines drops exponentially from k x = 0 due to the e −k 2 x σ 2 /4 pre-factor in Eqs. 2a-2b.The separation of lines is inversely proportional to the spacing between beamlets ∆x. Figure 2 represents normalized values (within each panel) of E(ω, k x , z) for various cases. In Fig.2(a), the case of a single beamlet of pulse duration τ = 50 fs and e −2 beamlet radius σ = 50 µm is depicted. The small size of the beamlet produces a large transverse momentum spread while the large bandwidth produces a significant spread in frequency. Figure 2(b) shows the case of N = 20 beamlets separated by a distance ∆x = 2σ = 100 µm. The separation between location of maxima in the ω − k x plane is inversely proportional to the separation of beamlets. Therefore, E(ω, k x , z) appears as a series of oblique lines.(c)The Figure 3: ω − k x plots of the field from Eq.4b for a τ = 50 fs,w = 0.5 mm pulse for different values of group velocity due to angular dispersion (a)k T = 0(b)k T = −50 × 10 −23 s 2 (c) −25 × 10 −23 s 2 .The first case represents a distortion free tilted-pulse-front. The subsequent cases do not satisfy phase-matching conditions over the entire bandwidth of the incident pulse.
spacing between beamlets is reduced to 25 µm. When the spacing between beamlets is sufficiently small, exactly a single oblique line falls within the transverse momentum spread defined by the e −k 2 x σ 2 m /4 pre-factor in Eqs. 2a-2b.Since, ∆z = −∆xtanγ, the line for l =0 is identical to a distortion-free tilted-pulse-front (See Fig.3(a)).
The condition for concentrating ≈ 90% of the spectral intensity (or a fraction of 1 − e −2 ) in the central oblique line is obtained by setting the separation between lines π/∆x according to e −(σπ/∆x) 2 /4 = e −2 . This provides the condition for preferred spacing of beamlets according to Eq.3.

Comparison to tilted-pulse-fronts from diffraction gratings
Equations 4a and 4b represent the fields due to a diffraction grating in the ω−x and ω−k x planes respectively.
Inspecting Eq.4b, it is evident that the distribution of the field in the the ω − k x plane is centered about k 0 β∆ω + k T ∆ω 2 , where β = dθ/dω is the angular dispersion term and k T defines the group-velocity dispersion due to angular dispersion (GVD-AD). For k T = 0, or in the absence of GVD-AD, the pattern formed in the ω − k x plane would be perfectly linear.
In Fig.3(a), a distortion free tilted-pulse-front for k T = 0 is plotted . For a tilt-angle of γ, it can be shown from Eq.4b that k 0 βv g = tanγ. Thus,a series of delayed beamlets (Fig.2(c))form a distortion free tilted pulse front in the ideal limit. Furthermore, notice that the distributions have a narrower spread due to the large value of beam size w in contrast to the beamlet radius σ.This is perhaps one advantage of tilted-pulse-fronts formed by gratings, where all the energy is concentrated at the desired transverse momentum value.
However, due to the non-zero values of k T (Figs.3(b)-(c)), one obtains a nonlinear distribution in ω − k x space for large detuning ∆ω from the central frequency of the incident optical field ω 0 . In other words, for equal increments in frequency, the angular separation is unequal for a field with finite GVD-AD.
A consequence of this nonlinear distribution is that phase-matching conditions are not satisfied across the bandwidth of the pulse.
for different values of pulse duration and beamlet separation ∆x/σ. The peak pump electric field is kept constant for all ∆x/σ for a given value of τ . When v g τ is large as is the case for τ = 500fs, the consideration of interaction between beamlets of the type depicted in Fig.4(b) become important for small ∆x/σ. On the other hand, for τ = 50fs, interaction within beamlets dominate the overall nonlinear polarization due to small values of v g τ .

Nonlinear polarization
In order to obtain expressions for the generated terahertz spectra and transients, we first obtain the nonlinear polarization term P T Hz (Ω, x, y, z) = ∞ −∞ E(∆ω + Ω, x, y, z)E * (∆ω, x, y, z)d∆ω at angular frequency Ω, which drives the terahertz generation process. After much algebraic manipulation, the expression for the most general case is provided in Eqs.5a-5c.

Bandwidth effects
In the equation above, χ (2) is the effective second order susceptibility of the nonlinear material. The first exponential term describes the terhertz spectral content. The effective duration is τ 2 /2 + 8β" 2 /τ 2 . Here, β", defined in Eq.5b represents the effects of any temporal chirp and material dispersion which would tend to spread the optical pump pulse in time and thereby reduce the generated terahertz frequency. Thus, their effects manifest as an increase in effective pulse duration, proportional to the initial bandwidth as is evident in the β" 2 /τ 2 dependence in Eq.5a. The second exponential term, denotes that the nonlinear polarization propagates at the group velocity cn −1 g in the z direction.

Interaction between beamlets
The summation over indices m, k delineates contributions due to interaction between beamlets m, k. Note that each of these polarization terms contains an average delay proportional to the average longitudinal position z m,k as delineated by the first exponential term within the summation. The second and third terms suggest that the nonlinear polarization terms decay exponentially with transverse separation ∆x m,k in relation to beamlet size σ and longitudinal separation ∆z m,k in relation to pulse length v g τ .This is intuitive in that the interaction between beamlets becomes exceedingly weak when their longitudinal separation significantly exceeds the pulse duration.
The third and fourth terms correspond to transverse spatial variations of the nonlinear polarization terms. Note that each of these terms is centered at the average transverse position x m,k of beamlets m, k.
The relative strength of interactions are depicted schematically in Fig.4. Interactions within beamlets are the strongest as shown in Fig.4(a). Moderately strong interactions occur when transverse separation is large but longitudinal overlap is still significant. This is the case when v g τ ≈ ∆z m,k . In Fig.4(c), the situation when interaction between beamlets becomes increasingly weak with spatial as well as temporal separation is shown. If only interactions of the type in Fig.4(a) are relevant, Eq. 5a can be reduced to a form which is purely a super position of radiation generated by each individual beamlet, i.e for m = k.
In Fig.4 , while keeping the peak pump electric for various values of ∆x, ∆z and τ = 500, 50 fs. Since ∆z ≈ v g τ in the former case, the full expansion over indices m, k (black,dashed) is necessary as the beamlets get closer or for ∆x/σ < 1. As ∆x ≈ σ, only nearest neighbour couplings(i.e. k − 1 ≤ m ≤ k + 1) become important (red,dashed). This region contains interactions of the form depicted in Fig.4(b) due to the large value of v g τ in relation to ∆z m,k while for even larger separations, only interaction within beamlets (blue,dashed) are required as seen by the convergence of the three curves.
For τ = 50f s, only interactions within beamlets are important, except at very small separations of ∆x/σ ≈ 0.1. This is due to the small value of v g τ which obviates consideration of interactions of the type depicted in Fig.4(b).
However, in general all interaction terms must be considered and interactions of the form represented by Eq.6 cannot be assumed.
What is particularly distinct from the case of conventional tilted-pulse-fronts is the absence of spatial dependence in the first exponential term in Eqs.5a and 6, which indicates relatively homogeneous spectral properties compared to tilted-pulse-fronts generated by diffraction gratings.

Terahertz spectra
Using the nonlinear polarization term in Eq.5a, one may obtain the following closed form expression for E T Hz (Ω, k x , y, z) in Ω − k x space by Fourier decomposition (See Appendix).
Note, the above expression can be utilized to evaluate the behavior of systems in general, in the undepleted limit.
Beamlet superposition The coherent superposition of beamlets becomes apparent upon examining the complex exponential terms within the summation. For ∆z = −∆xtanγ, the terms Ωn g ∆ z m,k c −1 + k x ∆ x m,k = 0 at k x = k(Ω)sinγ, leading to constructive interference of all beamlet interactions, for terahertz radiation propagating at an angle γ with respect to the pump beamlets. Phase-matching COHERENT SUPER POSITION Figure 5: (a) Larger terahertz frequencies are phase-matched at larger values of k z and hence require greater angular spread or equivalently, smaller beamlet radii σ. (b) Terahertz radiation from beamlet m − 1 only superposes with terahertz generated by beamlet m, due to its rapid diffraction. This requires the spacing between beamlets ∆x ≤ σ, consistent with Eq.3. The diffraction of terahertz radiation and optical beamlets places a lower bound on the value of σ The final term in Eq.7 represents the phase-matching condition and is maximized for k z = Ωv −1 g . In the absence of absorption, this reduces to the well known sinc(∆kz/2)z function, where ∆k = k z (Ω, k x ) − Ωv −1 g . Beam size dependencies From Eq.7, it is evident that to optimally utilize larger bandwidths to generate larger terahertz frequencies, a larger angular spread in beamlets is necessary . This is because, larger terahertz frequencies will be phase-matched at larger values of k z (See Fig.5(a)), which necessitates a greater spread in transverse momentum to overlap with the sinc∆kz/2 function. If the k x distribution in Eq.7 is set to be equal to e −1 at k(Ω)sinγ, it translates to the following condition for the maximum permissible beamlet radius σ.
σ ≤ λ T Hz πn T Hz sinγ (8) Terahertz and beamlet diffraction The small values of σ in relation to the wavelength of the terahertz radiation causes it diffract rapidly . However, if the beamlets are spaced by a distance ∆x ≤ σ as posited in Eq.3, the emitted terahertz radiation will be able to experience coherent growth by superposing with radiation emitted by adjacent beamlets. However, they do not coherently interfere with radiation emitted by beamlets farther away, as they fall outside the terahertz Rayleigh range as depicted in Fig.5(b).
Another salient feature, evident in Eq.7, is the longitudinal dependence on distance of the beamlet size σ(z) and radius of curvature R(z). The smaller beamlet radius σ results in smaller Rayleigh distances in comparison to conventional tilted-pulse-fronts produced by diffraction gratings, which results in nonnegligible diffraction. As shall be discussed later, the diffraction of terahertz radiation and pump beamlets place a lower bound on beamlet sizes.
Closed-form expressions While Eq.7 is useful in obtaining the distribution of the generated terahertz field in Ω − k x space, it does not directly lend itself to providing intuition in physical space or temporal domains. Complex functions of k x prevent obtaining closed form expressions in the general case. However, a situation when this is possible is in the large absorption limit, i.e. for α ∆k or when the distribtion in k x space is strongly localized around ∆k = 0. For tilted-pulse-fronts obtained from diffraction gratings, the latter condition is well satisfied in comparison to the case of beamlet superposition. However, qualitative behavior can still be explained by closed-form expressions in this limit. For quantitative predictions on efficiency and field, one would do better by turning to directly evaluating Eq.7. In the remainder of this section, we illustrate the general physics of this system, using the derived closed-form expressions.
Following the procedure outlined in the appendix, we obtain the following equation for the terahertz spectrum E T Hz (Ω, x, z) in Ω − x space.
In Eq.9, only the role of the final two terms have not been already addressed. The first corresponds to the polarization term which drives terahertz generation between beamlets m, k. The polarization is centered about the average position x m,k . The second term corresponds to the propagating wave that is generated by these beamlets. In contrast to the polarization, this term is centered about x − ztanγ − x m,k to signify that the terahertz wave propagates at an angle γ with respect to the pump. As it propagates, it is attenuated as evident in the exponentially varying loss term e −αz/2cosγ . Furthermore, the terahertz beam size increases due to terahertz diffraction as evident by the factor 16z 2 /(σ 4 (z)k 2 (Ω)cos 2 γ. Evident from the inverse dependence on terahertz wavenumber in this enhancement factor,is the fact that larger terahertz frequencies diffract lesser. Accompanying the diffraction is a reduction in strength of the beam by an amount specified by the factor 1/(1 − 4jz/kcosγσ 2 (z)) 1/2 before the propagation term.
The fast diffraction of terahertz beams implies that only terahertz radiated by beamlet m − 1 superposes on top of that produced by beamlet m (See Fig.5(b)).

Temporal profiles
In order to obtain the spatio-temporal field profiles, we need to calculate the inverse Fourier transform in the time-domain of Eq.9. However,the presence of frequency dependent terms in the denominator of Eq.9 makes this challenging. To eliminate these dependencies while retaining essential features, we set k(Ω)cosγ/2 ≈ σ before performing the inverse Fourier operation. This signifies the fact that the terahertz wavelengths are comparable to the size of individual beamlets. The resulting equation is presented in Eqs.10a-10c.
The spatio-temporal profiles of the generated terahertz pulses are evident from the above expressions. Firstly, there a series of terahertz pulses centered about the mean positions of beamlets m, k, with beam size σ(z)/ √ 2 and pulse duration τ 1 (z)/ √ 2. The √ 2 factor arises from the nature of second-order interactions. The z dependence of pulse duration arises from the balance between material dispersion and any external group-delay dispersion that may have been imparted to the pulse. Each waveform is a summation of a near-field (nonlinear polarization term) and a far-field (or propagating term) [27,28,29]. The field is of the form ∼ te −2t 2 /τ 2 , where t is time in the pulse frame. The near-field term is roughly z propagating as can be seen by inspecting the t term in Eq.10b. From Eq.10b, one can also discern that points along the radius of curvature of the beamlets are located at the same phase. Examining Eq.10c, one notices from the Figure 6: Temporal evolution of terahertz transients obtained from Eqs.10a-10c are plotted (a)A discrete set of beamlets produces a series of transients (b) These grow in amplitude and diffraction of terahertz waves makes them gradually fall into the form of (c) smooth tilted-pulse fronts. Diffraction at the exit surface of the crystal would further abet in washing out the discreteness of the tilted-pulse-front.
second term that close to the source location,i.e. z = 0, the field is virtually z propagating (the second term vanishes). However, within a distance of z ∼ σ, diffraction effects produce a propagating wave predominantly moving at an angle γ w.r.t to the pump beamlets. As the terahertz wave propagates, its radius of curvature continues to increase as evident from the final term in Eq.10c.
In Fig.10, we plot snapshots of the the normalized (to the peak value) terahertz electric field obtained from Eq.10a for lithium niobate with γ = 62 • , n g = 2.25, n T Hz = 4.73 at T= 100 K. In Fig.10(a), the field is weak and both near-field and far-field terms contain a plane of constant phase along t ∼ t − z/v g . However, as the field propagates, its strength increases and the second term in Eq.10c begins to grow due to terahertz diffraction. Therefore, the field gradually becomes aligned along the tilt-plane , i.e. t = t−xsinγ/v T Hz −z/v g as can be seen in Fig.10(c). Any discreteness in the tilted-pulse-front will further wash out as the beam exits the crystal surface.

Comparisons to TPFs from diffraction gratings
In this section, we present results that demonstrate the advantages of using a super position of beamlets as opposed to tilted-pulse-fronts obtained from diffraction gratings. At the outset, we first present key equations for tilted-pulse-fronts that were derived in a recent analysis [30], which account for spatio-temporal distortions that characterize tilted-pulse-fronts generated by diffraction gratings.
First, the counterpart of Eq.7 for gratings is presented below in Eq.11.
The above equation is generally valid and shall therefore be utilized for quantitative comparisons to beamlet superposition. Certain key differences in relation to Eq.7 are worth mentioning. Firstly, notice that the distribution in transverse momentum k x is centered about k(Ω)βv g and has a much narrower spread in k x space owing to the typically larger beam radius w in relation to σ. Secondly, notice a reduction in the amplitude of E T HZ (Ω, k x , z) by a factor proportional to (1 + k T Ω 2 w 2 /τ 2 ) −1/2 , which shall be proven in a subsequent section to most negatively affect the performance of tilted-pulse-fronts produced by diffraction gratings.
Proceeding along similar lines to that of beamlet superposition, one obtains the following closed form expression for terahertz transients generated by tilted-pulse-fronts from diffraction gratings.
Similar to Eq.10a, two principal terms corresponding to the forced oscillation(or near-field term) and propagating wave (or far-field term) are present. The key discrepancy in relation to Eq.10a however is the presence of spatially dependent durations τ 1 and τ 2 as depicted in Eq.12b. These terms produce a spatially chirped terahertz spectrum, commensurate to the amount of GVD-AD k T and bandwidth τ −1 . From Eq.10c, it is evident that both forced and propagating transients, lie along the line x + ztanγ. This is different from the case of beamlet superposition, where the field gradually forms the desired tilted-pulse-front as illustrated in Fig.10.
To illustrate the key differences in terahertz transients formed by beamlet superposition and diffraction gratings, we plot a representative snapshot of the field for each case for a pump pulse duration τ = 50 fs in Fig.7. A tilt angle of γ = 62 • , corresponding to lithium niobate was assumed. This tilt-angle is characterized by a GVD-AD value of k T = −1 × 10 −23 s 2 m −1 . The beamlet radius was σ = 25 µm, with a separation ∆x = 12.5 µm.
In Fig.7(a), the variation of terahertz duration or frequency along the transverse direction is clearly evident. In contrast, the case of beamlet superposition produces a uniform pulse duration across the transverse direction.

Efficiency and spectra
Having illustrated a salient difference between tilted-pulse-fronts produced by gratings and beamlet superposition in Fig.7, we proceed to highlighting situations when beamlet superposition generates terahertz radiation more efficiently compared to tilted-pulse-fronts produced by diffraction gratings. All calculations in this section are performed without approximations using the expressions from Eqs.7 and 11.
In Fig.8(a), we plot the conversion efficiency ratios η S /η DG of beamlet superposition (S) relative to diffraction gratings (DG) for various values of beamlet radii σ, for a pump duration τ = 500 fs. In accordance with threshold conditions outlined in Eq.8, we set ∆x = σ/2. The beam radius for the grating case is fixed at w = 2mm and the total number of beamlets is varied with ∆x as N = 2w/∆x. The peak pump electric field Figure 8: (a) Efficiency ratios of beamlet superpostion to that produced by diffraction gratings for τ = 500fs. As beam sizes get small enough to supply the necessary bandwidth, efficiencies of beamlet superposition is larger compared to grating based tilted-pulse-fronts. (b) Efficiency ratios for τ = 50fs : The threshold beamlet radius is smaller since a larger transverse momentum spread is required to extract the bandwidth of the incident pulse. The ratios increase with larger GVD-AD, which deteriorates terahertz generation by diffraction gratings.(c) Spectra for various cases of τ = 500 fs : The reducing beamlet radius produces an increase in terahertz frequency. (d) Terahertz spectra for the τ = 50fs case. Beamlet superposition produces higher frequencies in relation to gratings, with particularly improved performance for larger GVD-AD values.
for grating and beamlet superposition cases are maintained constant. We rely on efficiency ratios rather than absolute values since absolute efficiency values are only meaningful when depleted calculations (i.e cascading effects are accounted for) are performed.
The efficiency of terahertz radiation generated by beamlet superposition is actually smaller than that obtained from gratings for σ = 100µm (red). This is because the radius σ is too large to contain a sufficient spread in momentum to efficaciously utilize the pump bandwidth. Only for σ = 50µm (blue), is parity in performance obtained. Eq.8 can be re-written as being σ ≤ v T Hz τ /sinγ since the optimal terahertz angular frequency Ω T Hz = 2τ −1 (Obtained by maximizing Ω 2 e −Ω 2 τ 2 /4 ). As σ is reduced to 25µm, the efficiency ratio, significantly exceeds unity. In contrast to the cases of σ = 100, 50µm, which exhibit monotonic growth of efficiency over propagation length z, for σ = 25µm, the influence of terahertz diffraction begins to play a role(albeit only slightly). The minimum in efficiency is observed close to z = 4 mm, which was the location of the beam waist z 0 . Naturally, the smallest beam size has the most adverse effect on terahertz diffraction.
In Fig.8(b), efficiency ratios are plotted for τ = 50 fs. In this case, the larger bandwidth necessitates a smaller beam size to achieve efficiency parity. For instance, in Fig.8(a), efficiency parity is already achieved for σ = 50 µm but this is not the case in Fig.8(b). The reason for this is that the larger initial bandwidth of τ = 50 fs contributes to higher terahertz frequency components compared to τ = 500 fs in the grating case (see black dashed curves in Figs.8(c),(d)). However, the beamlet radius is not large enough to produce the same high frequency components as can be seen by comparing the terahertz spectra produced by the σ = 50 µm (blue, Fig.8(d)) and grating case (black dashed, Fig.8(d)). However, as the beamlet radius is reduced to 25 µm (black, Fig.8(b)), the efficiency ratio is much larger compared to that of the τ = 500 fs case, since the the effect of GVD-AD is more adverse for larger bandwidths. Larger values of GVD-AD deteriorate the performance of grating based tilted-pulse-fronts even further(green, Fig.8(b)). This is because,larger k T values push the generated terahertz frequency by gratings to a lower value (green dashed, Fig.8(d)).
The effects of terahertz diffraction for the τ = 50 fs case is not evident owing to the fact that a larger Figure 9: Effect of beam radius (a) Efficiency ratios for a larger total beam radius w = 5 mm. The performance of τ = 50 fs for beamlet super position is even further enhanced due to a greater impact of GVD-AD for larger beam sizes in the grating case (b) Terahertz frequencies for larger beam sizes are smaller for the grating case, while beamlet superposition shows scalable performance.
terahertz frequency is generated for the same beamlet size of 25 µm in the case of the τ = 50 fs pump pulse (black, Fig.8(d)) compared to the τ = 500 fs case (black, Fig.8(c)). The reduction in efficiency in the Fig.8(b) for longer values of z is because of diffraction of the optical beam and increased absorption at larger terahertz frequencies.
Thus from Fig.8, we can see that with the right beam size, conversion efficiencies of the beamlet superposition case can significantly exceed that from gratings.In addition, beamlet superposition is more scalable with beam size. Fig.9 plots the efficiency ratios for w = 5 mm for τ = 50 fs and τ = 500 fs. In this case, the efficiency ratios are even larger compared to the corresponding situations for w = 2 mm. While the effect is less pronounced for τ = 500 fs, it is significant for τ = 50 fs. This is because the effect of GVD-AD was seen to be more adverse for larger bandwidths . The average frequency also reduces (blue dashed, Fig.9(b)) compared to the w = 2 mm case (black dashed, Fig.8(d)).
In summary, in this section we showed that upon having beamlet radii (σ) which satisfy the threshold conditions (Eq.8), conversion efficiencies for beamlet super position can be far greater, particularly for larger bandwidths and total beam sizes (w).

Analytic proof of efficiency enhancement
The differences in efficiency trends between tilted-pulse-fronts and beamlet superposition cases may be discerned by approximating the sinc functions in Eqs.7 and 4b by a Dirac-Delta function δ((k x − k(Ω)sinγ)z/2). We then obtain the following closed form expressions for efficiency for the titled-pulse-front case.
In Eq.13a, we consider the limit when group velocity dispersion due to angular dispersion |k T |w/τ 2 1.In this regime, η, shows an optimal duration τ opt = 1.62(|k T |w) 1/2 . For w = 5 mm, this corresponds to ≈ 413 fs, which is not far from earlier predictions. Note that the value increases for larger beam sizes.
However, in the total absence of GVD-AD, the efficiency scales inversely as the pulse duration cubed or τ −3 . In contrast, in the regime where |k T | is very large, the efficiency shows the contrasting trend of reducing with τ as depicted in Eq.(13b). Furthermore, the efficiency drops with w 2 in this high |k T | regime.
Therefore, the above Eqs.13a and 13b vindicate earlier calculations and prove that diffraction gratings deteriorate performance for larger beam sizes and bandwidths.
In contrast, the beamlet superposition case is closer to the first term of the undistorted tilted-pulse-front (TPF) in Eq.(13a), with a minor modification as depicted below in Eq. (14). (14) Equation (14) suggests that the effective pulse duration is given by τ 2 + (σsinγ/v T Hz ) 2 , which increases with larger σ. In essence, this suggests that σ has to be sufficiently small to allow for sufficient momentum spread to generate the necessary terahertz radiation, identical to the threshold condition deduced in Eq.8.
If one takes the ratio of Eq.14 and Eq.13a at the threshold condition for σ given by Eq.8, then one one obtains the following conditions for efficiency parity, corresponding to the two limiting cases for gratings outlined in Eqs.13a,13b respectively.
Thus, we have been able to show quite generally that beamlet superposition offers a much better performance for large bandwidths and beam sizes.

Optimizing Efficiency
Obtaining global optima requires detailed numerical simulations, particularly accounting for cascading effects. However, some general bounds for parameters can be established analytically. Firstly, the pulse duration τ is bound by the constraint of containing sufficient bandwidth to generate the desired frequency according to f T Hz = 1/(πτ ). Secondly, the optimal length will either be constrained by the absorption length or the length over which material dispersion spreads the pulse to a value very different from its transform limit. This is encompassed in Eq.16a. Finally, the beamlet radius σ has to be small enough to provide the necessary angular spread to generate f T Hz but large enough to circumvent limitations of either pump diffraction or terahertz diffraction. For pump diffraction, the Rayleigh length has to be larger than the optimal interaction length outlined in Eq.16a. However, for terahertz diffraction, the Rayleigh range only need be larger than ∆x/sinγ, which is the distance the terahertz radiation propagates before encountering the adjacent beamlet as was illustrated in Fig.5(b). The constraint on beamlet radius is thus presented in Fig.16b.
We plot the bounds for beamlet radius in Fig.10 for cryogenically cooled lithium niobate as a function of f T Hz . In this case, absorption is the primary limitation at all values of f T Hz . The maximum value of σ is given by Eq.8 while assuming ∆x = σ/2, while the minimum is given by Eq.16b. The optimum lies somewhere in the shaded region and as stated previously, will require detailed simulations. At lower terahertz frequencies, terahertz diffraction is the main limitation, while for larger terahertz frequencies, one is limited by diffraction of the pump beam.

Cascading effects
While here we dealt with purely undepleted models, the important finding was that beamlet superposition performs better in relation to grating based tilted-pulse-fronts when beam sizes or bandwidths are larger. In effect, cascading is just an increase in the bandwidth of the pulse and hence the improvement in performance at larger bandwidths signals alleviation of cascading limitations. While initial simultations appear to attest this hypothesis [26], future work will address the ultimate limits of beamlet superposition. Figure 10: Optimal parameters for beamlet superposition for lithium niobate. The optimal pulse duration is related to the central terahertz frequency by f T Hz = (πτ ) −1 . The maximum length is absorption limited and given by L = 2cosγ/α. The upper limit of beamlet radius σ is given by Eq.8, while the lower limits are determined either by terahertz or pump beamlet diffraction as shown in Eq.16b

Conclusion
We were able to show by a combination of analytic and semi-analytic methods, the advantages of beamlet superposition over tilted-pulse-fronts generated by gratings for large beam sizes and bandwidths. Closed form expressions in space and time enabled us to shed light on the general physics of terahertz generation using this approach. Conditions for obtaining parity in efficiency were furnished and bounds for optimal parameters were provided. Since cascading is essentially an increase in bandwidth, this analysis signals alleviation of cascading based limitations using beamlet superposition.

Appendix
The following conventions are adopted for Fourier and inverse Fourier operations in the time-domain. The following conventions are adopted for Fourier and inverse Fourier operations in the spatial domain.
The representation of Parseval's theorem or energy conservation is as follows: The general three dimensional scalar wave equation after Fourier decomposition can be presented as follows: ∇ 2 E T Hz (Ω, k x , k y , z) + k 2 z (Ω)E T Hz (Ω, k x , k y , z) = − Ω 2 ε 0 c 2 P T Hz (Ω, k x , k y , z) The solution to the above equation maybe expressed by the ansatz E T Hz (Ω, k x , k y , z) = A T Hz (Ω, k x , k y , z)e −jkz.z . Using the above ansatz and integrating over z , one obtains the following expression.
A T Hz (Ω, k x , k y , z) = − jΩ 2 χ (2) z 0 P T Hz (Ω, k x , k y , z)e j∆kz+ αz 2cosγ dz 2k z (Ω)c 2 One may deal with the many z dependent terms in the polarization term as follows. If one has an integral of the form e Σifi(z) e j∆kz+ α 2 z dz, it may be reduced to the following form e Σifi(z) e j∆kz+ α 2 z dz = Π i Σ n f i (z) n n! e j∆kz+ α 2 z dz = P (z)e j∆kz+ α 2 z dz = P (z)e j∆kz+ α 2 z α/2 + j∆k − P (z)e j∆kz+ α 2 z (α/2 + j∆k) 2 + P "(z)e j∆kz+ α 2 z (α/2 + j∆k) 3 For typical absorption coefficients and ∆k values , the second terms and beyond are much smaller and can be ignored. Therefore, we have e Σifi(z) e j∆kz+ α 2 z dz = e Σifi(z) e j∆kz+ α 2 z (α/2 + j∆k) −1 . Using the above result, expanding k z is expanded about k(Ω)cosγ(since that is where phase-matched terahertz radiation propagates) and performing an inverse Fourier operation on Eq.21, one obtains Eq.23 under the condition that α ∆k. The first term denotes the nonlinear polarization term and the second term represents the field that it radiates. The radiation term is nothing but a Fresnel integral with modified kernel to account for oblique propagation.