1. Introduction

Over the last decades, ultrafast laser-driven terahertz time domain spectroscopy (THz-TDS) has developed into a mature technique [1]. Both the development of the driving lasers as well as advancements in the THz generation process itself have led to a large variety of available sources, which sparked diverse applications in biology and medicine [2], physics [3,4], quality control [5], as well as chemical sensing and the security sector [6], to only mention a few examples. One area that has seen strong development is the generation of strong-field THz pulses, which nowadays commonly reach several GV/m [7], and have enabled the field of THz nonlinear spectroscopy to emerge. However, this development has often been achieved at the expense of the source repetition rate, i.e. the average power of energetic THz sources has not seen comparable progress. This resulted in unwanted compromises for the experiments making use of these sources such as impractically long measurement times and/or reduced dynamic range (DR).

However, higher average power ultrafast driving lasers based on diode-pumped Yb-doped gain media have recently seen spectacular progress, making high-average power THz generation a recent area of active research [8–10]. In fact, nowadays ultrafast pump powers on the kW-level [11–13] are becoming available, rather than the commonly used Ti:Sapphire technology, which is typically limited to less than 10 W of driving power. Among the many available methods to make use of this new laser technology, optical rectification (OR) in lithium niobate (LN) using tilted pulse fronts (TPFs) [14] is well established for efficient high-field THz generation in the 1 THz range [15] and offers many advantages for average power scaling: a high nonlinear coefficient, a high damage threshold and a low impact of multi-photon absorption effects [16]. On the other hand, the pump pulse front needs to be tilted in order to achieve velocity matching between the pump and the THz pulses [17]. The THz radiation is emitted perpendicularly to the pulse front, which results in a large non-collinear angle. This complicates the experimental setup and the THz generation process as compared to collinear generation schemes. Nevertheless, high efficiencies on the order of a percent have been achieved with this technique, making it attractive for both energy and average power scaling. Using a 1 kHz Yb:KYW regenerative amplifier with 1.2 mJ pulse energy, an efficiency on the percent-level was obtained at room temperature, and up to 3.7 % using a cryogenically cooled LN crystal. Focusing on highest energies, 436 µJ THz pulse energy were obtained with a 10 Hz Yb:YAG amplifier system and 58 mJ pump pulse energy [18]. In another experiment 200 µJ were generated employing a 10 Hz Ti:Sapphire amplifier system delivering 70 mJ pump pulse energy [19]. Both experiments also achieved efficiencies on the percent-level. With respect to average power scaling of OR in LN using TPFs, Kramer et al. [9] recently reported on the use of an Yb:YAG innoslab amplifier with a repetition rate of 100 kHz and 3.4 mJ pulse energy (i.e. 340 W of average power) to generate 144 mW at 0.6 THz, which is currently the highest average power of a single-cycle laser-driven THz source. Whereas this is an impressive realization, proving the potential of OR in LN using TPFs for average power scaling, the nonlinear generation mechanism is unaltered compared to the usual mJ-class excitation except for potential thermal effects due to the high average power, which have not been explored in detail. In fact, the excitation regime using multi-mJ pulse energies and above is well explored experimentally as illustrated by the above-mentioned state-of-the-art and in theory [20,21].

In contrast, at much higher repetition rates of several MHz, the currently available average power levels still result in rather moderate pulse energies in the few to tens of µJ, setting additional unexplored constraints to this generation mechanism. Only few attempts have been made to use the TPF method in LN with such high (MHz) repetition rates and correspondingly lower pulse energies. Early on, a 1 MHz Yb-fiber amplifier system with 14 µJ pump pulse energy was used to generate 250 µW of THz average power (250 pJ) with a moderate efficiency of 2.8×10^-5 [22]. More recently, we have demonstrated the generation of 66 mW THz average power (5 nJ), driven by an amplifier-free compressed mode-locked thin-disk oscillator with a repetition rate of 13.3 MHz and up to 8.9 µJ pulse energy [10]. While this is currently the highest laser-driven THz average power at MHz repetition rates, the conversion efficiency was moderate, with only 5.6×10^-4. Furthermore, under optimal conditions, the highest pump peak intensity only reached 8 GW cm⁻², while typical intensities can reach 100 GW cm⁻² when higher pump pulse energies are used [15].

