# Universal properties of locally generated terahertz waveforms from
polarization-controlled two- and multi-color ionizing fields

H. Alirezaee  Center for Research on Laser and Plasma, Shahid Chamran
University of Ahvaz, Iran  S. Skupin  Institut Lumière Matière, UMR 5306 -
CNRS, Université de Lyon 1, 69622 Villeurbanne, France  V. Vaicaitis  Laser
Research Center, Vilnius University, Saulėtekio 10, Vilnius LT-10223,
Lithuania  A. Demircan  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  I. Babushkin  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  Luc
Bergé  Centre Lasers Intenses et Applications, Université de
Bordeaux–CNRS–CEA, 33405 Talence Cedex, France  U. Morgner  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

###### Abstract

The polarization states of terahertz (THz) radiation generated in a photo-
ionized gas driven by strong 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 the frequency. This “linear chirp of ellipticity” makes plasma-
based THz generation unique among other THz sources. However, it also puts
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.

##  I Introduction

Terahertz (THz) radiation has various applications, such as characterization,
spectroscopy and imaging of materials, nanostructures and plasmas Chan  _et
al._ (2007); Tonouchi (2007); Marx (2007); Dhillon  _et al._ (2017); Ma  _et
al._ (2022); Cocker  _et al._ (2013); Katletz  _et al._ (2012); Babushkin  _et
al._ (2022); Vaičaitis  _et al._ (2023), remote sensing Liu  _et al._ (2010);
Wang  _et al._ (2010); Brown (2003), biomedical applications Pickwell and
Wallace (2006); Marx (2007); Chan  _et al._ (2007); Tonouchi (2007); Dhillon
_et al._ (2017), 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.
Kress  _et al._ (2004); Bartel  _et al._ (2005); Kim  _et al._ (2007, 2008);
Koulouklidis  _et al._ (2020). 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 produced by two-color laser plasmas. In
particular, the efficiency of THz generation was shown to be influenced by the
relative phases between laser harmonics Kim  _et al._ (2007); Alirezaee and
Sharifian (2018); Babushkin  _et al._ (2011). Other important factors
affecting the efficiency of the produced THz radiation are pump wavelengths
Clerici  _et al._ (2013); Vvedenskii  _et al._ (2014); Nguyen  _et al._
(2017); Babushkin  _et al._ (2017), polarization of the pulses Meng  _et al._
(2016a); Zhang  _et al._ (2018); You  _et al._ (2013); Dai  _et al._ (2009);
Fedorov  _et al._ (2016); Tulsky  _et al._ (2018); Wen and Lindenberg (2009);
Kosareva  _et al._ (2018); Tailliez  _et al._ (2020); Stathopulos  _et al._
(2021); Meng  _et al._ (2016b); Song  _et al._ (2020), their durations
Babushkin  _et al._ (2011); Lu  _et al._ (2015); Borodin  _et al._ (2013) and,
more generally, their waveshapes Babushkin  _et al._ (2011); Bai  _et al._
(2012); de Alaiza Martínez  _et al._ (2015); Zhang  _et al._ (2017, 2018);
Wang  _et al._ (2023). In addition to the conventional two-color pump
configuration, multi-color pump fields have attracted considerable interest in
recent years Babushkin  _et al._ (2011); de Alaiza Martínez  _et al._ (2015);
Balčiūnas  _et al._ (2015); Zhang  _et al._ (2016, 2017); Lu  _et al._ (2017);
Vaičaitis  _et al._ (2019); Alirezaee  _et al._ (2020); Liu  _et al._ (2020);
Stathopulos  _et al._ (2021); Ma  _et al._ (2021); Wang  _et al._ (2023).

Many THz applications require polarization control Watanabe (2018); Cocker
_et al._ (2013); Baierl  _et al._ (2016); Katletz  _et al._ (2012); Hoshina
_et al._ (2011). The most common techniques for controlling THz polarization
are based on wire-grid polarizers Watanabe (2018); Amer  _et al._ (2005).
However, even if they might work in the whole THz region, they accept only
narrowband THz pulses. So far, most studies on THz radiation in two-color
schemes have been done by considering linear polarization of the pump. Apart
from specific dedicated works Meng  _et al._ (2016a); Zhang  _et al._ (2018);
You  _et al._ (2013); Dai  _et al._ (2009); Fedorov  _et al._ (2016); Tulsky
_et al._ (2018); Wen and Lindenberg (2009); Kosareva  _et al._ (2018);
Tailliez  _et al._ (2020); Stathopulos  _et al._ (2021); Meng  _et al._
(2016b); Song  _et al._ (2020) 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
superimposed, resulting in complex THz waveshapes Zhang  _et al._ (2018); You
_et al._ (2013); Köhler  _et al._ (2011); Cabrera-Granado  _et al._ (2015). In
these propagation scenarios, the THz field generated downstream is also
influenced by the previously generated THz components Cabrera-Granado  _et
al._ (2015).

In this paper, we consider the most important building block of ionization-
based THz generation – the locally generated THz waveforms in small spatial
volumes. We address a rich variety of possible THz waveforms, which can appear
in two- or multi-color schemes with arbitrary polarization of the pump pulses,
and particularly focus on the polarization state of such THz radiation.
Starting from general principles, we show how the 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 waveshape on the pump
configuration for a two-color scheme with controlled polarization and phase
difference between the two pump components.

##  II System and Model

###  II.1 Vectorial formulation of local current model

Figure 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. In this article, we assume that the polarization
ellipses and relative phases of harmonics remain constant throughout the pump
pulse. (b) Schematic spectrum showing polarization angle, ellipticity, and
intensity of the resulting THz radiation, which exhibits an intrinstic
“ellipticity chirp”.

The scheme of our setup is illustrated in Fig. 1a. 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 Brunel (1990);
Balčiūnas  _et al._ (2015); Babushkin  _et al._ (2017, 2022) (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 You  _et al._ (2013); Zhang  _et al._
(2018); Kosareva  _et al._ (2018); Köhler  _et al._ (2011); Cabrera-Granado
_et al._ (2015), here we focus on the response of an infinitely small emitting
volume of plasma. This corresponds to analyzing the elementary source term one
would have to introduce into, e.g., Maxwell equations to address later a more
complex propagation geometry.

The electric field of the pump pulse in the multi-color case, consisting of
the fundamental harmonic at frequency $\omega_{0}$ and their harmonics
$m\omega_{0}$, $m=2,3,\ldots$ is given by:

$\begin{split}\mathrm{\mathbf{E}}(t)=\sum_{m\geq
1}f_{m\omega_{0}}(t)e_{m\omega_{0}}\Big{[}&{\cos{\left(m\omega_{0}t+\phi_{x,m\omega_{0}}\right)}\mathrm{\mathbf{x}}}\\\
+{\epsilon}_{m\omega_{0}}&{\sin{\left(m\omega_{0}t+\phi_{y,m\omega_{0}}\right)}\mathrm{\mathbf{y}}}\Big{]}.\end{split}$
(1)

Here we assume laser field components polarized in $x$ and $y$ directions,
$\mathrm{\mathbf{x}}$ and $\mathrm{\mathbf{y}}$ being the corresponding unit
vectors, $e_{m\omega_{0}}$ define the amplitudes of the $m$th harmonics,
$f_{m\omega_{0}}(t)$ define the slow envelopes of the pulses,
${\epsilon}_{m\omega_{0}}$ determine their 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}\left(t\right)=W\left(t\right)\left[{\rho}_{\rm
at}-{\rho}_{e}\left(t\right)\right],$ (2)

with ${\rho}_{\rm at}$ being the neutral density. For simplicity, the
ionization rate $W\left(t\right)$ can be evaluated from the quasi-static
tunneling model Landau and Lifshitz (1981),

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

