# Quantum effects in the interaction of low-energy electrons with light

Adamantios P. Synanidis  ICFO-Institut de Ciencies Fotoniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
P. A. D. Gonçalves  ICFO-Institut de Ciencies Fotoniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
Claus Ropers  Department of Ultrafast Dynamics, Max Planck Institute for
Multidisciplinary Sciences, 37077 Göttingen, Germany  4th Physical Institute -
Solids and Nanostructures, University of Göttingen, 37077 Göttingen, Germany
F. Javier García de Abajo<EMAIL_ADDRESS>ICFO-Institut
de Ciencies Fotoniques, The Barcelona Institute of Science and Technology,
08860 Castelldefels (Barcelona), Spain  ICREA-Institució Catalana de Recerca i
Estudis Avançats, Passeig Lluís Companys 23, 08010 Barcelona, Spain

###### Abstract

The interaction between free electrons and nanoscale optical fields has
emerged as a unique platform to investigate ultrafast processes in matter and
explore fundamental quantum phenomena. In particular, optically modulated
electrons are employed in ultrafast electron microscopy as noninvasive probes
that push the limits of spatiotemporal and spectral resolution down to the
picometer–attosecond–microelectronvolt range. Electron kinetic energies well
above the involved photon energies are commonly employed, rendering the
electron–light coupling efficiency low and, thus, only providing limited
access to the wealth of quantum nonlinear phenomena underlying the dynamical
response of nanostructures. Here, we theoretically investigate electron–light
interactions when photons and electrons have comparable energies, revealing
strong quantum and recoil effects that include a nonvanishing coupling of
surface-scattered electrons to plane waves of light, inelastic electron
backscattering from localized optical fields, and strong electron–light
coupling under grazing electron diffraction by an illuminated crystal surface.
Our results open new vistas in electron–light–matter interactions with
promising applications in ultrafast electron microscopy.

## Introduction

The synergetic relation between short light pulses and free electron beams
(e-beams) underlies several recent advances in ultrafast electron microscopy
toward a combined sub-nm–sub-fs–sub-meV spatiotemporal and spectral
resolution, rapidly progressing toward the goal of mapping atomic-scale
spatial features and their evolution over unprecedentedly small time scales
Barwick  _et al._ (2009); Flannigan and Zewail (2012); Feist  _et al._ (2015);
Barwick and Zewail (2015); Polman  _et al._ (2019); García de Abajo and Di
Giulio (2021); Roques-Carmes  _et al._ (2023); García de Abajo and Ropers
(2023). A prominent example is photon-induced near-field electron microscopy
Barwick  _et al._ (2009); García de Abajo  _et al._ (2010); Park  _et al._
(2010) (PINEM), which is based on the synchronous arrival of laser and
electron femtosecond pulses at a sampled nanostructure, thus enabling optical-
pump/electron-probe spectroscopy to be performed with a nanoscale spatial
resolution inherited from the use of state-of-the-art electron optics setups.
This approach has been applied to image optical near fields in nanophotonics
Piazza  _et al._ (2015); Lummen  _et al._ (2016); Vanacore  _et al._ (2020);
Kurman  _et al._ (2021), the subcycle evolution of those fields Nabben  _et
al._ (2023); Gaida  _et al._ (2024); Bucher  _et al._ (2023), and the
nanoscale-resolved fluctuations of the light with which the electron has
interacted Di Giulio  _et al._ (2019); Dahan  _et al._ (2021); Yang  _et al._
(2024). PINEM can be regarded as a specific application of the more general
concept of stimulated inelastic electron–light scattering (SIELS). The latter
has been leveraged to gain control over the free-electron wave function by
customizing its interaction with light, including the generation of trains of
attosecond electron pulses Feist  _et al._ (2015); Priebe  _et al._ (2017);
Kozák  _et al._ (2018); Morimoto and Baum (2018) and the shaping of the
transverse e-beam profile Vanacore  _et al._ (2019); Konečná and García de
Abajo (2020); García de Abajo and Konečná (2021); Madan  _et al._ (2022);
Mihaila  _et al._ (2022). Laser-assisted photoemission Glover  _et al._
(1996); Saathoff  _et al._ (2008) constitutes another example of SIELS that
can be used to probe the ultrafast dynamics of condensed matter systems Arrell
_et al._ (2016) and produces interesting effects in the strong-field limit
Yalunin  _et al._ (2011); Dombi  _et al._ (2020).

Many of these advances have been accomplished in transmission electron
microscopes (TEMs) operating with relatively high e-beam energies ($\gtrsim
30\,$keV), orders of magnitude larger than those of the employed photons
(typically in the eV range) and, consequently, rendering the probability that
a single electron interacts with a single photon (e.g., in an individual
confined optical mode) much smaller than unity. Such a weak electron–light
interaction limits applications in metrology, imaging of atomic-scale
excitations, and the study of nonlinear phenomena. For robust structures, the
problem is circumvented in SIELS by using intense laser pulses Piazza  _et
al._ (2015); Lummen  _et al._ (2016); Vanacore  _et al._ (2020), but this
approach cannot be extended to sensitive specimens such as biological
materials. Phase matching between the electron excitation and the light field
can also boost the interaction Bendaña  _et al._ (2011); Kfir  _et al._
(2020); Wang  _et al._ (2020); Dahan  _et al._ (2020); Henke  _et al._ (2021),
although this strategy is only practical in specialized structures that host
modes with evanescent optical fields in vacuum.

Surprisingly, in the linear electron–photon interaction regime, describing
electrons as classical point-like charges produces the same results as a
quantum-mechanical treatment in which electron recoil is ignored García de
Abajo and Di Giulio (2021). Likewise, the wave function of energetic electrons
in PINEM, and more generally in SIELS, is modified by a global factor that
encapsulates the interaction with light through a single complex parameter
that depends linearly on the optical field Varshalovich and D’Yakonov (1971);
Weingartshofer  _et al._ (1977); García de Abajo  _et al._ (2010); Park  _et
al._ (2010); Feist  _et al._ (2015); García de Abajo and Di Giulio (2021) and
is therefore unsuited to access the subcycle nonlinear dynamics of a specimen
in general. Consequently, it is highly desirable to achieve strong interaction
between single electrons and atomic-scale excitations as a way to access their
ultrafast nonlinear dynamics. In this context, the use of low-energy electrons
with kinetic energies comparable to those of the quanta associated with the
employed optical fields opens a plausible avenue to overcome these challenges.
Indeed, the electron–light coupling strength increases when the electron
energy is reduced, while a deviation from a classical-probe behavior is
introduced by electron recoil Talebi (2018, 2020); García de Abajo  _et al._
(2022). Low-energy electrons thus hold the potential to reveal new physical
phenomena during their interaction with localized optical fields, a
possibility that demands the exploration of physically relevant
configurations.

Figure 1: Electron–light–matter interaction and recoil effects with low-energy
electrons. (a) Illustration of inelastic electron scattering by the evanescent
optical field associated with propagating surface polaritons. (b, c)
Transmission (b) and reflection (c) electron spectra in the configuration of
(a) for electrons of incident energy $\hbar\varepsilon_{0}=10\,$eV and angle
$\theta_{\rm e}=45\degree$ combined with polaritons of electric-field
amplitude $E_{0}=6\times 10^{8}\,$V/m, effective refractive index $n_{\rm
eff}=50$, and energy $\hbar\omega=1\,$eV. For readability, spectral peaks are
broadened with a Lorentzian of 0.1 FWHM in $\Delta\varepsilon/\omega$, with
$\Delta\varepsilon=\varepsilon-\varepsilon_{0}$. (d) Nonvanishing interaction
between a plane wave of light and an electron reflected on a light-transparent
surface. (e) Despite the small coupling in (d) ($P_{1}<2\%$), strong recoil
effects are observed for $\varepsilon_{0}\sim\omega$ in the asymmetry factor
$(P_{1}-P_{-1})/(P_{1}+P_{-1})$, which vanishes in the classical regime
($\varepsilon_{0}\gg\omega$); we take $E_{0}=8\times 10^{7}\,$V/m,
$\theta_{\rm e}=45\degree$, $\theta_{\rm l}=90\degree$, and
$\hbar\omega=1\,$eV.

Here, we reveal a wealth of previously unexplored quantum and recoil effects
taking place during electron–light–matter interactions when the electron and
photon energies are comparable. Considering realistic frameworks that admit
rigorous semi-analytical treatments, we theoretically explore electron–light
interactions mediated by scattering at planar interfaces of either the light,
the electrons, or both (Fig. 1). We adopt optical fields in the form of either
externally incident plane waves or surface polaritons. Importantly, surface
scattering leads to symmetry breaking that enhances the electron–light
coupling, even allowing otherwise forbidden electron–photon interactions. In
particular, we show that low-energy electrons can be inelastically scattered
solely due to an evanescent optical field propagating along an electron-
transparent surface (Fig. 1a–c), including the emergence of a back-reflected
electron signal (Fig. 1b). In addition, we demonstrate that a plane-wave
electron reflected on a light-transparent surface can produce inelastically
reflected electrons because of the nonvanishing electron–photon coupling
originating in translational symmetry breaking of the out-of-plane electron
wave function (Fig. 1d–e). The resulting reflected-electron spectrum exhibits
substantial recoil effects (Fig. 1e). We further report on the possibility of
reaching the strong electron–light coupling regime with moderate light
intensities through Bragg scattering at planar atomic lattices, whereby the
interaction is boosted under Rayleigh anomaly conditionsLord Rayleigh (1907a).
Besides its interest from a fundamental viewpoint, the present study unveils
exciting opportunities for the exploration of improved microscopy and
metrology in the regime of low-energy electrons exposed to optical fields of
comparable photon energy. Our results are particularly relevant to the
exploration of the rich phase and electronic-structure phenomenology exhibited
by material surfaces, which are cornerstones in many technological
applications.

## Results

### Theoretical framework for the interaction of low-energy electrons with
light

We first consider a planar material acting through a one-dimensional (1D)
potential $V(z)$ on the electron, while a generalization to laterally
corrugated atomic lattices is presented further below. To study electron
scattering by the planar structure in the presence of a classical optical
field (e.g., light plane waves or surface polaritons), we write the
Hamiltonian $\hat{\mathcal{H}}_{0}({\bf r})+\hat{\mathcal{H}}_{1}({\bf r},t)$,
where $\hat{\mathcal{H}}_{0}({\bf r})=-\hbar^{2}\nabla^{2}/2m_{\rm e}+V(z)$
describes the electron–material system and $\hat{\mathcal{H}}_{1}({\bf
r},t)=-({\rm i}e\hbar/m_{\rm e}c)\,{\bf A}({\bf r},t)\cdot\bm{\nabla}$
accounts for the electron–light interaction. The latter arises from the
minimal coupling prescription applied to a classical vector potential ${\bf
A}({\bf r},t)$ after neglecting $A^{2}$ terms and adopting a gauge in which
the scalar potential is zero.

We focus on monochromatic fields of frequency $\omega$ and in-plane wave
vector ${\bf k}_{\parallel}=(k_{x},k_{y})$, characterized by a vector
potential ${\bf A}({\bf r},t)={\bf A}(z)\,{\rm e}^{{\rm i}{\bf
k}_{\parallel}\cdot{\bf R}-{\rm i}\omega t}+{\rm c.c.}$, where the notation
${\bf R}=(x,y)$ is adopted. Also, in the absence of illumination, the electron
wave function is taken to be a solution of the $\hat{\mathcal{H}}_{0}({\bf
r})$ Hamiltonian, namely, $\psi_{0}({\bf r},t)=\varphi_{00}(z)\,{\rm e}^{{\rm
i}{\bf q}_{0\parallel}\cdot{\bf R}-{\rm i}\varepsilon_{0}t}$, having a well-
defined energy $\hbar\varepsilon_{0}$ and in-plane wave vector ${\bf
q}_{0\parallel}$ ($\perp\hat{\bf z}$). Now, in-plane translation symmetry and
energy conservation allow us to write the perturbation series

$\displaystyle\psi({\bf
r},t)=\sum_{n=0}^{\infty}\sum_{\ell=-n}^{n}\varphi_{n\ell}(z)\,{\rm e}^{{\rm
i}({\bf q}_{0\parallel}+\ell\,{\bf k}_{\parallel})\cdot{\bf R}-{\rm
i}(\varepsilon_{0}+\ell\omega)t},$ (1)

where $n$ runs over scattering orders, while $\ell$ denotes the net number of
exchanged photons (i.e., absorbed or emitted by the electron for $\ell>0$ and
$\ell<0$, respectively). In addition, $\psi({\bf r},t)$ satisfies the
Lippmann–Schwinger equation $\psi({\bf r},t)=\psi_{0}({\bf r},t)+\int{\rm
d}^{3}{\bf r}^{\prime}\int{\rm d}t^{\prime}\,\mathcal{G}_{0}({\bf r},{\bf
r}^{\prime},t-t^{\prime})\,\hat{\mathcal{H}}_{1}({\bf
r}^{\prime},t)\,\psi({\bf r}^{\prime},t^{\prime})$, where the Green function
$\mathcal{G}_{0}({\bf r},{\bf r}^{\prime},t-t^{\prime})$ is defined by
$[\hat{\mathcal{H}}_{0}({\bf r})-{\rm i}\hbar\partial_{t}]\mathcal{G}_{0}({\bf
r},{\bf r}^{\prime},t-t^{\prime})=-\delta({\bf r}-{\bf
r}^{\prime})\delta(t-t^{\prime})$. Combining these elements, we obtain the
recurrence relation

$\displaystyle\varphi_{n\ell}(z)=\frac{\hbar e}{m_{\rm e}c}\int{\rm
d}{z^{\prime}}\,\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp}_{\ell})$ (2)
$\displaystyle\times\Big{\\{}{\bf A}(z^{\prime})\cdot\big{[}{\bf
q}_{0\parallel}+(\ell-1)\,{\bf k}_{\parallel}-{\rm i}\hat{\bf
z}\,\partial_{z^{\prime}}\big{]}\,\varphi_{n-1,\ell-1}(z^{\prime})$
$\displaystyle\;\;\,+{\bf A}^{*}(z^{\prime})\cdot\big{[}{\bf
q}_{0\parallel}+(\ell+1)\,{\bf k}_{\parallel}-{\rm i}\hat{\bf
z}\,\partial_{z^{\prime}}\big{]}\,\varphi_{n-1,\ell+1}(z^{\prime})\Big{\\}}$