In our experiment [10], we attributed this reduction in achievable conversion efficiency to spatial walk-off due to the strongly non-collinear geometry and our small pump beam sizes, which were required to reach sufficient peak intensity in combination with our low pulse energies. This empirical observation was not confirmed with dedicated simulations, which prove critical in this rather complex nonlinear generation scheme. The effect of spatial walk-off and beam size on the THz conversion efficiency has been discussed before briefly [23–26] , but not investigated in detail for our particular excitation regime of small pulse energies and beam sizes. Furthermore, details of the THz generation process as well as the corresponding scaling laws are so far not understood. In this paper, we employ a 2+1-D model to study the details of the generation process under these circumstances. We reproduce our previous experimental findings and show that indeed, a combination of spatial walk-off and pump beam depletion is responsible for this limitation in conversion efficiency. We discuss strategies to overcome these limitations and predict that watt-level THz sources with MHz repetition rates will become available in the near future. Such sources will open the door for linear and nonlinear THz-TDS experiments with high signal-to-noise ratio (SNR) and DR at MHz repetition rates.

2. Basic considerations and previous results

We consider the traditional geometry for THz generation [17] in the TPF geometry, consisting of a diffraction grating and an aspherical lens, which produces an image of the grating in the LN crystal as depicted in Fig. 1(a). Note, that the calculations performed below do not depend on the type of imaging and would also be valid for a telescopic setup [27]. We choose this configuration, as the main target of our work was to reproduce trends in our experimental setup and draw conclusions about further steps to upscale such sources to higher power. The THz radiation is emitted from the front surface of the trapezoidal LN crystal, perpendicularly to the pulse front. Details of this noncollinear geometry and the used coordinate systems are shown in Fig. 1(b). The pump beam propagates along the $z$-direction, with an $\textrm {e}^{-2}$ beam diameter $2w_0$ along the transverse $x$-direction and is positioned at a distance $\Delta x_0 = h-x_0$ away from the crystal tip. Note that the beam size in the crystal is smaller by the demagnification factor $M=0.66$ compared to the beam size at the position of the grating. The image plane, where the pulse is expected to be shortest, is placed at a distance $z=s$ away from the crystal entrance surface. The $y$-direction is neglected for our purposes. The THz radiation on the other hand propagates along the $z^\prime$-direction, which is rotated counter-clockwise by the pulse front tilt (PFT) angle of $\gamma = 63.4^{\circ}.$

 

Fig. 1. Layout of the THz generation geometry. (a) TPF setup consisting of diffraction grating, aspherical lens and LN crystal. (b) Close-up of the generation crystal including the used coordinate system.

Download Full Size | PDF

Figure 1(b) provides an intuitive understanding of when spatial walk-off becomes significant. Efficient THz generation only takes place over an effective generation length $L_{\textrm {eff}}$, which is determined by the impact of group velocity dispersion due to angular dispersion (GVD-AD) on the pulse duration [17] as well as spatio-temporal break-up of the pump pulse, and typically is on the order of a few mm [20]. THz radiation will therefore be generated in the red, shaded area $A_{\textrm {gen}}.$ Due to the non-collinear geometry, the THz beam then propagates through unpumped parts of the crystal, depicted as the orange, shaded area. This will compromise the conversion efficiency, because (a) the THz radiation is not amplified further, and (b) a significant fraction of the radiation will be absorbed due to the high THz absorption coefficient of LN [28,29]. This walk-off effect will be increasingly significant if the beam size is considerably smaller than the effective generation length $L_{\textrm {eff}}.$ This is typically not the case for experiments with mJ-level pump pulse energies and above, with typical beam diameters of several mm. However, for µJ pulse energies, such as in our experimental study [10] with beam diameters on the order of 1 mm and below, this effect can be expected to be significant.

2.1 Summary of experimental results

Here, we briefly review our previous experimental observations [10], since the numerical results obtained below will be largely discussed in the context of these findings. Figure 2(a) shows the THz average power measured as a function of the pump average power for different pulse durations and beam sizes. The most important parameters used in these measurements are given in the table in Fig. 2(b). Our findings can be summarized as follows:

We will now proceed to reproduce these observations using our 2+1-D simulation tool and shed more light on the origin of the observed efficiency limitations.

 

Fig. 2. Summary of our previous experimental results published in [10]. (a) THz average power as a function of pump power for different beam diameters on the grating and for a pulse duration of 550 fs (top) and for different pulse durations at optimized beam diameters (bottom). (b) Summary of several used and obtained parameters in the experiment.

Download Full Size | PDF

2.2 Numerical model

In order to investigate effects related to such small beam sizes, we numerically solve the coupled wave equations for OR using a 2+1-D model ($x$, $z$, and $t$-dimensions) similarly as in [21]. Details about the model can be found in Supplement 1. The model takes into account pump beam depletion, the Kerr effect, THz absorption, phase matching, walk-off, and diffraction effects. The reflection of the pump beam on the front surface of the LN crystal and possible interference effects are not taken into account. Effects such as multi-photon absorption (MPA) or stimulated Raman scattering (SRS) can be added in a straight-forward fashion, but are neglected here due to our comparably low intensities.

