# Modeling terahertz emissions from energetic electrons and ions in foil
targets irradiated by ultraintense femtosecond laser pulses

E. Denoual<EMAIL_ADDRESS>CEA, DAM, DIF, F-91297 Arpajon, France
Université Paris-Saclay, CEA, LMCE, F-91680 Bruyères-le-Châtel, France  L.
Bergé<EMAIL_ADDRESS>CEA, DAM, DIF, F-91297 Arpajon, France  Université
Paris-Saclay, CEA, LMCE, F-91680 Bruyères-le-Châtel, France  Centre des Lasers
Intenses et Applications, Université de Bordeaux-CNRS-CEA, F-33405 Talence
Cedex, France  X. Davoine<EMAIL_ADDRESS>CEA, DAM, DIF, F-91297
Arpajon, France  Université Paris-Saclay, CEA, LMCE, F-91680 Bruyères-le-
Châtel, France  L. Gremillet<EMAIL_ADDRESS>CEA, DAM, DIF, F-91297
Arpajon, France  Université Paris-Saclay, CEA, LMCE, F-91680 Bruyères-le-
Châtel, France

(August 28, 2024)

###### Abstract

Terahertz (THz) emissions from fast electron and ion currents driven in
relativistic, femtosecond laser-foil interactions are examined theoretically.
We first consider the radiation from the energetic electrons exiting the
backside of the target. Our kinetic model takes account of the coherent
transition radiation due to these electrons crossing the plasma-vacuum
interface as well as of the synchrotron radiation due to their deflection and
deceleration in the sheath field they set up in vacuum. After showing that
both mechanisms tend to largely compensate each other when all the electrons
are pulled back into the target, we investigate the scaling of the net
radiation with the sheath field strength. We then demonstrate the sensitivity
of this radiation to a percent-level fraction of escaping electrons. We also
study the influence of the target thickness and laser focusing. The same
sheath field that confines most of the fast electrons around the target
rapidly sets into motion the surface ions. We describe the THz emission from
these accelerated ions and their accompanying hot electrons by means of a
plasma expansion model that allows for finite foil size and multidimensional
effects. Again, we explore the dependencies of this radiation mechanism on the
laser-target parameters. Under conditions typical of current ultrashort laser-
solid experiments, we find that the THz radiation from the expanding plasma is
much less energetic – by one to three orders of magnitude – than that due to
the early-time motion of the fast electrons.

##  I Introduction

Intense sources of terahertz (THz) radiation are drawing growing interest as
their oscillation period, of the order of a few picoseconds, makes them
ideally suited for the study of numerous phenomena evolving on similar time
scales [1]. Their main direct applications include medical imaging [2],
molecular spectroscopy [3, 4], tomography [5], and modification of condensed
matter properties [6, 7, 8], to cite only a few. While intense lasers offer
promising prospects for developing compact, ultrashort THz sources, the main
challenge nowadays is to produce broadband THz pulses with $\rm mJ$-level
energies as is required for various uses [9, 10, 8, 11, 12]. This is a
nontrivial task as the most widely explored THz generation mechanisms, namely,
optical rectification in asymmetric crystals [13, 14, 15] or photoionization
of gases by two-color, moderate-intensity ($\lesssim 10^{15}\,\rm
W\,cm^{-2}$), femtosecond laser pulses [16, 17, 18], are to date limited to
tens of $\rm\mu J$ THz pulse energies and $\sim 1\,\rm GV\,m^{-1}$ field
strengths.

A more auspicious approach is to irradiate gaseous targets at relativistic
laser intensities ($I_{L}>10^{18}\,\rm W\,cm^{-2}$). In this regime, it has
been demonstrated that coherent transition radiation (CTR) from wakefield-
accelerated relativistic electron bunches at the rear plasma boundary can lead
to intense THz emissions, characterized by a few $100\,\rm\mu J$ energy yield
and $>10\,\rm GVm^{-1}$ field strength [19, 20]. Such a radiation is coherent
because the typical dimensions of the electron bunches ($\sim 1-5\,\rm\mu m$)
are smaller than the THz radiation wavelengths ($>10-100\,\rm\mu m$).
Consequently, the THz pulse energy essentially scales as the square number of
fast electrons, which makes it a potentially very efficient mechanism [21].

CTR also operates in relativistic laser-solid interactions, whereby, compared
to gas targets, it benefits from a stronger absorption of the laser energy
into MeV-range electrons, and hence from an increased number of radiating
particles [22, 23, 24, 25, 26, 27, 28, 29, 30]. However, because different
acceleration mechanisms are at play [31], these energetic electrons are
generally characterized by a much larger ($\sim 100\times$) angular divergence
than those generated by laser wakefields, which translates into a broader CTR
emission cone. Yet, owing to its high density ($\sim 10^{19-21}\,\rm
cm^{-3}$), the hot-electron population does not only radiate via CTR when
exiting a solid foil.

The latter mechanism indeed assumes that the fast electrons propagate
ballistically across the plasma-vacuum interface whilst most of them actually
get reflected in the strong charge-separation field that they set up in vacuum
[32, 33, 34]. This results in an additional coherent, synchrotron-type
radiation (CSR) of polarity opposite to that of CTR [28, 29]. An additional
complication follows from the fraction of fast electrons that are able to
escape the target, and thus just emit a single burst of CTR. The net THz
radiation resulting from those combined processes, CTR and CSR, will be
referred to as CTSR [Fig. 1(top)].

The sheath electric field induced by the hot electrons on both sides of the
target subsequently sets into motion the surface ions, a process widely known
as target normal sheath acceleration (TNSA) in the context of relativistic
laser-plasma interactions [35, 36, 37]. Because of their highest charge-to-
mass ratio, the protons, generally present in the form of contaminants, react
the fastest to that field, reaching velocities $\sim 0.1c$ ($c$ is the speed
of light) on $\sim 1\,\rm ps$ timescales. The resultant expanding plasma at
the target backside, composed of accelerated ions and electrostatically
trapped hot electrons, comprises two charge-separation regions: one negatively
charged at its outer boundary, and one positively charged around its inner
boundary [36]. Their time-varying properties lead to a dipole-type, low-
frequency radiation [38, 39, 26, 27], henceforth labeled plasma expansion
radiation (PER) [Fig. 1(bottom)].

Here, we develop two models to estimate the THz radiated spectra and energy
yields from the two aforementioned mechanisms, CTSR and PER, based,
respectively, on the fast-electron dynamics alone at the unperturbed target
backside and, on slower time scales, fast-electron-induced ion acceleration
into vacuum. Compared to previous related works [21, 27, 28, 29], we propose a
unified kinetic treatment of the destructively interfering CTR and CSR as
resulting from the beam electrons’ trajectories in vacuum. Moreover, we model
PER using a refined description of ion acceleration, notably allowing for the
time-decreasing surface charge density in the expanding sheath and for the
hot-electron cooling in thin foil targets. We expect our modeling to be mostly
valid in the case of micrometer-range foil targets driven by relativistic
femtosecond laser pulses. Our main prediction is that, under such conditions,
the energy radiated by the sole fast electrons via CTR and CSR should exceed
that due to plasma expansion by at least tenfold, and more in the likely case
of a percent-level fraction of escaping electrons. Our work therefore settles
a lingering debate on the dominant THz radiation process in ultrashort-pulse
laser experiments [38, 39, 25, 24, 26, 27, 29].

Figure 1: The two major THz radiation mechanisms in relativistic laser-foil
interactions considered in this work. (top) First, coherent transition and
synchrotron radiations are generated by the laser-accelerated electrons
crossing the target backside and being reflected (or not) in the sheath field
they set up in vacuum. (bottom) Subsequently, the nonneutral layers at the
edges of the rear-side expanding plasma emit a dipole-type radiation.

We start this paper by outlining the framework of our study in Sec. II.1. In
Sec. II.2, we recall the basics of far-field radiation from a charged particle
moving in vacuum near a perfect conductor, using the method of the image
charge. Section II.3 then characterizes the radiation from a single electron
exiting a perfect conductor and experiencing a constant electric field in
vacuum. In Sec. II.4, this problem is generalized to the case of an energy-
angle distributed electron beam originating from the laser-irradiated side of
the target. The respective integral expressions of CTR and CSR are then
detailed. In Sec. II.5, we specify the various parameters of the model, of
relevance to femtosecond laser-foil interactions.

Section III.1 presents the main spectral features of CTR, CSR and CTSR in a
typical ultrashort laser-foil configuration. A major finding is that CTR and
CSR closely compensate each other in the THz domain, when all fast electrons
are made to reflect back into the target. The scaling of the CTSR yield with
the sheath field strength is examined in Sec. III.2, while the possibly
dominant contribution of even a small ($\sim 1\,\%$) fraction of escaping
electrons is discussed in Sec. III.3. The sensitivity of the net THz radiation
to target thickness and laser focusing is addressed in Secs. III.4 and III.5.

Section IV next considers the radiation arising from the plasma expansion. The
general formalism of our approach is presented in Sec. IV.1. After detailing
the underlying model of TNSA in Sec. IV.2, we derive the integral form of the
PER energy spectra in Sec. IV.3. The dependencies of PER on the system’s
parameters are investigated in Secs. IV.4 and IV.5. Finally, Sec. V gathers
our conclusions.

##  II Coherent transition and synchrotron radiations from fast electrons

###  II.1 Framework of the study

The system investigated consists of a thin ($d\sim 1\,\rm\mu m$) solid foil
impacted by an ultraintense ($I_{L}\sim 10^{20}\,\rm W\,cm^{-2}$) and
ultrashort ($\tau_{L}\sim 30\leavevmode\nobreak\ \rm fs$) laser pulse of
wavelength $\lambda_{L}\sim 1\leavevmode\nobreak\ \rm\mu m$, focused to a
few-$\mu\rm m$ spot size ($w_{L}$). We suppose that the target is fully
ionized within a few optical cycles, hence turning into a plasma of
overcritical electron density $n_{e0}\gg n_{c}$ (where
$n_{c}=4\pi^{2}\epsilon_{0}m_{e}c^{2}/e^{2}\lambda_{L}^{2}$ is the critical
density, $m_{e}$ the electron mass, $e$ the elementary charge and
$\epsilon_{0}$ the vacuum permittivity), treated hereafter as a perfect
conductor with sharp boundaries.

The relativistic electrons driven by the laser pulse at the target front side
form a bunch of typical length $c\tau_{L}\sim 10\,\rm\mu m$, width $w_{L}$,
number density $n_{h}\sim n_{c}$ and average Lorentz factor
$\langle\gamma_{h}\rangle$. They are assumed to propagate ballistically on
their first pass through the material. Upon crossing the backside of the
target, they undergo an abrupt change in permittivity which causes transition
radiation [40, 41]. Wavelengths larger than the bunch size, i.e., lying in the
THz domain, are emitted coherently [42, 21, 43].

When exiting the target, the hot electrons induce a strong electric sheath
field, parallel to the surface normal and of typical strength
$E_{0}\sim(m_{e}c\omega_{0}/e)\sqrt{\langle\gamma_{h}\rangle(n_{h}/n_{c})}\sim
10^{12-13}\,\rm V\,m^{-1}$ [32, 36, 44], which pulls the vast majority of the
accelerated electrons back into the foil [34]. The synchrotron emission that
they generate while experiencing the sheath field is another source of
coherent THz radiation. One can anticipate, as will be examined in detail
below, that this deceleration-induced radiation will interfere destructively
with CTR, hereafter interpreted as resulting from the apparent sudden
acceleration of the electrons at the conductor’s surface [41].

The problem of CTR from laser-generated electron beams crossing at constant
velocity a plasma-vacuum boundary has already been widely addressed
theoretically in multiple frequency domains [42, 21, 43]. However, with the
exception of Refs. 28, 29, the simultaneous modeling of CTR and CSR due to the
non-ballistic motion of beam electrons in vacuum has received little attention
so far. As we will see further on, when almost all of these electrons are made
to reflux into the target, CTR and CSR produce THz fields with very similar
spectra, yet of opposite polarity, which therefore tend to cancel each other
out. The net THz radiation is thus determined by an integral over the electron
trajectories in vacuum, which depends on the initial electron energy and
propagation angle.

This radiation, as will be shown in Sec. III.3, is also highly sensitive to
the non-compensated CTR from the fraction of high-energy electrons able to
escape the target permanently. These, believed to make up at most a few
percent of the whole hot-electron population [34], can escape the target at
near the speed of light and thus contribute to a single uncompensated CTR
flash [28, 29]. The interplay of those mechanisms, and their variations with
the laser-target parameters, will be thoroughly examined in Secs. III.4 and
III.5.

Before proceeding, two limitations of our modeling are already worth
mentioning. First, we will consider only the first excursion of the fast
electrons into vacuum, and thus neglect their subsequent THz emissions while
they bounce back and forth across the target [23, 25, 30]. Accordingly, the
following estimates of the energy radiated from the sole fast electrons should
be considered as lower values. Second, by assuming a target of infinite
transverse size, we will discard the multiple THz emissions that are expected
to arise from the laser axis and the target edges when the fast electrons
recirculate transversely across a finite-width target [45, 30].

###  II.2 Radiation from an electron accelerated in vacuum near a perfect
conductor

The electromagnetic field radiated in vacuum by an accelerated charged
particle can be split into two components [46, 47]: a velocity field that
rapidly decays in space as $R^{-2}$ and an acceleration field that decays as
$R^{-1}$, where $R=\left|\bm{r}-\bm{r}_{p}(t)\right|$ is the distance from the
particle, located at $\bm{r}_{p}(t)$, to the detector, located at $\bm{r}$. In
the far-field limit, the detected field reduces to the acceleration field [46,
47],

$\bm{E}_{\rm acc}(\bm{r},t)=\frac{q}{4\pi
c\varepsilon_{0}}\left[\frac{\bm{\hat{R}}\times\big{\\{}(\bm{\hat{R}}-\bm{\beta}_{p})\times\dot{\bm{\beta}}_{p}\big{\\}}}{R\big{(}1-\bm{\beta}_{p}\cdot\bm{\hat{R}}\big{)}^{3}}\right]_{t_{\rm
ret}}\,,$ (1)

where $\bm{\beta}_{p}(t)=\dot{\bm{r}}_{p}(t)/c$ is the normalized velocity of
the particle, $\dot{\bm{\beta}}_{p}(t)$ its normalized acceleration, $q$ its
charge, $\varepsilon_{0}$ the permittivity of free space and
$\bm{\hat{R}}=\left(\bm{r}-\bm{r}_{p}(t)\right)/R$ the unit direction of
observation. The operator $\left[\,.\,\right]_{t_{\rm ret}}$ evaluates its
argument at the retarded time $t_{\rm ret}$, implicitly defined by $t_{\rm
ret}=t-R(t_{\rm ret})/c$.

When the particle moves in vacuum in the vicinity of a conductor, Eq. (1) must
be corrected so that the tangential component of the electric field vanishes
at the conductor surface (assumed planar). Microscopically, this arises
because of the additional field generated by polarization surface currents
[48]. In the perfect conductor approximation, the field induced in vacuum by
those currents is identical to that from a virtual particle of charge $-q$ and
trajectory symmetric to that of the real particle with respect to the
conductor surface[41, 49] (see Fig. 2). The total radiation from the particle
in the presence of the conductor is therefore

$\bm{E}_{\rm rad}(\bm{r},t)=\bm{E}_{\rm acc}^{-}(\bm{r},t)+\bm{E}_{\rm
acc}^{+}(\bm{r},t)\,,$ (2)

where $\bm{E}_{\rm acc}^{-}(\bm{r},t)$ and $\bm{E}_{\rm acc}^{+}(\bm{r},t)$
represent the acceleration fields generated by, respectively, the real
(superscript -) and image (superscript +) particles. Both fields are evaluated
at their own retarded times, $t_{\rm ret}^{-}$ or $t_{\rm ret}^{+}$, depending
on the positions of the particle and observer.

The power radiated by the particle and its image charge per unit solid angle
reads [47, 46]

$\frac{dP_{\rm
rad}(\bm{\hat{r}},t)}{d\Omega}=c\varepsilon_{0}\left|\left[R^{-}\bm{E}_{\rm
acc}^{-}\right]_{t_{\rm ret}^{-}}+\left[R^{+}\bm{E}_{\rm
acc}^{+}\right]_{t_{\rm ret}^{+}}\right|^{2}\,,$ (3)

where we introduce the direction of observation
$\bm{\hat{r}}=(\sin\theta\cos\Psi,\sin\theta\sin\Psi,\cos\theta)=\bm{r}/r$,
the polar angle $\theta$, the azimuthal angle $\Psi$, and
$\mathrm{d}\Omega=\sin\theta\mathrm{d}\theta\mathrm{d}\Psi$. The total energy
radiated per unit solid angle can be expressed as [46, 47]:

$\frac{\partial\mathcal{E}_{\rm
rad}(\bm{\hat{r}})}{\partial\Omega}=\int_{-\infty}^{\infty}\frac{\partial
P_{\rm
rad}(\bm{\hat{r}},t)}{\partial\Omega}\mathrm{d}t\equiv\int_{0}^{\infty}\frac{\partial^{2}\mathcal{I}_{\rm
rad}(\bm{\hat{r}},\nu)}{\partial\nu\partial\Omega}\mathrm{d}\nu\,,$ (4)