for $n>0$, where $\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp}_{\ell})$ is
the frequency-domain 1D Green function satisfying $\mathcal{G}_{0}({\bf
r},{\bf r}^{\prime},t-t^{\prime})=(2\pi)^{-3}\int{\rm d}^{2}{\bf
q}_{\parallel}\int{\rm
d}\varepsilon\,\mathcal{G}_{0}(z,z^{\prime},\varepsilon)\,{\rm e}^{{\rm i}{\bf
q}_{\parallel}\cdot({\bf R}-{\bf R}^{\prime})-{\rm
i}\varepsilon(t-t^{\prime})}$, and we define
$\varepsilon^{\perp}_{\ell}=\varepsilon_{0}+\ell\omega-\hbar|{\bf
q}_{0\parallel}+\ell\,{\bf k}_{\parallel}|^{2}/2m_{\rm e}$ such that
$\hbar\varepsilon^{\perp}_{\ell}$ is the out-of-plane electron energy after
the exchange of a net number of photons $\ell$ (see Appendix for a self-
contained derivation).

We consider planar structures of negligible thickness and expand the electron
wave function components in Eq. (1) using the ansatz

$\displaystyle\varphi_{n\ell}(z)=\sum_{s=\pm}\sum_{j}\alpha_{n\ell}^{js}\,{\rm
e}^{\zeta_{n\ell}^{js}z}\,\Theta(sz)$ (3)

within the regions above ($s=+$) and below ($s=-$) the material. Here, $j$
labels contributions that can be either evanescent (${\rm
Re}\\{\zeta_{n\ell}^{j\pm}\\}\neq 0$) or propagating (imaginary
$\zeta_{n\ell}^{j\pm}$). Inserting Eq. (3) into Eq. (2), we find a recursive
expression with a unique solution for the coefficients $\alpha_{n\ell}^{j\pm}$
and $\zeta_{n\ell}^{j\pm}$, which automatically satisfy the physical
conditions ${\rm Re}\\{\zeta_{n\ell}^{j+}\\}\leq 0$ and ${\rm
Re}\\{\zeta_{n\ell}^{j-}\\}\geq 0$ (see SI). In addition, the purely
propagating components have exponential coefficients
$\zeta_{n\ell}^{j\pm}=\pm{\rm i}q_{\ell z}$ (i.e., imaginary and independent
of $n$ and $j$), where $q_{\ell z}=\sqrt{2m_{\rm
e}\varepsilon^{\perp}_{\ell}/\hbar}$ is determined by energy conservation for
a net number of photon exchanges $\ell$. Finally, the fractions of
$\ell$-resolved electrons scattered along the upward ($+$) and downward ($-$)
directions are given by

$\displaystyle P_{\ell}^{\pm}=\frac{q_{\ell
z}}{q_{0z}}\,\Big{|}\sideset{}{{}^{\prime}}{\sum}_{nj}\alpha_{n\ell}^{j\pm}\Big{|}^{2},$
(4)

where the primed sum indicates that it is restricted to purely propagating
waves.

Figure 2: Inelastic scattering of low-energy electrons upon total reflection
at a polariton-supporting surface. (a–c,e–g) Probabilities corresponding to a
net exchange of $\ell={0,\pm 1}$ quanta calculated without (a–c) and with
(e–g) inclusion of quantum recoil as a function of electron incidence angle
$\theta_{\rm e}$. (d, h) Electron spectra for $\theta_{\rm e}=86.2\degree$
obtained without (d) and with (h) recoil. We consider a polaritonic electric-
field amplitude $E_{0}=5\times 10^{6}\,$V/m, effective refractive index
$n_{\rm eff}=250$, and energy $\hbar\omega=0.2\,$eV.

### Recoil and quantum effects in surface-scattered electrons

Fully electron-reflecting surface.—As a tutorial configuration, we first
consider an electron with energy $\hbar\varepsilon_{0}$ and incident angle
$\theta_{\rm e}$ (with respect to the surface normal) that undergoes total
reflection at a planar surface supporting a surface polariton with in-plane
wave vector $k_{\parallel}=|{\bf k}_{\parallel}|=n_{\rm eff}k$, where
$k=\omega/c$ is the light wave vector and $n_{\rm eff}>1$ is an effective
index of refraction. For a thin metallic film of permittivity $\epsilon<0$ and
thickness $d$ embedded in a medium of permittivity $\epsilon_{s}>1$, we have
$n_{\rm eff}=\epsilon_{s}\lambda/[\pi(1-\epsilon)d]$ at a light wavelength
$\lambda$, so thin films favor large values of $n_{\rm eff}$ (e.g., $n_{\rm
eff}=40$ for $\lambda=500$ nm in currently available Abd El-Fattah  _et al._
(2019) 10-monolayer crystalline Ag(111) films deposited on Si). Higher values
of $n_{\rm eff}\sim 100$’s are displayed by infrared graphene plasmons
Woessner  _et al._ (2015); Ni  _et al._ (2018) and phonon-polaritons in few-
layer hBN Giles  _et al._ (2018); Li  _et al._ (2021).

For simplicity, we take the electron and the surface polariton to share the
same in-plane direction of incidence with ${\bf q}_{0\parallel}\parallel{\bf
k}_{\parallel}\parallel\hat{\bf x}$, such that the associated vector potential
can be written as ${\bf A}({\bf
r})=(E_{0}/k)\,\big{(}\kappa^{2}+k_{\parallel}^{2}\big{)}^{-1/2}\,\big{(}\kappa\,\hat{\bf
x}+{\rm i}\,k_{\parallel}\,\operatorname{sign}\\{z\\}\,\hat{\bf
z}\big{)}\,{\rm e}^{{\rm i}k_{\parallel}x-\kappa|z|}$, where $E_{0}$ is a
global electric-field amplitude (see SI). In addition, the electron–surface
interaction is assumed to be elastic, and thus, any inelastic electron signal
stems from surface-polariton emission and absorption processes by the
electron. In this scenario, we can set the $z$ component of the zeroth-order
electron wave function as $\varphi_{00}(z)=\big{[}{\rm e}^{-{\rm
i}q_{0z}z}-e^{{\rm i}q_{0z}z}\big{]}\Theta(z)$, while the Green function in
Eq. (2) reduces to
$\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp}_{\ell})=({\rm i}m_{\rm
e}/\hbar^{2}q_{\ell z})\big{[}{\rm e}^{{\rm i}q_{\ell z}(z+z^{\prime})}-{\rm
e}^{{\rm i}q_{\ell z}|z-z^{\prime}|}\big{]}\Theta(z)\Theta(z^{\prime})$ (see
SI). Inserting these elements into Eq. (2) and noticing that only reflected
components need to be considered, we find a set of analytical coefficients
$\alpha_{n,\ell}^{j+}$ and $\zeta_{n,\ell}^{j+}$, from which the reflection
probability $P_{\ell}\equiv P^{+}_{\ell}$ for a given $\ell$ channel is
obtained via Eq. (4).

It is instructive to examine the $\varepsilon_{0}\gg\omega$ limit, where
recoil effects should play a minor role. As a direct generalization of the
result obtained for an electron moving with constant velocity along a
straight-line trajectory García de Abajo and Di Giulio (2021), we approximate
the $\ell$-dependent inelastic probability as
$P_{\ell}=J_{\ell}^{2}(2|\beta|)$, where $\beta=({\rm i}e/\hbar c)\int{\rm
d}{t}\;\dot{{\bf r}}_{\rm e}(t)\cdot{\bf A}[{\bf r}_{\rm e}(t)]\,{\rm
e}^{-{\rm i}\omega t}$ is an electron–light coupling parameter obtained by
integrating over time the vector potential component parallel to the velocity
$\dot{{\bf r}}_{\rm e}(t)$ and evaluated at the electron position ${\bf
r}_{\rm e}(t)$. Taking a specularly reflected trajectory with the surface-
polariton field given above, we obtain

$\displaystyle\beta=\frac{2\,{\rm i}}{\sqrt{2n_{\rm
eff}^{2}-1}}\;\frac{eE_{0}c}{\hbar\omega^{2}}\frac{(n_{\rm
eff}c/v-\sin\theta_{\rm e})\cos\theta_{\rm e}}{(c/v-n_{\rm eff}\sin\theta_{\rm
e})^{2}+(n_{\rm eff}^{2}-1)\cos^{2}\theta_{\rm e}}.$ (5)

For reference, we note that the prefactor $eE_{0}c/\hbar\omega^{2}$ takes a
value of $\approx 2$ for $E_{0}=10^{7}$ V/m and $\hbar\omega=1$ eV. This mode
energy is characteristic of surface plasmons in ultrathin metal films Abd El-
Fattah  _et al._ (2019) and exciton-polaritons in transition-metal
dichalcogenides Li  _et al._ (2014); Epstein  _et al._ (2020), while
polaritons of lower energy (e.g., $\hbar\omega\sim 0.1$ eV in graphene
Woessner  _et al._ (2015); Ni  _et al._ (2018) and hBN Giles  _et al._ (2018);
Li  _et al._ (2021)) should produce larger coupling for the same field
amplitude in accordance to the scaling $\beta\propto 1/\omega^{2}$. Equation
(5) reveals the important role of confinement in enhancing the
electron–polariton coupling: given a certain electron velocity $v$, the
denominator reaches its minimum value under the condition

$\displaystyle c/n_{\rm eff}=v\sin\theta_{\rm e}$ (6)

(i.e., when the surface-polariton phase velocity matches the in-plane
projection of the electron velocity). In addition, Eq. (5) illustrates the
well-known linear scaling of the coupling coefficient with the applied
electric field amplitude $E_{0}$.

Figure 2 highlights the importance of recoil effects in the interaction
between a low-energy electron and a strongly confined surface polariton by
contrasting the classical nonrecoil PINEM theory [Fig. 2a–d, based on Eq. (5)]
with the rigorous quantum formalism introduced above [Fig. 2e–h, Eq. (4)]. As
a first observation, the classical treatment in Eq. (5) provides the necessary
conditions for strong coupling and leads to substantial inelastic
probabilities (Fig. 2a-d). Focusing for concreteness on the $\ell=\pm 1$
channels (see SI for more sidebands), the parameter space for which
electron–light coupling is maximized is well captured by the classical
framework, and it corresponds to the phase-matching condition in Eq. (6)
(i.e., $\hbar\varepsilon_{0}\approx m_{\rm e}c^{2}/2n_{\rm eff}^{2}\sim 4$ eV
with $\theta_{\rm e}=90\degree$), for which the coupling diverges as
$\beta\propto 1/\cos{\theta_{\mathrm{e}}}$ near $\theta_{\rm e}=90\degree$.
However, both the intensity profile in the $(\varepsilon_{0},\theta_{\rm e})$
phase space and the magnitude of the electron–light coupling strength are
markedly different in the classical and quantum theories. Specifically, the
incorporation of recoil in the latter leads to asymmetric loss–gain spectra
(i.e., $P_{-\ell}\neq P_{\ell}$) as well as abrupt thresholds in $P_{\ell}$
(Fig. 2e–g). These observations can be interpreted by noting that
energy–momentum conservation imposes the condition $q_{\ell
z}=q_{0}\sqrt{1+\ell\omega/\varepsilon_{0}-|\sin\theta_{\rm e}+\ell
k_{\parallel}/q_{0}|^{2}}$, and thus, not only the emission ($\ell<0$) and
absorption ($\ell>0$) probabilities are rendered different, but also the
inelastic signal vanishes whenever
$1+\ell\omega/\varepsilon_{0}<|\sin\theta_{\rm e}+\ell
k_{\parallel}/q_{0}|^{2}$, since $q_{\ell z}$ becomes purely imaginary (i.e.,
the corresponding electron wave is evanescent). Another dramatic consequence
of recoil is the redistribution of probability to neighboring $\ell$’s near
the aforementioned thresholds (see Fig. 2h), in contrast to the symmetric
spectrum produced by the classical description (Fig. 2d).

Figure 3: Inelastic scattering of low-energy electrons upon partial reflection
at a polariton-supporting thin film. We show transmitted (a–d) and reflected
(e-h) electron spectra as a function of polaritonic electric-field amplitude
$E_{0}$ for fixed effective refractive index ($n_{\rm eff}=50$) and energy
($\hbar\omega=1\,$eV). Selected incident electron energies
$\hbar\varepsilon_{0}$ are considered, while the incidence angle is fixed at
$\theta_{\rm e}=45\degree$. Electron–surface scattering is modeled through a
$\delta$-function potential of amplitude $U_{0}=1\,$eV nm. Spectral features
are broadened by a Lorentzian of 0.1 eV FWHM.

Partially electron-reflecting surface.—We expect partial transmission and
reflection when the electron is scattered by an atomically thin 2D material,
which we describe through a surface potential $V(z)=U_{0}\,\delta(z)$. The
parameter $U_{0}$ has units of energy times length, and arguing that an atomic
monolayer can be described by a barrier of finite thickness $d\lesssim 1$ nm
and internal potential $V_{0}$ in the eV range, we expect $U_{0}\approx
V_{0}d$ in the eV$\times$nm range. We also assume the material to support
long-lived polaritonic modes (e.g., phonon-polaritons in hBN Giles  _et al._
(2018) or plasmons in doped graphene Ni  _et al._ (2018)), so that they are
characterized by a real effective index $n_{\rm eff}$. The calculation of the
probabilities associated with the electron wave functions involved in the net
exchange of $\ell$ polariton quanta follows the same steps as in the above
scenario of full reflection, but now both reflected and transmitted electron
components are produced. For an electron prepared with an incident $-q_{0z}$
wave vector in the out-of-plane direction, the wave function in the absence of
illumination is given by $\varphi_{00}(z)=\big{[}{\rm e}^{-{\rm
i}q_{0z}z}+r_{0}\,{\rm e}^{{\rm i}q_{0z}z}\big{]}\Theta(z)+t_{0}\,{\rm
e}^{-{\rm i}q_{0z}z}\,\Theta(-z)$, where we use the transmission and
reflection coefficients $t_{\ell}=\big{(}1+{\rm i}m_{\rm
e}U_{0}/\hbar^{2}q_{\ell z}\big{)}^{-1}$ and $r_{\ell}=t_{\ell}-1$,
respectively. We use this result together with the Green function
$\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp}_{\ell})=-({\rm i}m_{\rm
e}/\hbar^{2}q_{\ell z})\Big{[}{\rm e}^{{\rm i}q_{\ell
z}|z-z^{\prime}|}+r_{\ell}\,{\rm e}^{{\rm i}q_{\ell
z}(|z|+|z^{\prime}|)}\Big{]}$ (see SI) to obtain the reflected and transmitted
$\ell$-resolved probabilities $P^{\pm}_{\ell}$ from Eq. (4) following the
formalism developed above. The results are plotted as a function of electric-
field amplitude for different electron energies in Fig. 3. Again, we find
recoil effects emerging through strong asymmetries in the electron spectra,
which, as anticipated, become more symmetric as the electron energy is
increased toward the $\varepsilon_{0}\gg\omega$ regime. These results are
qualitatively correct even when more involved $z$-dependent potentials are
considered (e.g., finite-thickness films) because the Green function outside
the material retains the same expression as above de Aguiar (1993), with
reflection and transmission coefficients depending on the details of the
potential.

### Nonvanishing interaction of surface-scattered electrons and unscattered
light plane waves

