Experimental three-dimensional beam profiling and modeling of a terahertz beam generated from a two-color air plasma

Enabled by recent developments in high-intensity table-top terahertz (THz) sources, there is currently a strong upsurge of interest in nonlinear interactions between THz waves and matter. The highest peak intensities are reached with ultrafast, pulsed THz sources. A second-order nonlinearity driven by intense femtosecond laser pulses in appropriate crystals results in difference-frequency generation of intense THz pulses, with the most prominent materials currently in use being the inorganic crystals LiNbO₃ [1, 2], ZnTe [3] and GaSe [4, 5], and the organic salt DAST [6] and related structures such as DASC [7] and OH1 [8]. With these crystals, pulse energies up to tens of μJ and peak field strengths in the MV cm⁻¹ range have been reported. Such field strengths are by far sufficient for nonlinear studies, and a range of fundamental physical phenomena such as impact ionization [9], side valley scattering [10], THz-induced optical Kerr effect in liquids [11], vibrational anharmonicities in molecular crystals [12], the dynamic Franz–Keldysh effect in multiple quantum wells [13] and quantum dots [14], and even THz-induced phase transitions [15] have been investigated in the low (< 3 THz)-frequency range. Nonlinear two-dimensional spectroscopy on bulk InSb [16], dual quantum wells and multilayer graphene [17] have been demonstrated using THz transients in the 10–40 THz range. In principle, the bandwidth of the generated THz pulse can be as high as the bandwidth of the femtosecond pump source. However, in practice, the generated bandwidth is always limited by the phase matching conditions offered by the specific crystal at the pump wavelength.

Nonlinear THz pulse generation in femtosecond two-color air plasmas has recently become an attractive alternative to crystal-based THz pulse generation mechanisms, in part due to the intriguing physical mechanisms involved in the interaction between the intense femtosecond laser pulse and the generated plasma, and in part due to a THz output limited only by the duration and bandwidth of the pump laser pulse [18, 19]. THz emission from a two-color plasma was first reported by Cook and Hochstrasser [20], who identified the THz generation mechanism as four-wave mixing due to the third-order nonlinearity of air. It was then shown that ionization of the air is required for efficient THz generation [21], and the coherence of the four-wave mixing mechanism in the plasma was then reported [22]. A quantum-mechanical treatment of the electron trajectories in the plasma is required for a full understanding of the polarization properties of the emitted THz radiation [23, 24]. Further investigations showed the strong influence of the transient currents in the plasma on the THz generation mechanism [25, 26], and the finite length of the plasma was shown to lead to conical emission of THz radiation in a preferred angle from the plasma due to constructive interference caused by the phase matching conditions between the fundamental and the second harmonic pump beams [27–29]. It has been predicted [30] and recently also experimentally verified [31] that excitation of the plasma with longer-wavelength pump pulses leads to strong enhancement of the generated THz energy, with resulting peak field strengths in the MV cm⁻¹ range, and pulse energies approaching the μJ range.

With offset in the theoretical and experimental findings on the generation mechanism in a two-color plasma discussed above, we investigate the full three-dimensional (3D) beam profile emitted from a two-color femtosecond plasma. In any nonlinear experiment utilizing THz pulses either from two-color plasmas or other sources, it is of paramount importance to have detailed information about the spatial distribution of the THz field in order to estimate the intensity incident on the sample under investigation. The THz pulse energy is always limited by the available pump energy and saturation effects in the generating medium, while the THz pulse duration is limited mainly by phase matching conditions or the duration of the excitation laser pulse itself. The spot size of the focused THz beam, on the other hand, is to a large extent determined by the geometrical arrangement of focusing elements in the experimental setup. The conical emission observed by You et al [29] must be carefully considered in a nonlinear experiment where the spatial distribution of the THz energy must be known, especially if e.g. the z-scan technique is employed [32]. We discuss the observed 3D beam profile in terms of measurements of the Gouy phase shift in the focal region, and simulations based on the transient photocurrent model put forward by You et al [29] show that the focusing of the beam profile is well described by propagation of a Bessel–Gauss (BG) beam.