where the second equality defines $\partial^{2}\mathcal{I}_{\rm
rad}/\partial\nu\partial\Omega$, the energy radiated per unit solid angle and
frequency ($\nu$) interval. Upon combining Eqs. (1)-(4), performing the change
of variable $t\mapsto t_{\rm ret}$ and assuming that the accelerated particle
remain far away from the observer, i.e., $R(t)\simeq
r-\bm{\hat{r}}\cdot\bm{r}_{p}(t)$, one obtains [46, 47]

$\frac{\partial^{2}\mathcal{I}_{\rm
rad}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)=\frac{q^{2}}{8\pi^{2}\varepsilon_{0}c}\left|\int\frac{\bm{\hat{r}}\times[(\bm{\hat{r}}-\bm{\beta}_{p}^{-})\times\dot{\bm{\beta}}_{p}^{-}]}{(1-\bm{\beta}_{p}^{-}\cdot\bm{\hat{r}})^{2}}e^{-2i\pi\nu(t-\bm{\hat{r}}\cdot\bm{r}_{p}^{-}(t)/c)}-\frac{\bm{\hat{r}}\times[(\bm{\hat{r}}-\bm{\beta}_{p}^{+})\times\dot{\bm{\beta}}_{p}^{+}]}{(1-\bm{\beta}_{p}^{+}\cdot\bm{\hat{r}})^{2}}e^{-2i\pi\nu(t-\bm{\hat{r}}\cdot\bm{r}_{p}^{+}(t)/c)}\,\mathrm{d}t\right|^{2}\,.$
(5)

This expression describes the full energy spectrum radiated by a charged
particle moving in vacuum in the vicinity of a perfect conductor. In the
following section, we will decompose this radiation into two components,
namely, transition radiation and synchrotron radiation.

###  II.3 Transition and synchrotron radiations from an electron exiting and
coming back into a perfect conductor

The generic scenario addressed in our study is sketched in Fig. 2. It consists
of an electron (i) laser-accelerated at the front side of a foil target, (ii)
traveling ballistically through it, (iii) crossing its rear side, (iv) being
reflected in vacuum by the sheath field established by the hot electrons and
(v) re-crossing the target’s rear surface. Because of plasma shielding, stages
(i) and (ii) of the electron motion do not lead to outgoing radiation from the
backside of the target. By contrast, stages (iii) and (v), corresponding to
apparent sudden accelerations of the electron and its image charge, give rise
to two consecutive flashes of transition radiation [40, 41], at the exit
($t=t_{\mathrm{e}}$) and return ($t=t_{\mathrm{r}}$) times, while stage (iv)
generates synchrotron-type radiation over the time interval
$t_{\mathrm{e}}<t<t_{\mathrm{r}}$.

To identify the transition- and synchrotron-type components in Eq. (5), where
the time integral is performed over the particle trajectory in vacuum, it is
convenient to express the normalized velocities ($\bm{\beta}_{p}^{\pm}$) of
the electron and its image charge as

$\bm{\beta}_{p}^{\pm}(t)=\begin{cases}\bm{\beta}_{\mathrm{e}}^{\pm}\left[\frac{t-t_{\mathrm{e}}}{\delta
t}+1\right]&\,t_{\mathrm{e}}-\delta t<t\leq t_{\mathrm{e}}\,,\\\
\bm{\beta}^{\pm}(t)&\,t_{\mathrm{e}}<t<t_{\mathrm{r}}\,,\\\
\bm{\beta}_{\mathrm{r}}^{\pm}\left[\frac{t_{\mathrm{r}}-t}{\delta
t}+1\right]&\,t_{\mathrm{r}}\leq t<t_{\mathrm{r}}+\delta t\,,\\\
0&\,\text{otherwise}\,,\end{cases}\vskip 6.0pt plus 2.0pt minus 2.0pt$ (6)

where $\bm{\beta}_{\mathrm{e}}^{\pm}\equiv\bm{\beta}^{\pm}(t_{\mathrm{e}})$
and $\bm{\beta}_{\mathrm{r}}^{\pm}\equiv\bm{\beta}^{\pm}(t_{\mathrm{r}})$
denote the velocities at the exit and return times, respectively, and $\delta
t$ is an infinitesimal time interval.

Figure 2: Transition and synchrotron radiations from a relativistic electron
confined at the target-vacuum interface. The electron trajectory starts at the
left (laser-irradiated) side of the target, at time $t_{\mathrm{f}}$ and
position $\bm{r}_{\mathrm{f}}$. The initial electron momentum is
$\bm{p}_{\mathrm{e}}=p_{\mathrm{e}}(\sin\psi\cos\varphi,\sin\psi\sin\varphi,\cos\psi)$
where $\psi$ and $\varphi$ are the polar and azimuthal angles, respectively.
After traveling ballistically through the target (red dashed curve), the
electron crosses its rear side at time $t_{\mathrm{e}}$ and position
$\bm{r}_{\mathrm{e}}$. In vacuum, it is then reflected by the longitudinal
sheath field $\bm{E}_{0}$, assumed homogeneous and stationary, and re-enters
the target at time $t_{\mathrm{r}}$ and position $\bm{r}_{\mathrm{r}}$ (blue
solid curve). The color gradient represents the electrostatic potential
profile $\Phi(z)$, see Sec. II.3. The trajectory of the associated image
charge is plotted as a grey dotted curve. The observation direction is
$\bm{\hat{r}}=(\sin\theta\cos\Psi,\sin\theta\sin\Psi,\cos\theta)=\bm{r}/r$,
where $\theta$ and $\Psi$ are the corresponding polar and azimuthal angles. In
this 2D representation, for simplicity, the observer is taken to lie in the
plane of particle motion (i.e. $\Psi=\varphi$).

Noting that [46]

$\frac{\bm{\hat{r}}\times[(\bm{\hat{r}}-\bm{\beta}_{p}^{\pm})\times\dot{\bm{\beta}}_{p}^{\pm}]}{\left(1-\bm{\beta}_{p}^{\pm}\cdot\bm{\hat{r}}\right)^{2}}=\frac{\mathrm{d}}{\mathrm{d}t}\left[\frac{\bm{\hat{r}}\times(\bm{\hat{r}}\times\bm{\beta}_{p}^{\pm})}{1-\bm{\beta}_{p}^{\pm}\cdot\bm{\hat{r}}}\right]\,,$
(7)

we can perform an integration by parts in Eq. (5) and simplify the resulting
expression using the vector identity
$[\bm{\hat{r}}\times(\bm{\hat{r}}\times\bm{a})]\cdot[\bm{\hat{r}}\times(\bm{\hat{r}}\times\bm{b})]=(\bm{\hat{r}}\times\bm{a})\cdot(\bm{\hat{r}}\times\bm{b})$.
We then obtain

$\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
rad}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)$
$\displaystyle=\frac{q^{2}}{8\pi^{2}\varepsilon_{0}c}\Bigg{|}2i\pi\nu\int_{t_{\mathrm{e}}-\delta
t}^{t_{\mathrm{r}}+\delta
t}\bm{\hat{r}}\times\left[\bm{\beta}_{p}^{-}(t)e^{-i\Theta^{-}(t)}-\bm{\beta}_{p}^{+}(t)e^{-i\Theta^{+}(t)}\right]\mathrm{d}t$
$\displaystyle+\left[\frac{\bm{\hat{r}}\times\bm{\beta}_{p}^{-}(t)}{1-\bm{\beta}_{p}^{-}(t)\cdot\bm{\hat{r}}}e^{-i\Theta^{-}(t)}-\frac{\bm{\hat{r}}\times\bm{\beta}_{p}^{+}(t)}{1-\bm{\beta}_{p}^{+}(t)\cdot\bm{\hat{r}}}e^{-i\Theta^{+}(t)}\right]_{t_{\mathrm{e}-\delta
t}}^{t_{\mathrm{r}+\delta t}}\Bigg{|}^{2}\,,$ (8)

where $\Theta^{\pm}(t)=2\pi\nu[t-\bm{\hat{r}}\cdot\bm{r}_{p}^{\pm}(t)/c]$.
Since, from Eqs. (6), $\bm{\beta}_{p}(t_{\mathrm{r}}+\delta
t)=\bm{\beta}_{p}(t_{\mathrm{e}}-\delta t)=0$, the boundary value terms vanish
and the full intensity distribution radiated by the electron is therefore
given by

$\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
TSR}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)=\frac{q^{2}\nu^{2}}{2\varepsilon_{0}c}$
$\displaystyle\times\left|\int_{t_{\mathrm{e}}}^{t_{\mathrm{r}}}\bm{\hat{r}}\times\left[\bm{\beta}_{p}^{-}(t)e^{-i\Theta^{-}(t)}-\bm{\beta}_{p}^{+}(t)e^{-i\Theta^{+}(t)}\right]\,\mathrm{d}t\right|^{2}\,,$
(9)

where the acronym TSR stands for transition and synchrotron radiation. Its
coherent version (CTSR) in the case of a compact electron bunch will be
addressed below.

The intensity distribution of transition radiation,
$\partial^{2}\mathcal{I}_{\rm TR}/\partial\nu\partial\Omega$, is obtained from
Eq. (5) by integrating by parts over the intervals $t_{\mathrm{e}}-\delta
t<t<t_{\mathrm{e}}$ and $t_{\mathrm{r}}<t<t_{\mathrm{r}}+\delta t$ and taking
the limit $\delta t\to 0$. Only the boundary value term then remains, leading
to a variant of the well-known Ginzburg formula [41]

$\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
TR}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)=\frac{q^{2}}{8\pi^{2}\varepsilon_{0}c}$
$\displaystyle\times\left|\left[\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{e}}^{-}}{1-\bm{\beta}_{\mathrm{e}}^{-}\cdot\bm{\hat{r}}}-\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{e}}^{+}}{1-\bm{\beta}_{\mathrm{e}}^{+}\cdot\bm{\hat{r}}}\right]e^{-i\Theta_{\mathrm{e}}}\right.$
$\displaystyle\left.-\left[\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{r}}^{-}}{1-\bm{\beta}_{\mathrm{r}}^{-}\cdot\bm{\hat{r}}}-\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{r}}^{+}}{1-\bm{\beta}_{\mathrm{r}}^{+}\cdot\bm{\hat{r}}}\right]e^{-i\Theta_{\mathrm{r}}}\right|^{2}\,,$
(10)

where $\Theta_{\mathrm{e},\mathrm{r}}\equiv
2\pi\nu(t_{\mathrm{e},\mathrm{r}}-\bm{\hat{r}}\cdot\bm{r}_{\mathrm{e},\mathrm{r}}/c)$.

The synchrotron component of the radiation is obtained by integrating Eq. (5)
by parts over the time interval $t_{\mathrm{e}}<t<t_{\mathrm{r}}$:

$\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
SR}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)=\frac{q^{2}}{8\pi^{2}\varepsilon_{0}c}$
$\displaystyle\times\Bigg{|}2i\pi\nu\int_{t_{\mathrm{e}}}^{t_{\mathrm{r}}}\bm{\hat{r}}\times\left[\bm{\beta}_{p}^{-}(t)e^{-i\Theta^{-}(t)}-\bm{\beta}_{p}^{+}(t)e^{-i\Theta^{+}(t)}\right]\,\mathrm{d}t$
$\displaystyle+\left[\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{r}}^{-}}{1-\bm{\beta}_{\mathrm{r}}^{-}\cdot\bm{\hat{r}}}-\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{r}}^{+}}{1-\bm{\beta}_{\mathrm{r}}^{+}\cdot\bm{\hat{r}}}\right]e^{-i\Theta_{\mathrm{r}}}$
$\displaystyle-\left[\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{e}}^{-}}{1-\bm{\beta}_{\mathrm{e}}^{-}\cdot\bm{\hat{r}}}-\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{e}}^{+}}{1-\bm{\beta}_{\mathrm{e}}^{+}\cdot\bm{\hat{r}}}\right]e^{-i\Theta_{\mathrm{e}}}\Bigg{|}^{2}\,,$
(11)

Of course, one recovers the full spectrum (II.3) by summing the terms within
the squared vertical bars in Eqs. (10) and (11). The above expression is
interesting in showing that whatever the evolution of $\bm{\beta}^{\pm}(t)$
during $t_{\mathrm{e}}<t<t_{\mathrm{r}}$, the synchrotron radiation generates
a spectrum increasingly resembling that due to transition radiation when
$t_{\mathrm{r}}-t_{\mathrm{e}}\equiv\Delta t_{r}\to 0$. This result will be
demonstrated numerically in the next sections.

In the following, for simplicity, the electric sheath field acting on the
electron will be assumed uniform, stationary and parallel to the target
surface normal, $\bm{E}_{0}=E_{0}\hat{\bm{z}}$ (with $E_{0}>0$). The electron
trajectory can be easily calculated using the Hamiltonian

$\mathcal{H}=m_{e}c^{2}\gamma(t)-e\Phi(z(t))\,,$ (12)

where $\gamma\equiv\sqrt{1+p^{2}/(m_{e}c)^{2}}$ is the electron Lorentz factor
and $\Phi=-E_{0}z$ the electrostatic potential.

The electron crosses the target backside ($z=0$) at time $t=t_{\mathrm{e}}$
and transverse position $\bm{r}_{\mathrm{e},\perp}$ (neglecting the $\delta t$
interval) with momentum
$\bm{p}_{\mathrm{e}}=p_{\mathrm{e}}(\sin\psi\cos\varphi,\sin\psi\sin\varphi,\cos\psi)$,
where $\psi$ and $\varphi$ correspond to the polar and azimuthal angles,
respectively. Since $\Phi$ only depends on $z$, $\mathcal{H}$ is a constant of
motion. Moreover,
$\gamma(t_{\mathrm{r}})=\gamma(t_{\mathrm{e}})\equiv\gamma_{\mathrm{e}}$ and
$\bm{p}_{\perp}(t)=\bm{p}_{\mathrm{e},\perp}$ (the subscript ⟂ refers to the
surface normal plane).

Introducing the normalized momentum
$\bm{u}\equiv\bm{p}/m_{e}c\equiv\gamma\bm{\beta}$, normalized position
$\bar{\bm{r}}\equiv eE_{0}\bm{r}/m_{e}c^{2}$, relative position
$\bar{\bm{\xi}}=\bar{\bm{r}}-\bar{\bm{r}}_{\mathrm{e},\perp}$, normalized time
$\bar{t}\equiv eE_{0}t/m_{e}c$ and relative time
$\bar{\tau}=\bar{t}-\bar{t}_{\mathrm{e}}$, one obtains

$\displaystyle\gamma(\bar{\tau})=\sqrt{\left(u_{\mathrm{e},z}-\bar{\tau}\,\right)^{2}+\gamma_{\perp}^{2}}\,,$
(13)
$\displaystyle\beta_{z}(\bar{\tau})=\left(u_{\mathrm{e},z}-\bar{\tau}\,\right)/\gamma(\bar{\tau})\,,$
(14)
$\displaystyle\bm{\beta}_{\perp}(\bar{\tau})=\bm{u}_{\mathrm{e},\perp}/\gamma(\bar{\tau})\,,$
(15)
$\displaystyle\bar{\xi}_{z}(\bar{\tau})=\gamma_{\mathrm{e}}-\gamma(\bar{\tau})\,,$
(16)
$\displaystyle\bar{\bm{\xi}}_{\perp}(\bar{\tau})=\frac{\bm{u}_{\mathrm{e},\perp}}{2}\ln\left(\frac{\left(\gamma_{\mathrm{e}}+u_{\mathrm{e},z}\right)\left[\gamma(\bar{\tau})-u_{\mathrm{e},z}+\bar{\tau}\right]}{\left(\gamma_{\mathrm{e}}-u_{\mathrm{e},z}\right)\left[\gamma(\bar{\tau})+u_{\mathrm{e},z}-\bar{\tau}\right]}\right)\,,$
(17)

for $0<\bar{\tau}<\overline{\Delta t}_{\mathrm{r}}$, where $\overline{\Delta
t}_{\mathrm{r}}\equiv\bar{t}_{\mathrm{r}}-\bar{t}_{\mathrm{e}}=2u_{\mathrm{e},z}$
is the (normalized) time spent by the electron in vacuum, and
$\gamma_{\perp}=\sqrt{1+u_{\mathrm{e},\perp}^{2}}$.

The full (TSR) and synchrotron-only (SR) radiation spectra are computed by
substituting the above expressions into Eq. (II.3) and Eq. (11), respectively.
Furthermore, given our choice of $\bm{E}_{0}$, we have the relations
$\bm{\beta}_{\mathrm{r}}^{\pm}=\bm{\beta}_{\mathrm{e}}^{\mp}$ and
$\beta_{\mathrm{r},z}^{\pm}=-\beta_{\mathrm{e},z}^{\pm}$, allowing the TR
spectrum (10) to be recast as