The results for a single simulation run after propagation through the crystal are depicted in Fig. 3. The used parameters match the experimental ones of our 66 mW result [10] assuming Gaussian spatial and temporal envelopes for the pump beam. The pump beam was positioned similarly as in the experiment at $\Delta x_0 = {2.5}\;\textrm{mm}$ away from the crystal tip and the image plane was located at $z=s= {4}\;\textrm{mm}.$

 

Fig. 3. Results for a single simulation run of the 2+1-D code for a 236 fs pulse with both Gaussian spatial and temporal profiles. A pump pulse energy of 8 µJ and a pump beam diameter of $0.66\times {1.4}\;\textrm{mm}$ was used. The crystal was positioned such that $\Delta x_0= {2.5}\;\textrm{mm}$ and $s= {4}\;\textrm{mm}.$ The frame of reference was moving at the group velocity of the pump (co-moving with the pump pulse) and tilted with the pulse front. (a) Pump intensity profile in real space. (b) THz electric field on the crystal exit surface. The data has been multiplied with the Fresnel coefficient and therefore represents the field in free space. (c) Spectral intensity of the pump. $\Delta \omega$ is the radial frequency offset from the central frequency. (d) Spectral intensity profile of the THz pulse with absolute radial frequency $\Omega .$

Download Full Size | PDF

Figure 3(a) shows the spatio-temporal intensity profile of the depleted pump beam. The frame of reference has been chosen to be co-moving with the pump and tilted with the pulse front, i.e. an undepleted pump beam without spatio-temporal couplings would appear as a straight vertical line in this plot. The pump beam does not exhibit a strong break-up behavior [21], due to the relatively low peak intensity of ~8 GW cm⁻². Depletion manifests itself in terms of a shift of pulse energy towards later times in the upper part of the beam. This is a direct signature of the spatial walk-off effect. Since the THz radiation propagates in the positive $x$-direction, the overlap with the upper part of the pump beam is better and therefore causes stronger depletion. Likewise, this effect can be seen in the spectral domain shown in Fig. 3(c). Here, depletion causes a shift towards lower frequencies due to cascading [20]. Again, this effect is more pronounced in the upper part of the pump beam as a result of better overlap with the generated THz radiation.

The spatio-temporal electric field distribution of the THz radiation on the LN crystal surface is depicted in Fig. 3(b). Only few spatio-temporal distortions are observed here, essentially resulting in a single-cycle THz pulse. The corresponding intensity spectrum centered around ~0.7 THz is depicted in Fig. 3(d), also showing only few spatio-spectral distortions.

The conversion efficiency is calculated by dividing the final THz pulse energy by the initial pump pulse energy. Fresnel losses at the crystal surface are accounted for by multiplication with the intensity transmission coefficient ($\approx 0.56$ at 0.5 THz). Furthermore, intensity averaging over the missing $y$-direction is approximated by multiplying with a factor of $ {1}/{\sqrt {2}}.$ In the configuration depicted here, we obtain an efficiency of 12×10^-4, which is a factor 2 higher than observed in the experiment. This can be explained with the more complex pulse envelope shape due to the used pulse compressor (see section 3.4).

3. Results and discussion

Having established the model, we will now proceed to investigate several properties of THz generation in the TPF geometry in our specific excitation regime of small (µJ) pulse energies and beam sizes below 1 mm.

3.1 Propagation dynamics

First, we study the evolution of the pump and THz beams as a function of the propagation distance $z$ inside the crystal. Figure 4(a) shows the obtained conversion efficiency as a function of the propagation coordinate for the three different pulse durations we had available in our experiment. As we can see, for shorter pulses, THz generation takes place within a relatively short distance, while for longer pulses it is spread out over a longer distance, i.e., the effective generation length $L_{\textrm {eff}}$ is bigger. At $z= {4.9}\;\textrm{mm},$ the pump beam is reflected off the crystal surface, which results in the rapid saturation of the efficiency around this value, as the THz radiation is coupled out of the crystal. The change in the effective generation length for different pulse durations can be understood when looking at the corresponding peak intensity change, displayed in Fig. 4(b). The dashed lines show the nominal peak intensity under the undepleted pump approximation. For longer pulses, the effect of GVD-AD is small, such that the pulse duration changes weakly as a function of $z$ and so does the peak intensity. Therefore, the overlap between THz and pump persists over a long propagation distance leading to strong depletion, which eventually destroys the temporal structure of the pulse and causes a drop in intensity. For intermediate values of $z,$ the spectral broadening due to cascading can constitute a self-compression mechanism, that can increase the intensity above its nominal value. This effect is observed for the 550 fs pulse in Fig. 4(b). For 236 fs, where the effect of GVD-AD is stronger, the nominal peak intensity is never reached. The same is true for the 97 fs pulse, but here the deviation as compared with the undepleted case is only small. The reason is, that the effect of GVD-AD is very strong due to the large bandwidth of the pulse. This leads to considerable intensity values only in the vicinity of the image plane at $s= {4}\;\textrm{mm},$ where the pulse is shortest, giving the pump pulse less time to deplete significantly. It becomes clear from these observations that the optimal pulse duration for highest THz generation efficiency depends on a complex interplay between several parameters and cannot be simply extrapolated from other excitation geometries; i.e., for example those observed at higher energies and larger spot sizes, further justifying the need for detailed numerical investigations.

 