with constants $\alpha$ and $\beta$ characterizing the ionization potential.
In the numerical simulations below we assume argon (ionization energy of 15.6
eV).

The electrons released acquire a drift velocity
$\mathrm{\mathbf{v}}(t,t^{\prime})=(q/m_{e})\int_{t^{\prime}}^{t}\vec{E}(t^{\prime\prime})dt^{\prime\prime}$,
where $t^{\prime}$ denotes the time when they were born and $q$ denotes the
electron charge. For simplicity, we assume that electrons are born with zero
drift velocity, and collisions are neglected in our analytical approach. The
macroscopic current density ${\mathrm{\mathbf{J}}}$ is then given by

${\mathrm{\mathbf{J}}}(t)=q\int^{t}_{-\infty}\mathrm{\mathbf{v}}(t,t^{\prime})\;\left[\partial_{t^{\prime}}\rho_{e}(t^{\prime})\right]dt^{\prime}.$
(4)

This net current density is responsible for the emission of the so-called
Brunel radiation ${\mathrm{\mathbf{E}}}_{\mathrm{Br}}$, which can be estimated
as

${{\mathrm{\mathbf{E}}}}_{\mathrm{Br}}(t)\propto{\partial}_{t}\mathrm{\mathbf{J}}=\frac{q^{2}}{m_{e}}{\rho_{e}{\mathrm{\mathbf{E}}}}.$
(5)

This key equation is often referred to as the local current (LC)
approximation. Note that this 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. Propagation effects in more complex
geometries can then be explained by superimposing these local contributions
You  _et al._ (2013); Zhang  _et al._ (2018); Kosareva  _et al._ (2018);
Köhler  _et al._ (2011); Cabrera-Granado  _et al._ (2015), provided that
charge separation effects in the plasma can be neglected Thiele  _et al._
(2018). We mention that the LC model does not take into account any effects
related to the electron deflection by the tail of the atomic Coulomb potential
Babushkin  _et al._ (2022).

###  II.2 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. 1b). To
describe this property, instead of considering the pump field
$\mathrm{\mathbf{E}}(t)$ and the Brunel radiation
$\mathrm{\mathbf{E}}_{\mathrm{Br}}$ in time, we operate with their Fourier
components
$\hat{\mathrm{\mathbf{E}}}(\omega)\propto\int\mathrm{\mathbf{E}}(t)\exp{-i\omega
t)}\mathrm{d}t$, and analogous to this
$\hat{\mathrm{\mathbf{E}}}_{\mathrm{Br}}(\omega)$.

Furthermore, we shall use the Jones formalism to describe our THz waveforms.
To any (real-valued) monochromatic field $\mathrm{\mathbf{E}}(t)$ oscillating
at frequency $\omega$, we assign a complex vector $\mathrm{\mathbf{j}}$,
describing its polarization state. For example, for the field
$\mathrm{\mathbf{E}}(t)=(\cos{\omega t},\epsilon\sin{\omega t})$ we consider
its complex representation $\propto(1e^{-i\omega t},i\epsilon e^{-i\omega
t})$, and $\mathrm{\mathbf{j}}$ can be defined as the pre-factors before the
oscillating term, that is, for this example,
$\mathrm{\mathbf{j}}=(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\leq\alpha$. The sign of $\epsilon$ specifies the direction of rotation
of the electric field. This direction is given by assigning a positive
(clockwise) or negative (anti-clockwise) sign to $\epsilon$. In the example
above, we chose anti-clockwise polarization. In this way, we can assign a
polarization state to $\hat{\mathrm{\mathbf{E}}}(\omega)$ and
$\hat{\mathrm{\mathbf{E}}}_{\mathrm{Br}}(\omega)$ for each frequency $\omega$,
thus obtaining a frequency-dependent ellipticity $\epsilon(\omega)$ and
polarization angle $\Psi(\omega)$.

Figure 2: Polarization ellipse in the laboratory coordinate system ($x,y$) and
the principal axes coordinate system ($\xi,\eta$). The angle between $x$-axis
of the laboratory corrdinate system and the major axis of the ellipse is
denoted by $\Psi$. The ellipticity of the polarization state $\epsilon$ is
defined as the ratio between the lengths of the minor and major axes. The sign
of $\epsilon$ specifies the direction of rotation of the electric field (see
text for details).

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\limits_{0}^{\omega_{\mathrm{co}}}\epsilon(\omega)I(\omega)\mathrm{d}\omega}{\int\limits_{0}^{\omega_{\mathrm{co}}}I(\omega)\mathrm{d}\omega},$
(6)

where $\epsilon(\omega)$ denotes the frequency-dependent ellipticity and
$I(\omega)\propto|\hat{\mathrm{\mathbf{E}}}_{\mathrm{Br}}(\omega)|^{2}$ is the
THz spectral intensity, later extracted from the LC model. The cutoff
frequency $\omega_{\mathrm{co}}$ (in the upcoming simulations we took
$\omega_{\mathrm{co}}=\omega_{0}/4$) selects the frequency range of interest
and thus gives a physical meaning to the above quantity, assuming that
$\epsilon(\omega)$ does not change sign. 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-
nonresolving THz detection.

##  III Two-color case: Exemplary waveshapes

Figure 3: A typical scheme of THz generation using a two-color pump,
consisting of funamental harmonic (FH) and second harmonic (SH) pump pulses.
The ellipticities of the initially linearly polarized FH and SH pulses are
tunable by rotating the quarter wave plates (QWP) by the angles $\theta_{1}$
and $\theta_{2}$, respectively. The angle $\phi$ denotes the phase difference
between FH and SH.

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, Eq. (1) reduces to

$\displaystyle{\mathrm{\mathbf{E}}}$
$\displaystyle={\mathrm{\mathbf{E}}}_{\omega_{0}}+{\mathrm{\mathbf{E}}}_{2\omega_{0}},$
(7)
$\displaystyle{{\mathrm{\mathbf{E}}}}_{\omega_{0}}=e_{\omega_{0}}e^{{-t}^{2}/{\tau}^{2}}\Big{[}$
$\displaystyle{\cos{\left(\omega_{0}t\right)}\mathrm{\mathbf{x}}}+{\epsilon}_{\omega_{0}}{\sin{\left(\omega_{0}t\right)}\mathrm{\mathbf{y}}}\Big{]},$
(8)
$\displaystyle{{\mathrm{\mathbf{E}}}}_{2\omega_{0}}=e_{2\omega_{0}}e^{{-t}^{2}/{\tau}^{2}}\Big{[}$
$\displaystyle{\cos{\left(2\omega_{0}t+\phi\right)}\mathrm{\mathbf{x}}}$
$\displaystyle+{\epsilon}_{2\omega_{0}}{\sin{\left(2\omega_{0}t+\phi\right)}\mathrm{\mathbf{y}}}\Big{]}.$
(9)

Here ${\mathrm{\mathbf{E}}}_{\omega_{0}}$ and
${\mathrm{\mathbf{E}}}_{2\omega_{0}}$ are electric fields of fundamental
harmonic (FH) and second harmonic (SH) pulses, respectively, with $\phi$ being
their relative phase. A typical setup is shown in Fig. 3. Before the pump is
focused, the ellipticities of FH and SH can be controlled by rotating the
quarter wave plates (QWP) by angles $\theta_{1}$ and $\theta_{2}$,
respectively. This affects the $E_{x}$ and $E_{y}$ components of each pump
harmonic by $\cos{\theta_{j}}$ and $\sin{\theta_{j}}$ (i.e.,
$\epsilon_{j\omega_{0}}=\tan{\theta_{j}}$, $j=1,2$). 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
waveshapes. Our baseline pump configurations when $\theta_{1}$ is fixed to 0
($x$-linearly polarized FH pulse) or $\pi/4$ (circularly polarized FH pulse)
are displayed in Fig. 4 with blue lines for zero relative phase $\phi$. For
the sake of clarity, in all cases the same field amplitudes for FH and SH are
taken, so that their extrema coincide. These Lissajou figures describe the
electric field vector in time for different values of $\theta_{2}$ over one
oscillation period and for equal amplitudes of FH and SH.