$\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
TR}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)=\frac{q^{2}}{2\pi^{2}\varepsilon_{0}c}\cos^{2}\left(\frac{\Theta_{\mathrm{e}}-\Theta_{\mathrm{r}}}{2}\right)$
$\displaystyle\times\left|\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{e}}^{-}}{1-\bm{\beta}_{\mathrm{e}}^{-}\cdot\bm{\hat{r}}}-\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{r}}^{-}}{1-\bm{\beta}_{\mathrm{r}}^{-}\cdot\bm{\hat{r}}}\right|^{2}\,.$
(18)

###  II.4 Full radiation from an energy-angle-distributed electron beam

The frequency-angle spectrum radiated by an electron bunch exiting a perfect
conductor and experiencing a uniform electric field can be readily evaluated
by summing Eq. (II.3) over $N_{h}\gg 1$ particles (labeled by the subscript
$l$) and taking an ensemble average:

$\frac{\partial^{2}\mathcal{I}_{\rm
CTSR}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)=\left\langle\frac{q^{2}\nu^{2}}{2\varepsilon_{0}c}\bigg{|}\sum_{l=1}^{N_{h}}\big{[}\bm{A}_{l}^{-}-\bm{A}_{l}^{+}\vphantom{\int}\big{]}(\bm{\hat{r}},\nu)\bigg{|}^{2}\right\rangle\,,$
(19)

where

$\bm{A}_{l}^{\pm}(\bm{\hat{r}},\nu)=\int_{t_{\mathrm{e},l}}^{t_{\mathrm{r},l}}\left[\bm{\hat{r}}\times\bm{\beta}_{p,l}^{\pm}(t)\right]e^{-i\Theta_{l}^{\pm}(t)}\,\mathrm{d}t\,.$
(20)

We recall that the acronym CTSR is for _coherent_ transition and synchrotron
radiation.

In the realistic case of a finite-size, energy-angle distributed electron
bunch, the full radiation pattern is determined by the relative phase shifts
between the fields emitted along the particle trajectories in vacuum. In the
following, each electron trajectory is taken to originate from the _front
side_ ($z=-d$) of the target, where (leaving out the subscript $l$) it is
parameterized by the injection time $t_{\mathrm{f}}$, initial transverse
position $\bm{r}_{\mathrm{f},\perp}$ and initial normalized momentum
$\bm{u}_{\mathrm{f}}=\bm{u}_{\mathrm{e}}=\bm{p}_{\mathrm{e}}/m_{e}c$ of the
electron. After traveling ballistically through the target, the electron
reaches its rear side ($z=0$) at time
$t_{\mathrm{e}}=t_{\mathrm{f}}+d/c\beta_{\mathrm{f},z}$ and transverse
position
$\bm{r}_{\mathrm{e},\perp}=\bm{r}_{\mathrm{f},\perp}+\bm{\beta}_{\mathrm{f},\perp}d/\beta_{\mathrm{f},z}$.

In order to pass to the continuous limit, we introduce $F(\bm{r},t,\bm{u})$,
the flux probability function characterizing the electron beam at the front
side of the target, normalized as

$\iiint
F(\bm{r}_{\mathrm{f},\perp},t_{\mathrm{f}},\bm{u}_{\mathrm{f}})\,\mathrm{d}^{2}\bm{r}_{\mathrm{f},\perp}\,\mathrm{d}t_{\mathrm{f}}\,\mathrm{d}^{3}\bm{u}_{\mathrm{f}}=1\,.$
(21)

This probability function is taken in the form of a separable function,
$F(\bm{r}_{\mathrm{f},\perp},t_{\mathrm{f}},\bm{u}_{\mathrm{f}})=f(\bm{r}_{\mathrm{f},\perp})h(t_{\mathrm{f}})g(\bm{u}_{\mathrm{f}})$,
where the distribution functions $f$, $h$ and $g$ will be specified in Sec.
II.5. The discrete sum in Eq. (19) can be replaced by integrals over time,
transverse position and momentum (see Appendix A):

$\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
CTSR}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)=\frac{q^{2}\nu^{2}}{2\varepsilon_{0}c}N_{h}^{2}\Big{|}\iiint
f(\bm{r}_{\mathrm{f},\perp})h(t_{\rm f})g(\bm{u}_{\mathrm{f}})$
$\displaystyle\times\big{[}\bm{A}^{-}(\bm{\hat{r}},\nu)-\bm{A}^{+}(\bm{\hat{r}},\nu)\big{]}\,\mathrm{d}^{2}\bm{r}_{\mathrm{f},\perp}\,\mathrm{d}t_{\mathrm{f}}\,\mathrm{d}^{3}\bm{u}_{\mathrm{f}}\Big{|}^{2}\,.$
(22)

This formula can further be simplified by rewriting the phase term in Eq. (20)
as

$\displaystyle\Theta^{\pm}(t)$
$\displaystyle=2\pi\nu\left[t-\bm{\hat{r}}\cdot\bm{r}^{\pm}(t)/c\right]$
$\displaystyle=2\pi\bar{\nu}\left[(\bar{\tau}+\bar{t}_{\mathrm{e}})-\bm{\hat{r}}\cdot(\bar{\bm{\xi}}^{\pm}(\bar{\tau})+\bar{\bm{r}}_{\mathrm{e},\perp})\right]$
$\displaystyle=\Theta_{\mathrm{e}}+2\pi\bar{\nu}\left[\bar{\tau}-\bm{\hat{r}}\cdot\bar{\bm{\xi}}^{\pm}(\bar{t})\right]\,,$
(23)

where $\bar{\nu}\equiv\nu m_{e}c/eE_{0}$ is the normalized frequency and where

$\displaystyle\Theta_{\mathrm{e}}$
$\displaystyle=2\pi\bar{\nu}\left[\bar{t}_{\mathrm{e}}-\bm{\hat{r}}\cdot\bar{\bm{r}}_{\mathrm{e},\perp}\right]$
$\displaystyle=2\pi\bar{\nu}\left(\bar{t}_{\mathrm{f}}-\bm{\hat{r}}\cdot\bar{\bm{r}}_{\mathrm{f}}\right)+2\pi\bar{\nu}\bar{d}\left(1-\bm{\hat{r}}\cdot\bm{\beta}_{\mathrm{f},\perp}\right)/\beta_{\mathrm{f},z}$
(24)

characterizes the propagation of the electron through the foil. Performing the
change of variable $t\mapsto\bar{\tau}$ and factoring out
$\Theta_{\mathrm{e}}$ from Eq. (20) yields

$\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
CTSR}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)=\frac{q^{2}\bar{\nu}^{2}}{2\varepsilon_{0}c}N_{h}^{2}\left|\int
g(\bm{u})F(\bm{\hat{r}},\nu,\bm{u})\right.$
$\displaystyle\times\left.\Big{[}\widebar{\bm{A}}^{-}(\bm{\hat{r}},\nu)-\widebar{\bm{A}}^{+}(\bm{\hat{r}},\nu)\vphantom{\int}\Big{]}\,\mathrm{d}^{3}\bm{u}\right|^{2}\,,$
(25)

with $F(\bm{\hat{r}},\nu,\bm{u})$ the form factor of the electron beam defined
by

$F(\bm{\hat{r}},\nu,\bm{u})=\iint f(\bm{r}_{\mathrm{f},\perp})h(t_{\rm
f})e^{-i\Theta_{\mathrm{e}}}\mathrm{d}t_{\rm
f}\,\mathrm{d}^{2}\bm{r}_{\mathrm{f},\perp}\,,$ (26)

and where

$\widebar{\bm{A}}^{\pm}(\bm{\hat{r}},\nu)=\int_{0}^{\overline{\Delta
t}_{\mathrm{r}}}\left[\bm{\hat{r}}\times\bm{\beta}^{\pm}(\bar{\tau})\right]e^{-2i\pi\bar{\nu}\left[\bar{\tau}-\bm{\hat{r}}\cdot\bar{\bm{\xi}}^{\pm}(\bar{\tau})\right]}\,\mathrm{d}\bar{\tau}\,.$
(27)

Equations (25)-(27), supplemented with Eqs. (13)-(17), give the energy-angle
spectrum of radiation from an electron beam exiting a perfect conductor and
being reflected by a stationary and homogeneous electric sheath field. This
spectrum depends not only on the distribution functions $f(\bm{r}_{\perp})$,
$h(t)$ and $g(\bm{u})$ characterizing the beam source at the target front
side, but also on the finite target thickness, which determines the
longitudinal and transverse spreading of the beam after its passage through
the foil.

As explained in Sec. II.3, the respective contributions of CTR
($\partial^{2}\mathcal{I}_{\rm CTR}/\partial\nu\partial\Omega$) and CSR
($\partial^{2}\mathcal{I}_{\rm CSR}/\partial\nu\partial\Omega$) to the total
spectrum can be distinguished in the time integral involved in
$\bar{\bm{A}}^{\pm}(\bm{r},\nu)$ [Eq. (27)], namely,

$\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
CTR}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)$
$\displaystyle=\frac{N_{h}^{2}q^{2}}{2\pi^{2}\varepsilon_{0}c}\left|\int
g(\bm{u})F(\theta,\nu,\bm{u})\bar{\bm{A}}_{\rm
TR}\left(\bm{\hat{r}},\nu\right)\,\mathrm{d}^{3}\bm{u}\right|^{2}\,,$ (28)
$\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
CSR}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)$
$\displaystyle=\frac{N_{h}^{2}q^{2}}{2\pi^{2}\varepsilon_{0}c}\left|\int
g(\bm{u})F(\theta,\nu,\bm{u})\Biggl{\\{}\vphantom{\int}i\pi\bar{\nu}\Bigl{[}\widebar{\bm{A}}^{-}(\bm{\hat{r}},\nu)-\widebar{\bm{A}}^{+}(\bm{\hat{r}},\nu)\Bigr{]}-\bar{\bm{A}}_{\rm
TR}\left(\bm{\hat{r}},\nu\right)\Biggr{\\}}\,\mathrm{d}^{3}\bm{u}\right|^{2}\,,$
(29)

where $\bar{\bm{A}}_{\rm TR}\left(\bm{\hat{r}},\nu\right)$ characterizes the
full transition radiation of a particle in the vicinity of the perfect
conductor [stages $\rm(iii)$ and $\rm(v)$ in Fig. 2] and is defined by:

$\bar{\bm{A}}_{\rm
TR}\left(\bm{\hat{r}},\nu\right)=\cos\left(\frac{\Theta_{\mathrm{e}}-\Theta_{\mathrm{r}}}{2}\right)\\\
\times\left[\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{e}}^{-}}{1-\bm{\beta}_{\mathrm{e}}^{-}\cdot\bm{\hat{r}}}-\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{r}}^{-}}{1-\bm{\beta}_{\mathrm{r}}^{-}\cdot\bm{\hat{r}}}\right]e^{-i\Theta_{\mathrm{e}}}\,,$
(30)

where
$\Theta_{\mathrm{r}}=\Theta_{\mathrm{e}}+2\pi\bar{\nu}\left[\widebar{\Delta
t_{\mathrm{r}}}-\bm{\hat{r}}\cdot\bar{\bm{\xi}}\left(\widebar{\Delta
t_{\mathrm{r}}}\right)\right]$.

The numerical evaluation of Eqs. (25), (28) and (29) is challenging because of
the four nested integrals involved, or even seven if one wishes to compute the
total radiated energy (integrated over observation angles $(\theta,\Psi)$ and
frequencies $\nu$). To simplify, we will assume hereafter that the system is
axisymmetric with respect to the target normal ($\bm{\hat{z}}$ axis). The
integration over $(\bm{\hat{r}},\nu)$ then reduces to an integration over
$(\theta,\nu)$ due to invariance of
$\partial^{2}\mathcal{I}/\partial\nu\partial\Omega$ over $\Psi$. In this work,
the calculation over the angle of observation, frequency and momentum is
parallelized on multiple GPUs, using a 2D kernel
$(\theta\times\nu,u\times\psi\times\varphi)$ to efficiently distribute the
workload. The innermost time integration [Eq. (27)] is performed over 1000
time steps through a type-2 nonuniform Fourier transform, as defined in Refs.
50, 51, 52. Numerical integrals over $(u,\psi,\varphi)$ are computed using
$512\times 128\times 64$ points for each $(\nu,\theta)$ pair. Radiated spectra
are discretized over $160\times 90$ points in $(\nu,\theta)$ space, amounting
to a total of $\sim 6\times 10^{10}$ nonuniform Fourier transforms. This makes
parametric scans quite computationally expensive despite the GPU
parallelization.

###  II.5 Model parameters

To solve the previous formulas, the distribution functions defining the
(axisymmetric) electron beam source need to be specified. The spatial and
temporal profiles of the beam are assumed to be Gaussian,

$\displaystyle f(r)$ $\displaystyle=\frac{4\ln 2}{\pi w_{L}^{2}}e^{-4\ln
2\,(r/w_{L})^{2}}\,,$ (31) $\displaystyle h(t)$
$\displaystyle=\frac{2\sqrt{\ln 2}}{\sqrt{\pi}\tau_{L}}e^{-4\ln
2\,(t/\tau_{L})^{2}}\,,$ (32)

where the full-width-at-half-maximum diameter ($w_{L}$) and duration
($\tau_{L}$) are taken equal to those characterizing the intensity profile of
the laser drive. The above profiles lead to the form factor

$\displaystyle F(\theta,\nu,\bm{u})=$ $\displaystyle\exp\Bigl{(}-(\pi^{2}/4\ln
2)\nu^{2}\left[\tau_{L}^{2}+(w_{L}/c)^{2}\sin^{2}\theta\right]\Bigr{)}$
$\displaystyle\times$ $\displaystyle\exp\Bigl{(}-2i\pi\nu
d\left(1-\bm{\hat{r}}_{\perp}\cdot\bm{\beta}_{\perp}\right)/c\beta_{z}\Bigr{)}\,.$
(33)

We further assume that the momentum distribution function of the beam is a
separable function of absolute momentum $u$ and direction angles
$(\psi,\varphi)$, i.e.,
$g(\bm{u})\mathrm{d}^{3}\bm{u}=g_{u}(u)g_{\psi}(\psi)g_{\varphi}(\varphi)u^{2}\sin\psi\,\mathrm{d}u\,\mathrm{d}\psi\,\mathrm{d}\varphi$.
Inspired by Ref. 21, we use

$\displaystyle g_{u}(u)$ $\displaystyle=\frac{1}{2\Delta u^{3}}\exp(-u/\Delta
u)\,,$ (34) $\displaystyle g_{\psi}(\psi)$
$\displaystyle=\frac{1}{\Delta\psi^{2}}\frac{e^{-\frac{1}{2}\left(\sin\psi/\Delta\psi\right)^{2}}}{1-e^{-1/2\Delta\psi^{2}}}\cos\psi\,,$
(35) $\displaystyle g_{\varphi}(\varphi)$ $\displaystyle=1/2\pi\,,$ (36)

where $\Delta u$ and $\Delta\psi\equiv\sin\Psi_{h}$ represent characteristic
spreads in momentum and polar angle, respectively. The above distributions are
normalized such that

$\int_{0}^{\infty}g_{u}(u)\,u^{2}\,\mathrm{d}u=\int_{0}^{\pi/2}g_{\psi}(\psi)\,\sin\psi\,\mathrm{d}\psi=1\,.$
(37)

The momentum spread $\Delta u$ is chosen so that the average Lorentz factor

$\langle\gamma_{h}\rangle=\int_{0}^{\infty}g_{u}(u)u^{2}\sqrt{u^{2}+1}\,\mathrm{d}u\sim\sqrt{9\Delta
u^{2}+1}$ (38)

equals the ponderomotive potential of the (linearly polarized) laser pulse,
$\gamma_{L}=\sqrt{1+a_{L}^{2}/2}$, where $a_{L}$ is the dimensionless laser
field strength [53].

The number of hot electrons carried by the beam is estimated from

$N_{h}\simeq\frac{\eta_{h}\mathcal{E}_{L}}{m_{e}c^{2}\langle\gamma_{h}-1\rangle}=\frac{\mathcal{E}_{\rm
beam}}{\langle\mathcal{E}_{h}\rangle}\,,$ (39)

where $\mathcal{E}_{L}=I_{L}(w_{L}^{2}\tau_{L}/8)(\pi/\ln 2)^{3/2}$ is the
laser pulse energy,
$\langle\mathcal{E}_{h}\rangle=m_{e}c^{2}\langle\gamma_{h}-1\rangle$ the mean
electron kinetic energy, $\eta_{h}$ the laser-to-hot-electron energy
conversion efficiency, and $\mathcal{E}_{\rm beam}=\eta_{h}\mathcal{E}_{L}$
the total kinetic beam energy. Based on experimental measurements [54, 55], a
constant value $\eta_{h}=0.2$ will be assumed in the following.

Introducing $n_{\rm hr0}$ as the hot-electron density at the target backside,
the strength of the sheath field is expected to be [32, 36]

$E_{0}\simeq\alpha_{E_{0}}\sqrt{\frac{n_{\rm
hr0}\langle\mathcal{E}_{h}\rangle}{\varepsilon_{0}}}\,,$ (40)