Fig. 4. Simulation of the evolution of the conversion efficiency (a) and pump peak intensity (b) as a function of propagation direction $z$ for different pulse durations. The dashed lines show the same calculation in the undepleted pump approximation. Gaussian temporal profiles and a nominal peak intensity of 8 GW cm⁻², as well as a 1 mm beam diameter on the grating were assumed in all cases. A pump beam position with $\Delta x_0 = {2.5}\;\textrm{mm}$ and $s= {4}\;\textrm{mm}$ was used.

Download Full Size | PDF

3.2 Influence of pump beam positioning

We have seen above, that the THz generation dynamics can vary significantly for different pulse durations. As clearly illustrated in Fig. 4, since different pulse durations result in very different propagation dynamics along the effective length of the crystal, the position and effective length over which THz radiation is emitted efficiently varies largely. This illustrates the importance of the geometry and alignment in such setups. In fact, the position of the pump beam on the crystal $\Delta x_0$ as well as the position of the image plane $s$ need to be optimized for highest efficiency for each pulse duration, beam geometry and intensity. Figure 5 displays the conversion efficiency as a function of these critical geometrical parameters. Figures 5(a)–5(c) show the calculation for three different pulse durations for the same parameters as in the experiment. Only the pulse energy was reduced from 9 µJ to 8 µJ for the 550 fs pulse to make these cases comparable. The investigation shows that in order to maximize the conversion efficiency, the pump beam has to be positioned with an accuracy on the mm-scale. The reason for this sensitivity to the position is that the point of maximum THz generation needs to be brought close to the crystal surface. Otherwise, the generated THz radiation will be reduced due to a combination of THz absorption, spatial walk-off and a reduction in pump beam intensity resulting from GVD-AD or pump pulse break-up as discussed in the last section. This procedure can also be understood as choosing the "effective" crystal thickness close to the effective generation length for a particular set of pump beam parameters.

 

Fig. 5. Simulated conversion efficiency as a function of the crystal position with respect to the pump beam for different pulse durations and pulse energies. $\Delta x_0 = h-x_0$ is the distance of the beam center from the LN trapezoid tip and $s$ is the distance of the image plane from the crystal entrance at $z=0.$ The white dots indicate positions, that will be referred to in the text. (a)-(c) Pump parameters chosen as available in the experiment: pulse durations of 97 fs, 236 fs, and 550 fs with beam diameters on the grating of 1.4 mm, 1.4 mm, and 1 mm, respectively. A pulse energy of 8 µJ was used in each case. The pulse shapes for the two shorter pulses were obtained from a frequency resolved optical gating (FROG) measurement, while that of the 550 fs pulse was assumed to be Gaussian. (d) Efficiency for a laser with 550 fs pulse duration, a beam diameter of 1 mm on the grating and a $10\times$ higher pulse energy of 80 µJ.

Download Full Size | PDF

For the short 97 fs pulse, the effect of GVD-AD dominates due to the larger bandwidth, causing the pulse duration to change significantly during propagation. For this reason, also the peak intensity changes quickly, such that it is only high over a relatively short distance. Therefore, pulse break-up as a consequence of pump depletion is less important. This results in the straight line formed by the efficiency maximum, which also allows high efficiency for rather high values of $\Delta x_0$, where the pump beam propagates through a thick part of the crystal.

The situation is reversed for the 550 fs pulse, which is less susceptible to GVD-AD so that both the pulse duration and intensity do not vary much during propagation. Therefore, the efficiency is relatively insensitive to the position of the image plane $s,$ but the pulse is more prone to break up, because the intensity is high over a longer distance than for the 97 fs pulse. As a result, the efficiency decreases for high values of $\Delta x_0$, i.e. thicker parts of the crystal, because the pulse starts to break up already at the beginning of the crystal.

The situation for the 236 fs pulse is a combination of both effects. It shows a higher sensitivity to the position of the image plane than the 550 fs pulse, because of its higher bandwidth, but the efficiency also decreases for high values of $\Delta x_0$, because the pulse is more prone to break up than the 97 fs pulse.