Figure 4: Pump waveshapes over one optical cycle for different configurations
of QWPs, indicated by numbers in Fig. 5 for (a) circular polarization with
clockwise rotation ($\theta_{1}=\pi/4$) and (b) linear polarization in $x$
direction ($\theta_{1}=0$) of the FH, assuming $\phi=0$. FH and SH components
have equal amplitudes. The red dots represent instants $t_{n}$ of the maxima
of $|\mathrm{\mathbf{E}}(t)|$ (irrelevant $t_{n}$ are omitted), and the arrows
indicate the direction of the product
$\mathrm{\mathbf{A}}(t_{n})\delta\rho_{n}$, involving the number of generated
free electrons ($\delta\rho_{n}$) and the vector potential
($\mathrm{\mathbf{A}}(t_{n})$) at these instants. See text in Section IV.

Using these baseline pump fields, we numerically integrated the LC model Eq.
(5) recalled in Section II.1 for a number of cases labeled 1-5 in Fig. 4.These
cases refer to different polarization states of the SH, namely

  1. 1.

linear polarization in $x$ direction ($\theta_{2}=0$),

  2. 2.

circular polarization (clockwise, $\theta_{2}=\pi/4$),

  3. 3.

linear polarization in $y$ direction ($\theta_{2}=\pi/2$),

  4. 4.

circular polarization (anti-clockwise, $\theta_{2}=3\pi/4$),

  5. 5.

linear polarization in $x$ direction ($\theta_{2}=\pi$).

The resulting properties of the generated THz radiation are summarized in Fig.
5. In Figs. 5(a,b) the THz pulse energies defined by the integral of
$I(\omega)$ over frequencies below the cutoff value
$\omega_{\mathrm{co}}=\omega_{0}/4$ are shown. The frequency-averaged
ellipticities defined by Eq. (6) in the same frequency range are detailed in
Figs. 5(c,d). The corresponding pump polarization is illustrated in Fig. 5 by
red ellipses (for FH) and blue ellipses (for SH), indicating polarization
states of SH and FH. The THz energy is maximal when SH and FH are circularly
polarized and co-rotating (in the following we will reffer to such
configuration as ”CP-S”), in agreement with recent experimental and
theoretical investigations Meng  _et al._ (2016b); Tulsky  _et al._ (2018);
Tailliez  _et al._ (2020). Although the frequency-averaged ellipticity is
generally not very large, it seems to behave inversely to the THz energy, i.e.
the ellipticity is maximum when the THz energy is low, and vice versa.

Figure 5: THz energy (a,b) and frequency-averaged ellipticity (c,d) as a
function of the rotation angle $\theta_{2}$ for SH considering (a,c) circular
polarization and (b,d) linear polarization for FH. The vertical lines indicate
selected values of $\theta_{1}$, $\theta_{2}$, [marked by numbers in (c,d)]
for which the respective pump waveshapes are shown in Fig. 4. Ellipses
illustrate the polarization states of FH and SH.

Furthermore, typical THz waveshapes obtained for 50 fs long pump pulses, also
calculated using the LC model, are shown in Fig. 6 for different valued of
$\theta_{2}$, assuming equal amplitudes of FH and SH, $\phi=0$, and
$\theta_{1}$ to be either 0 or $\pi/4$. We can observe that the THz waveshapes
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.

To explain these observations and gain a general understanding of how THz
waveshapes 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
waveshapes).

Figure 6: Selected THz waveshapes for circular (a,b) and linear (c,d)
polarization of FH for different values of $\theta_{2}\equiv\theta$.  Figure
7: Frequency-averaged THz ellipticities $\langle\epsilon\rangle$ according to
Eq. (6) as a function of the pump pulse duration $\tau$ [see Eqs. (7)-(9)] for
$\phi=0$: (a) co-rotating circularly polarized FH and SH with equal amplitudes
(CP-S), and (b) linearly polarized FH and circularly polarized SH.

##  IV THz waveshapes in the zero frequency limit

Assuming ionization far below saturation (that is, $\rho_{e}\ll{\rho}_{\rm
at}$), and neglecting the effects of the field envelope, we can expand the
pump electric field into harmonics of the fundamental frequency $\omega_{0}$:

$\mathrm{\mathbf{E}}(t)=\sum_{m}{\mathrm{\mathbf{E}}_{m}e^{-im{\omega}_{0}t}},$
(10)

where we assume $\mathrm{\mathbf{E}}_{0}=0$ in what follows (no DC component
in the pump). Note that because $\mathrm{\mathbf{E}}(t)$ is real-valued,
$\mathrm{\mathbf{E}}_{m}=\mathrm{\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}}.$ (11)

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

${\partial}_{t}\rho_{e}\approx{\rho}_{\rm at}W,$ (12)

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}_{\rm
at}W_{0}t,$ (13)

and the coefficients ${\rho}_{l}$ are given by

${\rho}_{l}=i\frac{{\rho}_{\rm at}}{\omega_{0}l}W_{l},\quad l\neq 0.$ (14)

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 non-
periodic term in Eq. (13) conflicts with our initial assumption
$\rho_{e}\ll\rho_{\rm at}$. This inconsistency comes from neglecting the
finite pulse duration and using Fourier series expansion. However, the non-
periodic term and the coefficient $\rho_{0}$ do not play any role in the
following considerations.

With Eqs. (5), (10) and (13) we have

$\mathrm{\mathbf{E}}_{\mathrm{Br}}\propto\rho_{e}\mathrm{\mathbf{E}}=\sum_{m,l}{{\rho}_{l}e^{-il{\omega_{0}}t}}{\mathrm{\mathbf{E}}_{m}e^{-im{\omega}_{0}t}}+{\rho}_{\rm
at}W_{0}t\mathrm{\mathbf{E}}.$ (15)

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

$\mathrm{\mathbf{E}}_{\mathrm{Br},0}\propto\sum_{m\neq
0}{\mathrm{\mathbf{E}}_{m}{\rho}_{-m}}.$ (16)

Here, we can exclude $m=0$ from the summation because
$\mathrm{\mathbf{E}}_{0}=0$, and therefore the coefficient $\rho_{0}$ does not
contribute to $\mathrm{\mathbf{E}}_{\mathrm{Br},0}$. Only the harmonics of the
fundamental frequency $\omega_{0}$ contribute to
$\mathrm{\mathbf{E}}_{\mathrm{Br},0}$, which is why we call Eq. (16)
”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 growths in
the form of ”sharp steps,” Babushkin  _et al._ (2011):

$\rho_{e}(t)\approx\sum_{n}\delta\rho_{n}\Theta(t-t_{n}),$ (17)

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}),$ (18)

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})$ (19)

as a sum of ${N_{c}}$ Dirac combs with period $2\pi/\omega_{0}$. The Fourier
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})}\,.$
(20)

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}\,.$
(21)

With Eq. (11) and Eq. (12) we can conclude that

$W_{l}=\frac{\omega_{0}}{2\pi{\rho}_{\rm
at}}\sum_{n=1}^{N_{c}}\delta\rho_{n}e^{il\omega_{0}t_{n}},$ (22)

and Eq. (14) 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.$ (23)

Thus, with Eq. (16) we arrive to:

$\mathrm{\mathbf{E}}_{\mathrm{Br},0}\propto\sum_{m\neq
0}\sum_{n=1}^{N_{c}}\frac{\delta\rho_{n}\mathrm{\mathbf{E}}_{m}}{m}e^{-im\omega_{0}t_{n}}.$
(24)