where $\alpha_{E_{0}}\leq 1$ is a scaling factor expressing that the effective
average field strength experienced by the hot electrons in vacuum is lower
than its maximum value at $z=0$, $E_{\rm max}=\sqrt{n_{\rm
hr0}\langle\mathcal{E}_{h}\rangle/\varepsilon_{0}}$. We will take
$\alpha_{E_{0}}=0.5$ as a fiducial parameter. Moreover, $n_{\rm hr0}$ can be
related to the hot-electron density at the front side, $n_{\rm hf0}$, defined
by

$n_{\rm
hf0}\simeq\frac{\eta_{h}I_{L}}{m_{e}c^{3}\langle\gamma_{h}-1\rangle}\,.$ (41)

through

$n_{\rm hr0}=n_{\rm hf0}\left(\frac{w_{L}}{w_{h}(d)}\right)^{2}\,.$ (42)

The transverse size of the beam, $w_{h}(z)$, after a ballistic propagation of
a depth $z$ through the target is estimated [56] as

$w_{h}(z)\simeq
w_{L}\sqrt{1+\left(\frac{2z\tan\Psi_{h}}{w_{L}}\right)^{2}}\,.$ (43)

The sheath field $E_{0}$ is obtained from Eqs. (40)-(43).

##  III Phenomenology of CTSR and parametric dependencies

###  III.1 Interplay of transition and synchrotron radiations

Figure 3: (a) Frequency-angle spectrum $\partial^{2}\mathcal{I}_{\rm
CTSR}/\partial\nu\partial\Omega$ (in $\rm\mu J\,ps\,sr^{-1}$ units) radiated
by hot electrons in a conducting foil of $d=2\,\rm\mu m$ thickness, exposed to
a $a_{L}=15$, $w_{L}=5\,\rm\mu m$, $\tau_{L}=30\,\rm fs$,
$\mathcal{E}_{L}=2.8\,\rm J$ laser pulse. The laser-to-hot-electron conversion
efficiency is set to $\eta_{h}=0.2$, leading to an electron beam energy of
$\mathcal{E}_{\rm beam}=0.56\,\rm J$. The white dashed line plots Eq. (49)
with $\eta_{\rm coh}=3\times 10^{-3}$. Panels (b) and (c) show the spectra (in
$\rm mJ\,ps\,sr^{-1}$ units) due to solely CTR and CSR, respectively.

As a first illustration of our modeling, we consider the case of a $2\,\rm\mu
m$ thick target exposed to an intense laser pulse characterized by $a_{L}=15$,
$\lambda_{L}=1\,\rm\mu m$, $\tau_{L}=30\,\rm fs$ and $w_{L}=5\,\rm\mu m$.
These parameters yield a peak pulse intensity $I_{L}=3.1\times 10^{20}\,\rm
W\,cm^{-2}$ and pulse energy $\mathcal{E}_{L}=2.8\,\rm J$. The beam divergence
is taken to be $\Psi_{h}=30^{\circ}$. For these parameters (together with
$\eta_{h}=0.2$), one obtains $\mathcal{E}_{\rm beam}\simeq 0.56\,\rm J$,
$\mathcal{E}_{h}\simeq 4.9\,\rm MeV$, $N_{h}\simeq 7.03\times 10^{11}$ and
$E_{0}\simeq 6.93\times 10^{12}\,\rm V\,m^{-1}$.

Figure 3(a) shows the angular distribution of the full energy spectrum (in
$\rm\mu J\,ps\,sr^{-1}$ units) for $0.1\leq\nu\leq 40\,\rm THz$, while Figs.
3(b) and (c) display the spectra (in $\rm mJ\,ps\,sr^{-1}$ units) due to only
CTR and CSR, respectively. One can see that the full spectrum is peaked around
$\nu\simeq 10-25\,\rm THz$ and $\theta\simeq 60-90^{\circ}$. By contrast, CTR
and CSR give rise separately to almost identical – but much more intense – THz
spectra concentrated at lower frequencies and angles.

The latter result can be readily understood upon noting that the first
(second) half of the electron trajectory in vacuum, before (resp. after) its
turning point, can be viewed, given its femtosecond timescale ($t\sim
m_{e}\gamma_{\mathrm{e}}c/eE_{0}\sim 1\,\rm fs$), as a sudden deceleration
along $\bm{\hat{z}}$ (resp. acceleration towards $-\bm{\hat{z}}$) with respect
to THz-range synchrotron radiation. Since, by contrast, transition radiation
at the exit (return) time corresponds to a sudden apparent acceleration (resp.
deceleration), the two radiation mechanisms generate very similar THz fields
but of opposite polarity and so almost exactly cancel out, especially at low
frequencies ($\lesssim 10\,\rm THz$ here) and angles ($\lesssim 40^{\circ}$).
Mathematically speaking, this description merely amounts to taking $\Delta
t_{\mathrm{r}}\to 0$ in Eq. (27) so that Eq. (29) converges to Eq. (28). It is
to be noted that the average time spent by the electrons in vacuum,
$\langle\Delta t_{\mathrm{r}}\rangle$, can be exactly solved as

$\displaystyle\langle\overline{\Delta t}_{\mathrm{r}}\rangle$
$\displaystyle=\iint g(u,\psi)\overline{\Delta
t}_{\mathrm{r}}u^{2}\sin\psi\,\mathrm{d}^{2}u\mathrm{d}\psi$ (44)
$\displaystyle=-3\Delta
u\left(1+\coth\frac{X_{\psi}^{2}}{2}\right)\left(X_{\psi}^{-1}D\left(X_{\psi}\right)-1\right)\,,$
(45)

where $D(x)=e^{-x^{2}}\int_{0}^{x}e^{t^{2}}\,\mathrm{d}t$ is the Dawson
function, $X_{\psi}\equiv 1/\sqrt{2}\Delta\psi$ and $\overline{\Delta
t}_{\mathrm{r}}=2u_{z}=2u\cos\psi$. For the parameters of Fig. 3, one finds
$\langle\Delta t_{\mathrm{r}}\rangle\simeq 4\,\rm fs$.

Such a strong interplay of CTR and CSR therefore precludes their separate
treatment, otherwise the low-frequency part of the spectrum would be
inaccurately modeled and, worse, the radiated energy would be unphysically
overestimated. Indeed, while the integration of $\partial^{2}\mathcal{I}_{\rm
CTSR}/\partial\nu\partial\Omega$ over $0.1<\nu<40\,\rm THz$ and $\Omega$ gives
a total radiated energy of $\mathcal{E}_{\rm CTCR}\simeq 4.86\times
10^{-3}\,\mathcal{E}_{\rm beam}\simeq 2.7\,\rm mJ$, the CTR and CSR spectra,
both scaling as $N_{h}^{2}$, contain alone an energy much exceeding the total
beam energy, i.e., $\mathcal{E}_{\rm CTR}\simeq\mathcal{E}_{\rm CSR}\simeq
88\,\mathcal{E}_{\rm beam}$.

Besides the high-pass filtering effect caused by the ultrashort reflection
timescale of the electrons in vacuum, the form factor $F(\theta,\nu,\bm{u})$
that characterizes the spatiotemporal coherence of the entire beam acts as a
low-pass filter. The high-frequency shape of the spectrum can be approximately
captured by the attenuation factor $\eta_{\rm coh}$ defined by

$\displaystyle\eta_{\rm coh}(\theta,\nu)$ $\displaystyle=\left|\int
g(\bm{u})F(\theta,\nu,\bm{u})\mathrm{d}^{3}\bm{u}\right|^{2}$
$\displaystyle=C_{\rm coh}\,e^{-(\pi^{2}/2\ln
2)\bar{\nu}^{2}\left(\bar{\tau}_{L}^{2}+\bar{w}_{L}^{2}\sin^{2}\theta\right)}\,,$
(46)

where

$C_{\rm coh}=\left|\int
g(\bm{u})e^{-2i\pi\bar{\nu}\bar{d}\left[1-\bm{\hat{r}}\cdot\bm{\beta}_{\perp}\right]/\beta_{z}}\mathrm{d}^{3}\bm{u}\right|^{2}$
(47)

quantifies the coherence of the beam after its propagation through the target
of thickness $d$. If the target is sufficiently thin, the propagation effects
can be neglected and since $\int g(\bm{u})\mathrm{d}^{3}\bm{u}=1$, one has
$C_{\rm coh}\simeq 1$. In this case, the contour lines of the low-pass filter
can be computed as a function of $\theta$ and $\nu$ through

$\eta_{\rm coh}(\theta,\nu)\simeq e^{-(\pi^{2}/2\ln
2)\bar{\nu}^{2}\left(\bar{\tau}_{L}^{2}+\bar{w}_{L}^{2}\sin^{2}\theta\right)}\,.$
(48)

As seen in Fig. 3(a), the radiated spectrum assumes the shape of a band
located, for our parameters, between $\nu\simeq 10\,\rm THz$ and the
approximate upper bound

$\bar{\nu}_{c}(\theta)\simeq\sqrt{\frac{-2\ln 2\ln\eta_{\rm
coh}}{\pi^{2}\left(\bar{\tau}_{L}^{2}+\bar{w}_{L}^{2}\sin^{2}\theta\right)}}\,.$
(49)

This expression is plotted as a white dashed line in Fig. 3(a) for $\eta_{\rm
coh}(\theta,\bar{\nu}_{c})=3\times 10^{-3}$, leading to
$25\leq\nu_{c}(\theta)\leq 30\,\rm THz$.

###  III.2 Influence of the sheath field

Equation (40) represents a crude expression of the sheath field strength
which, in reality, varies both with space and time [32, 36]. To examine the
dependency of the THz radiation spectrum on $E_{0}$, we perform a parametric
study over $\alpha_{E_{0}}$.

Increasing the sheath field strength shortens the electron excursion time in
vacuum: the CTR and CSR fields then compensate each other even better so that
the energy radiated by CTSR is further reduced. Conversely, decreasing $E_{0}$
allows the particles to propagate over larger distances and radiate more
efficiently.

Figure 4: Radiated THz energy (integrated over solid angles and $0.1-40\,\rm
THz$ frequencies) as a function of the sheath-field parameter
$\alpha_{E_{0}}$, as defined by Eq. (40). The blue solid curve, in $\rm
log_{10}$ scale, represents the full radiated energy while the red dashed and
dotted curves plot the energy yields of CTR and CSR, respectively. The laser-
plasma parameters are those used in Fig. 3 ($a_{L}=15$). The reference case
$\alpha_{E_{0}}=0.5$ corresponds to a field $E_{0}\simeq 6.9\times
10^{12}\,\rm V\,m^{-1}$.

Figure 4 shows the variation in the radiated energy (integrated over solid
angles and the $0.1-40$ THz frequency range) with $\alpha_{E_{0}}\leq 1$, for
the same laser-plasma parameters as in Fig. 3. The energies radiated through
CTR and CSR, taken separately, decrease when lowering $\alpha_{E_{0}}$ whilst
the total radiated energy rapidly rises (by $\sim 10\times$ from
$\alpha_{E_{0}}=1$ to $\alpha_{E_{0}}=0.5$, and by $\sim 20\times$ from
$\alpha_{E_{0}}=0.5$ to $\alpha_{E_{0}}=0.2$). As expected, CSR tends to
vanish when $\alpha_{E_{0}}\to 0$ so that CTR then accounts for all of the
radiated energy. Yet, the electrons then only radiate coherently while exiting
the target because they are greatly spread out upon re-entering the foil. When
$\alpha_{E_{0}}\to 1$, by contrast, they are reflected so rapidly into the
target that they remain packed enough to emit two almost coincident CTR bursts
[stages (iii) and (v) in Fig. 2] of same polarity (the first from the
accelerating real particles, the second from the decelerating image
particles), the total energy of which is four times that of the initial CTR
burst.

Of course, $\alpha_{E_{0}}$ is just a fitting parameter, introduced for
simplicity to capture the overall effect of the highly nonstationary sheath
field [32, 44]. Precautions must therefore be taken when analyzing the
radiated energy for very weak sheath fields. Notably, below
$\alpha_{E_{0}}\simeq 0.2$, the hot electrons spend so much time ($\gtrsim
10\,\rm fs$) in vacuum before returning to the target that a sizable ($>0.1$)
fraction of the beam energy is predicted to be radiated away through CTSR.
Even worse, below $\alpha_{E_{0}}\simeq 0.1$, the total emitted energy
approaches that potentially radiated by all accelerated electrons through CTR
alone, entailing non-physically high ($>\mathcal{E}_{\rm beam}$) radiated
energies. A proper treatment of the problem for such weak sheath fields would
therefore require a self-consistent description of the backreaction of
radiation on the electron dynamics, as is done, _e.g._ , in accelerator
physics [57]. This is a challenging task, well exceeding the scope of the
present semianalytical model.

###  III.3 Contribution of escaping ballistic electrons

Figure 5: Radiated THz energy (integrated over solid angles and $0.1-40\,\rm
THz$ frequencies) as a function of the fraction of escaping electrons
$\eta_{\rm esc}$. The blue solid curve, in $\rm log_{10}$ scale, represents
the full radiated energy while the red dashed and dotted curves plot the
energy yields of CTR and CSR, respectively. The laser-plasma parameters are
those used in Fig. 3 ($a_{L}=15$). The radiated spectra of cases (a)-(c) are
plotted in Fig. 6.

As seen previously, very efficient compensation of the radiated THz fields
takes place when the fast electrons are rapidly pulled back into the target,
as happens under the considered interaction conditions. We have so far assumed
that all of the fast electrons are drawn back into the target, yet it is well
known that a higher-energy fraction of them can escape the sheath potential
[33, 34]. If these escaping electrons are numerous enough, their uncompensated
transition radiation may dominate the total radiation yield. Here, we will
assess the impact on the overall THz radiation of an increasing fraction
($\eta_{\rm esc}$, considered as a free parameter) of escaping electrons.

Figure 6: Frequency-angle spectra for various fractions of escaping electrons:
(a) $\eta_{\rm esc}=1\,\%$, (b) $\eta_{\rm esc}=5\,\%$ and (c) $\eta_{\rm
esc}=10\,\%$, corresponding to cases (a)-(c) in Fig. 5. The other laser-plasma
parameters are those used in Fig. 3 ($a_{L}=15$). The cutoff frequency
$\nu_{c}(\theta)$ is obtained from Eqs. (49) and (47) with $\eta_{\rm
coh}=3\times 10^{-3}$.

To this purpose, we introduce the cutoff energy $u_{c}(\eta_{\rm esc})\equiv
p_{c}(\eta_{\rm esc})/m_{e}c$ beyond which the electrons escape the target,
which fulfills

$\displaystyle\eta_{\rm esc}$ $\displaystyle=\int_{u_{\rm
c}}^{\infty}g_{u}(u)u^{2}\mathrm{d}u$ $\displaystyle=e^{-u_{\rm c}/\Delta
u}(2\Delta u^{2}+u_{\rm c}^{2}+2u_{\rm c}\Delta u)/2\Delta u^{2}\,,$ (50)

where $\eta_{\rm esc}$ is usually estimated to be in the percent range [34].
To simplify, we assume that the escaping electrons (characterized by $u\geq
u_{\rm c}$) keep on propagating ballistically after exiting the target
backside, and hence only emit a single burst of CTR [stage (iii) in Fig. 2].
Equations (25), (28) and (29) can thus be recast as ($\bar{\nu}\equiv\nu
m_{e}c/eE_{0}$):