Finally, Fig. 5(d) shows the conversion efficiency for the 550 fs pulse with a $10 \times$ higher pulse energy of 80 µJ. Using the same laser technology as in this experiment, 80 µJ is the highest pulse energy demonstrated so far [30], thus representing a realistic possibility for future improvement. The efficiency in this case is found to be equally insensitive to the position of the image plane, but much more sensitive to the position of the laser beam $\Delta x_0$. While the maximum efficiency grows by about a factor of 3.3, the peak shifts to lower values of $\Delta x_0$, where the crystal is thinner. Following the arguments given above, this can again be understood in terms of pump pulse break-up. Since the intensity is higher, more THz radiation is generated over a shorter distance, leading to an increase in cascading and break-up of the pulse with a reduction in peak intensity. Since this happens already close to the crystal entrance surface, the generated THz radiation can only escape the crystal without being absorbed, if the used section of the crystal is thin.

The above considerations show, that for every simulated parameter set (pump beam size, pulse duration, peak intensity, PFT angle,…) an optimization of the beam position and image plane would be necessary. However, a full scan of $\Delta x_0$ and $s$ as shown above takes about 12 h on our GPU (see Supplement 1), which makes such a procedure impractical. For a global optimization of the efficiency, more sophisticated optimization procedures are currently under development. This will allow to determine universal scaling laws in this excitation regime. Nevertheless, the developed model is a powerful tool to evaluate and understand the trends observed in the experiment, and provides indications on how to optimize a given experimental setup which will be done in the next sections.

3.3 Comparison with experimental spectra

Having described basic features of the excitation regime with low pulse energies and beam sizes, we will now compare the results with some of our experimental results in order to verify the developed code. Note that the experiment was not designed to precisely measure the used values of $\Delta x_0$ or the position of the image plane $z=s.$ Consequently, reasonable estimates of $\Delta x_0= {2.5}\;\textrm{mm}$ and $s= {4}\;\textrm{mm}$ are used. Figure 6 shows normalized calculated power spectra for three different pulse durations, compared with the experimental spectra (dashed lines) [10]. The experimental results agree reasonably well with the simulations. While there is a slight disagreement in the exact peak position and bandwidth of each spectrum, the spectral shape and the shift of the peak position to higher frequencies for decreasing pump pulse duration is reproduced. In the 97 fs case, the flat shape of the spectrum was replicated, although in this case the bandwidth is bigger than for the experimental spectrum, followed by a steeper drop-off. This can be explained with the acceptance angle $\delta$ (see Supplement 1) of the parabolic mirrors, which becomes more significant for this bandwidth as the higher frequency components are not phase matched on-axis but under a small angle. The assumption of such an acceptance angle is a rough approximation, which leads to a steeper drop-off of the spectrum and an onset of this drop-off at higher frequencies.

 

Fig. 6. Comparison of normalized simulated (solid lines) and experimental (dashed lines) power spectra for all three available pulse durations at full pump power. For the simulations $\Delta x_0= {2.5}\;\textrm{mm}$ and $s= {4}\;\textrm{mm}$ was used.

Download Full Size | PDF

Another experimental observable is the depleted pump spectrum. Simulated and experimental spectra are compared in Fig. 7 for the 550 fs pulse. Note, that clean experimental spectra for the shorter pump pulses are not available, due to geometric clipping after the LN crystal resulting from a combination of the larger bandwidth and angular dispersion. The simulated spectrum with a peak intensity of $I_0= {8}\;\textrm{GW}\;\textrm{cm}^{-2}$ (purple line) shows a stronger broadening due to cascading than observed in the experiment. Since intensity averaging over the missing $y$-direction is not captured in the simulations, the same simulation is carried out with an averaged intensity $ {I_0}/{\sqrt {2}}$ (green line), which agrees better in terms of the bandwidth. Still, there is a certain discrepancy regarding the shape of the depleted spectra. However, the propagation of the pump beam after being reflected off the front crystal surface is not taken into account in the simulations, which can be responsible for a change of the depleted pump spectrum.

 

Fig. 7. Simulated (solid lines) and experimental (dashed line) depleted pump spectra for the 550 fs pulse. The peak intensity is $I_0= {8}\;\textrm{GW}\;\textrm{cm}^{-2}$ for the purple lines. The green line shows the spectrum for a peak intensity reduced to $ {I_0}/{\sqrt {2}}$ in order to approximate intensity averaging along the missing $y$-direction. The spectrum of the Gaussian 550 fs input pulse (simulation) is shown for comparison in grey.

Download Full Size | PDF

3.4 Scaling with beam size

One critical observation made during our experimental campaign was that in all explored configurations, the conversion efficiency saturates for small beam sizes rendering the optimal peak intensity clamped at only 8 GW cm⁻² with small efficiencies on the 1×10^-4 scale, thus not allowing us to achieve state-of-the-art percent-level efficiencies as achieved in other experimental realizations. Taking into account the previous considerations, we use our numerical tool to attempt to reproduce this trend. First, we consider the case of the 236 fs pulse duration and calculate the efficiency as a function of beam size for our experimental conditions. We take into account different effects, allowing us to visualize which effects are mostly responsible for the observed trend. The result is presented in Fig. 8(a) which depicts the calculated efficiency as a function of the beam diameter on the grating. Cases (i-iv) depict the different effects included.

 

