# Numerical study on femtosecond electro-optical spatial decoding of
transition radiation from laser wakefield accelerated electron bunches

K. Huang<EMAIL_ADDRESS>Kansai Institute for Photon Science (KPSI),
National Institutes for Quantum Science and Technology (QST), 8-1-7 Umemidai,
Kizugawa-city, Kyoto, 619-0215, Japan  Laser Accelerator R&D, Innovative Light
Sources Division, RIKEN SPring-8 Center, 1-1-1 Kouto, Sayo, Hyogo, 679-5148,
Japan  Z. Jin  SANKEN, Osaka University, 8-1 Mihogaoka, Ibaraki, Osaka,
567-0047, Japan  Laser Accelerator R&D, Innovative Light Sources Division,
RIKEN SPring-8 Center, 1-1-1 Kouto, Sayo, Hyogo, 679-5148, Japan  N. Nakanii
Kansai Institute for Photon Science (KPSI), National Institutes for Quantum
Science and Technology (QST), 8-1-7 Umemidai, Kizugawa-city, Kyoto, 619-0215,
Japan  Laser Accelerator R&D, Innovative Light Sources Division, RIKEN
SPring-8 Center, 1-1-1 Kouto, Sayo, Hyogo, 679-5148, Japan  T. Hosokai
SANKEN, Osaka University, 8-1 Mihogaoka, Ibaraki, Osaka, 567-0047, Japan
Laser Accelerator R&D, Innovative Light Sources Division, RIKEN SPring-8
Center, 1-1-1 Kouto, Sayo, Hyogo, 679-5148, Japan  M. Kando  Kansai Institute
for Photon Science (KPSI), National Institutes for Quantum Science and
Technology (QST), 8-1-7 Umemidai, Kizugawa-city, Kyoto, 619-0215, Japan  Laser
Accelerator R&D, Innovative Light Sources Division, RIKEN SPring-8 Center,
1-1-1 Kouto, Sayo, Hyogo, 679-5148, Japan

###### Abstract

This numerical study is focused on electro-optic (EO) spatial decoding of
transition radiation (TR) produced by a relativistic electron bunch passing
through a metal foil. The calculations included the imaging of polychromatic
transition radiation from an electron bunch and the process of EO spatial
decoding. From an experimental perspective, a careful examination of the
calculation approach of the data analysis is essential. Therefore, to
thoroughly understand the process of signal generation and examine the
possibility of adopting a less time-consuming treatment, comparative studies
were conducted on detailed and simplified models of both transition radiation
imaging and EO signal generation. All calculations are defined in SI units for
the convenience of experimental measurements. For TR imaging, the results
suggest that the simplified analytical model is sufficient to perform
polychromatic calculations with considerable accuracy. For EO spatial
decoding, we discussed the process of EO signal generation using 1D and 2D
models. We found that the 1D model was sufficient for rapid data analysis.
Furthermore, the temporal energy chirp was demonstrated to have a minimal
impact on the shape of the EO signal. Because both the transverse and
longitudinal profiles can be calculated with arbitrary distributions, this
numerical study can facilitate measurements of 3D electron charge density
profiles in both laser wakefield acceleration and conventional accelerator
research.

##  I Introduction

Laser wakefield acceleration[1, 2, 3, 4], with an inherent ultrashort
timescale and ultrahigh acceleration gradient, has attracted great interest
worldwide. With the incidence of a high-power femtosecond (fs) laser beam on
an underdense plasma, electron plasma waves with a wavelength of a few tens of
micrometers ($\mu$m) and an electric field strength of $>$ 100 GeV/m can be
stimulated. Trapped electrons achieve energies of over GeV [5, 6, 7, 8, 9, 10]
within a centimeter. In 2021, free-electron lasing at 27 nm driven by an laser
wakefield accelerator (LWFA) was demonstrated via the generation of high-
quality electron bunches [11, 12].

The three-dimensional (3D) charge density distribution of a relativistic
electron bunch is an important parameter for secondary radiation sources of an
LWFA, such as betatron X-rays [13, 14], inverse Compton scattering [15], and
X-ray free electron lasers (XFELs) [16]. The duration of the electron bunch
partially determines the temporal resolution of a pump-probe experiment. The
gain length of lasing processes in an undulator [17] is affected by the 3D
charge density.

Electro-optic (EO) sampling [18] [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
30] was applied to accelerator studies 20 years ago owing to its single-shot
and non-destructive capabilities. Pioneering numerical studies focused on this
convenient method can be found in [29, 30]. When the electric field carrying
the temporal information of the electron bunch propagates through an EO
crystal, the effective principal axis of the crystal rotates because of the
Pockels effect. With the incidence of a probe laser beam on a crystal, time
information of the electron bunch can be reconstructed by analyzing the phase
retardation. The longitudinal distribution of electrons is encoded
transversely into the probe laser beam by setting a relative angle between the
probe laser and the signal field. This is referred to as the “EO spatial
decoding technique” [26].

EO sampling of transition radiation (TR) of an electron bunch has been
conducted to measure the relative electron longitudinal profile [31, 32] or
the field strength of a terahertz pulse [33]. However, simultaneous
measurement of the absolute current and transverse profiles has not been
attempted. By the incidence of the TR pulse onto the EO crystal and the
performing the spatial decoding, the TR field produced by the electron bunch
is measured both temporally and spatially. We named the detection method as
the “TR-EO” method in this study.

In this article, we performed numerical calculations of EO spatial decoding of
TR produced by a relativistic electron bunch passing through a metal foil. The
calculations involved the imaging of polychromatic TR from electron bunches
with arbitrary shapes and energies and the EO spatial decoding process,
including TR absorption by the crystal, smearing, and phase mismatch, which
are affected by the relative angle and speed difference between the TR and the
probe laser. All calculations in this study are defined in SI units for ease
of implementation. Detailed and simplified models were discussed for the
calculations of both the TR imaging and EO signal generation. Through the
comparison of results, the feasibility of the simplified models for a quicker
analysis was demonstrated. In addition, we calculated the impact of temporal
energy chirps of electrons on the EO signals and concluded that the current
profiles dominated the shapes of the EO signals.

This paper is organized as follows: the experimental setup of this numerical
study is described in Sec. II. The calculation of the polychromatic TR field
via a detailed Huygens–Fresnel diffraction calculation and a simplified
analytic model based on Fraunhofer approximation is described in Sec. III. The
EO spatial decoding process is described in Sec. IV using a simplified model
and a detailed 2D interpolation method. The impact of the temporal energy
chirp of electrons on the EO signal is discussed in Sec. IV. Thereafter, noise
in transition radiation generated by the optical system is calculated in Sec.
V. Finally, the discussion and conclusions of this study are presented in Sec.
VI.

##  II Experimental concept of this numerical study

Figure 1: Conceptual scenario of the TR-EO spatial decoding set-up for an
LWFA. e-: electron bunch; OAP: Off-axis-parabolic gold mirror; S: polarizer
before EO crystal; P: polarizer after EO crystal; $\lambda$/2: half waveplate;
$\lambda$/4: quarter waveplate; and BP: bandpass filter. The inset shows the
coordinates in the imaging plate. The [-1,1,0] axis of the crystal is along
the $Y_{D}$ direction. The polarization of the incident probe laser is also
parallel to the $Y_{D}$ direction.

Numerical calculations were conducted to quantitatively explain the aspects
that can be achieved using the experimental concept illustrated in Fig. 1. In
the LWFA, a drive laser beam was incident on a gas target to generate an
electron beam. Therefore, electron bunches propagated with the outgoing drive
laser beam. However, an intense drive laser can result in damage and noise
during particle parameter measurements in the LWFA. Therefore, a metal foil
(e.g., aluminum, or stainless steel) was placed in the electron beam path to
eliminate noise generated by the drive laser beam and produce TR.
Subsequently, the TR emission was imaged using two off-axis parabolic (OAP)
mirrors and was incident onto an EO crystal. Typically, crystals with a zinc-
blende structure are used. An ultrashort probe laser beam with a small
incident angle $\theta_{p}$ to the EO crystal was used to perform spatial
decoding. The [-1, 1, 0] axis of the crystal was parallel to the polarization
direction of the probe laser and orthogonal to the plane formed by the probe
laser beam and the TR propagation direction. Further, half-wave and quarter-
wave plates were inserted to perform the near-cross-polarization [29]
detection. The TR passing through the crystal was then imaged again to measure
the optical transition radiation (OTR). Because the TR field had a wavelength
range from visible to infrared, it is recommended that the EO crystal be
placed inside a vacuum chamber to avoid unwanted absorption from the vacuum
window and air. For the same reason, the gold-coated OAP was used instead of a
lens to deliver the TR pulse beam.

Figure 2: Calculation geometry of the imaging of TR of an electron bunch. OAP1
and OAP2 have effective focal lengths of $F_{1}$ and $F_{2}$, respectively.
The distance between the OAPs is $L$.

##  III Three-dimensional TR electric field in the coordinates ($X_{D}$,
$Y_{D}$, $T$)

In this section, we present a 3D calculation of the TR field imaged on the EO
crystal. The calculation geometry is shown in Fig. 2. An electron bunch passed
through a radiator (source plane), and TR was produced. The produced TR was
collected by the first OAP (OAP1) with a focal length of $F_{1}$. The source
point was placed at the focus of OAP1. Consequently, the collimated TR light
was focused by OAP2 on the EO crystal. Here, ($X_{s},Y_{s}$) and
($X_{D},Y_{D}$) denote the mesh grids in the source and detector planes (EO
crystals), respectively, and ($X_{l1},Y_{l1}$) and ($X_{l2},Y_{l2}$) are the
mesh grids in the planes of OAP1 and OAP2, respectively, with the OAPs
considered as thin lenses. The mesh grids were arranged symmetrically with the
centers set at (0, 0). The distance between OAP1 and OAP2 was denoted as $L$.
Further, $\bm{r_{b}}=(X_{b},Y_{b},Z_{b})$ denotes the relative position of a
single electron inside the electron bunch. The wave number of the TR was
$k=\sqrt{k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}$ and it had a relative angle of
$\theta$ relative to the z-axis. The direction of the [-1,1,0] axis of the
crystal was along the “$Y$” axes in Fig. 2. For EO spatial decoding, the
polarization of the probe laser beam was set along the [-1,1,0] axis of the EO
crystal to achieve maximum phase retardation.

In this study, the third dimension “$T$” or “$\omega$” was assumed to be
decoupled from the transverse beam profile. First, we calculated the 2D
spatial distribution of the TR field $E_{y}(X_{D},Y_{D},\omega)$ for a single-
frequency component. Subsequently, the 2D image was multiplied by the
longitudinal form factor $F_{z}(\omega)$. Finally, we performed a Fourier
transform of $E_{y}(X_{D},Y_{D},\omega)$ in the $\omega$ dimension to obtain
$E_{y}(X_{D},Y_{D},T)$ distribution.

By treating the metal foil as an ideal conductor, the TR of a single electron
can be calculated as the total reflection of the self-field (superposition of
pseudo-photons). Here, without loss of generality, when passing through the
foil, a single electron at position (0, 0) of the source plane is considered.
The self-field expression (further details presented in Appendix A) of an
electron in the frequency and wavenumber domains is expressed as,

$E^{s}(X_{s},Y_{s})=i\dfrac{e}{(2\pi)^{2}\epsilon_{0}v}\iint{dk_{x}dk_{y}}\exp[i(k_{x}X_{s}+k_{y}Y_{s})]\dfrac{(k_{x},k_{y})}{k_{x}^{2}+k_{y}^{2}+\alpha^{2}}$
(1a)
$E^{s}(X_{s},Y_{s})=-\dfrac{e\alpha}{2\pi\epsilon_{0}v}\dfrac{(X_{s},Y_{s})}{\sqrt{X_{s}^{2}+Y_{s}^{2}}}K_{1}(\alpha\sqrt{X_{s}^{2}+Y_{s}^{2}})$
(1b)

where $\alpha=\omega/\gamma v$, where $v$ is the electron velocity, $K_{1}$ is
a Bessel function of the second type, and $\epsilon_{0}$ denotes the electric
constant. Eq. (1a) and Eq. (1b) are equivalent and used as the source field of
the TR from an electron.

###  III.1 2D spatial imaging of the TR field

Calculations related to OTR imaging have been investigated in accelerator
research to measure the transverse size of electron bunches. Both incoherent
and coherent OTR imaging [34, 35, 36, 37, 38] have been focused upon. In these
studies, a far-field model based on Fraunhofer assumption was widely used
because the TR was specially treated in the optical range. However, because
the EO sampling of the TR covers the wavelength range of zero to a few hundred
THz, the validity of the Fraunhofer model should be discussed.

####  III.1.1 Detailed numerical diffraction calculation based on the
Huygens–Fresnel principle

We calculated the 2D TR field $E_{y}(X_{D},Y_{D})$ based on the
Huygens–Fresnel principle[39] as follows:

$U(x_{2},y_{2})=\frac{z}{i\lambda}\iint_{\Sigma_{1}}U(x_{1},y_{1})\frac{\exp(ikr_{12})}{r^{2}_{12}}dx_{1}dy_{1}$
(2)

