# Waveshape of THz radiation produced by two-color laser-induced air plasmas

A. Stathopulos  CEA, DAM, DIF, 91297 Arpajon, France  Université Paris-Saclay,
CEA, LMCE, 91680 Bruyères-le-Châtel, France  S. Skupin  Institut Lumière
Matière, UMR 5306 - CNRS, Université de Lyon 1, 69622 Villeurbanne, France  B.
Zhou  DTU Electro - Department of Electrical and Photonics Engineering,
Technical University of Denmark, 2800 Kongens Lyngby, Denmark  P. U. Jepsen
DTU Electro - Department of Electrical and Photonics Engineering, Technical
University of Denmark, 2800 Kongens Lyngby, Denmark  L. Bergé
<EMAIL_ADDRESS>Centre des Lasers Intenses et Applications, Université
de Bordeaux-CNRS-CEA, 33405 Talence Cedex, France

(August 28, 2024)

###### Abstract

The spatial and spectral distributions of terahertz (THz) pulses emitted by
two-color air plasmas are theoretically investigated for focused laser pulses
and in the filamentation regime. We derive a so-called ”augmented” conical
emission model, which, similarly to the one originally proposed by You et al.
[Phys. Rev. Lett. 109, 183902 (2012)], involves phase matching between laser
harmonics along the plasma channel, the plasma density and length, and the
emitted frequency as key parameters. Our augmented model, however, treats
envelope effects and accounts for transverse variations of both plasma- and
Kerr-driven potential THz emitters. We highlight the importance of the
characteristic spatio-spectral distributions of these two conversion
mechanisms in the expression of the angular radiated power. The results of our
model are successfully compared with data provided by a comprehensive, fully
space and time-resolved unidirectional solver. Importantly, these numerical
simulations clear up the effective plasma length along which THz emission
develops, compared with the dephasing length along which the laser fundamental
and second harmonic become out-of-phase. The impact of common optical
aberrations, such as sphericity, astigmatism, and coma, on the THz generation
is also investigated. Aberrations are shown to generally decrease the laser-
to-THz conversion efficiency and potentially induce spatial asymmetries and
narrowing in the THz spectra.

##  I Introduction

The production of terahertz (THz) radiation is receiving increasing interest
in a wide range of domains, such as homeland security [1, 2, 3], medical
diagnosis [4, 5, 6], art conservation [7, 8, 9, 10] or archeology [11, 12].
THz waves can propagate through a few millimeters of optically opaque
materials, provided the water concentration is low [13, 14], and can be used
to identify complex molecules by coherent THz spectroscopy [15, 16, 17, 18].
Other notable applications of THz waves include control of leaf water content
[19, 20, 21], plasma excitation [22], electron acceleration [23], high
harmonic (sideband) generation [24, 25, 26], or the fabrication of
ferromagnetic structures [27]. Several emitters can nowadays produce intense
pulsed or continuous-wave THz radiation, including photoconductive antennas
[28, 29], molecular (CH3OH) gas sources [30], quantum cascade lasers [31, 32,
33] or $\chi^{(2)}$ crystals serving for optical rectification [34]. However,
THz time-domain spectroscopy (THz-TDS) requires broadband THz radiation to
access extended absorption spectra, allowing the identification of molecules
or precise determination of the physics of broadband processes such as
relaxation and conductivity. This need can be fulfilled by ionizing a gas with
an ultrashort multi-color laser pulse. This method generates strong ($\sim
0.1-1$ GV/m) and ultra-broadband (exceeding 100 THz) terahertz fields with
significant conversion efficiencies ($>10^{-4}$) [35, 36]. Another benefit of
this technique is the non-destructive nature of the mechanisms at play, e.g.,
in air plasmas. The laser driver usually comprises only two harmonics.
However, recent theoretical and experimental results showed that adding more
harmonics to the optical beam can dramatically increase the THz yield [37, 38,
39, 40].

When a temporally asymmetric femtosecond light pulse composed of a fundamental
component and its second harmonic reaches the $\sim 10^{14}$ W/cm2 intensity
range, sufficient to photoionize a gas (here, ambient air) in the tunnel
regime [41] an underdense plasma is created. The electrons freed by this
process are then accelerated by the asymmetric laser field and induce a
photocurrent that contains a low-frequency component. This quasi-DC
contribution causes the plasma to radiate in the THz domain. Although this
process is usually performed in the so-called focused regime, where the plasma
spot is formed a few tens of cm after a focusing element, recent works have
confirmed its efficiency also in the long-range filamentation regime [42] that
results from a sequence of local equilibria between Kerr self-focusing and
plasma defocusing. This regime can lead to meter-long filaments [43] and could
be advantageous for THz-TDS of explosives for safety reasons [44, 45].

One specificity of THz wave generation by two-color air plasmas is the
typically conical shape of the THz radiation [46]. In their seminal work, You
et al. [47, 48] advocated that conical emission is due to interferences
between the THz waves produced at different places along a generated
cylindrical homogeneous plasma, as their polarity changes with the phase angle
between the two laser harmonics caused by optical and plasma dispersion.
Importantly, the key parameters driving THz emission in this scenario are the
plasma length $L$, over which THz emission takes place, and the dephasing
length $l_{d}$, over which the two pump harmonics get out of phase. The latter
length depends on the plasma density $N_{e}$. Other models describing the
angular THz power spectrum have been elaborated in the literature; most of
them also assume a homogeneous plasma cylinder [49, 50]. Although these
assumptions are generally met, in practice, optical beam distortions are
always present. They can affect the emitter’s stability and symmetry, thereby
degrading the generated THz radiation. An alternative theory [51] tried to
account for the role of plasma opacity, relevant for THz frequencies lower
than the local plasma frequency, as well as oblique diffraction of THz waves
through a transversally varying plasma channel and the experimental detector’s
spectral sensitivity. Although showing reasonable agreements between semi-
analytical results and experimental measurements, this work does not include
enough information about key parameters in THz emission (selected THz
frequencies, plasma profiles, or the experimental detector frequency
responses) to provide an accurate and convincing model for conical emission.
In particular, this reference claims the emission is always conical at small
THz frequencies. At the same time, one of the remaining problems to date is
the seemingly generic coexistence of on-axis emission at low frequencies
($\lesssim 5$ THz) with a donut-shaped emission at large THz frequencies ($>5$
THz) in focused propagation [46, 52].

This paper addresses this open question and is organized into three parts.
Section II introduces the underlying physics of THz generation, mainly
dominated by photocurrents. It also recalls the so-called ”Local Current” (LC)
model, extended to include the instantaneous Kerr response. In section III, we
establish a revised (”augmented”) THz emission model that keeps the key
ingredients of the theory developed by You et al. [47] but accounts for (i)
the spectral distributions characterizing the propagation regimes through
three baseline configurations marked by different interplays between
photocurrents and Kerr effect, (ii) transverse variations in the emitting
source induced by laser spatial inhomogeneities or a simple linear focusing
phase, and (iii) optical envelope profiles and changes in the laser
polarization state, e.g., from linear to circular. Results obtained from this
augmented model are compared with the THz spatio-spectral distributions
provided by vectorial full space and time-resolved unidirectional pulse
propagation (UPPE [53]) simulations. Section IV is dedicated to the influence
of optical aberrations on laser propagation, plasma geometry, and THz
radiation. We model optical aberrations using Zernike polynomials that affect
the input pulse phase [54]. We then compare the resulting angular THz spectra
evaluated by our semi-analytical model with UPPE simulation results for the
most common optical aberrations in laser experiments: sphericity, astigmatism,
and coma. Section V concludes this work.

##  II Physical mechanism for THz generation and laser propagation model

We study the THz emission by two-color femtosecond (fs) light pulses that
trigger air plasmas. The laser field can be either linearly polarized along
the $x$-axis when the two harmonics are oriented along the same transverse
direction (LP-P for ”Linear Polarization - Parallel”) or circularly polarized
when the colors are circularly polarized and both rotate in the same direction
(CP-S for ”Circular Polarization - Same”). Previous work [55, 56] reported
that CP-S laser pulses can increase by a factor $\sim 4-6$ the laser-to-THz
conversion efficiency compared to their LP-P counterpart. The temporal
profiles of the input pulse at $z=0$ are chosen as

$\displaystyle\vec{E}_{\rm LP-P}$
$\displaystyle=\left[\mathcal{E}_{1}\cos(\omega_{0}t)+\mathcal{E}_{2}\cos(2\omega_{0}t+\phi)\right]\mathrm{\mathbf{e}}_{x}$
(1) $\displaystyle\vec{E}_{\rm CP-S}$
$\displaystyle=\frac{1}{\sqrt{2}}\binom{\mathcal{E}_{1}\cos(\omega_{0}t)+\mathcal{E}_{2}\cos(2\omega_{0}t+\phi)}{\mathcal{E}_{1}\sin(\omega_{0}t)+\mathcal{E}_{2}\sin(2\omega_{0}t+\phi)}$
(2)

where $\omega_{0}$ is the laser fundamental frequency, $\phi$ the phase angle
between the two colors and $\mathcal{E}_{1,2}$ their envelopes being initially
Gaussian in time and space:

$\mathcal{E}_{1,2}(\mathrm{\mathbf{r}}_{\perp},t)=\sqrt{r_{1,2}}E_{0}\exp\\!\left[-\frac{t^{2}}{t_{1,2}^{2}}-\frac{r_{\perp}^{2}}{w_{1,2}^{2}}\right].$
(3)

Here $E_{0}$ is the laser amplitude, $t_{1,2}$ and $w_{1,2}$ are the $1/e^{2}$
duration and width of the two harmonics, and $k_{1,2}$ their corresponding
wavenumbers. The peak intensity ratio of fundamental over second harmonic is
determined by $r_{2}=r$ and $r_{1}=\sqrt{1-r}$, with $0\leq r\leq 1$. A
focusing lens with focal length $f$ is taken into account by adding a phase to
the input pulse in the temporal Fourier domain,

$\widehat{\mathrm{\mathbf{E}}}(r_{\perp},\omega)\quad\rightarrow\quad\widehat{\mathrm{\mathbf{E}}}(r_{\perp},\omega)\exp\\!\left(-i\frac{\omega}{c}\frac{r_{\perp}^{2}}{2f}\right).$
(4)

Later in Sect. IV, the transverse spatial dependency of this phase term
$\propto r_{\perp}^{2}/2f$ describing a perfect focusing lens will be
generalized to an $(x,y)$-dependent function $Z_{\rm
OA}(\mathrm{\mathbf{r}}_{\perp})$, accounting for optical aberrations.

Photoionization caused by the laser pulse with intensity $\sim 10^{13-15}$
W/cm2 takes place in the tunneling regime [41, 57, 58], leading to a stepwise
increase in the free electron density described by the rate equation

$\partial_{t}N_{e}=W(E)\left(N_{a}-N_{e}\right),$ (5)

where $W(E)$ is the photo-ionization rate, and $N_{a}$ is the initial neutral
density (here written for one species for simplicity). The free electrons are
then accelerated by the laser field and generate photocurrents governed for
non-relativistic plasmas by

$\left(\partial_{t}+\nu_{c}\right)\vec{J}_{e}=\frac{e^{2}}{m_{e}}N_{e}\vec{E},$
(6)

where $\nu_{c}$ denotes an electron-neutral collision rate taken as
$\nu_{c}=2.86$ THz. Although low-frequency radiation is mainly generated from
photocurrents, for weaker laser intensities ($<10^{13}$ W/cm2), an additional
contribution to the overall THz yield may be provided by the Kerr
polarization, that is, the third-order nonlinear response of the bound
electrons [59, 46, 60] reading

$\mathrm{\mathbf{P}}_{K}(\mathrm{\mathbf{r}},t)=\varepsilon_{0}\chi^{(3)}{E}^{2}(\mathrm{\mathbf{r}},t)\mathrm{\mathbf{E}}(\mathrm{\mathbf{r}},t),$
(7)

with $\chi^{(3)}$ being the third-order electric susceptibility [61] assumed
to be non-dispersive, i.e.,
$\chi^{(3)}\equiv\chi^{(3)}(-\omega_{0},\omega_{0},\omega_{0})$. A second
well-known action of the Kerr response is to transversally compress the laser
beam and maintain it in a self-focused state through optical self-focusing and
plasma defocusing sequences as its peak power remains higher than the critical
power for self-focusing [43],

$P_{\rm cr}\approx\frac{3.72\lambda_{0}^{2}}{8\pi n_{0}n_{2}},$ (8)

with $n_{0}=n_{\rm opt}(\omega_{0})$, $n_{\rm opt}(\omega)$ denoting the
linear optical refractive index of air taken from Ref. [62], and
$n_{2}=3\chi^{(3)}/(4n_{0}^{2}c\varepsilon_{0})$ is the nonlinear refractive
index.

During propagation in air, the laser pulse can also excite rotational
transitions, mostly in $N_{2}$ molecules, leading to stimulated Raman
rotational scattering (SRRS) with polarization [63]

$\mathrm{\mathbf{P}}_{R}(t)=\frac{3}{2}x_{K}\varepsilon_{0}\chi^{(3)}\int_{0}^{+\infty}\mathrm{\mathbf{G}}(\mathrm{\mathbf{E}},t,\tau)R(\tau)d\tau,$
(9)

where $x_{K}$ is the ratio of SRRS to the total third-order response.
$\mathrm{\mathbf{G}}(\mathrm{\mathbf{E}},t,\tau)$ is a function describing the
third-order interaction:

$\mathrm{\mathbf{G}}(\mathrm{\mathbf{E}},t,\tau)=E^{2}(t-\tau)\mathrm{\mathbf{E}}(t)-\frac{2}{3}\binom{E^{2}_{y}(t-\tau)E_{x}(t)}{E^{2}_{x}(t-\tau)E_{y}(t)},$
(10)

and
$R(\tau)=\tau_{1}\left(\tau_{1}^{-2}+\tau_{2}^{-2}\right)\sin(\tau/\tau_{1})e^{-\tau/\tau_{2}}$
is the delayed SRRS response with inverse rotational frequency
$\tau_{1}\approx 62.5$ fs, and $\tau_{2}\approx 77$ fs refers to the
characteristic dipole dephasing time [64, 65]. The ratio $x_{K}$ is estimated
to $\approx 80$ % according to Zahedpour et al. [66]. Because both the
instantaneous Kerr response and the delayed SRRS response contribute to the
third-order nonlinear polarization, a factor $(1-x_{K})$ must be applied in
Eq. (7) when accounting for SRRS. As shown in Ref. [63], the SRRS contribution
to the overall THz yield is negligible compared to the instantaneous Kerr
response and, above all, to the photocurrents.

In the following, we shall perform a systematic comparison between semi-
analytically computed angular spectra and far-field THz spectra calculated
using the unidirectional pulse propagation equation (UPPE) [67, 53],

$\frac{\partial\overline{\mathrm{\mathbf{E}}}}{\partial
z}=i\sqrt{k^{2}(\omega)-k_{\perp}^{2}}\ \overline{\mathrm{\mathbf{E}}}\\\
+\frac{i\mu_{0}\omega^{2}}{2k(\omega)}\left[\frac{i}{\omega}\left(\overline{\mathrm{\mathbf{J}}}_{e}+\overline{\mathrm{\mathbf{J}}}_{\rm
loss}\right)+\overline{\mathrm{\mathbf{P}}}_{K}+\overline{\mathrm{\mathbf{P}}}_{R}\right].$
(11)

Here, the notation

$\overline{f}(\mathrm{\mathbf{k}}_{\perp},\omega)=\iiint_{-\infty}^{+\infty}f(\mathrm{\mathbf{r}}_{\perp},t)e^{-i\mathrm{\mathbf{k}}_{\perp}\cdot\vec{r}_{\perp}+i\omega
t}d\mathrm{\mathbf{r}}_{\perp}dt$

refers to the Fourier transform of any function
$f(\mathrm{\mathbf{r}}_{\perp},t)$ in time and transverse variables, and
$\mathrm{\mathbf{J}}_{\rm loss}$ is the loss current density associated with
photoionization. In our analysis and simulations, the propagation medium is
ambient air, composed of 20 % O2 and 80 % N2 with ionization energies
$U_{i}=12.1$ eV [68] and $U_{i}=15.6$ eV [68], respectively. The initial
neutral density is $N_{a}=2.7\times 10^{19}$ cm-3 and photoionization is
modeled using the ionization rate of Popov, Perelomov and Terent’ev (PPT) [58]
with effective charge numbers $Z_{O_{2}}^{*}=0.53$ and $Z_{N_{2}}^{*}=0.9$
according to Talebpour et al. [69].

Two representative $1/e^{2}$ beam widths will be considered, namely
$w_{1,2}=250\,\mu$m for spatially narrow beams and $w_{1,2}=2.5$ mm for broad
beams. To simulate the latter in the long-range filamentation regime and save
computational cost, a pre-processor was used, which integrates the linear part
of Eq. (11) before reaching the electric field threshold of $\sim 1$ GV/m
($\equiv 10^{11}$ W/cm2) above which the (Kerr) nonlinearity kicks in. Because
the pre-processor solves a linear propagation equation, the solution can be
obtained analytically:

$\overline{\mathrm{\mathbf{E}}}(\omega,z+\Delta
z)=\overline{\mathrm{\mathbf{E}}}(\omega,z)e^{i\Delta
z\sqrt{k(\omega)^{2}-k_{\perp}^{2}}}.$ (12)