Fig. 8. Simulated conversion efficiency as a function of the beam diameter on the grating. A crystal geometry of $\Delta x_0= {2.5}\;\textrm{mm}$ and $s= {4}\;\textrm{mm}$ was assumed. (a) Influence of different effects on the conversion efficiency for the 236 fs pulse, 8 µJ and 1.4 mm beam diameter on the grating. (i) undepleted pump approximation for a Gaussian temporal profile. (ii) depletion with Gaussian temporal profile. (iii) depletion plus the actual temporal profile of the compressed pulse obtained from FROG. (iv) all former effects plus a wavefront radius of curvature (ROC) of −50 mm. (b) Full simulation taking into account all effects for 97 fs, 236 fs, and 550 fs pulse durations with beam diameters (on the grating) of 1.4 mm, 1.4 mm, and 1 mm, respectively.

Download Full Size | PDF

For case (i), a Gaussian temporal pulse shape was considered under the undepleted pump approximation. As expected, a monotonic increase in efficiency is observed for decreasing beam size (increasing peak intensity).

Case (ii) takes into account depletion with a Gaussian temporal pulse shape. This leads to a drastically different behavior of the efficiency, which now shows a maximum near beam diameters of 1 mm and thus qualitatively reproduces the experimental trend, albeit at higher efficiency. The maximum is observed, because a reduction of the beam size increases the peak intensity, which then leads to a stronger pump pulse break-up. When the beam size is reduced further, the efficiency drops, because the relative contribution of THz absorption increases as a consequence of walk-off. This trend is also consistent with recent calculations based on a $3$D model [25].

Case (iii) uses the actual pulse profile (in amplitude and phase) after our pulse compressor, which we obtained from a FROG measurement. This maintains the maximum of the efficiency curve below 1 mm, but causes a drop to ~6×10^-4, which is the efficiency level observed in the experiment. The reason for this finding is twofold. On the one hand the peak power for the compressed pulse is about 20 % lower than for an ideal Gaussian pulse of this width. On the other hand its spectrum is broader and more structured, which increases the influence of GVD-AD and pulse break-up.

Case (iv) additionally includes a wavefront ROC of $-50\,\textrm {mm}$. The ROC causes a small change in efficiency, which will be further discussed in section 3.5. The Kerr effect was found to have a negligible effect on the efficiency and is therefore not shown here separately.

Full simulation results taking into account all effects for different pulse durations are presented in Fig. 8(b) where conversion efficiency as a function of beam diameter on the grating is shown. The observed maximum for beam diameters in the vicinity of 1 mm are reproduced in all cases. Moreover, a slight advantage of using 236 fs over 550 fs pulses is also reproduced. A strong increase is predicted for the 97 fs pulse. Note, that this difference of about a factor three, is due to the specific choice of the values for $\Delta x_0$ and $s.$ For fully optimized conditions, only a factor of ~1.5 is expected (see Fig. 5). Nevertheless, this does not agree with the experimental results, where a strong drop to about 1×10^-4 conversion efficiency was observed.

There are several reasons that could explain these discrepancies, such as imaging errors or higher-order angular dispersion effects, that are currently not taken into account in the simulation. Nevertheless, the simulation shows, that in principle short pulses could be advantageous at low pulse energies, to increase the efficiency with higher peak intensities without having to compromise on the beam size. This will be investigated further in future work.

A remarkable conclusion that can be made from these confirmed experimental trends is that the reduction in efficiency at small pump spots cannot be directly attributed to potential thermal effects from operating at small pump spot sizes and high average power, that could also have a similar effect. The fact that the simulations, that do not consider thermal effects, reproduce the observed trend reasonably well seems to suggest that thermal effects are not the main limitation here — however a combined effect of both geometrical and thermal effects cannot be excluded at this point. In a further exploration which is out of the scope of this publication, we are performing specific experiments to disentangle the influence of geometrical and thermal aspects in this generation geometry.

3.5 Influence of phase matching and phase front curvature

In Fig. 8(a) it was shown that a phase front ROC $R$ can lead to a slight change of the efficiency. In this section we will discuss the origin and significance of this effect. The occurrence of small ROC arises from our specific setup and the need to generate small spot sizes. A value of $R=- {50}\;\textrm{mm}$ was obtained by tracing a Gaussian beam through our imaging setup using ABCD matrices. The ROC is typically negligible in more conventional setups with higher pulse energies and large spot sizes. We will therefore further investigate its influence for our excitation regime in the following.