where “1” and “2” are the indices of planes “1” and “2”, respectively, with a
distance of $z$, and $U$ is the scalar field distribution. Further, ($x_{1}$,
$y_{1}$) and ($x_{2}$, $y_{2}$) are points on the two planes with the distance
between them being $r_{12}=\sqrt{z^{2}+(x_{2}-x_{1})^{2}+(y_{2}-y_{1})^{2}}$.
The reasons for choosing the Huygens–Fresnel principle instead of Fresnel
approximation are as follows: (i) the para-axial condition might not be
fulfilled if the electrons have low energy because the TR has an angular peak
at $\theta_{peak}\sim 1/\gamma$; and (ii) the acceptance angle of OAP1 is
occasionally large in the experiment.

Figure 3: Calculation results based on the Huygens–Fresnel principle. The
$|E_{y}(X,Y,f)|$ at the positions of OAP1 (a, d, g, and j), OAP2 (b, e, h, and
k), and the EO crystal (c, f, i, and l) are plotted. The first row (a–c) shows
the results of {electron energy, metal foil size, wavelength} = {250 MeV, 10
cm, 100 $\mu$m}.The second row (d–f) shows the results of {250 MeV, 1 cm, 100
$\mu$m}. The third row (g–i) shows the results of {250 MeV, 10 cm, 10 $\mu$m}.
The fourth row (j–l) shows the results of {10 MeV, 10 cm, 10 $\mu$m}. Figure
4: Calculation results of $|E_{y}(X,Y,f)|$ at the EO crystal based on an
analytic model based on Fraunhofer assumption. Here, (a) shows the result of
{electron energy, wavelength} = {250 MeV, 100 $\mu$m}, (b) shows the result of
{250 MeV, 10 $\mu$m}, and (c) shows the result of {10 MeV, 10 $\mu$m}.

Equation (2) can be treated easily as a convolution. We defined the transfer
function as:

$H(x,y)=\frac{z}{i\lambda}\frac{\exp(ik\sqrt{z^{2}+x^{2}+y^{2}})}{z^{2}+x^{2}+y^{2}}$
(3)

The Huygens–Fresnel principle expressed in Eq.(2) can be rewritten as follows:

$U(x_{2},y_{2})=\iint_{\Sigma_{1}}U(x_{1},y_{1})H(x_{2}-x_{1},y_{2}-y_{1})dx_{1}dy_{1}$
(4)

By setting both planes with the same mesh grid, the diffraction can be
calculated using the following Fourier transformation:

$U_{2}=U_{1}\ast H=\mathcal{F}^{-1}[\mathcal{F}(U_{1})\times\mathcal{F}(H)]$
(5)

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote the Fast Fourier
Transfomation (FFT) and Inverse FFT, respectively. Three transfer functions
$\\{H_{1},H_{L},H_{2}\\}$ existed in the system with different longitudinal
distances as $\\{F_{1},L,F_{2}\\}$.

There were three optical apertures with sharp edges (step functions) in the
imaging system: (i) a finite size of the metal foil $D_{s}$, (ii) a finite
size of OAP1; and (iii) a finite size of OAP2. In this study, we set up two
OAPs with the same radius $R_{l}$ and focal length $F$. These apertures can be
described by the following matrices:

$\displaystyle P_{s}(X_{s},Y_{s})$
$\displaystyle=\\{|X_{s}|,|Y_{s}|\\}<D_{s}/2$ (6) $\displaystyle
P_{l}(X_{l},Y_{l})$ $\displaystyle=\sqrt{X_{l}^{2}+Y_{l}^{2}}<R_{l}$

where $P^{s}(X_{s},Y_{s})$ is the square pupil (metal foils often have square
shapes) at the source plane and $P^{l}(X_{l},Y_{l})$ is the round pupil
created by the OAPs.

An OAP with an effective focal length $F$ can be considered an ideal lens
without spherical aberrations (explanations are provided in Appendix B). The
converging effect of a lens results in a relationship between the field at the
left side of the lens $E^{ll}$ and the right side of the lens $E^{lr}$ as
follows:

$E^{lr}(X_{l},Y_{l})=E^{ll}\exp(-ik\dfrac{X_{l}^{2}+Y_{l}^{2}}{2F})$ (7)

where $F$ is the focal length of the lens. The wavefront difference before and
after the OAP is expressed as,

$PT_{l}=\exp(-ik\dfrac{X_{l}^{2}+Y_{l}^{2}}{2F})$ (8)

Considering the transverse distribution of $g_{\perp}(X_{s},Y_{s})$, the
Coulomb field of one slice in the electron bunch at the source plane is a
convolution, which is expressed as follows:

$Cou^{s}(X_{s},Y_{s})=g_{\perp}\ast E^{s}$ (9)

Here, we used the expression of the field of a single electron $E^{s}$, as
expressed in Eq. (1b). The transverse distribution was arbitrary. The TR field
on the left side of the OAP1 was calculated as follows:

$E^{ll1}(X_{l1},Y_{l1})=(P_{s}\times Cou^{s})\ast H_{1}$ (10)

The diffractive propagation of the TR from the left side of OAP1 to the left
side of OAP2, including the phase shift and aperture of OAP1, was calculated
as follows:

$E^{ll2}(X_{l2},Y_{l2})=(P_{l}\times PT_{l}\times E^{ll1})\ast H_{L}$ (11)

Similarly, TR field propagation from the left side of OAP2 to the EO crystal
was calculated as follows:

$E^{D}(X_{D},Y_{D})=(P_{l}\times PT_{l}\times E^{ll2})\ast H_{2}$ (12)

Equations (9–12) show the quantitative calculation procedures for the
diffractive propagation of the TR field. We set $R_{l}$ = 25.4 mm and $F$ =
190 mm for the current and subsequent calculations. The distance between both
the OAPs was set to $L$ = 40 cm. The TR field profiles of $E^{ll1}$,
$E^{ll2}$, and $E^{D}$ were calculated to check the diffraction patterns, as
shown in Fig. 3. We first performed calculations using an electron energy of
250 MeV with a TR wavelength of 100 $\mu$m. In this case, the effective
electron radius[36, 34] was $\gamma\lambda\sim$ 5 cm. A transverse beam of
size (50, 50) $\mu$m was assumed. The results of metal foil sizes of $D_{S}$ =
{1, 10} cm were compared, as shown in Figs. 3(a–f). Although the diffraction
patterns at the positions of the OAPs were slightly different, the TR field
distributions in the EO crystal were almost identical. This calculation
facilitated an understanding of the diffraction of the TR field at
intermediate positions inside the optical system. Furthermore, this indicates
that, in most cases, the variation in the metal foil size should merely result
in a difference in the field strength at the imaging point.

In addition, diffraction calculations were conducted in the short-wavelength
region. The calculations for {250 MeV, 10 $\mu$m} and {10 MeV, 10 $\mu$m} are
shown in Figs. 3(g–i) and 3(j–l), respectively. The diffraction effect was
less noticeable at a wavelength of 10 $\mu$m, and the focus sizes at the EO
crystal were smaller. With an electron energy of 250 MeV, certain diffraction
patterns remained. The diffraction patterns were invisible for an electron
energy of 10 MeV because the effective electron radius was only $\sim$ 200
$\mu$m.

####  III.1.2 Analytic model based on Fraunhofer assumption

Although a detailed calculation based on the Huygens–Fresnel principle can
provide relatively accurate diffraction results for an arbitrary setup, the
process becomes time-consuming when the number of frequency components exceeds
a few hundred and the electron energy spectra are broad. Thus, we examined the
feasibility of using a simplified analytical model based on the Fraunhofer
assumption. The derivation process is described in [34, 35]. To perform a
quantitative calculation of the absolute field strength, we derived the TR
field produced by a single electron at the imaging point in the SI unit as
follows:

$E_{0}^{D}(X_{D},Y_{D})=-\frac{e}{\epsilon_{0}}\frac{1}{\lambda
Mv}\frac{(X_{D},Y_{D})}{\rho_{D}}f(\theta_{m},\gamma,\zeta)$ (13)

where $\rho_{D}=\sqrt{X_{D}^{2}+Y_{D}^{2}}$. $\theta_{m}=R_{l}/F$ is the
acceptance angle of the OAP and $\zeta=k\rho_{D}/M$ is the normalized distance
from the center of the TR field to the calculation point in the detector
plane. Further, $M=F_{2}/F_{1}$ is the magnification of the imaging system. In
this study, $M$ = 1.
$f(\theta_{m},\gamma,\zeta)=\int_{0}^{\theta_{m}}\theta^{2}/(\theta^{2}+\gamma^{-2})J_{1}(\zeta\theta)d\theta$
is a diffraction function defined by the optical system and the electron
energy, where $J$ is a Bessel function of the first kind. This integration can
be performed numerically. When $1/\gamma\ll\theta_{m}$, it is further
simplified as follows [34]:

$f(\theta_{m},\gamma,\zeta)\approx\zeta^{-1}[\gamma^{-1}\zeta
K_{1}(\gamma^{-1}\zeta)-J_{0}(\zeta\theta_{m})]$ (14)

The electric field formed by an electron slice was then calculated as follows:

$E^{D}(X_{D},Y_{D})=g_{\perp}\ast E_{0}^{D}$ (15)

The following assumptions were made to obtain Eqs. (13,14): (i) Fraunhofer
approximation. The quadratic phase term in the source plane
$\exp[ik(X_{s}^{2}+Y_{s}^{2})/2F]$ was neglected, and the metal foil size was
infinite. (ii) The distance between OAP1 and OAP2 is zero. An ideal thin lens
was assumed, and intermediate diffraction was ignored. (An explanation for why
the 2-OAP system can be geometrically treated as a thin-lens imaging system
has been provided in Appendix B). (iii) The diameter of OAP1 acts as the only
pupil in the system. (iv) The electrons possess high energy to satisfy the
condition $1/\gamma\ll\theta_{m}$. It is expected that the field strength
values from the detailed Huygens–Fresnel calculation should be smaller than
those from the simplified model, wherein several diffraction processes are not
considered. Therefore, the errors generated in this approach should be
assessed.

The calculation results for the same electron parameters as those in the
previous section are shown in Fig. 4. Comparisons of Figs. 3(c), 3(i), and
3(l) with Figs. 4(a), 4(b), and 4(c), respectively, revealed that the focus
shapes are almost identical over a large range of electron energies and
wavelengths, though small errors exist. As our prediction, the peak values
from the Huygens–Fresnel calculation were smaller than those of the analytical
model, with errors of (7%, 0.47%, and 5.5 %) in the cases of ({250 MeV, 100
$\mu$m}, {250 MeV, 10 $\mu$m}, and {10 MeV, 10 $\mu$m}), respectively. The
error in the case of {250 MeV, 100 $\mu$m} primarily originated from the
diffraction loss in the source plane and propagation. Whereas, the error for
{10 MeV, 10 $\mu$m} resulted from the assumption $1/\gamma\ll\theta_{m}$. In
the high-energy and short-wavelength regions, the errors were ignored.

###  III.2 Calculation of the OTR image

Before proceeding with the 3D calculations, we explain the calculation method
for OTR images. OTR was used to determine the electron transverse bunch size.
In the optical range, the wavelength was excessively short, which resulted in
the diffraction difference between the Huygens–Fresnel principle and the
Fraunhofer assumption being unnoticeable. Therefore, a simplified analytical
model was used. The OTR image is an integrated TR intensity with a dimension
“$z$”. Therefore, the issue of coherence expressed by the form factor must be
considered.

For a certain frequency, the TR field in the EO crystal is the sum of all
electrons in the bunch, that is, the integration of the TR produced by all
slices, as illustrated in Fig. 2. The TR field of one frequency component in
the EO crystal is expressed as,

$E^{D}=N\iint
dx_{b}dy_{b}g_{\perp}(x_{b},y_{b})E_{0}^{D}(X_{D}-x_{b},Y_{D}-y_{b})\\\
\times\int dz_{b}g_{\parallel}(z_{b})\exp(i\omega z_{b}/v)$ (16)

Figure 5: Calculation results of the OTR intensity distribution (arbitrary
unit) at a wavelength of 800 nm. Initial transverse sizes of the Gaussian
bunch are ($\sigma_{x}$, $\sigma_{y}$) = (50, 50) $\mu$m. The electron energy
is 100 MeV and the charge is 32 pC. (a) and (b) show $I_{icoh}$ and $I_{coh}$
at an electron bunch duration of 1 fs. (c) and (d) show the line-outs along
the $X$ axis with $Y$ = 0 in (a) and (b), respectively.

Here, $E^{D}$ represents any {x, y} component of the TR field and $N$ is the
number of electrons. The integration along the $z$ direction is the
longitudinal form factor $F_{z}=\int dz_{b}g_{\parallel}(z_{b})\exp(i\omega
z_{b}/v)$, where $g_{\parallel}(z_{b})$ is the longitudinal distribution of
the bunch. The term $\exp(i\omega z_{b}/v)$ is given in Eq. (57) in Appendix
A. The intensity of the OTR is
$I_{OTR}=\sqrt{\epsilon_{0}/\mu_{0}}(|E_{x}^{D}|^{2}+|E_{y}^{D}|^{2})/2$,
which can be separated into an incoherent component $I_{icoh}$ and a coherent
component $I_{coh}$ [40] as follows:

$I_{icoh}=\sqrt{\epsilon_{0}/\mu_{0}}/2\times
N\times[g_{\perp}\ast(|E_{0,x}^{D}|^{2}+|E_{0,y}^{D}|^{2})]$ (17)

$I_{coh}=\sqrt{\epsilon_{0}/\mu_{0}}/2\times N(N-1)\times|F_{z}|^{2}\\\
\times(|g_{\perp}\ast E_{0,x}^{D}|^{2}+|g_{\perp}\ast E_{0,y}^{D}|^{2})$ (18)

where $\mu_{0}$ is the magnetic constant. The ratio of $I_{coh}/I_{icoh}$ is
strongly dependent on the electron number and form factor. In the optical
wavelength region, the ratio was sensitive when the electron bunch had a
femtosecond or sub-femtosecond bunch duration. The calculation results for
Gaussian bunches with a duration of 1 fs are shown in Fig. 5. As evident,
$I_{icoh}$ exhibited a shape similar to the original electron bunch profile
(Figs. 5(a) and 5(c)), whereas $I_{icoh}$ exhibited a “donut” shape because of
radial polarization and the symmetry of the electron transverse profile (Figs.
5(b) and 5(d)). In another calculation with a bunch duration of 2 fs,
$I_{coh}$ decreased by a factor of $10^{7}$, whereas the change in $I_{icoh}$
was irrelevant with respect to the bunch duration.

With the number of electrons set as $2\times 10^{8}$ in this calculation, the
peak value of $I_{coh}$ was $2.6\times 10^{3}$ times higher than that of
$I_{icoh}$. This indicates that, in the bunch, if there exists a femtosecond
spike with 0.1% of the total charge (tens to hundreds of fC), the overall OTR
intensity profile exhibits a coherent pattern, as shown in Fig. 5(b). Such
coherent OTR caused by femtosecond or sub-femtosecond spikes have been
reported in both conventional accelerators [37] and LWFA [41].

When the OTR was incoherent, deconvolution was conducted to retrieve the
electron transverse profile by directly using the point spread function
$|E_{0,x}^{D}|^{2}+|E_{0,y}^{D}|^{2}$, as explained in Eq. (17) and Ref. [34].
However, when the OTR exhibits a coherent pattern, straightforward
deconvolution is not possible. Thus, to estimate transverse electron sizes, an
empirical linear relationship $X_{max}=A+B\sigma_{x}$ [38] was discovered by
evaluating the peak positions of the lineout, as shown in Fig. 5(d), where $A$
and $B$ are constants, $X_{max}$ is the peak position in the OTR signal, and
$\sigma_{x}$ is the Gaussian bunch size. Therefore, to a certain extent, this
method can be used to determine the electron bunch size. The reconstruction
method for the fine structure of the original electron transverse profile from
$I_{coh}$ will be investigated elsewhere.

Figure 6: Calculation results of the 3D TR field $E_{y}(X,Y,T)$. (a) shows the
($X$, $Y$) distribution of the maximum value of $|E_{y}|$ in the third
dimension $T$. (b) shows the ($T$,$Y$) distribution of $E_{y}$ at $X=0$. (c)
shows the ($X$,$T$) distribution of $E_{y}$ at $Y=-200\mu$m.

###  III.3 3D TR field at the EO crystal

With a TR field of single frequency denoted by Eq. (16), the 3D TR field can
be obtained using the Fourier transform as follows:

$E(X_{D},Y_{D},t)=\frac{1}{2\pi}\int d\omega
E^{D}(X_{D},Y_{D},\omega)\exp(-i\omega t)$ (19)

Calculations were performed with an electron energy of 100 MeV, a charge of 32
pC, and a bunch duration of 20 fs (rms). The size of the transverse bunch was
(50, 50) $\mu$m. The calculations dealt with the time range of (-0.5–0.5) ps.
The frequency range was zero to a few hundred THz.

To better explain the 3D (X,Y,T) profile, we plotted three 2D contours, as
shown in Fig. 6. The spatial distribution of the maximum temporal TR value
$|E_{y}|_{max}$ is shown in Fig. 6(a). The strong area of the TR field was
concentrated within a few hundred micrometers. The size of the concentrated
area was related to the electron energy and transverse bunch sizes. The
$(T,Y)$ distribution of $E_{y}$ is shown in Fig. 6(b). The opposite values on
both sides of $Y=0$ resulted from the radial polarization of TR. The EO
signals were primarily determined by this distribution. The $(X,T)$
distribution of $E_{y}$ at $y=-200\mu$m is shown in Fig. 6(c). The temporal
profile of the TR field was not uniform along the $X$ axis. Previously, when
experiments of EO spatial decoding were conducted, the transverse field
distributions were not considered. However, when the magnification is not
sufficiently large and the electron bunch has a small transverse size, the
spatial strength profile of the field should be counted when uncompromising
accuracy is required.

Figure 7: Lineouts of field profiles. (a) shows the temporal profiles at two
vertical positions $Y=-200$ $\mu$m and $Y=-800$ $\mu$m in Fig. 6(b). (b) shows
the temporal profiles at two horizontal positions $X=0$ $\mu$m and $X=200$
$\mu$m in Fig. 6(c). (c) shows the frequency domain profile of the TR field at
three points (0, -200) $\mu$m (blue), (0, -800) $\mu$m (red), and (200, -200)
$\mu$m (green).

For a better elaboration of the spatial broadening of the TR field in the
imaging plane, temporal and frequency domain lineouts were plotted, as shown
in Fig. 7. The lineout plots for $Y=-200$ $\mu$m and $Y=-800$ $\mu$m
illustrated in Fig. 6(b) are shown in Fig. 7(a). The lineouts at positions
$x=0$ and $x=200$ $\mu$m are illustrated in Fig. 7(b). As evident, the field
strength weakened and the duration was prolonged with an increase in the
distance from the center of the field, both horizontally and vertically. The
frequency distributions of the three points are shown in Fig. 7(c). The
frequency distribution was narrower and shifted toward the long-wavelength end
at greater distances from the center. The results shown in Fig. 7 indicate
clear differences in the frequency compositions when measuring the field at
different positions. Thus, the spatial distribution of the TR field should be
considered for all forms of measurements, including temporal [23] and
frequency [42] domain analyses.

Figure 8: Comparison of the 3D calculations of the $|E_{y}|_{max}$ using the
two approaches with different electron energies. All the figures share the
same colormap limit. (a) and (c) show the calculations of Huygens–Fresnel
diffractions. (b) and (d) show the calculations using the analytic model based
on Fraunhofer assumption. (a) and (b) show the results with an electron energy
of 20 MeV. (c) and (d) show the results with an electron energy of 200 MeV.

###  III.4 Comparison of 3D TR field calculations made using the two
approaches

In Sections III.A and III.B, two approaches for calculating the image of the
TR field: detailed diffraction based on the Huygens–Fresnel principle and an
analytic model based on the Fraunhofer assumption for calculating, were
elaborated. We observed slight differences in the monochromatic imaging
calculations for various electron energies and wavelengths, as shown in Figs.
3 and 4. To confirm the possibility of adopting a less time-consuming approach
in 3D calculations, the results of both approaches were compared.

The calculations from both approaches with different electron energies are
plotted in Fig. 8. The electron bunch parameters are as follows: charge = 32
pC, transverse beam size = (50, 50) $\mu$m, and bunch duration = 20 fs (rms,
gaussian). The TR field distributions obtained from both approaches were
almost identical. With electron energies of both 20 and 200 MeV, the peak
values of $|E_{y}|_{max}$ from the detailed calculations were only 3 % lower
than those of the simplified analytic model. The minute difference
necessitates the use of an analytical model when performing parameter fitting.
Based on the calculations above and in the previous subsections, we conclude
that the analytical model is appropriate in most cases. In the lateral
sections, without special announcements, an analytic model based on the
Fraunhofer approximation will be applied in the calculation of the 3D TR
field.

##  IV Numerical study of EO spatial decoding

The fundamental inline EO sampling process has been studied numerically in
detail in [30]. Maximum phase retardation is achieved when the polarization of
the probe laser, optical axis [-1,1,0] of the crystal, and polarization of the
external field are in the same direction.

$\Gamma=\frac{2\pi n_{0}^{3}d}{\lambda_{0}}r_{41}E_{TR}$ (20)

where $n_{0}$ is the refractive index of the central wavelength $\lambda_{0}$
of the probe laser, $r_{41}$ is the EO coefficient, and $d$ denotes the
crystal thickness. Notably, Eq. (20) can only be used to approximately
estimate the peak strength of the external field because the frequency-
dependent response of the crystal is not included. Steffen et al. [29]
proposed a near cross-polarization setup for measuring the shape of a TR field
with good linearity when $\Gamma$ was small. As shown in Fig. 1, a $\lambda/4$
plate was used to eliminate the residual birefringence caused by the EO
crystal. The polarizer pair “S” and “P” are cross polarized to each other. The
$\lambda/2$ plate was slightly rotated by an angle $\theta_{2}$, and the
detected EO signal is calculated as,

$I_{det}(\theta_{2},\Gamma)=\dfrac{I_{p}}{2}[1-\cos(\Gamma+4\theta_{2})]+\delta_{ext}I_{p}+B_{0}$
(21)

where $B_{0}$ is the dark image of the camera, $\delta_{ext}$ is the
extinction ratio of the polarizer pair, and $I_{p}$ is the original intensity
profile of the probe. The background without an external TR field is expressed
as follows:

$B_{1}=\dfrac{I_{p}}{2}[1-\cos(4\theta_{2})]+\delta_{ext}I_{p}+B_{0}$ (22)

We conducted data processing of $I_{sig}=(I_{det}-B_{1})/(B_{1}-B_{0})$ to
eliminate the impact of nonuniformity of the original probe transverse profile
as follows:

$I_{sig}=\dfrac{\cos(4\theta_{2})-\cos(\Gamma+4\theta_{2})}{1-\cos(4\theta_{2})+2\delta_{ext}}$
(23)

For EO spatial decoding, the detailed EO signal generation process has rarely
been investigated outside the geometric temporal mapping relationship:

$c\Delta\tau=\tan\theta_{p}\Delta\xi$ (24)

Figure 9: Geometric explanation of the EO spatial decoding process. (a) shows
the normalized probe intensity profiles at propagation depths of 0 (dotted
black), 20 $\mu$m (red), 100 $\mu$m (green), and 200 $\mu$m (blue). (b)
illustrates the geometry of the walk-off between the TR field and the probe
laser beam. The green and red solid lines show a timeline in the TR field and
probe pulse when incident on a slice of the EO crystal, respectively. The
dashed lines show those at the exit of the slice. The probe laser beam has
incident and propagation angles of $\theta_{p}$ and $\theta_{p}^{\prime}$,
respectively. (c) shows how the existing probe duration should be considered
in the calculation.

where $\Delta\tau$ is the temporal gap, $\xi$ is the displacement observed by
the camera, and $\theta_{p}$ denotes the incident angle of the probe laser
beam on the surface of the EO crystal. When the dispersion of the TR field
inside the crystal and the spatial distribution of the TR field are included,
such a temporal mapping relationship lacks sufficient information to calculate
the resulting EO signal shape. Thus, we conducted a detailed numerical study
to apply EO spatial decoding to various situations. We arranged the
elaboration from a simplified to a complicated case. The complexity of the
calculation is dependent on the thickness of the crystal and spot size of the
TR field in the image plane. Two main processes were involved in the
calculation:

(i) Temporal elongation of the probe laser beam: Here, we considered the use
of gallium phosphide (GaP) crystals. The broadening effect can be calculated
by considering the group delay dispersion (GDD) with the knowledge of the
refractive index curve in the optical range
$n(\lambda)=\sqrt{a_{1}+a_{2}\lambda^{2}/(\lambda^{2}-a_{3})}$, where $a_{1}$,
$a_{2}$, and $a_{3}$ are constants[29]. The central wavelength of the probe
laser beam was $\lambda_{0}$ = 800 nm. Assuming a Gaussian probe pulse with a
duration of 10 fs (rms) and a flat phase profile in the frequency domain, the
pulse shapes for various propagation depths are shown in Fig. 9(a). At a depth
of 20 $\mu$m inside the GaP crystal, the shape exhibited minimal change. The
probe pulse elongated upon further propagation inside the crystal. The curves
in Fig. 9(a) were calculated with the integrated energy conserved.

(ii) Walk-off between the probe laser beam and the TR field inside the EO
crystal: The TR field was normally incident on a slice of the EO crystal.
Owing to the relative angle and velocity differences between the TR field and
the probe laser beam, a 2D geometry should be considered, as shown in Fig.
9(b). A timeline in the TR field propagated from the entrance to the exit of
the crystal slice (solid green to dashed green). Within the same period, a
timeline in the probe pulse propagated from the solid line to the dashed red
line. As EO spatial decoding examines the transverse profile of the probe, the
points “A” and “A1” shared the same timing on the camera. However, the cross
point (coding point) shifted from “A” to “B”, which introduced walk-off in
both the probe laser beam (denoted by the region $\overline{BA1}$) and the TR
field (denoted by the region $\overline{BA2}$). $\overline{BA1}$ caused
smearing of the EO signal in the transverse probe profile, and
$\overline{BA2}$ is a factor that must be considered when the strength of the
TR field changes rapidly in the transverse direction.