The setup for THz generation with a laser-induced two-color air plasma [29] is shown in figure 1. The incident laser beam is from a regenerative amplifier system (SpectraPhysics Spitfire) giving ∼ 100 fs, 2.9 mJ pulses centered at 800 nm with a repetition rate of 1 kHz. The 800 nm pulses generate their second harmonic (400 nm) in a 100 μm thick β-barium borate crystal, and both beams are focused collinearly using a biconvex lens with a 300 mm focal length. From this, an air plasma with a length of approximately 20 mm is generated. A photo of the air plasma is shown in figure 1. The THz light emitted from the plasma is collimated and focused by two off-axis parabolic mirrors, OPM1 and OPM2, with effective focal lengths of 4'' and 3'', respectively. A silicon beam splitter transmits the THz light but reflects the residual optical beams so that only the generated THz light is detected on the camera. The frequency range of the THz light generated with this setup spans from ∼0.2 to 12 THz (10% bandwidth), as seen in figure 5(b). Using a translation stage, a THz camera (NEC, model IRV-T0831) is scanned through the focus in steps of 0.2 mm. The total scanning distance is 52.8 mm, and a THz radiation profile image is recorded for each step by integrating over 128 frames (camera video rate: 30 frames s⁻¹). The camera has a lock-in function chopping the incident laser beam at 15 Hz to reduce noise. The images recorded are 240 × 320 pixels with a pixel size of 23.5 μm, resulting in an image size of 5.64 mm × 7.52 mm [33].

Zoom In Zoom Out Reset image size

Two examples of recorded intensity images are shown in figure 2. A representative image of the unfocused beam is shown in figure 2(a) whereas figure 2(b) shows a measurement of the focused beam. In the unfocused intensity profile, we observe the conical emission pattern appearing as a donut-shaped beam profile, which is focused down to a single spot as seen in figure 2(b). We note a slightly uneven intensity distribution of the emission pattern through the whole beam scan, which was also observed in [29] for higher frequencies up to 10 THz. This is present in figure 2(a) as the brighter emission pattern in the left half of the image. Presently, we are unsure if this is the result of slight optical misalignment or some unhypothesized physical mechanism. The generated spectrum (see figure 5(b)) spans up to approximately 12 THz, but the camera is only sensitive from 1 to 7 THz, so we expect the recorded intensity images to be representative for this frequency span. The THz beam is linearly polarized, both in the focal plane and in the unfocused parts of the beam, as determined by insertion of wiregrid polarizers.

Zoom In Zoom Out Reset image size

To investigate the evolution of the intensity distribution over the entire beam scan in a manner which is independent of the beam shape, the standard deviation ($\sigma _x = \sqrt {\int_x {(x - \langle x\rangle )^2 I(x)\,{\rm{d}}x} }$) is calculated for the X (horizontal) and Y (vertical) directions for each recorded image. The calculated standard deviations in the two directions are plotted as a function of the camera stage position in figure 3. The intensity of the beam is distributed in a single spot when stage position is close to the focus position and the standard deviation is minimal in both directions, i.e. the stage position is between 25 and 32 mm. It is seen that the standard deviation for the two directions X and Y do not minimize at the exact same stage position, which we attribute to astigmatism in the imaging system. To get an exact measure of the beam size, the full-width at half-maximum (FWHM) was measured on the images recorded within this region. The minimum FWHM widths were measured to be 127.9 and 140.3 μm for the Y- and X-direction, respectively. It should be noted that only the region of the beam scan where the intensity pattern is a single spot is desired for a potential nonlinear spectroscopy experiment.

Zoom In Zoom Out Reset image size

The series of intensity profile images recorded with the THz camera can be combined into a 3D profile for full visualization of the beam properties. Figure 4 and the associated animation show the measured 3D profile of the THz beam. The intensity profile is normalized at each z position in order to visualize also the weaker portions of the beam. In the animation the beam is rotated around its focus point, and figure 4 shows two screen shots illustrating the 3D beam shape. The extent of the focal region is clearly seen, indicating a very low divergence of the THz beam in this region. The rings in the outer regions of this 3D beam profile clearly show the conical shape of the beam in the far field.