A wavefront ROC in combination with spatial chirp $\zeta$ changes the PFT, which for a beam with Gaussian spatial and temporal profiles can be calculated according to [23,31]

The first term describes the nominal PFT due to angular dispersion $\beta$ for a pump pulse with wave number $k_0$ at the central frequency $\omega _0$ and group velocity $v_{\mathrm {g}}.$ The two last terms cause a deviation from the nominal PFT angle as a result of spatial chirp, which in combination with a ROC $R$ and group dispersion delay $D$ causes an $x$-dependent pulse delay that translates to a change in the tilt angle (for details see Supplement 1). Here, $w_0$ is the beam radius and $t_{\mathrm {p}}$ is the transform limited pulse duration. The results for our experimental situation are shown in Fig. 9(a) for $R= {-50}\;\textrm{mm},$ which is the minimum value for our experiment, and $R=\infty .$ The curvature of $\gamma$ is caused by the second term in Eq. (1), while the third term causes a tilt for finite ROC.

 

Fig. 9. Influence of PFT and ROC on conversion efficiency for different pulse durations. (a) Calculated change of PFT angle as a function of propagation distance $z$ for a phase front ROC of −50 mm. The transverse positon is $x=0$ and the beam diameter $ {0.66}\times {1}\;\textrm{mm}.$ The dashed lines show the tilt angle for infinite ROC. (b) Efficiency $\eta$ as a function of the PFT angle $\gamma$. Pump beam parameters are the same as used in Fig. 8(b). The coherence length of LN for a frequency of 0.5 THz is plotted as the dashed line on the second $y$-axis for comparison.

Download Full Size | PDF

In order to estimate the significance of this effect, it is necessary to know how sensitive the efficiency is with respect to the PFT angle. For all simulations above, the PFT angle was set to its nominal value of $\gamma = {63.4}^{\circ}.$ In Fig. 9(b) the efficiency was simulated as a function of $\gamma$ for different pulse durations. It exhibits a plateau over a range of ~1° before dropping off steeply. This can be understood when considering the coherence length for the OR process [32]. If the coherence length is larger than the effective generation length in LN (~5 mm for the case of 550 fs pulse duration, see Fig. 4(a)), phase matching is less relevant, leading to a plateau in the efficiency. For angles, where the coherence length becomes comparable to the effective generation length, phase matching is relevant and the efficiency drops.

For the 97 fs pulse, the distribution is slightly shifted to higher angles. This is a consequence of the increased THz bandwidth and higher velocity matching angles for larger frequencies. The width of these distributions is on the same order as the angle change calculated in Fig. 9(a) of about 1° over the relevant propagation length (cf. Figure 4), which explains the small, but noticeable impact of the ROC on the conversion efficiency. Based on these results, we can conclude that the phase front curvature in our experimental setup did not pose a major limitation for the parameter space investigated in this work. However, one can expect a stronger negative impact on the overall conversion efficiency for smaller values than −50 mm.

3.6 Discussion on absolute efficiency values

The simulated values for the conversion efficiency of ~6×10^-4 are in good agreement with our experimental findings for the 236 fs and 550 fs pulses. Nevertheless, these values should be treated with care as there is a number of uncertainties in the numerical model compared to the real experimental conditions. We have already discussed the influence of the pump beam position, which could not be measured with high precision in the experiment, thus resulting in a certain degree of uncertainty. Furthermore, it remains unknown whether imaging errors play a role in our setup. Ray tracing simulations of the experimental geometry will provide additional elements to this question, and will be explored in further work.

As we mentioned in section 3.4, an important aspect that these simulations do not allow to elucidate is whether thermal effects play a role in the observed THz generation efficiency, and what the interplay between thermal effects and the effects simulated here is. While we have not observed strong temperature changes in the crystal with a thermal camera or seen any clear indications of thermal damage, further investigations are needed to rule thermal effects out entirely. In future work, we will target to disentangle thermal effects from generation effects, which can be done for example by varying the repetition rate at constant pulse energy for different crystal temperatures. In general, it is expected that heating (resp. cooling, e.g., with a cryostat) of the crystal would lead to an increase (resp. decrease) of the THz absorption coefficient and thereby generally negatively (resp. positively) influence the conversion efficiency. However, as analyzed in detail above, the influence of THz absorption in this geometry and particularly with small beam sizes strongly depends on alignment and input parameters. Therefore, it becomes clear that a general statement on what quantitative benefit would result from a decrease in THz absorption (for example via cooling) is difficult to make. In this context it should also be noted that reported values for the absorption coefficient vary in the literature by about a factor of two [28,29]. Therefore, the exact amount of THz absorption is subject to a significant error bar for our samples.