###  IV.1 EO spatial decoding in a simple 1D case

Figure 10: Calculated EO signals at points ($X$, $Y$) = (0, -200) $\mu$m
(blue) and (0, -800) $\mu$m (red) using the near cross-polarization setup with
a 20-$\mu$m crystal. (a) and (b) show the results with electron bunch
durations of 20 fs (rms). (c) and (d) show the results with electron bunch
durations of 50 fs (rms). (a) and (c) correspond to an electron energy of 100
MeV. (b) and (d) show the results with an electron energy of 20 MeV.

As shown in Fig. 7, the TR fields in the temporal and frequency domains can be
calculated at any spatial point on the imaging plane. The information
regarding the frequency-dependent refractive index
$N(\omega)=n(\omega)+i\kappa(\omega)$ and EO coefficient $r_{41}(\omega)$ can
be found in [30]. To further elaborate the process of EO signal generation, we
separated the EO crystals into slices. Each slice $j$ had a thickness of
$\Delta z$ and a propagation depth of $z_{j}$.

The smearing process is described in [43] in the general case considering that
the TR field has an incident angle of $\theta_{s}$. The smearing factor was
calculated as follows:

$S=\dfrac{\sin\theta_{p}+\sin\theta_{s}}{\dfrac{v_{p}\sin\alpha\cos^{2}\theta_{p}^{\prime}}{v_{s}-v_{p}\cos\alpha}+\sin\theta_{p}^{\prime}\cos\theta_{p}^{\prime}}$
(25)

where $\theta_{p}^{\prime}$ and $\theta_{s}^{\prime}$ are the propagation
angles of the probe and the TR field inside the crystal, respectively,
$\alpha=\theta_{s}^{\prime}+\theta_{p}^{\prime}$, $v_{p}$ denotes the group
velocity of the probe, and $v_{s}=c/n(\omega)$ is the phase velocity of the TR
field. Because $\theta_{s}=\theta_{s}^{\prime}=0$, based on the relationship
$\sin\theta_{p}=n_{0}\sin\theta_{p}^{\prime}$, Eq. (25) is reduced to

$S(\omega)=c(\frac{1}{v_{\parallel}}-\frac{1}{v_{s}})$ (26)

where $v_{\parallel}=v_{p}\cos\theta_{p}^{\prime}$ denotes the probe speed
component along the propagation direction of the TR field. For EO spatial
decoding, the initial transverse profile of the probe laser beam was mapped as
a time array $\tau$. In slice $j$, the time table changed owing to the walk-
off $\overline{BA1}$ as follows:

$t_{j}=\tau+S\times z_{j}/c$ (27)

Without considering the TR’s transverse distribution, we calculated only
$E_{TR}^{0}(\omega)$ at a certain spatial point ($X$, $Y$) at the entrance of
the crystal. At a propagation depth of $z_{j}$ inside the crystal, the TR
field in the frequency domain was calculated as:

$E_{TR,j}(\omega)=\frac{2}{1+N(\omega)}\exp[-\kappa(\omega)\omega
z_{j}/c]E_{TR}^{0}(\omega)$ (28)

as their own time references. Here, $2/(1+N(\omega))$ is the amplitude
transmission coefficient based on the Fresnel’s law and
$\kappa(\omega)\omega/c$ is the absorption rate. The effective TR field was
calculated by calculating the EO coefficient $r_{41}(\omega)$ as follows:

$E_{TR,j}^{eff}(\omega)=r_{41}(\omega)E_{TR}(\omega)$ (29)

The EO effect was calculated in the time domain. The effective TR field
$E_{TR,j}^{eff}(\tau)$ was obtained using Fourier transformation
$E_{TR,j}^{eff}(\tau)=\frac{1}{2\pi}\int\exp(-i\omega
t_{j})E_{TR,j}^{eff}(\omega)d\omega$.

$E_{TR,j}^{eff}(\tau)=\frac{1}{2\pi}\int\exp(-i\omega\tau)\exp[-i\omega
z_{j}(\frac{1}{v_{\parallel}}-\frac{1}{v_{s}})]\frac{2}{1+N(\omega)}\exp[-\kappa(\omega)\omega
z_{j}/c]E_{TR}^{0}(\omega)r_{41}(\omega)d\omega$ (30)

The phase retardation in slice $j$ was then calculated as,

$\Delta\Gamma_{j}(\tau)=\frac{2\pi n_{0}^{3}}{\lambda_{0}}\Delta
zE_{TR,j}^{eff}(\tau)$ (31)

By setting the slice number to infinity, the overall phase retardation is the
sum of the contributions from all slices
$\Gamma(\tau)=\sum_{j}\Delta\Gamma_{j}(\tau)$ and can be transformed into an
integral as follows:

$\Gamma(\tau)=\frac{n_{0}^{3}d}{\lambda_{0}}\int
d\omega\frac{2}{1+N(\omega)}E_{TR}^{0}(\omega)r_{41}(\omega)\times\frac{1}{d}\int_{0}^{d}\exp[i\omega(\frac{N(\omega)}{c}-\frac{1}{v_{\parallel}})z]dz$
(32)

Figure 11: Averaged geometric response function with different number of
slices. The results of $|G|_{integral}$ and $|G|_{sum}$ are denoted by blue
dots and red curves, respectively. (a) and (b) correspond to cases with a
small number of slices. (c) and (d) show the results from a large number of
slices. (a) and (c) show the results from a 30-$\mu$m crystal. (b) and (d)
show the results from a 100-$\mu$m crystal. Figure 12: Comparison of the
calculation results between the “1D integral” and “1D sum” models. The red
curves and black dots denote the calculation results from the “1D integral”
and “1D sum” models, respectively. (a) shows the results with a bunch duration
of 20 fs. (b) shows the results with a bunch duration of 50 fs.

Equation(32) expresses the phase retardation obtained using an infinitesimal
slice of the probe laser beam. However, the probe laser beam had a longer
pulse duration. We considered that the probe laser beam contained the maximum
number of slices possible and calculated the overall phase retardation by
summing all slices. This process is illustrated in Fig. 9(c). Two slices,
denoted as “P1” and “P2”, in the probe pulse had a time gap of $\Delta t_{p}$.
In the probe laser beam with an incident angle of $\theta_{p}$, the same TR
wavefront was encoded to the subsequent probe slice with a displacement of
$\Delta\xi_{p}=c\Delta t_{p}/\tan\theta_{p}$ on the camera chip. When using
the temporal mapping relationship in Eq. (24), signal blurring occurred with a
decoded time gap of $\Delta\tau_{p}=\tan\theta_{p}/c\times c\Delta
t_{p}/\tan\theta_{p}=\Delta t_{p}$.

Similarly, if the probe had temporal slices with an index $j$ and an intensity
profile at time $I(t_{p}^{j})$ in the time region between ($-T_{p},T_{p}$),
the total EO phase retardation is expressed as:
$\Gamma_{total}(\tau)=\sum_{j}\Gamma(\tau+t_{p}^{j})I(t_{p}^{j})/\sum_{j}I(t_{p}^{j})$.
This was rewritten as an integral as follows:

$\Gamma_{total}(\tau)=\int_{-T_{p}}^{T_{p}}dt_{p}f(t_{p})\Gamma(\tau+t_{p})$
(33)

where $f(t_{p})$ denotes the normalized temporal intensity profile of the
probe laser beam. Equation (33) is in the form of a cross-correlation between
the temporal profile of the probe laser beam and phase retardation of one
probe slice, as denoted in Eq. (32). Mathematically, when the kernel profile
is symmetrical around zero, similar to that in the case of a Gaussian-like
profile, Eq.(33) is equivalent to a convolution.

Examples of the calculations are shown in Fig. 10 for electron energies of
{100, 20} MeV and bunch durations of {20, 50} fs (rms). Frequency-domain TR
fields at two points ($X$, $Y$) = (0, -200) $\mu$m and (0, -800) $\mu$m were
used as inputs for EO signal generation. The electron transverse bunch sizes
were (50, 50) $\mu$m, and Eq. (23) was used to calculate $I_{sig}$. The probe
laser beam was assumed to exhibit a Gaussian shape with a duration of 10 fs
(rms), and the incident angle of the probe laser beam was $\theta_{p}$ =
25$\degree$. The oscillations in the $I_{sig}$ curves resulted from frequency-
dependent phase mismatching.

The following observations were obtained from Fig. 10. (i) The visibility of
oscillations is strongly dependent on the electron bunch duration. Fierce
oscillations are observed for electron bunching duration of 20 fs. For
electrons with a bunch duration of $>$ 50 fs, the oscillations are
indistinctive and the shapes of $I_{sig}$ are almost identical to the temporal
shape of the TR field. (ii) The signals are broadened and weaker when the
signal was measured at $Y$ = -800 $\mu$m. With a lower electron energy, the
signal intensity further decreases. (iii) The signal intensity ratio of
$I_{-200}/I_{-800}$ varies for different electron bunch durations and electron
energies. The features above indicate that, to monitor the electron bunch
duration or peak current, the spatial point at which the EO signal is measured
must be identified.

Figure 13: Schematic showing the calculation procedure for 2D EO spatial
decoding. (a) shows the 2D probe intensity distribution at coordinates
($X_{p}$, $Z_{p}$). (b) shows the effective TR field inside the crystal at
coordinates ($X_{TR}$, $Z_{TR}$). (c) shows the interpolated TR field at
coordinates ($X_{TR}$, $Z_{TR}$). (d) shows the pattern of the interpolated
effective TR field by coordinate transformation from ($X_{TR}$, $Z_{TR}$) to
($X_{p}$, $Z_{p}$). Figure 14: Schematic showing coordinate transformation
when performing the 2D EO calculation. “O” represents the origin of the
coordinates. “A” ($X_{p},Z_{p}$) is a point in the probe laser. “B”
($X_{interp},Z_{interp}$) is the point in the effective TR field meeting with
point “A” in the crystal slice $z_{j}$.  Figure 15: EO signal generation in
one slice. (a) shows the 2D phase retardation distribution at coordinates
($X_{p}$,$Z_{p}$). (b) shows the accumulated EO signal in one slice in the
near cross-polarization setup. Figure 16: Overall EO signals generated from a
GaP crystal with a thickness of 30 $\mu$m determined using different
calculation models. The field information is obtained at $Y$ = -200 $\mu$m.
(a) shows the calculation results with an electron bunch duration of 20 fs
(rms). (b) shows the calculation results with an electron bunch duration of 50
fs (rms). The results from the “2D”, “TU” and “1D” models are illustrated by
blue solid curves, red dashed, and black dotted curves, respectively.

####  IV.1.1 Reminder on the adequate number of slices in EO sampling
calculations

To manage the temporal elongation of the probe laser beam in a thick crystal,
the EO crystal was separated into slices, as described above. A similar
procedure was followed in [29, 30]. Here, the number of slices must be focused
upon. We proceeded with this explanation by introducing a parameter that
measures the phase mismatching process: the averaged geometric response
function $|G|$. The $|G|$ value of the two approaches is calculated as

$\displaystyle|G|_{integral}=\frac{1}{d}|\int_{0}^{d}\exp(i\Delta kz)dz|$ (34)

$\displaystyle|G|_{sum}=\frac{1}{d}|\sum_{j}\exp(i\Delta
kz_{j})\times(d/N_{slice})|$ (35)

Equation (34) is the final integral of Eq. (32) for a thin crystal and $\Delta
k(\omega)=\omega(\frac{N(\omega)}{c}-\frac{1}{v_{\parallel}})$. Notably, Eq.
(34) is a modification of the inline version reported in Ref. [30], which was
obtained by changing the group velocity of the probe $v_{p}$ to
$v_{\parallel}$, considering the relative angle between the probe and the TR
field. However, for thick crystals, it is inadequate to consider the entire
crystal as a single slice. The average geometric response function is the sum
of phase mismatches in each slice, as described in Eq. (35), where $N_{slice}$
denotes the number of slices used in the calculations and $d$ denotes overall
crystal thickness.

Figure 11 shows the results of the calculations for GaP crystals with
thicknesses of 30 $\mu$m (Figs. 11(a) and 11(c)) and 100 $\mu$m (Figs. 11(b)
and 11(d)). Cases with small $N_{slice}$ had a slice thickness of 10 $\mu$m,
whereas cases with large $N_{slice}$ had slice thickness of 0.25 $\mu$m. When
the number of slices was insufficient, $|G|_{sum}$ differed significantly from
$|G|_{integral}$, which resulted in fake responses in the high-frequency
region. This comparison indicates that, to achieve an appropriate shape of the
EO signal, the number of slices must be increased.