This expression can be simplified using the vector potential
$\partial_{t}\mathrm{\mathbf{A}}(t)=-\mathrm{\mathbf{E}}(t).$ With

$\mathrm{\mathbf{A}}(t)=\sum_{m}\mathrm{\mathbf{A}}_{m}e^{-im\omega_{0}t}$
(25)

we find

$\mathrm{\mathbf{A}}_{m}=-\frac{i}{m\omega_{0}}\mathrm{\mathbf{E}}_{m},\quad
m\neq 0$ (26)

for the $m$th harmonic. If we further set $\mathrm{\mathbf{A}}_{0}=0$, we get

$\mathrm{\mathbf{A}}(t)=-\int_{-\infty}^{t}\mathrm{\mathbf{E}}(t^{\prime})dt^{\prime},$
(27)

and Eq. (24) can be expressed as

$\mathrm{\mathbf{E}}_{\mathrm{Br},0}\propto\sum_{m}\sum_{n=1}^{N_{c}}\delta\rho_{n}\mathrm{\mathbf{A}}_{m}e^{-im\omega_{0}t_{n}}.$
(28)

Eq. (28) as well as Eq. (24) represent the ”mixed representation”: they
contain the harmonics of the pump, $\mathrm{\mathbf{E}}_{m}$ or
$\mathrm{\mathbf{A}}_{m}$, respectively, and, at the same time, the step sizes
of ionization $\delta\rho_{n}$ in time; see Fig. 8.

Figure 8: Different representations of the zero frequency component
$\mathrm{\mathbf{E}}_{\mathrm{Br},0}$. (a) Time representation – corresponds
to Eq. (29): $\mathrm{\mathbf{E}}_{\mathrm{Br},0}$ can be determined using the
information on the vector potential $\mathrm{\mathbf{A}}(t_{n})$ in the
vicinity of the ionization events $t_{n}$ and the corresponding ionization
steps $\delta\rho_{n}$. Note that $A(t_{n})$ represents the velocity
$v(t=\infty)$ that the electrons born at $t_{n}$ have acquired after the pulse
has passed. (b) Frequency representation – corresponds to Eqs. (16) and (26)
and utilizes the polarization states $\mathrm{\mathbf{A}}_{m}$ and the
electron density Fourier coefficients $\rho_{m}$ of the $m$th harmonic. (c)
Mixed representation – corresponds to Eq. (28), using
$\mathrm{\mathbf{A}}_{m}$ in frequency space and $\delta\rho_{n}$ in the time
domain, together with the phase information on the ionization event
$\exp\\!\left(imt_{n}\omega_{0}\right)$.

To further understand the physical meaning of
$\mathrm{\mathbf{E}}_{\mathrm{Br},0}$, we use Eq. (25) to obtain from Eq. (28)
the ”time representation”, cf. Ref Debayle  _et al._ (2014a),

$\mathrm{\mathbf{E}}_{\mathrm{Br},0}\propto\sum_{n=1}^{N_{c}}\delta\rho_{n}\mathrm{\mathbf{A}}(t_{n}).$
(29)

This equation allows to predict which pump waveshapes generate THz radiation
effectively Babushkin  _et al._ (2011): for each ionization event, one must
optimize $\delta\rho_{n}$ and $\mathrm{\mathbf{A}}(t_{n})$ simultaneously,
which is not trivial: Since the $t_{n}$ are located at the extrema of
$\mathrm{\mathbf{E}}(t)$, for simple (e.g., single-color) pump waveshapes,
$|\mathrm{\mathbf{A}}(t_{n})|\approx 0$, rendering
$\delta\rho_{n}\mathrm{\mathbf{A}}(t_{n})$ negligible. Examples of two-color
pump waveshapes with $\delta\rho_{n}\mathrm{\mathbf{A}}(t_{n})$ for each
ionization event (red arrows) are shown in Fig. 4. Comparing, for instance,
the cases Fig. 4(a).2 and Fig. 4(b).2, one can see that even though the
amplitude of $\mathrm{\mathbf{E}}(t_{n})$ is the same and thus
$\delta\rho_{n}$ is similar, $\delta\rho_{n}\mathrm{\mathbf{A}}(t_{n})$ is
larger in case (a).2 due to the larger amplitude of
$\mathrm{\mathbf{A}}(t_{n})$. For certain pump configurations, there is more
than one significant ionization event per cycle, e.g., for case 4 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
$\mathrm{\mathbf{A}}(t_{n})$ are pointing in different directions. As a
result, the maximum THz yield is achieved when one ionization event dominates.
For a circularly polarized two-color pump pulse with corotating SH and FH [see
Fig. 4(a).2], this property clearly appears and the large THz yield is
confirmed in Fig. 5(a) Meng  _et al._ (2016b); Tulsky  _et al._ (2018);
Tailliez  _et al._ (2020).

The physics behind Eq. (29) is easy to understand, as depicted in Fig. 8: The
electron, leaving the atom with negligible initial velocity, experiences
acceleration in the field $d\mathrm{\mathbf{v}}/dt\propto\mathrm{\mathbf{E}}$.
Integrating this expression from the electron birth time $t_{n}$ to
$t=\infty$, we obtain the electron velocity after the pulse,

$\mathrm{\mathbf{v}}(t\to\infty,t_{n})\propto\int_{t_{n}}^{\infty}\mathrm{\mathbf{E}}(t^{\prime})dt^{\prime}.$
(30)

Because

$0=\mathrm{\mathbf{E}}_{0}=\int_{-\infty}^{\infty}\mathrm{\mathbf{E}}(t^{\prime})dt^{\prime}=\int_{-\infty}^{t_{n}}\mathrm{\mathbf{E}}(t^{\prime})dt^{\prime}+\int_{t_{n}}^{\infty}\mathrm{\mathbf{E}}(t^{\prime})dt^{\prime}$

we find with Eq. (27) that

$\mathrm{\mathbf{v}}(t\to\infty,t_{n})\propto\mathrm{\mathbf{A}}(t_{n}).$ (31)

The fact that the electron velocity is proportional to the vector potential at
the time of ionization is well known Mori and Katsouleas (1992). It is the
heart of the so-called atto-streaking technique Goulielmakis  _et al._ (2004,
2008), 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 by the pulse reads

$\Delta\mathrm{\mathbf{J}}=\sum_{n=1}^{N_{c}}\delta\rho_{n}\mathrm{\mathbf{v}}(t\to\infty,t_{n})\propto\sum_{n=1}^{N_{c}}\delta\rho_{n}\mathrm{\mathbf{A}}(t_{n}),$
(32)

which was already exploited in the context of THz generation Babushkin  _et
al._ (2011); Debayle  _et al._ (2014b); de Alaiza Martínez  _et al._ (2015);
Debayle  _et al._ (2015)111The relevant quantity in Babushkin  _et al._ (2011)
is
$v_{f}(t)=\int_{-\infty}^{t}E(t^{\prime})e^{\gamma(t-t^{\prime})}dt^{\prime}$,
where $\gamma$ describes collision-induced decay of the current. Neglecting
$\gamma$ restores the vector potential $A(t)$ on the right side of this
expression.. Since we assume that there is no net current density before the
pulse arrives, $\mathrm{\mathbf{J}}(t\to-\infty)=0$, we find

$\mathrm{\mathbf{E}}_{\mathrm{Br},0}\propto\Delta\mathrm{\mathbf{J}}=\mathrm{\mathbf{J}}(t\to\infty)-\mathrm{\mathbf{J}}(t\to-\infty)=\mathrm{\mathbf{J}}(t\to\infty).$