Further Fourier interpolation was then carried out to refine the spatial
resolution of the resulting data, which were used as initial conditions for
the UPPE simulations to properly compute the subsequent nonlinear propagation.
The spatial resolution used in the UPPE simulations was $\Delta x=\Delta
y=4\,\mu$m for narrow beams and $\Delta x=\Delta y=5-7\,\mu$m for broad beams
and up to $15\,\mu$m for the filamentation of beam profiles distorted by coma.
The temporal resolution was $\Delta t=100$ as, and a (co-moving) time window
of $1.5$ ps guaranteed a spectral resolution of $\Delta\nu=0.66$ THz.

Besides performing costly UPPE simulations, we also integrated the ”Local
Current” (LC) model [70, 49, 71], which neglects propagation effects and
approximates the electric field emitted by the plasma and Kerr source terms as

$\widehat{\mathrm{\mathbf{E}}}(\omega)\propto\left(-i\omega\widehat{\mathrm{\mathbf{J}}}_{e}(\omega)-\omega^{2}\widehat{\mathrm{\mathbf{P}}}_{K}(\omega)\right),$
(13)

$\widehat{f}(\omega)$ denoting the Fourier transform of $f(t)$.

##  III THz emission model

Analytical results from the LC model provide the dependence on the phase angle
$\phi$ of the THz amplitudes generated by photocurrents and Kerr response [56]
as

$\displaystyle\mathrm{\mathbf{E}}_{\rm THz,LP-P}^{\rm J}$
$\displaystyle\propto\sin\phi\ \mathrm{\mathbf{e_{x}}},$
$\displaystyle\mathrm{\mathbf{E}}_{\rm THz,LP-P}^{\rm K}$
$\displaystyle\propto\cos\phi\ \mathrm{\mathbf{e_{x}}},$ (14)
$\displaystyle\mathrm{\mathbf{E}}_{\rm THz,CP-S}^{\rm J}$
$\displaystyle\propto\binom{\sin\phi}{\cos\phi},$
$\displaystyle\mathrm{\mathbf{E}}_{\rm THz,CP-S}^{\rm K}$
$\displaystyle\propto\binom{\cos\phi}{-\sin\phi}.$ (15)

Here, $E_{\rm THz}^{\rm J}$ and $E_{\rm THz}^{\rm K}$ are the THz electric
fields produced by the photocurrents (index J) and Kerr response (resp. K) of
LP-P and CP-S pump pulses, respectively. In practice, plasma channels extend
over millimeter to centimeter lengths in the focused regime and up to several
meters in the (collimated or loosely focused) filamentation regime. Over
sufficiently long propagation ranges, the optical and plasma dispersion induce
an important variation of the phase $\phi$ along the plasma channel, i.e., for
an axially homogeneous plasma:

$\phi(z)=\phi_{0}+\pi\frac{z}{l_{d}},$ (16)

where $\phi_{0}$ is the phase angle taken at the plasma onset distance. The
”dephasing length” $l_{d}$ represents the propagation distance over which
$\phi$ changes by $\pi$ due to dispersion, that is,

$l_{d}=\frac{\lambda_{0}}{4\left[n(2\omega_{0})-n(\omega_{0})\right]}$ (17)

with $\lambda_{0}=2\pi c/\omega_{0}$. Here, $n(\omega)$ is the medium
refractive index that accounts for the contributions of Kerr and plasma

$n^{2}(\omega)=n^{2}_{\rm
opt}(\omega)+2n_{0}n_{2}I_{0}-\frac{\omega_{p}^{2}}{\omega^{2}(1+i\nu_{c}/\omega)},$
(18)

where $I_{0}=c\varepsilon_{0}n_{0}E_{0}^{2}/2$ is the incident pulse intensity
and $\omega_{p}=\sqrt{e^{2}N_{e}/\varepsilon_{0}m_{e}}$ is the electron plasma
frequency. Here, we recall that $l_{d}$ depends on $N_{e}$ through the plasma
frequency in Eq. (18). Note the difference by a factor of 2 with You et al.’s
definition of the dephasing length, which, we believe, proceeds from a
misprint in the original Ref. [47].

###  III.1 You et al.’s conical emission model

Before dwelling upon our revised model, we find it instructive to recall the
main ingredients of You et al.’s [47, 48] theory that describes the THz
intensity radiated in the far-field by a homogeneous ($N_{e}=$ const.)
cylindrical plasma of length $L$ and radius $a$. Here, the only mechanism
generating THz waves is the photocurrent density excited by an LP-P laser
pulse, and thus, according to Eq. (14), the local THz yield evolves as
$\sin{\phi(z)}$. Derived from the radiated electric field in the far-field
approximation, You et al. obtained the angular THz intensity spectrum (see
Appendix for a detailed derivation):

$|\widehat{E}(\omega_{\rm THz},\Theta)|^{2}\propto|A(\omega_{\rm
THz})|^{2}\left(\frac{J_{1}(\beta)}{\beta}\right)^{2}\\\ \times
L^{2}\left[\kappa_{+}^{2}+\kappa_{-}^{2}-2\kappa_{+}\kappa_{-}\cos(2\phi_{0})\right],$
(19)

where $A(\omega)$ is the spectral amplitude of the THz source, and $J_{1}$ is
the Bessel function of first order. The quantities $\kappa_{\pm}$ and $\beta$
are given by

$\displaystyle\kappa_{\pm}$ $\displaystyle=\mathrm{sinc}\ \alpha_{\pm}$ (20)
$\displaystyle\alpha_{\pm}$ $\displaystyle=\frac{\pi L}{\lambda_{\rm
THz}}\left(1-\cos\Theta\pm\frac{\lambda_{\rm THz}}{2l_{d}}\right)$ (21)
$\displaystyle\beta$ $\displaystyle=\frac{2\pi}{\lambda_{\rm
THz}}a\sin\Theta,$ (22)

with $\mathrm{sinc}(x)=\sin x/x$ and $\lambda_{\rm THz}$ denotes the
wavelength associated with the THz frequency $\omega_{\rm THz}$. The functions
$\kappa_{\pm}$, linked to the phase matching condition for THz generation,
depend on the extent $(L)$ of the axial filament. The quantity $\beta$
accounts for the radial dimensions of the plasma ($\beta\propto a$) entering
the term $J_{1}^{2}(\beta)/\beta^{2}$ that describes the circular diffraction
of the THz waves emitted by the plasma cylinder. For common filaments, the
radius $a$ is limited to a few tens of $\mu$m and $\Theta<10$°, i.e.,
$\beta\ll 1$ in the main part of the spectrum, so this diffraction term has
limited influence. The maximum emission is reached for $\alpha_{-}=0$, thus
providing the theoretical emission angle $\Theta_{p}$:

$\cos\Theta_{p}=1-\frac{\lambda_{\rm THz}}{2l_{d}}.$ (23)

The phase mismatch between the THz waves emitted at different places in the
plasma filament causes destructive interference for the on-axis THz component
for plasmas of length $L>l_{d}$, and the THz radiation propagates as a cone in
space, far from the plasma volume. The emission angles are usually narrow,
i.e., $\Theta$ is in the 2° – 10° range [46, 72, 51, 73].

As summarized in Fig. 1, the angular THz spectrum described by this model
depends mainly on three parameters. The first is the length of the plasma $L$.
When $L<l_{d}$, the THz radiation is maximal on the optical axis ($\Theta=0$),
while for longer plasma lengths ($L>l_{d}$), the on-axis THz amplitude
vanishes and the THz beam becomes a cone with angle $\Theta_{p}$ [see Fig.
1(a)]. The second key parameter is the electron density $N_{e}$. Figure 1(b)
indicates that, at a fixed frequency, $\Theta_{p}$ increases for higher
$N_{e}>2\times 10^{17}$ cm-3 while it remains constant for lower density
levels. The dephasing length $l_{d}$ [Eq. (18)] decreases significantly as the
value of $N_{e}$ increases. For larger density levels, conical emission,
therefore, appears for shorter plasma lengths. The third parameter is the
radiated THz frequency $\nu$. In this regard, Fig. 1(c) highlights that the
emission angle $\Theta_{p}$ [Eq. (23)] strongly decreases with increasing
$\nu$ (white curve), consistent with Eq. (19).

Figure 1:  (a,b,c) Mapping of the angular THz spectral intensity
$|\widehat{E}(\omega_{\rm THz},\Theta)|^{2}$ given by Eq. (19) as a function
of the angle $\Theta$ and (a) the plasma length $L$, (b) the electron density
$N_{e}$, and (c) the frequency $\nu=\omega_{\rm THz}/2\pi$. The maps are
normalized to their respective maximum values, and we set $A(\omega)=1$ for
simplicity (as in Ref. [47]). The white curves report the theoretical conical
emission angle $\Theta_{p}$ predicted by Eq. (23).

In addition to the above well-known properties, additional information can be
extracted from the function describing the phase-matching condition in Eq.
(19), namely, $P(L,\Theta,\lambda_{\rm THz})\equiv
L^{2}[\kappa_{+}^{2}+\kappa_{-}^{2}-2\kappa_{+}\kappa_{-}\cos(2\phi_{0})]$.
First, on-axis components associated with the interference function
$P(L,\Theta=0,\lambda_{\rm THz})=(2l_{d}/\pi)^{2}[1-\cos{(\pi
L/l_{d})}][1-\cos{(2\phi_{0})]}$ develop as a $2l_{d}$-periodic function of
the plasma length $L$ and reach maxima with constant amplitude at odd-numbered
multiples of $l_{d}$, provided that the initial phase angle $\phi_{0}$ differs
from 0 [$\pi$]. Second, phase-matched off-axis components parabolically
increase like $P(L,\Theta=\Theta_{p},\lambda_{\rm THz})\simeq
L^{2}[1+\mathrm{sinc}^{2}(\pi L/l_{d})-2\cos(2\phi_{0})\mathrm{sinc}(\pi
L/l_{d})]$. They exceed the on-axis components from $L\approx l_{d}$ and
increase monotonously to higher amplitude levels for longer plasma lengths.
The value of $\phi_{0}$ does not change this behavior and simply introduces
secondary oscillations. Importantly, these properties result from interference
in the far-field THz intensity spectrum.

In the following, we revisit this THz emission model by making it as general
as possible. This includes (i) treating envelope effects, (ii) allowing for
transverse variations in the plasma density, (iii) accounting for the
characteristic spectral signatures of $A(\omega)$ depending on the nonlinear
propagation regime at play, and (iv) passing from LP-P to CP-S driver pulses.

###  III.2 THz conical emission: An augmented model

The electric field $\mathrm{\mathbf{E}}(\mathrm{\mathbf{r}},t)$ emitted by a
conductive current density $\mathrm{\mathbf{J}}(\mathrm{\mathbf{r}},t)$ in a
neutral plasma volume $V$ reads [74]

${\mathrm{\mathbf{E}}}({\mathrm{\mathbf{r}}},t)=-\frac{\mu_{0}}{4\pi}\int_{V}d^{3}{\mathrm{\mathbf{r}}}^{\prime}\frac{\partial_{t}{\mathrm{\mathbf{J}}}({\mathrm{\mathbf{r}}}^{\prime},t^{\prime})}{|{\mathrm{\mathbf{r}}}-{\mathrm{\mathbf{r}}}^{\prime}|},$
(24)

where $t^{\prime}=t-|{\mathrm{\mathbf{r}}}-{\mathrm{\mathbf{r}}}^{\prime}|/c$
is the retarded time at which the signal the detector receives at time $t$ was
generated. We assume that the refractive index of the medium is $n_{\rm
opt}(\omega_{\rm THz})\equiv 1$ in the THz frequency range. The notation
$\mathrm{\mathbf{r}}^{\prime}$ applies to the local coordinates inside the
plasma volume. Cartesian, cylindrical, and spherical coordinates are written
as ($x$, $y$, $z$), ($r$, $\varphi$, $z$) and ($r$, $\Theta$, $\Phi$),
respectively. We are interested in the far-field distribution, so the distance
between the plasma and the detector is large compared to the dimensions of the
filament, that is, $|\mathrm{\mathbf{r}}|\gg|\mathrm{\mathbf{r^{\prime}}}|$,
which leads to the simplification of the Fourier-transformed field (see
Appendix A):

$\widehat{\mathrm{\mathbf{E}}}(r,\Theta,\Phi,\omega)=\frac{i\omega\mu_{0}}{4\pi
r}\mbox{e}^{i\frac{\omega}{c}r}\int_{-L/2}^{+L/2}dz^{\prime}\mbox{e}^{-i\frac{\omega}{c}\cos\Theta
z^{\prime}}\\\ \times\iint
d^{2}\mathrm{\mathbf{r}}^{\prime}_{\perp}\widehat{{\mathrm{\mathbf{J}}}}({\mathrm{\mathbf{r}}}^{\prime},\omega)\mbox{e}^{-i{\mathrm{\mathbf{k}}}_{\perp}\cdot{\mathrm{\mathbf{r}}}^{\prime}_{\perp}},$
(25)

with
${\mathrm{\mathbf{k}}}_{\perp}=(\omega/c)(\sin\Theta\cos\Phi,\sin\Theta\sin\Phi)$.
The double integral formally applies to the entire transverse plane.
$\widehat{{\mathrm{\mathbf{J}}}}({\mathrm{\mathbf{r}}}^{\prime},\omega)$ is
zero for $|\mathrm{\mathbf{r^{\prime}}}|\rightarrow\infty$, in particular when
considering a finite plasma cylinder of radius $a$.

The second step consists in taking the Kerr polarization [Eq. (7)] into
account through the substitution

$\mathrm{\mathbf{J}}(\mathrm{\mathbf{r}},t)\rightarrow\mathrm{\mathbf{J}}_{e}(\mathrm{\mathbf{r}},t)+\partial_{t}\mathrm{\mathbf{P}}_{K}(\mathrm{\mathbf{r}},t),$
(26)

and calculating the free-electron contribution
$\mathrm{\mathbf{J}}_{e}(\mathrm{\mathbf{r}},t)$ using a toy ionization rate
satisfying $W(E)\propto E^{2}(t)$. This toy ionization rate was already
employed in Refs. [75, 76] and allows a semi-analytical treatment that
provides the major characteristics of multi-color-driven THz pulses, including
their polarization state. Most importantly, it captures the asymmetries in the
temporal location of the ionization peaks introduced by the second pump color,
which are crucial for the photocurrent mechanism. In this respect, approaches
based on Fourier series using more realistic, highly nonlinear ionization
rates but involving only a single pump color [77, 78] cannot account for these
asymmetries.

The third step is to consider a pump field that admits transverse envelopes
and different delays for the two pump colors. Introducing the retarded times
$t_{j=\\{1,2\\}}^{\prime}=t-z/v_{pj}$, with $v_{pj}=c/n(j\omega_{0})$ being
the respective phase velocities, our LP-P two-color laser field takes the
general form

${\mathrm{\mathbf{E}}}_{\rm
L}({\mathrm{\mathbf{r}}},t)={\mathrm{\mathbf{\cal{E}}}}_{1}({\mathrm{\mathbf{r}}},t_{1}^{\prime})\cos{[\omega_{0}t_{1}^{\prime}+\phi_{1}(\mathrm{\mathbf{r}})]}\mathrm{\mathbf{e_{x}}}\\\
+{\mathrm{\mathbf{\cal{E}}}}_{2}({\mathrm{\mathbf{r}}},t_{2}^{\prime})\cos{[2\omega_{0}t_{2}^{\prime}+\phi_{2}(\mathrm{\mathbf{r}})]}\mathrm{\mathbf{e_{x}}},$
(27)

where $\phi_{j}(\mathrm{\mathbf{r}})$ is the spatial phase of the $j$th color.
We assume that the envelopes $\mathcal{E}_{j}$ and their phases $\phi_{j}$
vary little along the filament, that is, they do not explicitly depend on $z$,
and we can write
$\mathcal{E}_{j}({\mathrm{\mathbf{r}}},t_{1}^{\prime})\approx\mathcal{E}_{j}({\mathrm{\mathbf{r}}}_{\perp},t_{1}^{\prime})$
and
$\phi_{j}(\mathrm{\mathbf{r}})\approx\phi_{j}(\mathrm{\mathbf{r}}_{\perp})$.
We note that this is a strong assumption, in particular for the phases, but is
necessary for the upcoming semianalytical treatment. In contrast, our UPPE
simulation results comprehensively account for the complete spatiotemporal
dynamics of the pump field. Furthermore, we assume that temporal and radial
variations of the laser envelopes are separable as

$\mathcal{E}_{j}({\mathrm{\mathbf{r}}}_{\perp},t_{j}^{\prime})=\mathcal{E}_{tr,j}({\mathrm{\mathbf{r}}}_{\perp})\mathcal{E}_{ax,j}(t_{j}^{\prime}),$
(28)

where the maximum of $\mathcal{E}_{tr,j=1,2}$ is normalized to unity and
$\mathcal{E}_{ax,j}$ contains the electric field amplitude expressed in V/m.

Straightforward computations detailed in Appendix A then enable us to
establish the following generic dependencies