Using the information described above, we carefully examined the validity of
“1D integral” calculations in a slightly thicker crystal with a thickness of
30 $\mu$m. In this crystal, the probe laser beam was slightly elongated in the
last 10 $\mu$m. Figure 12 shows the results of the comparison between the “1D
integral” and “1D sum” calculations. A total of 120 $N_{slice}$ were utilized
in the “1D sum” calculation and the probe elongation in each slice was
calculated. In the “1D integral” calculation, the entire crystal was
considered as a single slice. An unchanged probe shape was obtained by
calculating the average probe profile of 120 slices. At electron bunch
durations of 20 and 50 fs, no differences were observed between the two
models. This conclusion is essential in the analysis of the experimental EO
signal because the calculation with many slices is two orders of magnitude
slower than the “1D integral” calculation. However, as this study utilized the
numerical methodology, the notation “1D model” implies the “1D sum” model in
lateral context for accuracy.

###  IV.2 EO spatial decoding in a general 2D case

In this study, the transverse TR strength distribution was ignored, as
described. However, when the temporal scale of interest is longer, a larger
transverse region must be calculated, where the transverse TR field
distribution may have an impact. Here, the TR field distribution plotted in
Fig. 6(c) was used as the input.

The EO crystal was divided into thin slices. For example, a GaP crystal with a
thickness of 30 $\mu$m was separated into 120 slices. Each slice had a
thickness of $\Delta z$ = 0.25 $\mu$m. The propagation depth at each slice was
$z_{j}=j\Delta z$. The effective 2D TR fields $E_{TR}^{eff,j}(X_{TR},Z_{TR})$
in the time domain and the 2D probe intensities $I_{p}^{j}(X_{p},Z_{p})$ were
calculated for each slice. Here, $Z$ and $Z_{p}$ are $Z_{TR}=-ct_{TR}$ and
$Z=-ct_{p}$, respectively. further, $(X_{TR},Z_{TR})$, and $(X_{p},Z_{p})$ are
the coordinates that co-propagated with the TR and probe laser beams. Notably,
the calculation of $E_{TR}^{eff,j}(X_{TR},Z_{TR})$ already included the EO
effect and longitudinal phase mismatch.

The calculation procedure for one slice of the EO crystal is illustrated in
Fig. 13. The propagation depth was $z_{j}$ = 25 $\mu$m. The 2D intensity
distribution of the probe laser beam with coordinates ($X_{p}$, $Z_{p}$) is
shown in Fig. 13(a). Because the transverse intensity distribution of the
probe was excluded using Eq. (23), the transverse profile of the probe along
$X_{p}$ was set to be uniform. Further, the probe intensity was normalized in
the $Z_{p}$ ($t_{p}$) direction and calculated, including the GDD. The
effective 2D TR field in the same slice is shown in Fig. 13(b).

A mesh grid for interpolation was calculated for each slice. We wish to
clarify that coordinate transformation is not simply a rotation of ($X_{p}$,
$Z_{p}$). The relative propagation between the TR pulse and the probe laser
beam should be considered for the EO effect occurring in a certain slice
$z_{j}$. The geometry is shown in Fig. 14. By including the transverse shift
and longitudinal phase mismatch, the two coordinates ($X_{p}$, $Z_{p}$) and
($X_{TR}$, $Z_{TR}$) were set to have the same origin, O (0, 0). Point “H”
denotes the point where the EO effect occurred at an earlier timing,
contributed by point “A” in the probe laser beam and point “B” in the
effective TR field. Thus, we obtained the relationship
$\overline{HA}=\overline{HB}$. When $\overline{OC}=X_{p}$ and
$\overline{AC}=Z_{p}$, we showed that $\overline{OH}=X_{p}/\cos\theta_{p}$ and
$\overline{HA}=\overline{HB}=X_{p}\cdot\tan\theta_{p}+Z_{p}$. Furthermore, by
subtracting the relative transverse shift $d\cdot\tan\theta_{p}^{\prime}$, the
interpolation coordinates were expressed using Eq. 36.