$\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
CTSR}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)$
$\displaystyle=\frac{N_{h}^{2}q^{2}}{2\pi^{2}\varepsilon_{0}c}\left|\int_{0}^{u_{c}}\sideset{}{{}_{\psi,\varphi}}{\iint}g(\bm{u})F(\theta,\nu,\bm{u})i\pi\bar{\nu}\Bigl{[}\widebar{\bm{A}}^{-}(\bm{\hat{r}},\nu)-\widebar{\bm{A}}^{+}(\bm{\hat{r}},\nu)\Bigr{]}\mathrm{d}^{3}\bm{u}\right.$
$\displaystyle\hskip
28.0pt+\frac{1}{2}\left.\int_{u_{c}}^{\infty}\sideset{}{{}_{\psi,\varphi}}{\iint}g(\bm{u})F(\theta,\nu,\bm{u})\left[\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{e}}^{-}}{1-\bm{\beta}_{\mathrm{e}}^{-}\cdot\bm{\hat{r}}}-\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{e}}^{+}}{1-\bm{\beta}_{\mathrm{e}}^{+}\cdot\bm{\hat{r}}}\right]e^{-i\Theta_{\mathrm{e}}(\bm{u})}\,\mathrm{d}^{3}\bm{u}\right|^{2}\,,$
(51) $\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
CTR}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)$
$\displaystyle=\frac{N_{h}^{2}q^{2}}{2\pi^{2}\varepsilon_{0}c}\left|\int_{0}^{u_{c}}\sideset{}{{}_{\psi,\varphi}}{\iint}g(\bm{u})F(\theta,\nu,\bm{u})\bar{\bm{A}}_{\rm
TR}\left(\bm{\hat{r}},\nu\right)\,\mathrm{d}^{3}\bm{u}\right.$
$\displaystyle\hskip
28.0pt+\frac{1}{2}\left.\int_{u_{c}}^{\infty}\sideset{}{{}_{\psi,\varphi}}{\iint}g(\bm{u})F(\theta,\nu,\bm{u})\left[\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{e}}^{-}}{1-\bm{\beta}_{\mathrm{e}}^{-}\cdot\bm{\hat{r}}}-\frac{\bm{\hat{r}}\times\bm{\beta}_{\mathrm{e}}^{+}}{1-\bm{\beta}_{\mathrm{e}}^{+}\cdot\bm{\hat{r}}}\right]e^{-i\Theta_{\mathrm{e}}(\bm{u})}\,\mathrm{d}^{3}\bm{u}\right|^{2}\,,$
(52) $\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
CSR}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)$
$\displaystyle=\frac{N_{h}^{2}q^{2}}{2\pi^{2}\varepsilon_{0}c}\left|\int_{0}^{u_{c}}\sideset{}{{}_{\psi,\varphi}}{\iint}g(\bm{u})F(\theta,\nu,\bm{u})\Biggl{\\{}\vphantom{\int}i\pi\bar{\nu}\Bigl{[}\widebar{\bm{A}}^{-}(\bm{\hat{r}},\nu)-\widebar{\bm{A}}^{+}(\bm{\hat{r}},\nu)\Bigr{]}-\bar{\bm{A}}_{\rm
TR}\left(\bm{\hat{r}},\nu\right)\Biggr{\\}}\,\mathrm{d}^{3}\bm{u}\right|^{2}\,.$
(53)

Figure 5 shows the evolution of the radiation yield in the $0.1-40\,\rm THz$
frequency range when $\eta_{\rm esc}$ increases from $0\,\%$ to $10\,\%$. We
observe that the CTR yield progressively diminishes with $\eta_{\rm esc}$
while the CSR yield decreases twice as fast.

The decreasing CTR yield is due to the fact that, as $\eta_{\rm esc}$ rises,
fewer and fewer electrons emit CTR at time $t=t_{\mathrm{r}}$ while the number
of electrons emitting CTR at $t=t_{\mathrm{e}}$ remains unchanged. Similarly,
the CSR yield drops because fewer particles are decelerated. Therefore, when
increasing $\eta_{\rm esc}$, one source of CTR is preserved [stage (iii) in
Fig. 2] whereas the whole source of CSR is degraded [stage (iv) of Fig. 2],
which explains the trends seen in Fig. 5.

In parallel, the total (CTSR) THz yield weakly varies (reaching a slight
minimum around $\eta_{\rm esc}\lesssim 1\,\%$) for $\eta_{\rm esc}\lesssim
2\,\%$ but rapidly rises beyond this value, by $\sim 65\times$ between
$\eta_{\rm esc}=1\,\%$ and $\eta_{\rm esc}=10\,\%$. Note that our neglect of
radiative losses becomes questionable at $\eta_{\rm esc}=10\,\%$, for which
the yield attains $\mathcal{E}_{\rm CTSR}\simeq 0.2\,\mathcal{E}_{\rm beam}$.

Figure 6 displays the frequency-angle spectra associated with $\eta_{\rm
esc}=1\,\%$, $5\,\%$ and $10\,\%$, corresponding, respectively, to markers
(a), (b) and (c) in Fig. 5. When $\eta_{\rm esc}$ rises from zero, the
spectrum develops an increasingly prominent low-frequency structure, with a
large-angle tail around $\theta\simeq 40-70^{\circ}$, quite above the electron
beam’s divergence ($\Psi_{h}=30^{\circ}$). For thin ($d=2\,\rm\mu m$) targets,
the hierarchy between CTSR from the confined electrons and CTR from the
escaping electrons is reversed at around $\eta_{\rm esc}\simeq 1\,\%$ [Fig.
6(a)], leading to a doubly peaked spectrum, yet containing a comparable
integrated energy ($\simeq 3.5\times 10^{-3}\,\mathcal{E}_{\rm beam}\simeq
2.0\,\rm mJ$) as for $\eta_{\rm esc}=0\%$ (Fig. 3). When $\eta_{\rm esc}>1\%$
[Fig. 6(b,c)], the low-frequency signal from the escaping electrons largely
prevails.

Figure 7: Frequency-angle spectra for a $20\,\rm\mu m$ thick foil and various
fractions of escaping electrons: (a) $\eta_{\rm esc}=1\,\%$, (b) $\eta_{\rm
esc}=5\,\%$ and (c) $\eta_{\rm esc}=10\,\%$. The other laser-plasma parameters
are those used in Fig. 3 ($\mathcal{E}_{L}=2.8\,\rm J$, $\mathcal{E}_{\rm
beam}=0.56\,\rm J$). The cutoff frequency $\nu_{c}(\theta,d)$, plotted as a
white dashed line, is obtained from Eqs. (III.1) and (47) with $\eta_{\rm
coh}=3\times 10^{-3}$.

###  III.4 Influence of the target thickness

As the electron beam propagates ballistically through the target, it spreads
transversely and longitudinally depending on its momentum distribution
function. The increased electron dilution at the backside of a thicker target
also entails a weaker sheath field, thus allowing the electrons to propagate
in vacuum over larger distances and longer times. As seen in Sec. III.2, this
tends to reduce the respective yields and frequency ranges of CTR and CSR but
also, and above all, to hamper the field cancellation of these two emissions,
hence causing a net increase in the CTSR yield.

To illustrate this behavior, we plot in Fig. 7 the total THz spectra obtained
for $d=20\,\rm\mu m$ and an escaping electron fraction ranging from $1\,\%$ to
$10\,\%$. The main observations are that the radiated energy is $\sim
16\times$ higher than at $d=2\,\rm\mu m$ and that the spectrum remains
essentially unchanged up to $\eta_{\rm esc}\simeq 2-3\,\%$. This means that
the THz radiation is less sensitive to the (unknown) fraction of escaping
electrons in thicker foils. When $\eta_{\rm esc}=10\,\%$, the signal from the
recirculating electrons is admittedly less intense than the (lower-frequency)
one due to the escaping electrons, but still visible. Moreover, compared to
the case of $d=2\,\rm\mu m$ and $\eta_{\rm esc}=0\,\%$, the spectra associated
with $d=20\,\rm\mu m$ have shifted to lower angles ($\sim 25-45^{\circ}$) and
frequencies ($\sim 10\,\rm THz$). This behavior results from the lower
coherence of the electron beam (characterized by $C_{\rm coh}$, see Sec.
III.1) at the rear side of the target. To illustrate the effect of coherence
loss, we overlay in Fig. 7 the cutoff frequency $\nu_{c}(\theta,d)$ (white
dashed line) as obtained by numerically solving Eqs. (III.1) and (47) for
$\eta_{\rm coh}=3\times 10^{-3}$. This cutoff frequency captures fairly well
the high-$\nu$ shape of the CTSR spectrum.

The prediction of a THz yield increasing with the target thickness between
$d=2\,\rm\mu m$ and $d=20\,\rm\mu m$ and for $\eta_{\rm esc}=1-5\,\%$ [Figs.
6(a,b) and 7(a,b)] may seem dubious given the absence of supporting
experimental evidence, but should be considered along with the $\eta_{\rm
esc}=10\,\%$ case. The latter indeed suggests that, for a given fraction of
escaping electrons $\eta_{\rm esc}$, the THz yield will eventually drop with
the target thickness as CTSR, which is increasingly dominated by CTR when
$\eta_{\rm esc}>1\,\%$, weakens due to degraded beam coherence.

###  III.5 Influence of the laser focusing

We now examine how the THz radiation depends on the laser pulse intensity
($\propto\,a_{L}^{2}$) at fixed laser energy ($\propto\,a_{L}^{2}w_{L}^{2}$).
To do so, we vary the laser field strength $a_{L}$ while adjusting the laser
spot size $w_{L}\,\propto\,1/a_{L}$ and keeping the laser-to-hot-electron
conversion efficiency fixed to $\eta_{h}=0.2$. Since
$\langle\gamma_{h}\rangle\,\propto\,a_{L}$, lowering $a_{L}$ causes the number
of hot electrons $N_{h}\,\propto\,\mathcal{E}_{L}/\langle\gamma_{h}\rangle$ to
rise, leading to stronger separate CTR and CSR yields. Predicting the net
effect on CTSR, though, is not straightforward as it is extremely sensitive to
the normalized return time $\langle\widebar{\Delta
t_{r}}\rangle\,\propto\,\Delta u\,\propto\,a_{L}$, as shown by Eqs. (25) -
(27) and Secs. III.1 and III.2.

Figure 8: Frequency-angle spectra (in $\rm\mu J\,ps\,sr^{-1}$ units) for (top)
$d=2\,\rm\mu m$ and (bottom) $d=20\,\rm\mu m$ at $a_{L}=5$ and fixed laser
pulse energy $\mathcal{E}_{L}=2.8\,\rm J$ ($\mathcal{E}_{\rm beam}=0.56\,\rm
J$). The fraction of escaping electrons varies from left to right: (a,d)
$\eta_{\rm esc}=0\,\%$, (b,e) $\eta_{\rm esc}=1\,\%$ and (c,f) $\eta_{\rm
esc}=5\,\%$. Apart from $d$, $a_{L}$ and $w_{L}\propto a_{L}^{-1}$, the other
parameters are as in Fig. 3. Panels (b,c) (resp. (e,f)) are to be compared
directly to their $a_{L}=15$ equivalent [Figs. 6(a,b) and 7(a,b)]. The cutoff
frequencies (a) $\nu_{c}(\theta)$ and (d,e) $\nu_{c}(\theta,d)$ are computed
from, respectively, Eq. (49) and Eqs. (III.1)-(47) with $\eta_{\rm
coh}=1\times 10^{-3}$.