Zoom In Zoom Out Reset image size

The 3D visualization of the focused beam profile reveals several unusual features compared to a standard Gaussian beam that are worthy of a detailed discussion. The first unusual feature is the observation that the intensity profile is close to Lorentzian at the beam waist, and the second feature is the observation of the transition from the far field donut pattern to the spot-like profile in the focus region, with a prolonged, pencil-like focal region which is distinctly different from the focusing of a Gaussian beam. These features will be discussed in the following two sections. Additionally, we observe a small off-axis pencil-like structure reminiscent of the focal region at large z values (see figure 4(b)). This feature is weak (approximately 5%) compared to the focused intensity, but is rendered visible due to the normalization procedure. We believe that this feature is due to a spurious reflection in the optical system.

For most THz applications, and especially for nonlinear spectroscopy, Gaussian-shaped beam profiles are desired in order to achieve localized, well-controlled high electric field values. However, we do not observe a Gaussian beam shape for our measurements in the focal plane, and instead achieve a much better experimental fit using a Lorentzian profile (see the intensity image in figure 2(b) and its horizontal projection in figure 6(a)). For this reason, we perform a second experiment in order to investigate the beam waist at the focal plane.

When focusing a beam as in the experiment described in section 2, the wave experiences an additional phase shift—the so-called Gouy phase shift [34–36]. Kužel et al [37] defined how the Gouy phase shift affects the measurement of the refractive index of a sample placed in the focus of the THz beam compared to a sample placed in the collimated beam. By comparing those two extracted refractive indices for the same sample, the Gouy phase shift can be determined. Assuming a Gaussian beam in the focal plane of the beam, we applied the procedure described by Kužel et al to extract information about the beam waist size over the entire experimental frequency range. The setup used for the Gouy phase experiment is shown in figure 5(a). The THz beam is generated and focused in the same way as described previously, and now the THz spectrum is detected using the air-biased coherent detection (ABCD) technique [38]. A typical spectrum measured with this setup without a sample is shown in the inset of figure 5(b). The Gouy phase shift is extracted by measuring the difference in refractive index from a 525 μm thick high resistivity silicon wafer at two distinct positions, namely in the focal plane of the THz beam, and in the collimated part of the beam, as indicated in figure 5(a).

Zoom In Zoom Out Reset image size

Using these scans, we extract the experimental beam waist for all frequencies, as shown in figure 5(b). Assuming that the intensity profile at the beam waist is a superposition of frequency-dependent Gaussian beams with width w₀ (f) and with amplitudes weighted according to the spectrum in figure 5(b), we obtain the intensity profile displayed in figure 6(b). Comparing this intensity profile derived from the Gouy phase shift with the experimental measurement, we see that both profiles have Lorentzian shapes, and are in good agreement with each other. The low-intensity part of the profile measured with the THz camera is slightly broader than the profile extracted from the Gouy phase measurement. However, the exact frequency response function of the camera is unknown, and the discrepancy could easily be due to additional broadening from the response function of the camera.

Zoom In Zoom Out Reset image size

The anomalous behavior of the beam profile as it transforms from the far field to the beam waist and back suggests that we look at the existing literature for Helmholtz solutions that mimic this behavior. Just such a solution, the BG beam, was derived independently by both Sheppard and Wilson [39] and later Gori et al [40]. The BG beam is a linear superposition of traditional Hermite–Gaussian beams. While higher order solutions exist [41], for simplicity we will only refer to the fundamental mode of the ordinary BG beam. To conceptually describe a BG beam, imagine a traditional Gaussian TEM₀₀ mode propagating at a small angle, ε, from the z-axis. If ε is large relative to the rate of beam waist expansion, w(z), then the Gaussian profile traverses away from the z-axis with increasing distance. The on axis amplitude will go to zero except near z = 0, where the beam profile and the z-axis intersect. To construct the full BG beam, this TEM₀₀ mode is then swept a full 2π radians around the z-axis. This yields the full solution of [39–41]