Because the current density $\mathrm{\mathbf{J}}$ is a real-valued vector, the
0th harmonic $\mathrm{\mathbf{E}}_{\mathrm{Br},0}$ must be linearly polarized.
This follows directly from the Jones formalism, in which any elliptically (or
circularly) polarized field requires a complex-valued description. 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.

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.
(24) we find

$\begin{split}\mathrm{\mathbf{E}}_{\mathrm{Br},0}\propto\sum_{n=1}^{N_{c}}\delta\rho_{n}&\Big{(}\mathrm{\mathbf{E}}_{1}e^{-i\omega_{0}t_{n}}-\mathrm{\mathbf{E}}_{-1}e^{i\omega_{0}t_{n}}\\\
&+\frac{\mathrm{\mathbf{E}}_{2}}{2}e^{-2i\omega_{0}t_{n}}-\frac{\mathrm{\mathbf{E}}_{-2}}{2}e^{2i\omega_{0}t_{n}}\Big{)}.\end{split}$
(33)

To see how this ”mixed representation” can be used, let us consider the
particular waveshape Eqs. (7)-(9) with
$\epsilon_{\omega_{0}}=\epsilon_{2\omega_{0}}=1$ and $\phi=0$. This is the
well-known CP-S waveshape, with both FH and SH being circularly polarized and
co-rotating Tailliez  _et al._ (2020); Meng  _et al._ (2016b); Song  _et al._
(2020); Babushkin  _et al._ (2022). Here, we consider this pump waveshape
through our analytical approach. In Fig. 9, 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\rightarrow+\infty$. Then, the
pump harmonics amplitudes read

$\displaystyle\mathrm{\mathbf{E}}_{1}=\frac{e_{\omega_{0}}}{2}\left(1,+i\right),\,\mathrm{\mathbf{E}}_{-1}=\frac{e_{\omega_{0}}}{2}\left(1,-i\right),$
(34)
$\displaystyle\mathrm{\mathbf{E}}_{2}=\frac{e_{2\omega_{0}}}{2}\left(1,+i\right),\,\mathrm{\mathbf{E}}_{-2}=\frac{e_{2\omega_{0}}}{2}\left(1,-i\right),$
(35)

where $(a,b)$ on the right-hand side represents a Jones vector with the
corresponding $x$ and $y$ components. Substituting Eqs. (34) and (35) into Eq.
(24) gives

$\mathrm{\mathbf{E}}_{\mathrm{Br},0}\propto(0,1).$ (36)

Thus, the 0th harmonic for this pump configuration is linearly polarized in
the $y$-direction. More generally, for an arbitrary phase $\phi$ between the
co-rotating 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.
(35) has to be modified into

$\mathrm{\mathbf{E}}_{2}=\frac{e_{2\omega_{0}}}{2}\left(1,+i\right)e^{-i\phi},\,\mathrm{\mathbf{E}}_{-2}=\frac{e_{2\omega_{0}}}{2}\left(1,-i\right)e^{i\phi}.$
(37)

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

$\mathrm{\mathbf{E}}_{br,0}\propto\left(\sin{\phi},\cos{\phi}\right).$ (38)

This well-known property Dai  _et al._ (2009); Tailliez  _et al._ (2020); Meng
_et al._ (2016b); Song  _et al._ (2020) is illustrated in Fig. 9(a-c). In the
next section, more general examples are addressed.

Figure 9:  (a,b,d-f) Pump waveshapes (black lines) for the CP-S configuration
(a,b) and linearly polarized FH + circularly polarized SH (d-f) for an
exemplary pulse duration of 50 fs and different phases $\phi$ between FH and
SH. The associated THz waveshapes are shown by red lines (their amplitudes are
rescaled to make them comparable to the pump fields). (c) Numerically-computed
THz waveshapes for the CP-S case for various values of $\phi$ (solid lines)
compared to their corresponding analytical solution of Eq. (38) (dashed
lines).

##  V THz waveshapes 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 non-zero frequencies, the situation is more involved. In Fig.
10, 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”. As we show later, this 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 waveshape for (low) non-zero
frequencies, and present closed analytical expressions for the linear
ellipticity chirp for some selected pump waveshapes.

###  V.1 THz polarization state for generic multi-color pump pulses

For simplicity, we assume $\rho_{e}\ll\rho_{\rm at}$ as in the previous
section, but now we have to consider the finite duration of the pump pulse. In
this situation, the pump waveshape 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. Since
$\mathrm{\mathbf{E}}_{\mathrm{Br}}(t)\propto\rho_{e}(t)\mathrm{\mathbf{E}}(t)$,
we find a convolution integral in the Fourier domain,

$\hat{\mathrm{\mathbf{E}}}_{\mathrm{Br}}(\omega)\propto\int_{-\infty}^{\infty}\hat{\rho}_{e}(\omega-\Omega)\hat{\mathrm{\mathbf{E}}}(\Omega)d\Omega,$
(39)

where $\,\widehat{}\,$ symbol refers to the Fourier transform in time. From
Eq. (18) we immediately get

$-i\omega\hat{\rho}_{e}(\omega)\approx\sum_{n=1}^{N}\delta\rho_{n}e^{i\omega
t_{n}},$ (40)

which implies that $\hat{\rho}_{e}(\omega)$ is rather broad in frequency
space. In contrast, we assume that the multi-color pump field
$\hat{\mathrm{\mathbf{E}}}(\omega)$ is given by narrow peaks at the harmonic
frequencies. These narrow peaks will be approximated by Dirac $\delta$
functions in the following, and we write

$\hat{\mathrm{\mathbf{E}}}(\omega)\approx\sum_{m}\mathrm{\mathbf{E}}_{m}\delta(\omega-m\omega_{0}),$
(41)

setting $\mathrm{\mathbf{E}}_{0}=0$. We note that the approximation Eq. (41)
renders the pump field similar to Eq. (10) used in the previous section.
However, we must remember that Eq. (41) is only an approximation to facilitate
further computations. In this section, we will not use any discrete Fourier
series expansion of periodic waveforms but only continuous Fourier transforms
of finite pulses.

Figure 10: (a) Spectra of THz radiation, (b) THz spectral ellipticity, and (c)
orientation of the THz polarization ellipse 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 amplitudes of the two pump harmonics are equal, their durations are
50 fs, and $\phi=0$.

With Eqs. (39)-(41), the Brunel radiation in the THz range can be expressed as

$\begin{split}\hat{\mathrm{\mathbf{E}}}_{\mathrm{Br}}(\omega)&\propto\sum_{n=1}^{N}\delta\rho_{n}\int\hat{\mathrm{\mathbf{E}}}(\Omega)\frac{e^{i(\omega-\Omega)t_{n}}}{\omega-\Omega}d\Omega\\\
&=\sum_{m}\sum_{n=1}^{N}\delta\rho_{n}\mathrm{\mathbf{E}}_{m}\frac{e^{i(\omega-m\omega_{0})t_{n}}}{\omega-m\omega_{0}}.\end{split}$
(42)

This is, in general, a complex function of $\omega$. However, for small
$\omega$ one can approximate it as

$\hat{\mathrm{\mathbf{E}}}_{\mathrm{Br}}(\omega)\propto\boldsymbol{\mu}_{0}+\boldsymbol{\mu}_{1}\frac{\omega}{\omega_{0}}+\ldots,$
(43)

where the vector-valued coefficients $\boldsymbol{\mu}_{0}$ and
$\boldsymbol{\mu}_{1}$ are

$\displaystyle\boldsymbol{\mu}_{0}$
$\displaystyle=-\sum_{m}\sum_{n=1}^{N}\delta\rho_{n}\mathrm{\mathbf{E}}_{m}\frac{e^{-im\omega_{0}t_{n}}}{m\omega_{0}},$
(44) $\displaystyle\boldsymbol{\mu}_{1}$
$\displaystyle=-\sum_{m}\sum_{n=1}^{N}\delta\rho_{n}\mathrm{\mathbf{E}}_{m}\frac{e^{-im\omega_{0}t_{n}}\left(1+im\omega_{0}t_{n}\right)}{m^{2}\omega_{0}}.$
(45)