$\left\\{\begin{aligned}
X_{interp}&=X_{p}/\cos\theta_{p}-d\cdot\tan\theta_{p}^{\prime}\\\
Z_{interp}&=X_{p}\tan\theta_{p}+Z_{p}\end{aligned}\right.$ (36)

where $\theta_{p}^{\prime}$ denotes the propagation angle of the probe beam
inside the crystal. Subsequently, a 2D interpolation of
$E_{TR}^{eff,j}(X_{interp},Z_{interp})$ was conducted in the region determined
using Eq. (36). The interpolated TR field distribution is illustrated in Fig.
13(c). Next, we transformed $E_{TR}^{eff,j}(X_{interp},Z_{interp})$ into the
probe laser coordinates ($X_{p}$, $Z_{p}$), as shown in Fig. 13(d).

Subsequently, we calculated the phase retardation in each pixel as
$\Gamma_{j}(X_{p},Z_{p})=2\pi n_{0}^{3}\Delta d/\lambda_{0}\times
I_{p}^{j}(X_{p},Z_{p})E_{TR}^{eff,j}(X_{p},Z_{p})$, which is the product of
Fig. 13(a) and 13(d). A 2D distribution of $\Gamma$ is shown in Fig. 15(a).
The phase retardation profile $\Gamma_{j}(X_{p})$ along $X_{p}$ was then
calculated by accumulating the results shown in Fig. 15(a) in the $Z_{p}$
direction. The EO signal contributed by this slice was calculated using Eq.
(23). Earlier timing corresponded to larger $Z_{TR}$ and $X_{p}$ in Fig. 14.
The timetable was arranged as $t=-X_{p}\tan\theta_{p}/c$. The EO signal from
one slice is shown in Fig. 15(b). The peak position of the signal was at a
timing earlier than zero because the group velocity of the probe laser beam
was smaller than the phase velocity of the TR field inside the EO crystal. In
this article, without special clarification, the peak of the EO signal is
shifted to “zero” timing.

The information described above indicates the way EO signal generation is
observed within one slice in a 2D geometry. The overall EO signal was
calculated in three steps: (i) calculation of the accumulated phase
retardation $\Gamma_{j}(X_{p})$ in each slice, (ii) calculation of the overall
phase retardation as $\Gamma(X_{p})=\sum_{j}\Gamma_{j}(X_{p})$, and (iii)
calculation of the EO signal using Eq. (23).

Figure 16 shows the 2D EO spatial decoding results at the vertical position
$Y$ = -200 $\mu$m. We performed the calculations with three different models:
(i) “2D” (blue curve) implies the 2D EO calculation by considering the
transverse TR field strength distribution; (ii) “TU” (red dot) implies the 2D
EO calculation assuming that the TR field is transversely uniform; and (iii)
“1D” (black dashed line) denotes the analytical model described in Sec. IV. A.
$\Gamma_{j}$ was calculated for each slice using Eqs. (32) and (33).

The calculations revealed interesting results. (i) The results from “TU” and
“1D” were identical, indicating that our approach for the 2D EO calculation is
adequate. (ii) For a slightly longer bunch duration of 50 fs, the difference
between the “2D” and “TU” was unnoticeable. (iii) For a very short bunch
duration of 20 fs, determining a difference around the main peak of the signal
was challenging. However, the amplitudes of the oscillations at later times
were slightly smaller. Whereas, for a very long electron bunch, although no
calculations were performed, we expected visible differences between the “2D”
and “1D” models. In the case of LWFA, as the electron bunches have duration of
a few fs to tens of fs, with such small errors, the application of the “1D”
model should provide sufficient results when analyzing the EO signal around
the main peak.

###  IV.3 EO spatial decoding signal of temporally chirped electron bunches

Using the calculation method established above, we investigated electron
bunches with multiple energy components. Typically, electron bunches with
energy spreads exhibit chirps in the longitudinal phase space (LPS). We
addressed this issue by separating the electrons into energy slices
corresponding to different timings. For simplicity, a linear energy chirp with
an energy range of 20–300 MeV was assumed. For electron slice j, the timing
was $\delta t_{j}=chirp\cdot(E_{j}-160)$ and the charge weight was $\eta_{j}$.
The longitudinal form factor was calculated using the Fourier transform
$F_{z}^{j}(\omega)=\int\exp(i\omega t)\times\exp[-(t-\delta
t_{j})^{2}/2\sigma_{t0}^{2}]dt$. The result of this integral was obtained as

$F_{z}^{j}(\omega)=\exp(-\frac{\omega^{2}\sigma_{z0}^{2}}{2v^{2}}+\frac{i\omega\delta
z_{j}}{v})$ (37)

for relativistic electrons $v\approx c$. Each slice was assumed to have the
same slicing bunch duration of $\sigma_{t0}$ = 10 fs, corresponding to a
length of $\sigma_{z0}=c\sigma_{t0}$. With the relative center timing of each
slice set to $\delta t_{j}$, we defined $\delta z_{j}=c\delta t_{j}$. We
assessed the signals from “2D” EO calculations at a vertical distance of
$Y_{D}=-200$ $\mu$m. The electron energy spectrum of a two-temperature (“2T”)
distribution $dN/dE\propto 1/T_{1}\exp(-E/T_{1})+1/T_{2}\exp(-E/T_{2})$ was
evaluated, with $T_{1}$ = 20 MeV and $T_{2}$ = 300 MeV. Examples of temporally
chirped electron bunches are shown in Fig. 17. Figure 17(a) shows the energy
distributions (blue) of the “2T” model and a temporal chirp (red) of 0.2
fs/MeV. The 2D plot of the LPS and the current of a two-temperature electron
bunch with a chirp of 0.2 fs /MeV is illustrated in Fig. 17(b). Using the
field expressions presented in Eq. (13), the frequency-domain field component
is the sum of contributions from each electron energy component and the
corresponding $F_{z}^{j}(\omega)$, and is calculated as follows:

$E^{D}(X_{D},Y_{D},\omega)=\sum_{j}\eta_{j}E_{0}^{D}(X_{D},Y_{D},\gamma_{j},\omega)F_{z}^{j}(\omega)$
(38)

Figure 17: Longitudinal energy chirp and current profiles. (a) The blue line
illustrates the two-temperature energy distribution. The red line shows the
slice timings of a positive chirp of 0.2 fs/MeV. (b) shows the 2D plot of the
LPS of the two-temperature linear-chirped electron bunch. The red curve shows
the current profile. Figure 18: Impact of temporal energy chirp on the EO
signals. The blue and red curves illustrate the results of chirped and
nonchirped cases, respectively. (a) shows the field strength profiles at
($X_{D}$, $Y_{D}$) = (0, -200) $\mu$m. (b) and (c) show the EO signals from
current profiles with a chirp of 0.1 and 0.2 fs/MeV, respectively.

To verify the impact of LPS on the EO signal, we compared the chirped and non-
chirped electron bunches with a same current profile and overall charge. The
normalized longitudinal profile was retrieved from the current profile, as
shown in Fig. 17(b). Subsequently, we calculated the longitudinal form factor
$F_{z}(\omega)$ by performing Fourier transformation. In the non-chirped case,
$F_{z}(\omega)$ was identical for all electron energy components. The field
strengths $E_{y}$ are plotted in Fig. 18(a). The calculated EO signals from
the 30-$\mu$m GaP crystal are shown in Figs. 18 (b) and 18(c). Chirps of 0.1
and 0.2 fs/MeV corresponded to temporal spans of 28 and 56 fs, respectively.
The field from the chirped electron bunch exhibited a gentle rising edge
because the high-energy electrons resided at subsequent timings. However, this
small difference did not affect the shape of the EO signals. As shown in Figs.
18(b–c), despite the slight differences in the amplitudes, the shapes of the
EO signals were almost identical.

These calculations demonstrate that the current profile plays a more important
role than the temporal energy chirp in the shapes of the EO signals. When the
LPS is unclear, a non-chirped electron bunch is assumed to be sufficient for
calculating the longitudinal current distribution.

##  V Noise introduced by the TR when the electrons hit OAP1

Figure 19: (a) Calculation geometry of the TR from OAP1. TR generated by OAP1
has a plane wave front after OAP1 and is transported to OAP2 with diffraction.
(b) Beam size of the scattered electrons at the position of OAP1. Figure 20:
TR field strengths at the EO crystal. (a) TR generated from OAP1 with electron
energy of 20 MeV. (b) TR generated from metal foil with electron energy of 20
MeV. (c) TR generated from OAP1 with electron energy of 200 MeV. (d) TR
generated from metal foil with electron energy of 200 MeV.

In the experimental setup illustrated in Fig. 1, electrons passing through the
metal foil hit the gold surface of OAP1, causing secondary TR noise. Thus, to
remove this noise, the insertion of a bending magnet between the metal foil
and OAP1 is recommended. Occasionally, to achieve better spatial resolution,
the focal length of OAP1 is set as short, which renders the insertion of a
magnet difficult. Here, we evaluated the strength of the TR field produced by
OAP1. Figure 19(a) illustrates the calculation process. The scattering caused
by the metal foil was calculated using the following equation that is
dependent on the composition of the metal, thickness, and electron energy
[44]:

$\theta_{0}=\frac{13.6MeV}{\beta cp}z\sqrt{x/X_{0}}[1+0.038\ln(x/X_{0})]$ (39)

where $\theta_{0}$ is the $rms$ divergence of the scattered beam, $\beta c$
and $p$ are the velocity and momentum of the electrons, respectively, $X_{0}$
denotes the radiation length of the material, and $x$ denotes the thickness of
the material. With a propagation distance of $F$ to OAP1, the electron beam
size was $r_{1}=\sqrt{r_{0}^{2}+(\theta_{0}F)^{2}}$, where $r_{0}$ is the beam
size at the position of the metal foil. For alloys such as stainless steel,
the overall radiation length can be estimated as
$1/X_{0}=\sum_{j}w_{j}/X_{j}$, where the $j$th element has a weight of $w_{j}$
and a radiation length of $X_{j}$. The radiation length of each element can be
found in the NIST database [45]. The weights of the elements of a material are
also obtained from various databases. For example, we selected SUS304
stainless steel as the material for the metal foil [46] with a thickness of
100 $\mu$m, an initial beam size of $r_{0}$ = 50 $\mu$m, and $F$ = 190 mm. The
$r_{1}$ for various electron energies at OAP1 are plotted in Fig. 19(b). The
size of the electron beams increased to several millimeters when the electrons
had lower energies. For electrons with energies of 20 and 200 MeV, the beam
sizes at OAP1 were 7623 $\mu$m and 780 $\mu$m, respectively, resulting in
decreased charge densities.

The electron bunch on the metal foil was considered as the point source. The
scattered electrons traversed OAP1 with phases forming a spherical wavefront.
Thus, the Coulomb field of the electrons also exhibited a spherical phase. The
phase difference on a plane perpendicular to the direction of propagation was
described as $\Delta\phi=k(\sqrt{F^{2}+\rho^{2}}-F)\approx k\rho^{2}/2F$.
After reflection from OAP1, the spherical phase was corrected using the
transform function in Eq. (8). The TR emitted from OAP1 was calculated as
follows:

$E^{lr1}_{OAP1}=(g_{\perp}^{OAP1}\ast E^{S})\times P_{l}$ (40)

where $g_{\perp}^{OAP1}$ is the transverse distribution of the scattered
electron beam at OAP1 and $P_{l}$ is the pupil defined by the OAP. $E^{S}$ was
calculated using Eq. (1). The procedure described in Section III.A was
followed for calculating lateral Huygens–Fresnel diffractions in the optical
system. For a fair comparison, we calculated the TR fields at the position of
the EO crystal originating from two sources: the metal foil and OAP1. The TR
generated from OAP1 with electron energies of 20 and 200 MeV are shown in
Figs. 20(a) and 20(c), respectively. The corresponding results of the TR
generated from the metal foil are shown in Figs. 20(b) and 20(d),
respectively. With an electron energy of 20 MeV, the TR field generated from
OAP1 was three orders of magnitude smaller than that generated from the metal
foil. Whereas, at a higher electron energy of 200 MeV, the TR field generated
from OAP1 was six times smaller than that generated from the metal foil.

Apparently, a smaller energy corresponds to a larger scattering angle. This
results in a smaller charge density at the position of OAP1. If most of the
electron energy is concentrated at approximately 100 MeV or a few tens of MeV,
the noise TR source from OAP1 can be ignored. In the experiment, by measuring
the electron energy spectrum, the impact of such noise was examined by
performing this calculation, including the contributions of all energy
components. Notably, temporal elongation of the electron bunches owing to
energy spread was not considered. Thus, the noise should be even weaker when
calculating the bunch elongation.

##  VI Discussion and Summary

In this study, the shapes of the signals in the EO spatial decoding setup were
calculated at arbitrary positions in the image plane. We found that the
duration and oscillation behavior of the EO signals varied with parameters
such as the electron energy, bunch duration, and transverse beam size. Even
with the same electron bunch parameters, the EO signal was broadened further
from the center of the TR field. This behavior is similar to cases wherein
Coulomb fields were measured [26]. Thus, we did not intend to provide a quick
one-to-one correspondence between the duration of the main peak of the EO
signal and the original electron bunch duration. Using the methods described
herein, such correlations can be identified with specific electron bunch
parameters. Spatially resolved detection is strongly recommended when
conducting temporal or spectral measurements of a TR field using an imaging
system. Furthermore, it is necessary to determine the position of the spatial
point at which the experimental signals are measured.

In summary, we conducted a systematic numerical study of EO spatial decoding
of TR from relativistic electron bunches. 3D TR field was determined using
both a detailed calculation based on the Huygens–Fresnel principle and a
simplified analytical model based on Fraunhofer assumption. The results
suggest that the simplified analytical model is sufficient to perform such
polychromatic calculations with considerable accuracy. For EO spatial
decoding, we discussed the process of EO signal generation using both 1D and
2D models. The 1D model was concluded to be sufficient for quick data
analysis. In addition, we demonstrated the minimal impact of the temporal
energy chirp on the shape of the EO signal. Further, we proposed a method to
estimate the noise level for unwanted TRs created by reflection optics.
Therefore, this study will be useful for investigating the 3D charge-density
profiles of ultrafast electron bunches in both LWFA and conventional
accelerators.

###### Acknowledgements.

We are grateful for the encouragement provided by Dr. Y. Sano, and Dr. N.
Kumagai. We thank Dr. J. K. Koga for providing valuable insights into the
current research topic. This work was funded by the JST-Mirai Program, Japan
(Grant No. JPMJMI17A1); the Grant-in-Aid for Early-Career Scientists (Grant
Nos. JP21K17998 and JP23K17152); Grant-in-Aid for Scientific Research (C) from
JSPS KAKENHI, Japan (Grant No. JP22K12665); and the QST President’s Strategic
Grant (Exploratory Research), Japan.

##  Appendix A Formulation of the electric field of transition radiation in SI
units

Herein, the self-field expression for a single electron is derived. Similar
derivations can be found in several textbooks and literature on Gaussian units
[47, 48, 49]. We rewrote the formalization in SI units. Maxwell’s equations in
vacuum are as follows:

$\displaystyle\nabla\cdot\bm{E}$ $\displaystyle=\dfrac{\rho}{\epsilon_{0}}$
(41) $\displaystyle\nabla\cdot\bm{B}$ $\displaystyle=0$
$\displaystyle\nabla\times\bm{E}$
$\displaystyle=-\dfrac{\partial\bm{B}}{\partial t}$
$\displaystyle\nabla\times\bm{B}$
$\displaystyle=\mu_{0}(\bm{J}+\epsilon_{0}\dfrac{\partial\bm{E}}{\partial t})$

The scalar potential $\varphi$ and vector potential $\bm{A}$ are expressed as,

$\displaystyle\bm{B}$ $\displaystyle=\nabla\times\bm{A}$ (42)
$\displaystyle\bm{E}$
$\displaystyle=-\nabla\varphi-\dfrac{\partial\bm{A}}{\partial t}$

The Lorentz Gauge equation is expressed as,

$\nabla\cdot\bm{A}+\dfrac{1}{c^{2}}\dfrac{\partial\varphi}{\partial t}=0$ (43)

Via Eq. (41) and Eq.(42) and Eq. (43), with $c=1/\sqrt{\epsilon_{0}\mu_{0}}$,
the wave equations can be derived as:

$\displaystyle(\nabla^{2}-\dfrac{1}{c^{2}}\dfrac{\partial^{2}}{\partial
t^{2}})\varphi(\bm{r},t)$
$\displaystyle=-\dfrac{\rho(\bm{r},t)}{\epsilon_{0}}$ (44)
$\displaystyle(\nabla^{2}-\dfrac{1}{c^{2}}\dfrac{\partial^{2}}{\partial
t^{2}})\bm{A}(\bm{r},t)$ $\displaystyle=-\mu_{0}\bm{J}(\bm{r},t)$

Here, the variables ($\varphi$, $\bm{A}$, $\rho$, and $\bm{J}$) are expressed
in space and time coordinates ($\bm{r}$ and $t$), where the instant charge,
current, and displacement of a single electron are expressed as follows:

$\displaystyle\rho(\bm{r},t)$ $\displaystyle=-e\delta(\bm{r}-\bm{r}(t))$ (45)
$\displaystyle\bm{J}(\bm{r},t)$
$\displaystyle=-e\bm{v}\delta(\bm{r}-\bm{r}(t))$ $\displaystyle\bm{r}(t)$
$\displaystyle=\bm{v}t$

To avoid discrepancies caused by different conventions of Fourier
transformation in the angular frequency, we used the standard Fourier
transformation as follows:

$\displaystyle g_{1}(\bm{r},t)$
$\displaystyle=\int{df}\int{d^{3}\bm{\tilde{\nu}}}g_{2}(\bm{\tilde{\nu}},f)\exp[i2\pi(\bm{\tilde{\nu}}\cdot\bm{r}-ft)]$
(46) $\displaystyle g_{2}(\bm{\tilde{\nu}},f)$
$\displaystyle=\int{dt}\int{d^{3}\bm{r}}g_{1}(\bm{r},t)\exp[-i2\pi(\bm{\tilde{\nu}}\cdot\bm{r}-ft)]$

where ($\bm{\tilde{\nu}}$, $f$) and ($\bm{r}$, $t$) are the reciprocal
$\bm{\tilde{\nu}}=\bm{k}/2\pi$ and $f=\omega/2\pi$. By performing Fourier
transformations on both sides of Eq. (44), using the relationships in Eq.
(46), the wave equations in the frequency domain are derived as,

$\displaystyle(\bm{\tilde{\nu}}^{2}-\dfrac{f^{2}}{c^{2}})\varphi(\bm{\tilde{\nu}},f)$
$\displaystyle=\dfrac{\rho(\bm{\tilde{\nu}},f)}{4\pi^{2}\epsilon_{0}}$ (47)
$\displaystyle(\bm{\tilde{\nu}}^{2}-\dfrac{f^{2}}{c^{2}})\bm{A}(\bm{\tilde{\nu}},f)$
$\displaystyle=\dfrac{\bm{J}(\bm{\tilde{\nu}},f)}{4\pi^{2}\epsilon_{0}c^{2}}$

where $\rho(\bm{\tilde{\nu}},f)$ and $\bm{J}(\bm{\tilde{\nu}},f)$ are derived
using Eqs. (45) and (46) as

$\displaystyle\rho(\bm{\tilde{\nu}},f)$
$\displaystyle=-e\delta(f-\bm{\tilde{\nu}}\cdot\bm{v})$ (48)
$\displaystyle\bm{J}(\bm{\tilde{\nu}},f)$
$\displaystyle=-e\bm{v}\delta(f-\bm{\tilde{\nu}}\cdot\bm{v})$

The scalar and vector potentials of the field in the $(\bm{k},\omega)$ domain
are obtained as follows:

$\displaystyle\varphi(\bm{\tilde{\nu}},f)$
$\displaystyle=-\dfrac{e}{(2\pi)^{2}\epsilon_{0}}\dfrac{\delta(f-\bm{\tilde{\nu}}\cdot\bm{v})}{\bm{\tilde{\nu}}^{2}-f^{2}/c^{2}}$
(49) $\displaystyle\bm{A}(\bm{\tilde{\nu}},f)$
$\displaystyle=-\dfrac{e\bm{v}}{(2\pi)^{2}\epsilon_{0}c^{2}}\dfrac{\delta(f-\bm{\tilde{\nu}}\cdot\bm{v})}{\bm{\tilde{\nu}}^{2}-f^{2}/c^{2}}$

Using the relationships in Eq.(42), we obtain

$\bm{E}(\bm{\tilde{\nu}},f)=i2\pi
f\bm{A}(\bm{\tilde{\nu}},f)-i2\pi\bm{\tilde{\nu}}\varphi(\bm{\tilde{\nu}},f)$
(50)

The formula for the electric field $E(\bm{k},\omega)$ is derived as follows:

$\bm{E}(\bm{\tilde{\nu}},f)=-i\dfrac{e}{2\pi\epsilon_{0}}\dfrac{f\bm{v}/c^{2}-\bm{\tilde{\nu}}}{\bm{\tilde{\nu}}^{2}-f^{2}/c^{2}}\delta(f-\bm{\tilde{\nu}}\cdot\bm{v})$
(51)

For simplicity, we considered electron propagation along the z-axis. We have
$\bm{v}=(0,0,v)$ and $f=\bm{\tilde{\nu}}\cdot\bm{v}=\tilde{\nu}_{z}v$, where
$\bm{\tilde{\nu}}=(\tilde{\nu}_{x},\tilde{\nu}_{y},\tilde{\nu}_{z})$. With the
property of the delta function $\delta(ax)=\delta(x)/|a|$, Eq.(51) is
rewritten as follows:

$\bm{E}(\bm{\tilde{\nu}},f)=i\dfrac{e}{2\pi\epsilon_{0}v}\delta(f/v-\tilde{\nu}_{z})\dfrac{(\tilde{\nu}_{x},\tilde{\nu}_{y},\tilde{\nu}_{z}/\gamma^{2})}{\tilde{\nu}_{x}^{2}+\tilde{\nu}_{y}^{2}+\tilde{\nu}_{z}^{2}/\gamma^{2}}$
(52)

Thus, the frequency-domain electric field in space can be obtained using the
integral in Eq. (53), with $\bm{r}=(x,y,z)$. As the longitudinal component of
the field was weaker than the transverse component by a factor of
$1/\gamma^{2}$, we focused on analyzing the transverse component of the
electric field in this study. The $(E_{x},E_{y})$ components were calculated
using Eq. ( 54).

$\displaystyle\bm{E}(\bm{r},f)=i\dfrac{e}{2\pi\epsilon_{0}v}\int{d^{3}\bm{\tilde{\nu}}}\exp[i2\pi(\tilde{\nu}_{x}x+\tilde{\nu}_{y}y+\tilde{\nu}_{z}z)]\delta(f/v-\tilde{\nu}_{z})\dfrac{(\tilde{\nu}_{x},\tilde{\nu}_{y},\tilde{\nu}_{z}/\gamma^{2})}{\tilde{\nu}_{x}^{2}+\tilde{\nu}_{y}^{2}+\tilde{\nu}_{z}^{2}/\gamma^{2}}$
(53)

$\displaystyle E_{x,y}(\bm{r},f)=i\dfrac{e}{2\pi\epsilon_{0}v}\exp(i2\pi
fz/v)\iint{d\tilde{\nu}_{x}d\tilde{\nu}_{y}}\exp[i2\pi(\tilde{\nu}_{x}x+\tilde{\nu}_{y}y)]\dfrac{(\tilde{\nu}_{x},\tilde{\nu}_{y})}{\tilde{\nu}_{x}^{2}+\tilde{\nu}_{y}^{2}+f^{2}/v^{2}\gamma^{2}}$
(54)

$\displaystyle E_{x,y}(\bm{r},f)=-\frac{ef/\gamma v}{\epsilon_{0}v}\exp(i2\pi
fz/v)\frac{(x,y)}{\sqrt{x^{2}+y^{2}}}K_{1}(2\pi f/\gamma v\sqrt{x^{2}+y^{2}})$
(55)

$\displaystyle
E_{x,y}(\bm{r},\omega)=i\dfrac{e}{(2\pi)^{2}\epsilon_{0}v}\exp(i\omega
z/v)\iint{dk_{x}dk_{y}}\exp[i(k_{x}x+k_{y}y)]\dfrac{(k_{x},k_{y})}{k_{x}^{2}+k_{y}^{2}+\alpha^{2}}$
(56) $\displaystyle
E_{x,y}(\bm{r},\omega)=-\frac{e\alpha}{2\pi\epsilon_{0}v}\exp(i\omega
z/v)\frac{(x,y)}{\sqrt{x^{2}+y^{2}}}K_{1}(\alpha\sqrt{x^{2}+y^{2}})$ (57)

By calculating the integral in Eq. (54), the transverse field strength is
expressed as Eq. (55), where $K_{1}$ is a Bessel function of the second kind.
By introducing a parameter $\alpha=\omega/\gamma v=2\pi f/\gamma v$, Eqs. (54)
and (55) were changed to the form described by the angular frequency in Eqs.
(56) and (57). Notably, the factor $\dfrac{e\alpha}{2\pi\epsilon_{0}v}$ in
Eqs. (57) differs from $\dfrac{e\alpha}{(2\pi)^{3/2}\epsilon_{0}v}$ in [36,
50], where the unitary convention of Fourier transformation
$g(t)=\frac{1}{(2\pi)^{1/2}}\int g(\omega)\exp(-i\omega t)dt$ was applied.

The derivation of Eq. (55) required the integrals [51] listed as follows:

$\int_{-\infty}^{+\infty}\dfrac{x\sin{ax}}{b^{2}+x^{2}}dx=\pi e^{-ab}$ (58a)
$\int_{0}^{+\infty}e^{-a\sqrt{c^{2}+x^{2}}}\cos(bx)dx=\dfrac{ac}{\sqrt{a^{2}+b^{2}}}K_{1}(c\sqrt{a^{2}+b^{2}})$
(58b)

Figure 21: (a) Explanation of the phase shift introduced by a parabolic
mirror. (b) The overall transmission process from the source to the detector
in a 2-OAP imaging system.

##  Appendix B Imaging system with two OAPs

Herein, a simple explanation of the imaging system composed of two OAPs is
presented. We observed that, geometrically, such a system can be considered as
a thin lens with additional phase delay. The phase-shift effect of the
parabolic mirror is shown in Fig. 21(a). The parabolic surface, plotted as a
red curve, is defined by the function $y^{2}=4f_{p}x$, where $f_{p}$ is the
parent focal length of the parabola. The light indicated by $l_{1}$ was
incident at the center of a 90$\degree$ parabolic mirror. The sideline $l_{2}$
has a transverse distance $\rho$ to $\l_{1}$. The tangent of the parabola at
point $(f_{p},2f_{p})$ has a slope of 1 and a relative angle of 45$\degree$ on
the x-axis. Without a parabolic mirror, the light was reflected by a flat
mirror, as illustrated by the blue line, and the overall optical paths were
the same. The red curve indicates that $l_{2}$ is incident at point
$(x_{1},2f_{p}-\rho)$, where $x_{1}=(2f_{p}-\rho)^{2}/4f_{p}$. Thus, the
introduced path difference is $f_{p}-x_{1}-\rho$, which is calculated as
follows:

$\Delta L(\rho)=-\dfrac{\rho^{2}}{4f_{p}}$ (59)

For the 90$\degree$ OAP, the effective focal length was $f_{e}=2f_{p}$. For
the 2-OAP system shown in Fig. 21(b), the overall transformation function is
calculated as $T(\rho)=\exp[ik(F_{1}+L+F_{2})]\exp[ik(\Delta
L_{1}(\rho)+\Delta L_{2}(\rho))]$. Neglecting the constant phase delay, the
transformation can be written as

$\displaystyle T(\rho)$
$\displaystyle=\exp[-\dfrac{ik\rho^{2}}{2}(\dfrac{1}{F_{1}}+\dfrac{1}{F_{2}})]$
(60) $\displaystyle=\exp[-\dfrac{ik\rho^{2}}{2}\dfrac{1}{f}]$

where $1/f=1/F_{1}+1/F_{2}$ and $f$ is the effective focal length of the
overall system. By placing the source at the focal point of the first OAP,
this imaging system had an effective focal length of $f$ and a magnification
of $F_{2}/F_{1}$. This calculation was performed based on geometrical optics.

##  Appendix C Derivation of the formula for TR imaging using Fraunhofer
approximation

The derivation of an analytical formula in the frequency domain for OTR
imaging has been reported in existing literature [34, 35]. Here, we listed the
details used to confirm the factors in SI units. The imaging system comprised
only one thin lens. The propagation of TR was separated into two sessions: (i)
source to lens and (ii) lens to detector. Assuming a metal foil of infinite
size and neglecting the spherical phase in both the source and detector
planes, the overall diffractive propagation can be expressed as Eq. (61) using
the Fraunhofer approximation.

$E^{D}(X_{D},Y_{D})=\frac{\exp(ikF_{1})}{i\lambda
F_{1}}\frac{\exp(ikF_{2})}{i\lambda
F_{2}}\int_{X_{l,min}}^{X_{l,max}}dX_{l}\int_{Y_{l,min}}^{Y_{l,max}}dY_{l}\exp[-\frac{ik}{F_{2}}(X_{D}X_{l}+Y_{D}Y_{l})]\\\
\times\exp[-\frac{ik}{2F_{2}}(X_{l}^{2}+Y_{l}^{2})]\exp[-\frac{ik}{2F_{1}}(X_{l}^{2}+Y_{l}^{2})]\exp[\frac{ik}{2F_{2}}(X_{l}^{2}+Y_{l}^{2})]\exp[\frac{ik}{2F_{1}}(X_{l}^{2}+Y_{l}^{2})]\\\
\times\int_{-\infty}^{+\infty}dX_{s}\int_{-\infty}^{+\infty}dY_{s}E^{s}(X_{s},Y_{s})\exp[-\frac{ik}{F_{1}}(X_{s}X_{l}+Y_{s}Y_{l})]$
(61)

$E^{D}(\rho_{D},\phi_{D})=\frac{e\alpha}{2\pi\epsilon_{0}v}\frac{1}{\lambda^{2}M}\int_{0}^{\theta_{m}}\theta
d\theta\int_{0}^{2\pi}d\phi_{l}\exp[-\frac{ik}{M}\rho_{D}\theta\cos(\phi_{D}-\phi_{l})]\\\
\times\int_{0}^{+\infty}K_{1}(\alpha\rho_{s})\rho_{s}d\rho_{s}\int_{0}^{2\pi}d\phi_{s}(\cos\phi_{s},\sin\phi_{s})\exp[-ik\rho_{s}\theta\cos(\phi_{s}-\phi_{l})]$
(62)

This system exhibits cylindrical symmetry because the pupil defined by the OAP
has a circular shape. We change the notations of the coordinates as:
$(X_{s},Y_{s})=\rho_{s}(\cos\phi_{s},\sin\phi_{s})$,
$(X_{l},Y_{l})=\rho_{l}(\cos\phi_{l},\sin\phi_{l})$,
$(X_{D},Y_{D})=\rho_{D}(\cos\phi_{D},\sin\phi_{D})$. The Eq. (57) is rewritten
as $E^{s}(\rho_{s})=-\frac{e\alpha}{2\pi\epsilon_{0}v}\exp(i\omega
z/v)K_{1}(\alpha\rho_{s})$. With $\theta=\rho_{l}/F_{1}$ and $M=F_{2}/F_{1}$,
Eq.(61) can be simplified as Eq. (62) by omitting all longitudinal phase
factors. $\theta_{m}=R_{l}/F_{1}$ denotes the acceptance angle of the first
OAP. Using the integral
$\int_{0}^{2\pi}d\beta_{1}(\cos\beta_{1},\sin\beta_{1})\exp[-iA\cos(\beta_{1}-\beta_{2})]=-2\pi
i(\cos\beta_{2},\sin\beta_{2})J_{1}(A)$ and the definition of the normalized
transverse distance in the image plane $\zeta=k\rho_{D}/M$, Eq. (62) can be
simplified as

$E^{D}(\rho_{D},\phi_{D})=\frac{-ek\alpha}{\epsilon_{0}v\lambda
M}(\cos\phi_{D},\sin\phi_{D})\int_{0}^{\theta_{m}}\theta J_{1}(\zeta\theta)\\\
\times\int_{0}^{+\infty}\rho_{s}d\rho_{s}J_{1}(k\theta\rho_{s})K_{1}(\alpha\rho_{s})$
(63)

By using the integral
$\int_{0}^{+\infty}xJ_{1}(ax)K_{1}(bx)dx=a/(a^{2}b+b^{3})$ [51], we have

$E^{D}(\rho_{D},\phi_{D})=-\frac{e}{\epsilon_{0}}\frac{(\cos\phi_{D},\sin\phi_{D})}{\lambda
Mv}\\\
\times\int_{0}^{\theta_{m}}\frac{\theta^{2}}{\theta^{2}+(\gamma\beta)^{-2}}J_{1}(\zeta\theta)d\theta$
(64)

In most cases, $\beta\simeq 1$. Equation(64) can be expressed in Cartesian
coordinates as follows:

$E^{D}(X_{D},Y_{D})=-\frac{e}{\epsilon_{0}}\frac{1}{\lambda
Mv}\frac{(X_{D},Y_{D})}{\sqrt{X_{D}^{2}+Y_{D}^{2}}}\\\
\times\int_{0}^{\theta_{m}}\frac{\theta^{2}}{\theta^{2}+\gamma^{-2}}J_{1}(\zeta\theta)d\theta$
(65)

## References

  * Tajima and Dawson [1979] T. Tajima and J. M. Dawson, Physical Review Letters 43, 267 (1979).
  * Geddes  _et al._ [2004] C. Geddes, C. Toth, J. Van Tilborg, E. Esarey, C. Schroeder, D. Bruhwiler, C. Nieter, J. Cary, and W. Leemans, Nature 431, 538 (2004).
  * Faure  _et al._ [2004] J. Faure, Y. Glinec, A. Pukhov, S. Kiselev, S. Gordienko, E. Lefebvre, J.-P. Rousseau, F. Burgy, and V. Malka, Nature 431, 541 (2004).
  * Mangles  _et al._ [2004] S. P. Mangles, C. Murphy, Z. Najmudin, A. G. R. Thomas, J. Collier, A. E. Dangor, E. Divall, P. Foster, J. Gallacher, C. Hooker, _et al._ , Nature 431, 535 (2004).
  * Leemans  _et al._ [2006] W. P. Leemans, B. Nagler, A. J. Gonsalves, C. Toth, K. Nakamura, C. G. Geddes, E. Esarey, C. Schroeder, and S. Hooker, Nature physics 2, 696 (2006).
  * Leemans  _et al._ [2014] W. Leemans, A. Gonsalves, H.-S. Mao, K. Nakamura, C. Benedetti, C. Schroeder, C. Tóth, J. Daniels, D. Mittelberger, S. Bulanov, _et al._ , Physical review letters 113, 245002 (2014).
  * Gonsalves  _et al._ [2019] A. Gonsalves, K. Nakamura, J. Daniels, C. Benedetti, C. Pieronek, T. de Raadt, S. Steinke, J. Bin, S. Bulanov, J. van Tilborg, _et al._ , Physical Review Letters 122, 084801 (2019).
  * Wang  _et al._ [2013] X. Wang, R. Zgadzaj, N. Fazel, Z. Li, S. Yi, X. Zhang, W. Henderson, Y.-Y. Chang, R. Korzekwa, H.-E. Tsai, _et al._ , Nature communications 4, 1988 (2013).
  * Kim  _et al._ [2013] H. T. Kim, K. H. Pae, H. J. Cha, I. J. Kim, T. J. Yu, J. H. Sung, S. K. Lee, T. M. Jeong, and J. Lee, Physical review letters 111, 165002 (2013).
  * Mirzaie  _et al._ [2015] M. Mirzaie, S. Li, M. Zeng, N. Hafz, M. Chen, G. Li, Q. Zhu, H. Liao, T. Sokollik, F. Liu, _et al._ , Scientific reports 5, 14659 (2015).
  * Wang  _et al._ [2016] W. Wang, W. Li, J. Liu, Z. Zhang, R. Qi, C. Yu, J. Liu, M. Fang, Z. Qin, C. Wang, _et al._ , Physical review letters 117, 124801 (2016).
  * Wang  _et al._ [2021] W. Wang, K. Feng, L. Ke, C. Yu, Y. Xu, R. Qi, Y. Chen, Z. Qin, Z. Zhang, M. Fang, _et al._ , Nature 595, 516 (2021).
  * Albert  _et al._ [2014] F. Albert, A. Thomas, S. Mangles, S. Banerjee, S. Corde, A. Flacco, M. Litos, D. Neely, J. Vieira, Z. Najmudin, _et al._ , Plasma Physics and Controlled Fusion 56, 084015 (2014).
  * Huang  _et al._ [2016] K. Huang, Y. Li, D. Li, L. Chen, M. Tao, Y. Ma, J. Zhao, M. Li, M. Chen, M. Mirzaie, _et al._ , Scientific reports 6, 27633 (2016).
  * Yan  _et al._ [2017] W. Yan, C. Fruhling, G. Golovin, D. Haden, J. Luo, P. Zhang, B. Zhao, J. Zhang, C. Liu, M. Chen, _et al._ , Nature Photonics 11, 514 (2017).
  * Nakajima [2008] K. Nakajima, Nature physics 4, 92 (2008).
  * Huang and Kim [2007] Z. Huang and K.-J. Kim, Phys. Rev. ST Accel. Beams 10, 034801 (2007).
  * Shan  _et al._ [2000] J. Shan, A. S. Weling, E. Knoesel, L. Bartels, M. Bonn, A. Nahata, G. A. Reider, and T. F. Heinz, Optics letters 25, 426 (2000).
  * Yan  _et al._ [2000] X. Yan, A. MacLeod, W. Gillespie, G. Knippels, D. Oepts, A. van der Meer, and W. Seidel, Physical Review Letters 85, 3404 (2000).
  * Wilke  _et al._ [2002] I. Wilke, A. M. MacLeod, W. A. Gillespie, G. Berden, G. Knippels, and A. Van Der Meer, Physical review letters 88, 124801 (2002).
  * Berden  _et al._ [2004] G. Berden, S. P. Jamison, A. M. MacLeod, W. A. Gillespie, B. Redlich, and A. Van der Meer, Physical review letters 93, 114802 (2004).
  * Berden  _et al._ [2007] G. Berden, W. A. Gillespie, S. P. Jamison, E.-A. Knabbe, A. M. MacLeod, A. Van Der Meer, P. J. Phillips, H. Schlarb, B. Schmidt, P. Schmüser, _et al._ , Physical review letters 99, 164801 (2007).
  * Steffen  _et al._ [2009] B. Steffen, V. Arsov, G. Berden, W. A. Gillespie, S. Jamison, A. M. MacLeod, A. Van Der Meer, P. Phillips, H. Schlarb, B. Schmidt, _et al._ , Physical Review Special Topics-Accelerators and Beams 12, 032802 (2009).
  * Cavalieri  _et al._ [2005] A. L. Cavalieri, D. Fritz, S. Lee, P. Bucksbaum, D. Reis, J. Rudati, D. Mills, P. Fuoss, G. Stephenson, C. Kao, _et al._ , Physical review letters 94, 114801 (2005).
  * Yang  _et al._ [2009] X. Yang, T. Tsang, T. Rao, J. Murphy, Y. Shen, and X. Wang, Applied Physics Letters 95, 231106 (2009).
  * Wang  _et al._ [2017] W. Wang, Y. Du, L. Yan, Z. Chi, Z. Zhang, J. Hua, W. Huang, C. Tang, and M. Li, Physical Review Accelerators and Beams 20, 112801 (2017).
  * Scoby  _et al._ [2010] C. Scoby, P. Musumeci, J. Moody, and M. Gutierrez, Physical Review Special Topics-Accelerators and Beams 13, 022801 (2010).
  * Jamison  _et al._ [2006] S. Jamison, G. Berden, A. M. MacLeod, D. Jaroszynski, B. Redlich, A. Van der Meer, and W. A. Gillespie, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 557, 305 (2006).
  * Steffen [2007] B. R. Steffen,  _Electro-optic methods for longitudinal bunch diagnostics at FLASH_ , Tech. Rep. (2007).
  * Casalbuoni  _et al._ [2008] S. Casalbuoni, H. Schlarb, B. Schmidt, P. Schmüser, B. Steffen, and A. Winter, Physical Review Special Topics-Accelerators and Beams 11, 072802 (2008).
  * van Tilborg  _et al._ [2006] J. van Tilborg, C. B. Schroeder, C. V. Filip, C. Tóth, C. G. R. Geddes, G. Fubiani, R. Huber, R. A. Kaindl, E. Esarey, and W. P. Leemans, Phys. Rev. Lett. 96, 014801 (2006).
  * Debus  _et al._ [2010] A. D. Debus, M. Bussmann, U. Schramm, R. Sauerbrey, C. D. Murphy, Z. Major, R. Hörlein, L. Veisz, K. Schmid, J. Schreiber, K. Witte, S. P. Jamison, J. G. Gallacher, D. A. Jaroszynski, M. C. Kaluza, B. Hidding, S. Kiselev, R. Heathcote, P. S. Foster, D. Neely, E. J. Divall, C. J. Hooker, J. M. Smith, K. Ertel, A. J. Langley, P. Norreys, J. L. Collier, and S. Karsch, Phys. Rev. Lett. 104, 084802 (2010).
  * Tian  _et al._ [2020] Q. Tian, Y. Du, H. Xu, Y. Liang, Y. Wang, Y. Tan, L. Yan, R. Li, C. Cheng, W. Huang, and C. Tang, Phys. Rev. Accel. Beams 23, 102802 (2020).
  * Xiang and Huang [2007] D. Xiang and W.-H. Huang, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 570, 357 (2007).
  * Castellano and Verzilov [1998] M. Castellano and V. A. Verzilov, Phys. Rev. ST Accel. Beams 1, 062801 (1998).
  * Casalbuoni  _et al._ [2009] S. Casalbuoni, B. Schmidt, P. Schmüser, V. Arsov, and S. Wesch, Physical Review Special Topics-Accelerators and Beams 12, 030705 (2009).
  * Loos  _et al._ [2008] H. Loos, R. Akre, A. Brachmann, F. Decker, Y. Ding, D. Dowell, P. Emma, J. Frisch, S. Gilevich, G. Hays, _et al._ , in  _Proceedings of the 2008 Free-Electron Laser Conference (Gyeongju, Korea)_ (2008).
  * Potylitsyn  _et al._ [2020] A. Potylitsyn, L. Sukhikh, T. Gusvitskii, G. Kube, and A. Novokshonov, Physical Review Accelerators and Beams 23, 042804 (2020).
  * Goodman [1988] J. Goodman, Introduction to fourier optics, classic textbook reissue (1988).
  * Schroeder  _et al._ [2004] C. B. Schroeder, E. Esarey, J. van Tilborg, and W. P. Leemans, Phys. Rev. E 69, 016501 (2004).
  * Lumpkin  _et al._ [2020] A. H. Lumpkin, M. LaBerge, D. W. Rule, R. Zgadzaj, A. Hannasch, O. Zarini, B. Bowers, A. Irman, J. P. Couperus Cabadağ, A. Debus, A. Köhler, U. Schramm, and M. C. Downer, Phys. Rev. Lett. 125, 014801 (2020).
  * Schmidt  _et al._ [2020] B. Schmidt, N. M. Lockmann, P. Schmüser, and S. Wesch, Phys. Rev. Accel. Beams 23, 062801 (2020).
  * Huang  _et al._ [2019] K. Huang, H. Kotaki, M. Mori, T. Esirkepov, J. K. Koga, Y. Hayashi, N. Nakanii, S. V. Bulanov, and M. Kando, Phys. Rev. Accel. Beams 22, 121301 (2019).
  * Beringer  _et al._ [2012] J. Beringer, J. F. Arguin, R. M. Barnett, K. Copic, O. Dahl, D. E. Groom, C. J. Lin, J. Lys, H. Murayama, C. G. Wohl, W. M. Yao, P. A. Zyla, C. Amsler, M. Antonelli, D. M. Asner, H. Baer, H. R. Band, T. Basaglia, C. W. Bauer, J. J. Beatty, V. I. Belousov, E. Bergren, G. Bernardi, W. Bertl, S. Bethke, H. Bichsel, O. Biebel, E. Blucher, S. Blusk, G. Brooijmans, O. Buchmueller, R. N. Cahn, M. Carena, A. Ceccucci, D. Chakraborty, M. C. Chen, R. S. Chivukula, G. Cowan, G. D’Ambrosio, T. Damour, D. de Florian, A. de Gouvêa, T. DeGrand, P. de Jong, G. Dissertori, B. Dobrescu, M. Doser, M. Drees, D. A. Edwards, S. Eidelman, J. Erler, V. V. Ezhela, W. Fetscher, B. D. Fields, B. Foster, T. K. Gaisser, L. Garren, H. J. Gerber, G. Gerbier, T. Gherghetta, S. Golwala, M. Goodman, C. Grab, A. V. Gritsan, J. F. Grivaz, M. Grünewald, A. Gurtu, T. Gutsche, H. E. Haber, K. Hagiwara, C. Hagmann, C. Hanhart, S. Hashimoto, K. G. Hayes, M. Heffner, B. Heltsley, J. J. Hernández-Rey, K. Hikasa, A. Höcker, J. Holder, A. Holtkamp, J. Huston, J. D. Jackson, K. F. Johnson, T. Junk, D. Karlen, D. Kirkby, S. R. Klein, E. Klempt, R. V. Kowalewski, F. Krauss, M. Kreps, B. Krusche, Y. V. Kuyanov, Y. Kwon, O. Lahav, J. Laiho, P. Langacker, A. Liddle, Z. Ligeti, T. M. Liss, L. Littenberg, K. S. Lugovsky, S. B. Lugovsky, T. Mannel, A. V. Manohar, W. J. Marciano, A. D. Martin, A. Masoni, J. Matthews, D. Milstead, R. Miquel, K. Mönig, F. Moortgat, K. Nakamura, M. Narain, P. Nason, S. Navas, M. Neubert, P. Nevski, Y. Nir, K. A. Olive, L. Pape, J. Parsons, C. Patrignani, J. A. Peacock, S. T. Petcov, A. Piepke, A. Pomarol, G. Punzi, A. Quadt, S. Raby, G. Raffelt, B. N. Ratcliff, P. Richardson, S. Roesler, S. Rolli, A. Romaniouk, L. J. Rosenberg, J. L. Rosner, C. T. Sachrajda, Y. Sakai, G. P. Salam, S. Sarkar, F. Sauli, O. Schneider, K. Scholberg, D. Scott, W. G. Seligman, M. H. Shaevitz, S. R. Sharpe, M. Silari, T. Sjöstrand, P. Skands, J. G. Smith, G. F. Smoot, S. Spanier, H. Spieler, A. Stahl, T. Stanev, S. L. Stone, T. Sumiyoshi, M. J. Syphers, F. Takahashi, M. Tanabashi, J. Terning, M. Titov, N. P. Tkachenko, N. A. Törnqvist, D. Tovey, G. Valencia, K. van Bibber, G. Venanzoni, M. G. Vincter, P. Vogel, A. Vogt, W. Walkowiak, C. W. Walter, D. R. Ward, T. Watari, G. Weiglein, E. J. Weinberg, L. R. Wiencke, L. Wolfenstein, J. Womersley, C. L. Woody, R. L. Workman, A. Yamamoto, G. P. Zeller, O. V. Zenin, J. Zhang, R. Y. Zhu, G. Harper, V. S. Lugovsky, and P. Schaffner (Particle Data Group), Phys. Rev. D 86, 010001 (2012).
  * [45] https://pdg.lbl.gov/2009/AtomicNuclearProperties/.
  * [46] https://en.wikipedia.org/wiki/SAE_304_stainless_steel/.
  * Jackson [1999] J. D. Jackson, Classical electrodynamics (1999).
  * Kube [2008] G. Kube, TESLA-FEL Report 1, 2008 (2008).
  * Ter-Mikaelian [1972] M. L. Ter-Mikaelian,  _High-energy electromagnetic processes in condensed media_ , 29 (John Wiley & Sons, 1972).
  * Brau [2004] C. A. Brau,  _Modern problems in classical electrodynamics_ (Oxford univ. press, 2004).
  * Gradshteyn and Ryzhik [2014] I. S. Gradshteyn and I. M. Ryzhik,  _Table of integrals, series, and products_ (Academic press, 2014).