An interesting scenario is presented when a thin film is illuminated from the
far field and the scattered optical components are comparatively negligible,
so the electron mainly sees an external light plane wave. Electron–light
coupling can still take place under this configuration because the free-space
energy–momentum mismatch is broken by the fact that the electron is scattered
by the material. We regard this situation as the complementary of electron
shaping mediated by PINEM interaction (i.e., SIELS assisted by electron-
transparent, light-reflecting plates, in which the kinematic electron–photon
free-space coupling mismatch is circumvented by having light half-plane waves
instead of full plane waves Vanacore  _et al._ (2018)). We note that 2D
monolayers (e.g., graphene and hBN) are nearly transparent to light (e.g.,
$\sim 2.3\%$ absorption by graphene over a wide spectral range Mak  _et al._
(2008); Nair  _et al._ (2008)) and can thus be regarded as good candidates to
explore the interaction of light plane waves with surface-reflected low-energy
electrons.

Figure 4: Inelastic interaction between a light plane wave and surface-
scattered electrons. We plot the probability associated with $\ell=\pm 1$ net
photon exchanges as a function of the electron ($\theta_{\rm e}$) and photon
($\theta_{\rm l}$) incidence angles without (a,b) and with (c,d) inclusion of
recoil for an electron-opaque, light-transparent surface. We assume total-
electron reflection and take $\hbar\varepsilon_{0}=2\,$eV, $E_{0}=8\times
10^{7}\,$V/m, and $\hbar\omega=1\,$eV. In panel (d), kinematically forbidden
regions are shaded in gray.

We explore this idea in Fig. 4, where a fully reflected electron is considered
to be interacting with a freely propagating p-polarized light plane wave. The
analysis is analogous to that in Fig. 2, but using a plane-wave optical field
instead of a surface mode. We present the probabilities $P_{\ell=\pm 1}^{\pm}$
calculated to first-order in the light intensity as a function of light and
electron incidence angles, comparing classical (nonrecoil) and quantum (Eq.
(4), with recoil) descriptions. Incidentally, the coupling parameter in the
former is given by

$\displaystyle\beta=\frac{2eE_{0}c}{\hbar\omega^{2}}\frac{\cos\theta_{\rm
e}(\sin\theta_{\rm l}\;c/v-\sin{\theta_{\rm e}})}{(c/v-\sin\theta_{\rm
l}\sin\theta_{\rm e})^{2}-\cos^{2}\theta_{\rm l}\cos^{2}\theta_{\rm e}},$

where $E_{0}$ is the light-plane-wave amplitude. We find again that recoil
leads to asymmetric inelastic electron signals, as well as regions of the
$(\theta_{\rm l},\theta_{\rm e})$ parameter space in which electron–light
coupling becomes kinematically allowed or forbidden, accompanied by a transfer
of probability to the symmetric ($\ell\to-\ell$) channel. In contrast to the
interaction with surface polaritons, where phase-matching at grazing incidence
produced the strongest interaction, now phase-matching is forbidden (i.e., the
in-plane optical wave vector is always smaller than the in-plane electron wave
vector), rendering the coupling smaller, so we reduce the analysis to first-
order processes. The maximum coupling is observed at grazing light incidence
($\theta_{\rm l}=\pm 90\degree$) and normal electron incidence ($\theta_{\rm
e}=0$), which is consistent with the angular scaling of the classical coupling
coefficient as $\beta\propto\sin\theta_{l}\cos\theta_{\rm e}$ for $v\ll c$.
These conditions guarantee maximum overlap of the light electric field along
the electron trajectory.

Figure 5: Recoil effects in the interaction of light with lattice-diffracted
electrons. (a) We consider a normally incident electron undergoing diffraction
by an Au(111) monolayer (Au–Au distance $d\approx 0.288$ nm, with atomic bonds
along $y$) under grazing plane-wave light irradiation along $x$ (${\bf
k}_{\parallel}\parallel\hat{\bf x}$) with linear polarization in the $y-z$
plane as indicated by the angle $\phi_{\rm p}$. Each of the electron Bragg
diffraction orders (see yellow arrows for one of them) splits in energy and
direction of reflection/transmission upon exchange of a net number of photons
$\ell$. (b) Electron isoenergy contours after exchanging $\ell=-1$, 0, and 1
photons (blue, black, and red circles), superimposed on the reciprocal lattice
of the atomic monolayer ($4\pi/\sqrt{3}d$ distance between sites). We consider
71.4 eV incident electrons and 1 eV photons. (c-h) Intensity of transmitted
(c-e) and reflected (f-h) Bragg peaks for three different incident electron
energies (see labels above (c,f), (d,g), and (e,h)) upon exchange of
$\ell=-1$, 0, or 1 photons (1 eV energy). The light wave vector ${\bf
k}_{\parallel}$ is indicated in (c) and the electric field amplitude is
$E_{0}=2.5\times 10^{8}$ V/m with polarization set by $\phi_{\rm
p}=45\degree$. The area of the circles gives the fraction of electrons in each
Bragg peak (see log-scale legend). (i-k) Same as (e), but for varying
polarization angles $\phi_{\rm p}$ (see top labels).

### Light-assisted low-energy electron diffraction

Low-energy electrons with energies $\sim 10-500\,$eV are commonly used to
study the atomic structure of crystal surfaces in low-energy electron
diffraction Rocca (1995); Pendry (1974, 1984) (LEED) because they penetrate
only a few atomic layers and have de Broglie wavelengths commensurate with the
atomic spacings. In a related context, energy-resolved inelastic low-energy
electron surface scattering is also used to probe surface modes Claus  _et
al._ (1992); Nagao  _et al._ (2001, 2006), while ultrafast LEED grants one
access to time-resolved structural dynamics Gulde  _et al._ (2014);
Vogelgesang  _et al._ (2018). Here, we theoretically study the interaction of
surface-diffracted electrons with light plane waves and show the important
role played by recoil in the underlying electron–light coupling, including the
presence of lattice resonances that boost the interaction under Rayleigh
anomaly conditions Lord Rayleigh (1907a).

To illustrate electron–light–matter interactions in the presence of Bragg
diffraction, we consider a low-energy electron normally impinging on an
illuminated monolayer of gold atoms arranged in a (111) triangular lattice
with an Au–Au bond distance of $0.288\,$nm (Fig. 5a). In the absence of
external illumination, the diffracted electron wave function consists of
components with wave vectors given by ${\bf q}^{\pm}_{{\bf G}0}={\bf
G}\pm\sqrt{q_{0}^{2}-G^{2}}\,\hat{\bf z}$, where ${\bf G}$ are 2D reciprocal
lattice vectors, while the $+$ ($-$) sign corresponds to upward (downward)
electron motion relative to the atomic plane. Lattice scattering is elastic,
so all of these wave vectors have a magnitude $q_{0}=\sqrt{2m_{\rm
e}\varepsilon_{0}/\hbar}$ determined by the incident electron energy
$\hbar\varepsilon_{0}$. Diffracted electron plane waves with $G<q_{0}$
generate observable LEED spots, as determined by an Ewald sphere construction
(see Fig. 5), whereas waves with $G>q_{0}$ are evanescent. The latter do not
explicitly contribute to the far-field electron scattering probability, but
they have to be retained in the description of dynamical electron diffraction
by the atomic layer Rocca (1995); Pendry (1974, 1984). Starting from a
normally incident electron, the wave function associated with the incident and
scattered waves in the absence of illumination takes the form

$\displaystyle\psi_{0}({\bf r},t)=\Big{[}{\rm e}^{-{\rm
i}q_{0}z}+\sum_{\pm}B^{\pm}_{{\bf G}}{\rm e}^{{\rm i}{\bf q}^{\pm}_{{\bf
G}0}\cdot{\bf r}}\Theta(\pm z)\Big{]}\,{\rm e}^{-{\rm i}\varepsilon_{0}t},$
(7a) where $B^{\pm}_{{\bf G}}$ are scattering amplitudes.

Upon interaction with an incident light plane wave of wave vector ${\bf k}$,
every diffraction order (either propagating or evanescent) can exchange energy
with the light field in multiples of the photon energy and in-plane wave
vector (i.e., $\ell\hbar\omega$, and $\ell{\bf k}_{\parallel}$, respectively),
giving rise to diffracted components with wave vectors ${\bf q}^{\pm}_{{\bf
G}\ell}={\bf G}+\ell{\bf k}_{\parallel}\pm q_{{\bf G}\ell z}$, where $q_{{\bf
G}\ell z}=\sqrt{2m_{\rm e}(\varepsilon_{0}+\ell\omega)^{2}/\hbar-|{\bf
G}+\ell{\bf k}_{\parallel}|^{2}}\,\hat{\bf z}$, labeled by the direction of
motion ($+/-$ for upward/downward scattering) and the net lattice and photon
momentum exchanges ($\hbar{\bf G}$ and $\ell\hbar{\bf k}_{\parallel}$). The
total electron wave function takes the form

$\displaystyle\psi({\bf r},t)=\Big{[}{\rm e}^{-{\rm
i}q_{0}z}+\sum_{\pm}C^{\pm}_{{\bf G}\ell}{\rm e}^{{\rm i}{\bf q}^{\pm}_{{\bf
G}\ell}\cdot{\bf r}-{\rm i}\ell\omega t}\Theta(\pm z)\Big{]}\,{\rm e}^{-{\rm
i}\varepsilon_{0}t},$ (7b)

where the amplitudes $C_{{\bf G}\ell}^{\pm}$ are self-consistently determined
from an extension of LEED theory to incorporate the interaction with both the
atomic lattice and the optical field [see SI for details, including analytical
expressions for the coefficients $B^{\pm}_{{\bf G}}$ and $C^{\pm}_{{\bf
G}\ell}$ in Eqs. (7a) and (7b), respectively.] For simplicity, we limit our
analysis to first order in the electron–light interaction, but we incorporate
the interaction with the lattice to all orders.

As the electron energy increases, the number of diffracted spots also
increases because more points are inside the Ewald sphere. In the elastic part
($\ell=0$), a given diffraction order ${\bf G}$ is observed in the LEED
pattern when the electron energy exceeds a threshold energy
$\hbar^{2}G^{2}/2m_{\rm e}$. However, in the inelastic components, the changes
in electron energy and momentum enter the condition for far-field propagation.
For example, after exchanging $\ell$ photons, a previously evanescent order
may become propagating if $\hbar|{\bf G}+\ell{\bf k}_{\parallel}|^{2}/2m_{\rm
e}<\varepsilon_{0}+\ell\omega$, and likewise, a propagating Bragg-diffracted
beam may become evanescent (Fig 5b). These effects are more relevant when
$\varepsilon_{0}\sim\hbar G^{2}/2m_{\rm e}$. In addition, electrons scattered
at the onset of a diffraction order move under grazing conditions, so they
spend more time near the surface and, therefore, undergo a stronger
interaction with light. Consequently, by choosing the electron energy close to
the threshold of one of the ${\bf G}$ beams, we expect to increase the
coupling of diffracted electrons to light, emphasizing the importance of
recoil (see below).

This phenomenology is illustrated in the calculations presented in Fig. 5c-h
for three different incident electron energies (in separate columns) near the
threshold of the $(1,1)$ and its symmetry-equivalent diffraction spots. We
show the intensities of transmitted (Fig. 5c-e) and reflected (Fig. 5f-h)
beams. Following the conditions under which a maximum electron–light
interaction was observed in Fig. 4, we take the electron to be normally
impinging on the atomic plane, while a light plane wave is incident parallel
to the $x$ surface direction with linear polarization as indicated in Fig. 5a.
The atomic lattice is oriented with Au–Au bonds along $y$. Upon inspection of
our numerical results, we estimate an overall inelastic scattering probability
of $\sim 10-20\%$ for incident electron energies up to 50 eV and optical
electric-field strengths $E_{0}\sim 2.5\times 10^{8}\,$V/m. Additionally, we
find that the waves corresponding to $\ell=\pm 1$ typically have comparable
intensities to the transmitted $\ell=0$ beam. In contrast to the amplitude
found when adopting the nonrecoil approximation, which only depends on the
amplitude of each ${\bf G}$-dependent LEED spot and its multiplexing into
different energy sidebands according to the corresponding electron–light
coupling coefficient $\beta$ (i.e., considering the time integral of the field
along the classical incoming and Bragg-reflected electron paths), a full
quantum treatment including recoil reveals that all diffraction orders
(propagating and evanescent) can contribute.

Figure 6: Strong coupling between light and diffracted electrons. (a)
Probability (blue curve, to first order) and azimuthal outgoing angle (orange
curves) of the $(2,1)$ diffraction spot [${\bf
G}=(2\pi/d)\,\big{(}\sqrt{3}\hat{\bf x}-\hat{\bf y}\big{)}$] after emitting
one photon ($\ell=-1$) as a function of electron energy $\hbar\varepsilon_{0}$
under the configuration of Fig. 5a for a light plane wave impinging along $x$
with polarization along $y$, amplitude $E_{0}=2.5\times 10^{7}$ V/m, and
photon energy $\hbar\omega=1$ eV. Rigorous first-order theory (solid-blue
curve) is compared with the analytical approximation in Eq. (8) (dashed
curve). (b) Same as (a), but shown as a function of photon energy for an
incident electron energy $\hbar\varepsilon_{0}=73.3$ eV. A divergent
probability is observed when the inelastically diffracted beam becomes grazing
at electron and photon energies $\hbar\varepsilon_{0{\rm r}}$ and
$\hbar\omega_{\rm r}$ in (a) and (b), respectively. The probability also
diverges at low frequencies with the asymptotic behavior $\propto
1/\omega^{4}$ indicated by the dotted-blue curve in (b).

The onset of a new diffraction order during the scattering of waves by a
periodic structure causes an anomaly consisting of the depletion of the
specularly reflected and directly transmitted beams, as pointed out by Lord
Rayleigh in the context of light diffraction by periodic gratings Lord
Rayleigh (1907a, b); García de Abajo (2007). In light-assisted inelastic
electron diffraction, a related anomaly takes place when a scattered beam
becomes grazing (i.e., a vanishing out-of-plane wave vector component $q_{{\bf
G}\varepsilon_{\ell}z}=0$ for a combination of reciprocal lattice vector ${\bf
G}$ and sideband order $\ell$, or equivalently,
$\hbar\varepsilon_{0}+\ell\hbar\omega=(\hbar^{2}/2m_{\rm e})|{\bf G}+\ell{\bf
k}_{\parallel}|^{2}$ under normal electron-incidence conditions). Upon
examination of the corresponding coefficients $C_{{\bf G}\ell}^{\pm}$ in Eq.
(7b) (see SI), considering $q_{{\bf G}\varepsilon_{\ell}z}\approx 0$ and
taking the field amplitude ${\bf E}_{0}\perp\hat{\bf z}$ for simplicity, we
can approximate