Figure 8 plots the THz radiated spectra for (top) $d=2\,\rm\mu m$ and (bottom)
$d=20\,\rm\mu m$, the laser pulse energy being set to
$\mathcal{E}_{L}=2.8\,\rm J$, the laser amplitude to $a_{L}=5$, and thus the
waist to $w_{L}=15\,\rm\mu m$. Panels (b,c) [(e,f)] are to be compared
directly to their $a_{L}=15,w_{L}=5\,\rm\mu m$ [Fig. 6(a,b) [resp. Fig.
7(a,b)] equivalent. For $d=2\,\rm\mu m$, and assuming full refluxing
($\eta_{\rm esc}=0\,\%$), the total radiated energy ($\mathcal{E}_{\rm
CTSR}\simeq 5.2\times 10^{-3}\,\mathcal{E}_{\rm beam}$) is almost identical to
that achieved at $a_{L}=15$, yet the radiation is more sensitive to the
escaping electrons: their contribution becomes overwhelmingly dominant even
when their fraction is as low as $1\,\%$, in which case the radiation yield is
more than doubled ($\mathcal{E}_{\rm CTSR}\simeq 1.3\times
10^{-2}\,\rm\mathcal{E}_{\rm beam}$). Again, when $\eta_{\rm esc}$ is
sufficiently low, increasing the target thickness from $d=2\,\rm\mu m$ to
$d=20\,\rm\mu m$ significantly enhances the radiated energy, albeit to a
lesser degree than at $a_{L}=15$. The contribution of the recirculating
electrons is then still significant at $\eta_{\rm esc}=1\,\%$, but no longer
so at $\eta_{\rm esc}=5\,\%$. Note, however, that for both $d=2\,\rm\mu m$ and
$d=20\,\rm\mu m$, the model’s predictions become questionable when $\eta_{\rm
esc}\gtrsim 5\,\%$ since $\mathcal{E}_{\rm CTSR}$ then exceeds $\sim
0.3\,\mathcal{E}_{\rm beam}$.

Interestingly, the $\sim 10^{-2}$ beam-to-THz conversion efficiency predicted
for $a_{L}=5$ and $\eta_{\rm esc}=1\,\%$ appears to be roughly consistent with
the $\sim 1\,\rm mJ$ yield reported in Ref. 26 under comparable interaction
conditions ($a_{L}\simeq 5$, $\tau_{L}\simeq 30\,\rm fs$, $w_{L}\simeq
10\,\rm\mu m$), but with a thicker foil ($d=5\,\rm\mu m$) and a lower cutoff
frequency ($\nu<9\,\rm THz$).

Overall, our model pinpoints the critical influence of some of its parameters,
notably the sheath field strength and the fraction of escaping electrons. Yet,
to our knowledge, most experiments to date have characterized THz emissions
from laser-foil interactions over restricted frequency and angular ranges[25,
26], with the notable exception of Ref. 58 where the full THz distribution was
diagnosed, but with an obliquely incident laser pulse and no variation in the
target thickness was reported.

##  IV Radiation from the expanding plasma

We now address the THz radiation that is subsequently emitted by the electron-
proton plasma expanding due to the fast-electron-induced sheath field.
Compared to previous efforts [59, 26, 27, 28, 29], our modeling will hinge on
a more realistic description of the plasma acceleration, taking account of the
time-decreasing areal charge at the accelerated ion front and of the adiabatic
electron cooling taking place in thin foils.

###  IV.1 General formalism

A convenient formula to describe the far-field plasma expansion radiation
(PER) is [60, 61]

$\frac{\partial^{2}\mathcal{I}_{\rm
PER}}{\partial\nu\partial\Omega}(\bm{\hat{r}},\nu)=\frac{1}{8\pi^{2}\varepsilon_{0}c^{3}}\\\
\times\left|\iint\bm{\hat{r}}\times\left[\frac{\partial\bm{j}(\bm{r}^{\prime},t)}{\partial
t}\right]_{t_{\rm ret}}e^{-2i\pi\nu
t}\mathrm{d}t\mathrm{d}^{3}r^{\prime}\right|^{2}\,,$ (54)

where $\bm{j}$ is the plasma current density which will be estimated below. We
will assume, for simplicity, that the plasma is accelerated along the target
normal only. Thus, $\bm{j}(\bm{r},t)=j(\bm{r},t)\bm{\hat{z}}$ and the radiated
spectrum associated with Eq. (54) reads

$\frac{\partial^{2}\mathcal{I}_{\rm
PER}}{\partial\nu\partial\Omega}(\theta,\nu)=\frac{\sin^{2}\theta}{8\pi^{2}\varepsilon_{0}c^{3}}\\\
\times\left|\iint\left[\frac{\partial j}{\partial
t}(\bm{r}^{\prime},t)\right]_{t_{\rm ret}}e^{-2i\pi\nu
t}\,\mathrm{d}t\,\mathrm{d}^{3}\bm{r}^{\prime}\right|^{2}\,.$ (55)

Performing the change of variable $t\to t_{\rm ret}$ leads to

$\frac{\partial^{2}\mathcal{I}_{\rm
PER}}{\partial\nu\partial\Omega}(\theta,\nu)=\frac{\sin^{2}\theta}{8\pi^{2}\varepsilon_{0}c^{3}}\\\
\times\left|\iint\frac{\partial j}{\partial
t}(\bm{r}^{\prime},t^{\prime})e^{-2i\pi\nu\left(t^{\prime}-\frac{\bm{\hat{r}}\cdot\bm{r}^{\prime}}{c}\right)}\,\mathrm{d}t^{\prime}\,\mathrm{d}^{3}\bm{r}^{\prime}\right|^{2}\,.$
(56)

This formula is the basis of the PER model constructed below.

###  IV.2 Plasma expansion in the isothermal and adiabatic regimes

To evaluate the current density $j(\bm{r},t)$, we need to describe the
dynamics of the expanding plasma. To do so, we start by considering the one-
dimensional (1D) model proposed in Ref. 36, which applies to a collisionless
plasma accelerated towards vacuum by the sheath field created by an isothermal
electron population. According to this model, the plasma ion profile,
initialized as a step function, retains a sharp front located at $z_{\rm
f,iso}(t)$ and moving at the velocity

$\dot{z}_{\rm f,iso}(t)\equiv\beta_{\rm f,iso}(t)c\simeq 2c_{\rm
s0}\ln\left(\tau+\sqrt{\tau^{2}+1}\right)\,,$ (57)

where $c_{\rm s0}\simeq\sqrt{\langle\mathcal{E}_{h}\rangle/m_{p}}$ is the ion
acoustic velocity, $m_{p}$ the proton mass, $\tau=\omega_{pi}t/\sqrt{2e_{N}}$
the normalized time, $\omega_{pi}=\sqrt{n_{\rm
hr0}e^{2}/m_{p}\varepsilon_{0}}$ the ion plasma frequency,
$e_{N}\equiv\exp(1)$, and $n_{\rm hr0}$ the initial electron density at the
target rear side. We estimate $n_{\rm hr0}$ as in Eq. (42) but consider in
addition a folding term $(1+c\tau_{L}/2d)$ to describe the electron
accumulation in thin foils during the laser irradiation [62]:

$n_{\rm
hr0}=\frac{\eta_{h}I_{L}}{m_{e}c^{3}\langle(\gamma_{h}-1)\beta_{z}\rangle}\left(\frac{w_{h}(0)}{w_{h}(d)}\right)^{2}\left(1+\frac{c\tau_{L}}{2d}\right)\,.$
(58)

When the longitudinal extent $\sim c\tau_{L}$ of the fast electron bunch is
short compared to the target thickness, i.e., when $c\tau_{L}\ll 2d$, the
folding term approaches unity and the electron density is that given by Eq.
(42). Conversely, when $c\tau_{L}\gtrsim 2d$, the fast electrons recirculate
several times before the pulse ends and the folding term increases the
electron density accordingly.

The approximation of isothermal hot electrons ceases to be valid when the
rarefaction waves coming from the two sides of the foil have converged to its
center, which occurs at a time $t_{\rm ad}\simeq d/2c_{\rm s0}$ [63]. The
adiabatic cooling experienced by the fast electrons at $t>t_{\rm ad}$
precipitates the decay of the sheath field, so that the ions eventually reach
a maximum velocity. In a 1D geometry, this maximum velocity is predicted to be
[63]

$\dot{z}_{\rm f,ad}\simeq 2c_{\rm
s0}\ln\left(0.32d/\lambda_{D0}+4.2\right)\,,$ (59)

where

$\lambda_{D0}=c_{\rm
s0}/\omega_{pi}=\sqrt{\varepsilon_{0}\langle\mathcal{E}_{h}\rangle/n_{\rm
hr0}e^{2}}$ (60)

is the initial Debye length of the hot electrons.

To model the dynamic transition between the isothermal and adiabatic (cooling)
expansion regimes, we use the simple interpolation formula proposed in Ref.
64:

$\dot{z}_{\rm f,1D}(t)\simeq\left[\dot{z}^{-2}_{\rm f,iso}(t)+\dot{z}_{\rm
f,ad}^{-2}\right]^{-1/2}.$ (61)

As the time-dependent sheath field (as seen by the front ions) fulfills
$\ddot{z}_{\rm f,1D}(t)=eE_{\rm f,\,1D}(t)/m_{p}$, we deduce

$E_{\rm f,1D}(t)=\frac{m_{p}}{e}\frac{\ddot{z}_{\rm
f,iso}(t)}{{\left[1+\dot{z}^{2}_{\rm f,iso}(t)/\dot{z}^{2}_{\rm
f,ad}\right]^{3/2}}}\,.$ (62)

Next, in order to consider the effect of the transverse spreading of the hot
electrons, we note that the electric field initially scales as $\sqrt{n_{\rm
hr0}}$ in the purely 1D case [36], and so apply for $t>\tau_{L}$ a correction
factor $\sqrt{n_{\rm hr}(t)/n_{\rm hr0}}=w_{h}(d)/w_{h}(d+c(t-\tau_{L}))$ to
$E_{\rm f,1D}$ above. Finally, following Refs. 65, 64, we take into account
the expected weakening of the sheath field when the ion front has moved a
distance $z_{\rm f}$ greater than its transverse extent $w_{h}(d)$ [$w_{h}(z)$
being the transverse electron beam size given by Eq. (43)]. The sheath field
is then expected to decay in time as $z_{\rm f}(t)^{-2}$ in a realistic 3D
geometry [65, 64].

All these considerations motivate the following approximate expression of the
sheath field

$\displaystyle E_{\rm f,3D}(t)$ $\displaystyle=\begin{cases}E_{\rm
f,1D}(t)\left[1\\!+\\!\frac{z^{2}_{\rm
f,3D}(t)}{w_{h}^{2}(d)}\right]^{-1}&,t<\tau_{L}\\\\[6.0pt] \frac{E_{\rm
f,1D}(t)w_{h}(d)}{w_{h}(d+c(t-\tau_{L}))}\left[1\\!+\\!\frac{z^{2}_{\rm
f,3D}(t)}{w_{h}^{2}(d)}\right]^{-1}&,t\geq\tau_{L}\end{cases}$ (63)
$\displaystyle\ddot{z}_{\rm f,3D}(t)$ $\displaystyle=\frac{e}{m_{p}}E_{\rm
f,3D}(t)\,.$ (64)

The above equations are solved numerically to obtain the time-varying
position, velocity and acceleration of the (fastest) front ions.

As an example, Fig. 9 compares the normalized proton front velocities,
$\beta_{\rm f}=\dot{z}_{\rm f}/c$, as obtained from the above models. The
laser-plasma parameters ($a_{L}=15$, $w_{L}=5\,\rm\mu m$, $\tau_{L}=30\,\rm
fs$, $d=2\,\rm\mu m$) are the same as in Fig. 3. Plotted are the maximum
adiabatic velocity [63] $\dot{z}_{\rm f,ad}/c$, the unbounded ion velocity in
the isothermal [36] regime, $\dot{z}_{\rm f,iso}/c$, and the ion velocity in
the mixed isothermal/adiabatic regime without ($\dot{z}_{\rm f,1D}/c$) or with
($\dot{z}_{\rm f,3D}/c$) 3D corrections [64]. The rapid saturation of the ion
velocity (reaching a maximum value $\beta_{\rm f,3D}\simeq 0.18$ by $t\simeq
0.1\,\rm ps$) that is predicted by the 3D adiabatic model is particularly
manifest.

Figure 9: Time evolution of the ion front velocity as predicted by various
models. Green dotted dashed curve: isothermal model [$\beta_{\rm iso,1D}$ from
Eq. (57)]. Blue dashed curve: mixed 1D isothermal/adiabatic model [$\beta_{\rm
f,1D}$ from Eq. (61)]. Orange solid curve: mixed 3D isothermal/adiabatic
regime [$\beta_{\rm f,3D}$ from Eq. (64)]. Black dashed curve: maximum
velocity in the adiabatic regime [$\beta_{\rm f,ad}$ from Eq. (59)]. The
laser-plasma parameters are those used in Fig. 3.

###  IV.3 Evaluation of the plasma current density

The plasma current density $j(\bm{r},t)$ (resulting from both electron and ion
contributions) involved in Eq. (56) is inferred from the charge density
$\rho(\bm{r},t)$ via the 1D charge conservation equation
$\partial_{t}\rho+\partial_{z}j=0$. As described in Ref. 36, the density
profile of the expanding plasma exhibits two charge-separation regions: one
around the mobile ion front ($z\simeq z_{\rm f}$) with negative areal charge
density $-\sigma_{\rm f}(t)$ and one around the initial plasma surface
($z\simeq 0$) with positive areal charge density $+\sigma_{\rm f}(t)$. We will
neglect the longitudinal extent of those nonneutral layers and model
$\rho(\bm{r},t)$ as a sum of two charged planes, centered at $z\simeq z_{\rm
f}(t)$ and $z=0$,

$\rho(\bm{r},t)=\sigma_{\rm f}(t)\left[\delta(z)-\delta(z-z_{\rm
f}(t))\right]f_{\perp,\textsc{PE}}(\bm{r}_{\perp})\,.$ (65)

Here $f_{\perp,\textsc{PE}}(\bm{r}_{\perp})$ describes the transverse density
profile of the expanding plasma. In principle, this profile should vary in
time since the plasma expands both in the longitudinal and transverse
directions. For simplicity, we consider a constant transverse distribution of
Gaussian shape and same width as that of the recirculating electron beam:

$f_{\perp,\textsc{PE}}(\bm{r}_{\perp})=e^{-4\ln 2(r_{\perp}/w_{h}(d))^{2}}\,.$
(66)

In a 1D geometry, the time-varying areal charge density is linked to the
sheath field through $\sigma_{\rm f}(t)=\varepsilon_{0}E_{\rm f,1D}(t)$. For
simplicity, we assume that this relation remains approximately valid in the 3D
expansion regime, i.e.,

$\sigma_{\rm f}(t)=\varepsilon_{0}E_{\rm f,3D}(t)\,,$ (67)

with $E_{\rm f,3D}(t)$ as defined by Eq. (63). This leads to the current
density

$j(\bm{r},t)=j_{\parallel}(z,t)f_{\perp,\textsc{PE}}(\bm{r}_{\perp})\,,$ (68)

where

$\displaystyle j_{\parallel}(z,t)$
$\displaystyle=-\int_{-\infty}^{z}\frac{\partial\rho(\bm{r},t)}{\partial
t}\mathrm{d}z$ $\displaystyle=-\dot{\sigma}_{\rm
f}(t)\left[\mathcal{H}_{0}-\mathcal{H}_{z_{\rm f}(t)}\right]-\sigma_{\rm
f}(t)\dot{z}_{\rm f}(t)\delta_{z_{\rm f}(t)}\,.$ (69)

Here, $\mathcal{H}_{z_{\rm f}(t)}\equiv\mathcal{H}\left(z-z_{\rm f}(t)\right)$
and $\delta_{z_{\rm f}(t)}\equiv\delta\left(z-z_{\rm f}(t)\right)$ denote,
respectively, the Heaviside and Dirac delta functions centered on $z=z_{\rm
f}(t)$. Finally, introducing the transverse spatial Fourier transform

$\displaystyle\mathcal{F}[f_{\perp,\textsc{PE}}](k_{x},k_{y})$
$\displaystyle=\iint_{-\infty}^{\infty}\mathrm{d}x\mathrm{d}y\,f_{\perp,\textsc{PE}}(x,y)$
$\displaystyle\times e^{-2i\pi(k_{x}x+k_{y}y)}$
$\displaystyle=\pi\frac{w_{h}^{2}(d)}{4\ln
2}e^{-\pi^{2}\frac{w_{h}^{2}(d)}{4\ln 2}(k_{x}^{2}+k_{y}^{2})}\,,$ (70)

we can rewrite Eq. (56) as

$\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
PER}}{\partial\nu\partial\Omega}(\theta,\nu)=\frac{w_{h}^{4}(d)\sin^{2}\theta}{2(8\ln
2)^{2}\varepsilon_{0}c^{3}}e^{-\left(w_{h}(d)\sin\theta\frac{\pi\nu}{c}\right)^{2}/2\ln
2}$ $\displaystyle\times\left|\int\mathcal{F}\left[\frac{\partial
j_{\parallel}}{\partial
t^{\prime}}\right]\left(\frac{-\nu\cos\theta}{c}\right)e^{-2i\pi\nu
t^{\prime}}\mathrm{d}t^{\prime}\right|^{2}\,.$ (71)

where $\mathcal{F}\left[\frac{\partial j_{\parallel}}{\partial
t^{\prime}}\right]$ is the spatial (along $z$) Fourier transform of the time
derivative $\partial j_{\parallel}/\partial t$, evaluated analytically at
$k_{z}=-\nu\cos\theta/c$ as a function of $z_{\rm f}(t)$ and $\sigma_{\rm
f}(t)$ (See Appendix B). The time integral in the above expression is then
computed numerically.

In the following sections, we study the variations in the energy yield and
spectra of PER with the main laser-plasma parameters.

###  IV.4 Influence of the laser focusing

We first perform a scan over the laser pulse intensity (or, equivalently,
degree of laser focusing) at fixed laser pulse energy
($\mathcal{E}_{L}=2.8\,\rm J$) and duration ($\tau_{L}=30\,\rm fs$), by
keeping $a_{L}w_{L}$ constant. Figure 10 plots the evolution of the energy
radiated via PER ($\mathcal{E}_{\rm PER}$, integrated over $0.1\leq\nu\leq
40\,\rm THz$) as a function of $a_{L}$ (in the range $5\leq a_{L}\leq 15$) and
for different target thicknesses ($2\leq d\leq 20\,\rm\mu m$). As in Sec. III,
the laser-to-hot-electron conversion efficiency is set to $\eta_{\rm h}=0.2$.

First, one observes that the predicted yields lie in the $3\times
10^{-5}-1.7\times 10^{-4}\,\mathcal{E}_{\rm beam}$ range, _i.e._ , at least
one order of magnitude below the previous estimates of the energy radiated by
the fast electrons alone. All curves show a monotonic increase with $a_{L}$.
In detail, for $d=2\,\rm\mu m$, $\mathcal{E}_{\rm PER}$ is enhanced by a
factor of $\sim 2.5$ when $a_{L}$ is increased from 5 to 15. When $d=5\,\rm\mu
m$ ($d=10\,\rm\mu m$), the enhancement factor is of $\sim 2.2$ (resp. $\sim
1.7$).

Figure 10: THz energy yield (integrated over $0.1-40\,\rm THz$) of PER as a
function of the laser field strength $a_{L}$ (right axis), for various target
thicknesses $d$. The product $a_{L}w_{L}$ is kept fixed to ensure a constant
laser pulse energy $\mathcal{E}_{L}=2.8\,\rm J$ (or, equivalently, a constant
electron beam energy $\mathcal{E}_{\rm beam}=0.56\,\rm J$). The blue dashed
curve plots $w_{L}$ (left axis). Apart from $a_{L}$ and $w_{L}$, the laser-
plasma parameters are those used in Fig. 3.

As displayed in Fig. 11, intensifying the laser field (from $a_{L}=5$ to
$a_{L}=15$) by narrowing its spot size (from $w_{L}=15\,\rm\mu m$ to
$5\,\rm\mu m$) broadens and shifts the THz spectrum to higher THz frequencies
(from $\sim 5-15\,\rm THz$ to $\sim 10-25\,\rm THz$) because of a faster
plasma expansion. This effect is reinforced by the accompanying decrease in
transverse width of the sheath field $w_{h}(d)$ [see Eq. (43)], which weakens
the attenuation factor associated with the Fourier transform of $f_{\perp,\rm
PE}$ in Eq. (IV.3). The two spectra, however, are maximized at large angles
$\theta\sim 70-90^{\circ}$, as is typical for nonrelativistic dipole-like
radiation.

Figure 11: THz energy yield (integrated over $0.1-40\,\rm THz$) of PER for
laser field strengths (a) $a_{L}=5$ and (b) $a_{L}=15$. The product
$a_{L}w_{L}$ is kept fixed to ensure a constant laser pulse energy
$\mathcal{E}_{L}=2.8\,\rm J$ (or, equivalently, a constant electron beam
energy $\mathcal{E}_{\rm beam}=0.56\,\rm J$). Apart from $a_{L}$ and $w_{L}$,
the laser-plasma parameters are those used in Fig. 3.

###  IV.5 Influence of the target thickness

We now vary the target thickness $d$ at fixed laser parameters ($a_{L}=15$,
$w_{L}=5\,\rm\mu m$, $\tau_{L}=30\,\rm fs$, $\mathcal{E}_{L}=2.8\,\rm J$).
Figure 12 plots, as a function of $d\geq 2\,\rm\mu m$, the maximum kinetic
energy reached by the fastest protons as well as the corresponding THz energy
yield. The general trend is that of a monotonic decrease in the ion kinetic
energy with the target thickness, from $\mathcal{E}_{\rm f}\simeq 15\,\rm MeV$
at $d=2\,\rm\mu m$ down to $\simeq 5\,\rm MeV$ at $d=20\,\rm\mu m$. These
results are consistent with the measurements reported in Refs. 62, 66, 67.

Figure 12: Kinetic energy of the fastest ions (dashed red curve) and energy
yield of PER integrated over $0.1-40\,\rm THz$ (solid blue curve) as a
function of the target thickness $d\geq 2\,\rm\mu m$. Apart from $d$, the
laser-plasma parameters are those used in Fig. 3 ($a_{L}=15$).  Figure 13:
Energy-angle spectra of PER for (a) $5\,\rm\mu m$, (b) $10\,\rm\mu m$ and (c)
$20\,\rm\mu m$ thick targets. The laser parameters are those of Fig. 3.

Closely following that trend, the model predicts a continuous drop in the
radiated energy from $\mathcal{E}_{\rm PER}\simeq 1.7\times
10^{-4}\,\mathcal{E}_{\rm beam}\simeq 95\,\rm\mu J$ at $d=2\,\rm\mu m$ down to
$\leavevmode\nobreak\ 3.4\times 10^{-5}\,\mathcal{E}_{\rm beam}\simeq
19\,\rm\mu J$ at $d\simeq 20\,\rm\mu m$. Figure 10 further shows that this
decrease in $\mathcal{E}_{\rm PER}$ at larger $d$ is less pronounced when
$a_{L}$ is reduced: $\mathcal{E}_{\rm PER}$ indeed drops by $\sim 5\times$ at
$a_{L}=15$ and by $\sim 3\times$ at $a_{L}=5$.

Modifying the thickness of the foil also strongly impacts the frequency-angle
spectrum of PER. Figures 13(a)-(c) display the spectra obtained for
$d=5\leavevmode\nobreak\ \rm\mu m$, $10\leavevmode\nobreak\ \rm\mu m$ and
$20\leavevmode\nobreak\ \rm\mu m$, respectively. The $d=2\,\rm\mu m$ case is
depicted in Fig. 11(b). We notice that, as the target is made thicker, the
spectrum, initially distributed in a wide frequency range ($5\,\rm
THz\lesssim\nu\lesssim 30\,\rm THz$), shrinks to a much narrower bandwidth
while shifting to smaller frequencies ($1\,\rm THz\lesssim\nu\lesssim 6\,\rm
THz$). As expected, though, its angular profile remains unchanged, with a
broad maximum around $\sim 70-90^{\circ}$.