Alternatively, in terms of $\mathrm{\mathbf{A}}$ we can write, cf. Eq. (26):

$\displaystyle\boldsymbol{\mu}_{0}$
$\displaystyle=-i\sum_{m}\sum_{n=1}^{N}\delta\rho_{n}\mathrm{\mathbf{A}}_{m}e^{-im\omega_{0}t_{n}},$
(46) $\displaystyle\boldsymbol{\mu}_{1}$
$\displaystyle=-i\sum_{m}\sum_{n=1}^{N}\delta\rho_{n}\mathrm{\mathbf{A}}_{m}\frac{e^{-im\omega_{0}t_{n}}\left(1+im\omega_{0}t_{n}\right)}{m}.$
(47)

With these expressions, it is possible to compute the ellipticity
$\epsilon(\omega)$ of the THz polarization ellipse at each frequency. We note
that the zero-order term Eq. (46) coincides with Eq. (28), except that here we
sum over all ionization events of the pulse. Thus, for $\omega=0$, we obtain a
linearly polarized THz field, consistent with the previous section.

###  V.2 THz ellipticity $\epsilon(\omega)$ for the multi-color case

When taking into account zero ellipticity at $\omega=0$ and neglecting higher
order terms in $\omega$, Eq. (43) predicts a linearly increasing ellipticity
with $\omega$:

$\epsilon(\omega)\approx\mathcal{B}\frac{\omega}{\omega_{0}},\quad\omega\ll\omega_{0},$
(48)

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 first recall that since the electric field satisfies
$\mathrm{\mathbf{E}}_{m}=\mathrm{\mathbf{E}}_{-m}^{*}$ and
$\mathrm{\mathbf{E}}_{0}=0$, we can rewrite Eqs. (44) and (45) as

$\displaystyle\boldsymbol{\mu}_{0}$
$\displaystyle=-2i\sum_{n=1}^{N}\delta\rho_{n}\sum_{m>0}\Im\left[\frac{\mathrm{\mathbf{E}}_{m}}{m\omega_{0}}e^{-im\omega_{0}t_{n}}\right],$
(49) $\displaystyle\boldsymbol{\mu}_{1}$
$\displaystyle=-2\sum_{n=1}^{N}\delta\rho_{n}\sum_{m>0}\Re\left[\frac{\mathrm{\mathbf{E}}_{m}(1+im\omega_{0}t_{n})e^{-im\omega_{0}t_{n}}}{m^{2}\omega_{0}}\right].$
(50)

Therefore, $\boldsymbol{\mu}_{0}$ is pure 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,

${\cal
B}=\frac{{|\boldsymbol{\mu}}_{0}\times{\boldsymbol{\mu}}_{1}|}{|{\boldsymbol{\mu}}_{0}|^{2}},$
(51)

or equivalently:

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

For elliptically polarized multi-color pump fields, $|\mathcal{B}|>0$ and
therefore non-zero frequency components of the generated THz pulse are
generally elliptically polarized.

###  V.3 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. (7)–(9). 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. As illustrative examples,
we will address two polarization geometries: A CP-S polarization state (co-
rotating FH and SH) where $\epsilon_{\omega_{0}}=1$, and a
$\mathrm{\mathbf{x}}$ linearly polarized FH coupled to a circularly polarized
SH where $\epsilon_{\omega_{0}}=0$. Examples of pump and numerically
calculated THz waveshapes for different values of the phase $\phi$ between FH
and FH are shown in Fig. 9(a-c) for the former (CP-S) and Fig. 9(d-f) for the
latter case.

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

$\displaystyle\mathrm{\mathbf{E}}_{1}=\frac{e_{\omega_{0}}}{2}\left(1,+i\epsilon_{\omega_{0}}\right),\,\mathrm{\mathbf{E}}_{-1}=\frac{e_{\omega_{0}}}{2}\left(1,-i\epsilon_{\omega_{0}}\right),$
(53)
$\displaystyle\mathrm{\mathbf{E}}_{2}=\frac{e_{2\omega_{0}}}{2}\left(1,+i\right)\mbox{e}^{-i\phi},\,\mathrm{\mathbf{E}}_{-2}=\frac{e_{2\omega_{0}}}{2}\left(1,-i\right)\mbox{e}^{+i\phi}.$
(54)

Evaluating Eqs. (44) and (45) yields

$\displaystyle\boldsymbol{\mu}_{0}$
$\displaystyle=-i\sum_{n=1}^{N}\frac{\delta\rho_{n}}{\omega_{0}}\left[e_{\omega_{0}}\begin{pmatrix}-\sin{(\omega_{0}t_{n})}\\\
\epsilon_{\omega_{0}}\cos{(\omega_{0}t_{n})}\end{pmatrix}\right.$
$\displaystyle+\left.\frac{e_{2\omega_{0}}}{2}\begin{pmatrix}-\sin{(2\omega_{0}t_{n}+\phi)}\\\
\cos{(2\omega_{0}t_{n}+\phi)}\end{pmatrix}\right],$ (55)
$\displaystyle\boldsymbol{\mu}_{1}$
$\displaystyle=-\sum_{n=1}^{N}\frac{\delta\rho_{n}}{\omega_{0}}\left[{e_{\omega_{0}}}\begin{pmatrix}\cos{(\omega_{0}t_{n})}+\omega_{0}t_{n}\sin{(\omega_{0}t_{n})}\\\
\epsilon_{\omega_{0}}\left[\sin{(\omega_{0}t_{n})}-\omega_{0}t_{n}\cos{(\omega_{0}t_{n})}\right]\end{pmatrix}\right.$
$\displaystyle+\left.\frac{e_{2\omega_{0}}}{4}\begin{pmatrix}\cos{(2\omega_{0}t_{n}+\phi)}+2\omega_{0}t_{n}\sin{(2\omega_{0}t_{n}+\phi)}\\\
\sin{(2\omega_{0}t_{n}+\phi)}-2\omega_{0}t_{n}\cos{(2\omega_{0}t_{n}+\phi)}\end{pmatrix}\right].$
(56)

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 $|\mathrm{\mathbf{E}}|$. For our two-color pulse, one has

$\displaystyle|\mathrm{\mathbf{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)}\right.$
$\displaystyle+\left.\frac{1-\epsilon_{\omega_{0}}^{2}}{2}\cos{(2\omega_{0}t)}+r(1-\epsilon_{\omega_{0}})\cos{(3\omega_{0}t+\phi)}\right),$
(57)