$\displaystyle\begin{split}&\quad(\partial_{t}+\nu_{c})\mathrm{\mathbf{J}}_{e}^{\rm
THz}\\\
&\propto\mathcal{E}_{1}^{2}({\mathrm{\mathbf{r}}}_{\perp},t_{1}^{\prime})\mathcal{E}_{2}({\mathrm{\mathbf{r}}}_{\perp},t_{2}^{\prime})\sin{\left[\Delta
kz+\phi(\mathrm{\mathbf{r}}_{\perp})\right]}\mathrm{\mathbf{e_{x}}},\end{split}$
(29) $\displaystyle{\mathrm{\mathbf{P}}}_{K}^{\rm THz}$
$\displaystyle\propto\mathcal{E}_{1}^{2}({\mathrm{\mathbf{r}}}_{\perp},t_{1}^{\prime})\mathcal{E}_{2}({\mathrm{\mathbf{r}}}_{\perp},t_{2}^{\prime})\cos{\left[\Delta
kz+\phi(\mathrm{\mathbf{r}}_{\perp})\right]}\mathrm{\mathbf{e_{x}}},$ (30)

where $\Delta k$ is the phase mismatch in wavenumber between the harmonics and
$\phi(\mathrm{\mathbf{r}}_{\perp})$ their phase angle:

$\displaystyle\Delta
k=\frac{2\omega_{0}}{c}\left[n(\omega_{0})-n(2\omega_{0})\right]=\pi/l_{d},$
(31)
$\displaystyle\phi(\mathrm{\mathbf{r}}_{\perp})\equiv\phi_{2}(\mathrm{\mathbf{r}}_{\perp})-2\phi_{1}(\mathrm{\mathbf{r}}_{\perp}).$
(32)

Further computations lead to

$\left|\widehat{\mathrm{\mathbf{E}}}({\mathrm{\mathbf{r}}},\omega)\right|^{2}=\frac{\mu_{0}^{2}}{16\pi^{2}r^{2}}\left|\mathrm{\mathbf{C}}_{J}(\Theta,\Phi,\omega)+\mathrm{\mathbf{C}}_{K}(\Theta,\Phi,\omega)\right|^{2},$
(33)

where the functions $\mathrm{\mathbf{C}}_{J}(\Theta,\Phi,\omega)$ and
$\mathrm{\mathbf{C}}_{K}(\Theta,\Phi,\omega)$ represent the contributions to
the THz yield from the photocurrents (index $J$) and the Kerr effect (index
$K$), respectively. These general functions are expressed as

$\mathrm{\mathbf{C}}_{J,K}(\Theta,\Phi,\omega)=A_{J,K}(\omega)\mathrm{\mathbf{I}}_{J,K}(\Theta,\Phi,\omega),$
(34)

$A_{J,K}(\omega)$ containing the frequency dependency of the considered source
term and $\mathrm{\mathbf{I}}_{J,K}(\Theta,\Phi,\omega)$ being the associated
angle-dependent longitudinal phase matching function, namely,

$\displaystyle\mathrm{\mathbf{I}}_{J}^{\rm LP-P}(\Theta,\Phi,\omega)$
$\displaystyle=\frac{L}{2}\left[I_{S}^{+}\kappa_{+}-I_{S}^{-}\kappa_{-}\right]\mathrm{\mathbf{e_{x}}},$
(35) $\displaystyle\mathrm{\mathbf{I}}_{K}^{\rm LP-P}(\Theta,\Phi,\omega)$
$\displaystyle=i\frac{L}{2}\left[I_{S}^{+}\kappa_{+}+I_{S}^{-}\kappa_{-}\right]\mathrm{\mathbf{e_{x}}}.$
(36)

Here, $\kappa_{\pm}$ are those given in Eq. (20).

The quantities $I_{S}^{\pm}$ describe the transverse effects of the plasma on
the far-field THz radiation,

$I_{S}^{\pm}(\Theta,\Phi,\omega)=\mathcal{F}_{\perp}\left[J_{\rm
tr}(\mathrm{\mathbf{r}}_{\perp}^{\prime})e^{\pm
i\phi(\mathrm{\mathbf{r}}_{\perp}^{\prime})}\right](\Theta,\Phi,\omega),$ (37)

where notation $\mathcal{F}_{\perp}$ refers to the Fourier transform with
respect to the transverse variables $\mathrm{\mathbf{r}}_{\perp}^{\prime}$ and

$J_{\rm tr}(\mathrm{\mathbf{r_{\perp}}})=\mathcal{E}_{\rm
tr,1}^{2}(\mathrm{\mathbf{r}}_{\perp})\mathcal{E}_{\rm
tr,2}(\mathrm{\mathbf{r}}_{\perp})$ (38)

represents the transverse envelope of the source terms. $I_{S}$ reduces to the
Bessel-like circular diffraction function recalled in Eq. (19) for a
cylindrical homogeneous plasma of radius $a$ (see Appendix B). We also mention
that CP-S pulses can be treated similarly, cf. Appendix A, and one finds

$\displaystyle\mathrm{\mathbf{I}}_{J}^{\rm CP-S}(\Theta,\Phi,\omega)$
$\displaystyle=\frac{L}{2}\binom{I_{S}^{+}\kappa_{+}-I_{S}^{-}\kappa_{-}}{I_{S}^{+}\kappa_{+}+I_{S}^{-}\kappa_{-}},$
(39) $\displaystyle\mathrm{\mathbf{I}}_{K}^{\rm CP-S}(\Theta,\Phi,\omega)$
$\displaystyle=i\frac{L}{2}\binom{I_{S}^{+}\kappa_{+}+I_{S}^{-}\kappa_{-}}{-I_{S}^{+}\kappa_{+}+I_{S}^{-}\kappa_{-}}.$
(40)

By construction, our model does not account for variations in the plasma
density along $z$, and plasma absorption and opacity are also neglected for
reasons addressed in the appendix.

Integration of our model requires some prior knowledge of key parameters of
the emitting filament. First, we need to model the transverse shape of the THz
source $J_{\rm tr}(\mathrm{\mathbf{r}}_{\perp}^{\prime})$ to calculate
$I_{S}^{\pm}$. To this end, we use our linear laser propagation code up to a
point $z=z_{p}$ where the width of the laser beam is small enough to reach
high intensities corresponding to efficient plasma generation. At this
distance, assuming that the envelopes of the two colors remain superimposed in
time, we can extract the quantities $\phi({\mathrm{\mathbf{r}}}_{\perp})$ and
the transverse field profiles, allowing us to reconstruct $J_{\rm
tr}(\mathrm{\mathbf{r_{\perp}}})$ at this distance and thus evaluate
$I_{S}^{\pm}$. Second, as addressed in the appendix for illustrative purposes,
Gaussian temporal envelopes impose a radiated Gaussian spectrum. This is
unrealistic because the temporal pump pulse profiles strongly evolve along
propagation due to the Kerr and plasma response and typical THz spectra
obtained from experiments or simulations look very different. We therefore
replace the profiles $A_{J,K}(\omega=2\pi\nu)$ in Eq. (34) by typical spectral
shapes delivered by the respective source term in the specific propagation
regime under consideration. In what follows, we use

$\displaystyle A_{J}(\nu)$ $\displaystyle\propto\left(1-\mbox{e}^{-\nu/f_{\rm
max}}\right)\mbox{e}^{-(\nu/f_{0})^{2}}$ (41) $\displaystyle A_{K}(\nu)$
$\displaystyle\propto\left(1-\mbox{e}^{-\nu/f_{\rm
max}}\right)\mbox{e}^{-(\nu/f_{0})^{4}}$ (42)

yielding the synthetic profiles shown in Fig. 2 for typical THz spectra
produced by photocurrents alone (blue curve), or Kerr effect discarding SRRS
(red curve) in either short-range focused (solid) and long-range filamentation
(dashed curves) regimes. They have been inspired by Figs. 3(c,d) of Ref. [60]
and Fig. 2 of Ref. [59] and are comparable to our transversely integrated UPPE
spectra (not shown). Typically, plasma-driven THz spectra go up to $\sim
70-100$ THz, whereas bound electron responses only produce THz waves below
$\sim 40$ THz in the focused regime. Their contribution may be extended to 80
THz in the long-range filamentation regime. In this self-guiding regime,
nonlinearities develop a supercontinuum through which the THz spectrum can
merge with the broadened optical pump spectrum [79]. Three typical interaction
regimes will be explored in the following three subsections: purely plasma-
driven emission in focused geometry, plasma regime interplaying with a strong
Kerr response, and long-range filamentation regime.

Figure 2:  Characteristic shapes of THz spectra due to photocurrents (blue
curves) and Kerr effect (red curve) in focused (solid curves) and
filamentation (dotted curves) regimes inspired from Refs. [60] and [59]. The
fraction of Kerr response can increase up to $10\%$ in amplitude in the
filamentation regime or up to $77\%$ in the strong Kerr regime. Parameters in
Eq. (41) are $f_{\rm max}=5$ THz and $f_{0}=60$ THz in focused regime and
$f_{\rm max}=5$ THz and $f_{0}=120$ THz in filamentation regime. Parameters in
Eq. (42) are $f_{\rm max}=10$ THz and $f_{0}=30$ THz in focused regime and
$f_{\rm max}=5$ THz and $f_{0}=60$ THz in filamentation regime.

Lastly, the ”conical” functions $\kappa_{\pm}$ are obtained by estimating the
plasma length $L$ and the on-axis electron density $N_{e}$, which provides the
dephasing length $l_{d}$. In estimating $N_{e}$ one has to remember that the
THz waves are produced before the air plasma is fully established [72], which
means that the peak electron density (henceforth denoted by $N_{e}^{\rm max}$)
is considerably higher than the effective electron density $N_{e}^{\rm eff}$
that should be used in the model. Illustrative examples obtained from UPPE
simulations are given in Fig. 3 detailing the maps of THz fields emerging from
the interaction between the two-color laser field and the generated free
electron density created in focused or extended propagation regimes. We can
observe that in the focused regime [Figs. 3(a,b)] the THz field is always
created at the front of the plasma, which justifies that the effective density
level $N_{e}^{\rm eff}$ associated with the emitted THz wave is always lower
than the peak plasma density. Another important observation is that the
effective plasma length $L$ corresponding to the longitudinal distance along
which the THz field is emitted remains approximately on the order of the
distance along which the relative phase between the two laser harmonics
undergoes a $\pi$ phase shift, that is, $L\gtrsim l_{d}$. By contrast, in the
filamentation regime [Figs. 3(c,d)], we can see that although the THz pulse
still emerges at the plasma front along Kerr self-focusing - plasma defocusing
sequences, the plasma length generally satisfies $L\geq 2\,l_{d}$.

Figure 3:  Left-hand side column: on-axis THz fields (blue areas) generated by
two-color laser fields (green curves) triggering electron plasma densities
(red areas) in the $(z,t)$ plane in (a) focused plasma regime and (c) meter-
range filamentation regime. Right-hand side column: (b,d) corresponding on-
axis relative phase between the fundamental and second harmonic
$\phi(z,t)=\phi_{2}(z,t)-2\phi_{1}(z,t)$ at $t=0$ together with the created
THz field and plasma density. The simulation parameters are detailed in the
respective subsections below.

###  III.3 Focused regime without Kerr response

In this subsection, we compare the results of our augmented model with You et
al.’s predictions and confront them with UPPE simulation results. We simulate
the nonlinear propagation of a two-color LP-P Gaussian laser pulse with 0.8 µm
fundamental wavelength, 500 $\mu$J energy, 34 fs duration, and 250 µm initial
radius. 35% of the laser energy is allocated to the second harmonic with zero
relative phase. First, the pulse is focused in air by a lens of 2.6 cm focal
length. The ionization is modeled with the PPT rate [58]
($Z_{O_{2}}^{*}=0.53,Z_{N_{2}}^{*}=0.9$). To make photocurrents the only
source of THz radiation in the simulation, we neglect the Kerr response and
SRRS by setting the nonlinear index of air $n_{2}$ to zero.

Figures 4(a,b) plot the evolutions of the maximum electron density along $z$
and the THz energy contained in the simulation window, with the associated
effective density $N_{e}^{\rm eff}$ at the time instants when the THz field is
maximum. The peak density $N_{e}^{\rm max}$ is rather constant $\approx
10^{18}$ cm-3 (resp. $N_{e}^{\rm eff}\approx 2\times 10^{17}$ cm-3) during the
increase in THz energy, taking place between 1.5 and 2.25 cm. Figures 4(c,d,e)
illustrate the angular THz intensity spectra calculated by Fourier transform
in time and transverse space at $z=1.75$ cm (plasma onset), $z=2$ cm (middle)
and $z=2.5$ cm (end of the emitting zone), respectively. Upon laser
propagation, the initial on-axis THz radiation progressively transforms into a
cone, starting from the upper part of the spectrum. This spectral pattern is
often met in the literature [46, 72]. The generic coexistence of on-axis
emission at low frequencies with a donut-shaped emission at large THz
frequencies in such focused propagation regimes [46, 52] is associated with an
effective plasma length $L$ that remains close to the dephasing length
$l_{d}$, as justified by Fig. 3(a,b). The spectrum plotted in Fig. 4(e)
contains all the THz energy produced during the simulation. It presents a
notable contribution on the optical axis for $\nu<10-15$ THz, and conical
emission at $\Theta\approx 3$° for $\nu>15$ THz.

Figure 4:  Results of UPPE simulations for laser-air parameters described in
the text. (a,b) Evolutions along $z$ of (a) maximum electron density and (b)
THz energy ($\nu<90$ THz) in the simulation box (blue curve, left axis) and
plasma density at the instants when the THz field is maximum (red curve, right
axis). (c,d,e) Angular THz spectra in intensity calculated by transverse
Fourier transform of the THz field at three different propagation distances.

We now compare Figs. 4(c,d,e) to the angular intensity spectra predicted by
You et al.’s model and by our augmented model in Fig. 5. We assume that the
only source of THz waves is photocurrents [$A_{\rm K}(\omega)=0$]. Because we
operate in the focused regime, the spectrum $A_{J}(\omega)$ is taken as the
solid blue curve of Fig. 2. When integrating You et al.’s model, the plasma
radius is chosen as $a=30$ µm and the effective electron density is
$N_{e}^{\rm eff}=2\times 10^{17}$ cm-3, yielding the dephasing length
$l_{d}=4$ mm. The maximum emitter length is determined by the growth of THz
energy, thus $L=0.75$ cm as justified above. When integrating our augmented
model, we employ the same input laser-air parameters in our linear propagator
as those used in the UPPE simulation. The relative phase $\phi_{0}=\pi/4$
between the two colors is chosen at the input ($z=0$). We propagate the laser
pulse linearly up to $z_{p}=1.85$ cm ($f=2.6$ cm), where the laser fluency
takes a shape similar to that calculated by UPPE, i.e., it reaches a
transverse Gaussian profile being 75 µm wide (not shown). The corresponding
laser amplitude was checked to provide the same peak plasma density value
($N_{e}^{\rm max}\approx 10^{18}$ cm-3) as in Fig. 4(a), and a 50 µm wide
transverse plasma density profile.

Figure 5:  THz angular spectra evaluated by (a,b,c) You et al.’s model and
(d,e,f) our augmented model at (a,d) $z=1.75$ cm $(L=2.5$ mm), (b,e) $z=2$ cm
$(L=5$ mm) and (c,f) $z=2.5$ cm $(L=10$ mm). The laser-gas parameters are
specified in the text; $z_{p}=1.85$ cm.

Figure 5 illustrates the results given by You et al.’s model (a,b,c) and by
our augmented model (d,e,f) at $z=1.75$ cm, $2$ cm, and $2.5$ cm. At these
distances, the effective plasma length $L$ is 2.5 mm, 5 mm and 10 mm,
respectively. We notice that the spectra predicted by the two models in the
cases where $L<l_{d}$ [Figs. 5(a,d)] and $L>l_{d}$ [Figs. 5(c,f)] are very
similar, confirming that You and al.’s model is well suited for radially
symmetric plasma channels. However, for $L=5\leavevmode\nobreak\
\textrm{mm}\approx l_{d}$, the strong on-axis contribution predicted by the
UPPE simulation is better restored by the augmented model [Fig. 5(e) to be
compared with Fig. 5(b)]. This difference comes from the phase modulation
caused by the phase term $e^{\pm i\phi(\mathrm{\mathbf{r}}_{\perp}^{\prime})}$
in Eq. (37), which accounts for additional interferences between the THz waves
formed in the overall plasma volume. This effect is absent in You et al.’s
model. Note that for $L=10\leavevmode\nobreak\ \mathrm{mm}\approx 3\,l_{d}$,
revivals of on-axis emission occur in this spectral dynamics – in agreement
with the interference function $P(L,0,\lambda_{\rm THz})$ discussed in the
previous subsection – and they are more pronounced in the augmented model.

###  III.4 Focused regime with strong Kerr response

We now examine what happens when a relevant Kerr contribution affects the
former focused propagation. A second UPPE simulation was carried out with a
strong instantaneous Kerr polarization by setting $x_{K}=0$ and
$n_{2}=3.8\times 10^{-19}$ cm2/W. This choice of $x_{K}$ implies that SRRS is
ignored, which maximizes self-focusing and related four-wave mixing.

