Universal properties of locally generated terahertz waveforms from polarization-controlled
two- and multicolor ionizing fields

H. Alirezaee,¹ S. Skupin,² V. Vaicaitis,³ A. Demircan,⁴,⁵ I. Babushkin,⁴,⁶,⁵ Luc Bergé,⁷ and U. Morgner⁴,⁵

¹Center for Research on Laser and Plasma, Shahid Chamran University of Ahvaz, Iran

²Université Claude Bernard Lyon 1, CNRS, Institut Lumière Matière, UMR5306, F-69100 Villeurbanne, France

³Laser Research Center, Vilnius University, Saulėtekio 10, Vilnius LT-10223, Lithuania

⁴Institute of Quantum Optics, Leibniz University Hannover, Welfengarten 1, 30167 Hannover, Germany

⁵Cluster of Excellence PhoenixD (Photonics, Optics, and Engineering - Innovation Across Disciplines),
Welfengarten 1, 30167 Hannover, Germany

⁶Max Born Institute, Max-Born Straße 2A, D-12489 Berlin, Germany

⁷Centre Lasers Intenses et Applications, Université de Bordeaux-CNRS-CEA, 33405 Talence Cedex, France

<figure><img src="image_1.png" /></figure>

(Received 7 June 2024; accepted 6 June 2025; published 14 July 2025)

The polarization states of terahertz (THz) radiation generated in a photoionized gas driven by intense two- or multi-frequency fields with locally controlled polarization are studied. We reveal a universal property of the resulting THz waveforms: the ellipticity of their polarization state increases linearly with frequency. This “linear chirp of ellipticity” makes plasma-based THz generation unique among other THz sources. However, it also imposes some constraints on the polarization properties of the generated THz radiation. We derive a general expression for the THz ellipticity and demonstrate how the polarization states of the generated THz waveforms can be manipulated and controlled by the polarization of the pump pulses.

DOI: 10.1103/dj4k-q4x2

I. INTRODUCTION

Terahertz (THz) radiation has various applications, such as characterization, spectroscopy, and imaging of materials, nanostructures and plasmas [1–8], remote sensing [9–11], and biomedical applications [1–3,12], to name a few. Generation of strong broadband THz fields by two-color laser-induced plasmas has attracted significant attention because of the high field strength and broad spectral bandwidth of the emitted radiation [13–17]. In this method, THz radiation is generated by photocurrents produced by ionization in an asymmetric driving field consisting of the fundamental pulse (FH) and its second harmonic (SH). Recently, considerable efforts have been devoted to optimizing and controlling various characteristics of THz radiation from two-color laser plasmas. In particular, the efficiency of THz generation was shown to be influenced by the relative phases between laser harmonics [15,18,19]. Other important factors affecting the efficiency of the produced THz radiation are pump wavelengths [20–23], polarization of the pulses [24–34], their durations [19,35,36] and, more generally, their waveshapes [19,25,37–40]. In addition to the conventional two-color pump configuration, multicolor pump fields have attracted considerable interest in recent years [19,33,38–47].

Many THz applications require polarization control
[5,6,48–50]. The most common techniques for controlling
THz polarization are based on wire-grid polarizers [48,51].
However, even if they might work in the whole THz re-
gion, they accept only narrowband THz pulses. So far, most
studies on THz radiation with two-color schemes have con-
sidered linear polarization of the pump. Apart from specific
dedicated works [24–34], little attention has been paid to
pump pulses with elliptical and circular polarizations. In
filaments or other extended propagation geometries of the
laser driver, local THz fields, produced at different points
in space with different polarization and phases, are superim-
posed, resulting in complex THz wave shapes [25,26,52,53].
In these propagation scenarios, the THz field generated down-
stream is also influenced by the previously generated THz
components [53].

In this paper, we consider the most important build-
ing block of ionization-based THz generation – the locally
generated THz wave forms in small spatial volumes. We
address a rich variety of possible THz fields, which can
appear in two- or multicolor schemes with arbitrary po-
larization of the pump pulses, and particularly focus on
the polarization state of such THz radiation. Starting from
general principles, we show how elliptically polarized THz
radiation arises, what states of polarization are possible,
and to what extent these can be controlled by the pump
configuration. We illustrate our findings by systematically
studying the dependence of the THz wave shape on the
pump configuration for a two-color scheme with controlled
polarization and phase difference between the two pump
components.

<footnote>Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International license. Further
distribution of this work must maintain attribution to the author(s)
and the published article's title, journal citation, and DOI.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) General setting. An intense complex pump field, characterized by multiple harmonics of some fundamental frequency $\omega_0$ together with their initial polarization states (visualized as ellipses) and amplitudes, is focused into a gas, producing plasma and thereby generating THz radiation with its own polarization state. (b) Schematic spectrum showing polarization angle $\Psi$, see Fig. 2, ellipticity $\epsilon$, and intensity of the resulting THz radiation, which exhibits an intrinsic “ellipticity chirp”.</figcaption></figure>

## II. SYSTEM AND MODEL

### A. Vectorial formulation of the local current model

The scheme of our setup is illustrated in Fig. 1(a). A laser pulse with two and more frequency components is focused into a gas, producing plasma and thereby low-frequency THz radiation via the so-called Brunel mechanism [7,23,41,54] (see also below). Although in general incommensurate combinations of frequencies can be used, here we assume that all frequencies are harmonics of the fundamental carrier frequency $\omega_0$. Whereas evolution of the resulting polarization can be very complex due to propagation effects [25,26,31,52,53], here we focus on the response of an infinitely small emitting volume of plasma at a remote detector.

The electric field of the pump pulse in the multicolor case, consisting of the fundamental at frequency $\omega_0$ and its harmonics $m\omega_0$, $m = 2, 3, ...$ is given by

$$ \mathbf{E}(t) = \sum_{m \ge 1} f_{m\omega_0}(t) e_{m\omega_0} [\cos(m\omega_0 t + \phi_{x,m\omega_0})\mathbf{x} + \epsilon_{m\omega_0}\sin(m\omega_0 t + \phi_{y,m\omega_0})\mathbf{y}]. \quad (1) $$

Here we assume laser field components polarized in x and y directions, x and y being the corresponding unit vectors, $e_{m\omega_0}$ define the amplitudes of the mth harmonics, $f_{m\omega_0}(t)$ define the slowly varying envelopes<fnref n="1" /> of the pulses, $0 \le \epsilon_{m\omega_0} \le 1$ determine the relative amplitudes and $\phi_{x,m\omega_0}$, $\phi_{y,m\omega_0}$ the carrier envelope phases of x and y components.

The generation of the free electron density is described by the rate equation

$$ \partial_t \rho_e(t) = W(t)[\rho_{at} - \rho_e(t)], \quad (2) $$

with $\rho_{at}$ being the neutral density. For simplicity, the ionization rate $W(t)$ can be evaluated from the quasi-static tunneling model [55],

$$ W(t) = \frac{\alpha}{|\mathbf{E}(t)|} e^{-\frac{\beta}{|\mathbf{E}(t)|}}, \quad (3) $$

with constants

$$ \alpha = 4 \frac{m_e^3 q^9}{(4\pi\epsilon_0)^5 \hbar^7} \left(\frac{U_1}{U_H}\right)^{\frac{5}{2}} \quad \text{and} \quad \beta = \frac{2}{3} \frac{m_e^2 q^5}{(4\pi\epsilon_0)^3 \hbar^4} \left(\frac{U_1}{U_H}\right)^{\frac{3}{2}} $$

characterizing the ionization potential. Here, $m_e$ and $q$ denote the electron mass and charge, respectively, and $U_H = 13.6$ eV is the Rydberg energy. In the numerical simulations below we assume helium (ionization potential $U_1 = 24.6$ eV) at atmospheric pressure and apply the limit $\rho_e \ll \rho_{at}$ in our analytics for technical convenience. For the parameter range used in our simulations, this is a rather good approximation since the typical ionization ratio $\rho_e/\rho_{at}$ is of the order of $10^{-4}$. We also note that, as soon as $\rho_e \ll \rho_{at}$, the polarization state of the THz radiation is nearly independent of the ionization rate (in contrast to the THz amplitude).

The electrons released acquire the drift velocity $\mathbf{v}(t, t') = (q/m_e) \int_{t'}^t \mathbf{E}(t'')dt''$, where $t'$ denotes their time of birth. For simplicity, we assume that electrons are born with zero drift velocity, and collisions are neglected in our analytical approach. The macroscopic current density $\mathbf{J}$ is then given by