$\displaystyle C_{{\bf G},\ell}^{\pm}\approx$ $\displaystyle\frac{{\rm
i}^{\ell}e}{\hbar\omega}\;\frac{{\bf E}_{0}\cdot{\bf G}}{q_{{\bf
G}\varepsilon_{\ell}z}\;q_{{\bf G}\varepsilon_{0}z}}\Big{(}B^{+}_{{\bf
G}}+B^{-}_{{\bf G}}\Big{)},$ (8)

and consequently, the probability $\big{|}C_{{\bf G},\ell}^{\pm}\big{|}^{2}$
diverges as $1/q_{{\bf G}\varepsilon_{\ell}z}^{2}$. An illustrative example is
presented in Fig. 6 when varying either the electron energy (Fig. 6a) or the
photon energy (Fig. 6b) around the conditions for the grazing emission of an
inelastic electron beam. A divergence in the probability calculated to the
first order of interaction with the light is observed, leading to unphysical
values above unity. This indicates that the system enters into the
nonpertubative regime, requiring higher orders of interaction to describe the
system, but also revealing that strong electron–light coupling can be reached
even for small light intensities. Such a strong interaction results from in-
phase scattering by a large number of atoms in the planar lattice, which
demands the use of sufficiently broad electron and light beams that can be
regarded as plane waves over a wide surface area. This divergence is well
captured by Eq. (8) [Fig. 6, dashed-blue curves], in reasonable agreement with
our rigorous first-order results [Fig. 6, solid-blue curves, obtained from Eq.
(B34) in Appendix]. We note that, in addition to the explicit $1/\omega$
factor in Eq. (8), the inelastic scattering coefficient $C_{{\bf
G},\ell}^{\pm}$ is dominated by a $1/(q_{{\bf G}\varepsilon_{l}z}-q_{{\bf
G}\varepsilon_{0}z})$ term in Eq. (B34) (see SI), thus resulting in an overall
$\propto 1/\omega^{4}$ scaling of the probability with decreasing photon
frequency $\omega$, as indicated in Fig. 6b.

## Discussion

In summary, based on a comprehensive theoretical treatment of the quantum-
mechanical interaction between low-energy free electrons and optical fields
with comparable photon energies, we have identified a plethora of recoil and
quantum effects emerging in the form of dramatic modifications in the energy
and angular distribution of electrons undergoing elastic surface scattering
and inelastic interaction with optical fields associated with either surface
polaritons or propagating light. In particular, we have shown that free
electrons interacting with evanescent optical fields can undergo classically
forbidden backscattering. Furthermore, the interaction between surface-
scattered electrons and unscattered light plane waves renders a nonzero
electron–light coupling due to the breaking of translational symmetry in the
electron wave function; we propose that suspended atomically thin layers may
provide suitable conditions (high transparency to light and large electron
scattering) to observe this effect.

As a common element in the interaction between free electrons and crystal
surfaces, we have incorporated Bragg diffraction, which leads to an interplay
between scattering by the atomic lattice and inelastic photon exchanges. The
latter can transform propagating diffraction orders into evanescent or the
other way around. Importantly, strong electron–light coupling is predicted at
the onset of an inelastically scattering electron beam, capitalizing on the
in-phase interaction with many atoms in the structure. Such a strong coupling
regime could be leveraged to optically shape free electrons using moderate
light intensities, potentially operating in the continuous-wave regime without
damaging the scattering material.

Although we have focused on planar surfaces, the concept of combining elastic
scattering of low-energy electrons by a material structure and the inelastic
interaction with light of comparably low photon energy is more general and
could involve the use of non-periodic nanostructures such as holes, tips, and
other curved elements to guide and reshape the electron wave function and also
increase or spatially modulate its interaction with specific optical modes. In
a related context, electron interaction with illuminated atoms in the gas
phase has a long tradition Weingartshofer  _et al._ (1977, 1983); Francken and
C.J.Joachain (1990); E. Arqué López and García de Abajo (2022), which could be
revisited as a platform to optically modulate electrons and explore new
physics. The novel directions opened by the presented theory and simulations
could be experimentally explored in currently available low-energy electron-
microscope setups.

## Acknowledgments

The authors acknowledge insightful discussions with V. Di Giulio, S. V.
Yalunin, and J. Otto. This work has been supported in part by the European
Research Council (Advanced Grants 789104-eNANO and 101055435-ULEEM), the
European Commission (Horizon 2020 Grants No. 101017720 FET-Proactive EBEAM and
No. 964591-SMART-electron), the Spanish MICINN (PID2020-112625GB-I00 and
Severo Ochoa CEX2019-000910-S), the Catalan AGAUR (Grant No. 2023 FI-1 00052)
and CERCA Programs, and Fundaciós Cellex and Mir-Puig.

##  Appendix A Low-energy electron scattering by an illuminated homogeneous
surface

We present a self-contained theory of low-energy electron scattering at planar
interfaces that are either supporting surface polaritons or subject to
external illumination.

###  A1.1 Electron scattering by a planar surface

We consider an electron of energy $\hbar\varepsilon_{0}$ elastically scattered
by a homogeneous planar structure represented through a one-dimensional (1D)
potential energy $V(z)$. Different types of surfaces are discussed below, with
the electron impinging from the vacuum region ($z>0$) as a plane wave of well-
defined in-plane wave vector ${\bf q}_{0\parallel}$ and total wave vector
$q_{0}=\sqrt{2m_{\rm e}\varepsilon_{0}/\hbar}$. The incident out-of plane wave
vector is then $-q_{0z}$ with $q_{0z}=\sqrt{q_{0}^{2}-q_{0\parallel}^{2}}$.
The Hamiltonian of the system can be written as

$\displaystyle\hat{\mathcal{H}}_{0}({\bf
r})=-\frac{\hbar^{2}\nabla^{2}}{2m_{\rm e}}+V(z),$

with $V(z)=0$ in the $z>0$ region. Translational invariance in the $x-y$ plane
allows us to factorize the electron wave function as

$\displaystyle\psi_{0}({\bf r},t)={\rm e}^{{\rm i}{\bf
q}_{0\parallel}\cdot{\bf R}-{\rm i}\varepsilon_{0}t}\,\psi_{0}(z),$ (A9)

where ${\bf R}=(x,y)$ and the subscript in $\psi_{0}$ indicates that the
interaction with light is still not included.

In the formalism that follows, $\hat{\mathcal{H}}_{0}({\bf r})$ determines
$\psi_{0}(z)$ starting from an incident electron plane wave, as well as the
noninteracting Green function $\mathcal{G}_{0}({\bf r},{\bf
r}^{\prime},t-t^{\prime})$ implicitly defined by Economou (2006)

$\displaystyle\big{[}\hat{\mathcal{H}}_{0}({\bf r})-{\rm
i}\hbar\partial_{t}\big{]}\,\mathcal{G}_{0}({\bf r},{\bf
r}^{\prime},t-t^{\prime})=-\delta({\bf r}-{\bf
r}^{\prime})\delta(t-t^{\prime}).$ (A10)

Time invariance allows us to write the Green function in frequency space as
$\mathcal{G}_{0}({\bf r},{\bf r}^{\prime},t-t^{\prime})=(2\pi)^{-1}\int{\rm
d}\varepsilon\;{\rm e}^{{\rm
i}\varepsilon(t^{\prime}-t)}\,\mathcal{G}_{0}({\bf r},{\bf
r}^{\prime},\varepsilon)$. Likewise, surface translation symmetry reduces the
dependence on in-plane coordinates to the difference ${\bf R}-{\bf
R}^{\prime}$. Moving to Fourier space, this dependence can be expressed as a
combination of ${\rm e}^{{\rm i}{\bf q}_{\parallel}\cdot({\bf R}-{\bf
R}^{\prime})}$ waves with in-plane wave vectors ${\bf q}_{\parallel}$. In
addition, the Hamiltonian $\hat{\mathcal{H}}_{0}({\bf r})$ depends on ${\bf
R}$ only through the $\nabla^{2}_{\bf R}$ term, so we can separate the
frequency in parallel and perpendicular components as $\varepsilon=\hbar
q_{\parallel}^{2}/2m_{\rm e}+\varepsilon^{\perp}$. Putting these elements
together, we can write

$\displaystyle\mathcal{G}_{0}({\bf r},{\bf
r}^{\prime},t-t^{\prime})=\int\frac{{\rm d}^{2}{\bf
q}_{\parallel}}{(2\pi)^{2}}\;{\rm e}^{{\rm i}{\bf q}_{\parallel}\cdot({\bf
R}-{\bf R}^{\prime})}\int\frac{{\rm d}\varepsilon}{2\pi}\;{\rm e}^{{\rm
i}\varepsilon(t^{\prime}-t)}\;\mathcal{G}_{0}\big{(}z,z^{\prime},\varepsilon-\hbar
q_{\parallel}^{2}/2m_{\rm e}\big{)},$ (A11)

which, upon insertion into Eq. (A10), leads to the equation

$\displaystyle\bigg{[}-\frac{\hbar^{2}}{2m_{\rm
e}}\partial_{z}^{2}+V(z)-\hbar\varepsilon^{\perp}\bigg{]}\,\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp})=-\delta(z-z^{\prime})$
(A12)

for the 1D Green function $\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp})$.

We consider three different types of potentials, for which the Green function
can be calculated analytically as follows:

  * •

Finite potential barrier. To describe atomically thin two-dimensional (2D)
materials, we set $V(z)=U_{0}\delta(z)$. The parameter $U_{0}$ has units of
energy times length. Arguing that an atomic monolayer can be described by a
barrier of finite thickness $d\lesssim 1$ nm and internal potential $V_{0}$ in
the eV range, we expect $U_{0}\approx V_{0}d$ in the eV$\times$nm range when
adopting the zero-thickness approximation. The corresponding Green function
reads

$\displaystyle\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp})$
$\displaystyle=-\frac{{\rm i}m_{\rm e}}{\hbar^{2}q_{z}}\,\Big{[}{\rm e}^{{\rm
i}q_{z}|z-z^{\prime}|}+r_{q_{z}}\,{\rm e}^{{\rm
i}q_{z}(|z|+|z^{\prime}|)}\Big{]},$ (A13a)

where $q_{z}=\sqrt{2m_{\rm e}\varepsilon^{\perp}/\hbar}$ is the out-of-plane
wave vector, the first term accounts for free electron propagation, and the
second term is proportional to the energy-dependent reflection coefficient of
the potential barrier $r_{q_{z}}=\big{(}{\rm i}\hbar^{2}|q_{z}|/m_{\rm
e}U_{0}-1\big{)}^{-1}$. Upon inspection, Eq. (A13a) can be readily verified to
satisfy Eq. (A12) for the chosen potential.

  * •

Infinite potential barrier. As an instructive limit, we also consider complete
surface electron reflection, which can be realized through the potential
barrier introduced above by setting $U_{0}=\infty$ (i.e., with
$r_{q_{z}}=-1$). From Eq. (A13a), the Green function now reduces to

$\displaystyle\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp})$
$\displaystyle=-\frac{{\rm i}m_{\rm e}}{\hbar^{2}q_{z}}\big{(}{\rm e}^{{\rm
i}q_{z}|z-z^{\prime}|}-\,{\rm e}^{{\rm i}q_{z}(|z|+|z^{\prime}|)}\big{)}$
(A13b) $\displaystyle=-\frac{{\rm i}m_{\rm e}}{\hbar^{2}q_{z}}\left({\rm
e}^{{\rm i}q_{z}|z-z^{\prime}|}-\,{\rm e}^{{\rm
i}q_{z}(|z|+|z^{\prime}|)}\right)\,\big{[}\Theta(z)\,\Theta(z^{\prime})+\Theta(-z)\,\Theta(-z^{\prime})\big{]},$

where the second line shows that positive and negative $z$ regions do not mix
[i.e., an electron coming from $z>0$ stays fully inside that region after
scattering, as one can verify from Eq. (A22) below].

  * •

Free electron. When the interaction with the material can be neglected, we set
$V(z)=0$ (i.e., $U_{0}=0$), so we have $r_{q_{z}}=0$, and therefore, Eq.
(A13a) reduces to

$\displaystyle\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp})$
$\displaystyle=-\frac{{\rm i}m_{\rm e}}{\hbar^{2}q_{z}}{\rm e}^{{\rm
i}q_{z}|z-z^{\prime}|}.$ (A13c)

In addition, the initial wave function $\psi_{0}(z)$ produced by an incident
wave ${\rm e}^{-{\rm i}q_{0z}}$ also depends on the type of potential and
takes the explicit form