$$\begin{eqnarray} &&U_{{\rm{BG}}} (r,z) = A_0 \frac{{q_0 }}{{q(z)}}{\rm exp}\,\left( {\frac{{{\rm{i}}k}}{{2q(z)}}(r^2 + r_c^2 (z))} \right)J_0 \left( {\beta r - \frac{{kr_c (z)}}{{q(z)}}r} \right), \nonumber\\\\ &&\nonumber q(z) = q_0 + z,\quad r_c (z) = \left( {\beta /k} \right)z,\quad \varepsilon = \beta /k. \end{eqnarray} \tag{ 1 }$$

In particular, it is also worth examining this solution at the beam waist. At z = 0, the profile is Gaussian with the addition of an interference term, a zeroth-order Bessel function, which arises from the integration of the slight angular deviations from each independent Gaussian profile [42],

$$\begin{equation} U_{{\rm{BG}}} (r,z = 0) = A_0 \,{\rm exp}\,\left( {\frac{{ - r^2 }}{{w_0^2 }}} \right)J_0 \left( {\beta r} \right). \end{equation} \tag{ 2 }$$

To investigate how this beam compares with our experimental data, we implemented the transient photocurrent model popularized by Kim et al [25, 43]. The physical distribution of the induced plasma is approximated as a simple one-dimensional line. In order to maximize the donut-like behavior of the photocurrent model, a plasma length of 2l_d ≈ 40 mm [29] was chosen, where l_d is the dephasing length. This results in a full 2π phase shift between the 400 and 800 nm beams for a realistic electron density of 10¹⁶ cm⁻³. In practice, the experimental plasma length is closer to 20–25 mm. However, the lengths of the plasmas in simulation and experiment are not directly comparable with each other. Changes in the incident laser intensity alter the plasma density, which will in turn alter the phase shift between the 400 and 800 nm beams along the filament, resulting in distinct changes to l_d and the emitted THz radiation profile. A more detailed treatment, and the resulting direct comparisons, could be achieved with a full 3D finite difference time domain (FDTD) simulation.

There are three small differences between [29] and the model implemented here. Firstly, the FWHM of the laser pulse was changed to 100 fs to better match our experimental conditions. Secondly, we used the Keldysh–Faisal–Reiss [44] equation for ionization instead of Ammosov–Delone–Krainov, and thirdly, we calculated the far field electric pattern in the time domain, according to [45]

$$\begin{equation} E(t,\vec r) \approx \frac{1}{{4\pi \varepsilon _0 c^2 }}\int_{V'} {{\rm{d}}\vec r'\frac{{[\dot j(t'_{{\rm{ret}}} ,\vec r) \times (\vec r - \vec r')] \times (\vec r - \vec r')}}{{\left| {\vec r - \vec r'} \right|^3 }}} , \end{equation} \tag{ 3 }$$

where

