Spatio-spectral characteristics of ultra-broadband THz emission from two-colour photoexcited gas plasmas and their impact for nonlinear spectroscopy

We present a characterization of the combined spatial and spectral properties of the terahertz (THz) and mid-infrared emission from gas plasmas generated and driven by two-colour femtosecond optical pulses. For its use in nonlinear spectroscopy, the impact of the relatively complex spatial profile for both broadband (∼ 10 THz) and ultra-broadband (> 100 THz) emission needs to be considered, in particular for experiments based on z-scan techniques. Here we apply spatially resolved measurements based on both field autocorrelation and sum-frequency (up-conversion) detection. Based on these results, we present simulations of the ultra-broadband profile during its passage through a focal region. In addition to the inherent features of the emission profile due to the generation mechanism in the plasma filament, we also analyse the role of the semconductor (silicon) wafer typically placed after the plasma to discard the optical pump beams, whose photoexcitation also can play a role in the resultant THz profile.


Introduction
Among the various high-field table-top-scale terahertz (THz) sources for emerging nonlinear THz spectroscopies, the emission from gas plasmas generated and driven by two-colour femtosecond photoexcitation [1][2][3][4][5][6][7][8][9] has certain attributes which make it attractive for a range of experiments. While optical rectification using a tilted optical pump pulse front in LiNbO 3 can achieve high fields (∼1 MV cm −1 [10]) and very high pulse energies ( 10 µJ [11]), the spectral bandwidth is limited to 3 THz due to phase-matching and phonon absorption. For the nonlinear excitation of spectral features at higher frequencies, an alternative approach is difference-frequency generation in GaSe crystals [12]. In terms of experimental apparatus, such sources require two multi-stage optical parametric amplifiers to provide the input optical pump and signal pulses. Peak fields of >10 MV cm −1 with spectra covering the range from ∼10-60 THz are achieved (the latter determined by the phase-matching conditions in GaSe). For experiments with even shorter THz pulses, two-colour plasma-THz emission offers bandwidths well in excess of 100 THz [9,13], with a continuous bandwidth extending from <1 THz allowing one to realize true single-cycle transients and consider their use for the study of nonequilibrium dynamics on time scales <10 fs. Based on pulse energies approaching 1 µJ [9], peak fields >10 MV cm −1 can be expected for near-transform-limited pulses.
An important requirement for nonlinear spectroscopy is a means to vary the peak field/intensity incident at the sample, in order to characterize the nonlinear response in detail. For the ultra-broadband THz-mid-IR pulses from a plasma emitter, however, this is not a trivial matter. While for second-order nonlinear generation methods, one can simply control the optical pump energy (below saturation thresholds), such a method is not robust for the plasma emitter due to the large degree of self-action of the optical pulses on the plasma medium-in particular, the plasma spatial profile, defocusing and phase modulation [5,8,9] which depend in a highly nonlinear fashion on the pump pulse energy-such that the spectral/temporal properties of the THz emission will also be significantly influenced. Commonly used wire-grid polarizers [14,15] become ineffective above ∼5 THz, while mid-IR polarizers suffer from a low-frequency cutoff and their substrates can introduce unacceptable temporal pulse broadening. In practice, an ultra-broadband THz-mid-IR beam should be guided to the sample with only reflective optics, as no practical means for pulse chirp compensation over such enormous bandwidths are currently available. Moreover, simple means such as placing a variable iris in the THz beam has a significant effect on the THz pulse properties (as will be demonstrated in this paper).
One of the remaining options to control the intensity at the sample is to employ a z-scan technique, i.e. by translating the sample along the beam axis through the focal region, as is commonly employed e.g. for nonlinear measurements in the optical spectral region [16]. In order to interpret the results of such measurements, one must apply certain assumptions concerning the focusing behaviour of the beam (or perform a complete characterization using media with a known response). For super-octave-spanning pulses, such as the ultra-broadband pulses considered here, even for an ideal Gaussian beam profile one already has the complication of a strong spatio-temporal transformation in traversing the focus, due to the effects of a frequency-dependent Gouy phase [17] and a frequency-dependent focal beam radius (whose iω-dependence manifests as a temporal derivative relative to the collimated beam for the onaxis field).
In addition to these effects, for plasma-THz emission it is well established that the beam profile deviates strongly from an ideal Gaussian shape [18,19], which has been reconciled in previous reports using simplified one-dimensional (1D) models where a phase-mismatch between the two-colour optical beat and THz wave during propagation results in off-axis phasematching and a conical emission profile. However, the THz spectral dependence of this beam profile has only been studied at low frequencies ( 2 THz, [18]), or qualitatively using two spectral-range filters (Teflon and Ge [19]).
In this paper, we present spatially and spectrally resolved measurements covering the full spectral range, both for (i) broadband THz emission (bandwidth ∼10 THz) using optical pump pulses with duration (full-width-at-half-maximum) T ω = 150 fs, and (ii) ultra-broadband THz-mid-IR emission (>150 THz) with T ω < 20 fs. The ultra-broadband detection is achieved using both field autocorrelation (Fourier-transform infrared spectrometry (FTIR)) with a thermal detector, as well as femtosecond time-gated detection using optical-THz sum-frequency generation (up-conversion) in a χ (2) -medium. Besides a discussion of the inherent emission profile from the plasma filament, we also analyse the role of the semiconductor wafer (here high-resistivity Si) placed after the plasma to discard the optical pump beams. It will be shown that the unavoidable photoexcitation of charge-carriers in this wafer has a significant effect on the THz radiation close to the beam axis. Besides providing an experimental basis for deeper investigations of the physical plasma-THz generation mechanism, such results allow one to assess and account for the effects of such a frequency-dependent beam profile at a subsequent focus. In the final section, we provide simulation results for the predicted development of the spatio-temporal field in a focal region and discuss its impact for z-scan-based nonlinear spectroscopy experiments.