$\displaystyle\psi_{0}(z)=\left\\{\begin{array}[]{ll}\big{[}{\rm e}^{-{\rm
i}q_{0z}z}+r_{q_{0z}}\,{\rm e}^{{\rm
i}q_{0z}z}\big{]}\,\Theta(z)+t_{q_{0z}}\,{\rm e}^{-{\rm
i}q_{0z}z}\,\Theta(-z),&\quad\quad\quad\text{(finite potential barrier)}\\\
\big{[}{\rm e}^{-{\rm i}q_{0z}}-{\rm e}^{{\rm
i}q_{0z}}\big{]}\,\Theta(z),&\quad\quad\quad\text{(infinite potential
barrier)}\\\ {\rm e}^{-{\rm i}q_{0z}},&\quad\quad\quad\text{(free
electron)}\end{array}\right.$ (A17)

where $t_{q_{z}}=1+r_{q_{z}}$ is the transmission coefficient (evaluated at
the incident out-of-plane wave vector $q_{z}=-q_{0z}$), which guarantees the
continuity of $\psi_{0}(z)$. In the analysis presented below, we use the 1D
Green function and initial wave function for the finite potential barrier,
from which the free-electron and infinite-potential-barrier configurations are
directly obtained by setting $r_{q_{z}}=0$ and $r_{q_{z}}=-1$, respectively.

###  A1.2 Interaction of light

The optical field is introduced through a classical vector potential ${\bf
A}({\bf r},t)={\bf A}(z)\,{\rm e}^{{\rm i}{\bf k}_{\parallel}\cdot{\bf R}-{\rm
i}\omega t}+{\rm c.c.}$ with well-defined in-plane wave vector ${\bf
k}_{\parallel}$ and frequency $\omega$. By working in a gauge with vanishing
scalar potential, adopting the minimal-coupling prescription, and neglecting
$A^{2}$ terms, the interaction Hamiltonian reduces to

$\displaystyle\hat{\mathcal{H}}_{1}({\bf r},t)=\frac{-{\rm i}e\hbar}{m_{\rm
e}c}{\bf A}({\bf r},t)\cdot\bm{\nabla},$ (A18)

so the total Hamiltonian becomes $\hat{\mathcal{H}}({\bf
r},t)=\hat{\mathcal{H}}_{0}({\bf r})+\hat{\mathcal{H}}_{1}({\bf r},t)$. The
electron wave function $\psi({\bf r},t)$ is then modified relative to
$\psi_{0}({\bf r},t)$ as described by the Lippmann–Schwinger (LS) equation
Messiah (1966); Sakurai (1994) $\psi({\bf r},t)=\psi_{0}({\bf r},t)+\int{\rm
d}^{3}{\bf r}^{\prime}{\rm d}t^{\prime}\;\mathcal{G}_{0}({\bf r},{\bf
r}^{\prime},t-t^{\prime})\,\hat{\mathcal{H}}_{1}({\bf
r}^{\prime},t^{\prime})\,\psi({\bf r}^{\prime},t^{\prime})$. At this point, it
is convenient to write the perturbation series

$\displaystyle\psi({\bf r},t)=\sum_{n=0}^{\infty}\psi^{(n)}({\bf r},t),$ (A19)

where the $n$-th term is order $n$ in the vector potential ${\bf A}$. By
inserting Eq. (A19) into the LS equation, we find the recursion formula

$\displaystyle\psi^{(n)}({\bf r},t)=\int{\rm d}^{3}{\bf r}^{\prime}\int{\rm
d}t^{\prime}\,\mathcal{G}_{0}({\bf r},{\bf
r}^{\prime},t-t^{\prime})\,\hat{\mathcal{H}}_{1}({\bf
r}^{\prime},t^{\prime})\,\psi^{(n-1)}({\bf r}^{\prime},t^{\prime})$ (A20)

for $n>0$, starting with $\psi^{(0)}\equiv\psi_{0}$.

Because the vector potential depends on time through ${\rm e}^{\pm{\rm
i}\omega t}$, each subsequent perturbation order introduces the possibility of
absorbing or emitting one additional photon. Hence, at a given order $n$ the
electron energy can be modified to $\hbar(\varepsilon_{0}+\ell\omega)$ with
$|\ell|\leq n$. Using this and further imposing conservation of in-plane
momentum, we can decompose the post-interaction wave function as

$\displaystyle\psi({\bf
r},t)=\sum_{n=0}^{\infty}\sum_{\ell=-n}^{n}\varphi_{n\ell}(z)\,{\rm e}^{{\rm
i}({\bf q}_{0\parallel}+\ell\,{\bf k}_{\parallel})\cdot{\bf R}}\,{\rm
e}^{-{\rm i}(\varepsilon_{0}+\ell\omega)t}.$ (A21)

We now introduce Eqs. (A11) and (A21) into Eq. (A20), use the factorization in
Eq. (A9), and identify terms with the same $({\bf R},t)$ dependence on both
sides of the equation. This procedure shows that Eq. (A21) is indeed a
solution when the $z$-dependent coefficients are determined by iterating the
relation

$\displaystyle\varphi_{n\ell}(z)=\frac{\hbar e}{m_{\rm e}c}\int{\rm
d}{z^{\prime}}\,\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp}_{\ell})\;\Big{\\{}{\bf
A}(z^{\prime})$ $\displaystyle\cdot\big{[}{\bf q}_{0\parallel}+(\ell-1)\,{\bf
k}_{\parallel}-{\rm i}\hat{\bf
z}\,\partial_{z^{\prime}}\big{]}\,\varphi_{n-1,\ell-1}(z^{\prime})$
$\displaystyle+{\bf A}^{*}(z^{\prime})$ $\displaystyle\cdot\big{[}{\bf
q}_{0\parallel}+(\ell+1)\,{\bf k}_{\parallel}-{\rm i}\hat{\bf
z}\,\partial_{z^{\prime}}\big{]}\,\varphi_{n-1,\ell+1}(z^{\prime})\Big{\\}},$
(A22)

where $\varepsilon_{\ell}^{\perp}=\varepsilon_{0}+\ell\omega-\hbar|{\bf
q}_{\parallel}+\ell\,{\bf k}_{\parallel}|^{2}/2m_{\rm e}$. Equation (A22),
which is reproduced as Eq. (2) in the main text, applies to $n>0$ with a
seeding term $\varphi_{00}(z)=\psi_{0}(z)$ and imposing $\varphi_{n\ell}(z)=0$
for $|\ell|>n$.

For simplicity, we assume the in-plane wave vectors of light and electron to
be both along the surface direction $x$ (i.e., ${\bf
q}_{0\parallel}=q_{0\parallel}\hat{\bf x}$ and ${\bf
k}_{\parallel}=k_{\parallel}\hat{\bf x}$). In what follows, we consider two
types of optical fields:

  * •

Unscattered light plane waves. For sufficiently thin films (e.g., self-
standing atomic monolayers), light scattering by the material can be
approximately neglected. Then, taking an incident wave vector ${\bf
k}=k_{\parallel}\hat{\bf x}-k_{z}\hat{\bf z}$ and an electric field amplitude
${\bf E}_{0}$, the $z$ dependence of the vector potential reduces to

$\displaystyle{\bf A}(z)=\frac{1}{{\rm i}k}{\bf E}_{0}\,{\rm e}^{{\rm
i}k_{z}z}.$ (A23a)

We consider propagating light with $k_{\parallel}\leq k=\omega/c$ and
$k_{z}=\sqrt{k^{2}-k_{\parallel}^{2}}$.

  * •

Surface polaritons. For films of negligible thickness, the electric field
associated with a polariton of in-plane wave vector ${\bf
k}_{\parallel}=k_{\parallel}\hat{\bf x}$ and frequency $\omega$ can be written
as ${\bf E}({\bf r},t)={\bf E}({\bf r})\,{\rm e}^{-{\rm i}\omega t}+{\rm
c.c.}$ with ${\bf E}({\bf
r})=E_{0}\;\big{(}k_{\parallel}^{2}+\kappa^{2}\big{)}^{-1/2}\,\big{(}{\rm
i}\,\kappa\,\hat{\bf x}-\operatorname{sign}\\{z\\}\,k_{\parallel}\,\hat{\bf
z}\,\big{)}\,{\rm e}^{{\rm i}k_{\parallel}x-\kappa|z|}$, where $E_{0}$ is a
global amplitude and $\kappa=\sqrt{k_{\parallel}^{2}-k^{2}}$ determines an
out-of-plane exponential decay with $k_{\parallel}>k$. This expression is
constructed to satisfy $\nabla\cdot{\bf E}({\bf r},t)=0$ at $z\neq 0$. Then,
the vector potential reduces to ${\bf A}({\bf r})={\bf E}({\bf r})/{\rm
i}k={\bf A}(z)\,{\rm e}^{{\rm i}k_{\parallel}x}$, where the $z$ dependence is
given by

$\displaystyle{\bf
A}(z)=\frac{E_{0}}{k}\frac{1}{\sqrt{\kappa^{2}+k_{\parallel}^{2}}}\;\Big{(}\kappa\hat{\bf
x}+{\rm i}\,k_{\parallel}\operatorname{sign}\\{z\\}\,\hat{\bf z}\Big{)}\,{\rm
e}^{-\kappa|z|}.$ (A23b)

We find it convenient to express the in-plane wave vector as
$k_{\parallel}=n_{\rm eff}k$ in terms of an effective refractive index $n_{\rm
eff}$. In this work, we take $n_{\rm eff}\gg 1$ values, so have $\kappa\approx
k_{\parallel}$, and therefore, ${\bf
A}(z)\approx\big{(}E_{0}/\sqrt{2}\,k\big{)}\;\big{(}\hat{\bf x}+{\rm
i}\,\operatorname{sign}\\{z\\}\,\hat{\bf z}\big{)}\,{\rm
e}^{-k_{\parallel}|z|}$. However, we use the fully retarded expressions in our
calculations. Incidentally, these expressions can also be used for arbitrary
planar material systems contained within the $z<0$ region if we are only
interested in the outer $z>0$ half-space (i.e., without electron penetration
inside the material).

Next, we solve Eq. (A22) by plugging the expressions for
$\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp}_{\ell})$ and ${\bf A}(z)$
given in Eqs. (A13) and (A23).

###  A1.3 Recursive solution for the self-consistent electron wave function

We solve Eq. (A22) by expanding the wave function components
$\varphi_{n\ell}(z)$ following a recursive procedure to determine their
explicit dependence on $z$, the parameters involved in such dependence, and
how they are modified at each iteration step given the analytical expressions
for $\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp}_{\ell})$ and ${\bf
A}(z)$ in Eqs. (A13) and (A23). Because these quantities depend on $z$ through
exponential factors, and therefore, the $z^{\prime}$ integral at each
iteration $n$ in Eq. (A22) produces additional exponential factors consisting
of a combination of those at order $n-1$ plus those of
$\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp}_{\ell})$ and ${\bf A}(z)$,
the wave function components must take the general form

$\displaystyle\varphi_{n\ell}(z)=\sum_{s=\pm
1}\sum_{j=1}^{N_{n\ell}}\alpha_{n\ell}^{js}\,{\rm
e}^{\zeta_{n\ell}^{js}z}\,\Theta{(sz)}$ (A24)

[Eq. (3) in the main text], where the number of terms in this sum
($N_{n\ell}$) is obviously increasing with $n$ and also depends on $\ell$.
Specifically, the vector potentials under consideration can be written as

$\displaystyle{\bf A}(z)=\sum_{s=\pm 1}{\bf A}^{s}\,{\rm
e}^{\eta^{s}z}\,\Theta{(sz)},$ (A25)

where the coefficients ${\bf A}^{\pm 1}$ and $\eta^{\pm 1}$ can be directly
identified from Eqs. (A23). Likewise, the Green function
$\mathcal{G}_{0}(z,z^{\prime},\varepsilon^{\perp}_{\ell})$ is given by Eq.
(A13a) with the reflection coefficient $r_{q_{\ell z}}$ set to different
values depending on the choice for the potential $V(z)$. Note that in the
$\ell$ channel, the out-of-plane wave vector is given by $q_{\ell
z}=\sqrt{2m_{\rm e}\varepsilon^{\perp}_{\ell}/\hbar+{\rm
i}0^{+}}=\sqrt{2m_{\rm e}(\varepsilon_{0}+\ell\omega)/\hbar-|{\bf
q}_{\parallel}+\ell\,{\bf k}_{\parallel}|^{2}+{\rm i}0^{+}}$, where the square
root is taken to yield a positive imaginary part. By inserting Eqs. (A13a),
(A24), and (A25) into Eq. (A22), we obtain the equivalent recursion formula
(for $\pm z>0$)

$\displaystyle\sum_{j}\alpha_{n\ell}^{j\pm}\,\boxed{{\rm
e}^{\zeta_{n\ell}^{j\pm}z}}=$ $\displaystyle-\frac{{\rm i}e}{\hbar cq_{\ell
z}}\sum_{s=\pm 1}\sum_{j}\Bigg{\\{}2{\rm i}q_{\ell z}\delta_{s,\pm
1}\Bigg{[}\boxed{{\rm e}^{(\eta^{s}+\zeta_{n-1,\ell-1}^{js})z}}\;{\bf
A}^{s}\cdot\big{[}{\bf q}_{0\parallel}+(\ell-1)\,{\bf k}_{\parallel}-{\rm
i}\hat{\bf
z}\,\zeta_{n-1,\ell-1}^{js}\big{]}\,\frac{\alpha_{n-1,\ell-1}^{js}}{q_{\ell
z}^{2}+\big{(}\eta^{s}+\zeta_{n-1,\ell-1}^{js}\big{)}^{2}}$
$\displaystyle\quad\quad\quad+\boxed{{\rm
e}^{(\eta^{s*}+\zeta_{n-1,\ell+1}^{js})z}}\;{\bf A}^{s*}\cdot\big{[}{\bf
q}_{0\parallel}+(\ell+1)\,{\bf k}_{\parallel}-{\rm i}\hat{\bf
z}\,\zeta_{n-1,\ell+1}^{js}\big{]}\,\frac{\alpha_{n-1,\ell+1}^{js}}{q_{\ell
z}^{2}+\big{(}\eta^{s*}+\zeta_{n-1,\ell+1}^{js}\big{)}^{2}}\Bigg{]}$ (A26)
$\displaystyle+\boxed{{\rm e}^{\pm{\rm i}q_{\ell z}z}}\;\Bigg{[}$
$\displaystyle{\bf A}^{s}\cdot\big{[}{\bf q}_{0\parallel}+(\ell-1)\,{\bf
k}_{\parallel}-{\rm i}\hat{\bf
z}\,\zeta_{n-1,\ell-1}^{js}\big{]}\,\alpha_{n-1,\ell-1}^{js}\;\Bigg{(}\frac{-r_{q_{\ell
z}}}{{\rm i}q_{\ell
z}+s\;\big{(}\eta^{s}+\zeta_{n-1,\ell-1}^{js}\big{)}}\pm\frac{s}{{\rm
i}q_{\ell z}\mp\big{(}\eta^{s}+\zeta_{n-1,\ell-1}^{js}\big{)}}\Bigg{)}$
$\displaystyle+$ $\displaystyle{\bf A}^{s*}\cdot\big{[}{\bf
q}_{0\parallel}+(\ell+1)\,{\bf k}_{\parallel}-{\rm i}\hat{\bf
z}\,\zeta_{n-1,\ell+1}^{js}\big{]}\,\alpha_{n-1,\ell+1}^{js}\;\Bigg{(}\frac{-r_{q_{\ell
z}}}{{\rm i}q_{\ell
z}+s\;\big{(}\eta^{s*}+\zeta_{n-1,\ell+1}^{js}\big{)}}\pm\frac{s}{{\rm
i}q_{\ell
z}\mp\big{(}\eta^{s*}+\zeta_{n-1,\ell+1}^{js}\big{)}}\Bigg{)}\Bigg{]}\Bigg{\\}},$

where we have performed the $z^{\prime}$ integrals using the identities

$\displaystyle\int{\rm d}{z^{\prime}}\,\Theta{(\pm z^{\prime})}\,{\rm e}^{{\rm
i}q_{\ell z}|z^{\prime}|+\Delta z^{\prime}}=\frac{-1}{{\rm i}q_{\ell
z}\pm\Delta},$ (A27a) $\displaystyle\int{\rm d}{z^{\prime}}\,\Theta{(\pm
z^{\prime})}\,{\rm e}^{{\rm i}q_{\ell z}|z-z^{\prime}|+\Delta
z^{\prime}}=\frac{\pm\operatorname{sign}\\{z\\}}{{\rm i}q_{\ell
z}-\operatorname{sign}\\{z\\}\Delta}\,{\rm e}^{{\rm i}q_{\ell
z}|z|}+\frac{2{\rm i}q_{\ell z}}{q_{\ell z}^{2}+\Delta^{2}}\,{\rm e}^{\Delta
z}\;\Theta(\pm z).$ (A27b)