$$\begin{equation} \dot j(t'_{{\rm{ret}}} ,\vec r) \equiv \left( {\frac{{\partial j}}{{\partial t}}} \right)_{t = t_{{\rm{ret}}} } \end{equation} \tag{ 4 }$$

is the derivative of the current evaluated at the retarded time. Using this model, we calculated the electric field in two far field planes, at z₁ = 0.5 m and z₂ = 0.75 m.

The solution at these two planes allows us to analyze the beam profiles calculated from the transient photocurrent model as a BG beam. The two far field profiles are Fourier transformed in the time domain to obtain single frequency profiles. ε was determined from a simple geometric calculation ε = tan⁻¹((x₂–x₁)/(z₂–z₁)), where x₂ and x₁ are the physical locations of the maximum value of the donut-shaped profile. In this way, the only free parameter remaining is q₀.

Once the BG beam parameters are known, it is possible to propagate the beam parameters through optical systems according the ABCD matrix formalism [41, 46] (not to be confused with the ABCD—air biased coherent detection abbreviation mentioned previously). Here we have applied that formalism to account for OPM1 and OPM2 in figure 1.

In figures 7(a)–(c), we show three single-frequency slices that result from this full treatment. It is clearly seen in (a)–(c) that the beam profiles have the experimentally measured donut profile in the far field, and collapse to a spot-like profile at the beam waist as seen in (d)–(f). However, the beam waist profiles shown in (d)–(f) also display an extra ring, resulting from the Bessel function in equation (2), which is not present in experiment. It is important to recall that the experimental beam measurements themselves are broadband in nature, and that the extra ring radius in the BG beam is frequency dependent. When the broadband beam waist profile is reconstructed, using the full bandwidth of the model, we see that rings of different frequency interfere destructively with each other and gradually disappear with increasing bandwidth, resulting in the smooth, single spot shown in figure 7(g). To demonstrate the effect of added frequency components on the spot size, we show horizontal projections of the beam spot with different bandwidth contributions in figure 7(h). As the bandwidth increases, the wings of the projection are clearly diminished.

Zoom In Zoom Out Reset image size

In the same manner as for the experimentally measured beam profile in the previous section, a broadband (1–7 THz) 3D beam profile is reconstructed in figure 8. The associated animation shows a rotation of the 3D beam profile around its focal plane. In comparison with the measured 3D beam profile in figure 4, we observe the same overall behavior of the simulated beam profile. In both cases the donut-shaped beam is focused down to a narrow single-spot intensity profile with a low divergence around the focal region. It is worth noting that the focal region is not as narrow and extended as it appears. The intensity profile of each z slice is independently normalized so that the far field donut behavior can be seen in the same plot as the focus. This makes the focal region appear highly localized over an extended range.

Zoom In Zoom Out Reset image size

While our combined model of the transient photocurrent and BG beams qualitatively recreates the experimental data, there is an outstanding disagreement concerning the shape and size of the beam waist profile. As seen in the previous sections, we have a Lorentzian profile at the beam waist. The treatment here predicts a Gaussian profile, as seen in figure 7(h). Furthermore, the FWHM from this model is approximately twice that observed in our experiment. We believe that the discrepancy lies in the physical length of the generated plasma chosen for the photocurrent model. Here we have chosen the length to be 2l_d in order to maximize the donut-like behavior and minimize the radiation in the forward direction. In practice, our experimental plasma length is expected to be different from this value, which should result in a combination of donut-like emission as well as a radiation component in the forward direction. The faint shading in the interior of the donut in figure 2(a) is suggestive of this forward radiation, but is not conclusive. As shown in section 4, its presence could have a non-intuitive effect on the beam profile, which is completely neglected in our modeling. Furthermore, the physical width of the plasma filament has also been neglected in our one-dimensional model. This may also explain the discrepancy through interference from different physical locations in the plasma. The relative contribution of these effects is still under investigation.

In conclusion, we have demonstrated that a wealth of important information can be obtained from the full 3D characterization of the beam profile of THz radiation emitted from femtosecond two-color plasmas. While a simple knife-edge scan across the beam profile in the focal plane would have revealed a nicely focused THz spot, the donut-like full beam profile away from the focal region is far from simple, and a nonlinear experiment such as a z-scan will require precise information about the full beam profile in order to be interpreted correctly.

The intriguing combination of the nonlinear optical response of the plasma and the transient photocurrents in the plasma induced by the asymmetry of the optical field results in emission of a conical beam of THz radiation which can be focused back to a central spot, here determined to be of Lorentzian cross section. This cross section is due to a frequency-dependent focused spot size, inferred by considering the frequency-dependent Gouy phase shift in the focal region.

The transformation from a donut-shaped beam profile to the central peak can be explained in terms of the propagation of a BG beam. The combined transient photocurrent generation model and BG propagation model presented here, which qualitatively describes the 3D experimental profiles seen in section 3, is far from conclusive. To fully verify this interpretation, single-frequency profiles must be measured that match the predicted profiles in figure 7, and that validation is still on-going. However, the applicability of the ABCD matrix formalism to BG propagation through optical systems, and its potential utilization in system design, provides strong motivation for this continued study.

We acknowledge financial support from the Danish Council for Independent Research (FTP Project HI-TERA), the Carlsberg Foundation and the H C Ørsted Foundation. ACS thanks Dr Damian N Schimpf for fruitful discussions on BG beams.