Broadband terahertz (THz) emission
We begin by treating the 'conventional' case of a plasma-THz emitter with relatively long optical pump pulse (T ω = 150 fs). For the experiments here we employed a 1 kHz Ti:sapphire amplifier laser (Clark-MXR, CPA-2101) delivering pulses with an energy ∼800 µJ and centre wavelength of 775 nm. One attractive aspect of two-colour plasma-THz generation is that even with such an optical pulse duration one still can reach bandwidths exceeding 10 THz (10%maximum level). This significantly larger bandwidth compared to optical rectification is not only due to the absence of phase-matching and absorption effects present in χ (2) -crystals [20], but is also related to the high order of the effective nonlinearity in the plasma generation mechanism (whose first step is tunnel photoionization of the gas medium) [5,21]. The fundamental (ω) optical pump pulse (which is p-polarized in the laboratory axes) is focused to form a plasma in ambient air via a β-BBO crystal (150 µm thickness, cut angle 32 • , unless otherwise stated) where the second-harmonic optical pulse (SH, 2ω) is generated [5]. The azimuthal rotation of the BBO crystal is chosen to optimize the THz pulse energy, which corresponds to an angle of ∼45 • between the p-polarized fundamental and extraordinary axis of the BBO crystal. Following the plasma, a high-resistivity Si wafer (thickness 450 µm) at a distance d Si = 45 mm is used to discard the optical pump beams (which otherwise cause damage to the subsequent optics), where d Si is measured relative to the nominal focal plane of the optical pump focusing lens. The THz beam is subsequently collimated by an off-axis paraboloidal mirror (OAPM) with an effective focal length of f eff = 76.2 mm. The emitted THz pulse energies were typically in the range ∼20 nJ, as measured by a pyroelectric detector (Molectron P4-42) in a subsequent beam focus, accounting for the losses at two Si wafers in the intermediate beam path.
We begin by presenting the THz beam fluence profiles for different values of the focal length f for the optical beam before the plasma and input pulse energy J ω . For these measurements the OAPM directly after the plasma was removed and a Golay cell (QMC Instruments, OAD-7) on a dual-axis translation stage was used to acquire the beam profile via raster-scanning. Moreover, in order to gain information on the polarization state of the emission, a free-standing wire-grid polarizer was placed in the THz beam, and additional measurements were performed for both p-and s-polarized signal components (S p and S s , respectively) in addition to the signal S without polarizer. The extinction and peak transmission of the polarizer degrade gradually above ∼5 THz (with the extinction ratio rising to ∼40% at 10 THz). The measurement of all three signal scans allowed us to determine and correct the signals for this non-ideal power transmission for radiation polarised parallel (T = 0.95) and perpendicularly (T ⊥ = 0.05) to the wire axis.
In figures 1(a)-(d) we present the results of such measurements for f = 50, 100, 150 and 200 mm (all with J ω = 710 µJ). Here we plot the total fluence S versus the horizontal and vertical inclination angles to the beam axis (α and β, respectively). Also included in each figure are arrow fields indicating the relative magnitude of the p-and s-polarized components. The distance d BBO between the BBO crystal and the focal plane (which affects the relative ω-2ωphase entering the plasma due to the dispersion of the intermediate air path [5,22]) was held constant at d BBO = 75 mm for the case of f = 100, 150, 200 mm, while for f = 50 mm a value of d BBO = 40 mm was used. The distance d between the plasma focus and Golay cell detector was adjusted for each focal length in order to provide a suitable spatial scan range in each case, ranging from d = 80 mm ( f = 50 mm) to d = 190 mm ( f = 200 mm), corresponding to the far-field of the plasma emitter in all cases. Based on the 6 mm acceptance aperture of the Golay cell, this corresponded to a nominal angular resolution δα ranging from 4 • to 1.8 • . In all cases a clear minimum is observed on the beam axis, similar to that reported previously [19,23]. The profiles exhibit maxima at half-angles of α m = 10.7 • ( f = 50 mm), 6.4 • ( f = 100 mm), 4.7 • ( f = 150 mm) and 4.5 • ( f = 200 mm), respectively. Note that the present case with a two-colour excitation should not be confused with the case of the conical THz emission from a single-colour plasma filament [24] where the on-axis intensity must vanish due to lateral symmetry (which is broken in the two-colour case) and a radially polarized emission results.
The polarization vector fields shown in figures 1(a)-(d) show that the polarization orientation for each focal length is relatively constant over the whole beam profile, ruling  out any additional contribution from a radial polarization, which would not have been readily distinguished in previous analyses of the polarization state for the focused THz beam [22,25]. We note that the measurements here of the power along the two polarization directions does not yield the ellipticity of the polarization state in general, although the data for f = 200 mm (with a very small s-component) must correspond closely to a linear polarization. A trend with f is observed in the orientation of the polarization vector, which is approximately parallel to the 2ω-field polarization for f = 50 mm and becomes p-polarized for f = 200 mm. It is well-established that a variation in the ω-2ω-phase entering the plasma can give rise to a corresponding rotation in the THz polarization [22,25], e.g. by varying d BBO (where a 180 •rotation is induced by a change in d BBO of ∼26 mm). However, for the cases of f = 100-200 mm here, d BBO was kept constant (as was the orientation of the BBO crystal). One contribution to an ω-2ω-phase shift (and hence THz polarization rotation) may come from the fact that for the different focal lengths, the z-position for the onset of plasma formation shifts relative to the focal plane (changing the effective BBO-plasma distance). However, in going from f = 100 to 200 mm this amounts to a shift of only ∼1 mm, which is not sufficient to account for the ∼45 •rotation seen in the data. We note that a significant THz polarization rotation (32 • ) was also observed in changing only the optical pump energy (by a factor 2 [22]), which was attributed to ω-2ω phase walk-off effects inside the plasma due to plasma dispersion. Optical microscope measurements of the plasma profiles showed that the plasma length L p varied from L p = 1.5 mm for f = 50 mm to L p = 14 mm for f = 200 mm in the experiments here. Hence one can expect that the magnitude of the phase walk-off in the plasma depends strongly on focal length and lead to an analogous THz polarization rotation as in [22]. We note here that the Kerr cross-phase modulation between the ωand 2ω-fields has also been asserted as playing a role in addition to the plasma dispersion [8].
Within the previously reported 1D models to account for two-colour plasma-THz emission profiles [18,19], the beam cone angle α m is related to the optical-THz phase walk-off rate dϕ/dz in the plasma (leading to far-field constructive interference for the THz waves further from the axis with increasing dϕ/dz). Hence the dependence of α m on f is physically reasonable, as the peak achievable optical intensity (before plasma defocusing sets in [5]) is higher for shorter values of f . The plasma length L p , on the other hand, should affect the angular width of the cone, becoming sharper with increasing length. Such behaviour is indeed borne out in the data in figures 1(a)-(d). In order to investigate this further, we performed a set of 2D THz beam fluence measurements versus the optical pump energy J ω , as shown in figures 1(e)-(g) for f = 100 mm. Normalized 1D vertical slices of these profiles are shown in figure 1(h). Here it is evident that despite a strong reduction in the emitted THz pulse energy (a factor ∼5) in going from J ω = 715 to 310 µJ, the normalized profile remains relatively constant in shape. This is somewhat unexpected, as previous studies [5] have shown that such a variation of J ω should give rise to significant changes in the plasma length L p and hence the THz emission profile.
Up until now, we have only addressed the conical beam profile on the basis of the generation mechanism in the plasma and previously published models. However, as mentioned above, in these experiments a Si wafer is used to discard the optical beams after the plasma, raising the question as to what effect the resulting photoexcited charge-carrier profile should have on the THz beam. To address this, we used an approximate model to calculate the expected transmission profile of the Si wafer using the optical pump beam parameters used in the experiment (see appendix A for details), accounting for both Fresnel reflection and bulk absorption induced by the ωand 2ω-optical pulses. Such transmission profiles are shown in figures 1(i) and (j) for pump pulse energies of J ω = 310 and 715 µJ, respectively, and THz frequencies from ν = 1-10 THz. The on-axis optical fluence values F nω (r = 0) used for the simulation in each case are given in the caption. These model predictions indicate that for the majority of the spectral content of the beam, only a negligible or very low transmission should occur with a beam half-angle range of ∼2 • about the axis. The magnitude of these predicted photoinduced losses are supported by additional optical-pump THz-probe measurements (data not shown). Note that the finite on-axis fluence in figure 1(h) is likely influenced by the angular resolution (δα = 2.5 • ) used in the measurements. A comparison of figures 1(i) and (j) shows that the photoinduced loss profile is only weakly affected by the change in optical pump energy (as per the normalized experimental profiles in figure 1(h)), raising the question of whether this effect is actually predominantly responsible for the observed THz beam profiles. This issue can be resolved in the following, where we present frequency-resolved measurements of the THz beam profile. In the present case, this was achieved by guiding the collimated THz beam through a FTIR spectrometer arrangement (employing a 450 µm thick 100 mm diameter Si wafer as beamsplitter, as depicted in figure 2. By placing the Golay cell detector (equipped with an iris with diameter ∼1 mm) and scanning stage at the output of the interferometer, we could obtain a field autocorrelation for each vertical position in the beam. In figure 2(a) we plot a typical intensity spectrum of the full THz beam obtained from a Fourier transformation of the field autocorrelation measured at the output of the FTIR spectrometer (after refocusing the full beam onto the Golay detector), demonstrating the achievable bandwidth with T ω = 150 fs. The spatially resolved intensity spectra are shown in figure 2(b). (Note that in order to obtain a reasonable acquisition time for the full set of autocorrelation traces with acceptable signal-to-noise, the delay range of each scan was limited to 1.33 ps resulting in a spectral resolution of ν = 0.75 THz.) As can be seen, the axial intensity minimum becomes more pronounced with increasing THz frequency, while the spatial width of the corresponding off-axis peaks becomes sharper. Such a behaviour is indeed expected from the simple 1D emission model [18,19], as the off-axis phase-matching condition versus emission angle becomes more stringent as the THz frequency is increased. In figure 2(c) we apply this model to fit the experimental data, using an effective plasma length of L p = 2 mm and a phase walk-off rate based on a refractive index mismatch of n = n THz − n eff opt = 0.015. The total spectral intensity for each frequency was normalized to that of the experimental data. A good semi-quantitative agreement can seen between model prediction and experiment. We note, however, that in this analysis we have not accounted for the transmission profile of the photoexcited Si wafer (cf figures 1(i) and (j)), which should lead to a strong depression of the on-axis intensity especially for lower frequencies. That this is not observed in figure 2(b) is possibly due to the relatively long propagation distance (∼500 mm) to the output plane of the FTIR spectrometer, whereby diffraction effects following the Si wafer could lead to a partial 'healing' of the axial depression, especially for low frequencies. Hence, in the present case the application of a 1D emission model should be treated with caution, and in some sense the agreement with experiment is fortuitous. Nevertheless, one key aspect of the data is that while the Si loss profile is predicted to narrow and weaken as the THz frequency is increased, the experimental data reflect the opposite trend. Hence the photoinduced losses of the Si wafer cannot account qualitatively for the observed conical beam profile versus frequency and instead one must invoke a mechanism consistent with the predictions of the 1D emission model. The model frequency-dependent profile in figure 2(c) will be employed for the beam focusing simulations in section 4.

Ultra-broadband THz-mid-infrared emission
We now proceed to the case of ultra-broadband THz-mid-IR emission, that can be achieved using ultra-short optical pump pulses (T ω 20 fs). In this case, a continuous spectrum extending well beyond 100 THz can be achieved [9,13]. Previous modelling of this ultra-broadband generation process [9] indicated that even in a plane-wave treatment one must consider a relatively complex interplay of temporal and polarization effects during propagation in the plasma. Hence one could also anticipate that the spatial characteristics of the emission may exhibit a non-trivial structure, which may have an important impact on its use for spectroscopic applications.
For the experiments here, we employed a commercial hollow-core fibre compressor (Femtolasers KALEIDOSCOPE, containing Ar at 3.5 atm, length 1 m) to spectrally broaden the 150 fs pulses from the CPA-2101 laser, followed by a series of negative dispersion mirrors to obtain sub-20 fs duration [9]. As per the previous section, these pulses (J ω = 420 µJ) were focused to form a plasma in ambient air with a f = 200 mm lens via a 150 µm thick β-BBO crystal for generating the 2ω-pulse (placed at a distance d BBO = 40 mm before the focal plane). A BK7 wedge pair before the setup allows fine tuning of the ω-pulse chirp, which was adjusted to optimize the emitted THz-mid-IR pulse energy. The optimum pulse energy (J = 360 nJ, directly after the plasma) and intensity spectra were presented in a previous report [9].
After initial measurements of the 2D emission fluence profile, and considering the possible role of photoinduced losses in the Si wafer following the plasma, we adopted a measurement setup where the THz-mid-IR emission was allowed to diverge from the plasma and was measured using the Golay cell on the 2D raster-scan stage at a distance of 850 mm after the plasma, allowing us to vary the Si-wafer position d Si (and hence incident optical fluence F nω ) over a broad range. In figure 3(a) we plot the measured 2D fluence profile (for d Si = 30 mm), which, as per the broadband emission in the previous section, exhibits a pronounced ring structure with a maximum at a half-angle of α m = 3.2 • (compared to the value 4.5 • in the previous section for this optical focal length). In order to inspect the effect of the Si wafer, we THz (d Si = 800 mm). Note that F nω ∝ 1/d 2 Si and ω p ∝ 1/d Si hold for all intermediate values. (h) Measured FTIR field correlation and (i) corresponding intensity spectra for full THz beam, and with an on-axis aperture with an angular aperture with d Si = 130 mm limited to the axial peak region. The shaded region about each curve represents the 95% confidence interval for the signal (determined from a statistical analysis of multiple-scan data; 52-116 scans taken with lock-in time constant of 50 ms). repeated this measurement for d Si = 300 mm, as shown in figure 3(b). As can be seen, this gave rise to a new feature in the 2D profile, i.e. a new fluence peak on the beam axis. In order to study this axial peak in more detail, we performed high-spatial resolution line-scans of this central portion of the beam for a range of values of d Si , as shown in figure 3(c). Here one observes a smooth trend in the appearance of the axial peak with increasing d Si , approaching saturation for d Si = 300 mm, demonstrating that its absence for small values of d Si is due to the photoexcitation of the Si wafer. In order to model this effect, as per the last section, in figures 3(d)-(g) we plot the predicted transmission of the photoexcited Si wafer for the same values of d Si in figure 3(c) for frequencies of ν = 5, 10, 20 and 50 THz, respectively. A comparison with the experimental data indicates that the dependence of the axial peak on d Si is best reproduced by the curves for ν = 20 THz, suggesting that this spectral range contributes significantly to the appearance of the axial peak.
In order to investigate this experimentally, we collimated the THz-mid-IR beam using an OAPM with f eff = 152.4 mm and guided it to the FTIR spectrometer. This value of f eff allowed us to place the Si wafer at a distance d Si = 130 mm, where the appearance of the axial peak has already set in. In figure 3(h) we present measured FTIR intensity spectra for the full THz-mid-IR beam as well as with an intermediate iris centred on the collimated beam to transmit only the central portion corresponding to a maximal emission angle of α = 1.5 • and 1.0 • , respectively. A comparison of these intensity spectra reveals that the high-frequency spectral components (i.e. 60 THz) drop dramatically within the half-angle of α = 1.5 • , and hence correspond to a conical emission profile. As the aperture diameter is further reduced (α = 1.0 • ), the lower-frequency side of the spectrum is suppressed, indicating that the axial peak contains spectral components from low frequencies to ∼70 THz with a spectral concentration in the range ∼50-70 THz. Note that the case of α = 1.0 • corresponds to an aperture diameter of 2.6 mm-i.e. much larger than the wavelength of the spectral components affected in going from α = 1.5 • to 1.0 • , such that beam diffraction at the aperture should not play a significant role.
These measurements indicate a relatively complex spectral dependence for the emission pattern, and compel further detailed simulations of the ultra-broadband plasma emission process in the future, which will require one to go beyond simplified 1D models and treat the full spatiotemporal propagation in the plasma [5,26]. While the axial peak represents only ∼15% of the total emitted pulse energy, its distinct spectral and spatial properties suggest it may be connected to a second emission mechanism in the plasma in addition to a plasma photocurrent [21].

Sum-frequency detection
In order to characterize the ultra-broadband emission in more detail, we employ here coherent detection using the method of optical-THz sum-frequency generation (SFG, or 'up-conversion') in χ (2) -crystals [27] using 150-fs optical detection pulses. This has several advantages compared to FTIR detection (with typical thermal detectors), such as a higher sensitivity and signal-to-noise (especially at higher frequencies), as well as providing both spectral and temporal information through the measurement of SFG spectrograms. While our work on the optimization and deconvolution of such spectrograms is ongoing, one can already use this method to directly obtain estimates of the intensity spectrum of the THz-mid-IR radiation.
In figure 4(a) we present a spectrogram S(ν, τ ) of the ultra-broadband plasma emission, measured by focusing in a 110 -cut ZnTe crystal (Ingcrys, with a nominal centre thickness Corresponding frequency marginal M(ν) = dτ S(ν, τ ) shown at left. The pronounced spectral fringes are due to phase-matching oscillations. Also included (far left) is the model SFG phase-matching response function [28,29] (including response curve due to finite spectral width of optical detection pulse), generated using literature data [30] and assuming a crystal thickness of 68 µm (as required to match the periodicity of the phase-matching oscillations). of 25 µm on a 0.5 mm thick 100 -cut ZnTe supporting substrate) using an OAPM ( f eff = 50.8 mm). The optical detection beam was focused collinearly via a 2.5 mm diameter hole in the OAPM. A notch filter is used to block the central portion of the input optical spectrum before coupling the detected signal to a grating spectrometer. Also shown at the left in figure 4(a) is the spectral marginal M(ν), i.e. obtained by integrating S(ν, τ ) over delay τ , which provides an estimate for the intensity spectrum of pulse, albeit with an additional spectral modulation and envelope due to phase-matching oscillations (the latter resulting in a response roll-off with increasing frequency ν which goes as 1/ k ∼ 1/ν). As can be seen, a measurable SFG signal extends out to ∼150 THz, despite this spectral roll-off.
In order to gain information on the spatial beam distribution versus frequency, we placed an iris centred on the collimated beam axis, and recorded SFG spectra versus aperture radius R for a fixed detection delay (indicated by the vertical line in figure 4(a)). Note that this method was chosen over raster scanning an aperture or a knife-edge measurement in order to avoid any possible detection side-effects due to creating a transmitted beam with an off-axis centre. The measured SFG spectra versus aperture radius are shown in figure 4(c), after performing a local averaging over the phase-matching oscillations in order to produce a smoothly varying data set (at the expense of a small reduction in spectral resolution). As the aperture is placed in the collimated beam approximately a distance f eff before the OAPM, the aperture and focal plane constitute a Fourier optics pair [31]. As the optical detection beam essentially samples the THz-mid-IR beam on the focal beam axis, this results in a simple relation between the truncated collimated beam field profile E 0 (r, ν) and the detected SFG signal S(ν) ∝ |E(r = 0, ν)| 2 (see appendix B). Note that the Si wafer for these measurements was placed at d Si = 45 mm, i.e. such that the axial peak is suppressed and a conical beam profile is expected for all frequencies ( figure 3(c)). Hence, in order to fit the signal S(ν, R) versus R for each frequency, we had to adopt a model including axial suppression for the collimated beam field profile E 0 (r, ν). This was taken as a Gaussian beam profile modulated with an additional Gaussian loss profile to introduce the axial suppression, i.e.
which should provide a reasonable phenomenological description of the conical profile due to both the plasma emission mechanism and loss profile in the photoexcited Si wafer, and has the attractive mathematical property that the on-axis focal intensity versus R can be derived analytically (see appendix B) to allow rapid fitting of the experimental data. A plot of S(ν, R) for a representative frequency (ν = 40 THz) is shown in figure 4(b), along with the fitted curve, which is in good agreement. Also included is a model fit assuming an unmodulated Gaussian beam (i.e. without axial suppression), where it is apparent that the rate of decay as R → 0 is not sufficiently rapid to fit the data well. The model fits to the experimental data S(ν, R) for all ν are shown in figure 4(d), while the underlying fitted collimated beam profiles E 0 (r, ν) are shown in figure 4(e) (normalized for each frequency to have unity spectral intensity for the full beam). A smooth trend is observed in the collimated beam profile E 0 (r, ν) versus frequency, as the position of maximum fluence ranges from r = 5.3 mm for ν = 10 THz (corresponding to a half-angle α = 4.0 • before the OAPM, f eff = 76.2 mm) to r = 2.2 mm for ν = 150 THz (α = 2.2 • ). These results constitute a useful experimental basis for evaluating the predictions of future simulations of the plasma emission process, although referring to the results in figure 3, for the value of d Si = 45 mm used here the axial peak has been suppressed (larger values of d Si were not accessible with the SFG setup employed here). Nevertheless, the measured spectral dependence of E 0 (r, ν) clearly goes beyond the predictions of 1D emission models with a fixed optical-THz phase walk-off rate in the plasma. Moreover, the data do describe the actual experimental situation and hence serve as a valid input for the beam focusing simulations in the next section.

Simulation of spatio-temporal dynamics in focal region
In order to evaluate the consequences of the spatially modulated emission profiles presented in the previous sections for nonlinear spectroscopy experiments (in particular based on z-scan techniques), in this section we employ propagation simulations based on this data to predict the spatial and temporal behaviour of the fields traversing a focal region. These simulations are based on a fast Hankel transformation [32,33] to apply Helmholtz paraxial propagation for each frequency component (see appendix B) using the spectrally resolved intensity profiles established above as input, and assuming a linearly polarized (scalar) field with cylindrical symmetry. The model results are summarized in figure 5, for both broadband THz and ultrabroadband THz-mid-IR emission.
Beginning with the data for the broadband THz emission (T ω = 150 fs), as can be seen in figure 5(a), at the focal plane z = 0 the beam profile is approximately Gaussian in shape for all frequencies ν, with a beam radius going approximately as 1/ν. Here the modulation of the input beam profile (shown in the left panels for each frequency) manifests only as weak satellite features. In moving away from the focal plane, however, a more complex structure develops. This is to be expected, as over a length scale comparable to the effective Rayleigh range [34] one must have a transition from a scaled version of the input beam profile to its spatial Fourier transform (similar to the transition from Fresnel to Fraunhofer diffraction [31]). Due to the particularly strong modulation of the input profile for ν = 10 THz as well as the weak oscillations on the wings of the input profile, the simulations even predict the transient appearance of sharp local intensity maxima versus z on the beam axis.
In order to inspect how these features translate into the temporal properties of the pulse and achievable field amplitudes, in figure 5(c) we plot the corresponding simulated on-axis time-domain electric fields for a range of values of z. For these results, we propagated the beam profiles for a complete set of frequency components (N = 256) and employed the measured intensity spectrum of the full beam (figure 2(a)) to apply the correct amplitude for each component before Fourier transformation back to the time domain. The pulse energies used were J = 20 nJ (broadband) and J = 100 nJ (ultra-broadband), which are typical values obtained for the experimental conditions in sections 2 and 3 [9], respectively. To more closely model the experimental situation, we also applied a pulse chirp to the input field, i.e. by adding a groupdispersion delay (GDD) of ϕ = 30 fs 2 , which corresponds to the dispersion from a 100 µm thick Si wafer (although this only leads to significant effects for the ultra-broadband THzmid-IR pulse discussed below). At z = 0 the simulations predict a single-cycle sine-pulse (i.e. with a carrier-envelope phase of ϕ(t = 0) = π/2), which is as expected as the near-transformlimited input pulse corresponds to a cosine-pulse and focusing of the broadband pulse leads to a differentiation in the time domain for the on-axis field. For |z| 2 mm, the electric field waveform begins to broaden and develop a complex structure. This is due to the relatively complicated superposition of the various frequency components as they converge onto the beam axis (each with distinct behaviour). This effect is already significant at z-positions where one still has a peak field of ∼10% the maximum at z = 0 (which amounts to ∼1 MV cm −1 for the parameters used here). One also observes that the electric field profile undergoes both a polarity and temporal reversal in traversing the focus, which is due to the effect of the Gouy phase as discussed in other theoretical reports on the focusing of single-cycle pulses [35].
Hence, while the simulations predict a near-ideal beam and temporal profile at the focal plane for the broadband THz pulse, for experiments based on a z-scan technique one must bear in mind that significant spatio-temporal distortions may already be present within a relatively close z-range to the focus. These effects would be further exacerbated for experiments where the sample induces nonlinear spatial effects (e.g. Kerr focusing) where the spatial structure can be expected to affect the detected signals even more significantly.
We turn now to the corresponding focal simulation results for the ultra-broadband THzmid-IR pulses, as shown in figures 5(b) and (d). Again, an inspection of the spatial intensity profiles versus z reveal a near-Gaussian profile for z = 0 and the development of a non-trivial structure as one moves away from the focus. We note that due to the less drastic modulation and absence of satellite oscillations of the input profile for ν = 10 THz, the development of this structure for this frequency is less complex than that for the same frequency component in figure 5(a) (moreover, here the radial range shown is lower, such that a significant amount of the outer beam structure for ν = 10, 20 THz is not visible). The general qualitative differences between the spatial profiles in the two cases arise principally due to the greatly different ratio between the THz(-mid-IR) wavelengths involved and the focal length (which is f = 100 mm for both cases). In this sense, the focusing of the ultra-broadband radiation is actually more paraxial and well-behaved in nature than for the broadband case, despite covering a much broader absolute frequency range. (To avoid any confusion, we note here that the local axial peaks away from z = 0 in the focal profiles are not related to that observed experimentally for the collimated beam in figure 3(b).) The predicted on-axis temporal electric field profiles (figure 5(d)) also exhibit a distinct behaviour. Due to the GDD (ϕ = 30 fs 2 ) added in the simulation, a significant pulse broadening occurs with this ultra-broad bandwidth giving rise to a ∼two-cycle pulse with a discernible frequency modulation in the time-domain field. This pulse chirp also lifts the time-polarity reversal transformation for the field in traversing the profile, although a clear carrier-envelope phase change with z is still evident. Most notably, the relative degree of distortion on the temporal pulse envelope away from the focal plane is seen to be much milder for the ultrabroadband pulse, even though the z-range chosen here also corresponds to a distance where the on-axis peak field has dropped to ∼5% of the peak field at the focus, which is ∼10 MV cm −1 for the parameters used here (note that a peak field of ∼15 MV cm −1 results in the absence of any GDD in the simulation). Hence, in particular for experiments mostly dictated by the pulse intensity envelope (and not the sub-cycle field), one can expect only moderate temporal effects in performing z-scan measurements with such an ultra-broadband THz-mid-IR pulse, despite the strong modulation of the collimated beam profile and extreme broadband nature of the pulse.

Conclusions
The spatio-spectral characterization of both the broadband and ultra-broadband THz(-mid-IR) emission from two-colour optically excited plasmas has provided an extended picture with which to evaluate their use for nonlinear spectroscopy, where the resulting spatio-temporal features of the focused beam have a much larger impact on experiments than for linear spectroscopic measurements. Such information can also be used to validate the predictions of different theoretical models for the plasma emission, which in turn may be used to aid the rational optimization of the plasma emission for different experimental situations. These results extend the findings in previous reports [18,19,23] concerning the conical emission from gas plasmas with two-colour optical excitation, which we have also shown possesses a constant polarization state across the beam, in contrast to the emission from gas plasmas with single-colour excitation [36], where the conical emission and a radial polarization is imposed by the lateral symmetry and a different emission mechanism is operative. While the role of photoexcited losses in the semiconductor (Si) wafer, typically used directly after the plasma, has not received much attention in the past, we have shown that it should be taken into account if one wishes to use beam profile measurements to better understand the physical mechanism behind the plasma emission. From the results here, both the inherent conical emission mechanism and Si absorption should contribute to the spatial profile, especially at lower frequencies. For the ultra-broadband THz-mid-IR emission, the Si loss profile actually masks an on-axis emission feature whose distinct features suggest a second emission mechanism in the plasma. In terms of a medium which is sufficiently transparent across the whole THz-mid-IR range and opaque for the optical pump beams, Si and Ge are essentially the only practical candidates due to the effect of strong (single-)phonon resonances in other semiconductors, which not only lead to a spectral gap but more importantly give rise to a magnitude of dispersion which rapidly degrades the pulse duration (and hence, also peak field). Another option for discarding the optical beams lies in the use of a small circular obstruction on the optical axis after the plasma [6], i.e. blocking both THz and optical beams but without any dispersion/absorption for the THz beam periphery. In this sense, the conical emission characteristics are actually favourable for achieving higher throughput. Using the results here, we were able to perform simulations of the focal spatio-temporal fields, which require spectrally resolved profiles as input in order to bear out the detailed behaviour. Despite the spatial modulation of the emitted beam, the predicted field distribution at the focal plane has a near-ideal Gaussian spatial profile with little temporal distortion, such that peak fields of ∼10 MV cm −1 are predicted for the ultra-broadband THz-mid-IR emission with realistic experimental parameters. Nevertheless, for nonlinear z-scan measurements, care should be taken to consider the role of spatio-temporal distortions of the field away from the focal plane. In future investigations, we will aim at an experimental verification of the simulation results here, in order to aid the application of theoretical calculations (which are often carried out assuming ideal Gaussian temporal and spatial profiles) for the nonlinear response of various materials. by assuming a Gaussian temporal intensity profile for the optical pulses I ω,2ω ∝ e −2t 2 /T 2 0 and calculate the cumulative fluence as a function of time. The average transmission is then obtained by numerically integrating the time-dependent transmission over the time interval t ∈ [−T 0 , T 0 ]. (Note that a detailed co-propagation of the THz and optical pulses in the Si medium is not currently accessible, e.g. using finite-difference methods [38], as it would require a detailed knowledge of the temporal THz pulse profile.) Using literature values for the absorption coefficients α ω,2ω [39] the excitation density at radius r , depth z and time t can be written as N (r, z; t) = N ω (r, z; t) + N 2ω (r, z; t) where N nω (r, z; t) = (1 − R f (nω)) α nω hnω e −α nω z F nω (r ) The on-axis field is simply given by For the modulated field profile used in section 3 (equation (1)), E 0 (r ) = E 00 e −r 2 /w 2 (1 − β e −r 2 /a 2 ), with truncation by an aperture with radius R we then have the result Instead here we exploit the cylindrical symmetry and make use of the Hankel transformation H [32,33], which allows one to map the spatial and spectral domains to their corresponding radii r = x 2 + y 2 and ρ = √ u 2 + v 2 . Hence definingẼ(ρ, z) = H[ρ]{E(r, z)} we haveẼ(ρ, z) =Ẽ(ρ, 0) e i 2k 0 ρ 2 z , while the focal field profile is given in terms of the input collimated field by E(r, 0) = i 2π f k 0 H −1 [r ]{E 0 ( f k 0 ρ )}, where ρ is related to the input radius variable r by ρ = (k 0 / f )r . In order to implement the discrete forward and inverse Hankel transforms, we use the improved-accuracy quasi-fast Hankel transformation [33], which maps the radial variables onto domains whose sample points form a geometric series and reduces the discrete Hankel transformation to the application of three 1D FFTs, which is orders of magnitude faster than the corresponding 2D FFT calculation. Besides exploiting the dimensionality reduction of the problem, this numerical method has a distinct advantage here for the propagation of the frequency components of an ultra-broadband (single-cycle) pulse, i.e. the geometric sampling points allow one to realize simultaneously high sampling rates near the beam axis and large sampling intervals further from the axis, such that a single pair of r -and ρ-domains can be used for all frequencies, despite a wavelength (and hence focusing) ratio of ∼100 across the spectrum.