Incidentally, the exponentials inside the integrals vanish at
$|z^{\prime}|\to\infty$ because they contain either evanescent electron/light
components or propagating fields in which the electron wave vectors $q_{\ell
z}$ possess an arbitrarily small positive imaginary part that makes them
consistent with the use of retarded Green functions. At every iteration step
$n>0$, we solve Eq. (A26) by setting the coefficients $\alpha_{n\ell}^{j\pm}$
and $\zeta_{n\ell}^{j\pm}$ in the left-hand side to match the exponential
terms in the right-hand side (see boxed expressions), starting with the $n=0$
seeding values $\big{\\{}\alpha_{00}^{1,+1}=1,\,\zeta_{00}^{1,+1}=-{\rm
i}q_{0z}\big{\\}}$,
$\big{\\{}\alpha_{00}^{2,+1}=r_{q_{0z}},\,\zeta_{00}^{2,+1}={\rm
i}q_{0z}\big{\\}}$, and
$\big{\\{}\alpha_{00}^{1,-1}=t_{q_{0z}},\,\zeta_{00}^{1,-1}=-{\rm
i}q_{0z}\big{\\}}$, as determined from Eq. (A17), and noticing that
$|\ell|\leq n$.

When considering polaritons, for which $\eta^{s}$ in Eq. (A25) has a finite
real part that makes the optical field evanescent, upon inspection of the
numerical solution of Eq. (A26), we corroborate that the obtained exponential
coefficients of all propagating components in Eq. (A26) (i.e., with vanishing
${\rm Re}\\{\zeta_{n\ell}^{js}\\}$) satisfy $\zeta_{n\ell}^{js}={\rm
i}\,sq_{\ell z}$, as expected from energy conservation after a net exchange of
$\ell$ photons. The remaining coefficients are found to possess finite real
parts satisfying the condition $s\,{\rm Re}\\{\zeta_{n\ell}^{js}\\}<0$ (i.e.,
they are evanescent).

However, under external illumination, Eq. (A25) comprises propagating
components (i.e., imaginary coefficients $\eta^{s}$), so the interaction
region extends indefinitely in $z$, affecting the electron in the way
described by Volkov Wolkow (1935). A direct solution of Eq. (A26) is then
problematic because one needs to separate surface-mediated inelastic
transitions from far-field electron–light interactions. In actual experiments,
the range of the applied light fields is limited to a finite region near the
surface (e.g., when using a laser beam), so we find it more convenient to
introduce a finite real part in the exponential coefficients as ${\rm
Re}\\{\eta^{s}\\}=-s\gamma$ with $\gamma>0$ to prevent electron exposure to
far optical fields. Then, we iterate Eq. (A26) to the desired order $n$,
perform the $|z|\to\infty$ limit to calculate the scattered electron
intensity, and finally take $\gamma\to 0^{+}$. In practice, this amounts to
calculating the electron intensity after eliminating the terms of Eq. (A25) in
which $\zeta_{n\ell}^{js}$ depends explicitly on $\gamma$. The remaining far-
field terms also satisfy $\zeta_{n\ell}^{js}={\rm i}\,sq_{\ell z}$.

##  Appendix A2 Low-energy electron scattering by an illuminated atomic
monolayer including diffraction

While in Sec. A we discussed homogenous surfaces, the in-plane atomic
corrugation needs to be incorporated in a more realistic scenario. Next, we
extend the theory to incorporate diffraction by an illuminated monoatomic
layer.

###  B2.1 Low-energy electron diffraction by an atomic monolayer

Low-energy electron diffraction (LEED) is a well-established technique for the
determination of surface atomic structures based on the measurement and
comparison with the theory of the electron-energy dependence of the
intensities associated with different Bragg reflections by periodic atomic
lattices Pendry (1974); Van Hove  _et al._ (1986). Here, we borrow the
theoretical methods developed to simulate LEED and describe electron
diffraction by an atomic monolayer Pendry (1974). Following the notation of
Sec. A, we consider an electron beam of energy $\hbar\varepsilon_{0}$
impinging on a crystalline atomic monolayer with wave vector ${\bf q}_{0}={\bf
q}_{0\parallel}-q_{0z}\,\hat{\bf z}$, where
$q_{0z}=\sqrt{q_{0}^{2}-q_{0\parallel}^{2}}$. For simplicity, we assume a
simple lattice of the same type of atoms, periodically arranged at positions
${\bf R}_{\alpha}$ in the $z=0$ plane and with unit cell area $A$. The
Hamiltonian of the system can be approximated as

$\displaystyle\hat{\mathcal{H}}_{0}({\bf
r})=-\frac{\hbar^{2}\nabla^{2}}{2m_{\rm e}}+\sum_{\alpha}V_{\rm atom}(|{\bf
r}-{\bf R}_{\alpha}|),$

where $V_{\rm atom}(r)$ is the atomic potential, considered to be isotropic.
The electron is diffracted by the lattice and exchanges discrete amounts of
in-plane momentum $\hbar{\bf G}$ determined by the 2D reciprocal lattice
vectors ${\bf G}$, resulting in scattered beams with three-dimensional wave
vectors ${\bf q}_{0\parallel}+{\bf G}\pm q_{{\bf G}\varepsilon_{0}z}\hat{\bf
z}$ directed along upward ($+$) and downward ($-$) directions. The out-of-
plane wave vector component $q_{{\bf G}\varepsilon_{0}z}=\sqrt{q_{0}^{2}-|{\bf
q}_{0\parallel}+{\bf G}|^{2}+{\rm i}0^{+}}$ is determined by energy
conservation with the square root taken to yield a positive imaginary part.
For $|{\bf q}_{0\parallel}+{\bf G}|<q_{0}$, the diffracted waves propagate to
the far field, giving rise to observable LEED spots, whereas waves with $|{\bf
q}_{0\parallel}+{\bf G}|>q_{0}$ are evanescent. The total electron wave
function outside the layer should take the form

$\displaystyle\psi_{0}({\bf r},t)=\bigg{[}{\rm e}^{-{\rm i}q_{0z}z}+\sum_{{\bf
G},\pm}B^{\pm}_{{\bf G}}\;{\rm e}^{{\rm i}({\bf G}\cdot{\bf R}\pm q_{{\bf
G}\varepsilon_{0}z}z)}\;\Theta{(\pm z)}\bigg{]}\;{\rm e}^{{\rm i}({\bf
q}_{0\parallel}\cdot{\bf R}-\varepsilon_{0}t)},$ (B28)

where $B^{\pm}_{{\bf G}}$ are Bragg scattering coefficients.

To determine $B^{\pm}_{{\bf G}}$ in Eq. (B28), we start by expressing the
incident wave around each atomic position ${\bf R}_{\alpha}$ Pendry (1974) as

$\displaystyle{\rm e}^{{\rm i}{\bf q}_{0}\cdot{\bf r}}={\rm e}^{{\rm i}{\bf
q}_{0\parallel}\cdot{\bf
R}_{\alpha}}\sum_{L}\phi^{\mathrm{inc}}_{L}j_{L}\big{[}q_{0}({\bf r}-{\bf
R}_{\alpha})\big{]},$ (B29a) where
$\displaystyle\phi^{\mathrm{inc}}_{L}=4\pi\,Y_{L}^{*}\big{(}\Omega_{{\bf
q}_{0}}\big{)}$ (B29b)

are expansion coefficients involving spherical harmonics $Y_{L}$ labeled by
the quantum numbers $L=(l,m)$, while $j_{L}(q_{0}{\bf r})={\rm
i}^{l}j_{l}(q_{0}r)Y_{L}\big{(}\Omega_{\hat{\bf r}}\big{)}$ are regular
spherical waves with a radial dependence introduced through spherical Bessel
functions $j_{l}$. Note that we separate the phase associated with the
propagation of the incident wave to each atomic site as a global factor ${\rm
e}^{{\rm i}{\bf q}_{0\parallel}\cdot{\bf R}_{\alpha}}$ in Eq. (B29a).
Incidentally, under normal incidence (${\bf q}_{0}=-\hat{\bf z}$) we have
$\phi^{\mathrm{inc}}_{L}=\delta_{m0}\,(-1)^{l}\sqrt{4\pi\,(2l+1)}$. In the far
field ($q_{0}r\gg 1$), regular waves have the asymptotic behavior

$j_{l}(q_{0}r)=\frac{1}{2}\Big{[}h^{(1)}_{l}(q_{0}r)+h^{(2)}_{l}(q_{0}r)\Big{]}\approx\frac{1}{2q_{0}r}\Big{[}{\rm
i}^{-(l+1)}{\rm e}^{{\rm i}q_{0}r}+{\rm i}^{l+1}{\rm e}^{-{\rm
i}q_{0}r}\Big{]},$

which corresponds to a combination of outgoing and incoming waves expressed in
terms of the Hankel functions $h^{(1)}_{l}$ and $h^{(2)}_{l}$, respectively.
Scattering by an atom only affects the outgoing wave, which retains its
magnitude and quantum numbers $L$ because $V_{\rm atom}(r)$ is a real and
isotropic potential. Therefore, the waves become

$j_{l}(q_{0}r)\rightarrow\frac{1}{2q_{0}r}\Big{[}e^{2i\delta_{l}}h^{(1)}_{l}(q_{0}r)+h^{(2)}_{l}(q_{0}r)\Big{]},$