A major prediction of our modeling is that the THz energy radiated from the
plasma expansion should be about one to three orders of magnitude lower than
that from the fast electrons only. The dominance of the latter is expected to
be particularly marked for a fraction of escaping electrons $\eta_{\rm
esc}\gtrsim 5\%$ (compare Figs. 5 and 12).

##  V Conclusions

In summary, we have developed a novel theoretical model of far-field THz
emissions in relativistic ultrashort laser-foil interactions, assuming a two-
stage scenario. In the first stage, we consider the radiation resulting from
the fast electrons alone at the target backside, during their first and brief
(lasting only a few fs) excursion into vacuum. In this phase, the vast
majority of them are reflected into the target by the sheath field they have
themselves created. The THz emission then originates from the combination of
coherent transition/synchrotron radiations from the recirculating electrons,
as well as from coherent transition radiation from the escaping electrons. In
the second stage, all of the confined fast electrons are assumed to have
relaxed to a thermal distribution that drives, through the sheath field, the
acceleration of the target ions. The dynamics of the nonneutral layers at the
edges of the expanding plasma then gives rise to a dipole-type radiation.

To our knowlege, our unified kinetic model of CTR and CSR, based on the image
charge method and integrated over an ensemble of single-particle trajectories,
is the first of its kind. Its predictions highlight the critical importance of
treating simultaneously these two mechanisms, which tend to interfere
destructively in the THz domain. Their degree of imperfect cancellation, which
determines the energy yield and spectrum of the net (CTSR) radiation, depends
on the details of the electron trajectories, notably the extent of their
excursion in vacuum – a function of the electron momentum and sheath field
strength – and the amount of escaping (assumed ballistic) electrons.

We have examined the sensitivity of this radiation to the main laser-target
parameters, varying these around a reference setup characterized by
$a_{L}=15$, $\tau_{L}=30\,\rm fs$, $w_{L}=5\,\rm\mu m$ and $\eta_{h}=0.2$,
corresponding to a (fixed) electron beam energy of $\mathcal{E}_{\rm
beam}=0.56\,\rm J$. An important finding is that a fraction of escaping
electrons as low as $\eta_{\rm esc}\simeq 1\,\%$ – a value consistent with
expectations [34] – may suffice to shape the THz spectrum radiated from thin
($d=2\,\rm\mu m$) foils. At this threshold value, the spectrum is shifted to
relatively low frequencies ($\nu\lesssim 5\,\rm Thz$) compared to the case of
full electron refluxing, though containing about the same integrated energy
($\simeq 4\times 10^{-3}\,\mathcal{E}_{\rm beam}\simeq 2\,\rm mJ$). The
latter, however, is predicted to rise tenfold or more should the escaping
electron fraction reach $\sim 5\%$. The minimum fraction of escaping electrons
needed to govern the radiation is an increasing function of the target
thickness and degree of laser focusing. For a $20\,\rm\mu m$ thick foil and a
$5\,\rm\mu m$ laser spot size ($a_{L}=15$), this threshold fraction is found
to lie between $5\%$ and $10\,\%$, quite a bit above the expected amount
($\sim 1\,\%$) of those electrons [34]. Consequently, the contribution of the
confined electrons, characterized by higher frequencies ($\nu\simeq 5-20\,\rm
THz$), may then well dominate the total spectrum when using not-so-thin
targets and tightly focused laser pulses. Near the threshold between the two
regimes, the THz spectrum will exhibit two distinct bands at low and high
frequencies. More generally, the angular distribution of the radiation is
broad and tends to peak outside the emission cone of the fast electrons,
especially in thin foils or when CTR from escaping electrons prevails.

Of course, those trends should be considered qualitative as they depend on a
number of coupled parameters, which themselves depend on the detailed
experimental conditions. In particular, all throughout, we have set to
$20\,\%$ the fraction of the laser pulse energy carried by the fast electrons,
and to $30^{\circ}$ the angular spread of the latter, even though those
quantities are likely affected by the laser intensity and spot size. Moreover,
to ensure the tractability of our already computationally heavy model, we have
assumed that the fast electrons propagate ballistically across the solid foil
(which implies a small enough target thickness) and that the sheath field that
confines most of them is both uniform and stationary. Finally, our model being
restricted to the first excursion of the fast electrons into vacuum, it
discards the subsequent, possibly significant, CTSR-type bursts as the
electrons recirculate through the foil and, in so doing, progressively
decelerate and scatter.

These limitations notwithstanding, our modeling indicates that when CTR and
CSR insufficiently compensate each other, as occurs for a large enough
fraction of escaping electrons or a weak enough sheath field, the radiated
energy may approach, or even exceed, the electron beam energy. This result
suggests that, under such conditions, one should describe the self-consistent
effect of the collective radiation on the electron dynamics. This complex
problem can only be quantitatively addressed through three-dimensional, fully
kinetic numerical simulations, ideally equipped with a far-field radiation
diagnostic [68].

Our model for the THz radiation (PER) emitted at later times as a result of
ion acceleration no longer relies on a kinetic treatment of the accelerated
particles, because of the daunting complexity of the problem, but rather on a
simplified description of the space- and time-varying current density
distribution in the expanding plasma. Specifically, we approximate the charge-
separation regions to two infinitely thin disks of opposite areal charge,
located at the moving ion front and the initial target surface, and describe
their dynamics by adapting the plasma expansion model proposed in Ref. 64,
itself built upon several previous works [36, 63, 65]. Compared to previous
attempts at estimating the THz emission from accelerated ions [26, 27, 28,
29], our description takes account of the time-varying areal charge of the
nonneutral layers and goes beyond the customarily considered isothermal
electron limit by considering the transition to the adiabatic regime, as is
relevant to micron-thick targets.

Our computations, conducted over parameter ranges similar to those considered
for CTSR – notably assuming the same driving laser pulses –, predict that the
THz energy radiated during the plasma expansion phase should be at least one
order of magnitude below that emitted earlier by the fast electrons. Overall,
for the interaction conditions considered, the beam-to-THz conversion
efficiency of PER varies in the $\sim 10^{-5}-10^{-4}$ range and is a
decreasing function of the target thickness. It also rises with the degree of
laser focusing, albeit more and more slowly as the foil is made thicker (for
$2\leq d\leq 20\,\rm\mu m$). As expected, PER produces increasingly low
frequencies, and over an increasingly narrow range, when the ion acceleration
diminishes, that is, when the target thickness or the laser intensity are
reduced. This radiation is mainly emitted at large angles ($>50^{\circ}$), yet
this feature alone should not suffice to allow PER to be discerned over the
intense background due to CTSR.

Our results evidently call for numerical and experimental validations.
Regarding the latter, it should be noted that most experiments to date have
characterized the THz emissions over restricted frequency and angular ranges
[25, 26]. A notable exception is Ref. 58, where the full THz distribution was
diagnosed but with an obliquely incident laser pulse and where no variation in
the target thickness was reported. Recording detailed frequency-angle spectra
that can be compared with the theoretical spectra would be highly desirable in
order to better tune the model parameters and support its predictions.

##  Appendix A Discrete to continuous description of the energy spectrum

We detail here the calculation steps between Eq. (19) and Eq. (II.4) to obtain
a continuous description of the energy spectrum radiated by a set of charged
particles obeying a given distribution function. We start from the discrete
expression of the energy spectrum [see Eq. (19)] and introduce, for each
particle (labeled by the subscript $l$), the amplitude $\bm{A}_{l}$ and the
complex phase $\phi_{l}$ such that

$\displaystyle\frac{\partial^{2}\mathcal{I}_{\rm
rad}}{\partial\nu\partial\Omega}(\theta,\nu)$
$\displaystyle=\left\langle\frac{q^{2}\nu^{2}}{2\varepsilon_{0}c}\Biggl{|}\sum_{l=1}^{N_{h}}\left[\bm{A}_{l}^{-}(\bm{\hat{r}},\nu)-\bm{A}_{l}^{+}(\bm{\hat{r}},\nu)\vphantom{\int}\right]\Biggr{|}^{2}\right\rangle$
$\displaystyle=\left\langle\frac{q^{2}\nu^{2}}{2\varepsilon_{0}c}\Biggl{|}\sum_{l=1}^{N_{h}}\bm{A}_{l}e^{i\phi_{l}}\Biggr{|}^{2}\right\rangle\,,$
(72)

where $\langle\cdot\rangle$ indicates an ensemble average. We then introduce a
continuous description $\bm{A}_{\bm{r},t,\bm{u}}(\bm{\hat{r}},\nu)$ and
$\phi_{\bm{r},t,\bm{u}}(\bm{\hat{r}},\nu)$ of the terms $\bm{A}_{l}$ and
$\phi_{l}$ by imposing that for every particle,

$\bm{A}_{l}e^{i\phi_{l}}=\bm{A}_{\bm{r}_{l},t_{l},\bm{u}_{l}}(\bm{\hat{r}},\nu)\exp\left(i\phi_{\bm{r}_{l},t_{l},\bm{u}_{l}}(\bm{\hat{r}},\nu)\right)\,,$
(73)

where $\bm{r}_{l}$, $t_{l}$ and $\bm{u}_{l}$ represent the position, time and
initial momentum of the accelerated particle. The sum of this quantity over
$N_{h}$ particles sampling the $(\bm{r},t,\bm{u})$ space can be expressed as

$\displaystyle\Biggl{\langle}\Biggl{|}\sum_{l=1}^{N_{h}}\bm{A}_{l}e^{i\phi_{l}}\Biggr{|}^{2}\Biggr{\rangle}$
$\displaystyle=\Biggl{\langle}\Biggl{|}\,\iiint_{\bm{r},t,\bm{u}}\bm{A}_{\bm{r},t,\bm{u}}(\bm{\hat{r}},\nu)e^{i\phi_{\bm{r},t,\bm{u}}(\bm{\hat{r}},\nu)}\times\Biggl{(}\sum_{l=1}^{N_{h}}\delta(\bm{r}-\bm{r}_{l})\delta(t-t_{l})\delta(\bm{u}-\bm{u}_{l})\Biggr{)}\,\mathrm{d}^{2}\bm{r}\,\mathrm{d}t\,\mathrm{d}^{3}\bm{u}\Biggr{|}^{2}\Biggr{\rangle}\,.$
(74)

Next, we introduce the distribution function $F(\bm{r},t,\bm{u})$ that
fulfills

$\Biggl{\langle}\sum_{l=1}^{N_{h}}\delta(\bm{r}-\bm{r}_{l})\delta(t-t_{l})\delta(\bm{u}-\bm{u}_{l})\Biggr{\rangle}\to
N_{h}F(\bm{r},t,\bm{u})\,.$ (75)

We hereafter assume that the distribution function can be taken in the form of
a separable function, $F(\bm{r},t,\bm{u})=f(\bm{r})h(t)g(\bm{u})$.

Finally, we substitute Eq. (75) into Eq. (74) and make use of
$\bm{A}_{\bm{r},t,\bm{u}}(\bm{\hat{r}},\nu)\exp\left(i\phi_{\bm{r},t,\bm{u}}(\bm{\hat{r}},\nu)\right)=\bm{A}^{-}(\bm{\hat{r}},\nu)-\bm{A}^{+}(\bm{\hat{r}},\nu)\vphantom{\int}$
to obtain the expression of the energy radiated coherently by an ensemble of
$N_{h}\gg 1$ electrons,

$\frac{\partial^{2}\mathcal{I}_{\rm
rad}}{\partial\nu\partial\Omega}(\theta,\nu)=\frac{q^{2}\nu^{2}}{2\varepsilon_{0}c}N_{h}^{2}\Biggl{|}\,\iiint_{\bm{r},t,\bm{u}}f(\bm{r})h(t)g(\bm{u})\\\
\times\left[\bm{A}^{-}(\bm{\hat{r}},\nu)-\bm{A}^{+}(\bm{\hat{r}},\nu)\vphantom{\int}\right]\,\mathrm{d}^{2}\bm{r}\,\mathrm{d}t\,\mathrm{d}^{3}\bm{u}\Biggr{|}^{2}\,.$
(76)

Equation (76) is identical to Eq. (II.4).

##  Appendix B Spatial Fourier transform of $\partial_{t}j_{\parallel}$

The $\rm 1D$ current density defined by Eq. (IV.3) involves several Dirac
delta functions. Its time derivative can be written as a sum of three
components:

$\displaystyle\partial j_{\parallel}(z,t)/\partial t=$
$\displaystyle-\ddot{\sigma}(t)\left[\mathcal{H}_{0}-\mathcal{H}_{z_{\rm
f}(t)}\right]$ $\displaystyle-\left[2\dot{\sigma}(t)\dot{z}_{\rm
f}(t)+\sigma(t^{\prime})\ddot{z}_{\rm f}(t)\right]\delta_{z_{\rm f}(t)}$
$\displaystyle+\sigma(t)\dot{z}^{2}_{f}(t)\delta^{\prime}_{z_{\rm f}(t)}\,,$
(77)

where $\partial\delta(z-z_{\rm f}(t))/\partial t=-\dot{z}_{\rm
f}(t)\delta^{\prime}(z-z_{\rm f}(t))$.

The spatial Fourier transform of $\partial j_{\parallel}(z,t)/\partial t$ that
appears in the expression of the PER spectrum, Eq. (IV.3), can then be
evaluated analytically upon noting that each of the above components admits a
closed-form Fourier transform:

$\displaystyle\mathcal{F}\left[\partial j_{\parallel}(z,t)/\partial
t\right](k_{z})=$ $\displaystyle\hskip
30.0pt-\ddot{\sigma}(t)\mathcal{F}\left[\mathcal{H}_{0}-\mathcal{H}_{z_{\rm
f}(t)}\right](k_{z})$ $\displaystyle\hskip
30.0pt-\left[2\dot{\sigma}(t)\dot{z}_{\rm
f}(t)+\sigma(t^{\prime})\ddot{z}_{\rm
f}(t)\right]\mathcal{F}\left[\delta_{z_{\rm f}(t)}\right](k_{z})$
$\displaystyle\hskip 30.0pt+\sigma(t)\dot{z}^{2}_{\rm
f}(t)\mathcal{F}\left[\delta^{\prime}_{z_{\rm f}(t)}\right](k_{z})\,,$ (78)

where

$\displaystyle\mathcal{F}\left[\mathcal{H}_{0}-\mathcal{H}_{z_{\rm
f}(t)}\right](k_{z})=$ $\displaystyle\hskip 20.0ptz_{\rm
f}(t^{\prime})\,\mathrm{sinc}\left(k_{z}z_{\rm
f}(t^{\prime})\right)\exp\Bigl{(}-2i\pi k_{z}z_{\rm f}(t^{\prime})/2\Bigr{)},$
(79) $\displaystyle\mathcal{F}\left[{\delta_{z_{\rm
f}(t)}(z)}\right](k_{z})=\exp\Bigl{(}-2i\pi k_{z}z_{\rm
f}(t^{\prime})\Bigr{)},$ (80)
$\displaystyle\mathcal{F}\left[{\delta^{\prime}_{z_{\rm
f}(t)}(z)}\right](k_{z})=2i\pi k_{z}\exp\Bigl{(}-2i\pi k_{z}z_{\rm
f}(t^{\prime})\Bigr{)}\,.$ (81)

Equation (B) is then injected into Eq. (IV.3) where it is evaluated at the
wavenumber $k_{z}=-\nu\cos\theta/c$.