$$ \mathbf{J}(t) = q \int_{-\infty}^{t} \mathbf{v}(t, t') [\partial_{t'} \rho_e(t')] dt'. \quad (4) $$

This net current density is responsible for the emission of the so-called Brunel radiation $\mathbf{E}_{Br}$, which can be estimated by using Jefimenko's theory [56] (see Appendix A) as

$$ \mathbf{E}_{Br}(t) = g \partial_t \mathbf{J} = g \frac{q^2}{m_e} \rho_e \mathbf{E}. \quad (5) $$

This key equation is often referred to as the local current (LC) model. The proportionality factor $g$ depends on the observation point and thus on the experimental setup. As explained in more detail in Appendix A, assuming a sufficiently large distance $R$ from the observer this factor expresses as $g \approx -\delta V/(4\pi\epsilon_0 c^2 R)$, where $c$ is the speed of light, $\epsilon_0$ is the vacuum permittivity and $\delta V$ is the emission volume. The LC model does not take into account any propagation effects of the pump. However, it allows for a simple and vivid analysis of the locally generated THz waveforms and predicts relative THz emission yields. It was successfully verified against simulations of the unidirectional pulse propagation equation in full (3D + 1)

<footnote n="1">¹In accordance with the usual slowly varying envelope approximation, the envelope functions $f_{m\omega_0}(t)$ must change slowly on time scales defined by the optical cycle of the fundamental.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Polarization ellipse in the laboratory coordinate system (x, y) and the principal axes coordinate system (ξ, η). The angle between x axis of the laboratory coordinate system and the major axis of the ellipse is denoted by Ψ. The ellipticity ϵ of the polarization state is defined as the ratio between the lengths of the minor and major axes (see text for details).</figcaption></figure>

dimensional geometry, cf. Refs. [19,38], and used to explain experimental observations on laser-generated THz radiation in [15,16,20,21,24,27,30,36,41,44,46,57–59]. Propagation effects in more complex geometries can be explained by superimposing LC contributions [25,26,37,52,53,60–62], provided that charge separation effects in the plasma can be neglected [63]. In such integral formulation the knowledge of the time derivative of the local electron current density already provides good estimates of the direction of THz emission and far-field spatial emission profiles in a realistic experiment with extended plasma volume. It is well established that two-color laser-induced elongated plasmas emit THz radiation mainly in the forward direction within a narrow conical distribution close to the propagation axis [16,52,62,64,65]. However, these characteristics will not be addressed here as they need complementary information on the spatial geometry of the emitting plasma. Note also that the LC model does not take into account any effects related to electron deflection by the tail of the atomic Coulomb potential [7].

## B. Polarization analysis

As we shall see later, the polarization states of the emitted THz pulses are quite unusual, since they exhibit an “ellipticity chirp” [see Fig. 1(b)]. To describe this property, instead of considering the pump field $E(t)$ and the Brunel radiation $\mathbf{E}_{\text{Br}}(t)$ in time, we operate with their Fourier components $\hat{\mathbf{E}}(\omega) = \frac{1}{2\pi} \int \mathbf{E}(t) \exp(i\omega t) dt$ and $\hat{\mathbf{E}}_{\text{Br}}(\omega)$.

Furthermore, we shall use the Jones formalism, see for example Ref. [66], to describe our THz waveforms. To a (real-valued) monochromatic field $\mathbf{E}(t)$ oscillating at frequency $\omega$, we assign a complex vector $\mathbf{j}$, which describes its polarization state. For instance, for the exemplary field $\mathbf{E}(t) = e_0(\cos \omega t, \epsilon \sin \omega t)$ we consider its complex representation $\propto e_0(e^{-i\omega t}, i\epsilon e^{-i\omega t})$, and $\mathbf{j}$ can be defined as the prefactors before the oscillating term, that is, for this example, $\mathbf{j} = e_0(1, i\epsilon)$. This complex vector is commonly represented as an “oriented ellipse” [see Fig. 2 and also Fig. 1(a)] whose principal axis is generally rotated with respect to the x axis of the laboratory coordinate system by some angle $\Psi$. The ellipticity of the

polarization state $\epsilon$ is defined as the ratio between the lengths of the minor and major axes, that is, $\epsilon = \beta/\alpha$ for $\beta \le \alpha$. In some works (see for instance [32]), the sign of $\epsilon$ specifies the direction of rotation of the electric field. Here, the direction of rotation is treated as a separate quantity $d$, taking values $d = \pm 1$, independent of $\epsilon$, so $\epsilon \ge 0$. The sign of $d$ is positive (clockwise rotation from the point of view of the emitter) or negative (resp. counterclockwise rotation). In the example above, we chose counterclockwise polarization. In this way, we can assign a polarization state to $\hat{\mathbf{E}}(\omega)$ and $\hat{\mathbf{E}}_{\text{Br}}(\omega)$ for each frequency $\omega$, thus obtaining a frequency-dependent ellipticity $\epsilon(\omega)$ and polarization angle $\Psi(\omega)$.

To assign some well-defined ellipticity to a THz pulse as a whole, we will also use the frequency-averaged ellipticity

$$ \langle \epsilon(\omega) \rangle = \frac{\int_{0}^{\omega_{\text{co}}} \epsilon(\omega) I(\omega) d\omega}{\int_{0}^{\omega_{\text{co}}} I(\omega) d\omega}, \quad (6) $$

where $\epsilon(\omega)$ denotes the frequency-dependent ellipticity and $I(\omega) \propto |\hat{\mathbf{E}}_{\text{Br}}(\omega)|^2$ is the THz spectral intensity, later extracted from the LC model. The cutoff frequency $\omega_{\text{co}}$ (in the upcoming simulations we took $\omega_{\text{co}} = \omega_0/4$) selects the frequency range of interest and thus gives a physical meaning to the above quantity. The frequency-averaged ellipticity is what one would expect from a typical measurement setup, for instance, a conventional scheme employing a linear rotating polarizer and frequency-averaged THz detection. In the same way, the averaged polarization angle is defined as

$$ \langle \Psi(\omega) \rangle = \frac{\int_{0}^{\omega_{\text{co}}} \Psi(\omega) I(\omega) d\omega}{\int_{0}^{\omega_{\text{co}}} I(\omega) d\omega}, \quad (7) $$

where $\Psi(\omega)$ is the frequency dependent angle (see Fig. 2).

## III. TWO-COLOR CASE: EXEMPLARY WAVE SHAPES

In the following, we will frequently refer to the conventional two-color pump, that is, when $m = 1, 2$ in Eq. (1). In this case, when assuming two Gaussian components with identical duration $\tau$ and neglecting carrier-envelope-phase (CEP) effects [67] (that is, setting the CEP phase of the fundamental harmonic to zero), the pump electric field Eq. (1) reduces to

$$ \mathbf{E} = \mathbf{E}_{\omega_0} + \mathbf{E}_{2\omega_0}, \quad (8) $$

$$ \mathbf{E}_{\omega_0} = e_{\omega_0} e^{-t^2/\tau^2} [\cos \theta_1 \cos(\omega_0 t) \mathbf{x} + \sin \theta_1 \sin(\omega_0 t) \mathbf{y}], \quad (9) $$

$$ \mathbf{E}_{2\omega_0} = e_{2\omega_0} e^{-t^2/\tau^2} [\cos \theta_2 \cos(2\omega_0 t + \phi)\mathbf{x} + \sin \theta_2 \sin(2\omega_0 t + \phi)\mathbf{y}]. \quad (10) $$

Here $\mathbf{E}_{\omega_0}$ and $\mathbf{E}_{2\omega_0}$ are electric fields of FH and SH pulses, respectively, with $\phi$ being their relative phase. $\theta_1$ and $\theta_2$ determine the ellipticities of each harmonic as $\epsilon_{j\omega_0} = |\tan \theta_j|$, $j = 1, 2$. A possible experimental realization is described in more detail in the next paragraph. For the particular numerical simulations below we assume $\tau = 50$ fs, $\omega_0$ corresponding to 800 nm wavelength, and $e_{\omega_0} = e_{2\omega_0} = 0.053$ a.u. In this paper, we use Hartree atomic units for the electric field, 1 a.u. $\approx 514.22$ GV/m.

A practical setup capable of producing of such wave shapes is shown in Fig. 3. The beams are prepared with initial
<figure><img src="image_4.png" /><figcaption>FIG. 3. Controlling pump wave shapes in the two-color case, consisting of fundamental (FH) and second (SH) harmonics. (a) The polarization directions of the incoming linearly polarized FH and SH pulses are rotated by the angles $\theta_1$ and $\theta_2$ with respect to the x axis. The light and bold traces schematically represent the waveforms before and after this rotation, respectively, with orange representing the FH and blue the SH. Then, the harmonics pass through quarter wave plates (QWPs), which are alligned with their fast and slow axes along x and y direction. (b), (c) By varying $\theta_2$, the SH polarization state (after QWP) changes and follows the path on the Bloch sphere shown in (b). In (c), the SH ellipticity versus $\theta_2$ is drawn. The starting point of this path and its direction are indicated by the black dot with a black arrow in (b), (c). Black dots on the Bloch sphere (without arrows) indicate the positions of south and north pole, respectively, where the SH is circularly polarized. Background color in (b), (c) describes the field rotation direction for the corresponding area, whereas the color of the line segments denotes the polarization direction [see legend in (c)]. The dashed lines and numbers 1, 2, 3, 4a, and 4b indicate several characteristic pump configurations to be further treated.</figcaption></figure>

polarization angles $\theta_1$ and $\theta_2$ (with respect to the x axis) for FH and SH, correspondingly. Before the pump is focused, each component passes through a quarter wave plate (QWP) with fast and slow axis aligned along x and y direction. The relative phase $\phi$ is typically controlled by changing the delays of FH relative to SH. By superimposing FH and SH fields, we obtain a large diversity of possible pump waveforms. Our baseline pump configurations when $\theta_1$ is fixed to 0 (x-linearly polarized FH pulse) or $\pi/4$ (clockwise circularly polarized FH pulse) are displayed in Fig. 4 over one optical cycle for relative phase $\phi = 0$ and $\phi = \pi/2$. For these FH configurations, Fig. 4 shows Lissajou figures, describing the electric field vector trajectories for several exemplary values of $\theta_2$, when FH and SH have equal amplitudes. Note that for the circularly polarized FH pulse with $\theta_1 = \pi/4$ we show only the $\phi = 0$ configuration because changing the relative phase $\phi$ causes in this case only a trivial rotation of the Lissajou figure [see also Figs. 8(a) and 8(b)]. Exemplary SH configurations are marked by numbers 1, 2, 3, 4, 4a, and 4b, indicated also in Fig. 3(c), and refer to different polarization states of the SH, namely

(1) linear polarization in x direction ($\theta_2 = 0$),

(2) circular polarization (clockwise, $\theta_2 = \pi/4$),

(3) linear polarization in y direction ($\theta_2 = \pi/2$),

(4a) circular polarization (counterclockwise, $\theta_2 = 3\pi/4$),

(4b) elliptically polarized in y direction ($\theta_2 \approx 0.6\pi$).

Note that following Eq. (10), changing $\theta_2$ for given $\theta_1$ and $\phi$ accordingly modifies the effective relative phase between FH and SH.

These baseline pump field configurations of SH can be connected to form a single smooth curve on the Bloch sphere [see the corresponding curve in Figs. 3(b) and 3(c)]. The respective parameters of the Bloch sphere are $S_1 = \cos\tilde{\theta}\cos\tilde{\phi}$, $S_2 = \cos\tilde{\theta}\sin\tilde{\phi}$, $S_3 = -\sin\tilde{\theta}$, where $\tilde{\theta} = 2\theta_2$ and $\tilde{\phi} = \pi[1 - \sin(\cos\tilde{\theta})]/2$ describe the angles on the Bloch sphere. We numerically integrated the LC model Eq. (5) recalled in

Sec. II A for SH pump field configurations along the above mentioned curve. The resulting properties of the generated THz radiation are summarized in Fig. 5, where the first column corresponds to a circularly polarized FH pump field and the second and third columns to a linearly polarized FH and relative phase $\phi = 0$ and $\phi = \pi/2$, respectively. In Figs. 5(a)–5(c) the THz pulse energies defined by the integral of $I(\omega)$ over all frequencies below the cutoff value $\omega_{c0} = \omega_0/4$ are shown. The frequency-averaged THz polarization angles and ellipticities defined by Eqs. (7) and (6) in the same frequency range are detailed in Figs. 5(d)–5(f) and 5(g)–5(i), respectively. Some corresponding pump polarizations are shown in Fig. 4. The THz energy reaches its maximum when SH and FH are circularly polarized and corotating, and we will refer to such configuration as “CP-S”, in agreement with recent experimental and theoretical investigations [24,29,32]. The frequency-averaged ellipticity behaves inversely to the THz energy, i.e. the ellipticity is high when the THz energy is low and vice versa. Maximum THz ellipticity ($\sim 0.6$) is reached for linearly x-polarized FH, relative phase $\phi = 0$ or $\phi = \pi/2$, and SH departing from linear x or y polarization by an angle $\Delta\theta_2 < \pm\pi/10$, respectively. Comparable behaviors have been reported in Ref. [31]. However, note that in this work propagation effects inside the filament play an important role.

Furthermore, typical THz wave shapes obtained for 50 fs long pump pulses, also calculated using the LC model, are shown in Fig. 6 for different values of $\theta_2$, assuming equal amplitudes of FH and SH, $\phi = 0$, and $\theta_1$ to be 0 or $\pi/4$. We can observe that the THz wave shapes for these parameters develop within a single, nearly half-cycle spike, which is a typical signature for this photoionization-induced generation mechanism. For longer pump pulse durations, the ellipticity of the THz radiation decreases, as can be observed from Fig. 7.
<figure><img src="image_5.png" /><figcaption>FIG. 4. Pump wave shapes over one optical cycle of the FH for different configurations of QWPs, indicated by numbers in Fig. 5 for (a) circular polarization with clockwise rotation ($\theta_1 = \pi/4$) and (b), (c) linear polarization in x direction ($\theta_1 = 0$) of the FH, assuming (a), (b) $\phi = 0$ and (c) $\phi = \pi/2$. FH and SH components have equal amplitudes. The red dots represent instants $t_n$ of the maxima of $|\mathbf{E}(t)|$, the arrows indicate the product $\mathbf{A}(t_n)\delta\rho_n$, involving the $n$th density step of generated free electrons $\delta\rho_n$ and the vector potential $\mathbf{A}(t_n)$ at these instants. For better visibility we omitted some extrema and rescaled some arrows.</figcaption></figure>

To explain these observations and gain a general understanding of how THz wave shapes are formed, we approach the problem analytically in the following and initially focus on the limit $\omega \to 0$. Furthermore, we mainly discuss the ellipticity of the emitted THz radiation (rather than particular wave shapes).

#### IV. THZ WAVE SHAPES IN THE ZERO FREQUENCY LIMIT

We first consider the polarization properties of $\mathbf{E}_{Br}$ in the zero frequency limit $\omega \to 0$. We note that the exact zero frequency has, strictly speaking, no physical relevance for propagating electromagnetic fields, and the amplitude of the emitted radiation should vanish as $\omega \to 0$ (see Ref. [19] and Appendix A). However, its (formal) polarization properties are easily accessible and indicate what can be expected for finite frequencies close to zero.

Assuming ionization far below saturation and neglecting the effects of the field envelope, we can expand the pump electric field into harmonics of the fundamental frequency $\omega_0$:

$$ \mathbf{E}(t) = \sum_{m} \mathbf{E}_{m} e^{-im\omega_{0}t}, \quad (11) $$

where we assume $\mathbf{E}_0 = 0$ in what follows (no DC component in the pump). Note that because $\mathbf{E}(t)$ is real-valued, $\mathbf{E}_m = \mathbf{E}_{-m}^*$ holds, where * denotes the complex conjugate.

Similarly, the ionization rate $W(t)$ can be written as a sum over harmonics $l\omega_0$,

$$ W(t) = \sum_{l} W_{l} e^{-il\omega_{0}t}. \quad (12) $$

Then, according to Eq. (2), which simplifies to

$$ \partial_t \rho_e \approx \rho_{at} W, \quad (13) $$

the time derivative of the electron density $\partial_t \rho_e$ is periodic and can also be written as a Fourier series. Thus, the electron density $\rho_e(t)$ reads

$$ \rho_e(t) \approx \sum_l \rho_l e^{-il\omega_0 t} + \rho_{at} W_0 t, \quad (14) $$

and the coefficients $\rho_l$ are given by

$$ \rho_l = i \frac{\rho_{at}}{\omega_0 l} W_l, \quad l \neq 0. \quad (15) $$

The coefficient $\rho_0$ can be formally computed as $\rho_0 = \rho_e(t = 0) - \sum_{l \neq 0} \rho_l$. We note that the nonperiodic term in Eq. (14) conflicts with our initial assumption $\rho_e \ll \rho_{at}$. This inconsistency comes from neglecting the finite pulse duration and using Fourier series expansion. However, the nonperiodic term and the coefficient $\rho_0$ do not play any role in the following considerations.
<figure><img src="image_6.png" /><figcaption>FIG. 5. (a), (b), (c) THz energy, (d), (e), (f) frequency-averaged THz polarization angle, and (g), (h), (i) frequency-averaged THz ellipticity as a function of the rotation angle $\theta_2$ for SH considering (a), (d), (g) circular polarization for FH, and (b), (e), (h) linear polarization for FH with $\phi = 0$ and (c), (f), (i) linear polarization for FH with $\phi = \pi/2$. The vertical lines indicate selected values of $\theta_1$, $\theta_2$, [marked by numbers in (d), (e), (f)] for which the respective pump wave shapes are shown in Fig. 4. The other parameters of the pump pulse as well as the parameters of the medium are given below Eqs. (10) and (3) correspondingly.</figcaption></figure>

With Eqs. (5), (11), and (14) we have

$$ 
\begin{aligned}
\mathbf{E}_{\text{Br}} &= g \frac{q^2}{m_e} \rho_e \mathbf{E} \\
&= g \frac{q^2}{m_e} \left( \sum_{m,l} \rho_l e^{-il\omega_0 t} \mathbf{E}_m e^{-im\omega_0 t} + \rho_{\text{at}} W_0 t \mathbf{E} \right).
\end{aligned} 
\quad (16) $$

In this section, we are interested in the zero frequency limit. Therefore, only summands with $m = -l$ contribute:

$$ \mathbf{E}_{\text{Br},0} = g \frac{q^2}{m_e} \sum_{m \neq 0} \mathbf{E}_m \rho_{-m}. \quad (17) $$

Here, we can exclude $m = 0$ from the summation because $\mathbf{E}_0 = 0$, and therefore the coefficient $\rho_0$ does not contribute to $\mathbf{E}_{\text{Br},0}$. Only the harmonics of the fundamental frequency $\omega_0$ contribute to $\mathbf{E}_{\text{Br},0}$, which is why we call Eq. (17) “frequency representation” in what follows.

Since ionization takes place on the subcycle scale of the optical driver near the maxima of the pump electric field, the free electron density grows in the form of “sharp steps,” [19]:

$$ \rho_e(t) \approx \sum_n \delta \rho_n \Theta(t - t_n), \quad (18) $$

where $\Theta(t)$ is the Heaviside step function, $t_n$ are the positions of the ionization events, and $\delta \rho_n$ are the step sizes. The time

derivative of the electron density then reads

$$ \partial_t \rho_e(t) \approx \sum_n \delta \rho_n \delta(t - t_n), \quad (19) $$

where $\delta(t)$ denotes the Dirac $\delta$ function. For the assumed pump electric field configuration, the sequence of ionization steps $t_n$ is periodic with period $2\pi/\omega_0$. Assuming $N_c$ ionization events per optical cycle, we can write

$$ \partial_t \rho_e(t) \approx \sum_{n=1}^{N_c} \delta \rho_n \sum_j \delta(t - t_n + 2\pi j / \omega_0) \quad (20) $$

as a sum of $N_c$ Dirac combs with period $2\pi/\omega_0$. The Fourier series representation of one Dirac comb reads

$$ \sum_j \delta(t - t_n + 2\pi j / \omega_0) = \frac{\omega_0}{2\pi} \sum_l e^{-il\omega_0(t - t_n)}. \quad (21) $$

Therefore, we find the Fourier series representation for the time derivative of the electron density as

$$ \partial_t \rho_e(t) = \sum_l \left( \frac{\omega_0}{2\pi} \sum_{n=1}^{N_c} \delta \rho_n e^{il\omega_0 t_n} \right) e^{-il\omega_0 t_n}. \quad (22) $$

With Eqs. (12) and (13) we can conclude that

$$ W_l = \frac{\omega_0}{2\pi \rho_{\text{at}}} \sum_{n=1}^{N_c} \delta \rho_n e^{il\omega_0 t_n}, \quad (23) $$
<figure><img src="image_7.png" /><figcaption>FIG. 6. Selected THz wave shapes for (a), (b) circular and (c), (d) linear polarization of FH for different values of $\theta_2$. The left and right panels show the y and x components of the THz fields versus time, respectively. The bottom panels show the resulting Lissajou figures of the electric field vector. The other parameters of the pump pulse, as well as the parameters of the medium, are given after Eqs. (10) and (3) correspondingly.</figcaption></figure>

and Eq. (15) yields

$$ \rho_l = \frac{i}{2\pi l} \sum_{n=1}^{N_c} \delta\rho_n e^{il\omega_0 t_n}, \quad l \neq 0. \qquad (24) $$

Thus, with Eq. (17) we arrive to

$$ \mathbf{E}_{\text{Br},0} = -\frac{igq^2}{2\pi m_e} \sum_{m \neq 0} \sum_{n=1}^{N_c} \frac{\delta\rho_n \mathbf{E}_m}{m} e^{-im\omega_0 t_n}. \qquad (25) $$

In contrast to Eq. (17), this equation provides a “mixed representation” of the process, where both information in time domain $t_n$, $\delta\rho_n$ and in frequency domain $\mathbf{E}_m$ is used. We will use this representation extensively in our following analytics.

To further understand the physical meaning of the process, one can rewrite it in “time representation” (cf. Ref. [61]) as

$$ \mathbf{E}_{\text{Br},0} \propto \sum_{n=1}^{N_c} \delta\rho_n \mathbf{A}(t_n), \qquad (26) $$

where

$$ \mathbf{A}(t) = - \int_{-\infty}^{t} \mathbf{E}(t') dt' \qquad (27) $$

is the vector potential (see Appendix B for details). According to Eq. (26), $\mathbf{E}_{\text{Br},0}$ is linearly polarized and parallel to the net current $\mathbf{J}(t \to +\infty)$ formed by the ionized electrons after the pulse has interacted with the gas. Furthermore, Eq. (26) allows one to predict which pump wave shapes generate THz radiation effectively [19]: for each ionization event, one must optimize $\delta\rho_n$ and $\mathbf{A}(t_n)$ simultaneously, which is not trivial: Since the $t_n$ are located at the extrema of $\mathbf{E}(t)$, for simple (e.g., single-color) pump wave shapes, $|\mathbf{A}(t_n)| \approx 0$, rendering $\delta\rho_n \mathbf{A}(t_n)$ negligible. Examples of two-color pump wave shapes with $\delta\rho_n \mathbf{A}(t_n)$ for each ionization event (red arrows) are shown in Fig. 4. Comparing, for example, the cases Fig. 4(a).2 and Fig. 4(b).2, one can see that although the amplitude of $\mathbf{E}(t_n)$ is the same and thus $\delta\rho_n$ is similar, $\delta\rho_n \mathbf{A}(t_n)$ is larger in case 4(a).2 due to the larger amplitude of $\mathbf{A}(t_n)$. For certain pump configurations, there is more than one significant ionization event per cycle, e.g., for case 4a in Fig. 4(a). In general, if several ionization events with comparable amplitudes are present, they tend to cancel each other out, at least partially, since the vector potentials $\mathbf{A}(t_n)$ are pointing in different directions. As a result, the maximum THz yield is achieved when one ionization event per cycle dominates. For a circularly polarized two-color pump pulse with corotating SH and FH [see Fig. 4(a).2], this property clearly
<figure><img src="image_8.png" /><figcaption>FIG. 7. Frequency-averaged THz ellipticities $\langle \epsilon \rangle$ according to Eq. (6) as a function of the pump pulse duration $\tau$ [see Eqs. (8)–(10)] for $\phi = 0$: (a) corotating circularly polarized FH and SH with equal amplitudes (CP-S), and (b) linearly polarized FH and circularly polarized SH. The other parameters of the pump pulse as well as the parameters of the medium are given after Eqs. (10) and (3) correspondingly.</figcaption></figure>

appears and the large THz yield is confirmed in Fig. 5(a) [24,29,32].

For the particular case of a two-color pump, the polarization of the THz field is determined by the harmonics $m = 1, 2$. In the “mixed representation” of Eq. (25) we find

$$ \mathbf{E}_{\text{Br},0} \propto \sum_{n=1}^{N_c} \delta\rho_n \left( \mathbf{E}_1 e^{-i\omega_0 t_n} - \mathbf{E}_{-1} e^{i\omega_0 t_n} + \frac{\mathbf{E}_2}{2} e^{-2i\omega_0 t_n} - \frac{\mathbf{E}_{-2}}{2} e^{2i\omega_0 t_n} \right). \quad (28) $$

To see how this “mixed representation” can be used, let us consider the particular waveshape Eqs. (8)–(10) with $\epsilon_{\omega_0} = \epsilon_{2\omega_0} = 1$ and $\phi = 0$. This is the well-known CP-S wave shape, with both FH and SH being circularly polarized and corotating [7,24,32,34]. Here, we consider this pump wave shape through our analytical approach. In Fig. 8, where FH and SH have identical amplitudes, there is only one “main” maximum per optical cycle [cf. Fig. 4(a).2], which appears at $t_1 = 0$ for our choice of $\phi = 0$. As before, we neglect the effects of the field envelope in the analytical treatment and consider

plane waves by taking the limit $\tau \to +\infty$. Then, the pump harmonics amplitudes read

$$ \mathbf{E}_1 = \frac{e_{\omega_0}}{2}(1, +i), \quad \mathbf{E}_{-1} = \frac{e_{\omega_0}}{2}(1, -i), \quad (29) $$

$$ \mathbf{E}_2 = \frac{e_{2\omega_0}}{2}(1, +i), \quad \mathbf{E}_{-2} = \frac{e_{2\omega_0}}{2}(1, -i), \quad (30) $$

where the right-hand sides represent Jones vectors with the corresponding x and y components. Substituting Eqs. (29) and (30) into Eq. (25) gives

$$ \mathbf{E}_{\text{Br},0} \propto (0, 1). \quad (31) $$

Thus, the 0th harmonic for this pump configuration is linearly polarized in the y direction. More generally, for an arbitrary phase $\phi$ between the corotating circularly polarized FH and SH pump fields, the “main” maximum of the electric field shifts from $t_1 = 0$ to $t_1 = -\phi/\omega_0$, and Eq. (30) has to be modified into

$$ \mathbf{E}_2 = \frac{e_{2\omega_0}}{2}(1, +i)e^{-i\phi}, \quad \mathbf{E}_{-2} = \frac{e_{2\omega_0}}{2}(1, -i)e^{i\phi}. \quad (32) $$

Then, using common trigonometric transformations, it is easy to see that the entire wave shape is rotated by $-\phi$, giving the polarization of the THz field as

$$ \mathbf{E}_{\text{br},0} \propto -(\sin \phi, \cos \phi). \quad (33) $$

This well-known property [24,27,32,34] is illustrated in Figs. 8(a)–8(c). In the next section, more general examples are addressed.

## V. THZ WAVE SHAPES AT NON ZERO FREQUENCIES

As shown in the previous section, the THz polarization at $\omega \to 0$ is always linear. This is of course not surprising, since no temporal dynamics is possible at zero frequency. The orientation of the polarization depends on the final direction of the net current density, which in turn depends on the pump waveshape. For nonzero frequencies, the situation is more complicated. In Fig. 9, we present the spectrum, ellipticity, and polarization angle of the THz radiation computed for an exemplary pump pulse of finite duration using the LC model Eq. (5). The remarkable feature is that, while the polarization angle remains constant with frequency in first-order approximation ($\omega/\omega_0 < 0.2$), the THz ellipticity increases linearly, revealing a “linear chirp of ellipticity”. Thus, the THz radiation produced by the two-color mechanism may have rather unusual polarization states. If each frequency component of the THz pulse has a different ellipticity, the resulting field in the time domain may become very complex. However, the ellipticity changes with frequency in a controlled, linear way. As we show later, this linear chirp of ellipticity is a generic feature of ionization-based THz generation.

In the next subsections, we derive the general expressions for the polarization state for arbitrary pump wave shape for (low) nonzero frequencies, and present closed analytical expressions for the linear ellipticity chirp for some selected pump configurations.

### A. THz polarization state for generic multicolor pump pulses

For simplicity, we again assume $\rho_e \ll \rho_{at}$ as in the previous section, but now we have to consider the finite duration of
<figure><img src="image_9.png" /><figcaption>FIG. 8. (a), (b), (d)–(f) Pump wave shapes (black lines) for (a), (b) the CP-S configuration and (d)–(f) linearly polarized FH with circularly polarized SH for an exemplary pulse duration of 50 fs and different phases $\phi$ between FH and SH. The associated THz wave shapes are shown by red lines (their amplitudes are rescaled to make them comparable to the pump fields). (c) Numerically-computed THz wave shapes for the CP-S case for various values of $\phi$ (solid lines) compared to their corresponding analytical solution of Eq. (33) (dashed lines). The other parameters of the pump pulse as well as the parameters of the medium are given below Eq. (10) and Eq. (3) correspondingly.</figcaption></figure>

the pump pulse. In this situation, the pump wave shape is no longer periodic. Therefore, in this section, $N$ denotes the (finite) total number of ionization events in the pulse and must not be confused with the number of ionization events per optical cycle $N_c$ used before. The radiated field $E_{Br}(t) \propto \rho_e(t)E(t)$ expresses as a convolution integral in the Fourier domain,

$$ \hat{E}_{Br}(\omega) = g \frac{q^2}{m_e} \int_{-\infty}^{\infty} \hat{\rho}_e(\omega - \Omega) \hat{E}(\Omega) d\Omega, \quad (34) $$

where the $\hat{\omega}$ symbol refers to the Fourier transform in time. From Eq. (19) we immediately get

$$ -i\omega\hat{\rho}_e(\omega) \approx \frac{1}{2\pi} \sum_{n=1}^{N} \delta\rho_n e^{i\omega t_n}, \quad (35) $$

which implies that $\hat{\rho}_e(\omega)$ is rather broad in frequency space. In contrast, we assume that the multicolor pump field $\hat{E}(\omega)$ is given by narrow peaks at the harmonic frequencies, that is, the corresponding pulse envelopes are slowly varying in time.

<figure><img src="image_10.png" /><figcaption>FIG. 9. (a) Spectral amplitude $|E_{Br}(\omega)|$ of THz radiation, (b) invariant orientation of the THz polarization ellipse and (c) THz spectral ellipticity versus frequency considering right-handed circular polarization for SH and linear or right-handed circular polarization for FH (see legend), numerically computed from the LC model Eq. (5). The other parameters of the pump pulse as well as the parameters of the medium are given below Eq. (10) and Eq. (3) correspondingly. Here $\phi = 0$.</figcaption></figure>
These narrow peaks will be approximated by Dirac δ functions in the following, and we write

$$ \hat{\mathbf{E}}(\omega) \approx \sum_m \mathbf{E}_m \delta(\omega - m\omega_0), \quad (36) $$

setting $\mathbf{E}_0 = 0$. We note that the approximation Eq. (36) renders the pump field similar to Eq. (11) used in the previous section. However, we must remember that Eq. (36) is only an approximation to facilitate further computations. We could replace in Eq. (36) all $\delta(\omega)$ with smooth narrow peaks $\tilde{\delta}_m(\omega)$ obeying the normalization $\int \tilde{\delta}_m(\omega)d\omega = 1$, provided that their widths are much smaller than $\omega_0$. In the case of a Gaussian pulse with duration $\tau$, one finds $\tilde{\delta}_m(\omega) \propto \tau \exp[-\tau^2(\omega - m\omega_0)^2/4]$ for the spectrum of the mth harmonic, whose amplitude and direction are defined by $\mathbf{E}_m$. However, the following integrals over the frequency would then have to be treated numerically, while the analytical results using $\delta(\omega)$ already provide quite good approximations.

With Eqs. (34)–(36), the Brunel radiation in the THz range can be expressed as

$$ \begin{aligned} \hat{\mathbf{E}}_{\text{Br}}(\omega) &= ig \frac{q^2}{2\pi m_e} \sum_{n=1}^{N} \delta\rho_n \int \hat{\mathbf{E}}(\Omega) \frac{e^{i(\omega-\Omega)t_n}}{\omega-\Omega} d\Omega \\ &= ig \frac{q^2}{2\pi m_e} \sum_{m} \sum_{n=1}^{N} \delta\rho_n \mathbf{E}_m \frac{e^{i(\omega-m\omega_0)t_n}}{\omega-m\omega_0}. \end{aligned} \quad (37) $$

This is, in general, a complicated function of $\omega$. However, for small $\omega$ one can approximate this expression as

$$ \hat{\mathbf{E}}_{\text{Br}}(\omega) = ig \frac{q^2}{2\pi m_e} \left( \boldsymbol{\mu}_0 + \boldsymbol{\mu}_1 \frac{\omega}{\omega_0} + \dots \right), \quad (38) $$

where the vector-valued coefficients $\boldsymbol{\mu}_0$ and $\boldsymbol{\mu}_1$ are

$$ \boldsymbol{\mu}_0 = - \sum_m \sum_{n=1}^{N} \delta \rho_n \mathbf{E}_m \frac{e^{-im\omega_0 t_n}}{m\omega_0}, \quad (39) $$

$$ \boldsymbol{\mu}_1 = - \sum_m \sum_{n=1}^{N} \delta \rho_n \mathbf{E}_m \frac{e^{-im\omega_0 t_n} (1 + im\omega_0 t_n)}{m^2 \omega_0}. \quad (40) $$

With these expressions, it is possible to compute the ellipticity $\epsilon(\omega)$ of the THz polarization ellipse at each frequency.

## B. THz ellipticity $\epsilon(\omega)$ for the multicolor case

With zero ellipticity at $\omega \to 0$ and when neglecting higher order terms in $\omega$, Eq. (38) predicts a linearly increasing ellipticity with $\omega$:

$$ \epsilon(\omega) \approx \mathcal{B} \frac{\omega}{\omega_0}, \quad \omega \ll \omega_0, \quad (41) $$

that is, the ellipticity is linearly chirped in frequency, and we call $\mathcal{B}$ the “ellipticity chirp”. To obtain an analytical expression for $\mathcal{B}$, let us recall that the electric field satisfies $\mathbf{E}_m = \mathbf{E}_{-m}^*$ and $\mathbf{E}_0 = 0$, so that we can rewrite Eqs. (39) and

(40) as

$$ \boldsymbol{\mu}_0 = -2i \sum_{n=1}^{N} \delta \rho_n \sum_{m>0} \Re \left[ \frac{\mathbf{E}_m}{m\omega_0} e^{-im\omega_0 t_n} \right], \quad (42) $$

$$ \boldsymbol{\mu}_1 = -2 \sum_{n=1}^{N} \delta \rho_n \sum_{m>0} \Re \left[ \frac{\mathbf{E}_m (1 + im\omega_0 t_n) e^{-im\omega_0 t_n}}{m^2 \omega_0} \right]. \quad (43) $$

This shows that, $\boldsymbol{\mu}_0$ is purely imaginary, while $\boldsymbol{\mu}_1$ always takes real values.

To compute the ellipticity $\epsilon(\omega)$, we have to evaluate the orthogonal component of $\boldsymbol{\mu}_1$ with respect to $\boldsymbol{\mu}_0$. This amounts to computing the cross-product of both vectors normalized to the reference vector $\boldsymbol{\mu}_0$, namely,

$$ \mathcal{B} = \frac{|\boldsymbol{\mu}_0 \times \boldsymbol{\mu}_1|}{|\boldsymbol{\mu}_0|^2}, \quad (44) $$

or equivalently:

$$ \mathcal{B} = \frac{|\mu_{1y}\mu_{0x} - \mu_{1x}\mu_{0y}|}{|\mu_{0x}|^2 + |\mu_{0y}|^2}. \quad (45) $$

For elliptically polarized multicolor pump fields, $\mathcal{B} > 0$ and therefore nonzero frequency components of the generated THz pulse are generally elliptically polarized.

## C. THz ellipticity $\epsilon(\omega)$ for the two-color case

Here we consider a vectorial electric field for the incident two-color pump pulse given by Eqs. (8)–(10). We simplify the SH component to circular polarization by setting $\epsilon_{2\omega_0} = 1$. Still, the FH component may have an arbitrary ellipticity $\epsilon_{\omega_0}$, and we allow for an arbitrary phase and ratio of amplitudes between SH and FH. Furthermore, we ignore the contributions of the Gaussian envelope in our analytical calculations. As illustrative examples, we will address two polarization geometries: A CP-S polarization state (corotating FH and SH) where $\epsilon_{\omega_0} = 1$, and a $\mathbf{x}$ linearly polarized FH coupled to a circularly polarized SH where $\epsilon_{\omega_0} = 0$. Examples of pump and numerically calculated THz wave shapes for different values of the phase $\phi$ between FH and SH are shown in Fig. 8(a-c) for the former (CP-S) and Figs. 8(d)-8(f) for the latter case.

In the Jones representation, such an input field is expressed as

$$ \mathbf{E}_1 = \frac{e_{\omega_0}}{2} (1, +i\epsilon_{\omega_0}), \quad \mathbf{E}_{-1} = \frac{e_{\omega_0}}{2} (1, -i\epsilon_{\omega_0}), \quad (46) $$

$$ \mathbf{E}_2 = \frac{e_{2\omega_0}}{2} (1, +i)e^{-i\phi}, \quad \mathbf{E}_{-2} = \frac{e_{2\omega_0}}{2} (1, -i)e^{+i\phi}. \quad (47) $$

Evaluating Eqs. (39) and (40) yields

$$ \boldsymbol{\mu}_0 = -i \sum_{n=1}^{N} \frac{\delta \rho_n}{\omega_0} \left[ e_{\omega_0} \left( -\sin(\omega_0 t_n) \right) + \frac{e_{2\omega_0}}{2} \left( -\sin(2\omega_0 t_n + \phi) \right) \right], \quad (48) $$

$$ \boldsymbol{\mu}_1 = - \sum_{n=1}^{N} \frac{\delta \rho_n}{\omega_0} \left[ e_{\omega_0} \left( \cos(\omega_0 t_n) + \omega_0 t_n \sin(\omega_0 t_n) \right) + \frac{e_{2\omega_0}}{4} \left( \cos(2\omega_0 t_n + \phi) + 2\omega_0 t_n \sin(2\omega_0 t_n + \phi) \right) \right]. \quad (49) $$
The instants of ionization $t_n$ depend on the polarization geometry of the input field. They correspond to maxima of the absolute value of the pump electric field $|E|$. For our two-color pulse, one has

$$|E(t)|^2 = e_{\omega_0}^2 \left( \frac{1 + \epsilon_{\omega_0}^2}{2} + r^2 + r(1 + \epsilon_{\omega_0}) \cos(\omega_0 t + \phi) + \frac{1 - \epsilon_{\omega_0}^2}{2} \cos(2\omega_0 t) + r(1 - \epsilon_{\omega_0}) \cos(3\omega_0 t + \phi) \right), \quad (50)$$

where $r \equiv e_{2\omega_0}/e_{\omega_0}$. The ionization instants $t_n$ are determined by the roots of $\partial_t |E(t)|^2 = 0$, from which we have to select those that satisfy $\partial_t^2 |E(t)|^2 < 0$.

In CP-S geometry ($\epsilon_{\omega_0} = 1$) one has to solve $\sin(\omega_0 t + \phi) = 0$. The ionization instants are given by

$$\omega_0 t_n = 2n\pi - \phi. \quad (51)$$

For a x linearly polarized FH ($\epsilon_{\omega_0} = 0$), the ionization instants proceed from the roots of

$$\frac{\partial}{\partial t} \{\cos(\omega_0 t)[\cos(\omega_0 t) + 2r \cos(2\omega_0 t + \phi)]\} = 0. \quad (52)$$

Solving this equation is difficult for an arbitrary $\phi$, and in this paper only the limit $r \to 1$ is considered. Yet, there exist exact solutions for particular values of $\phi$, namely,

$$\phi = 0 \quad \rightarrow \quad \omega_0 t_n = 2n\pi, \quad (53)$$

$$\phi = \frac{\pi}{4} \quad \rightarrow \quad \omega_0 t_n = \arctan \left( \frac{\sqrt{2} - \sqrt{6}}{\sqrt{2} + \sqrt{6}} \right) + 2n\pi, \quad (54)$$

$$\phi = \frac{\pi}{2} \quad \rightarrow \quad \omega_0 t_n = n\pi - (-1)^n \frac{\pi}{6}. \quad (55)$$

Whereas for $\phi = 0$ and $\phi = \pi/4$ there is only one strong ionization event over the cycle, for the case $\phi = \pi/2$ there are two events with equal amplitudes, see insets in Fig. 10. We note that the influence of the phase angle $\phi$ on the laser-to-THz conversion efficiency was studied extensively in Refs. [15,18,19,42,57]. Here, we complement these findings by its impact on the THz waveforms and ellipticities. In the following, we consider two distinct situations: zero or nonzero relative phase $\phi$.

1. Zero relative phase

For $\phi = 0$, Eqs. (51) and (53) show that ionization events occur at $t_n = 0$ modulo $2\pi/\omega_0$. In the CP-S geometry ($\epsilon_{\omega_0} = 1$, $\epsilon_{2\omega_0} = 1$), one finds from Eqs. (48) and (49):

$$\mu_0 = -\frac{1}{\omega_0} \sum_{n=1}^{N} \delta\rho_n \left( e_{\omega_0} + \frac{e_{2\omega_0}}{2} \right) \binom{0}{i}, \\ \mu_1 = -\frac{1}{\omega_0} \sum_{n=1}^{N} \delta\rho_n \left( e_{\omega_0} + \frac{e_{2\omega_0}}{4} - 2n\pi \left[ e_{\omega_0} + \frac{e_{2\omega_0}}{2} \right] \right), \quad (56)$$

and the slope $B$, cf. Eq. (45), is given by

$$B = \left| \frac{\mu_{1x}}{\mu_{0y}} \right| = \left| \frac{e_{\omega_0} + e_{2\omega_0}/4}{e_{\omega_0} + e_{2\omega_0}/2} \right|. \quad (57)$$

<figure><img src="image_11.png" /><figcaption>FIG. 10. Numerical results obtained from the LC model (lines and circles): (a) Relative THz intensities, (b) values of ellipticity chirp factor $B$ and (c) polarization angle $\Psi$ [see Fig. 2 and Eq. (7)], as functions of the relative phase $\phi$ between FH and SH, for FH being linearly and SH circularly polarized. Black stars show analytical predictions. Insets show the pump wave shapes for the same values of $\phi$ as in Figs. 8(d)–8(f) [$\phi = 0$ in (a), $\phi = \pi/4$ in (b) and $\phi = \pi/2$ in (c)], with arrows indicating the electric field vectors for the dominant ionization events $t_n$; shaded region shows the “extent” of the ionization events taken as the full-width-at-half-maximum of the ionization rate peaks. Note the asymmetry of the field shape around the center of ionization event in all cases except for $\phi = 0$. The other parameters of the pump pulse as well as the parameters of the medium for the numerical simulations are given after Eq. (10) and Eq. (3) correspondingly.</figcaption></figure>

For FH and SH with equal amplitudes this yields $B = 5/6$. That is, for circularly polarized FH and SH we have

$$\epsilon(\omega) = \frac{5}{6} \frac{\omega}{\omega_0} \quad \text{for} \quad e_{\omega_0} = e_{2\omega_0}, \quad (58)$$

which slightly corrects the result obtained in Ref. [7].
In the hybrid polarization geometry (FH is LP; SH is CP $\epsilon_{\omega_0} = 0, \epsilon_{2\omega_0} = 1$), we find

$$ 
\begin{aligned}
\boldsymbol{\mu}_0 &= -\frac{1}{\omega_0} \sum_{n=1}^{N} \delta\rho_n \frac{e_{\omega_0}}{2} \begin{pmatrix} 0 \\ i \end{pmatrix}, \\
\boldsymbol{\mu}_1 &= -\frac{1}{\omega_0} \sum_{n=1}^{N} \delta\rho_n \begin{pmatrix} e_{\omega_0} + \frac{e_{2\omega_0}}{4} \\ -n\pi e_{2\omega_0} \end{pmatrix}
\end{aligned} 
\quad (59) $$

and

$$ \mathcal{B} = \left| \frac{\mu_{1x}}{\mu_{0y}} \right| = \left| \frac{e_{\omega_0} + e_{2\omega_0}/4}{e_{\omega_0}/2} \right| \quad (60) $$

For FH and SH with equal amplitudes, we easily obtain $\mathcal{B} = 5/2$, that is,

$$ \epsilon(\omega) = \frac{5}{2} \frac{\omega}{\omega_0} \quad \text{for} \quad e_{\omega_0} = e_{2\omega_0}. \quad (61) $$

In Fig. 9, we compare the ellipticities given by Eq. (58) and Eq. (61) with the numerical simulations of the LC model Eq. (5). That is, the spectral intensity of the THz radiation, together with the frequency-dependent ellipticity and rotation angle, is shown for pump pulses with parameters given after Eq. (10), and circular polarization of SH, combined with linear or circular polarization of FH. As observed in Fig. 9(c), the analytical solutions given by Eqs. (58) and (61), which are indicated in green and red dashed lines for the cases of circular and linear polarization of FH, respectively, perfectly match the simulation results represented by blue and orange dots. In contrast to ellipticity, which changes with frequency, the polarization angle for small frequencies is constant. This is supported by Fig. 9(b), where the polarization angle $\Psi(\omega) \equiv -\pi/2$, given by the direction of $A(t_n)$, coincides with the analytical predictions.

## 2. Nonzero relative phase

When $\phi$ is nonzero, Eq. (51) shows that for the CP-S geometry ($\epsilon_{\omega_0} = 1$) the maximum of the electric field shifts from 0 to $-\phi/\omega_0$ modulo $2\pi/\omega_0$. Repeating the previous calculations and applying Eq. (45) leads to the $\phi$-independent expression of the THz ellipticity:

$$ \mathcal{B} = \frac{|\mu_{1y}(\mu_{0x}) - \mu_{1x}(\mu_{0y})|}{|\mu_0|^2} = \left| \frac{e_{\omega_0} + e_{2\omega_0}/4}{e_{\omega_0} + e_{2\omega_0}/2} \right| \quad (62) $$

Thus, the ellipticity does not change with varying $\phi$, which is in agreement with Figs. 8(a)–8(c).

By contrast, for a linearly polarized FH and a circularly polarized SH ($\epsilon_{\omega_0} = 0$), the THz ellipticity becomes dependent on $\phi$. For certain configurations, cf. Eq. (54) and Eq. (55), analytical results for FH and SH with equal amplitudes ($e_{\omega_0} = e_{2\omega_0}$) can be found and Eq. (45) yields

$$ \phi = \frac{\pi}{4} \quad \rightarrow \quad \mathcal{B} = 1 + \frac{3\sqrt{3}}{4}, \quad (63) $$

$$ \phi = \frac{\pi}{2} \quad \rightarrow \quad \mathcal{B} = \frac{\pi}{\sqrt{3}} - \frac{1}{2}. \quad (64) $$

Thus, from the slope $5/2$ achieved when $\phi = 0$, the THz ellipticity decreases for nonzero $\phi$, corresponding to a flattening of the polarization ellipse. The minimum slope is found for

$\phi = \pi/2$, where two ionization events per cycle contribute. This behavior is shown in Fig. 10. We note that the numerical solution of the LC model Eq. (5) yields a significantly lower value for $\mathcal{B}$ for $\phi = \pi/2$ than our analytical model Eq. (45). This discrepancy can be attributed to the finite duration of the ionization events and their asymmetry with respect to the field maxima (see the insets of Fig. 10), which are not taken into account in Eq. (45).

However, Eq. (45) predicts the right trend, namely a decrease of $\mathcal{B}$ when $\phi \rightarrow \pi/2$. This trend can be understood by the appearance of a second ionization event per cycle that renders $|\mu_0|$ and thus the denominator in Eq. (44) large. This decrease in $\mathcal{B}$ for pump configurations with more than one ionization events per optical cycle is, however, not generic because of the complex exponential factors in Eqs. (39) and (40). In particular, Eq. (38) implies that

$$ \mu_0 \propto E_{Br,0}, \quad (65) $$

where $E_{Br,0}$ is the THz field at $\omega \rightarrow 0$. Thus, contributions to $\mu_0$ from different ionization events within the optical cycle may cancel each other out, rendering $|\mu_0|$ very small [see also the discussion about possible vanishing of $E_{Br,0}$ in Appendix B]. In contrast, $\mu_1$ is not affected by this cancellation. Therefore, we observe a significantly larger frequency-averaged ellipticity for cases with small $\mu_0$ in the numerical solutions of the LC model Eq. (5). That is, when the THz energy is small, the ellipticity is large. In fact, this behavior can be observed in Fig. 5.

## D. Frequency-averaged THz ellipticity

The previous analysis suggests a straightforward physical explanation for the nonzero ellipticity of the THz pulses produced by pump fields with finite duration. We can expect that the larger the THz pulse bandwidth $\Delta$, the larger the frequency-averaged ellipticity, provided that $\omega_{co} > \Delta$. The dependence of the ellipticity of THz pulses on the duration of the pump pulse was already reported in Refs. [19,35,68]. Using Eq. (41), we can provide a rough estimate of this dependency. The bandwidth $\Delta$ can be estimated as $\Delta \approx 2/\tau$ (which is a suitable approximation for unchirped Gaussian pulses [36]). We may also assume that the THz spectrum spans from 0 to $\Delta$ and has a roughly constant spectral intensity in this range. Under these assumptions, we obtain

$$ \langle \epsilon(\omega) \rangle \approx \frac{1}{\Delta} \int_{0}^{\Delta} \epsilon(\omega) d\omega = \frac{\mathcal{B}}{\omega_0 \tau}. \quad (66) $$

This estimate is in good agreement with Fig. 7, where the dependency $\langle \epsilon(\omega) \rangle \propto \mathcal{B}/\tau$ predicted by the analytics is clearly visible.

## VI. CONCLUSIONS AND DISCUSSION

In conclusion, we analyzed the THz wave shapes generated via the ionization-based Brunel mechanism by two- and multicolor pump pulses with nontrivial polarization. By numerically solving the local current model and analytical
evaluations, we have shown that photoionization-driven THz wave shapes have in general a frequency-dependent ellipticity. In the limit of zero frequency, the generated radiation is linearly polarized, and the polarization vector points in the direction of the free current formed by the ionized electrons moving away from the parent ions after the pump pulse has passed. The ellipticity is in general linearly dependent on the frequency, i.e., the THz wave shapes possess a “linear ellipticity chirp,” whose magnitude B depends on the particular driving pulse. This behavior sets ionization-based THz emission apart from other schemes, which tend to produce, at least in the first approximation, pulses with constant ellipticity across the frequency bandwidth. This also means, that such pulses can not be characterized anymore by a well defined polarization ellipse. Instead, one must assign such an ellipse to each frequency separately [cf. Fig. 1(b)], or consider frequency-averaged ellipticities and polarization angles, which, however, characterize only partly the polarization state of the pulse.

The ellipticity chirp B tends to increase strongly for those pump wave shapes where free electron currents from several ionization events cancel each other out. Since the net THz energy for such configurations is low, as a rule of thumb, the ellipticity chirp is high for low THz energy and vice versa (for comparable free electron densities). Pulses with nonzero ellipticity chirp B do not have a well-defined polarization state. Yet, one can introduce a frequency-averaged THz ellipticity $\langle\epsilon\rangle$, which is proportional to the spectral width of the pump pulse. Decreasing, for instance, the pump pulse duration $\tau$ increases the THz bandwidth, and thereby increases its average ellipticity, so in general we have $\langle\epsilon\rangle \sim 1/\tau$. In contrast, the polarization angle is in first-order approximation independent of the frequency and coincides with the net free electron current which, in turn, is determined by the values of the vector potential $\mathbf{A}(t_n)$ at the ionization instants $t_n$ ($n = 1, 2, ...$).

Our results demonstrate the unique properties of THz waveforms produced in gas plasmas by polarization controlled multicolor pulses. If the characteristics of the pump (polarization states and relative phases of its components) are similar in the whole plasma volume, these properties should be directly observable in the generated THz radiation. If the pump characteristics change over the pump pulse duration, or in space, such as in the case of filaments, our results can be considered as the first building block for more complex THz wave shapes by superimposing different locally generated waveforms [25,52], giving rise to an even larger diversity of THz wave shapes.

## ACKNOWLEDGMENTS

I.B., A.D. and U.M. are thankful for funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy within the Cluster of Excellence PhoenixD (EXC 2122, Project ID 390833453). V.V. acknowledges support from the “Universities’ Excellence Initiative” programme. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under Grant Agreement No. 871124 Laserlab-Europe. S.S. and L.B. acknowledge

support from the Agence National de la Recherche (ANR-24-CE30-0621).

## APPENDIX A: LOCAL CURRENT (LC) MODEL

In this Appendix, we discuss the derivation of the LC model in detail. The LC model was successfully tested with experiments, either directly [15,16,20,21,24,27,30,36,41,44,46,57–59] or as the basic building block of more advanced models that involve the superposition of LC emitters [25,26,37,60]. It is also often used as a stand-alone theoretical tool [19,22,23,28,29,33,35,38–40,42,43,47] or in the framework of more complex investigations that include, for example, comparisons with simulations describing nonlinear propagation [18,19,22,31–33,38,45,52,53,61,62,69].

We start with the well-known Jefimenko equation [56]

$$ \mathbf{E}_{\text{Br}}(\mathbf{r}, t) = \frac{-1}{4\pi\varepsilon_0 c^2} \int \frac{1}{|\mathbf{r} - \mathbf{r}'|} \frac{\partial \mathbf{J}(\mathbf{r}', t - \frac{|\mathbf{r} - \mathbf{r}'|}{c})}{\partial t} dV', \quad (\text{A1}) $$

which allows us to calculate the electric field $\mathbf{E}_{\text{Br}}$ at the point $\mathbf{r}$ and time $t$, created by the current $\mathbf{J}$ at the point $\mathbf{r}'$ (retarded in time), integrated over all the points belonging to the emitter. In our case, the current $\mathbf{J}$ is produced by the free electrons created through the ionization events, as discussed in the main text. The classical motion of a free electron after its birth is given by

$$ m_e \frac{d\mathbf{v}(t, t')}{dt} = -q\mathbf{E}(t). \quad (\text{A2}) $$

Differentiating Eq. (4) of the main text, taking into account Eq. (A2) and imposing zero electron velocity at birth, we obtain

$$ \partial_t \mathbf{J} = \frac{q^2}{m_e} \rho_e \mathbf{E}, \quad (\text{A3}) $$

see also Ref. [70]. Here we neglected the current decay due to collisions, and our model is valid for nonrelativistic plasmas. We checked (not shown) that collisions only play a minor role in the resulting THz polarization. We note that the simple classical description of the current used in this work may be replaced by semiclassical [25] or fully quantum approaches [7,27,29].

Assuming furthermore a small generation volume $\delta V$, that is, a pointlike source (as seen from the detector) and a homogeneous current across the whole volume, we obtain from Eq. (A1)

$$ \mathbf{E}_{\text{Br}}(\mathbf{r}, t) = g \frac{\partial \mathbf{J}(\mathbf{r}', t - \frac{R}{c})}{\partial t}, \quad g = \frac{-\delta V}{4\pi \varepsilon_0 c^2} \frac{1}{R}, \quad (\text{A4}) $$

where $R = |\mathbf{r} - \mathbf{r}'|$ is the distance from the emitter to the detector. Note that the generated THz radiation is usually very weak, so nonlinear effects in the propagation of the THz radiation on the way to the detector do not need to be considered, while linear dispersion is also neglected. Thus, we assume a linear and dispersionless propagation of THz radiation from the source to the detector. Substituting Eq. (A3) into Eq. (A4) and ignoring the time delay that is not relevant to our purposes, we obtain Eq. (5).
APPENDIX B: FORMULATION IN TERMS OF VECTOR POTENTIAL AND PHYSICAL MEANING OF E<sub>Br</sub> IN THE LIMIT ω → 0

This Appendix is dedicated to a deeper analysis of the physical meaning of Eq. (25). With the Fourier series expansion of the vector potential

$$A(t) = \sum_{m} A_{m} e^{-im\omega_{0}t}, \quad (B1)$$

we find

$$A_m = -\frac{i}{m\omega_0} E_m, \quad m \neq 0 \quad (B2)$$

for the mth harmonic. If we further set A<sub>0</sub> = 0, we get

$$A(t) = -\int_{-\infty}^{t} E(t')dt', \quad (B3)$$

so that Eq. (25) can be expressed as

$$E_{Br,0} \propto \sum_{m} \sum_{n=1}^{N_c} \delta \rho_n A_m e^{-im\omega_0 t_n}. \quad (B4)$$

Equation (B4) as well as Eq. (25) provide a “mixed representation”: They contain the harmonics of the pump, E<sub>m</sub> or A<sub>m</sub>, respectively, and, at the same time, the step sizes of ionization δρ<sub>n</sub> in time.

The physics behind Eq. (26) is easy to understand: The electron, leaving the atom with negligible initial velocity, experiences acceleration in the field [cf. Eq. (A2)]. Integrating this expression from the electron birth time t<sub>n</sub> to t = +∞, we obtain the electron velocity after the interaction with the

electric field,

$$v(t \to +\infty, t_n) \propto \int_{t_n}^{+\infty} E(t')dt'. \quad (B5)$$

Because

$$0 = E_0 = \int_{-\infty}^{+\infty} E(t')dt' = \int_{-\infty}^{t_n} E(t')dt' + \int_{t_n}^{+\infty} E(t')dt',$$

we find with Eq. (B3) that

$$v(t \to +\infty, t_n) \propto A(t_n). \quad (B6)$$

The fact that the electron velocity is proportional to the vector potential at the time of ionization is well known [71]. It is the heart of the so-called atto-streaking technique [72,73], which allows one to measure the waveform of a strong pulse in a rather direct way – in contrast to other pulse-characterization techniques.

The net current density produced over one optical cycle reads

$$\Delta \mathbf{J} = \sum_{n=1}^{N_c} \delta \rho_n v(t \to +\infty, t_n) \propto \sum_{n=1}^{N_c} \delta \rho_n A(t_n), \quad (B7)$$

which was already exploited in the context of THz generation [19,38,61,74].<fnref n="2" /> Since we assume that there is no net current density before the pulse arrives, J(t → -∞) = 0, we find

$$E_{Br,0} \propto \Delta \mathbf{J} = \mathbf{J}(t \to +\infty) - \mathbf{J}(t \to -\infty) = \mathbf{J}(t \to +\infty).$$

Thus, the 0th Brunel harmonic is linearly polarized in the direction in which the electrons fly after the pulse has passed. This property is valid for arbitrary pump shapes.

[1] W. L. Chan, J. Deibel, and D. M. Mittleman, Imaging with terahertz radiation, Rep. Prog. Phys. 70, 1325 (2007).

[2] M. Tonouchi, Cutting-edge terahertz technology, Nat. Photon. 1, 97 (2007).

[3] S. S. Dhillon, M. S. Vitiello, E. H. Linfield, A. G. Davies, M. C. Hoffmann, J. Booske, C. Paoloni, M. Gensch, P. Weightman, G. P. Williams et al., The 2017 terahertz science and technology roadmap, J. Phys. D 50, 043001 (2017).

[4] Zh. Ma, P. Li, S. Chen, and X. Wu, Optical generation of strong-field terahertz radiation and its application in nonlinear terahertz metasurfaces, Nanophotonics 11, 1847 (2022).

[5] T. L. Cocker, V. Jelic, M. Gupta, S. J. Molesky, J. A. J. Burgess, G. De Los Reyes, L. V. Titova, Y. Y. Tsui, M. R. Freeman, and F. A. Hegmann, An ultrafast terahertz scanning tunnelling microscope, Nat. Photon. 7, 620 (2013).

[6] S. Katletz, M. Pfleger, H. Pühringer, M. Mikulics, N. Vieweg, O. Peters, B. Scherger, M. Scheller, M. Koch, and K. Wiesauer, Polarization sensitive terahertz imaging: Detection of birefringence and optical axis, Opt. Express 20, 23025 (2012).

[7] I. Babushkin, A. J. Galán, J. R. C. de Andrade, A. Husakou, F. Morales, M. Kretschmar, T. Nagy, V. Vaičaitis, L. Shi, D. Zuber et al., All-optical attoclock for imaging tunnelling wavepackets, Nat. Phys. 18, 417 (2022).

[8] V. Vaičaitis, O. Balachninaitė, A. Matijošius, I. Babushkin, and U. Morgner, Direct time-resolved plasma characterization with broadband terahertz light pulses, Phys. Rev. E 107, 015201 (2023).

[9] J. Liu, J. Dai, S. L. Chin, and X.-C. Zhang, Broadband terahertz wave remote sensing using coherent manipulation of fluorescence from asymmetrically ionized gases, Nat. Photon. 4, 627 (2010).

[10] T.-J. Wang, Sh. Yuan, Y. Chen, J.-F. Daigle, C. Marceau, F. Théberge, M. Châteauneuf, J. Dubois, and S. L. Chin, Toward remote high energy terahertz generation, Appl. Phys. Lett. 97, 111108 (2010).

[11] E. R. Brown, Fundamentals of terrestrial millimeter-wave and THz remote sensing, Int. J. Hi. Spe. Ele. Syst. 13, 995 (2003).

[12] E. Pickwell and V. P. Wallace, Biomedical applications of terahertz technology, J. Phys. D 39, R301 (2006).

<footnote n="2">²The relevant quantity in Ref. [19] is v<sub>f</sub>(t) = ∫<sub>-∞</sub><sup>t</sup> E(t')e<sup>γ(t-t')</sup>dt', where γ describes collision-induced decay of the current. Neglecting γ restores, up to a minus sign, the vector potential A(t) on the right side of this expression.</footnote>
[13] M. Kress, T. Löffler, S. Eden, M. Thomson, and H. G. Roskos, Terahertz-pulse generation by photoionization of air with laser pulses composed of both fundamental and second-harmonic waves, Opt. Lett. 29, 1120 (2004).

[14] T. Bartel, P. Gaal, K. Reimann, M. Woerner, and T. Elsaesser, Generation of single-cycle THz transients with high electric-field amplitudes, Opt. Lett. 30, 2805 (2005).

[15] K.-Y. Kim, J. H. Glownia, A. J. Taylor, and G. Rodriguez, Terahertz emission from ultrafast ionizing air in symmetry-broken laser fields, Opt. Express 15, 4577 (2007).

[16] K. Y. Kim, A. G. Taylor, A. G. Glownia, and G. Rodriguez, Coherent control of terahertz supercontinuum generation in ultrafast laser-gas interactions, Nat. Photon. 2, 605 (2008).

[17] A. D. Koulouklidis, C. Gollner, V. Shumakova, V. Yu. Fedorov, A. Pugžlys, A. Baltuška, and S. Tzortzakis, Observation of extremely efficient terahertz generation from mid-infrared two-color laser filaments, Nat. Commun. 11, 292 (2020).

[18] H. Alirezaee and M. Sharifian, Contribution of photocurrent mechanism and influence of plasma length in THz generation by two-color laser induced plasma, Phys. Plasmas 25, 043112 (2018).

[19] I. Babushkin, S. Skupin, A. Husakou, C. Köhler, E. Cabrera-Granado, L. Bergé, and J. Herrmann, Tailoring terahertz radiation by controlling tunnel photoionization events in gases, New J. Phys. 13, 123029 (2011).

[20] M. Clerici, M. Peccianti, B. E. Schmidt, L. Caspani, M. Shalaby, M. Giguère, A. Lotti, A Couairon, F. Légaré, T. Ozaki et al., Wavelength scaling of terahertz generation by gas ionization, Phys. Rev. Lett. 110, 253901 (2013).

[21] N. V. Vvedenskii, A. I. Korytin, V. A. Kostin, A. A. Murzanev, A. A. Silaev, and A. N. Stepanov, Two-color laser-plasma generation of terahertz radiation using a frequency-tunable half harmonic of a femtosecond pulse, Phys. Rev. Lett. 112, 055004 (2014).

[22] A. Nguyen, P. González de Alaiza Martínez, J. Déchard, I. Thiele, I. Babushkin, S. Skupin, and L. Bergé, Spectral dynamics of THz pulses generated by two-color laser filaments in air: the role of Kerr nonlinearities and pump wavelength, Opt. Express 25, 4720 (2017).

[23] I. Babushkin, C. Brée, Ch. M. Dietrich, A. Demircan, U. Morgner, and A. Husakou, Terahertz and higher-order Brunel harmonics: From tunnel to multiphoton ionization regime in tailored fields, J. Mod. Opt. 64, 1078 (2017).

[24] Ch. Meng, W. Chen, X. Wang, Zh. Lü, Y. Huang, J. Liu, D. Zhang, Z. Zhao, and J. Yuan, Enhancement of terahertz radiation by using circularly polarized two-color laser fields, Appl. Phys. Lett. 109, 131105 (2016).

[25] Zh. Zhang, Y. Chen, S. Cui, F. He, M. Chen, Z. Zhang, J. Yu, L. Chen, Zh. Sheng, and J. Zhang, Manipulation of polarizations for broadband terahertz waves emitted from laser plasma filaments, Nat. Photon. 12, 554 (2018).

[26] Y. S. You, T. I. Oh, and K.-Y. Kim, Mechanism of elliptically polarized terahertz generation in two-color laser filamentation, Opt. Lett. 38, 1034 (2013).

[27] J. Dai, N. Karpowicz, and X.-C. Zhang, Coherent polarization control of terahertz waves generated from two-color laser-induced gas plasma, Phys. Rev. Lett. 103, 023001 (2009).

[28] V. Yu. Fedorov, A. D. Koulouklidis, and S. Tzortzakis, THz generation by two-color femtosecond filaments with complex

polarization states: Four-wave mixing versus photocurrent contributions, Plasma Phys. Control. Fusion 59, 014025 (2017).

[29] V. A. Tulsky, M. Baghery, U. Saalmann, and S. V. Popruzhenko, Boosting terahertz-radiation power with two-color circularly polarized midinfrared laser pulses, Phys. Rev. A 98, 053415 (2018).

[30] H. Wen and A. M. Lindenberg, Coherent terahertz polarization control through manipulation of electron trajectories, Phys. Rev. Lett. 103, 023902 (2009).

[31] O. Kosareva, M. Esaulkov, N. Panov, V. Andreeva, D. Shipilo, P. Solyankin, A. Demircan, I. Babushkin, V. Makarov, U. Morgner et al., Polarization control of terahertz radiation from two-color femtosecond gas breakdown plasma, Opt. Lett. 43, 90 (2018).

[32] C. Tailliez, A. Stathopulos, S. Skupin, D. Buožius, I. Babushkin, V. Vaičaitis, and L. Bergé, Terahertz pulse generation by two-color laser fields with circular polarization, New J. Phys. 22, 103038 (2020).

[33] A. Stathopulos, S. Skupin, and L. Bergé, Terahertz pulse generation by multi-color laser fields with linear versus circular polarization, Opt. Lett. 46, 5906 (2021).

[34] Q. Song, Q. Lin, H. Wang, H. Zhong, Y. Cai, Sh. Zheng, Zh. Chen, X. Lu, X. Zeng, H. Shangguan et al., Efficient nearly-circularly-polarized terahertz generation from an air plasma pumped by collinear and circularly polarized two-color laser fields, Phys. Rev. A 102, 023506 (2020).

[35] Ch. Lu, T. He, L. Zhang, H. Zhang, Y. Yao, Sh. Li, and Sh. Zhang, Effect of two-color laser pulse duration on intense terahertz generation at different laser intensities, Phys. Rev. A 92, 063850 (2015).

[36] A. V. Borodin, N. A. Panov, O. G. Kosareva, V. A. Andreeva, M. N. Esaulkov, V. A. Makarov, A. P. Shkurinov, S. L. Chin, and X.-C. Zhang, Transformation of terahertz spectra emitted from dual-frequency femtosecond pulse interaction in gases, Opt. Lett. 38, 1906 (2013).

[37] Y. Bai, L. Song, R. Xu, Ch. Li, P. Liu, Zh. Zeng, Z. Zhang, H. Lu, R. Li, and Zh. Xu, Waveform-controlled terahertz radiation from the air filament produced by few-cycle laser pulses, Phys. Rev. Lett. 108, 255004 (2012).

[38] P. González de Alaiza Martínez, I. Babushkin, L. Bergé, S. Skupin, E. Cabrera-Granado, C. Köhler, U. Morgner, A. Husakou, and J. Herrmann, Boosting terahertz generation in laser-field ionized gases using a sawtooth wave shape, Phys. Rev. Lett. 114, 183901 (2015).

[39] L. Zhang, G.-L. Wang, S.-F. Zhao, and X.-X. Zhou, Controlling of strong tunable THz emission with optimal incommensurate multi-color laser field, Phys. Plasmas 24, 023116 (2017).

[40] H. Wang, W. Fan, X. Chen, and H. Yan, Polarization control of terahertz waves generated by a femtosecond three-color pulse scheme, J. Opt. Soc. Am. B 40, 1375 (2023).

[41] T. Balčiūnas, D. Lorenc, M. Ivanov, O. Smirnova, A. M. Zheltikov, D. Dietze, K. Unterrainer, T. Rathje, G. G. Paulus, A. Baltuška et al., CEP-stable tunable THz-emission originating from laser-waveform-controlled sub-cycle plasma-electron bursts, Opt. Express 23, 15278 (2015).

[42] L. Zhang, G.-L. Wang, and X.-X. Zhou, Optimized two-and three-colour laser pulses for the intense terahertz wave generation, J. Mod. Opt. 63, 2159 (2016).

[43] Ch. Lu, Ch. Zhang, L. Zhang, X. Wang, and Sh. Zhang, Modulation of terahertz-spectrum generation from an air plasma
by tunable three-color laser pulses, Phys. Rev. A 96, 053402 (2017).

[44] V Vaičaitis, O Balachninaitė, U Morgner, and I Babushkin, Terahertz radiation generation by three-color laser pulses in air filament, J. Appl. Phys. 125, 173103 (2019).

[45] H. Alirezaee, M. Sharifian, S. Mohammad R. Darbani, M. Saeed, A. Eslami Majd, and A. R. Niknam, Terahertz radiation emission from three-color laser-induced air plasma, Eur. Phys. J. Plus 135, 342 (2020).

[46] Sh. Liu, Zh. Fan, Ch. Lu, J. Gui, Ch. Luo, Sh. Wang, Q. Liang, B. Zhou, A. Houard, A. Mysyrowicz, V. Kostin, and Y. Liu, Coherent control of boosted terahertz radiation from air plasma pumped by a femtosecond three-color sawtooth field, Phys. Rev. A 102, 063522 (2020).

[47] D. Ma, L. Dong, R. Zhang, C. Zhang, Y. Zhao, and L. Zhang, Enhancement of terahertz wave emission from air plasma excited by harmonic three-color laser fields, Opt. Commun. 481, 126533 (2021).

[48] S. Watanabe, Terahertz polarization imaging and its applications, Photonics 5, 58 (2018).

[49] S. Baierl, M. Hohenleutner, T. Kampfrath, A. K. Zvezdin, A. V. Kimel, R. Huber, and R. V. Mikhaylovskiy, Nonlinear spin control by terahertz-driven anisotropy fields, Nat. Photon. 10, 715 (2016).

[50] H. Hoshina, Y. Morisawa, H. Sato, H. Minamide, I. Noda, Y. Ozaki, and C. Otani, Polarization and temperature dependent spectra of poly (3-hydroxyalkanoate)s measured at terahertz frequencies, Phys. Chem. Chem. Phys. 13, 9173 (2011).

[51] N. Amer, W. C. Hurlbut, B. J. Norton, Y.-Sh. Lee, and T. B. Norris, Generation of terahertz pulses with arbitrary elliptical polarization, Appl. Phys. Lett. 87, 221111 (2005).

[52] C. Köhler, C. Cabrera-Granado, I. Babushkin, L. Bergé, J. Herrmann, and S. Skupin, Directionality of terahertz emission from photoinduced gas plasmas, Opt. Lett. 36, 3166 (2011).

[53] E. Cabrera-Granado, Y. Chen, I. Babushkin, L. Bergé, and S. Skupin, Spectral self-action of THz emission from ionizing two-color laser pulses in gases, New J. Phys. 17, 023060 (2015).

[54] F. Brunel, Harmonic generation due to plasma effects in a gas undergoing multiphoton ionization in the high-intensity limit, J. Opt. Soc. Am. B 7, 521 (1990).

[55] L. D. Landau and L. M. Lifshitz, *Quantum Mechanics non-Relativistic Theory: Vol. 3*, 3rd ed. (Butterworth-Heinemann, London, 1981).

[56] O. D. Jefimenko, *Electricity and Magnetism: An Introduction to the Theory of Electric and Magnetic Fields* (Electret Scientific Company, Star City, WV, 1989).

[57] T. I. Oh, Y. S. You, and K. Y. Kim, Two-dimensional plasma current and optimized terahertz generation in two-color photoionization, Opt. Express 20, 19778 (2012).

[58] A. Martinez, J. Houard, A. Hideur, D. Paparo, and A. Vella, Controlling the time shape of terahertz pulses from two-color plasma by combining wavelength dispersion and laser chirp, Appl. Phys. Lett. 124, 021105 (2024).

[59] D. Paparo, A. Martinez, A. Rubano, J. Houard, A. Hideur, and A. Vella, THz generation by two-color plasma: Time shaping and ultra-broadband polarimetry, Sensors 24, 4265 (2024).

[60] A. Yousef-Zamanian and M. Neshat, Investigation of polarization state of terahertz radiation from compact laser-induced plasma in air, J. Mod. Opt. 64, 300 (2017).

[61] A. Debayle, L. Gremillet, L. Bergé, and Ch. Köhler, Analytical model for THz emissions induced by laser-gas interaction, Opt. Express 22, 13691 (2014).

[62] A. Gorodetsky, A. D. Koulouklidis, M. Massaouti, and S. Tzortzakis, Physics of the conical broadband terahertz emission from two-color laser-induced plasma filaments, Phys. Rev. A 89, 033838 (2014).

[63] I. Thiele, B. Zhou, A. Nguyen, E. Smetaninaa, R. Nuter, K. J. Kaltenecker, P. González de Alaiza Martínez, J. Déchard, L. Bergé, P. U. Jepsen, and S. Skupin, Terahertz emission from laser-driven gas plasmas: A plasmonic point of view, Optica 5, 1617 (2018).

[64] M. Rasmussen, O. Nagy, S. Skupin, A. Stathopulos, L. Bergé, P. Uhd Jepsen, and B. Zhou, Frequency-resolved characterization of broadband two-color air-plasma terahertz beam profiles, Opt. Express 31, 9287 (2023).

[65] A. Stathopulos, S. Skupin, B. Zhou, P. U. Jepsen, and L. Bergé, Waveshape of terahertz radiation produced by two-color laser-induced air plasmas, Phys. Rev. Res. 6, 043274 (2024).

[66] E. Hecht, *Optics*, 5th ed. (Pearson Education, Upper Saddle River, NJ, 2017).

[67] L. Bergé, C.-L. Soulez, C. Köhler, and S. Skupin, Role of the carrier-envelope phase in laser filamentation, Appl. Phys. B 103, 563 (2011).

[68] T.-J. Wang, Y. Chen, C. Marceau, F. Théberge, M. Châteauneuf, J. Dubois, and S. L. Chin, High energy terahertz emission from two-color laser-induced filamentation in air with pump pulse duration control, Appl. Phys. Lett. 95, 131108 (2009).

[69] I. Babushkin, W. Kuehn, C. Köhler, S. Skupin, L. Bergé, K. Reimann, M. Woerner, J. Herrmann, and T. Elsaesser, Ultrafast spatiotemporal dynamics of terahertz generation by ionizing two-color femtosecond pulses in gases, Phys. Rev. Lett. 105, 053903 (2010).

[70] I. Babushkin, S. Skupin, and J. Herrmann, Generation of terahertz radiation from ionizing two-color laser pulses in Ar filled metallic hollow waveguides, Opt. Express 18, 9658 (2010).

[71] W. B. Mori and T. Katsouleas, Ponderomotive force of a uniform electromagnetic wave in a time varying dielectric medium, Phys. Rev. Lett. 69, 3495 (1992).

[72] E. Goulielmakis, M. Uiberacker, R. Kienberger, A. Baltuska, V. Yakovlev, A. Scrinzi, Th. Westerwalbesloh, U. Kleineberg, U. Heinzmann, M. Drescher, and F. Krausz, Direct measurement of light waves, Science 305, 1267 (2004).

[73] E. Goulielmakis, M. Schultze, M. Hofstetter, V. S. Yakovlev, J. Gagnon, M. Uiberacker, A. L. Aquila, E. M. Gullikson, D. T. Attwood, R. Kienberger et al., Single-cycle nonlinear optics, Science 320, 1614 (2008).

[74] A. Debayle, P. González de Alaiza Martínez, L. Gremillet, and L. Bergé, Nonmonotonic increase in laser-driven THz emissions through multiple ionization events, Phys. Rev. A 91, 041801 (2015).