where $\delta_{l}$ are phase shifts that depend on the orbital quantum number
$l$ and the atomic potential (i.e., the type of atom). Here, we calculate
$\delta_{l}$ numerically assuming free atoms, for which the potential is
obtained by using a Dirac-Fock code Ankudinov  _et al._ (1996), while electron
scattering is computed by solving the Schrödinger equation following a well-
established procedure Salvat and Mayol (1991). For an individual atom, the
scattered components can thus be written as García de Abajo  _et al._ (2001)
$\sum_{L}\phi^{\mathrm{scat}}_{L}h^{(+)}_{L}\big{[}(q_{0}({\bf r}-{\bf
R}_{\alpha})\big{]}$, where $h^{(+)}_{L}(q_{0}{\bf r})={\rm
i}^{l+1}h^{(1)}_{l}(q_{0}r)Y_{L}\big{(}\Omega_{\hat{\bf r}}\big{)}$ are
spherical outgoing waves and
$\phi^{\mathrm{scat}}_{L}=t_{l}\,\phi^{\mathrm{inc}}_{L}$ are scattering
coefficients that involve the scattering matrix elements $t_{l}={\rm e}^{{\rm
i}\delta_{l}}\sin\delta_{l}$.

In the atomic monolayer, lattice symmetry allows us to write

$\displaystyle\psi_{0}({\bf r},t)=\bigg{[}{\rm e}^{{\rm i}{\bf q}_{0}\cdot{\bf
r}}+\sum_{\alpha}{\rm e}^{{\rm i}{\bf q}_{0\parallel}\cdot{\bf
R}_{\alpha}}\sum_{L}\phi^{\mathrm{scat}}_{L}\,h^{(+)}_{L}\big{[}(q_{0}({\bf
r}-{\bf R}_{\alpha})\bigg{]}\;{\rm e}^{-{\rm i}\varepsilon_{0}t},$ (B30)

where the propagation phase factor ${\rm e}^{{\rm i}{\bf
q}_{0\parallel}\cdot{\bf R}_{\alpha}}$ is retained in the contribution of
scattering by each lattice site $\alpha$. Now, we incorporate multiple
scattering in the monolayer by supplementing the incident wave at each atom
with the result of scattering from the rest of the atoms. More precisely, we
expand the waves outgoing from each position $\beta$ as a combination of
regular waves at each other position $\alpha$ using the identity Pendry (1974)
$h^{(+)}_{L^{\prime}}\big{[}(q_{0}({\bf r}-{\bf
R}_{\beta})\big{]}=\sum_{L}G_{\alpha\beta,LL^{\prime}}\,j_{L}\big{[}(q_{0}({\bf
r}-{\bf R}_{\alpha})\big{]}$, where
$G_{\alpha\beta,LL^{\prime}}=4\pi\sum_{L^{\prime\prime}}h^{(+)}_{L^{\prime\prime}}\big{[}(q_{0}({\bf
R}_{\alpha}-{\bf R}_{\beta})\big{]}\int{\rm
d}\Omega\;Y_{L}(\Omega)Y_{L^{\prime\prime}}(\Omega)Y_{L^{\prime}}^{*}(\Omega)$.
From these considerations, we find the self-consistent equation
$\phi^{\text{scat}}_{L}=t_{l}\,\big{[}\phi^{\text{inc}}_{L}+\sum_{L^{\prime}}\,G_{LL^{\prime}}({\bf
q}_{0\parallel})\,\phi^{\text{scat}}_{L^{\prime}}\big{]}$, or equivalently,

$\displaystyle\phi^{\text{scat}}_{L}=\sum_{L^{\prime}}\big{[}S^{-1}\big{]}_{LL^{\prime}}t_{l^{\prime}}\,\phi^{\text{inc}}_{L^{\prime}},$
(B31)

where we define the lattice sum $G_{LL^{\prime}}({\bf
q}_{0\parallel})=\sum_{\beta\neq\alpha}\,{\rm e}^{{\rm i}{\bf
q}_{0\parallel}\cdot({\bf R}_{\beta}-{\bf
R}_{\alpha})}\,G_{\alpha\beta,LL^{\prime}}$, which is trivially independent of
$\alpha$, and the matrix $S_{LL^{\prime}}=t_{l}\,G_{LL^{\prime}}({\bf
q}_{0\parallel})$. In our simulations, we use the efficient numerical methods
developed by Kambe Kambe (1967); Pendry (1974) to carry out the sum over
$\beta$ and find converged results by truncating the multipolar expansion to
$l\leq 5$.

To bring Eq. (B30) into the form of Eq. (B28), we use the identity (Eq. (A.39)
in Ref. Pendry (1974))

$h_{L}^{(+)}(q_{0}{\bf r})=-\frac{1}{2\pi^{2}q_{0}}\int{\rm d}^{3}{\bf
q}\;Y_{L}\big{(}\Omega_{\bf q}\big{)}\,\frac{{\rm e}^{{\rm i}{\bf q}\cdot{\bf
r}}}{q_{0}^{2}-q^{2}+{\rm i}0^{+}}$

and then convert the sum over $\alpha$ into a sum over reciprocal lattice
vectors as

$\sum_{\alpha}{\rm e}^{{\rm i}({\bf q}_{0\parallel}-{\bf
q}_{\parallel})\cdot{\bf R}_{\alpha}}=\frac{(2\pi)^{2}}{A}\sum_{\bf
G}\delta({\bf q}_{0\parallel}+{\bf G}-{\bf q}_{\parallel}).$

The ${\bf q}_{\parallel}$ integral is done using these $\delta$-functions,
while the integral over $q_{z}$ can be performed by closing the integration
contour in the upper (lower) complex $q_{z}$ plane for $z>0$ ($z<0$). Finally,
we recover an expression like Eq. (B28), from which we identify

$\displaystyle B^{\pm}_{\bf G}$ $\displaystyle=\frac{2\pi{\rm
i}}{A\,q_{0}q_{{\bf G}\varepsilon_{0}z}}\sum_{L}Y_{L}\big{(}\Omega_{{\bf
q}_{0\parallel}+{\bf G}\pm q_{{\bf G}\varepsilon_{0}z}\hat{\bf
z}}\big{)}\,\phi_{L}^{\mathrm{scat}}=\frac{8\pi^{2}{\rm i}}{A\,q_{0}q_{{\bf
G}\varepsilon_{0}z}}\sum_{LL^{\prime}}Y_{L}\big{(}\Omega_{{\bf
q}_{0\parallel}+{\bf G}\pm q_{{\bf G}\varepsilon_{0}z}\hat{\bf
z}}\big{)}\,\big{[}S^{-1}\big{]}_{LL^{\prime}}\,t_{l^{\prime}}\,Y_{L^{\prime}}^{*}\big{(}\Omega_{{\bf
q}_{0}}\big{)},$ (B32)

where the rightmost identity is obtained by inserting the scattering wave
coefficients $\phi_{L}^{\mathrm{scat}}$ obtained from Eqs. (B29b) and (B31).

From the above considerations, we can also obtain the Green function
$\mathcal{G}_{0}({\bf r},{\bf r}^{\prime},\varepsilon)$ corresponding to an
electron energy $\hbar\varepsilon$. In analogy to Eq. (A11), we separate it
into the contribution of different in-plane wave vector components ${\bf
q}_{\parallel}$, each of them connected to the whole set of in-plane wave
vectors $\big{\\{}{\bf q}_{\parallel}+{\bf G}\big{\\}}$ (differing by
reciprocal lattice vectors ${\bf G}$) via scattering by the atomic monolayer.
We thus restrict ${\bf q}_{\parallel}$ to the first Brillouin zone (1BZ) and
write the frequency-space Green function as

$\displaystyle\mathcal{G}_{0}({\bf r},{\bf
r}^{\prime},t-t^{\prime})=\sum_{{\bf G}{\bf G}^{\prime}}\int_{\rm
1BZ}\frac{{\rm d}^{2}{\bf q}_{\parallel}}{(2\pi)^{2}}\;{\rm e}^{{\rm i}({\bf
q}_{\parallel}+{\bf G})\cdot{\bf R}-{\rm i}({\bf q}_{\parallel}+{\bf
G}^{\prime})\cdot{\bf R}^{\prime}}\int\frac{{\rm d}\varepsilon}{2\pi}\;{\rm
e}^{{\rm i}\varepsilon(t^{\prime}-t)}\;\mathcal{G}_{{\bf G}{\bf
G}^{\prime}}\big{(}{\bf q}_{\parallel},z,z^{\prime},\varepsilon\big{)}.$
(B33a) Energy conservation in the regions outside the monolayer (i.e., above
and below) implies that $\mathcal{G}_{{\bf G}{\bf G}^{\prime}}\big{(}{\bf
q}_{\parallel},z,z^{\prime},\varepsilon\big{)}$ must be made of components
with a $(z,z^{\prime})$ dependence given by ${\rm e}^{\pm{\rm i}q_{{\bf
G}\varepsilon z}z\pm{\rm i}q_{{\bf G}^{\prime}\varepsilon z}z^{\prime}}$,
where $q_{{\bf G}\varepsilon z}=\sqrt{2m_{\rm e}\varepsilon/\hbar-|{\bf
q}_{\parallel}+{\bf G}|^{2}+{\rm i}0^{+}}$. Similarly to Eq. (A13a), we expect
the Green function to consist of a free-space term plus contributions arising
from wave scattering at the atomic layer (i.e., reflection and transmission),
that is, $\displaystyle\mathcal{G}_{{\bf G}{\bf G}^{\prime}}\big{(}{\bf
q}_{\parallel},z,z^{\prime},\varepsilon\big{)}=-\frac{{\rm i}m_{\rm
e}}{\hbar^{2}q_{{\bf G}\varepsilon z}}\,\Big{[}\delta_{{\bf G}{\bf
G}^{\prime}}{\rm e}^{{\rm i}q_{{\bf G}\varepsilon
z}|z-z^{\prime}|}+r^{ss^{\prime}}_{{\bf G}{\bf G}^{\prime}}({\bf
q}_{\parallel},\varepsilon)\,{\rm e}^{{\rm i}q_{{\bf G}\varepsilon z}|z|+{\rm
i}q_{{\bf G}^{\prime}\varepsilon z}|z^{\prime}|}\Big{]},$ (B33b) where
$r^{ss^{\prime}}_{{\bf G}{\bf G}^{\prime}}({\bf q}_{\parallel},\varepsilon)$
are Bragg scattering coefficients with $s=\operatorname{sign}\\{z\\}$ and
$s^{\prime}=\operatorname{sign}\\{z^{\prime}\\}$ denoting the direction of
wave propagation (i.e., upward and downward for $s,s^{\prime}=1$ and -1,
respectively). Following the same steps as in the derivation presented above
for Eqs. (B28) and (B32), but now applied to incident waves ${\rm e}^{{\rm
i}({\bf q}_{\parallel}+{\bf G}^{\prime})\cdot{\bf R}+{\rm i}s^{\prime}q_{{\bf
G}^{\prime}\varepsilon z}z}$ and outgoing waves ${\rm e}^{{\rm i}({\bf
q}_{\parallel}+{\bf G})\cdot{\bf R}+{\rm i}sq_{{\bf G}\varepsilon,z}z}$ (with
$s=\pm 1$ and $s^{\prime}=\pm 1$ indicating upward/downward propagation), we
can readily write $\displaystyle r^{ss^{\prime}}_{{\bf G}{\bf
G}^{\prime}}({\bf q}_{\parallel},\varepsilon)=\frac{8\pi^{2}{\rm
i}}{A\,q_{0}q_{{\bf G}\varepsilon
z}}\sum_{LL^{\prime}}Y_{L}\big{(}\Omega_{{\bf q}_{\parallel}+{\bf G}+sq_{{\bf
G}\varepsilon z}\hat{\bf
z}}\big{)}\,\big{[}S^{-1}\big{]}_{LL^{\prime}}\,t_{l^{\prime}}\,Y_{L^{\prime}}^{*}\big{(}\Omega_{{\bf
q}_{\parallel}+{\bf G}^{\prime}+s^{\prime}q_{{\bf G}^{\prime}\varepsilon
z}\hat{\bf z}}\big{)}.$ (B33c)

We note that Eqs. (B33) are only valid outside the range of the atomic
potentials $V_{\rm atom}(|{\bf r}-{\bf R}_{\alpha}|)$, which represents a thin
region and can therefore be safely neglected to describe the interaction with
light via the LS equation [Eq. (A20)].

###  B2.2 Interaction with light

The Green function in Eqs. (B33) accounts for electron–lattice interaction to
all orders of scattering. Inserting it into Eq. (A20) together with the
initial wave function in Eq. (B28), we can obtain a perturbative solution to
the wave function up to an arbitrary order $n$ in the interaction with light.
After a net exchange of $\ell$ photons (with $|\ell|\leq n$), the in-plane
wave vector and frequency of the electron become ${\bf
q}_{0\parallel}+\ell{\bf k}_{\parallel}$ and
$\varepsilon_{\ell}=\varepsilon_{0}+\ell\omega$, respectively, as deduced from
the LS equation by performing the ${\bf R}^{\prime}$ and $t^{\prime}$
integrals, which then select the $\mathcal{G}_{{\bf G}{\bf
G}^{\prime}}\big{(}{\bf q}_{0\parallel}+\ell{\bf
k}_{\parallel},z,z^{\prime},\varepsilon_{0}+\ell\omega\big{)}$ Green function
component.

A general analysis is involved and does not provide substantial insights.
Instead, we present results assuming normally incident electrons (${\bf
q}_{0\parallel}=0$ and ${\bf q}_{0}=-q_{0}\hat{\bf z}$) and limit our
discussion to first-order interaction with a grazingly incident light plane
wave of wave vector ${\bf k}_{\parallel}=k\,\hat{\bf x}$ and real electric
field amplitude ${\bf E}_{0}$ in the $y-z$ plane (i.e., ${\bf A}({\bf
r},t)=({\bf E}_{0}/{\rm i}k)[{\rm e}^{{\rm i}(kx-\omega t)}-{\rm e}^{-{\rm
i}(kx-\omega t)}]$). Then, the $\ell=0$ wave function component remains as in
Eq. (B28), while two sidebands $\ell=\pm 1$ emerge at order $n=1$. Inserting
the field into Eq. (A18) and combining the result with Eqs. (A20), (B28), and
(B33), we obtain

$\displaystyle\psi^{(1)}({\bf r},t)=\frac{e}{\hbar\omega}$
$\displaystyle\sum_{\ell=\pm 1}{\rm i}^{\ell-1}\;{\rm e}^{-{\rm
i}\varepsilon_{\ell}t}\sum_{{\bf G}{\bf G}^{\prime}}\frac{1}{q_{{\bf
G}\varepsilon_{\ell}z}}{\rm e}^{{\rm i}(\ell{\bf k}_{\parallel}+{\bf
G})\cdot{\bf R}}$ $\displaystyle\times{\bf E}_{0}\cdot\sum_{\pm}\int{\rm
d}z^{\prime}\;\Theta{(\pm z^{\prime})}\;{\rm
e}^{-\gamma|z^{\prime}|}\;\Big{[}\delta_{{\bf G}{\bf G}^{\prime}}{\rm e}^{{\rm
i}q_{{\bf
G}\varepsilon_{\ell}z}|z-z^{\prime}|}+r^{\operatorname{sign}\\{z\\},\pm}_{{\bf
G}{\bf G}^{\prime}}(\ell{\bf k}_{\parallel},\varepsilon_{\ell})\,{\rm e}^{{\rm
i}q_{{\bf G}\varepsilon_{\ell}z}|z|+{\rm i}q_{{\bf
G}^{\prime}\varepsilon_{\ell}z}|z^{\prime}|}\Big{]}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\;\,\times\Big{[}-\hat{\bf
z}\,q_{0}\,\delta_{{\bf G}^{\prime}0}\;{\rm e}^{-{\rm
i}q_{0}z^{\prime}}+B^{\pm}_{{\bf G}^{\prime}}\;({\bf G}^{\prime}\pm q_{{\bf
G}^{\prime}\varepsilon_{0}z}\hat{\bf z})\;{\rm e}^{{\rm i}q_{{\bf
G}^{\prime}\varepsilon_{0}z}|z^{\prime}|}\Big{]},$

where we have inserted the ${\rm e}^{-\gamma|z^{\prime}|}$ factor that we
discuss in Sec. A1.3. Finally, performing the $z^{\prime}$ integral with the
help of Eqs. (A27) and taking the $|z|\to\infty$ and $\gamma\to 0^{+}$ limits
in that order, we obtain

$\displaystyle\psi^{(1)}({\bf r},t)=\sum_{{\bf G},\ell,\pm}C_{{\bf
G},\ell}^{\pm}\;{\rm e}^{{\rm i}(\ell{\bf k}_{\parallel}+{\bf G})\cdot{\bf
R}\pm{\rm i}q_{{\bf G}\varepsilon_{\ell}z}z-{\rm
i}\varepsilon_{\ell}t}\;\Theta(\pm z),$

where

$\displaystyle C_{{\bf G},\ell}^{s}=$ $\displaystyle-\frac{{\rm
i}^{\ell}e}{\hbar\omega}\;\frac{1}{q_{{\bf G}\varepsilon_{\ell}z}}{\bf
E}_{0}\cdot\sum_{\pm}\Bigg{\\{}\hat{\bf z}\,r^{s\pm}_{{\bf G}0}(\ell{\bf
k}_{\parallel},\varepsilon_{\ell})\,\frac{q_{0}}{q_{{\bf
G}\varepsilon_{\ell}z}\mp q_{0}}$ (B34) $\displaystyle+\sum_{{\bf
G}^{\prime}}B^{\pm}_{{\bf G}^{\prime}}\;({\bf G}^{\prime}\pm q_{{\bf
G}^{\prime}\varepsilon_{0}z}\hat{\bf z})\bigg{[}\delta_{{\bf G}{\bf
G}^{\prime}}\frac{\pm s}{q_{{\bf G}\varepsilon_{\ell}z}\mp sq_{{\bf
G}^{\prime}\varepsilon_{0}z}}-r^{s\pm}_{{\bf G}{\bf G}^{\prime}}(\ell{\bf
k}_{\parallel},\varepsilon_{\ell})\,\frac{1}{q_{{\bf
G}\varepsilon_{\ell}z}+q_{{\bf
G}^{\prime}\varepsilon_{0}z}}\bigg{]}\Bigg{\\}}$

are the sought-after inelastic scattering coefficients.

## References

  * Barwick  _et al._ (2009) B. Barwick, D. J. Flannigan, and A. H. Zewail, Nature 462, 902 (2009).
  * Flannigan and Zewail (2012) D. J. Flannigan and A. H. Zewail, Acc. Chem. Res. 45, 1828 (2012).
  * Feist  _et al._ (2015) A. Feist, K. E. Echternkamp, J. Schauss, S. V. Yalunin, S. Schäfer, and C. Ropers, Nature 521, 200 (2015).
  * Barwick and Zewail (2015) B. Barwick and A. H. Zewail, ACS Photonics 2, 1391 (2015).
  * Polman  _et al._ (2019) A. Polman, M. Kociak, and F. J. García de Abajo, Nat. Mater. 18, 1158 (2019).
  * García de Abajo and Di Giulio (2021) F. J. García de Abajo and V. Di Giulio, ACS Photonics 8, 945 (2021).
  * Roques-Carmes  _et al._ (2023) C. Roques-Carmes, S. E. Kooi, Y. Yang, N. Rivera, P. D. Keathley, J. D. Joannopoulos, S. G. Johnson, I. Kaminer, K. K. Berggren, and M. Soljačić, Appl. Phys. Rev. 10, 011303 (2023).
  * García de Abajo and Ropers (2023) F. J. García de Abajo and C. Ropers, Phys. Rev. Lett. 130, 246901 (2023).
  * García de Abajo  _et al._ (2010) F. J. García de Abajo, A. Asenjo-Garcia, and M. Kociak, Nano Lett. 10, 1859 (2010).
  * Park  _et al._ (2010) S. T. Park, M. Lin, and A. H. Zewail, New J. Phys. 12, 123028 (2010).
  * Piazza  _et al._ (2015) L. Piazza, T. T. A. Lummen, E. Quiñonez, Y. Murooka, B. Reed, B. Barwick, and F. Carbone, Nat. Commun. 6, 6407 (2015).
  * Lummen  _et al._ (2016) T. T. A. Lummen, R. J. Lamb, G. Berruto, T. LaGrange, L. D. Negro, F. J. García de Abajo, D. McGrouther, B. Barwick, and F. Carbone, Nat. Commun. 7, 13156 (2016).
  * Vanacore  _et al._ (2020) G. M. Vanacore, I. Madan, and F. Carbone, Riv. Nuovo Cimento 43, 567 (2020).
  * Kurman  _et al._ (2021) Y. Kurman, R. Dahan, H. H. Sheinfux, K. Wang, M. Yannai, Y. Adiv, O. Reinhardt, L. H. G. Tizei, S. Y. Woo, J. Li, J. H. Edgar, M. Kociak, F. H. L. Koppens, and I. Kaminer, Science 372, 1181 (2021).
  * Nabben  _et al._ (2023) D. Nabben, J. Kuttruff, L. Stolz, A. Ryabov, and P. Baum, Nature 619, 63 (2023).
  * Gaida  _et al._ (2024) J. H. Gaida, H. Lourenço-Martins, M. Sivis, T. Rittmann, A. Feist, F. J. García de Abajo, and C. Ropers, “Attosecond electron microscopy by free-electron homodyne detection,”  (2024), Nat. Photon. https://doi.org/10.1038/s41566-024-01380-8, arXiv:2305.03005 .
  * Bucher  _et al._ (2023) T. Bucher, H. Nahari, H. H. Sheinfux, R. Ruimy, A. Niedermayr, R. Dahan, Q. Yan, Y. Adiv, M. Yannai, J. Chen, Y. Kurman, S. T. Park, D. J. Masiel, E. Janzen, J. H. Edgar, F. Carbone, G. Bartal, S. Tsesses, F. H. L. Koppens, G. M. Vanacore, and I. Kaminer, “Coherently amplified ultrafast imaging in a free-electron interferometer,”  (2023), arXiv:2305.04877 .
  * Di Giulio  _et al._ (2019) V. Di Giulio, M. Kociak, and F. J. García de Abajo, Optica 6, 1524 (2019).
  * Dahan  _et al._ (2021) R. Dahan, A. Gorlach, U. Haeusler, A. Karnieli, O. Eyal, P. Yousefi, M. Segev, A. Arie, G. Eisenstein, P. Hommelhoff, and I. Kaminer, Science 373, eabj7128 (2021).
  * Yang  _et al._ (2024) Y. Yang, J.-W. Henke, A. S. Raja, F. J. Kappert, G. Huang, G. Arend, Z. Qiu, A. Feist, R. N. Wang, A. Tusnin, A. Tikan, C. Ropers, and T. J. Kippenberg, Science 383, 168 (2024).
  * Priebe  _et al._ (2017) K. E. Priebe, C. Rathje, S. V. Yalunin, T. Hohage, A. Feist, S. Schäfer, and C. Ropers, Nat. Photon. 11, 793 (2017).
  * Kozák  _et al._ (2018) M. Kozák, N. Schönenberger, and P. Hommelhoff, Phys. Rev. Lett. 120, 103203 (2018).
  * Morimoto and Baum (2018) Y. Morimoto and P. Baum, Nat. Phys. 14, 252 (2018).
  * Vanacore  _et al._ (2019) G. M. Vanacore, G. Berruto, I. Madan, E. Pomarico, P. Biagioni, R. J. Lamb, D. McGrouther, O. Reinhardt, I. Kaminer, B. Barwick, H. Larocque, V. Grillo, E. Karimi, F. J. García de Abajo, and F. Carbone, Nat. Mater. 18, 573 (2019).
  * Konečná and García de Abajo (2020) A. Konečná and F. J. García de Abajo, Phys. Rev. Lett. 125, 030801 (2020).
  * García de Abajo and Konečná (2021) F. J. García de Abajo and A. Konečná, Phys. Rev. Lett. 126, 123901 (2021).
  * Madan  _et al._ (2022) I. Madan, V. Leccese, A. Mazur, F. Barantani, T. LaGrange, A. Sapozhnik, P. M. Tengdin, S. Gargiulo, E. Rotunno, J.-C. Olaya, I. Kaminer, V. Grillo, F. J. García de Abajo, F. Carbone, and G. M. Vanacore, ACS Photonics 9, 3215 (2022).
  * Mihaila  _et al._ (2022) M. C. C. Mihaila, P. Weber, M. Schneller, L. Grandits, S. Nimmrichter, and T. Juffmann, Phys. Rev. X 12, 031043 (2022).
  * Glover  _et al._ (1996) T. E. Glover, R. W. Schoenlein, A. H. Chin, and C. V. Shank, Phys. Rev. Lett. 76, 2468 (1996).
  * Saathoff  _et al._ (2008) G. Saathoff, L. Miaja-Avila, M. Aeschlimann, M. M. Murnane, and H. C. Kapteyn, Phys. Rev. A 77, 022903 (2008).
  * Arrell  _et al._ (2016) C. A. Arrell, J. Ojeda, L. Mewes, J. Grilj, F. Frassetto, L. Poletto, F. van Mourik, and M. Chergui, Phys. Rev. Lett. 117, 143001 (2016).
  * Yalunin  _et al._ (2011) S. V. Yalunin, M. Gulde, and C. Ropers, Phys. Rev. B 84, 195426 (2011).
  * Dombi  _et al._ (2020) P. Dombi, Z. Pápa, J. Vogelsang, S. V. Yalunin, M. Sivis, G. Herink, S. Schäfer, P. Groß, C. Ropers, and C. Lienau, Rev. Mod. Phys. 92, 025003 (2020).
  * Bendaña  _et al._ (2011) X. M. Bendaña, A. Polman, and F. J. García de Abajo, Nano Lett. 11, 5099 (2011).
  * Kfir  _et al._ (2020) O. Kfir, H. Lourenço-Martins, G. Storeck, M. Sivis, T. R. Harvey, T. J. Kippenberg, A. Feist, and C. Ropers, Nature 582, 46 (2020).
  * Wang  _et al._ (2020) K. Wang, R. Dahan, M. Shentcis, Y. Kauffmann, A. Ben Hayun, O. Reinhardt, S. Tsesses, and I. Kaminer, Nature 582, 50 (2020).
  * Dahan  _et al._ (2020) R. Dahan, S. Nehemia, M. Shentcis, O. Reinhardt, Y. Adiv, X. Shi, O. Be’er, M. H. Lynch, Y. Kurman, K. Wang, and I. Kaminer, Nat. Phys. 16, 1123 (2020).
  * Henke  _et al._ (2021) J.-W. Henke, A. S. Raja, A. Feist, G. Huang, G. Arend, Y. Yang, F. J. Kappert, R. N. Wang, M. Möller, J. Pan, J. Liu, O. Kfir, C. Ropers, and T. J. Kippenberg, Nature 600, 653 (2021).
  * Varshalovich and D’Yakonov (1971) D. A. Varshalovich and M. A. D’Yakonov, Sov. Phys. JETP 33, 51 (1971).
  * Weingartshofer  _et al._ (1977) A. Weingartshofer, J. K. Holmes, G. Caudle, E. M. Clarke, and H. Krüger, Phys. Rev. Lett. 39, 269 (1977).
  * Talebi (2018) N. Talebi, Adv. Phys. X 3, 1499438 (2018).
  * Talebi (2020) N. Talebi, Phys. Rev. Lett. 125, 080401 (2020).
  * García de Abajo  _et al._ (2022) F. J. García de Abajo, E. J. C. Dias, and V. Di Giulio, Phys. Rev. Lett. 129, 093401 (2022).
  * Lord Rayleigh (1907a) Lord Rayleigh, Proc. R. Soc. Lond. A 79, 399 (1907a).
  * Abd El-Fattah  _et al._ (2019) Z. M. Abd El-Fattah, V. Mkhitaryan, J. Brede, L. Fernández, C. Li, Q. Guo, A. Ghosh, A. Rodríguez Echarri, D. Naveh, F. Xia, J. E. Ortega, and F. J. García de Abajo, ACS Nano 13, 7771 (2019).
  * Woessner  _et al._ (2015) A. Woessner, M. B. Lundeberg, Y. Gao, A. Principi, P. Alonso-González, M. Carrega, K. Watanabe, T. Taniguchi, G. Vignale, M. Polini, J. Hone, R. Hillenbrand, and F. H. Koppens, Nat. Mater. 14, 421 (2015).
  * Ni  _et al._ (2018) G. X. Ni, A. S. McLeod, Z. Sun, L. Wang, L. Xiong, K. W. Post, S. S. Sunku, B.-Y. Jiang, J. Hone, C. R. Dean, M. M. Fogler, and D. N. Basov, Nature 557, 530 (2018).
  * Giles  _et al._ (2018) A. J. Giles, S. Dai, I. Vurgaftman, T. Hoffman, S. Liu, L. Lindsay, C. T. Ellis, N. Assefa, I. Chatzakis, T. L. Reinecke, J. G. Tischler, M. M. Fogler, J. H. Edgar, D. N. Basov, and J. D. Caldwell, Nat. Mater. 17, 134 (2018).
  * Li  _et al._ (2021) N. Li, X. Guo, X. Yang, R. Qi, T. Qiao, Y. Li, R. Shi, Y. Li, K. Liu, Z. Xu, L. Liu, F. J. García de Abajo, Q. Dai, E.-G. Wang, and P. Gao, Nat. Mater. 20, 43 (2021).
  * Li  _et al._ (2014) Y. Li, A. Chernikov, X. Zhang, A. Rigosi, H. M. Hill, A. M. van der Zande, D. A. Chenet, E.-M. Shih, J. Hone, and T. F. Heinz, Phys. Rev. B 90, 205422 (2014).
  * Epstein  _et al._ (2020) I. Epstein, B. Terrés, A. J. Chaves, V.-V. Pusapati, D. A. Rhodes, B. Frank, V. Zimmermann, Y. Qin, K. Watanabe, T. Taniguchi, H. Giessen, S. Tongay, J. C. Hone, N. M. R. Peres, and F. H. L. Koppens, Nano Lett. 20, 3545 (2020).
  * de Aguiar (1993) M. A. M. de Aguiar, Phys. Rev. A 48, 2567 (1993).
  * Vanacore  _et al._ (2018) G. M. Vanacore, I. Madan, G. Berruto, K. Wang, E. Pomarico, R. J. Lamb, D. McGrouther, I. Kaminer, B. Barwick, F. J. García de Abajo, and F. Carbone, Nat. Commun. 9, 2694 (2018).
  * Mak  _et al._ (2008) K. F. Mak, M. Y. Sfeir, Y. Wu, C. H. Lui, J. A. Misewich, and T. F. Heinz, Phys. Rev. Lett. 101, 196405 (2008).
  * Nair  _et al._ (2008) R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim, Science 320, 1308 (2008).
  * Rocca (1995) M. Rocca, Surf. Sci. Rep. 22, 1 (1995).
  * Pendry (1974) J. B. Pendry,  _Low Energy Electron Diffraction_ (Academic Press, London, 1974).
  * Pendry (1984) J. B. Pendry, in  _Determination of Surface Structure by LEED_ , edited by P. M. Marcus and F. Jona (Plenum Press, New York, 1984).
  * Claus  _et al._ (1992) H. Claus, A. Büssenschütt, and M. Henzler, Rev. Sci. Instrum. 63, 2195 (1992).
  * Nagao  _et al._ (2001) T. Nagao, T. Hildebrandt, M. Henzler, and S. Hasegawa, Phys. Rev. Lett. 86, 5747 (2001).
  * Nagao  _et al._ (2006) T. Nagao, S. Yaginuma, T. Inaoka, and T. Sakurai, Phys. Rev. Lett. 97, 116802 (2006).
  * Gulde  _et al._ (2014) M. Gulde, S. Schweda, G. Storeck, M. Maiti, H. K. Yu, A. M. Wodtke, S. Schäfer, and C. Ropers, Science 345, 200 (2014).
  * Vogelgesang  _et al._ (2018) S. Vogelgesang, G. Storeck, J. G. Horstmann, T. Diekmann, M. Sivis, S. Schramm, K. Rossnagel, S. Schäfer, and C. Ropers, Nat. Phys. 14, 184 (2018).
  * Lord Rayleigh (1907b) Lord Rayleigh, Philos. Mag. 14, 60 (1907b).
  * García de Abajo (2007) F. J. García de Abajo, Rev. Mod. Phys. 79, 1267 (2007).
  * Weingartshofer  _et al._ (1983) A. Weingartshofer, J. K. Holmes, J. Sabbagh, and S. L. Chin, J. Phys. B 16, 1805 (1983).
  * Francken and C.J.Joachain (1990) P. Francken and C.J.Joachain, J. Opt. Soc. Am. B 7, 554 (1990).
  * E. Arqué López and García de Abajo (2022) V. D. E. Arqué López and F. J. García de Abajo, Phys. Rev. Research 4, 013241 (2022).
  * Economou (2006) E. N. Economou,  _Green’s Functions in Quantum Physics_ (Springer, Heidelberg, 2006).
  * Messiah (1966) A. Messiah,  _Quantum Mechanics_ (North-Holland, New York, 1966).
  * Sakurai (1994) J. J. Sakurai,  _Modern Quantum Mechanics_ (Addison-Wesley, Boston, 1994).
  * Wolkow (1935) D. M. Wolkow, Z. Phys. 94, 250 (1935).
  * Van Hove  _et al._ (1986) M. A. Van Hove, W. H. Weinberg, and C.-M. Chan,  _Low Energy Electron Diffraction_ (Springer-Verlag, Heidelberg, 1986).
  * Ankudinov  _et al._ (1996) A. L. Ankudinov, S. I. Zabinsky, and J. J. Rehr, Comput. Phys. Commun. 98, 359 (1996).
  * Salvat and Mayol (1991) F. Salvat and R. Mayol, Comput. Phys. Commun. 62, 65 (1991).
  * García de Abajo  _et al._ (2001) F. J. García de Abajo, M. A. Van Hove, and C. S. Fadley, Phys. Rev. B 63, 075404 (2001).
  * Kambe (1967) K. Kambe, Z. Naturforsch. A 22, 322 (1967).

Figure S7: Inelastic scattering of low-energy electrons upon total reflection
at a polariton-supporting surface. Panels (b-d,g-i) are reproduced from Fig.
2a-c,e-g in the main text, but now including dashed curves indicating the
kinematic threshold beyond which energy–momentum conservation renders $q_{\ell
z}$ purely imaginary (i.e., $1+\ell\omega/\varepsilon_{0}=|\sin\theta_{\rm
e}+\ell k_{\parallel}/q_{0}|^{2}$). In addition, we present results for
$\ell=\pm 2$ in panels (a,e,f,j). Figure S8: Influence of electron incidence
angle on inelastic scattering under the conditions of Fig. 2 in the main text.
We present results analogous to Fig. 2d,h, but for different values of the
electron incidence angle $\theta_{\rm e}$.  Figure S9: Inelastic scattering of
low-energy electrons at an electron-transparent, polariton-supporting thin
film. Same as Fig. 3 in the main text, but with $U_{0}=0$ (i.e., the electron
does not see the material).