## References

  * Tonouchi [2007] M. Tonouchi, “Cutting-edge terahertz technology,” Nat. Photon. 1, 97–105 (2007).
  * Pickwell and Wallace [2006] E. Pickwell and V. P. Wallace, “Biomedical applications of terahertz technology,” J. Phys. D 39, R301–R310 (2006).
  * Mittleman  _et al._ [1999] D. Mittleman, M. Gupta, R. Neelamani, R. Baraniuk, J. Rudd, and M. Koch, “Recent advances in terahertz imaging,” Appl. Phys. B 68, 1085–1094 (1999).
  * 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,” Eur. Phys. Lett. 126, 24001 (2019).
  * Chan, Deibel, and Mittleman [2007] W. L. Chan, J. Deibel, and D. M. Mittleman, “Imaging with terahertz radiation,” Rep. Prog. Phys. 70, 1325 (2007).
  * Orenstein and Millis [2000] J. Orenstein and A. Millis, “Advances in the physics of high-temperature superconductivity,” Science 288, 468–474 (2000).
  * Kampfrath, Tanaka, and Nelson [2013] T. Kampfrath, K. Tanaka, and K. A. Nelson, “Resonant and nonresonant control over matter and light by intense terahertz transients,” Nat. Photon. 7, 680–690 (2013).
  * Li  _et al._ [2019] X. Li, T. Qiu, J. Zhang, E. Baldini, J. Lu, A. M. Rappe, and K. A. Nelson, “Terahertz field-induced ferroelectricity in quantum paraelectric SrTiO3,”  364, 1079–1082 (2019).
  * Andreasen and Kolesik [2012] J. Andreasen and M. Kolesik, “Nonlinear propagation of light in structured media: Generalized unidirectional pulse propagation equation,” Phys. Rev. E 86, 036706 (2012).
  * Nanni  _et al._ [2015] E. A. Nanni, W. R. Huang, K.-H. Hong, K. Ravi, A. Fallahi, G. Moriena, R. J. Dwayne Miller, and F. X. Kärtner, “Terahertz-driven linear electron acceleration,” Nat. Commun. 6, 8486 (2015).
  * Novelli, Guchhait, and Havenith [2020] F. Novelli, B. Guchhait, and M. Havenith, “Towards Intense THz Spectroscopy on Water: Characterization of Optical Rectification by GaP, OH1, and DSTMS at OPA Wavelengths,” Materials 13, 1311 (2020).
  * Vella  _et al._ [2021] A. Vella, J. Houard, L. Arnoldi, M. Tang, M. Boudant, A. Ayoub, A. Normand, G. Da Costa, and A. Hideur, “High-resolution terahertz-driven atom probe tomography,” Sci. Adv. 7, eabd7259 (2021).
  * Hebling  _et al._ [2002] J. Hebling, G. Almasi, I. Z. Kozma, and J. Kuhl, “Velocity matching by pulse front tilting for large-area THz-pulse generation,” Opt. Express 10, 1161–1166 (2002).
  * Vicario, Monoszlai, and Hauri [2014] 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 (2014).
  * Yeh  _et al._ [2007] K. L. Yeh, M. C. Hoffmann, J. Hebling, and K. A. Nelson, “Generation of 10 $\mu$J ultrashort THz pulses by optical rectification,” Appl. Phys. Lett. 90, 171121 (2007).
  * 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 interactions,” Nat. Phys. 2, 605 (2008).
  * Clerici  _et al._ [2013] M. Clerici, M. Peccianti, B. E. Schmidt, L. Caspani, M. Shalaby, M. Giguère, A. Lotti, A. Couairon, F. Légaré, T. Ozaki, D. Faccio, and R. Morandotti, “Wavelength scaling of terahertz generation by gas ionization,” Phys. Rev. Lett. 110, 253901 (2013).
  * A. Nguyen and P. González de Alaiza Martínez and I. Thiele and S. Skupin and L. Bergé [2018] A. Nguyen and P. González de Alaiza Martínez and I. Thiele and S. Skupin and L. Bergé, “Broadband terahertz radiation from two-color mid- and far-infrared laser filaments in air,” Phys. Rev. A 97, 063839 (2018).
  * Leemans  _et al._ [2003] W. P. Leemans, C. Geddes, J. Faure, C. Tóth, J. van Tilborg, C. B. Schroeder, E. Esarey, G. Fubiani, D. Auerbach, B. Marcelis, M. A. Carnahan, R. A. Kaindl, J. Byrd, and M. C. Martin, “Observation of terahertz emission from a laser-plasma accelerated electron bunch crossing a plasma-vacuum boundary,” Phys. Rev. Lett. 91, 074802 (2003).
  * Déchard  _et al._ [2018] J. Déchard, A. Debayle, X. Davoine, L. Gremillet, and L. Bergé, “Terahertz Pulse Generation in Underdense Relativistic Plasmas: From Photoionization-Induced Radiation to Coherent Transition Radiation,” Phys. Rev. Lett. 120, 144801 (2018).
  * Schroeder  _et al._ [2004] C. B. Schroeder, E. Esarey, J. van Tilborg, and W. P. Leemans, “Theory of coherent transition radiation generated at a plasma-vacuum interface,” Phys. Rev. E 69, 016501 (2004).
  * Ding, Sheng, and Koh [2013] W. J. Ding, Z. M. Sheng, and W. S. Koh, “High-field half-cycle terahertz radiation from relativistic laser interaction with thin solid targets,” Appl. Phys. Lett. 103, 204107 (2013).
  * Ding and Sheng [2016] W. J. Ding and Z. M. Sheng, “Sub GV/cm terahertz radiation from relativistic laser-solid interactions via coherent transition radiation,” Phys. Rev. E 93, 063204 (2016).
  * Liao  _et al._ [2016] G.-Q. Liao, Y.-T. Li, Y.-H. Zhang, H. Liu, X.-L. Ge, S. Yang, W.-Q. Wei, X.-H. Yuan, Y.-Q. Deng, B.-J. Zhu, Z. Zhang, W.-M. Wang, Z.-M. Sheng, L.-M. Chen, X. Lu, J.-L. Ma, X. Wang, and J. Zhang, “Demonstration of Coherent Terahertz Transition Radiation from Relativistic Laser-Solid Interactions,” Phys. Rev. Lett. 116, 205003 (2016).
  * Jin  _et al._ [2016] Z. Jin, H. B. Zhuo, T. Nakazawa, J. H. Shin, S. Wakamatsu, N. Yugami, T. Hosokai, D. B. Zou, M. Y. Yu, Z. M. Sheng, and R. Kodama, “Highly efficient terahertz radiation from a thin foil irradiated by a high-contrast laser pulse,” Phys. Rev. E 94, 033206 (2016).
  * Herzer  _et al._ [2018] S. Herzer, A. Woldegeorgis, J. Polz, A. Reinhard, M. Almassarani, B. Beleites, F. Ronneberger, R. Grosse, G. G. Paulus, U. Hübner, T. May, and A. Gopal, “An investigation on THz yield from laser-produced solid density plasmas at relativistic laser intensities,” New J. Phys. 20, 063019 (2018).
  * Woldegeorgis  _et al._ [2019] A. Woldegeorgis, S. Herzer, M. Almassarani, S. Marathapalli, and A. Gopal, “Modeling terahertz emission from the target rear side during intense laser-solid interactions,” Physical Review E 100, 053204 (2019).
  * Liao  _et al._ [2019] G. Liao, Y. Li, H. Liu, G. G. Scott, D. Neely, Y. Zhang, B. Zhu, Z. Zhang, C. Armstrong, E. Zemaityte, _et al._ , “Multimillijoule coherent terahertz bursts from picosecond laser-irradiated metal foils,” Proc. Natl. Acad. Sci 116, 3994–3999 (2019).
  * Liao  _et al._ [2020] G.-Q. Liao, H. Liu, G. G. Scott, Y.-H. Zhang, B.-J. Zhu, Z. Zhang, Y.-T. Li, C. Armstrong, E. Zemaityte, P. Bradford, D. R. Rusby, D. Neely, P. G. Huggard, P. McKenna, C. M. Brenner, N. C. Woolsey, W.-M. Wang, Z.-M. Sheng, and J. Zhang, “Towards Terawatt-Scale Spectrally Tunable Terahertz Pulses via Relativistic Laser-Foil Interactions,” Phys. Rev. X 10, 031062 (2020).
  * Déchard  _et al._ [2020] J. Déchard, X. Davoine, L. Gremillet, and L. Bergé, “Terahertz emission from submicron solid targets irradiated by ultraintense femtosecond laser pulses,” Phys. Plasmas 27, 093105 (2020).
  * Kemp  _et al._ [2014] A. J. Kemp, F. Fiuza, A. Debayle, T. Johzaki, W. B. Mori, P. K. Patel, Y. Sentoku, and L. O. Silva, “Laser-plasma interactions for fast ignition,” Nucl. Fusion 54, 054002 (2014).
  * Fill [2001] E. E. Fill, “Analytical theory of pulsed relativistic electron beams entering a vacuum,” Phys. Plasmas 8, 4613–4617 (2001).
  * Link  _et al._ [2011] A. Link, R. R. Freeman, D. W. Schumacher, and L. D. Van Woerkom, “Effects of target charging and ion emission on the energy spectrum of emitted electrons,” Phys. Plasmas 18, 053107 (2011).
  * Rusby  _et al._ [2019] D. R. Rusby, C. D. Armstrong, G. G. Scott, M. King, P. McKenna, and D. Neely, “Effect of rear surface fields on hot, refluxing and escaping electron populations via numerical simulations,” High Power Laser Sci. Eng. 7, e45 (2019).
  * Wilks  _et al._ [2001] S. C. Wilks, A. B. Langdon, T. E. Cowan, M. Roth, M. Singh, S. Hatchett, M. H. Key, D. Pennington, A. MacKinnon, and R. A. Snavely, “Energetic proton generation in ultra-intense laser-solid interactions,” Phys. Plasmas 8, 542 (2001).
  * Mora [2003] P. Mora, “Plasma Expansion into a Vacuum,” Phys. Rev. Lett. 90, 185002 (2003).
  * Macchi, Borghesi, and Passoni [2013] A. Macchi, M. Borghesi, and M. Passoni, “Ion acceleration by superintense laser-plasma interaction,” Rev. Mod. Phys. 85, 751–793 (2013).
  * Gopal  _et al._ [2013a] A. Gopal, S. Herzer, A. Schmidt, P. Singh, A. Reinhard, W. Ziegler, D. Brömmel, A. Karmakar, P. Gibbon, U. Dillner, _et al._ , “Observation of gigawatt-class thz pulses from a compact laser-driven particle accelerator,” Phys. Rev. Lett. 111, 074802 (2013a).
  * Gopal  _et al._ [2013b] A. Gopal, P. Singh, S. Herzer, A. Reinhard, A. Schmidt, U. Dillner, T. May, H.-G. Meyer, W. Ziegler, and G. G. Paulus, “Characterization of 700 $\mu$J T rays generated during high-power laser solid interaction,” Opt. Lett. 38, 4705 (2013b).
  * Ter-Mikaelian [1972] M. L. Ter-Mikaelian,  _High-energy electromagnetic processes in condensed media_, 29 (John Wiley & Sons, 1972).
  * Ginzburg [1982] V. L. Ginzburg, “Transition radiation and transition scattering,” Phys. Scr. T2A, 182 (1982).
  * Zheng  _et al._ [2003] J. Zheng, K. A. Tanaka, T. Miyakoshi, Y. Kitagawa, R. Kodama, T. Kurahashi, and T. Yamanaka, “Theoretical study of transition radiation from hot electrons generated in the laser-solid interaction,” Phys. Plasmas 10, 2994 (2003).
  * Bellei  _et al._ [2012] C. Bellei, J. R. Davies, P. K. Chauhan, and Z. Najmudin, “Coherent transition radiation in relativistic laser-solid interactions,” Plasma Phys. Control. Fusion 54, 035011 (2012).
  * Ridgers  _et al._ [2011] C. P. Ridgers, M. Sherlock, R. G. Evans, A. P. L. Robinson, and R. J. Kingham, “Superluminal sheath-field expansion and fast-electron-beam divergence measurements in laser-solid interactions,” pre 83, 036404 (2011).
  * Zhuo  _et al._ [2017] H. B. Zhuo, S. J. Zhang, X. H. Li, H. Y. Zhou, X. Z. Li, D. B. Zou, M. Y. Yu, H. C. Wu, Z. M. Sheng, and C. T. Zhou, “Terahertz generation from laser-driven ultrafast current propagation along a wire target,” Phys. Rev. E 95, 013201 (2017).
  * Jackson [1999] J. D. Jackson,  _Classical Electrodynamics_ (John Wiley & Sons, 1999).
  * Zangwill [2012] A. Zangwill,  _Modern Electrodynamics_ (Cambridge University Press, 2012).
  * Shkvarunets and Fiorito [2008] A. G. Shkvarunets and R. B. Fiorito, “Vector electromagnetic theory of transition and diffraction radiation with application to the measurement of longitudinal bunch size,” Physical Rev. Accel. Beams 11, 012801 (2008).
  * Bolotovskii and Serov [2009] B. M. Bolotovskii and A. V. Serov, “Features of the transition radiation field,” Phys. Usp. 52, 487–493 (2009).
  * Dutt and Rokhlin [1993] A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data,” SIAM Journal on Scientific computing 14, 1368–1393 (1993).
  * Greengard and Lee [2004] L. Greengard and J.-Y. Lee, “Accelerating the nonuniform fast fourier transform,” SIAM Review 46, 443–454 (2004).
  * Ruiz-Antolín and Townsend [2018] D. Ruiz-Antolín and A. Townsend, “A Nonuniform Fast Fourier Transform Based on Low Rank Approximation,” SIAM J. Sci. Comp. 40, A529–A547 (2018).
  * Macchi [2013] A. Macchi,  _A Superintense Laser–Plasma Interaction Theory Primer_ (Springer, 2013).
  * Nilson  _et al._ [2008] P. M. Nilson, W. Theobald, J. Myatt, C. Stoeckl, M. Storm, O. V. Gotchev, J. D. Zuegel, R. Betti, D. D. Meyerhofer, and T. C. Sangster, “High-intensity laser-plasma interactions in the refluxing limit),” Phys. Plasmas 15, 056308 (2008).
  * Carroll  _et al._ [2010] D. C. Carroll, O. Tresca, R. Prasad, L. Romagnani, P. S. Foster, P. Gallegos, S. Ter-Avetisyan, J. S. Green, M. J. V. Streeter, N. Dover, C. A. J. Palmer, C. M. Brenner, F. H. Cameron, K. E. Quinn, J. Schreiber, A. P. L. Robinson, T. Baeva, M. N. Quinn, X. H. Yuan, Z. Najmudin, M. Zepf, D. Neely, M. Borghesi, and P. McKenna, “Carbon ion acceleration from thin foil targets irradiated by ultrahigh-contrast, ultraintense laser pulses,” New J. Phys. 12, 045020 (2010).
  * Debayle  _et al._ [2010] A. Debayle, J. J. Honrubia, E. d’Humières, and V. T. Tikhonchuk, “Divergence of laser-driven relativistic electron beams,” Phys. Rev. E 82, 036405 (2010).
  * Stupakov and Penn [2018] G. Stupakov and G. Penn,  _Classical Mechanics and Electromagnetism in Accelerator Physics_ (Springer, 2018).
  * Gopal  _et al._ [2019] A. Gopal, A. Woldegeorgis, S. Herzer, and M. Almassarani, “Spatiotemporal visualization of the terahertz emission during high-power laser-matter interaction,” Phys. Rev. E 100, 053203 (2019).
  * Kuratov  _et al._ [2016] A. S. Kuratov, A. V. Brantov, Y. M. Aliev, and V. Y. Bychenkovba, “Terahertz radiation in laser-induced charge separation in the irradiated plasma target,” Quantum Electron. 46, 1023–1030 (2016).
  * McDonald [1997] K. T. McDonald, “The relation between expressions for time-dependent electromagnetic fields given by Jefimenko and by Panofsky and Phillips,” Am. J. Phys. 65, 1074–1075 (1997).
  * Panofsky and Phillips [2005] W. K. Panofsky and M. Phillips,  _Classical electricity and magnetism, Second edition_ (Dover Publications, 2005).
  * MacKinnon  _et al._ [2002] A. J. MacKinnon, Y. Sentoku, P. K. Patel, D. W. Price, S. Hatchett, M. H. Key, C. Andersen, R. Snavely, and R. R. Freeman, “Enhancement of Proton Acceleration by Hot-Electron Recirculation in Thin Foils Irradiated by Ultraintense Laser Pulses,” Phys. Rev. Lett. 88, 215006 (2002).
  * Mora [2005] P. Mora, “Thin-foil expansion into a vacuum,” Phys. Rev. E 72, 056401 (2005).
  * Ferri  _et al._ [2018] J. Ferri, L. Senje, M. Dalui, K. Svensson, B. Aurand , M. Hansson, A. Persson, O. Lundh, C. G. Wahlström, L. Gremillet, E. Siminos, T. C. DuBois, L. Yi, J. L. Martins, and T. Fülöp, “Proton acceleration by a pair of successive ultraintense femtosecond laser pulses,” Phys. Plasmas 25, 043115 (2018).
  * Brantov  _et al._ [2015] A. V. Brantov, E. A. Govras, V. Y. Bychenkov, and W. Rozmus, “Ion energy scaling under optimum conditions of laser plasma acceleration from solid density targets,” Phys. Rev. ST Accel. Beams 18, 021301 (2015).
  * Zeil  _et al._ [2010] K. Zeil, S. D. Kraft, S. Bock, M. Bussmann, T. E. Cowan, T. Kluge, J. Metzkes, T. Richter, R. Sauerbrey, and U. Schramm, “The scaling of proton energies in ultrashort pulse laser plasma acceleration,” New J. Phys. 12, 045015 (2010).
  * Poole  _et al._ [2018] P. L. Poole, L. Obst, G. E. Cochran, J. Metzkes, H. P. Schlenvoigt, I. Prencipe, T. Kluge, T. Cowan, U. Schramm, D. W. Schumacher, and K. Zeil, “Laser-driven ion acceleration via target normal sheath acceleration in the relativistic transparency regime,” New J. Phys. 20, 013019 (2018).
  * Pardal  _et al._ [2023] M. Pardal, A. Sainte-Marie, A. Reboul-Salze, R. A. Fonseca, and J. Vieira, “RaDiO: An efficient spatiotemporal radiation diagnostic for particle-in-cell codes,” Comp. Phys. Commun. 285, 108634 (2023).