Tailoring THz radiation by controlling tunnel photoionization events in gases

Applications ranging from nonlinear terahertz spectroscopy to remote sensing require broadband and intense THz radiation which can be generated by focusing two-color laser pulses into a gas. In this setup, THz radiation originates from the buildup of the electron density in sharp steps of attosecond duration due to tunnel ionization, and subsequent acceleration of free electrons in the laser field. We show that the spectral shape of the THz pulses generated by this mechanism is determined by superposition of contributions from individual ionization events. This provides a straightforward analogy with linear diffraction theory, where the ionization events play the role of slits in a grating. This analogy offers simple explanations for recent experimental observations and opens new avenues for THz pulse shaping based on temporal control of the ionization events. We illustrate this novel technique by tailoring the spectral width and position of the resulting radiation using multi-color pump pulses.


INTRODUCTION
Recently, THz technology has attracted significant attention, because it can provide unique analytical and imaging tools in nonlinear and time-domain spectroscopy, remote sensing, biology and medicine, security screening, information and communication technology (see e.g. [1][2][3][4]). Among the various THz sources, employing two-color ionizing femtosecond pulses in gases [5][6][7][8][9][10][11][12][13][14] stands out through absent damage threshold, low phonon absorption and zero interface reflection. The two-color scheme provides striking performance in terms of bandwidth that can exceed 100 THz and thus even covers farand mid-infrared, as well as high peak fields up to the MV/cm range.
In the standard realization of this method a femtosecond laser pulse and its second harmonic are collinearly focused into a gas, such as ambient air, yielding powerful THz radiation with well-defined electric field values. Initially, this process has been explained by four-wave mixing rectification via third-order nonlinearity [5,12]. However, a threshold for THz generation [6,10] as well as plasma current measurements [13] indicated that plasma formation plays an important role in the generation process. Besides, the Kerr nonlinearity is far too small [7,15] to explain the observed THz field strength in such a setup. In order to solve this discrepancy, Kim et al. [11,13] have put forward the so-called photocurrent mechanism. In this mechanism tunneling ionization and subsequent electron motion produce a quasi-DC photoinduced current, which in turn emits THz radiation. The quasi-DC plasma current can be efficiently produced in a gas irradiated by two-color laser pulses [5][6][7][8][9][10][11][12][13][14], in multicolor pump schemes [16] and in chirped [17] or few-cycle pulses [18,19].
Results published so far [10,13,15,16,[20][21][22][23][24] provide several important features of the THz spectrum. However, a general framework for analysis, control and design of THz radiation is still missing. Similar to pulse shaping techniques in the optical range, the broadband control of THz waveforms and its spectral properties is highly desirable for many applications, while the existing techniques are applicable only for a narrow bandwidth of few THz [23,25,26].
In this paper we advance the understanding of the fundamental physics in THz generation in gases as well as examine its basic mechanism associated with plasma formation. We demonstrate that THz generation in gases is intrinsically connected to optically-induced stepwise increase of the plasma density due to tunneling ionization. New spectral components are emitted in a discrete set of attosecond-scale, ultra-broadband bursts associated with these ionization events; in the spectral representation THz emission results from an interference of these radiation bursts.
Using appropriate field shapes of the pump pulses, e.g. tuning the frequencies and phases of two-or multi-color pulses, the temporal positions of the ionization events can be controlled. In turn, this allows to design the emitted radiation within a broad spectral range of more than 100 THz. This new concept of prescribing the phase and amplitude of each elementary contribution responsible for THz emission has a remarkable analogy in linear diffraction theory: The produced spectrum can be understood like the far-field diffraction pattern of a grating, where the ionization events play the role of the slits. A comprehensive (3+1)-dimensional non-envelope model of propagation is used to verify by means of numerical simulations that the simple picture drawn above captures the essential physics of THz generation.

STRUCTURE OF THE THZ SPECTRUM
We start our analysis in the so-called local current (LC) limit, that is, we consider a small volume of gas irradiated by the strong pump field E(t) [see Fig. 1(a)]. Electrons produced due to ionization build up a current J(t) which, in turn, produces an electromagnetic wave E J (t) ∝ ∂ t J(t) [27]. We focus here on the pump intensities where tunnel ionization is the dominant ionization process. In this regime, the free electron density ρ(t) increases stepwise in short attosecond-scale ionization events corresponding to maximum amplitudes of the pump field at times t n [see Fig. 1(c)]. This stepwise increase of the free electron density ρ(t) was recently confirmed experimentally [28,29]. The ionization events are well separated from each other and thus ρ(t) can be written as a sum over contributions from all ionization events: ρ(t) = n δρ n H n (t). Each ionization event has a well defined amplitude δρ n and shape H n (t). Since all ionization events have almost the same shape [see green filled curve in Fig. 1 where H is a "smoothed" step function. An approximate analytical formula for H(t) is given in the Methods.
In Fig. 1(b-e) an exemplary two-color pump field where ω 0 = 2π × 375 ps −1 is the fundamental frequency (corresponding to the wavelength of 800 nm), r = √ 0.2 and θ = π/2 are relative phase and amplitude of the second harmonic (SH) component. For such a two-color pulse (and assuming r < 1) the positions of ionization events are given by ω 0 t n ≈ πn − 2(−1) n r sin θ [15]. E(t) is a Gaussian pulse envelope shape with full-width at half-maximum (FWHM) pulse duration τ in = 40 fs.
With the assumption that electrons are born with zero initial velocity, the plasma current J(t) [blue curves in Fig. 1 where γ is the electron-ion collision rate (γ ∼ = 5 ps −1 at atmospheric pressure) and m, q are electron mass and charge [11,15]. Then, in Fourier domain the generated electromagnetic wave in the LC limit E J (ω) ∝ ωJ(ω) reads (see Methods) Here, dτ is the free electron velocity, and g is a constant.
Eq. (2) allows us to identify the impact of each ionization event. The interpretation of Eq. (2) is straightforward when we consider the spectral dependence of the amplitudes A n (ω) and B(ω), which can be calculated analytically in a reasonable approximation (see Methods): where σ ∼ = 10 fs −1 is the spectral width of a typical ionization peak [see thin black curves in Fig. 1(c)]. In contrast to B(ω), the amplitudes A n (ω) depend nontrivially on the pump pulse. Examples of A n (ω) and B(ω) are shown in Fig. 1(d,e). B(ω) features two clearly separated frequency scales determined by γ and σ. On the short frequency scale ω γ (up to 5 meV or ∼ 0 − 2 THz for atmospheric pressure), B(ω) ∝ ω, while for larger frequencies B(ω) depends only slowly on ω (on the scale of ∼ 10 eV).
If we consider the range from a few to approximately 100 THz, Eq. (2) can be significantly simplified. In this range B(ω) can be considered as a constant, and additionally the A n (ω) are negligible as seen from Fig. 1(d).
In other words, the resulting THz spectrum is a simple linear superposition of the contributions E J n = C n e iωtn and its structure is mainly governed by interference due to the spectral phases e iωtn . This interference is illustrated in Fig. 1(e) [cf. also Fig. 2(b)]. Note that in temporal domain Eq. (4) corresponds to a set of sharp peaks localized near t n , each of them having the amplitude C n and the shape ∂ t H(t − t n )e −γ(t−tn) . Thus, Eq. (4) describes sharp attosecond-long bursts of radiation localized near the ionization events.
Further on, the same simple superposition principle is also valid for higher frequencies in the more general Eq. (2). Indeed, Fig. 1(d) shows that constructive spectral interference in Eq.
(2) appears also for frequencies ω ∼ lω 0 , where l is integer. This finding is consistent with recent experimental results [29], where optical harmonics generated due to tunneling-ionization-induced modulations of the electron density were observed. Therefore, our superposition principle may also offer a complementary understanding of Brunel radiation [29,30].
In fact, Eq. (4) shows that THz generation can be understood completely analogous to far-field interference patterns produced by a diffraction grating [31], where the ionization events play the role of the slits. Each of such "slits" produces a broad "secondary wave", and their interference results in narrow lines, according to the Huygens-Fresnel principle [31]. THz radiation then corresponds to the zeroth diffraction order. In contrast to diffractive gratings, interference here takes place in the frequency domain and not in position space.
In a somewhat oversimplified picture, one can explain the spectrum shown in Fig. 1(d) analogous to interference patterns produced by diffraction gratings of equidistant two-color pump field E(t) (red curve) and current J(t) (blue curve) induced by tunneling-ionization, (c) step-wise modulation of the free electron density (filled green curve) and ionization rate W (t) (thin dotted black curve), besides, few cycles of the electric field E(t) (red curve) and the current J(t) (blue curve) are shown. (d) Spectrum |E J (ν)| 2 vs. frequency ν = ω/2π on a large frequency scale (red curve). Each ionization event at time tn induces a broadband radiation described by the form-factors B(ν) (black curve) and An(ν) (blue curve). The sum of all partial contributions in Eq. (2) yields the resulting spectrum E J (ν). (e) Partial contributions Re E J n (ν) from the nth ionization event for −16 < n < 16 (grey curves) and the total spectrum Re E J (ν) (filled red curve) in the THz to mid-infrared spectral region (Im E J (ν) is negligible here). The black curve shows B(ν). The decay of E J (ν) for ν 20 THz is solely due to destructive interference of different E J n (ν). In contrast, in the shaded region (ν < 5 THz) E J (ν) decays because of decrease of B(ν). (f) Spectrum |E J (ν)| 2 vs. input pulse duration τin. As the number of ionization events decreases with τin, the spectral width increases, as also confirmed by 3D simulations for τin = 20, 40, and 80 fs (black curves).
slits, E J (ω) ∼ sin(N ωδt/2)/ sin(ωδt), where N is the number of ionization events and δt ∼ π/ω 0 is the time interval between two subsequent ionization events (see Auxiliary Material). The width of the spectral line is inversely proportional to N , that is, to the input pulse duration [see Fig. 1(f)]. Such inverse dependence of the THz spectral width on the pump duration was already observed experimentally [7]. It is also compatible [black lines in Fig. 1(f)] with comprehensive 3D simulations, where propagation and transverse effects are accounted for (see Methods).

CONTROL OF THE THZ SPECTRUM
Understanding the spectrum of the radiation emitted by the plasma current J as an interference pattern provides a deeper insight into the mechanism responsible for the THz generation. Beyond that, this understanding provides us with the possibility to shape the emitted radiation by suitable engineering of the pump field. Here, we will restrict ourself to the THz and far-infrared do-main described by Eq. (4). The latter indicates that, in order to control the THz spectral shape, we have to target the times of the ionization events t n and the values of C n = qδρ n v f (t n ).  (1) and (2)], in addition to the known dependency of the THz yield [13]. This observation can be readily explained within our linear interference framework. For θ ∼ = 0 (case 1) the coefficients C n are not sign definite [red and blue circles in Fig. 2(a)]. Thus, for ω ∼ = 0 different contributions in Eq. (4) cancel each other, that is, destructive interference dominates. Nevertheless, with increasing frequency ω the destructive interference is (partially) compensated by the spectral phases e iωtn . This behavior is shown in Fig. 2(b), where contributions from the ionization events E J n are visualized in the complex plane for different frequencies ω. Consequently, the spectral maximum is shifted from zero [see Fig. 2(c)]. In contrast, for  (1) and (3) destructive interference occurs at ν ∼ = 0 and is partially compensated at larger frequencies. (c),(d) Corresponding THz spectral shapes and temporal profiles obtained from Eq. (4) ( n E J n , filled green curve) and from 3D simulations at the linear focus z = 5 mm (dashed yellow curves). Note that the amplitudes for case (1) are about one order of magnitude smaller in comparison to the other cases. θ = π/2 (case 2) all coefficients C n are positive. This obviously ensures that at ω ∼ = 0 constructive interference dominates, and the phase factor e iωtn can lead only to destructive interference when ω is increased. Therefore, the spectral maximum is located at zero frequency in this case. Moreover, the coefficients C n are much larger in amplitude, which explains the increased THz yield for θ = π/2. We verified this behavior by comprehensive (3+1)dimensional (3D) numerical modeling (see Methods), which was recently successfully compared to experiments [15]. We consider a Gaussian input beam with width w 0 = 100 µm and duration (FWHM) τ in = 40 fs for the 800 nm fundamental pump pulse (pulse energy ∼ 300 µJ). The energy of the second harmonic at 400 nm is chosen as 7.5 % of the fundamental, with duration and width smaller than the values for the fundamental by a factor √ 2 (similar to typical experimental conditions). These pulsed input beams are focused with f = 5 mm into argon at atmospheric pressure. In contrast to the simple LC model, the phase angle θ between the two components shifts during propagation. For our particular configuration, an initial value θ in ∼ = π/4 guarantees maximum THz energy yield of ∼ 1.8 µJ; θ in ∼ = −π/4 leads to a minimum yield of 0.25 µJ; the generic value θ in = 0 yields ∼ 1 µJ. THz spectra obtained for the two extremal cases are shown in Fig. 2(c), cases (2) and (1), respectively. They agree remarkably well with the predictions obtained from Eq. (4).
We also observe a high sensitivity of the THz spectra towards detuning of the relative frequency ratio between the two pump components [see case (3) and (4) in Fig. 2]. Such pulses with so-called incommensurate SH (frequency 2ν 0 is shifted by δν) offer a straightforward control of the position of the maximum THz spectral intensity, without decreasing the total THz yield. Recently, this shift of the THz spectrum with δν was predicted [20] and demonstrated experimentally [21]. Again, we can use our linear interference framework to explain this effect. The responsible shifts of the ionization events can be computed approximatively (see Auxiliary material) as ω 0 t n ≈ nπ − 2r(−1) n sin (θ ± πnδν/ν 0 ). The SH frequency shift ±δν makes the effective phase angle n dependent, and thus shifts the ionization events in time. Using our diffraction-grating analogy, the effect corresponds to increasing or decreasing the slit-distance along the grating. This frequency shift can be also calculated analytically under some further approximations (see the Auxiliary Material).
Because the coefficients C n are not sign definite [see Fig. 2(a,b), the "secondary waves" are not in-phase], such non-equidistant grating can lead to a spatial shift of the zeroth diffraction order, i.e., a frequency shift of the maximum THz spectral intensity. As it is seen from Fig. 2(c) and Fig. 3(a) the THz frequency shift is not symmetric with respect to the sign of δν. This effect occurs for pulses with finite duration (see Auxiliary Material).
Again, our predictions from Eq. (4) are in good agreement with 3D numerical simulations. Incommensurate two-color pulses clearly produce the predicted THz spectral shift for negative SH detuning [see Fig. 3(c)]. Moreover, neither linear nor nonlinear propagation effects of the pump pulse due to focusing or laser-induced plasma seem to jeopardize this effect. The energy in the shifted THz pulse is ∼ 1 µJ, like in the corresponding commensurate case [see Fig. 3(b)]. On the other hand, Fig. 3(d) reveals that, as expected, a positive SH detuning does not alter the THz spectral shape significantly, but the energy yield reduces by a factor of two to ∼ 0.5 µJ.
Let us finally consider a more complex example of THz spectral engineering. Namely, we attempt to significantly increase the spectral width by using a field shape which optimizes the "amplitudes" C n of just a few ionization events, suppressing the others. This leads to THz-tomid-infrared supercontinuum, which is of great importance for many applications. For this purpose, we resort to a three-color pump setup involving an optical parametric amplifier (OPA). We assume a fundamental frequency ν 0 = 375 THz (λ 0 = 800 nm) and signal and idler frequencies ν s and ν i obeying ν s + ν i = ν 0 (ν s = 0.55ν 0 , ν i = 0.45ν 0 ). From Fig. 4(b), one can see that in such a configuration only a few ionization events have large amplitudes |C n |. This immediately leads to a much broader THz spectrum, as it is obvious from the above discussed diffraction grating analogy. In Fig. 4(c,d) the results of 3D simulations for this case are shown, assuming that each of the three pump components contains one third of the total pulse energy of 300 µJ while all other parameters are kept unchanged.
Quite interestingly, the three-color configuration we use provides much larger values for the |C n | than the previous two-color scheme with comparable pump energy [cf. Fig. 4(a) and (b)]. In the LC limit, the resulting THz yield is about 40 times larger. In our simulations, this increased THz yield is indeed visible in the beginning of the propagation z < 2.5 mm. However, upon further propagation saturation effects, mainly coming from plasma defocusing of the pump and depletion of neutral atoms, limit the total THz energy to a few µJ.

CONCLUSION AND OUTLOOK
In summary, we have shown that THz emission in gases by multi-color laser fields is intrinsically connected with the attosecond jumps of the plasma density induced by tunneling ionization. The THz spectrum is described by a superposition of the spectrally ultrabroad contributions from separate ionization events.
This approach explains basic features of THz radiation such as the inverse dependence of the spectral width on the pulse duration and the shift of THz spectral maxima for non-integer frequency offsets of two-color pumps. Moreover it allows to manipulate the THz waveform and to shape the THz spectrum by controlling the tunneling photo-ionization. In the present article we demonstrated only basic features of such a control, using frequencies and phases of two-or three-colour pump fields. We note that other degrees of freedom of light such as polarization or spatial shape of the pump pulses can also be used for THz control. Noteworthy, our approach is not limited to THz, far-and mid-infrared range but can be potentially extended to frequencies much higher than the pump one, i.e., to optical harmonics.
We believe that this novel strategy to control fundamental properties of THz pulses could become a basic tool in optical technology: a ultra-broadband pulse shaper working in the range from THz to mid-infrared. Just as optical pulse shapers for femtosecond pulses are indispensable in ultrafast physics, the proposed THz control can find promising applications in many fields of research, including, among others, THz time-domain spectroscopy, chemical reaction dynamics, and optical signal processing.

Discrete ionization model
The electron density ρ(t) is described by the equation where ρ at denotes the neutral atomic density. We use a quasi-static tunneling ionization rate for hydrogen-like atoms [10,15]: where E a = m 2 q 5 /(4πǫ 0 ) 3 4 , ω a = mq 4 /(4πǫ 0 ) 2 3 , and r H = U Ar /U h . U h and U Ar are the ionization potentials of hydrogen and argon, respectively.
Supposing that the ionization process occurs only near maxima of the pump electric field at times t n , one finds E(t) ≈ E 0 + E 2 (t − t n ) 2 with E 2 = ∂ tt E(t n )/2 and E 0 = E(t n ). Then, assuming ρ at ≫ ρ, we can evalu- with Strictly speaking, σ is different for every ionization event.
The difference, however, is typically less than a few percent, so we assume here a "typical" value σ ∼ = 10 fs −1 .
The free electron current J is given by which is equivalent to Eq. (1) (see [11,15,22]). Using Eq. (6), the LC J can be written as a sum over all contributions from separate ionization events J(t) = n J n (t), and hence In the LC limit, assuming that the plasma source is contained in a small volume ∆V , the field produced by the plasma current J(ω) in frequency domain is E J = gωJ(ω) [27], where g = i ω∆V 4πǫ0c 2 r exp (iω r c ), c is the speed of light, ε 0 is the vacuum permittivity and r is the distance between the plasma source and the point where the field is measured.
Thus, in order to obtain Eq. (2), we have to multiply the Fourier transform of Eq. (7) by ωg. Further on, by performing the Fourier transform of H(t) exp (−γt) and neglecting a slowly varying phase term exp (−2iγω/σ 2 ) we obtain Eq. (3).

Comprehensive numerical model
We use the unidirectional pulse propagation equation [15,36] for linearly polarized pulses Here, E(k x , k y , z, ω) is the relevant electric field component in Fourier domain with respect to the transverse coordinates (x, y) and time, k = ωn(ω)/c is the wavenumber, c is the speed of light and n(ω) is the linear refractive index of argon [37]. The nonlinear polarization P NL (ω) = P Kerr (ω) + iJ(ω)/ω + iJ loss (ω)/ω accounts for third-order nonlinear polarization P Kerr (t), electron current J(t) and a loss term J loss (t) due to photon absorption during ionization. In our 3D numerical code, Eqs. (8), (1), and (5) are solved using a standard spectral operator splitting scheme. Particular attention is paid to resolve all relevant time scales, from tens of attosecond (ionization steps) to a few picosecond (THz radiation). Although, strictly speaking, Eq. (8) does not describe the propagation of waves below the plasma frequency properly, this plays only a minor role because the THz radiation is emitted from the leading plasma front and thus does not propagate in the plasma [38]. Spectra obtained using Eq. (8) were found to be in good agreement with experimental results [15].

Auxiliary Material
We consider for simplicity the case of a two-color pulse with duration T 0 having the square-shaped envelope E(t) = 0 for |t| > T 0 /2 and for |t| ≤ T 0 /2. Here δω is a detuning from the second harmonic, which is assumed to be small compared to the main frequency (δω ≪ ω 0 ). We also neglect the plasma decay γ and assume r ≡ A 2 /A 1 ≪ 1.
To determine the positions t n of the tunneling ionization events (maxima of the electric field) along the time axis we equate the time derivative of Eq. (9) to zero, set t n = t n0 + δt n and expand the resulting expression up to first orders in δω and δt n , yielding ω 0 t n ≈ nπ − 2r(−1) n (nπδω/ω 0 cos θ + sin θ) ≈ nπ − 2r(−1) n sin(θ + nπδω/ω 0 ), where it is assumed that πnδω/ω 0 ≪ 1. Using Eqs. (9), (10) and neglecting again higher-order terms, we find where a 1 = −3qA 1 r sin θ/2mω 0 , We further assume that the amplitude δρ n of each step in the modulation of the plasma density does not depend on the index n. For the square-shaped pulse envelope considered here this essentially means that saturation of the free electron density is neglected. Then, according to Eq. (4) of the main article where Here, the summation runs over all ionization events. Following Eq. (10), the exponent e iωtn can be expressed as e −2irω sin θ+ing+ω for even n and e 2irω sin θ−ing−ω for odd n, whereω = ω/ω 0 is the normalized frequency and g + = π(2rδω cos θ/ω 0 + 1), g − = π(2rδω cos θ/ω 0 − 1).
Considering an odd number of ionization events N , being much larger than unity, we calculate the expressions in Eq. (14) by performing the summation from −M to M , where M = (N − 1)/2. The term S 1 in Eq. (15) comprises simple geometric series and gives after separate summation for odd and even n: S 1 ≈ e 2irω sin θ sin (M g +ω ) sin (g +ω ) + e −2irω sin θ sin (M g −ω ) sin (g −ω ) , where M ≫ 1, so that M + 1 ≈ M . A similar analysis is possible for the term S 2 in Eq. (16). As before, one performs the sum for odd and even values of n separately. We are then left with an arithmetico-geometric series [39], i.e., for arbitrary q and integers m 1 , m 2 : m2 n=m1 nq n = m 1 (q − 1) + q + q m1+m2+1 [m 2 (q − 1) − 1] q m1 (q − 1) 2 .
Important features can be refound from the above expressions. In particular, for δω = 0 we have a 2 = 0 and hence E J ∝ S 1 . In addition, g + = −g − = π and thus, by introducing δt = π/ω 0 we obtain: The equation given at the end of Sec. II in the main article is obtained from Eq. (20) by assuming θ = 0 and N ≫ 1, that is, M = (N − 1)/2 ≈ N/2. For nonzero detuning δω we can at least qualitatively reproduce the behavior presented in Fig. 2 (case 3) and in Fig. 3(a) of the main article, i.e., the shift of the maxima of the THz spectra from zero frequency. The frequency shift obtained from Eq. (14), Eq. (17), Eq. (19) is, however, symmetric with respect to the sign of δω, in contrast to the full LC model. The asymmetry reported in the main text [see, e.g., Fig. 1(d)] disappears with increasing of the pulse duration and thus is induced by the presence of the Gaussian pulse envelope.