It is clear from the discussion above and from the general complexity of the generation mechanism that the absolute efficiency values obtained by the simulations in this work should be taken as an estimate for the order of magnitude. Nevertheless, the good agreement with the experimental results shows that the developed model can be used to predict the relative scaling of the conversion efficiency with the parameters of the pump laser.

3.7 Summary of simulation results

This section briefly summarizes the findings we obtained with the 2+1-D simulation tool.

4. Outlook: scaling to higher THz powers and pulse energies

The simulations presented in this work clearly show that small beam sizes represent a limiting factor for the conversion efficiency, thus setting constraints on the minimum pulse energy that makes sense to apply using this generation geometry. In our experiments, this trend was observed in that reducing the beam size in order to increase the peak intensity beyond a certain moderate value did not increase the obtained THz power further, an effect that was consistently observed in several scenarios. Reducing the pulse duration to increase the peak intensity only mitigated this effect moderately, and as we highlighted above, shorter pulses increase the influence of other effects that adversely affect the generation efficiency and is therefore not a straightforward path to significantly improve the current results. The most straightforward approach to circumvent the current limitations is to increase the intensity by using higher pulse energies. This is illustrated in Fig. 10 which shows the evolution of the efficiency as a function of the driving pulse energy for two different pump beam positions for the 550 fs pulse duration case. Compared to our experimental situation with ~9 µJ, we can expect an increase of the efficiency by about a factor of four, when using a laser system with 80 µJ pulse energy. This $9\times$ increase in pulse energy is within reach of mode-locked thin-disk technology, albeit at the moment at a slightly reduced repetition rate of 3 MHz [30]. Accounting for the increased average power as well as a factor of four in efficiency and considering our present THz average power of 66 mW, such a driving laser should provide a record THz power of ~500 mW in the near future. Another potential route, without sacrificing repetition rate, is to increase the driving average power at a constant, or even higher pulse energy; an attractive set of parameters could be 100 µJ at 10 MHz repetition rate, resulting in 1 kW of average power. Whereas this parameter range has not yet been directly achieved with oscillators, it is achievable with state-of-the-art amplified systems [11,13]. This will however certainly require a more thorough understanding of thermal effects in this geometry.

 

Fig. 10. THz conversion efficiency in LN with high pump pulse energies for two different positions as indicated in Fig. 5(d). The vertical dashed lines mark the pulse energy of the present driving laser and a possible future driving laser with an increased pulse energy of 80 µJ. The top $x$-axis shows the average power of this laser with a repetition rate of 3 MHz.

Download Full Size | PDF

In our best experimental realization so far, we reach an electric field strength of 16.7 kV cm⁻¹, which can be optimized using tighter focusing by at least a factor of two. Using a laser system such as the one mentioned above with 80 µJ and 3 MHz repetition rate, we can expect an increase in the achievable electric field strength by about a factor of 5 to 6. Depending on the THz imaging configuration, this would lead to a maximum field strength on the order of 100 kV cm⁻¹ to 200 kV cm⁻¹, which is well in the range where nonlinear THz spectroscopy experiments can be considered [34].

5. Conclusion

In this work, we presented a thorough numerical investigation of THz generation with TPFs in LN for the unusual excitation regime of small pulse energies on the order of 10 µJ and small pump beam sizes below 1 mm, such as those obtained by state-of-the-art high average power MHz repetition rate ultrafast Yb-lasers. This regime has become increasingly relevant for the THz community, where more and more applications call for highest possible repetition rates without sacrificing pulse strength. We employed a 2+1-D pulse propagation model to simulate the THz generation process for our previous experiments, in which we demonstrated record high powers at 13.3 MHz repetition rate [10], but also observed unexplained limitations in conversion efficiency. Our advanced simulation tool was able to reproduce most observed experimental trends, and in particular shed light on the conversion efficiency clamping observed when using small beam sizes. We conclude that the currently observed limitation is a direct consequence of the excitation geometry—i.e., a combination of spatial walk-off and spatio-temporal break-up of the pump pulse prevents an increase of the peak intensity by reducing the pump beam size—rather than thermal effects. However, further experiments are required to confirm a total absence of thermal effects, which could additionally adversely affect the generation efficiency at small spot sizes and high average powers. A detailed investigation to disentangle these effects is currently in progress.

The simulations show that using comparable average power levels, but scaling the pump pulse energy towards 100 µJ will significantly mitigate the observed limitations, which is possible with comapct state-of-the-art thin-disk oscillator technology, but also with Yb- amplifier schemes. Using such a laser, possibly in combination with crystal cooling, we predict average powers approaching the watt-level with MHz repetition rates and electric field strengths exceeding 100 kV cm⁻¹. This will enable novel linear and non-linear experiments at high repetition rates with unprecedented SNR and DR.