where $r\equiv e_{2\omega_{0}}/e_{\omega_{0}}$. The ionization instants
$t_{n}$ are determined by the roots of
$\partial_{t}|\mathrm{\mathbf{E}}(t)|^{2}=0$, from which we have to select
those that satisfy $\partial^{2}_{t}|\mathrm{\mathbf{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.$ (58)

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

$\frac{\partial}{\partial
t}\left\\{\cos{(\omega_{0}t)}[\cos{(\omega_{0}t)}+2r\cos{(2\omega_{0}t+\phi)}]\right)\\}=0.$
(59)

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

$\displaystyle\phi=0$
$\displaystyle\quad\rightarrow\quad\omega_{0}t_{n}=2n\pi,$ (60)
$\displaystyle\phi=\frac{\pi}{4}$
$\displaystyle\quad\rightarrow\quad\omega_{0}t_{n}=\arctan\left(\frac{\sqrt{2}-\sqrt{6}}{\sqrt{2}+\sqrt{6}}\right)+2n\pi,$
(61) $\displaystyle\phi=\frac{\pi}{2}$
$\displaystyle\quad\rightarrow\quad\omega_{0}t_{n}=n\pi-(-1)^{n}\frac{\pi}{6}.$
(62)

Whereas for $\phi=0$ and $\phi=\pi/4$ there are only one strong ionization
event over the cycle, for the case $\phi=\pi/2$ there are two events with
equal amplitudes, given by $n=2\pi n$ and $n=(2n+1)\pi$ (see insets in Fig.
11). We note that the influence of the phase angle $\phi$ on the laser-to-THz
conversion efficiency was studied extensively in Refs. Kim  _et al._ (2007);
Babushkin  _et al._ (2011); Oh  _et al._ (2012); Alirezaee and Sharifian
(2018); Zhang  _et al._ (2016). Here, we complement these findings by its
impact on the THz waveforms and ellipticities. In the following, we consider
two distinct situations: zero or non-zero relative phase $\phi$.

####  V.3.1 Zero relative phase

For $\phi=0$, Eqs. (58) and (60) show that ionization events occur at
$t_{n}=0$ modulo $2\pi/\omega_{0}$. In the CP-S geometry
($\epsilon_{\omega_{0}}=1$), one finds from Eqs. (55) and (56):

$\begin{split}\boldsymbol{\mu}_{0}&=-\frac{1}{\omega_{0}}\sum_{n=1}^{N}\delta\rho_{n}\left(e_{\omega_{0}}+\frac{e_{2\omega_{0}}}{2}\right)\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}\\\
-2n\pi\left[{e_{\omega_{0}}}+\frac{e_{2\omega_{0}}}{2}\right]\end{pmatrix},\end{split}$
(63)

and the slope ${\cal B}$, cf. Eq. (52), is given by