Figure 6(a) shows that the maximum plasma density value, $N_{e}^{\rm
max}\approx 10^{19}$ cm-3, is about 10 times higher than in the previous
configuration, which is caused by the strong Kerr self-focusing effect.
Noticeably, this increase in density does not induce a growth in the laser-to-
THz conversion efficiency [compare Figs. 4(b) and 6(b)]. The THz energy in the
simulation box increases between 1.5 and 2.5 cm, hence $L=10$ mm. Although the
density strongly varies with $z$, we select $N_{e}^{\rm eff}\approx 3\times
10^{17}$ cm-3 [consistently with the red curve in Fig. 6(b)], which
corresponds to a dephasing length $l_{d}=2.8$ mm. Figures 6(c,d,e) illustrate
the angular THz spectra calculated by the UPPE code at $z=1.75$ cm, $2$ cm,
and $2.5$ cm. Again, the corresponding plasma lengths $L$ are $2.5$ mm, $5$ mm
and $10$ mm. These spectra present strong on-axis components whose amplitudes
increase along propagation. On-axis emission occurs in the frequency range
$10-25$ THz, which we attribute to the strong Kerr self-focusing that affects
the pump pulse propagation dynamics [60]. Photocurrent-induced off-axis
emission ($\Theta=3-4$°) fully develops for high frequencies $\nu>20$ THz in
Figs. 6(d,e), where the criterion $L>l_{d}$ is fulfilled.

Figure 6:  UPPE results displaying the same quantities as in Fig. 4 for a
medium with strong Kerr response ($x_{K}=0$, $n_{2}=3.8\times 10^{-19}$
cm2/W).

Now, let us evaluate the THz radiation using our semi-analytical model. The
functions $A_{J}(\omega)$ and $A_{\rm K}(\omega)$ are taken as the solid blue
and red curves of Fig. 2, respectively, with the Kerr maximum spectral
amplitude being 77% of the photocurrent one. Here, we linearly propagated the
laser pulse up to $z_{p}=2$ cm (instead of $z_{p}=1.85$ cm) because stronger
Kerr self-focusing fosters a more narrow plasma channel whose equivalent
diameter is attained at longer distances $z$ upon linear propagation. The
amplitude of the optical field reached at $z=z_{p}$ then corresponds to a peak
electron density of $10^{19}$ cm-3, in agreement with Fig. 6(a).

Figure 7:  THz angular spectra evaluated by (a,b) You et al.’s model and
(c,d,e,f) our augmented model at (a,c) $z=1.75$ cm $(L=2.5$ mm), (d) $z=2$ cm
$(L=5$ mm), (b,d) $z=2.3$ cm $(L=8$ mm), and (f) $z=2.5$ cm $(L=10$ mm);
$z_{p}=2$ cm. See text for details.

Figure 7 shows the results of You et al.’s model (a,b) and of the augmented
model (c,d,e,f) for the same plasma lengths as analyzed in Fig. 6, completed
with additional information for $L=8$ mm. The two models predict similar
conical emission angles. However, the angular spectra predicted by our
augmented model include an extra contribution on the optical axis, peaked
around $\nu=15$ THz, which is shaped by the spectral amplitude
$|A_{J}(\omega)-i\omega A_{K}(\omega)|$ where the parabolic Kerr contribution,
here artificially increased, is absent in You et al.’s model. In the full
plasma regime ($z\approx 2.3$ cm), Figs. 7(b,e) both describe a net conical
emission due to photocurrents at angles between 3 and 5°. Yet, only the
augmented model produces off-axis emissions ($\Theta\lesssim 3.5$°) that
develop in the high-frequency range $>25$ THz, similarly to the UPPE
simulations. On-axis revivals are more pronounced for $L=8$ mm $\simeq
3\,l_{d}$ in the augmented model and reach a maximum at $\nu\approx 15$ THz,
due to our parabolic spectral profile. Radial phase variations
$\phi(\mathrm{\mathbf{r}}_{\perp})$ are also expected to contribute on-axis,
following Figs. 5(e,f) of the previous subsection. The difference with UPPE
data may be explained by the strong nonlinearities producing Kerr self-
modulation and plasma defocusing, which render the plasma density in the UPPE
simulation inhomogeneous, hence contradicting the model assumptions. On-axis
components decrease to lower amplitudes in the semi-analytical spectra whereas
they increase along $z$ in the UPPE simulation data. Finally, we underline
that the angular spectra computed from the augmented model become closer to
You et al.’s patterns when the Kerr spectral amplitude decreases below $50\%$
of the photocurrent spectral amplitude. This indicates that such large on-axis
components developed in the UPPE simulations cannot be caused directly by the
Kerr response alone.

###  III.5 Filamentation regime

We conclude this section by addressing the long-range filamentation regime.
Here, the UPPE code was used to simulate a two-color laser pulse ($r=35$%)
with fundamental operating at 0.8 $\mu$m, duration of 34 fs, initial radius of
2.5 mm, and energy of 1.5 mJ, loosely focused in the air (20% O2, 80% N2 \-
$n_{2}=3.8\times 10^{-19}$ cm2/W, $x_{K}=0.8$) by a lens of focal length
$f=2.6$ m. For such a beam, using the linear pre-processor was necessary to
perform a linear propagation of the laser pulse over 2.3 m until reaching
intensities close to the Kerr nonlinearity threshold. The resulting pulse was
then used as the input condition for the UPPE code at $z=0$.

Figure 8:  (a) Maximum electron density and (b) THz energy ($\nu<90$ THz)
contained in the UPPE simulation window as a function of $z$ for a two-color
long-range filament. (c) Angular THz spectrum calculated at $z\simeq 20.4$ cm.

No multiple filamentation occurs (not shown) as the input peak power remains
rather close to critical $P_{\rm cr}$ [80], i.e., $P_{\rm in}\simeq
2.6\,P_{\rm cr}^{\rm inst}$ where $P_{cr}^{\rm inst}\equiv P_{cr}/(1-x_{k})$
refers to the instantaneous critical power. Figure 8(a) shows that the
extended filament maintains a relatively constant maximum electron density
$\approx 10^{17}$ cm-3 between $z=10$ cm and $z=60$ cm, due to clamping
effect, then turns off afterward. Figure 8(b), plotting the THz energy in the
simulation box along $z$, shows a strong decrease beyond $z=18$ cm, due to the
THz radiation leaving the simulation box. Moreover, THz production becomes
less efficient at such long distances because the temporal overlap between the
two laser harmonics deteriorates through temporal walk-off. Figure 8(c)
describes the THz spectrum calculated at $z\simeq 20.4$ cm. We can evaluate
from Fig. 3(c) and Figs. 8(a) and 8(b) a minimum length $L_{\rm min}\simeq 5$
cm associated with an efficient plasma emitting zone in this configuration.
This length corresponds to the distance that a THz wave generated on-axis
($x=y=0$) at an angle $\Theta\approx 1.5$° must travel to leave the simulation
box with a transverse window extending from $-1.3$ mm to $1.3$ mm.

Setting $L=5$ cm, the effective electron density can be estimated as
$N_{e}^{\rm eff}\approx 10^{16}$ cm-3 from Figs. 8(a,b). Note that for such
density levels, the dephasing length depends little on $N_{e}^{\rm eff}$ so
that $l_{d}=2.6$ cm $<L$. We linearly propagated the same laser pulse over a
distance of $z_{p}=2.6$ m to integrate our augmented model. We chose
$A(\omega)$ and $A_{J,K}(\omega)$ as given by Fig. 2 (dashed curves), where
the Kerr component is about $10\%$ that of photocurrents in intensity. Figures
9(a,b) plot the angular spectra calculated from You et al.’s model and ours,
respectively. These spectra are almost identical: They present a conical
emission for $\Theta=1-2$° with a flat frequency profile, which agrees with
the UPPE simulation. The residual central emission $(<10\%$ in spectral power)
is consistent with the expected on-axis contribution for $L<2\,l_{d}$.

Figure 9:  THz angular spectra estimated using (a) You et al.’s model and (b)
our augmented model for the filamentation setup. Laser-medium parameters are
specified in the text.

##  IV Influence of optical aberrations

This section aims to characterize the impact of common optical aberrations
(OAs) on the nonlinear propagation of two-color pulses producing THz
radiation. OAs modify the propagation of a laser beam through distortion of
its transverse phase, thus of its optical wavefront [81]. They tend to alter
the position and shape of the focal spot and, thereby, the geometry of the
plasma channel. OAs are usually classified using two integers: the radial
degree $n$ and the azimuthal degree $m\leq n$ entering the Zernike polynomial
$Z_{n}^{m}(\rho,\Phi)$ [54]:

$Z_{n}^{m}(\rho,\Phi)=\left\\{\begin{array}[]{ll}R_{n}^{m}(\rho)\cos(m\Phi)&\mathrm{if\
}m\geq 0\\\ -R_{n}^{-m}(\rho)\sin(m\Phi)&\mathrm{if\ }m<0\end{array}\right.,$
(43)

where $\rho\equiv\sqrt{2}r_{\perp}/w_{0}$ with $r_{\perp}=\sqrt{x^{2}+y^{2}}$
and $w_{0}$ being the beam width.

$R_{n}^{m}(\rho)=\sum_{k=0}^{\frac{n-m}{2}}\frac{(-1)^{k}(n-k)!}{k!\left(\frac{n+m}{2}-k\right)!\left(\frac{n-m}{2}-k\right)!}\rho^{n-2k}$
(44)

is the radial distribution for an even number of $n-m$, and
$R_{n}^{m}(\rho)\equiv 0$ otherwise.

Below we shall concentrate on the basic configuration of focus ($n=2$, $m=0$)
- consisting of a simple focusing lens effect - and will add to this
elementary lensing three characteristic OAs of interest, namely, vertical
astigmatism ($n=m=2$), vertical coma aberration ($n=3$, $m=-1$) and sphericity
($n=4$, $n=0$). For these four configurations, the polynomials
$Z_{n}^{m}(\rho,\Phi)$ express as:

Focus $\displaystyle Z_{2}^{0}(\rho)=2\rho^{2}-1,$ (45) Astigmatism
$\displaystyle Z_{2}^{2}(\rho,\Phi)=\rho^{2}\cos(2\Phi),$ (46) Coma
$\displaystyle Z_{3}^{-1}(\rho,\Phi)=(3\rho^{3}-2\rho)\sin(\Phi),$ (47)
Sphericity $\displaystyle Z_{4}^{0}(\rho)=6\rho^{4}-6\rho^{2}+1.$ (48)

Each OA is associated with an amplitude coefficient, $c_{n}^{m}$, which has
the dimension of a distance and determines the impact of the respective OA on
the transverse phase of the beam. The total phase distortion is then
conveniently expressed as

$-\frac{\omega}{c}Z_{\rm OA}(r_{\perp},\Phi)=-\frac{\omega}{c}\sum_{n,m\leq
n}c_{n}^{m}Z_{n}^{m}\left(\frac{\sqrt{2}r_{\perp}}{w_{0}},\Phi\right).$ (49)

The polynomial $Z_{\rm OA}(r_{\perp},\Phi)$ henceforth substitutes the
transverse spatial dependency of the phase term $r_{\perp}^{2}/2f$ in Eq. (4).
The sum symbol means that when considering OAs such as sphericity,
astigmatism, or coma, we add the latter to the baseline focus configuration.

Figure 10(a) depicts linear laser propagation in the absence of OA (i.e., with
focus only). The initially 250-µm wide optical pulse is focused in air by a
lens with focal length $f=2.6$ cm. The inset shows the initial spatial phase
of the laser pulse, corresponding to the Zernike polynomial
$Z_{2}^{0}(r_{\perp})$ linked to the focusing effect. This reference case is
compared with sphericity ($c_{4}^{0}=5$ nm), astigmatism ($c_{2}^{2}=300$ nm)
and coma ($c_{3}^{-1}=300$ nm), further added to focus.

Figure 10:  Laser fluency maps in the linear regime (a,b,c,e) in the ($x$,
$z$) plane and (d,f) in the ($y$, $z$) plane for a two-color laser pulse with
800-nm fundamental ($r=35$ %) and initial radius $w_{0}=250$ µm focused in the
air by a lens of focal length $f=2.6$ cm (a) without additional optical
aberration and with added (b) sphericity ($c_{4}^{0}=5$ nm), (c,d) astigmatism
($c_{2}^{2}=300$ nm) and (e,f) coma ($c_{3}^{-1}=300$ nm). Insets show the
corresponding laser spatial phases at $z=0$.

Figures 10(b-f) evidence that laser propagation is strongly influenced by the
OAs. In particular, astigmatism and coma introduce transverse asymmetries. In
the following, we discuss nonlinear UPPE simulation results characterizing the
impact of astigmatism, coma, and sphericity in the focused regime. We begin by
addressing LP-P laser fields, which we will compare to circularly polarized
beams (CP-S) simulations. Finally, we present results in the filamentation
regime for OA-influenced LP-P laser pulses. UPPE far-field spectra are
compared with the results of our augmented model.

###  IV.1 LP-P Polarization

Our two-color laser driver again involves a $0.8$ µm fundamental pulse with
34-fs duration and 250 µm initial radius. The laser energy is 500 µJ, with
$r=35\%$ in the second harmonic; the pump beam is linearly polarized (LP-P)
along $\mathrm{\mathbf{e_{x}}}$ and focused ($f=2.6$ cm equivalent to
$c_{2}^{0}=300$ nm) in air ($n_{2}=3.8\times 10^{-19}$ cm2/W, $x_{k}=0.8$) to
provide us with the reference case (focus). We then carried out three other
simulations including sphericity ($c_{4}^{0}=5$ nm), astigmatism
($c_{2}^{2}=300$ nm) and coma ($c_{3}^{-1}=300$ nm) imposed at $z=0$.

Figure 11:  UPPE-computed maps of the laser fluency for a two-color LP-P pulse
propagating in nonlinear regime (a,b,c,e) in the ($x$, $z$) plane and (d,f)
($y$, $z$) plane (a) in the reference (focus) case, with added (b) sphericity,
(c,d) astigmatism and (e,f) coma aberration. The blue lines show density iso-
contours at $N_{e}^{\rm max}=5\times 10^{16}$ cm-3. Simulation parameters are
specified in the text.

To start with, Fig. 11 shows the evolution of the laser fluencies calculated
with the UPPE code in the nonlinear regime. The propagation dynamics look
similar to those observed in the linear regime (compare with Fig. 10).
Differences (reduction) in the maximum fluency values can be attributed to the
clamping of the maximum laser intensity by the peak plasma density, as
indicated by the blue iso-contour density levels $N_{e}^{\rm max}=5\times
10^{16}$ cm-3. For the reference case [Fig. 11(a)] and in the presence of
sphericity [Fig. 11(b)], the plasma remains axisymmetric. However, adding
astigmatism [Figs. 11(c,d)] and coma [Figs. 11(e,f)] breaks the transverse
symmetry. Note that in the configuration with coma, the plasma is, in fact,
composed of a filament extending along the laser propagation direction,
accompanied by small off-axis zones with non-zero electron density around
$z=1$ cm.

Figure 12(a) describes the axial evolutions of the maximum electron density
for the four considered configurations. Compared to focus (reference case,
blue curve), reaching a peak electron density of $N_{e}^{\rm max}\approx
2-3\times 10^{18}$ cm-3, astigmatism (green curve) and coma aberration
(magenta curve) reach lower density levels by a factor of $\approx 2$ and
$\approx 5$, respectively. Sphericity generates very intense density peaks
($N_{e}\approx 10^{19}$ cm-3) corresponding to the off-axis hot spots visible
in the laser fluence map of Fig. 11(b) at $z\simeq 1.8-2$ m. Sphericity
produces a longer plasma channel than is observed in the focus case, which is
in accordance with the references [82, 83]. Astigmatism leads to a shorter
plasma because this aberration induces an additional focusing effect. Coma
generates a longer and smoother plasma channel because it decreases the
maximum fluency level (and, therefore, the peak electron density). Plasma
generation starts at smaller $z$ distances, in accordance with Fig. 10(e,f).
Figure 12(b), plotting the THz energy ($\nu<90$ THz) accumulated along $z$ in
the simulation window, shows that OA hardly modifies the laser-to-THz
conversion efficiency, except for coma, for which the THz yield drops by a
factor $\approx 5$ due to the lower plasma density.

Figure 12:  Top: (a) Maximum electron density and (b) THz energy ($\nu<90$
THz) along $z$ in the reference (focus) case (blue curve), with sphericity
(red curve), astigmatism (green curve) and coma aberration (magenta curve).
Bottom: Angular THz spectra calculated with the UPPE code (c,d,e,g) in the
($x$, $z$) plane and (f,h) in the ($y$, $z$) plane (c) in the focus case at
$z=2.65$ cm, (d) in the sphericity case at $z=2.24$ cm, (e,f) for astigmatism
at $z=1.93$ cm and (g,h) for coma at $z=1.32$ cm.

Figures 12(c-h) illustrate the associated angular THz spectra computed at the
distances for which the growth of the THz energy plotted in Fig. 12(b) is
maximum. The reference case [Fig. 12(c)] exhibits a conical emission at angles
2.5° $<\Theta<$ 4° and extends up to 60 THz. The THz spectrum in the spherical
case [Fig. 12(d)] appears similar. However, it presents a broadening up to 75
THz, which we attribute to strong nonlinearities associated with the local
density peaks observed in Fig. 12(a). Astigmatism [Figs. 12(e,f)] and coma
[Figs. 12(g,h)] display a reduction of the THz spectral width ($\nu<45$ THz).
In the former case, we observe that the radiated wave appears rather
axisymmetric despite the asymmetric character of astigmatism. It presents a
strong on-axis contribution at low frequencies ($\nu<20-25$ THz) and a conical
emission at $\Theta\approx 4$° in the upper part of the spectrum ($\nu>25$
THz). In the coma case, the azimuthally-varying plasma shape renders the THz
spectrum asymmetrical with different frequency bandwidths for positive and
negative values of $\Theta$. This spectrum comprises on-axis emission for
$\nu<35-40$ THz and conical emission for $\nu>40$ THz at angles between 2° and
3° in absolute value.

For comparison, we integrated our augmented model to reproduce these spectra.
The parameters used for its integration are summarized in table 1. We have
selected our spectral profiles $A_{J}(\omega)$ so that they remain close to
the transversely integrated UPPE THz spectra [see, e.g., solid blue curve in
Fig. 2]. Here, only a plasma driver is accounted for in our test spectral
profile. Introducing OA only shortens the bandwidth for asymmetric (e.g.,
astigmatism) aberrations, as exemplified by Figs. 13(a,b).

| Focus | Sphericity | Astigmatism | Coma
---|---|---|---|---
$z_{p}$ [cm] | 2 | 1.5 | 1.8 | 0.8
$N_{e}^{\rm eff}$ [$10^{18}$ cm-3] | 0.2 | 0.4 | 0.3 | 0.01
$L$ [mm] | 7 | 5 | 3 | 10
$l_{d}$ [mm] | 3.9 | 2.1 | 2.8 | 20
Table 1: Parameters used to integrate our augmented model for conical emission
in the focused regime for LP-P laser pulses. Figure 13:  (a,b) Transversely
integrated THz spectra of Fig. 12 (black curves) and profiles $A_{J}(\omega)$
(blue curves) chosen to integrate our conical emission model (a) in the
sphericity case and (b) the astigmatism case. (c-h) Angular THz spectra along
(c,d,e,g) $k_{y}=0$ and (f,h) $k_{x}=0$ estimated from our augmented conical
emission model for a laser propagation affected by (c) no OA (focus only), (d)
sphericity, (e,f) astigmatism and (g,h) coma aberration. The parameters used
are summarized in table 1.

Figure 13(c) plots the angular THz spectrum estimated by our augmented model
for the focus (reference) case. The distance $z_{p}$ reached by our linear
pre-processor is chosen to match beam size and peak intensity similar to those
developed by the UPPE runs. The predicted emission angles $\sim 2.5-3$° are in
reasonable agreement with those of our UPPE simulations [see Fig. 12(c)]. The
sphericity case, shown in Fig. 13(d), describes a conical emission very
similar to that of the reference case and predicts angles of about $3-4$°. In
agreement with the UPPE calculations, these two figures show that sphericity
has little impact on the produced THz spectrum. For astigmatism [Figs.
13(e,f)], our model reproduces spectral narrowing and relevant on-axis
components. However, it slightly overestimates the spectral asymmetries by
amplifying more azimuthal variations (up to $\approx 9$% in amplitude). These
asymmetries are due to the terms $I_{S}^{\pm}(k_{x},k_{y})$, which ignore
nonlinear propagation effects attached to Kerr self-focusing and plasma
defocusing; the latter locally modify the spatiotemporal distributions of the
laser pulse and the plasma channel along $z$.

Nonetheless, the THz emission angles ($\Theta\approx 4$°) are correctly
restored. In the case of coma [see Figs. 13(g,h)], asymmetries emerging in the
UPPE spectra are partially reproduced. They are also associated with the
asymmetric diffraction caused by the term $J_{\rm
tr}(\mathrm{\mathbf{r}}_{\perp})$. Discrepancies with the UPPE simulation data
can again be attributed to the nonlinear distortion of the plasma channel
geometry omitted in our model. The overall spectral dynamics predicted by our
model are nevertheless consistent with the UPPE simulation results, namely, a
conical emission at angles $\Theta\approx 3$°, with a non-negligible on-axis
contribution as $L\lesssim l_{d}$.

###  IV.2 CP-S Polarization

We performed a similar study for circularly polarized laser pulses using the
same laser-air parameters, except for the pump polarization state. We did not
observe major modifications in the nonlinear laser fluency patterns calculated
by UPPE for the four configurations of interest. We observed that the
geometries of the laser beam and plasma remained similar to Fig. 11, i.e., in
the present focused geometry the laser polarization barely modifies pulse
propagation in the presence of optical aberrations.

Figure 14(a) illustrates the axial evolutions of the maximum electron
densities reached without OA (focus, blue curve), with sphericity (red curve),
astigmatism (green curve) and coma (magenta curve). Again, the peak plasma
dynamics are very similar to those observed in the LP-P configuration. Figure
14(b) shows the related THz energy. We refind the well-known property
following which CP-S laser pulses strongly increase the THz yield compared to
their LP-P counterpart [56]. Without OA (blue curve), the THz energy is
enhanced by a factor $\approx 6$. This property holds with sphericity (red
curve, $\approx\times 4$) and astigmatism (green curve, $\approx\times 3$),
even if these two configurations significantly reduce the laser-to-THz
conversion efficiency compared to the reference case. Coma aberration (magenta
curve) is again harmful for THz radiation due to its detrimental action on the
generated electron density.

Figure 14:  (a,b) Same quantities as in Fig. 12 but for CP-S configuration.
(c-h) Corresponding angular THz spectra calculated by UPPE (c,d,e,g) along
$\mathrm{\mathbf{e_{x}}}$ and (f,h) along $\mathrm{\mathbf{e_{y}}}$ (c) in the
reference case at $z=2.65$ cm, (d) for sphericity at $z=2.24$ cm, (e,f) for
astigmatism at $z=1.93$ cm and (g,h) for coma at $z=1.32$ cm.

Figures 14(c-h) illustrate the associated angular THz spectra. These are close
to the LP-P case, although with slightly narrower frequency bandwidth, which
we attribute to relatively weaker nonlinearities. The action of astigmatism
and coma on the THz distribution [see Figs. 14(e,f,g,h)] is similar to the
LP-P case, although with a less pronounced asymmetry at negative angles in the
coma case. Hence, LP-P and CP-S laser beams provide comparable THz emission
patterns.

Let us now integrate our conical emission model. We assume that the CP-S
polarization state barely modifies the parameters $L$ and $N_{e}^{\rm eff}$ so
that the values presented in Table 1 can be reused here. The $A_{J}(\omega)$
factors were checked to remain similar to those used in the LPP configuration.

Figure 15 plots the results of our model. Angular spectra appear very close to
those calculated in the LP-P case, confirming again that laser polarization
has a limited impact on conical emission. The overall angular spectra for
asymmetric aberrations reasonably agree with their UPPE counterparts in the
predicted spectral extent and production of on-axis components. We observe,
however, a few discrepancies in the on-axis contributions, which are
attributed to a strong dependence of the THz spectrum on the phase
$\phi(\mathrm{\mathbf{r}}_{\perp})$ when $L\approx l_{d}$. We again emphasize
that because the nonlinear responses of the medium alter the laser pulse
during propagation, our model can overestimate some asymmetries in the THz
spectra, e.g., for astigmatism [see Figs. 15(c,d)]. These asymmetries
disappear in the coma case [Figs. 15(e,f)] that amplifies on-axis contribution
and agrees better with the UPPE spectra.

Figure 15:  Same quantities as in Fig. 13(c-h), but for a CP-S laser pulse.

###  IV.3 Influence of optical aberrations in filamentation regime

To explore THz emission in the filamentation regime modified by adaptative
optics, we performed UPPE simulations of LP-P laser pulses without OA (focus),
with sphericity ($c_{4}^{0}=10$ nm), astigmatism ($c_{2}^{2}=150$ nm) and coma
($c_{3}^{-1}=100$ nm). Our 34-fs, two-color laser pulses now have a 2.5 mm
initial radius and 1.5 mJ energy, and they are loosely focused with $f=2.6$ m
(equivalent to $c_{2}^{0}=300$ nm) in air. These simulations required using
the linear pre-processor to propagate the laser pulses up to the threshold
intensity of the first nonlinearities before using them as input data for the
UPPE code. Laser beams without OA, with sphericity, astigmatism, and coma,
were propagated over 2.3 m, 1.7 m, 1.7 m, and 1.5 m, respectively. The value
of the distance $z$ where nonlinear propagation begins was always set to $z=0$
by convention.

Figure 16 shows the laser fluency calculated with UPPE. The reference case
[Fig. 16(a)] displays a plasma zone that extends over several tens of
centimeters. Kerr self-focusing/plasma defocusing cycles induce deviations
from the patterns shown in Fig. 11, particularly for pulses affected by
asymmetric aberrations. The hollow zones observed in the reference case and
with sphericity [Fig. 16(a,b)] at $z\approx 15$ cm and $z\approx 30$ cm,
respectively, indicate a replenishment dynamics [84]. Unlike the focused
regime, despite the presence of astigmatism over a long propagation range, a
rather axisymmetric filament develops over the propagation axis. In addition,
Figs. 16(e,f) for coma contain fewer hot spots compared to the fluency
patterns observed in Figs. 11(e,f). We conclude that in the filamentation
regime, the propagation dynamics is dominated by the nonlinear balance between
Kerr self-focusing and plasma defocusing, and linear effects due to OA are
less important.

Figure 16:  Laser fluency maps calculated in filamentation regime by the UPPE
code (a,b,c,e) in the ($x$, $z$) plane and (d,f) in the ($y$, $z$) plane (a)
in the reference (focus) case, (b) with sphericity, (c,d) astigmatism and
(e,f) coma. The blue iso-contours represent the iso-densities $N_{e}^{\rm
max}=5\times 10^{16}$ cm-3.

Figure 17(a) describes the longitudinal evolutions of $N_{e}^{\rm max}$ in the
reference case (blue curve), with sphericity (red curve), astigmatism (green
curve) and coma (magenta curve). The differences in the plasma onset distances
are mainly due to the selected final distances of the pre-processed linear
propagation runs, which vary from one case to another. Focus and sphericity
induce comparable maximum densities, although the latter reduces the plasma
length. On the other hand, asymmetric aberrations reduce plasma density levels
by a factor of $\approx 2$ and also decrease the length of the plasma zone.
Figure 17(b) plots the produced THz energy ($\nu<90$ THz). We observe that all
aberrations reduce the production of THz yield. OAs are thus generally harmful
for THz generation in the filamentation regime, similar to their impact in the
focused regime.

Figure 17:  Top: (a) Maximum electron density and (b) THz energy ($\nu<90$
THz) along $z$ in the reference case (focus, blue curve), with sphericity (red
curve), astigmatism (green curve) and coma (magenta curve). Bottom: Angular
THz spectra calculated with UPPE (c,d,e,g) along $\mathrm{\mathbf{e_{x}}}$ and
(f,h) along $\mathrm{\mathbf{e_{y}}}$ (c) in the reference case at $z=20.4$
cm, (d) with sphericity at $z=35.6$ cm, (e,f) astigmatism at $z=61.1$ cm and
(g,h) coma at $z=61.1$ cm.

UPPE angular THz spectra are presented in Fig. 17(c-h). In both the reference
case [Fig. 17(c)] and with sphericity [Fig. 17(d)], the THz radiation exhibits
a relatively flat spectrum up to 90 THz, which is typical of the filamentation
regime. Additionally, the emission angle remains nearly constant over a wide
frequency bandwidth. Astigmatism and coma illustrated by Figs. 17(e,f) and
17(g,h) modify the spectrum with an almost zero contribution at low
frequencies ($\nu<15-20$ THz). Like in the focused regime, astigmatism does
not induce any asymmetry in the THz intensity radiated in filamentation
regime. Coma tends to slightly unbalance THz emission between positive and
negative $\Theta$ because of the misalignment of the plasma channel with
respect to the optical axis. The observed THz emission angles remain similar
to those of the reference case in the filamentation regime; although affected
by OA, the effective plasma density is weaker. Moreover, the plasma length in
this regime always satisfies the condition $L\gg l_{d}$, ensuring perfectly
conical emission in the far field. The oblique lobes particularly visible at
low frequencies in the spectra of Figs. 17(c,d) are a numerical artifact due
to a coarser resolution in $\Theta$ at low frequency $\nu$.

| Focus | Sphericity | Astigmatism | Coma
---|---|---|---|---
$z_{p}$ [m] | 2.6 | 2 | 2.2 | 2
$N_{e}^{\rm eff}$ [$10^{18}$ cm-3] | 0.01 | 0.01 | 0.01 | 0.01
$L$ [cm] | 5 | 5 | 5 | 5
$l_{d}$ [cm] | 2 | 2 | 2 | 2
Table 2: Parameters used to integrate our augmented conical emission model in
filamentation regime.

To end with, we integrated our conical emission model using the parameters
summarized in Table 2. The spectral distribution $A_{J}(\omega)$ (dashed blue
curve in Fig. 2) was chosen in each case to approach the broadband UPPE
spectral profiles, which did not change significantly by adding OA. These
include $32\%$ of Kerr contribution in amplitude ($10\%$ in energy ratio)
extending till $\sim$70 THz. The results provided by our model are illustrated
in Figs. 18. As in the UPPE simulations, OAs have no significant impact on the
THz emission angles ($\Theta=1-2$°), which are correctly restored here.
Despite overestimated asymmetries in Figs. 18(c,d), our augmented model
restores the vanishing of the low-frequency region for astigmatism and coma,
together with the main trends of the UPPE spectrum.

Figure 18:  Angular THz spectra (a,b,c,e) along $\mathrm{\mathbf{e_{x}}}$ and
(d,f) along $\mathrm{\mathbf{e_{y}}}$ estimated using our augmented conical
emission model for a laser propagation affected by (a) no OA (lens of 2.6 m
focal length), (b) sphericity, (c,d) astigmatism and (e,f) coma. Key
parameters used for integration are detailed in table 2.

##  V Conclusion

In this paper, we have derived a so-called ”augmented” conical emission model,
able to satisfactorily reproduce angular THz spectra provided by 3D UPPE
simulations, both in focused propagation and filamentation regimes. This model
is called ”augmented” because it still relies on the key parameters found by
You et al. [47]. Yet, it includes complementary information on three other
crucial quantities, namely, (i) the spectral distribution of the THz source in
the nonlinear propagation regime of interest, (ii) the transverse profile and
phase of the beam when it undergoes lensing effect or more complex phase
distortions such as, e.g., those induced by deformable mirrors, and (iii) Kerr
contribution to the overall THz radiation when self-focusing may become
relevant, e.g., in the long-range filamentation regime. Our augmented model
also addressed the role of pulse envelopes, which has been shown to be
relatively limited, and it could be generalized to circularly polarized two-
color pulses.

We carefully and systematically compared the results obtained from our
augmented model with comprehensive UPPE simulation data. From this, we find an
alternative to running heavy simulations before evaluating the angular spectra
of the THz to interpret current experiments. Indeed, it appears sufficient to
model a class of spectral profiles appropriately with respect to the targetted
nonlinear propagation regime (tightly or loosely focused geometry) and
appreciate the weight of the Kerr contribution to get reliable angular
radiation patterns. Alternatively, LC-computed spectra could have been used,
or, if available, the 3D UPPE spectral profiles could directly be implemented
in the spectral amplitude $A(\omega)$. However, we have shown that simple test
functions for the spectral distributions are sufficient to obtain valuable
radiated intensities.

Importantly, instead of selecting one single radiated spectrum at a given
propagation distance, we probed, for the first time to our knowledge, the
whole spectral dynamics by examining changes in the radiated patterns for
various plasma lengths around the dephasing length $l_{d}$. It appears that a
good knowledge of the spectral dynamics along the propagation range covering
the THz emitting zone is much more instructive than only analyzing the
spectrum at the end of the plasma channel.

In connection with this property, our 3D UPPE simulations displayed evidence
that the choice of the actual plasma length $L$ must be given by the length
over which THz field emission takes place. This length is generally much
shorter than the overall plasma extent along the optical axis. By computing
the evolution of the relative phase between the two laser harmonics
numerically along this effective plasma length (Fig. 3), we have proven that
in focused propagation geometry $L$ remains close to $l_{d}$, which justifies
why most of the radiated spectra measured in such emission regimes [46, 52,
85] contain both on-axis emissions at small THz frequencies and net conical
emission at higher frequencies. In contrast, a collimated or loosely focused
propagation promotes plasma lengths $L\gg l_{d}$ and leads to conical emission
in the entire THz frequency range.

In conclusion, we think that the present study provides the terahertz
community with a deeper understanding of the THz spectral dynamics and
additional tools to predict angular THz radiation patterns delivered by two-
color air plasmas in future experiments.

##  Appendix A THz emission model

Here, we establish a general conical emission model that includes laser
envelope profiles in time and in the transverse plane. Let $t$ be the time
variable and $\mathrm{\mathbf{r}}$ be the position of the THz wave detector.
The frame’s origin is located at the center of the plasma with length $L$
assumed to be finite. In what follows, $\Theta$ corresponds to the polar angle
at which the radiation is emitted with respect to the optical axis; $\Phi$ is
the azimuthal angle.

Helmholtz’s theorem gives the vector potential
$\mathrm{\mathbf{A}}(\mathrm{\mathbf{r}},t)$ emitted by a plasma source of
volume $V$ with current density $\mathrm{\mathbf{J}}(\mathrm{\mathbf{r}},t)$
as

${\mathrm{\mathbf{A}}}({\mathrm{\mathbf{r}}},t)=\frac{\mu_{0}}{4\pi}\int_{V}d^{3}{\mathrm{\mathbf{r}}}^{\prime}\frac{{\mathrm{\mathbf{J}}}({\mathrm{\mathbf{r}}}^{\prime},t^{\prime})}{|{\mathrm{\mathbf{r}}}-{\mathrm{\mathbf{r}}}^{\prime}|},$
(50)

where $t^{\prime}=t-|{\mathrm{\mathbf{r}}}-{\mathrm{\mathbf{r}}}^{\prime}|/c$
corresponds to the delayed time at which the detector receives the signal
generated in $t$. The associated electric field
$\mathrm{\mathbf{E}}(\mathrm{\mathbf{r}},t)$, in the absence of a scalar
potential, follows from

${\mathrm{\mathbf{E}}}({\mathrm{\mathbf{r}}},t)=-\frac{\mu_{0}}{4\pi}\int_{V}d^{3}{\mathrm{\mathbf{r}}}^{\prime}\frac{\partial_{t}{\mathrm{\mathbf{J}}}({\mathrm{\mathbf{r}}}^{\prime},t^{\prime})}{|{\mathrm{\mathbf{r}}}-{\mathrm{\mathbf{r}}}^{\prime}|},$
(51)

and the magnetic field from

${\mathrm{\mathbf{B}}}({\mathrm{\mathbf{r}}},t)=\frac{\mu_{0}}{4\pi}\int_{V}d^{3}{\mathrm{\mathbf{r}}}^{\prime}\frac{\nabla\times{\mathrm{\mathbf{J}}}({\mathrm{\mathbf{r}}}^{\prime},t^{\prime})}{|{\mathrm{\mathbf{r}}}-{\mathrm{\mathbf{r}}}^{\prime}|}.$
(52)

For the photocurrent mechanism considered here, it is safe to assume that the
relevant current density is transverse and that no charge separation occurs
($\nabla\cdot\mathrm{\mathbf{J}}=0$), which justifies the use of Eq. (51) to
compute the radiated far-field power. We note that for the transition
Cherenkov (TC) mechanism [86, 87, 88, 89, 90, 50], where the relevant current
density is longitudinal and $\nabla\cdot\mathrm{\mathbf{J}}\neq 0$ produces a
charge density, one would need to account for a scalar potential in Eq. (51)
or resort to Eq. (52).

The temporal Fourier transform of Eq. (51) leads to

$\widehat{{\mathrm{\mathbf{E}}}}({\mathrm{\mathbf{r}}},\omega)=\frac{i\omega\mu_{0}}{4\pi}\int_{V}d^{3}{\mathrm{\mathbf{r}}}^{\prime}\frac{\mbox{e}^{i|{\mathrm{\mathbf{r}}}-{\mathrm{\mathbf{r}}}^{\prime}|\omega/c}}{|{\mathrm{\mathbf{r}}}-{\mathrm{\mathbf{r}}}^{\prime}|}\widehat{{\mathrm{\mathbf{J}}}}({\mathrm{\mathbf{r}}}^{\prime},\omega)\,.$
(53)

We assume that the distance between the plasma and the detector is large
compared to the filament dimensions, that is,
$|\mathrm{\mathbf{r}}|\gg|\mathrm{\mathbf{r^{\prime}}}|$. Cartesian,
cylindrical and spherical coordinates are written ($x$, $y$, $z$), ($r$,
$\varphi$, $z$) and ($r$, $\Theta$, $\Phi$). The notation ′ applies to the
local coordinates in the plasma. The far-field assumption leads to the
simplifications

$|{\mathrm{\mathbf{r}}}-{\mathrm{\mathbf{r}}}^{\prime}|\simeq
r-\frac{{\mathrm{\mathbf{r}}}^{\prime}\cdot{\mathrm{\mathbf{r}}}}{r}=r-\frac{c}{\omega}{\mathrm{\mathbf{k}}}_{\perp}\cdot{\mathrm{\mathbf{r}}^{\prime}}_{\perp}-z^{\prime}\cos\Theta$
(54)

with
${\mathrm{\mathbf{k}}}_{\perp}=(\omega/c)(\sin\Theta\cos\Phi,\sin\Theta\sin\Phi)$,
providing Eq. (25) of the main article.

To evaluate the source terms $\mathrm{\mathbf{J}}_{e}(\mathrm{\mathbf{r}},t)$
and $\mathrm{\mathbf{P}}_{K}(\mathrm{\mathbf{r}},t)$ in Eq. (26) for a two-
color laser pulse, we employ the general expression of the LP-P two-color
laser field [Eq. (27)] and apply the assumptions discussed in the main text.
The free-electron contribution
$\mathrm{\mathbf{J}}_{e}(\mathrm{\mathbf{r}},t)$ is calculated using the toy
ionization rate and assuming a small ionization yield,

$W(E)=CE^{2}(t),\quad
N_{e}(t)=CN_{a}\int_{-\infty}^{t}E^{2}(t^{\prime})dt^{\prime}.$ (55)

The square of the LP-P two-color laser field Eq. (27) with the assumptions
discussed in the main text reads

$E_{L}^{2}({\mathrm{\mathbf{r}}},t)=\mathcal{E}_{1}^{2}\cos^{2}{[\omega_{0}t^{\prime}_{1}+\phi_{1}}]+\mathcal{E}_{2}^{2}\cos^{2}{[2\omega_{0}t^{\prime}_{2}+\phi_{2}]}\\\
+2\mathcal{E}_{1}\mathcal{E}_{2}\cos{[2\omega_{0}t^{\prime}_{2}+\phi_{2}]}\cos{[\omega_{0}t^{\prime}_{1}+\phi_{1}]},$
(56)

where we omitted the arguments of
$\mathcal{E}_{1}(\mathrm{\mathbf{r}}_{\perp},t^{\prime}_{1})$,
$\mathcal{E}_{2}(\mathrm{\mathbf{r}}_{\perp},t^{\prime}_{2})$,
$\phi_{1}({\mathrm{\mathbf{r}}}_{\perp})$, and
$\phi_{2}({\mathrm{\mathbf{r}}}_{\perp})$ for conciseness. For the plasma
term, straightforward integration of Eq. (55) leads to linear and sine
temporal dependencies in time when slowly varying envelopes are assumed. The
product $N_{e}(t){\mathrm{\mathbf{E}}}_{L}(t)$ thus contains quasi-stationary
contributions in the terms

$\frac{CN_{a}}{4\omega_{0}}\mathcal{E}_{1}^{2}\sin[2\omega_{0}t^{\prime}_{1}+2\phi_{1}]\mathcal{E}_{2}\cos[2\omega_{0}t^{\prime}_{2}+\phi_{2}]\mathrm{\mathbf{e_{x}}}$

and

$\frac{CN_{a}}{\omega_{0}}\mathcal{E}_{1}\mathcal{E}_{2}\sin[\omega_{0}(2t^{\prime}_{2}-t^{\prime}_{1})+\phi_{2}-\phi_{1}]\mathcal{E}_{1}\cos[\omega_{0}t^{\prime}_{1}+\phi_{1}]\mathrm{\mathbf{e_{x}}}.$

The treatment of the Kerr term proceeds similarly from Eq. (7) without
integration in time and, therefore, only contains cosine contributions. After
retaining only the quasi-stationary terms, the two potential THz sources can
then be expressed as

$(\partial_{t}+\nu_{c})\mathrm{\mathbf{J}}_{e,{\rm
THz}}({\mathrm{\mathbf{r}}},t)={J}_{e,0}\
\mathcal{E}_{1}^{2}({\mathrm{\mathbf{r}}}_{\perp},t_{1}^{\prime})\mathcal{E}_{2}({\mathrm{\mathbf{r}}}_{\perp},t_{2}^{\prime})\\\
\times\sin{\left[\Delta
kz+\phi(\mathrm{\mathbf{r}}_{\perp})\right]}\mathrm{\mathbf{e_{x}}}$ (57)

using ${J}_{e,0}=3CN_{a}e^{2}/(8m_{e}\omega_{0})$, and

${\mathrm{\mathbf{P}}}_{K,{\rm THz}}({\mathrm{\mathbf{r}}},t)=P_{K,0}\
\mathcal{E}_{1}^{2}({\mathrm{\mathbf{r}}}_{\perp},t_{1}^{\prime})\mathcal{E}_{2}({\mathrm{\mathbf{r}}}_{\perp},t_{2}^{\prime})\\\
\times\cos{\left[\Delta
kz+\phi(\mathrm{\mathbf{r}}_{\perp})\right]}\mathrm{\mathbf{e_{x}}}$ (58)

with $P_{K,0}=3\varepsilon_{0}\chi^{(3)}/4$, $\Delta k=\pi/l_{d}$, and $l_{d}$
denoting the dephasing length given in Eq. (17). With the assumption of
separable laser envelopes in space and time [Eq. (28)] we can simplify the
relations (57) and (58), and the unified source term
$\widehat{\mathrm{\mathbf{J}}}(\mathrm{\mathbf{r}},\omega)$, cf. Eq. (26),
reads:

$\widehat{\mathrm{\mathbf{J}}}(\mathrm{\mathbf{r}},\omega)=\widehat{\mathrm{\mathbf{J}}_{0}}(\mathrm{\mathbf{r}},\omega)J_{\rm
tr}(\mathrm{\mathbf{r}}_{\perp})\mathcal{F}_{t}\\!\left[{\cal{E}}_{ax,1}^{2}(t_{1}^{\prime}){\cal{E}}_{ax,2}(t_{2}^{\prime})\right]\\!(\omega),$
(59)

where $J_{\rm tr}(\mathrm{\mathbf{r}}_{\perp})=\mathcal{E}_{\rm
tr,1}^{2}({\mathrm{\mathbf{r}}}_{\perp})\mathcal{E}_{\rm
tr,2}({\mathrm{\mathbf{r}}}_{\perp})$ [see Eq. (38)], and

$\widehat{\mathrm{\mathbf{J}}_{0}}(\mathrm{\mathbf{r}}_{\perp},z,\omega)=\frac{J_{e,0}}{\nu_{c}-i\omega}\sin{(\Delta
kz+\phi(\mathrm{\mathbf{r}}_{\perp}))}\mathrm{\mathbf{e_{x}}}\\\ -i\omega
P_{K,0}\cos{(\Delta
kz+\phi(\mathrm{\mathbf{r}}_{\perp}))}\mathrm{\mathbf{e_{x}}}.$ (60)

$\mathcal{F}_{t}$ refers to Fourier transform in time, and using
$t_{j=\\{1,2\\}}^{\prime}=t-zn(j\omega_{0})/c$ we can rewrite

$\displaystyle\quad{\cal
F}_{t}\\!\left[{\cal{E}}_{ax,1}^{2}(t_{1}^{\prime}){\cal{E}}_{ax,2}(t_{2}^{\prime})\right]\\!(\omega)=\int{\cal{E}}_{ax,1}^{2}(t_{1}^{\prime}){\cal{E}}_{ax,2}(t_{2}^{\prime})e^{i\omega
t}dt$
$\displaystyle=e^{ikn(\omega_{0})z}\int{\cal{E}}_{ax,1}^{2}(t_{1}^{\prime}){\cal{E}}_{ax,2}\\!\left(t_{1}^{\prime}-\frac{\pi
z}{2\omega_{0}l_{d}}\right)e^{i\omega t_{1}^{\prime}}dt_{1}^{\prime},$

with $k=\omega/c$ being the THz wavenumber in vacuum. We can use the
convolution theorem to rewrite the integral as

$\displaystyle\quad\int{\cal{E}}_{ax,1}^{2}(t){\cal{E}}_{ax,2}(t-\gamma)e^{i\omega
t}dt$
$\displaystyle=\frac{1}{2\pi}\iiint{\cal{E}}_{ax,1}^{2}(t)e^{i\omega^{\prime}t}{\cal{E}}_{ax,2}(t^{\prime}-\gamma)e^{i(\omega^{\prime}-\omega)t^{\prime}}dtdt^{\prime}d\omega^{\prime}$
$\displaystyle=\frac{1}{2\pi}\iiint{\cal{E}}_{ax,1}^{2}(t)e^{i\omega^{\prime}t}{\cal{E}}_{ax,2}(\tilde{t})e^{i(\omega^{\prime}-\omega)(\tilde{t}+\gamma)}dtd\tilde{t}d\omega^{\prime}$
$\displaystyle=\frac{e^{-i\omega\gamma}}{2\pi}\int{\cal
F}_{t}\\!\left[{\cal{E}}_{ax,1}^{2}(t)\right]\\!(\omega^{\prime}){\cal
F}_{t}\\!\left[{\cal{E}}_{ax,2}({t})\right]\\!(\omega^{\prime}-\omega)e^{i\omega^{\prime}\gamma}d\omega^{\prime}$

where $\gamma=\pi z/(2\omega_{0}l_{d})$. Provided that the spectra of the
temporal envelopes are sufficiently narrow (not too short pulses), we can
conclude that the integrand is nonzero for $\omega^{\prime}\ll\omega_{0}$
only. Then, for plasma lengths of at most a few $l_{d}$ we can neglect the
$z$-dependent term $\exp(i\omega^{\prime}\gamma)$ under the integral. Because
we are interested in THz frequencies only ($\omega\ll\omega_{0}$), with the
same argument we can neglect the factor $\exp(-i\omega\gamma)$ in front of the
integral. Thus, we have established that

$\quad{\cal
F}_{t}\\!\left[{\cal{E}}_{ax,1}^{2}(t_{1}^{\prime}){\cal{E}}_{ax,2}(t_{2}^{\prime})\right]\\!(\omega)\approx
e^{ikn(\omega_{0})z}C_{\rm env}(\omega)$ (61)

with

$C_{\rm
env}(\omega)=\int{\cal{E}}_{ax,1}^{2}(t_{1}^{\prime}){\cal{E}}_{ax,2}(t_{1}^{\prime})e^{i\omega
t_{1}^{\prime}}dt_{1}^{\prime}.$ (62)

While this expression is valid for arbitrary (slowly varying) pulse envelopes,
it is instructive to inspect it for Gaussian envelopes:

${\cal{E}}_{ax,j}(t)\propto\mbox{e}^{-(t-t_{c}^{j})^{2}/t_{j}^{2}}$ (63)

with $t_{j}$ being the duration of the $j^{th}$ color envelope and $t_{c}^{j}$
its temporal offset. After straightforward computations, we obtain a Gaussian
power spectrum

$C_{env}(\omega)\propto\frac{t_{1}t_{2}\sqrt{\pi}}{t_{eq}}e^{-\frac{2\Delta
t^{2}}{t_{eq}^{2}}}e^{i\omega(t_{c}^{1}-\Delta
t\frac{t_{1}^{2}}{t_{eq}^{2}})}e^{-\frac{t_{1}^{2}t_{2}^{2}}{4t_{eq}^{2}}\omega^{2}},$
(64)

where $t_{eq}=\sqrt{2t_{2}^{2}+t_{1}^{2}}$ and $\Delta t=t_{c}^{1}-t_{c}^{2}$
is the temporal separation of the two envelopes. Most notably, the amplitude
strongly decreases with increasing $\Delta t^{2}$. Interestingly, deviations
from the above slowly-varying assumptions and $\omega\ll\omega_{0}$ would
produce additional contributions of orders $\Delta
tL/(\omega_{0}t_{eq}^{2}l_{d})$ and $\omega\lambda_{\rm
THz}/(\omega_{0}l_{d})$, respectively, in the argument of the sinc function of
Eq. (20). Deviations from the above requirement $\omega\gamma\ll 1$ would
introduce a $z$-dependent correction $e^{-i\gamma\omega t_{1}^{2}/t_{eq}^{2}}$
into $C_{env}$. For a $0.8\,\mu$m fundamental pulse, this contribution remains
small for THz frequencies $<60$ THz (120 THz) as long as $L\leq 4\,l_{d}$
(resp. $L\leq 2\,l_{d}$).

The source term $\widehat{\mathrm{\mathbf{J}}}(\mathrm{\mathbf{r}},\omega)$
hence reads

$\widehat{\mathrm{\mathbf{J}}}(\mathrm{\mathbf{r}},\omega)\approx\widehat{\mathrm{\mathbf{J}}_{0}}(\mathrm{\mathbf{r}},\omega)J_{\rm
tr}(\mathrm{\mathbf{r}}_{\perp})e^{ikzn(\omega_{0})}C_{\rm env}(\omega).$ (65)

By injecting relations (65) into equation (25), the radiated electric field
takes the form:

$\widehat{\mathrm{\mathbf{E}}}(r,\Theta,\Phi,\omega)=\frac{\mu_{0}}{4\pi
r}\mbox{e}^{ikr}\sum_{j=J,K}C_{j}(\Theta,\Phi,\omega),$ (66)

where

$C_{J,K}(\Theta,\Phi,\omega)=A_{J,K}(\omega)\mathrm{\mathbf{I}}_{J,K}(\Theta,\Phi,\omega).$
(67)

Here, $A_{J,K}(\omega)$ are the THz spectral profiles associated with the
photocurrents and Kerr response,

$\displaystyle A_{J}(\omega)$ $\displaystyle=i\omega
C_{env}(\omega)\frac{J_{e,0}}{\nu_{c}-i\omega},$ (68) $\displaystyle
A_{K}(\omega)$ $\displaystyle=\omega^{2}C_{env}(\omega)P_{K,0},$ (69)

while $\mathrm{\mathbf{I}}_{J,K}(\Theta,\Phi,\omega)$ describes the
interferences between the THz waves generated along the plasma channel:

$\displaystyle\mathrm{\mathbf{I}}_{J}(\Theta,\Phi,\omega)$
$\displaystyle=\int_{-\frac{L}{2}}^{\frac{L}{2}}dz^{\prime}\mbox{e}^{ikz^{\prime}[n(\omega_{0})-\cos{\Theta}]}$
$\displaystyle\times\frac{\mbox{e}^{i\Delta
kz^{\prime}}I_{S}^{+}(\Theta,\Phi,\omega)-\mbox{e}^{-i\Delta
kz^{\prime}}I_{S}^{-}(\Theta,\Phi,\omega)}{2i}\mathrm{\mathbf{e_{x}}},$
$\displaystyle\mathrm{\mathbf{I}}_{K}(\Theta,\Phi,\omega)$
$\displaystyle=\int_{-\frac{L}{2}}^{\frac{L}{2}}dz^{\prime}\mbox{e}^{ikz^{\prime}[n(\omega_{0})-\cos{\Theta}]}$
$\displaystyle\times\frac{\mbox{e}^{i\Delta
kz^{\prime}}I_{S}^{+}(\Theta,\Phi,\omega)+\mbox{e}^{-i\Delta
kz^{\prime}}I_{S}^{-}(\Theta,\Phi,\omega)}{2}\mathrm{\mathbf{e_{x}}}.$

The quantities $I_{S}^{\pm}(\Theta,\Phi,\omega)$ refer to the transverse
integral contributions to the THz spectrum:

$I_{S}^{\pm}(\Theta,\Phi,\omega)=\iint
d^{2}\mathrm{\mathbf{r}}^{\prime}_{\perp}J_{\rm
tr}(\mathrm{\mathbf{r}}_{\perp}^{\prime})e^{\pm
i\phi(\mathrm{\mathbf{r}}_{\perp}^{\prime})}\mbox{e}^{-i{\mathrm{\mathbf{k}}}_{\perp}\cdot{\mathrm{\mathbf{r}}}_{\perp}^{\prime}},$
(70)

with
${\mathrm{\mathbf{k}}}_{\perp}=(\omega/c)(\sin\Theta\cos\Phi,\sin\Theta\sin\Phi)$.
The functions $I_{S}^{\pm}(\Theta,\Phi,\omega)$ can be rewritten as

$I_{S}^{\pm}(\Theta,\Phi,\omega)=\mathcal{F}_{\perp}\left[J_{\rm
tr}(\mathrm{\mathbf{r}}_{\perp}^{\prime})e^{\pm
i\phi(\mathrm{\mathbf{r}}_{\perp}^{\prime})}\right](\Theta,\Phi,\omega),$ (71)

where $\mathcal{F_{\perp}}$ represents the transverse Fourier transform.

It can be checked that the elementary transformations

$\displaystyle\sin[\Delta
kz+\phi(\mathrm{\mathbf{r}}_{\perp}^{\prime})]\mathrm{\mathbf{e_{x}}}$
$\displaystyle\to\binom{\sin[\Delta
kz+\phi(\mathrm{\mathbf{r}}_{\perp}^{\prime})]}{\cos[\Delta
kz+\phi(\mathrm{\mathbf{r}}_{\perp}^{\prime})]}$ (72)
$\displaystyle\cos[\Delta
kz+\phi(\mathrm{\mathbf{r}}_{\perp}^{\prime})]\mathrm{\mathbf{e_{x}}}$
$\displaystyle\to\binom{\cos[\Delta
kz+\phi(\mathrm{\mathbf{r}}_{\perp}^{\prime})]}{-\sin[\Delta
kz+\phi(\mathrm{\mathbf{r}}_{\perp}^{\prime})]}$ (73)

allow one to treat similarly CP-S laser pulses, cf. Eq. (2).

Straightforward calculation of the integrals
$\mathrm{\mathbf{I}}_{J,K}(\Theta,\Phi,\omega)$ gives, for a LP-P laser pulse:

$\displaystyle\mathrm{\mathbf{I}}_{J}^{\rm LP-P}(\Theta,\Phi,\omega)$
$\displaystyle=\frac{L}{2}\left[I_{S}^{+}\kappa_{+}-I_{S}^{-}\kappa_{-}\right]\mathrm{\mathbf{e_{x}}},$
(74) $\displaystyle\mathrm{\mathbf{I}}_{K}^{\rm LP-P}(\Theta,\Phi,\omega)$
$\displaystyle=i\frac{L}{2}\left[I_{S}^{+}\kappa_{+}+I_{S}^{-}\kappa_{-}\right]\mathrm{\mathbf{e_{x}}},$
(75)

and for a CP-S pulse:

$\displaystyle\mathrm{\mathbf{I}}_{J}^{\rm CP-S}(\Theta,\Phi,\omega)$
$\displaystyle=\frac{L}{2}\binom{I_{S}^{+}\kappa_{+}-I_{S}^{-}\kappa_{-}}{I_{S}^{+}\kappa_{+}+I_{S}^{-}\kappa_{-}},$
(76) $\displaystyle\mathrm{\mathbf{I}}_{K}^{\rm CP-S}(\Theta,\Phi,\omega)$
$\displaystyle=i\frac{L}{2}\binom{I_{S}^{+}\kappa_{+}+I_{S}^{-}\kappa_{-}}{-I_{S}^{+}\kappa_{+}+I_{S}^{-}\kappa_{-}}.$
(77)

The quantities $\kappa_{\pm}$ are defined by:

$\displaystyle\kappa_{\pm}$ $\displaystyle=\mathrm{sinc}(\alpha_{\pm})$ (78)
$\displaystyle\alpha_{\pm}$
$\displaystyle=\frac{1}{2}kL\left(n_{0}-\cos\Theta\pm\frac{\lambda_{\rm
THz}}{2l_{d}}\right),$ (79)

with $\mathrm{sinc}(x)=\sin x/x$ and $k=\omega/c=2\pi/\lambda_{\rm THz}$.

Equations (62), (74) and (75) [alternatively (76) and (77)] can finally be
injected into the relation (66) to obtain the far-field spectral intensity:

$\left|\widehat{\mathrm{\mathbf{E}}}(r,\Theta,\Phi,\omega)\right|^{2}=\frac{\mu_{0}^{2}}{16\pi^{2}r^{2}}\left|\sum_{j=J,K}C_{j}(\Theta,\Phi,\omega)\right|^{2}.$
(80)

We have implicitly assumed that the dephasing length $l_{d}$ is constant in
space, although the source terms of THz radiation include radial variations.
This means that similarly to You et al.’s model, we need to select an
effective electron density level $N_{e}^{\rm eff}$ to estimate $l_{d}$ and
thus evaluate our model.

Note that similarly to [47], we ignore plasma absorption and opacity in the
THz frequency range. These assumptions can be justified as follows. By
neglecting the Kerr contribution and focusing on plasma dispersion, Eq. (18)
reads

$\begin{split}n^{2}(\omega)&=n^{2}_{\rm
opt}(\omega)-\frac{\omega_{p}^{2}}{\omega^{2}(1+i\nu_{c}/\omega)}\\\
&=n^{2}_{\rm
opt}(\omega)-\frac{\omega_{p}^{2}}{\omega^{2}+\nu_{c}^{2}}+i\frac{\nu_{c}}{\omega}\frac{\omega_{p}^{2}}{\omega^{2}+\nu_{c}^{2}}.\end{split}$
(81)

The real and imaginary parts of the medium index are given by

$\displaystyle\Re[n(\omega)]$
$\displaystyle=\left(\frac{|n^{2}(\omega)|+\Re[n^{2}(\omega)]}{2}\right)^{1/2},$
(82) $\displaystyle\Im[n(\omega)]$
$\displaystyle=\left(\frac{|n^{2}(\omega)|-\Re[n^{2}(\omega)]}{2}\right)^{1/2}.$
(83)

For $\Re[n^{2}(\omega)]>0$, attenuation of THz waves through $\Im[n(\omega)]$
is due to plasma absorption ($\nu_{c}\approx 3$ THz), which recovers Kim et
al.’s[49] approximation
$\Im[n(\omega)]\approx\nu_{c}\omega_{p}^{2}/(2\omega(\omega^{2}+\nu_{c}^{2}))$
in the limits $\nu_{c}/\omega\ll 1,\,n^{2}_{\rm opt}(\omega)\approx 1$. At
density levels $N_{e}\lesssim 10^{17}$ cm-3, the absorption length
$L_{a}=c/(\omega\Im[n(\omega)])$ is of the order of $\sim 100$ µm at
$\omega=2\pi\nu_{c}$. As typical filaments have a radius of $a\approx 20$ µm
and THz radiation is emitted at angles $\Theta\approx 5$°, the propagation
distance of those waves to exit the plasma is $\approx 200$ µm $\sim L_{a}$,
meaning that plasma absorption should already have a limited effect on the
emitted THz field. At higher frequencies, $L_{a}$ quickly increases ($\sim
500$ µm at 5 THz), and absorption becomes negligible.

In the lower region of the THz spectrum, $\Re[n^{2}(\omega)]=n^{2}_{\rm
opt}(\omega)-\omega_{p}^{2}/(\omega^{2}+\nu_{c}^{2})$ may become negative.
Then, the attenuation of the THz waves through $\Im[n(\omega)]$ is due to the
plasma opacity. However, this effect barely impacts the forward emitted
spectrum, as the THz field is emitted at the front of the plasma, i.e., before
the plasma builds up [91, 92], so that the local $N_{e}$ and thus $\omega_{p}$
take smaller values.

##  Appendix B Derivation of You et al.’s model

The original You et al.’s model [Eq. (19)] [47, 48] can be easily derived from
our augmented emission model. We assume a homogeneous cylindrical plasma of
length $L$ and radius $a$. Thus, the source term $J_{\rm tr}$ contained in
$I_{S}^{\pm}(\Theta,\Phi,\omega)$, is constant on the disk of radius
$r_{\perp}<a$ and zero outside. Moreover, we introduce $\phi_{0}$ as the
constant value taken by the phase $\phi(\mathrm{\mathbf{r}}_{\perp})$. Hence,
we can write

$\displaystyle I_{S}^{\pm}$ $\displaystyle=e^{\pm i\phi_{0}}\iint
d^{2}\mathrm{\mathbf{r}}_{\perp}^{\prime}e^{-i\mathrm{\mathbf{k}}_{\perp}\cdot\mathrm{\mathbf{r}}_{\perp}^{\prime}}$
$\displaystyle=e^{\pm
i\phi_{0}}\int_{0}^{a}dr_{\perp}^{\prime}r_{\perp}^{\prime}\int_{0}^{2\pi}d\varphi^{\prime}e^{-ik_{\perp}r_{\perp}^{\prime}\cos\varphi^{\prime}}$
$\displaystyle=2\pi e^{\pm
i\phi_{0}}\int_{0}^{a}dr_{\perp}^{\prime}r_{\perp}^{\prime}J_{0}(k_{\perp}r_{\perp}^{\prime}),$

with $J_{0}(x)$ being the Bessel function of first kind and rank 0. The
simplifications continue by using the identity $\int dxxJ_{0}(x)=xJ_{1}(x)$
[93], where $J_{1}(x)$ describes the Bessel function of first kind and rank 1:

$I_{S}^{\pm}(\Theta,\omega)=2\pi e^{\pm
i\phi_{0}}a^{2}\frac{J_{1}(\beta)}{\beta},$ (84)

where we used $k_{\perp}=k\sin\Theta$ and $\beta=ak\sin\Theta$.

If we now neglect the contribution of the Kerr response to the THz radiation
($A_{K}(\omega)=0$), assume that the laser pulse is linearly polarized (LP-P)
and $n(\omega_{0})=1$, we can evaluate Eq. (80) as

$\left|\widehat{E}(r,\Theta,\omega)\right|^{2}=\frac{\mu_{0}^{2}}{16r^{2}}a^{4}L^{2}\left|A_{J}(\omega)\right|^{2}\left(\frac{J_{1}(\beta)}{\beta}\right)^{2}\\\
\times\left|\mbox{e}^{i\phi_{0}}\kappa_{+}-\mbox{e}^{-i\phi_{0}}\kappa_{-}\right|^{2},$

where the term depending on $\alpha_{\pm}$ develops as:

$\left|\mbox{e}^{i\phi_{0}}\kappa_{+}-\mbox{e}^{-i\phi_{0}}\kappa_{-}\right|^{2}=\kappa_{+}^{2}+\kappa_{-}^{2}-2\kappa_{+}\kappa_{-}\cos(2\phi_{0}).$
(85)

The far-field THz intensity radiated by the plasma thus restores You et al.’s
Eq. (19), where $A(\omega)\equiv A_{J}(\omega)$.

## Author’s Contributions

All authors contributed equally to this work.

## Acknowledgments

## Data availability

The data supporting this study’s findings are available within the present
article. Complementary data can be made available from the corresponding
authors upon reasonable request.

## References

  * Tonouchi [2007] M. Tonouchi, “Cutting-edge terahertz technology,” Nature Photonics 1, pp. 97 – 105 (2007).
  * Currie  _et al._ [1997] N. C. Currie, F. J. Demma, D. D. Ferris Jr., R. W. McMillan, M. C. Wicks, and K. Zyga, “Imaging sensor fusion for concealed weapon detection,” Investigative Image Processing 2942 (1997).
  * Federici [2005] J. F. Federici, “Thz imaging and sensing for security applications - explosives, weapons and drugs,” Semicond. Sci. Technol. 20, pp. S266 – S280 (2005).
  * Cole  _et al._ [2001] B. E. Cole, R. M. Woodward, D. A. Crowley, V. P. Wallace, D. D. Arnone, and M. Pepper, “Terahertz imaging and spectroscopy of human skin in vivo,” Commercial and Biomedical Applications of Ultrashort Pulse Lssers 4276 (2001).
  * Woodward  _et al._ [2003] R. M. Woodward, V. P. Wallace, D. D. Arnone, E. H. Linfield, and M. Pepper, “Terahertz pulsed imaging of skin cancer in the time and frequency domain,” Journal of Biological Physics 29, pp. 257 – 261 (2003).
  * Pickwell  _et al._ [2004] E. Pickwell, B. E. Cole, A. J. Fitzgerald, M. Pepper, and V. P. Wallace, “In vivo study of human skin using pulsed terahertz radiation,” Phys. Med. Biol. 49, pp. 1595 – 1607 (2004).
  * Fukunaga  _et al._ [2007] K. Fukunaga, Y. Ogawa, S. Hayashi, and I. Hosako, “Terahertz spectroscopy for art conservation,” IEICE Electronics Express 4, pp. 258 – 263 (2007).
  * Yasui  _et al._ [2005] T. Yasui, T. Tasuda, K. Sawanaka, and T. Araki, “Terahertz paintmeter for noncontact monitoring of thickness and drying progress in paint film,” Appl. Opt. 44, pp. 6849–6856 (2005).
  * Koch Dandolo and Jepsen [2015] C. L. Koch Dandolo and P. U. Jepsen, “Wall painting investigation by means of non-invasive terahertz time-domain imaging (THz-TDI): Inspection of subsurface structures buried in historical plasters,” Journal of Infrared, Millimeter, and Terahertz Waves 37, 198–208 (2015).
  * Koch Dandolo  _et al._ [2015] C. L. Koch Dandolo, T. Filtenborg, J. Skou-Hansen, and P. U. Jepsen, “Analysis of a seventeenth-century panel painting by reflection terahertz time-domain imaging (THz-TDI): contribution of ultrafast optics to museum collections inspection,” Applied Physics A 121, 981–986 (2015).
  * Bianca Jackson  _et al._ [2011] J. Bianca Jackson, J. Bowen, G. Walker, J. Labaune, G. Mourou, M. Menu, and K. Fukunaga, “A survey of terahertz applications in cultural heritage conservation science,” IEEE Transactions on Terahertz Science and Technology 1, pp. 220–231 (2011).
  * Yin  _et al._ [2007] X. Yin, S. Hadjiloucas, B. M. Fischer, B. W.-H. Ng, H. M. Paiva, R. K. H. Galvao, G. C. Walker, J. W. Bowen, and D. Abbott, “Classification of lactose and mandelic acid thz spectra using subspace and wavelet-packet algorithms,” SPIE 6798, 679814 (2007).
  * Thrane  _et al._ [1995] L. Thrane, R. H. Jacobsen, P. U. Jepsen, and S. R. Keiding, “Thz reflection spectroscopy of liquid water,” Chem. Phys. Lett. 240, pp. 330 – 333 (1995).
  * Rønne, Åstrand, and Keiding [1999] C. Rønne, P. Åstrand, and S. R. Keiding, “THz spectroscopy of Liquid H2O and D2O,” Phys. Rev. Lett. 84, pp. 2888 – 2891 (1999).
  * Fattinger and Grischkowsky [1988] C. Fattinger and D. Grischkowsky, “Point source terahertz optics,” Appl. Phys. Lett. 53, 1480 (1988).
  * van Exter, Fattinger, and Grischkowsky [1989] M. van Exter, C. Fattinger, and D. Grischkowsky, “Terahertz time-domain spectroscopy of water vapor,” Optics Letters 14, pp. 1128 – 1130 (1989).
  * Walther  _et al._ [2000] M. Walther, B. Fischer, M. Schall, H. Helm, and P. U. Jepsen, “Far-infrared vibrational spectra of all-trans, 9-cis and 13-cis retinal measured by thz time-domain spectroscopy,” Chemical Physics Letters 332, 389–395 (2000).
  * Neu and Schmuttenmaer [2018] J. Neu and C. A. Schmuttenmaer, “Tutorial: An introduction to terhaertz time domain spectroscopy (thz-tds),” J. Appl. Phys 124, 231101 (2018).
  * Hu and Nuss [1995] B. B. Hu and M. C. Nuss, “Imaging with terahertz waves,” Optics Letters 20, 1716–1718 (1995).
  * Gente and Koch [2015] R. Gente and M. Koch, “Monitoring leaf water content with thz and sub-thz waves,” Plant Methods 00, pp. 11 – 15 (2015).
  * Gente  _et al._ [2018] R. Gente, A. Rehn, T. Probst, E.-M. Stübling, E. C. Camus, A. A. Covarrubias, J. C. Balzer, and M. Koch, “Outdoor measurements of leaf water content using thz quasi time-domain spectroscopy,” J. Infrared Milli. Terahz Waves 39, pp. 943 – 948 (2018).
  * Curcio  _et al._ [2018] A. Curcio, A. Marocchino, V. Dolci, S. Lupi, and M. Petrarca, “Resonant plasma excitation by single-cycle THz pulses,” Sci. Rep. 8, 1052 (2018).
  * Zhang  _et al._ [2018] D. Zhang, A. Fallahi, M. Hemmer, X. Wu, M. Fakhari, Y. Hua, H. Cankaya, A.-L. Calendron, L. E. Zapata, N. H. Matlis, and F. X. Kärtner, “Segmented terahertz electron accelerator and manipulator (steam),” Nature Photonics 12, 336–342 (2018).
  * Zaks, Liu, and Sherwin [2012] B. Zaks, R. B. Liu, and M. S. Sherwin, “Experimental observation of electron-hole recollisions,” Nature 483, 580–3 (2012).
  * Schubert  _et al._ [2014] O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek, C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W. Koch, and R. Huber, “Sub-cycle control of terahertz high-harmonic generation by dynamical bloch oscillations,” Nature Photonics 8, 119–123 (2014).
  * Langer  _et al._ [2016] F. Langer, M. Hohenleutner, C. P. Schmid, C. Poellmann, P. Nagler, T. Korn, C. Schüller, M. S. Sherwin, U. Huttner, J. T. Steiner, S. W. Koch, M. Kira, and R. Huber, “Lightwave-driven quasiparticle collisions on a subcycle timescale,” Nature 533, 225–229 (2016).
  * Shalaby, Vicario, and Hauri [2016] M. Shalaby, C. Vicario, and C. P. Hauri, “Simultaneous electronic and the magnetic excitation of a ferromagnet by intense THz pulses,” New J. Phys. 18, 013019 (2016).
  * Burford and El-Shenawee [2017] N. M. Burford and M. O. El-Shenawee, “Review of terahertz photoconductive antenna technology,” Optical Engineering 56, 010901 (2017).
  * Park  _et al._ [2012] S. Park, Y. Choi, Y. Oh, and K. Jeong, “Terahertz photoconductive antenna with metal nanoislands,” Optics Express 20, pp. 25530 – 25535 (2012).
  * Zhong  _et al._ [2017] K. Zhong, W. Shi, D. Xu, P. Liu, Y. Wang, J. Mei, C. Yan, S. Fu, and J. Yao, “Optically pumped terahertz sources,” Sci. China Tech. Sci. 60, pp. 1801 – 1818 (2017).
  * Kohler  _et al._ [2002] R. Kohler, A. Tredicucci, F. Beltram, H. E. Beere, E. H. Linfield, A. G. Davies, D. A. Ritchie, R. C. Lotti, and F. Rossi, “Terahertz semiconductor-heterostructure laser,” Nature 417, pp. 156 – 159 (2002).
  * Williams [2007] B. S. Williams, “Terahertz quantum-cascade lasers,” Nature Photonics 1, pp. 517 – 525 (2007).
  * Li  _et al._ [2014] L. Li, L. Chen, J. Zhu, J. Freeman, P. Dean, A. Valavanis, A. G. Davies, and E. H. Linfield, “Terahertz quantum cascade lasers with $>1$ W output powers,” Electronic Letters 50, pp. 309 – 311 (2014).
  * Vicario, Monoszlai, and Hauri [2012] C. Vicario, B. Monoszlai, and C. P. Hauri, “GV/m single-cycle terahertz fields from a laser-driven large-size partitioned organic crystal,” Phys. Rev. Lett. 112, 213901 (2012).
  * Cook and Hochstrasser [2000] D. J. Cook and R. M. Hochstrasser, “Intense terahertz pulses by four-wave rectification in air,” Optics Letters 25, pp. 1210 – 1212 (2000).
  * Koulouklidis  _et al._ [2020] A. D. Koulouklidis, C. Gollner, V. Shumakova, V. Y. 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).
  * González 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).
  * Stathopulos, Skupin, and Bergé [2021] A. Stathopulos, S. Skupin, and L. Bergé, “Terahertz pulse generation by multi-color laser fields with linear versus circular polarization,” Optics Letters 46, pp. 5906 – 5909 (2021).
  * Liu  _et al._ [2020] S. Liu, Z. Fan, C. Lu, J. Gui, C. Luo, S. 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).
  * 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,” Journal of Applied Physics 125, 173103 (2019).
  * Keldysh [1965] L. V. Keldysh, “Ionization in the field of a strong electromagnetic wave,” Sov. Phys. JETP 20, 1307 (1965).
  * Daigle  _et al._ [2012] J.-F. Daigle, F. Théberge, M. Henriksson, T.-J. Wang, S. Yuan, M. Châteauneuf, J. Dubois, M. Piché, and S. L. Chin, “Remote thz generation from two-color filamentation: long distance dependence,” Opt.Express 20, 6825 (2012).
  * Bergé  _et al._ [2007] L. Bergé, S. Skupin, R. Nuter, J. Kasparian, and J.-P. Wolf, “Ultrashort filaments of light in weakly ionized, optically transparent media,” Reports on Progress in Physics 70, 1633 (2007).
  * Bergé  _et al._ [2019] L. Bergé, K. Kaltenecker, S. Engelbrecht, A. Nguyen, S. Skupin, L. Merlat, B. Fischer, B. Zhou, I. Thiele, and P. U. Jepsen, “Terahertz spectroscopy from air plasmas created by two-color femtosecond laser pulses: The ALTESSE project,” EPL 126, 24001 (2019).
  * Talbi  _et al._ [2023] A. Talbi, B. Zhou, P. U. Jepsen, S. Skupin, A. Courjaud, and L. Bergé, “Remote terahertz spectroscopy from extended two-color plasma filaments: The altesse 2 project,” EPL 143, 10001 (2023).
  * Andreeva  _et al._ [2016] V. A. Andreeva, O. G. Kosareva, N. A. Panov, D. E. Shipilo, P. M. Solyankin, M. N. Esaulkov, P. González de Alaiza Martínez, A. P. Shkurinov, V. A. Makarov, L. Bergé, and S. L. Chin, “Ultrashort terahertz spectrum generation from an air-based filament plasma,” Phys. Rev. Lett. 116, 063902 (2016).
  * You, Oh, and Kim [2012] Y. S. You, T. I. Oh, and K. Y. Kim, “Off-axis phase-matched terahertz emission from two-color laser-induced plasma filaments,” Phys. Rev. Lett. 109, 183902 (2012).
  * Oh  _et al._ [2013] T. I. Oh, Y. S. You, N. Jhajj, E. W. Rosenthal, H. M. Milchberg, and K. Y. Kim, “Intense terahertz generation in two-color laser filamentation: energy scaling with terawatt laser systems,” New J. Phys. 15, 075002 (2013).
  * Kim  _et al._ [2008] K. Y. Kim, A. J. Taylor, J. H. Glownia, and G. Rodriguez, “Coherent control of terahertz supercontinuum generation in ultrafast laser-gas interaction,” Nature Photonics 2, pp. 605 – 609 (2008).
  * Thiele  _et al._ [2017] I. Thiele, P. González de Alaiza Martínez, R. Nuter, A. Nguyen, L. Bergé, and S. Skupin, “Broadband terahertz emission from two-color femtosecond-laser-induced microplasmas,” Phys. Rev. A 96, 053814 (2017).
  * Gorodetsky  _et al._ [2014] A. Gorodetsky, A. D. Koulouklidis, M. Massaouti, and S. Tzortzakis, “Physics of the broadband terahertz emission from two-color laser-induced plasma filaments,” Phys. Rev. A 89, 033838 (2014).
  * Ushakov  _et al._ [2018] A. A. Ushakov, P. A. Chizhov, V. A. Andreeva, N. A. Panov, D. E. Shipilo, M. Matoba, N. Nemoto, N. Kanda, K. Konishi, V. V. Bukin, M. Kuwata-Gonokami, O. G. Kosareva, S. V. Garnov, and A. B. Savel’ev, “Ring and unimodal angular-frequency distribution of THz emission from two-color femtosecond plasma spark,” Opt. Express 26, pp. 18202 – 18213 (2018).
  * Kolesik and Moloney [2004] M. Kolesik and J. V. Moloney, “Nonlinear optical pulse propagation simulation: From Maxwell’s to unidirectional equations,” Phys. Rev. E 70, 036604 (2004).
  * Noll [1976] R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am. 66 (1976).
  * Meng  _et al._ [2016] C. Meng, W. Chen, X. Wang, Z. 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).
  * Tailliez  _et al._ [2020a] 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 (2020a).
  * Landau and Lifshitz [1965] L. D. Landau and E. M. Lifshitz,  _Quantum Mechanics_ (Pergamon, New York, 1965).
  * Perelomov, Popov, and Teret’ev [1966] A. M. Perelomov, V. S. Popov, and M. V. Teret’ev, “Ionization of atoms in an alternating electric field: I,” Sov. Phys. JETP 23, 924 (1966).
  * 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).
  * Nguyen  _et al._ [2017] 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,” Optics Express 25, pp. 4720 – 4740 (2017).
  * Boyd [2008] R. Boyd,  _Nonlinear Optics_ (Academic Press, San Diego, 2008).
  * Peck and Reeder [1972] E. R. Peck and K. Reeder, “Dispersion of air,” J. Opt. Soc. Am. A 62, 958 (1972).
  * Stathopulos, Skupin, and Bergé [2023] A. Stathopulos, S. Skupin, and L. Bergé, “A Raman-delayed nonlinearity for elliptically polarized ultrashort optical pulses,” Eur. Phys. J. Spec. Top. 232, 2285–2292 (2023).
  * Peñano  _et al._ [2003] J. R. Peñano, P. Sprangle, P. Serafim, B. Hafizi, and A. Ting, “Stimulated Raman scattering of intense laser pulses in air,” Phys. Rev. E 68, 056502 (2003).
  * Pitts  _et al._ [2004] T. A. Pitts, T. S. Luk, J. K. Gruetzner, T. R. Nelson, A. McPherson, S. M. Cameron, and A. C. Bernstein, “Propagation of self-focusing laser pulses in atmosphere: experiment versus numerical simulation,” J. Opt. Soc. Am. B 21, 2008 (2004).
  * Zahedpour, Wahlstrand, and Milchberg [2015] S. Zahedpour, J. K. Wahlstrand, and H. M. Milchberg, “Measurement of the nonlinear refractive index of air constituents at mid-infrared wavelengths,” Optics Letters 40, 5794 (2015).
  * Kolesik, Moloney, and Mlejnek [2002] M. Kolesik, J. V. Moloney, and M. Mlejnek, “Unidirectional optical pulse propagation equation,” Phys. Rev. Lett. 89, 283902 (2002).
  * NIS [2023] “NIST Atomic Spectra Database,” (2023).
  * Talebpour, Yang, and Chin [1999] A. Talebpour, J. Yang, and S. L. Chin, “Semi-empirical model for the rate of tunnel ionization in N2 and O2 molecule in an intense Ti:Saphir laser pulse,” Opt. Commun. 163, 29 (1999).
  * 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,” Optics Express 15, pp. 4577 – 4584 (2007).
  * Babushkin  _et al._ [2011] I. Babushkin, S. Skupin, A. Husakou, C. Köhler, E. Cabrera-Granado, L. Bergé, and J. Herrman, “Tailoring terahertz radiation by controling tunnel photoionization events in gases,” New J. Phys. 13, 123029 (2011).
  * Rasmussen  _et al._ [2023] M. Rasmussen, O. Nagy, S. Skupin, A. Stathopulos, L. Bergé, P. U. Jepsen, and B. Zhou, “Frequency-resolved characterization of broadband two-color air-plasma terahertz beam profiles,” Optics Express 31, pp. 9287 – 9298 (2023).
  * Zhong, Karpowicz, and Zhang [2006] H. Zhong, N. Karpowicz, and X.-C. Zhang, “Terahertz emission profile from laser-induced air plasma,” Appl. Phys. Lett. 88, 261103 (2006).
  * Jackson [1975] J. D. Jackson,  _Classical Electrodynamics_ (Wiley, New York, 1975).
  * Nguyen  _et al._ [2018] A. Nguyen, P. González de Alaiza Martínez, I. Thiele, S. Skupin, and L. Bergé, “THz field engineering in two-color femtosecond filaments using chirped and delayed laser pulses,” New. J. Phys. 20, 033026 (2018).
  * Tailliez  _et al._ [2020b] 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,” njp 22, 103038 (2020b).
  * Brunel [1990] F. Brunel, “Harmonic generation due to plasma effects in a gas undergoing multiphoton ionization in the high-intensity limit,” josab 7, 521 (1990).
  * Kim [2009] K.-Y. Kim, “Generation of coherent terahertz radiation in ultrafast laser-gas interactions,” pop 16, 056706 (2009).
  * Bergé  _et al._ [2013] L. Bergé, S. Skupin, C. Köhler, I. Babushkin, and J. Herrmann, “3D Numerical Simulations of THz Generation by Two-Color Laser Filaments,” Phys. Rev. Lett. 110, 073901 (2013).
  * Akturk  _et al._ [2007] S. Akturk, C. D’Amico, M. Franco, A. Couairon, and A. Mysyrowicz, “A simple method for determination of nonlinear propagation regimes in gases,” Opt. Express 15, pp. 15260 – 15267 (2007).
  * Born and Wolf [1959] M. Born and E. Wolf,  _Principles of Optics_ (Pergamon, New York, 1959).
  * Ionin  _et al._ [2013] A. A. Ionin, N. G. Iroshnikov, O. G. Kosareva, A. V. Larichev, D. V. Mokrousova, N. A. Panov, L. V. Seleznev, D. V. Sinitsyn, and E. S. Sunchugasheva, “Filamentation of femtosecond laser pulses governed by variable wavefront distortions via a deformable mirror,” J. Opt. Soc. Am. B 30 (2013).
  * Apeksimov  _et al._ [2018] D. V. Apeksimov, Y. E. Geints, A. A. Zemlyanov, A. M. Kabanov, V. K. Oshlanov, A. V. Petrov, and G. G. Matvienko, “Controlling tw-laser pulse long-range filamentation in air by a deformable mirror,” Applied Optics 57 (2018).
  * Mlejnek  _et al._ [1999] M. Mlejnek, M. Kolesik, J. V. Moloney, and E. M. Wright, “Optically turbulent femtosecond light guide in air,” Phys. Rev. Lett. 83, 2938 (1999).
  * Nikolaeva  _et al._ [2022] I. A. Nikolaeva, D. E. Shipilo, N. A. Panov, W. Liu, A. B. Savel’ev, and O. G. Kosareva, “Scaling Law of THz Yield from Two-Color Femtosecond Filament for Fixed Pump Power,” Photonics 9, 974 (2022).
  * D’Amico  _et al._ [2007] C. D’Amico, A. Houard, M. Franco, B. Prade, and A. Mysyrowicz, “Conical forward THz emission from femtosecond-laser-beam filamentation in air,” Phys. Rev. Lett. 98, 235002 (2007).
  * D’Amico  _et al._ [2008] C. D’Amico, A. Houard, S. Akturk, Y. Liu, J. Le Bloas, M. Franco, B. Prade, A. Couairon, V. T. Tikhonchuk, and A. Mysyrowicz, “Forward THz radiation emission by femtosecond filamentation in gases: theory and experiment,” New J. Phys. 10, 013015 (2008).
  * Balakin  _et al._ [2010] A. V. Balakin, A. V. Borodin, I. A. Kotelnikov, and A. P. Shkurinov, “Terahertz emission from a femtosecond laser focus in a two-color scheme,” J. Opt. Soc. Am. B 27, pp. 16 – 26 (2010).
  * Peñano  _et al._ [2010] J. R. Peñano, P. Sprangle, B. Hafizi, D. Gordon, and P. Serafim, “Terahertz generation in plasmas using two-color laser pulses,” Phys. Rev. E 81, 026407 (2010).
  * Clerici  _et al._ [2013] M. Clerici, M. Peccianti, B. E. Schmidt, L. Caspiani, M. Shalaby, M. Giguère, A. Lotti, A. Couairon, F. Légaré, T. Ozaki, D. Faccio, and R. Morandotti, “Wavelength scaling of terahertz generation by gas ionization,” Phys. Rev. Lett. 110, 253901 (2013).
  * Köhler  _et al._ [2011] C. Köhler, E. Cabrera-Granado, I. Babushkin, L. Bergé, J. Herrmann, and S. Skupin, “Directionality of terahertz emission from photoinduced gas plasmas,” Opt.Lett. 36, 3166 (2011).
  * Déchard  _et al._ [2017] J. Déchard, A. Nguyen, P. González de Alaiza Martínez, I. Thiele, S. Skupin, and L. Bergé, “Validity of the unidirectional propagation model: application to laser-driven terahertz emission,” J. Phys. Commun. 1, 055009 (2017).
  * Gradshteyn and Ryzhik [2014] I. S. Gradshteyn and I. M. Ryzhik,  _Table of Integrals, Series, and Products, Eighth Edition_ (Academic Press, San Diego, 2014).