${\cal
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|.$
(64)

For FH and SH with equal amplitudes this yields ${\cal B}=5/6$, that is,

$\epsilon(\omega)=\frac{5}{6}\frac{\omega}{\omega_{0}}\quad\textrm{for}\quad
e_{\omega_{0}}=e_{2\omega_{0}},$ (65)

which slightly corrects the result obtained in Babushkin  _et al._ (2022).

In the hybrid polarization geometry (FH is LP; SH is CP –
$\epsilon_{\omega_{0}}=0$), we find

$\begin{split}\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{split}$ (66)

and

${\cal
B}=\left|\frac{\mu_{1x}}{\mu_{0y}}\right|=\left|\frac{e_{\omega_{0}}+e_{2\omega_{0}}/4}{e_{\omega_{0}}/2}\right|.$
(67)

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

$\epsilon(\omega)=\frac{5}{2}\frac{\omega}{\omega_{0}}\quad\textrm{for}\quad
e_{\omega_{0}}=e_{2\omega_{0}}.$ (68)

In Fig. 10, we compare the ellipticities given by Eq. (65) and Eq. (68) 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 a pulse with 50 fs duration and
100 TW/cm2 intensity and circular polarization of SH, combined with linear or
circular polarization of FH. As observed in Fig. 10(b), the analytical
solutions given by Eq. (65) and Eq. (68), which are indicated in red and black
lines for the cases of circular and linear polarization of FH, respectively,
perfectly match the simulation results represented by the circle and square
symbols. In contrast to ellipticity, which changes with frequency, the
polarization direction for small frequencies is constant. This is supported by
Fig. 10(c), where the 90 degree polarization direction, given by the direction
of $\mathrm{\mathbf{A}}(t_{n})$ [cf. also Fig. 8(a)], coincides with the
analytical predictions.

####  V.3.2 Non-zero relative phase

When $\phi$ is non-zero, Eq. (58) 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. (52) leads to the $\phi$-independent expression
of the THz ellipticity:

${\cal
B}=\frac{|\mu_{1y}(\mu_{0x})-\mu_{1x}(\mu_{0y})|}{|\boldsymbol{\mu}_{0}|^{2}}=\left|\frac{e_{\omega_{0}}+e_{2\omega_{0}}/4}{e_{\omega_{0}}+e_{2\omega_{0}}/2}\right|.$
(69)

Thus, the ellipticity does not change with varying $\phi$, which is in
agreement with Fig. 9(a-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. (61) and Eq. (62), analytical results for
FH and SH with equal amplitudes ($e_{\omega_{0}}=e_{2\omega_{0}}$) can be
found and Eq. (52) yields

$\displaystyle\phi=\frac{\pi}{4}$ $\displaystyle\quad\rightarrow\quad{\cal
B}=1+\frac{3\sqrt{3}}{4},$ (70) $\displaystyle\phi=\frac{\pi}{2}$
$\displaystyle\quad\rightarrow\quad{\cal B}=\frac{\pi}{\sqrt{3}}-\frac{1}{2}.$
(71)

Thus, from the slope $5/2$ achieved when $\phi=0$, the THz ellipticity
decreases for non-zero $\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. 11. We
note that the numerical solution of the LC model Eq. (5) yields a
significantly lower value for $\cal{B}$ for $\phi=\pi/2$ than our analytical
model Eq. (52). 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. 11), which are not taken into account in Eq. (52).

However, Eq. (52) predicts the right trend, namely a decrease of $\cal{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. (51) large. This decrease in $\cal{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.
(46)-(47). In particular, Eq. (43) implies that

$\boldsymbol{\mu}_{0}\propto\mathrm{\mathbf{E}}_{\mathrm{Br,0}},$ (72)

where $\mathrm{\mathbf{E}}_{\mathrm{Br,0}}$ is the THz field at $\omega\to 0$.
Thus, contributions to $\boldsymbol{\mu}_{0}$ from different ionization events
within the optical cycle may cancel each other out, rendering
$|\mathrm{\mathbf{\mu}}_{0}|$ very small [see previous discussion about
possible vanishing of $\mathrm{\mathbf{E}}_{\mathrm{Br,0}}$ below Eq. (29)].
By contrast, $\boldsymbol{\mu}_{1}$ is not affected by this cancelation
because with

$\int_{-\infty}^{t}\mathrm{\mathbf{A}}(t^{\prime})dt^{\prime}=\frac{i}{\omega_{0}}\sum_{m}\frac{\mathrm{\mathbf{A}}_{m}}{m}e^{-im\omega_{0}t}$
(73)

we find

$\boldsymbol{\mu}_{1}=\omega_{0}\sum_{n=1}^{N}\delta\rho_{n}\left[t_{n}\mathrm{\mathbf{A}}(t_{n})-\int_{-\infty}^{t_{n}}\mathrm{\mathbf{A}}(t^{\prime})dt^{\prime}\right]\neq
0.$ (74)

Therefore, we observe a significantly larger frequency-averaged ellipticity
for cases with small $\boldsymbol{\mu}_{0}$ in the numerical solutions of the
LC model Eq. (5). That is, when the THz energy is small, the ellipticity is
large. This behaviour can indeed be observe in Fig. 5.

Figure 11: Numerical results obtained from the LC model (lines and circles):
(a) Relative THz intensities, (b) values of ellipticity chirp factor
$\mathcal{B}$ and (c) polarization angle, functions of the relative phase
$\phi$ between FH and SH, for FH being linearly and SH circularly polarized.
Stars show analytical predictions. Insets show the pump waveshapes for the
same values of $\phi$ as in Fig. 9(d-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 “extend”
of the ionization events (full-width-half-maximum of the ionization rate’s
peaks are indicated). Note the asymmetry of the field shape around the center
of ionization event in all cases except for $\phi=0$.

###  V.4 Frequency-averaged THz ellipticity

The previous analysis suggests a straightforward physical explanation for the
non-zero 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_{\mathrm{co}}>\Delta$. The dependence of the ellipticity of THz pulses
on the duration of the pump pulse was already reported in Babushkin  _et al._
(2011); Wang  _et al._ (2009); Lu  _et al._ (2015). Using the equations
derived above, 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 Borodin  _et al._
(2013)). 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\limits_{0}^{\Delta}\epsilon(\omega)d\omega=\frac{\mathcal{B}}{\omega_{0}\tau}.$
(75)

This estimate is in good agreement with Fig. 7, where the dependencies
$\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 waveshapes generated via the ionization-
based Brunel mechanism by two- and multi-color pump pulses with nontrivial
polarization. By numerically solving the local current model and analytical
evaluations, we have shown that photo-ionization-driven THz waveshapes 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
generally linearly dependent on the frequency, i.e. the THz waveshapes possess
a “linear ellipticity chirp”, whose magnitude $\mathcal{B}$ depends on the
particular driving pulse. This behavior sets ionization-based THz emission
appart from other schemes, which tend to produce, at least in the first
approximation, pulses with constant ellipticity across the frequency
bandwidth.

The ellipticity chirp $\mathcal{B}$ tends to strongly peak for the pump
waveshapes, where free electron currents from several ionization events cancel
each other. Since the net THz energy in such situations is low, as a rule of
thumb, the ellipticity chirp is high for low THz energy and vise versa (for
comparable free electron densities). Pulses with non-zero ellipticity chirp
$\mathcal{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 direction
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 $\mathrm{\mathbf{A}}(t_{n})$ at the ionization
instants $t_{n}$ ($n=1,2,\ldots$).

Our results demonstrate the unique properties of THz waveforms produced in gas
plasmas by polarization controled multi-color 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
waveshapes by superimposing different locally generated waveforms Köhler  _et
al._ (2011); Zhang  _et al._ (2018), giving rise to a even larger diversity of
THz waveshapes.

## Acknowledgement

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.

## References

  * Chan  _et al._ (2007) W. L. Chan, J. Deibel, and D. M. Mittleman, “Imaging with terahertz radiation,” Rep. Prog. Phys. 70, 1325 (2007).
  * Tonouchi (2007) M. Tonouchi, “Cutting-edge terahertz technology,” Nat. Photonics 1, 97 (2007).
  * Marx (2007) B. Marx, “Terahertz technology detects counterfeit drugs,” Laser Focus World 43, 44 (2007).
  * Dhillon  _et al._ (2017) 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: Appl. Phys. 50, 043001 (2017).
  * Ma  _et al._ (2022) 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).
  * Cocker  _et al._ (2013) 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. Photonics 7, 620 (2013).
  * Katletz  _et al._ (2012) 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).
  * Babushkin  _et al._ (2022) 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. Physics 18, 417 (2022).
  * Vaičaitis  _et al._ (2023) 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).
  * Liu  _et al._ (2010) 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. Photonics 4, 627 (2010).
  * Wang  _et al._ (2010) 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).
  * Brown (2003) E. R. Brown, “Fundamentals of terrestrial millimeter-wave and THz remote sensing,” International journal of high speed electronics and systems 13, 995 (2003).
  * Pickwell and Wallace (2006) E. Pickwell and V. P. Wallace, “Biomedical applications of terahertz technology,” J. Phys. D: Appl. Phys. 39, R301 (2006).
  * Kress  _et al._ (2004) M. Kress, T. Löffler, S. Eden, M. Thomson, and H. G. Roskos, “Terahertz-pulse generation by photoionization of airwith laser pulses composed of both fundamental and second-harmonicwaves,” Opt. Lett. 29, 1120 (2004).
  * Bartel  _et al._ (2005) 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).
  * Kim  _et al._ (2007) 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).
  * Kim  _et al._ (2008) K. Y. Kim, A. G. Taylor, A. G. Glownia, and G. Rodriguez, “Coherent control of terahertz supercontinuum generation in ultrafast laser-gas interactions,” Nat. Photonics 2, 605 (2008).
  * Koulouklidis  _et al._ (2020) 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).
  * Alirezaee and Sharifian (2018) 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).
  * Babushkin  _et al._ (2011) 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).
  * Clerici  _et al._ (2013) 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).
  * Vvedenskii  _et al._ (2014) 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).
  * Nguyen  _et al._ (2017) A Nguyen, P González de Alaiza Martínez, and J Déchard, “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).
  * Babushkin  _et al._ (2017) 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).
  * Meng  _et al._ (2016a) 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 (2016a).
  * Zhang  _et al._ (2018) 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. Photonics 12, 554 (2018).
  * You  _et al._ (2013) 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).
  * Dai  _et al._ (2009) 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).
  * Fedorov  _et al._ (2016) 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 (2016).
  * Tulsky  _et al._ (2018) 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).
  * Wen and Lindenberg (2009) H. Wen and A. M. Lindenberg, “Coherent terahertz polarization control through manipulation of electron trajectories,” Phys. Rev. Lett. 103, 023902 (2009).
  * Kosareva  _et al._ (2018) 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).
  * Tailliez  _et al._ (2020) 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).
  * Stathopulos  _et al._ (2021) 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).
  * Meng  _et al._ (2016b) 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 (2016b).
  * Song  _et al._ (2020) 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).
  * Lu  _et al._ (2015) 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).
  * Borodin  _et al._ (2013) 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).
  * Bai  _et al._ (2012) 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).
  * de Alaiza Martínez  _et al._ (2015) 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).
  * Zhang  _et al._ (2017) 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,” Physics of Plasmas 24, 023116 (2017).
  * Wang  _et al._ (2023) H. Wang, W. Fan, X. Chen, and H. Yan, “Polarization control of terahertz waves generated by a femtosecond three-color pulse scheme,” JOSA B 40, 1375 (2023).
  * Balčiūnas  _et al._ (2015) 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).
  * Zhang  _et al._ (2016) 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).
  * Lu  _et al._ (2017) 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).
  * Vaičaitis  _et al._ (2019) 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 (2019).
  * Alirezaee  _et al._ (2020) 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).
  * Liu  _et al._ (2020) 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).
  * Ma  _et al._ (2021) 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).
  * Watanabe (2018) S. Watanabe, “Terahertz polarization imaging and its applications,” Photonics 5, 58 (2018).
  * Baierl  _et al._ (2016) 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. Photonics 10, 715 (2016).
  * Hoshina  _et al._ (2011) 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).
  * Amer  _et al._ (2005) 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).
  * Köhler  _et al._ (2011) 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).
  * Cabrera-Granado  _et al._ (2015) Eduardo Cabrera-Granado, Yxing Chen, Ihar Babushkin, Luc Bergé, and Stefan Skupin, “Spectral self-action of thz emission from ionizing two-color laser pulses in gases,” New J. Phys. 17, 023060 (2015).
  * Brunel (1990) 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).
  * Landau and Lifshitz (1981) L. D. Landau and L. M. Lifshitz,  _Quantum Mechanics Non-Relativistic Theory, Third Edition: Volume 3_ , 3rd ed. (Butterworth-Heinemann, 1981).
  * Thiele  _et al._ (2018) 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).
  * Debayle  _et al._ (2014a) A. Debayle, L. Gremillet, L. Bergé, and Ch. Köhler, “Analytical model for THz emissions induced by laser-gas interaction,” Opt. Express 22, 13691 (2014a).
  * Mori and Katsouleas (1992) 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).
  * Goulielmakis  _et al._ (2004) 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).
  * Goulielmakis  _et al._ (2008) 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).
  * Debayle  _et al._ (2014b) A. Debayle, L. Gremillet, L. Bergé, and Ch. Köhler, “Analytical model for THz emissions induced by laser-gas interaction,” Opt. Express 22, 13691 (2014b).
  * Debayle  _et al._ (2015) 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).
  * Oh  _et al._ (2012) 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).
  * Wang  _et al._ (2009) 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).