# TocBERT: Medical Document Structure Extraction Using Bidirectional
Transformers

Sarra Baghdadi12, Majd Saleh12 and Stéphane Paquelet1  1Institute of Research
and Technology b-com, Rennes, France
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>2Equal contribution

###### Abstract

Text segmentation holds paramount importance in the field of Natural Language
Processing (NLP). It plays an important role in several NLP downstream tasks
like information retrieval and document summarization. In this work, we
propose a new solution, namely TocBERT, for segmenting texts using
bidirectional transformers. TocBERT represents a supervised solution trained
on the detection of titles and sub-titles from their semantic representations.
This task was formulated as a named entity recognition (NER) problem. The
solution has been applied on a medical text segmentation use-case where the
Bio-ClinicalBERT model is fine-tuned to segment discharge summaries of the
MIMIC-III dataset. The performance of TocBERT has been evaluated on a human-
labeled ground truth corpus of $250$ notes. It achieved an F1-score of
$84.6\%$ when evaluated on a linear text segmentation problem and $72.8\%$ on
a hierarchical text segmentation problem. It outperformed a carefully designed
rule-based solution, particularly in distinguishing titles from subtitles.

###### Keywords:

Title detection; text segmentation; NLP; language models; transformers, BERT,
information retrieval; medical text cleaning.

##  I Introduction

Text segmentation is the task of partitioning a document into topically
coherent segments [1, 2, 3, 4]. It represents an important pre-processing step
in several Natural Language Processing (NLP) downstream tasks inducing
information retrieval and summarization[1, 2, 3, 4, 5, 6]. Two types of text
segmentation can be distinguished: 1) linear segmentation where a document is
divided into sequential contiguous segments, and 2) hierarchical segmentation
where higher level segments are further segmented into smaller sub-segments
[3].

Existing solutions to text segmentation are applied either to unformatted,
i.e. plain-text, documents like [3, 4] or to formatted ones like [7, 8, 9]. In
the first family, recent text segmentation methods are based on the semantic
representation of words and/or of sentences. For instance, word embeddings and
sentence embeddings are used to detect topic changes in a sequence of
sentences in [3, 4]. In the second family, document layout, orthographic
features, geometric features and stylistic properties, e.g. font properties
and line spaces, are exploited to extract the document structure. Some of the
aforementioned features can be extracted using the Document Image Analysis
(DIA) [10, 9] while others can be extracted using regex matches of predefined
patterns [8]. In both families, titles and subtitles play an important role in
text segmentation. For example, in semantic-aspects-based methods, supervised
solutions like [3] and [4] are trained on a large corpus named WIKI-727K [4]
which is automatically labeled such that topic-change cutoff points are
completely dependent on titles. Similarly, visual-aspects-based methods like
[9] are mainly based on title detection.

Figure 1: An extract from a discharge summary from the MIMIC-III database
annotated with titles (yellow) and subtitles (blue)

In this article, the focus is set on the hierarchical segmentation of
unformatted medical documents by detecting titles and subtitles using their
semantic vector representations. Particularly, the considered use-case is the
segmentation of the discharge summaries of the MIMIC-III database [11] which
represents an unformatted free-text corpus. Figure 1 shows an extract of a
discharge summary where we highlighted titles in yellow and subtitles in blue.
The motivation underlying the aforementioned segmentation task is two-fold:

  * •

To clean the MIMIC-III discharge summaries corpus such that it can be
efficiently used to domain-adapt pretrained language models. For example,
noisy sections like ”Discharge medications” or ”Admission labs” can be
detected and removed since their content doesn’t represent a fluent language
and can harm the domain adaption task.

  * •

To enable the extraction of specialized sub-corpora and facilitate the
information retrieval tasks. For instance, sub-corpora about radiology reports
or cardiac diseases reports can be extracted and used to train relevant
information extraction systems.

While the aforementioned advantages have been discussed on the considered use-
case, the proposed method is general and can be extended to solve many similar
tasks as will be shown in Section IV.

The rest of the paper is organized as follows: we start by exploring related
works in Section II. Training and test corpora are described in Section III.
The proposed method is presented in Section IV. Experimental results are
discussed in Section V while Section VI concludes the paper and discuses
future work.

##  II Related work

Figure 2: Two-level neural networks for text segmentation: the first-level
network builds sentence embeddings from token embeddings while the second-
level network labels the sequence of sentence representations as: topic-change
$(1)$ or no-topic-change $(0)$

One remarkable family of solutions to topical text segmentation was proposed
in [4, 3]. As depicted in Figure 2, these methods are based on two-level
neural networks that divide the segmentation problem into two sub-problems: 1)
starting from token embeddings, build sentence embeddings, and 2) based on the
sentence representations, annotate the sentences’ sequence with binary labels
i.e. $0$ for no-topic-change and $1$ for topic-change. Specifically, these
methods start by dividing the input document into sentences (sentence
tokenization) using the NLTK PUNKT tokenizer [12]. Then, they use pretrained
token embeddings to represent input tokens $\textbf{w}_{i}$. Concretely,
word2vec embeddings [13, 14] are used in [4] while FastText embeddings [15]
are used in [3]. Token embeddings of each sentence are fed to a sentence
embedding model in order to create sentence representations i.e. each sentence
is represented by a dense vector $\textbf{s}_{j}$. This sentence embedding
model is a 2-layer bidirectional LSTM [16] in the method proposed in [4] while
it represents a transformer encoder [17] in the method proposed in [3]. The
sequence of sentence embeddings is then fed to a sentence labeling model which
assigns a binary label $l_{i}\in\\{0,1\\}$ to each sentence where $l_{i}=1$
encoding the existence of a topic change. This sentence labeling model is a
2-layer bidirectional LSTM in [4] while it is a transformer encoder in [3].
Both of the aforementioned solutions are trained in a supervised manner where
the labeled corpus WIKI-727K [4] is used for training. Note that the
transformer-based solution adds another feed-forward head besides the sentence
labeling model. This head is responsible for explicit coherence modeling [3].
A similar work was proposed in [18] where authors tested a solution based on a
transformer to encode sentences (first level) and a bidirectional LSTM to
detect topic changes (second level). For deeper details about how transformers
and LSTMs work, interested readers are referred to this tutorial [19].

While the aforementioned solutions are supervised, many solutions based on
unsupervised-learning have been proposed. For instance, the method proposed in
[20] uses RoBERTa embeddings [21] to calculate semantic similarity scores
between sentences. These scores are then used to perform segmentation.
Particularly, the variation in similarity scores over time is leveraged to
recognize changes in topics. The problem of representing sentences from token
embeddings was addressed in two ways: 1) to use sentence-BERT [22], a version
of BERT [23] fine-tuned on sentence representation tasks, and 2) to aggregate
token embeddings from the 2nd to the last layer of RoBERTa using max pooling.

Before the deep-learning era, several interesting statistical solutions were
proposed. Interested readers are referred for example to 1) TextTiling [24]
which captures semantic similarity between sentences using word frequency
vectors, 2) TopicTiling [25] an algorithm based on TextTiling [24] and Latent
Dirichlet Allocation (LDA) [26], 3) LDA-based topic modeling [27], and 4) the
statistical models for text segmentation proposed in [2].

The solutions discussed so far represent the semantic-aspects-based text
segmentation family. Let us explore now some visual-aspects-based text
segmentation methods i.e. segmentation applied on formatted documents. In
fact, in the financial domain, several interesting solutions to the document-
structure-extraction problem have been proposed thanks to the FinTOC shared
tasks [8, 7, 9]. These latter aim at extracting the structure of complex
layout formatted documents stored in pdf format. The objective is to detect
titles and to hierarchically organize them into a Table of Contents (ToC). For
instance, authors of [9] proposed a title detection algorithm based on
Document Image Analysis (DIA). Particularly, they started by Faster R-CNN [28]
model which is pertained on the PubLayNet dataset [29] and they fine tuned it
on the FinTOC-2022 training set. In order to recognize the hierarchical levels
of the detected titles, they extracted orthographic and layout features from
the detected titles and they proposed a supervised random forest model to
predict the title level.

##  III Data analysis and annotation

The experimental work has been performed on the discharge summaries corpus, a
sub-corpus from the MIMIC-III [11] database. MIMIC-III contains more than two
million reports including discharge summaries, EEG reports, radiology reports,
and many others. Figure 3-a shows the distribution of MIMIC-III reports’
families. The discharge summaries corpus contains around $60~{}000$ reports.
Figure 3-b shows the distribution of discharge summaries’ lengths given in
number of tokens. The average discharge summary length is $1435$ tokens
(roughly $3$ pages).

We divided the discharge summaries corpus to a test corpus of $250$ randomly
selected reports and a training corpus containing the rest of reports. The
test corpus has been manually annotated using the web application ACUITEE [30]
as illustrated in Figure 1. Note that two hierarchical level of titles are
considered: titles (highlighted in yellow) and sub-titles (highlighted in
blue).

In order to annotate the training corpus, we conducted a deep analysis of the
discharge summaries structure and we proposed a rule-based title detection
system, namely TocRegex. This latter has been implemented using regular
expressions. For example, the first considered pattern was any character
sequence that:

  * •

begins with a newline ”\n” followed by a valid title content and terminates
with a colon followed by newline ”:\n”; or

  * •

begins with a double newline ”\n\n” followed by a valid title content and
terminates with a colon ”:”; or

  * •

begins with a beginning-of-document character followed by a valid title
content and terminates with a colon ”:”

The definition of ”valid title content” as well as the full list of the $12$
considered regular expressions are discussed in the Appendix A.
Figure 4 shows the frequencies of the top 35 matches detected using the above-
described pattern. For instance, the title ”history of present illness” has
been detected $10e5$ times in the corpus. We note that several false positives
have also matched the considered pattern like e.g. ”tablet(s)*refills” which
appears frequently in medication lists. In order to filter-out such kind of
false positives, a manual curation has been applied to the detected titles. As
a result, the total remaining numbers of unique titles is $56~{}357$. They
have appeared $2~{}930~{}725$ times in the corpus i.e. the average number of
titles is around $49$ per discharge summary.

Figure 3: a) distribution of sub-corpora sizes of the MIMIC-III database. 2)
distribution of reports’ lengths of the discharge summaries corpus.  Figure 4:
Top $35$ ”candidate” titles extracted using the first pattern of the TocRegex
solution including true positives (blue) and false positives (light brown)

##  IV The proposed solution: TocBERT

In this section, we explore the proposed solution TocBERT (Table of Content
BERT). The hierarchical title detection task is formulated as a sequence
labeling problem i.e. token classification problem. Particularly, it is
considered as a named entity recognition (NER) problem with three entity
types: ”I-title” for titles, ”I-Stitle” for subtitles and ”O” for other
tokens. Note that we adopt the ”Inside–outside–beginning” (IOB) standard
labeling format. Figure 5 shows two examples of labeled sequences where tokens
are tagged with title, subtitle or outside labels. For instance, the second
sequence has one title ”Physical exam” and four subtitles ”HEENT”, ”Neck”,
”Lungs” and ”Extremities”.

Figure 5: Labeling token sequences in TocBERT: Example 1) two titles ”Past
Medical History” and ”Family History” are labeled with ”I-title” while other
tokens are labeled with ”O”; Example 2) one title ”Physical exam” is labeled
with ”I-title” and four subtitles ”HEENT”, ”Neck”, ”Lungs” and ”Extremities”
are labeled with ”I-Stitle”, while other (outside) tokens are labeled with
”O”.

TocBERT is based on fine-tuning the pretrained model Bio-ClinicalBERT [31] on
the above-mentioned NER task. Note that bidirectional transformers, like the
BERT family, are more convenient for tackling sequence labeling tasks, like
NER, compared to generative i.e. auto-regressive transformers like GPT [32].
Bio-ClinicalBERT represents a variant of BERT adapted to biomedical and
clinical domains. We fine-tuned this pretrained model using the discharge
summaries training corpus that has been semi-automatically labeled using
TocRegex followed by a manual curation as described in the previous section
III.
Bio-ClinicalBERT, like BERT, is pretrained using a fixed size vocabulary:
$28~{}996$ tokens. This vocabulary was created using the word-piece algorithm
[33]. Some of the resultant tokens represent sub-words. For this reason, the
first step of TocBERT is to project labels of the full-word tokens (extracted
by simple pre-tokenization procedure) to sub-word tokens (extracted using Bio-
ClinicalBERT tokenizer). The second step consists in preparing convenient
training windows from the training corpus. This is important because the
maximum window size of BERT is $512$ (sub-word) tokens while the average
length of a discharge summary is $1435$ (full-word) tokens. To this end,
discharge summaries have been segmented into windows of $384$ words where we
used the approximate formula: $1$ token $=0.75$ words. The total size of the
resultant training set is $144~{}000$ labeled windows. Finally, the last step
consists of training the TocBERT model.

##  V Experimental results

###  V-A Experimental configurations

TocBERT was trained on one NVIDIA A100 GPU equipped with 80 GB of RAM. It was
trained for $20$ epochs with a batch size of $16$ training samples. The
training time was around $17$ hours. The inference hardware was an NVIDIA RTX
A3000 laptop GPU equipped with 6 GB of dedicated RAM.

###  V-B Results

Table I shows the experimental results of the proposed text segmentation
solutions, TocBERT and TocRegex, in the hierarchical segmentation
configurations, i.e. detecting both titles and sub-titles. While TocBERT is
trained on a corpus labeled using TocRegex, the former considerably
outperforms the latter in all the considered performance criteria, i.e.
precision, recall and F1-score.

TABLE I: Experimental results: hierarchical text segmentation | Precision | Recall | F1-score
---|---|---|---
TocBERT | 0.714 | 0.754 | 0.728
TocRegex | 0.667 | 0.563 | 0.606

Table II shows the results of the aforementioned solutions in the linear text
segmentation configurations. TocBERT and TocRegex show comparable overall
performance measured by the F1-score metric. TocBERT shows higher sensitivity
(recall) while it is less specific. Note that TocBERT is completely based on
the semantic representation of titles and their context. It doesn’t depend on
the existence of visual aspects like newlines and colons. This explains why it
detects more titles (higher recall). On the other hand, TocRegex is more
specific since it is entirely based on patterns that almost always exist only
in titles.

Comparing the results of hierarchical and linear segmentation, the value of
TocBERT is clearly its strong capacity to exploit the context and the semantic
aspects to discriminate between titles and sub-titles.

TABLE II: Experimental results: linear text segmentation | Precision | Recall | F1-score
---|---|---|---
TocBERT | 0.829 | 0.878 | 0.846
TocRegex | 0.932 | 0.784 | 0.845
TABLE III: Regular expressions used to detect titles
---

In terms of inference execution time, TocRegex takes $87$ ms while TocBERT
takes $195$ ms, in average, to segment a discharge summary. While TocRegex is
faster, both solutions satisfy near-real-time requirements.

Figure 6: Grouping of some key concepts towards constructing topics’ ontology

Finally, we observed that the detected groups of titles and subtitles can be
exploited to automatically construct an ontology of topics. Figure 6 shows an
example of initial manual grouping of key concepts that help constructing such
ontology. The semantic vector representations of the detected titles can play
an important role in organizing the topics’ ontology. The interest of building
this ontology is to facilitate tasks like text segmentation and information
retrieval. This idea will be investigated in future works.

##  VI Conclusion

In this paper, we proposed a new solution, TocBERT, for the hierarchical
segmentation of medical reports. The segmentation task was formulated as a
named entity recognition problem. TocBERT was initialized by a pretrained
model, Bio-ClinicalBERT, and fine tuned on a the MIMIC-III discharge summaries
corpus. This latter was semi-automatically labeled with titles and sub-titles.
TocBERT showed very good results considerably outperforming a carefully-
designed rule-based system. Particularly, it showed a good performance in
discriminating between titles and subtitles by leveraging their semantic
representations and employing their context.

The semantic representations of the extracted titles can be exploited to
automatically construct an ontology of topics. Such an ontology can further
facilitate tasks like text segmentation and information retrieval. The
investigation of this idea is left for future work.

##  Appendix A TocRegex

Table III shows the regular expressions list used in the proposed rule-based
solution, TocRegex. In the first line, the four supported title-content
patterns are defined while the remaining $12$ lines list the full-title
patterns.

## References

  * [1] A. A. Alemi and P. Ginsparg, “Text segmentation based on semantic word embeddings,” arXiv, 2015.  https://arxiv.org/abs/1503.05543.
  * [2] D. Beeferman, A. Berger, and J. Lafferty, “Statistical models for text segmentation,” Machine Learning, vol. 34, pp. 177–210, 1999.  https://link.springer.com/article/10.1023/A:1007506220214.
  * [3] G. Glavas and S. Somasundaran, “Two-level transformer and auxiliary coherence modeling for improved text segmentation,” arXiv, 2020.  http://arxiv.org/abs/2001.00891.
  * [4] O. Koshorek, A. Cohen, N. Mor, M. Rotman, and J. Berant, “Text segmentation as a supervised learning task,” in Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 2 (Short Papers) (M. Walker, H. Ji, and A. Stent, eds.), (New Orleans, Louisiana), pp. 469–473, Association for Computational Linguistics, June 2018.  https://aclanthology.org/N18-2075.
  * [5] L. Miculicich and B. Han, “Document summarization with text segmentation,” arXiv, 2023.  https://arxiv.org/abs/2301.08817.
  * [6] Y. Zhu, H. Yuan, S. Wang, J. Liu, W. Liu, C. Deng, H. Chen, Z. Dou, and J.-R. Wen, “Large language models for information retrieval: A survey,” arXiv, 2024.  https://arxiv.org/abs/2308.07107.
  * [7] E. Giguet, G. Lejeune, and J.-B. Tanguy, “Daniel@FinTOC’2 shared task: Title detection and structure extraction,” in Proceedings of the 1st Joint Workshop on Financial Narrative Processing and MultiLing Financial Summarisation, 2020.  https://aclanthology.org/2020.fnp-1.30/.
  * [8] I. Maarouf, J. Kang, A. Azzi, S. Bellato, M. Gan, and M. El-Haj, “The financial document structure extraction shared task (FinTOC2021): Fintoc 2021,” in 3rd Financial Narrative Processing Workshop (FNP 2021), pp. 111–119, Oct. 2021.  https://aclanthology.org/2021.fnp-1.21.pdf.
  * [9] P. Cassotti, C. Musto, M. DeGemmis, G. Lekkas, and G. Semeraro, “swapuniba@fintoc2022: Fine-tuning pre-trained document image analysis model for title detection on the financial domain,” in Proceedings of the 4th Financial Narrative Processing Workshop, 2022.  http://www.lrec-conf.org/proceedings/lrec2022/workshops/FNP/index.html.
  * [10] L. O’Gorman and R. Kasturi, eds., Document image analysis.  Washington, DC, USA: IEEE Computer Society Press, 1997.  https://cse.usf.edu/~r1k/DocumentImageAnalysis/DIA.pdf.
  * [11] A. E. Johnson, T. J. Pollard, L. Shen, L.-w. H. Lehman, M. Feng, M. Ghassemi, B. Moody, P. Szolovits, L. Anthony Celi, and R. G. Mark, “Mimic-iii, a freely accessible critical care databases,” Scientific Data, vol. 3, no. 1, 2016.  https://doi.org/10.1038/sdata.2016.35.
  * [12] S. Bird, E. Klein, and E. Loper, Natural language processing with Python: analyzing text with the natural language toolkit.  O’Reilly Media, Inc., 2009.  https://www.nltk.org/book/.
  * [13] T. Mikolov, K. Chen, G. Corrado, and J. Dean, “Efficient estimation of word representations in vector space,” in Proceedings of the International Conference on Learning Representations ICLR 2013, (Scottsdale, Arizona, USA), Association for Computational Linguistics, July 2013.  https://doi.org/10.48550/arXiv.1301.3781.
  * [14] T. Mikolov, I. Sutskever, K. Chen, G. S. Corrado, and J. Dean, “Distributed representations of words and phrases and their compositionality,” in Proceedings of the Advances in Neural Information Processing Systems 26 (NIPS 2013), (Stateline, Nevada, USA), Association for Computational Linguistics, July 2013.  https://proceedings.neurips.cc/paper_files/paper/2013/file/9aa42b31882ec039965f3c4923ce901b-Paper.pdf.
  * [15] P. Bojanowski, E. Grave, A. Joulin, and T. Mikolov, “Enriching Word Vectors with Subword Information,” Transactions of the Association for Computational Linguistics, vol. 5, pp. 135–146, 06 2017.  https://doi.org/10.1162/tacl_a_00051.
  * [16] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural Comput., vol. 9, p. 1735–1780, nov 1997.  https://doi.org/10.1162/neco.1997.9.8.1735.
  * [17] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, “Attention is all you need,” in Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, (Red Hook, NY, USA), p. 6000–6010, Curran Associates Inc., 2017.  https://dl.acm.org/doi/pdf/10.5555/3295222.3295349.
  * [18] M. Lukasik, B. Dadachev, K. Papineni, and G. Simões, “Text segmentation by cross segment attention,” in Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP), 2020.  https://aclanthology.org/2020.emnlp-main.380.
  * [19] M. Saleh and S. Paquelet, “Anatomy of neural language models,” arXiv, 2024\.  https://arxiv.org/abs/2401.03797.
  * [20] A. Solbiati, K. Heffernan, G. Damaskinosa, S. Poddar, S. Modi, and J. Cali, “Unsupervised topic segmentation of meetings with bert embeddings,” arXiv, 2021.  https://arxiv.org/abs/2106.12978.
  * [21] Y. Liu, M. Ott, N. Goyal, J. Du, M. Joshi, D. Chen, O. Levy, M. Lewis, L. Zettlemoyer, and V. Stoyanov, “RoBERTa: A Robustly Optimized BERT Pretraining Approach,” 2020.  https://openreview.net/forum?id=SyxS0T4tvS.
  * [22] N. Reimers and I. Gurevych, “Sentence-BERT: Sentence embeddings using Siamese BERT-networks,” in Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP) (K. Inui, J. Jiang, V. Ng, and X. Wan, eds.), (Hong Kong, China), pp. 3982–3992, Association for Computational Linguistics, Nov. 2019.  https://aclanthology.org/D19-1410.
  * [23] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” arXiv, 2019.  https://arxiv.org/abs/1810.04805.
  * [24] M. A. Hearst, “TextTiling: segmenting text into multi-paragraph subtopic passages,” Comput. Linguist., vol. 23, p. 33–64, mar 1997.  https://dl.acm.org/doi/10.5555/972684.972687.
  * [25] M. Riedl and C. Biemann, “TopicTiling: A text segmentation algorithm based on LDA,” in Proceedings of ACL 2012 Student Research Workshop (J. C. K. Cheung, J. Hatori, C. Henriquez, and A. Irvine, eds.), (Jeju Island, Korea), pp. 37–42, Association for Computational Linguistics, July 2012\.  https://aclanthology.org/W12-3307.
  * [26] D. M. Blei, A. Y. Ng, and M. I. Jordan, “Latent dirichlet allocation,” J. Mach. Learn. Res., vol. 3, p. 993–1022, mar 2003.  https://dl.acm.org/doi/pdf/10.5555/944919.944937.
  * [27] M. Riedl and C. Biemann, “Text segmentation with topic models,” Journal for Language Technology and Computational Linguistics (JLCL), vol. 27, no. 47-69, pp. 13–24, 2012.  https://www.inf.uni-hamburg.de/en/inst/ab/lt/publications/2012-riedletal-jlcl.pdf.
  * [28] S. Ren, K. He, R. Girshick, and J. Sun, “Faster R-CNN: Towards real-time object detection with region proposal networks,” in Advances in Neural Information Processing Systems (C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett, eds.), vol. 28, Curran Associates, Inc., 2015.  https://proceedings.neurips.cc/paper_files/paper/2015/file/14bfa6bb14875e45bba028a21ed38046-Paper.pdf.
  * [29] X. Zhong, J. Tang, and A. J. Yepes, “PubLayNet: Largest dataset ever for document layout analysis,” in 2019 International Conference on Document Analysis and Recognition (ICDAR), (Los Alamitos, CA, USA), pp. 1015–1022, IEEE Computer Society, sep 2019.  https://doi.ieeecomputersociety.org/10.1109/ICDAR.2019.00166.
  * [30] M. Saleh, “ACUITEE : Annotation and curation user interface for terms extraction engines,” 2021.  https://acuitee.labs.b-com.com/.
  * [31] E. Alsentzer, J. Murphy, W. Boag, W.-H. Weng, D. Jindi, T. Naumann, and M. McDermott, “Publicly available clinical BERT embeddings,” in Proceedings of the 2nd Clinical Natural Language Processing Workshop (A. Rumshisky, K. Roberts, S. Bethard, and T. Naumann, eds.), (Minneapolis, Minnesota, USA), pp. 72–78, Association for Computational Linguistics, June 2019\.  https://aclanthology.org/W19-1909.
  * [32] A. Radford, K. Narasimhan, T. Salimans, and I. Sutskever, “Improving language understanding by generative pre-training,” Technical report, vol. 3, no. 1, 2018.  https://www.mikecaptain.com/resources/pdf/GPT-1.pdf.
  * [33] Y. Wu, M. Schuster, Z. Chen, Q. V. Le, M. Norouzi, W. Macherey, M. Krikun, Y. Cao, Q. Gao, K. Macherey, J. Klingner, A. Shah, M. Johnson, X. Liu, L. Kaiser, S. Gouws, Y. Kato, T. Kudo, H. Kazawa, K. Stevens, G. Kurian, N. Patil, W. Wang, C. Young, J. R. Smith, J. Riesa, A. Rudnick, O. Vinyals, G. S. Corrado, M. Hughes, and J. Dean, “Google’s neural machine translation system: Bridging the gap between human and machine translation,” ArXiv, vol. abs/1609.08144, 2016.  https://api.semanticscholar.org/CorpusID:3603249.


# Complete Excitation of Discrete Quantum Systems by Single Free Electrons

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  Eduardo
J. C. Dias  ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of
Science and Technology, 08860 Castelldefels (Barcelona), Spain  Valerio Di
Giulio  ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of
Science and Technology, 08860 Castelldefels (Barcelona), Spain

###### Abstract

We reveal a wealth of nonlinear and recoil effects in the interaction between
individual low-energy electrons and samples comprising a discrete number of
states. Adopting a quantum theoretical description of combined free-electron
and two-level systems, we find a maximum achievable excitation probability of
100%, which requires specific conditions relating to the coupling strength and
the transition symmetry, as we illustrate through calculations for dipolar and
quadrupolar modes. Strong recoil effects are observed when the kinetic energy
of the probe lies close to the transition threshold, although the associated
probability remains independent of the electron wave function even when fully
accounting for nonlinear interactions with arbitrarily complex multilevel
samples. Our work reveals the potential of free electrons to control localized
excitations and delineates the boundaries of such control.

##  I Introduction

Free electron beams (e-beams) allow us to image material nanostructures and
their excitations with an unsurpassed combination of space-energy resolution
in the subångstrom-meV domain thanks to a sustained series of advances in
electron microscope instrumentation over the last decades Nellist and
Pennycook (1998); Batson et al. (2002); Krivanek et al. (2014); Lagos et al.
(2017); Krivanek et al. (2019); Hage et al. (2020); Hachtel et al. (2019);
Mkhitaryan et al. (2021); Yan et al. (2021). In particular, electron energy-
loss spectroscopy (EELS) is widely used to identify localized excitations and
map their spatial distributions with atomic precision Egerton (1996, 2003);
Erni and Browning (2005); Brydson (2001); García de Abajo (2010); Krivanek et
al. (2014); Hage et al. (2018, 2020); Yan et al. (2021); García de Abajo and
Di Giulio (2020), as exemplified by recent studies of photon confinement in
optical cavities Kfir et al. (2020); Wang et al. (2020); Auad et al. (2022),
atomic vibrations in thin layers Hage et al. (2018, 2020); Yan et al. (2021)
and molecules Rez et al. (2016); Haiber and Crozier (2018); Hachtel et al.
(2019), and collective excitations such as phonon polaritons Krivanek et al.
(2014); Lagos et al. (2017); Li et al. (2020); Konečná et al. (2021) and
plasmons Bosman et al. (2007); Nelayah et al. (2007); Rossouw and Botton
(2013); Tan et al. (2014); Mkhitaryan et al. (2021).

At e-beam energies $>30\,$keV, typically employed in transmission electron
microscopes to perform EELS analyses, the per-electron excitation probability
of each individual mode in the specimen lies several orders of magnitude below
unity. While such weak interaction is beneficial to grant us clean access into
the nanoscale optical response over a wide spectral range
($10^{-3}-10^{3}\,$eV), a low excitation probability also implies that we
operate in the linear regime, which is useless to track the ultrafast dynamics
associated with a nonlinear behavior. This situation can be mitigated by
resorting to less energetic probes like those available in low-energy electron
microscopes Rocca (1995); Tromp (2019). Indeed, individual $\lesssim 100\,$eV
electrons are predicted to generate multiple excitations of a single optical
mode by appropriately adjusting the beam energy García de Abajo (2013), while
the onset of anharmonic response in this regime is expected to produce mode
saturation and spectral shifts Cox and García de Abajo (2020). In a different
approach, femtosecond resolution is achieved in ultrafast electron microscopy
by synchronizing laser and electron pulses in their arrival at the sampled
structure Barwick et al. (2009); García de Abajo et al. (2010); Feist et al.
(2015); Piazza et al. (2015); Kfir et al. (2020); Wang et al. (2020); Henke et
al. (2021), a method that potentially enables the determination of nonlinear
response functions with nanoscale resolution Konečná et al. (2019).

Many of the aforementioned studies focus on bosonic excitations (e.g., phonons
Krivanek et al. (2014); Lagos et al. (2017); Li et al. (2020); Konečná et al.
(2021) and plasmons Bosman et al. (2007); Nelayah et al. (2007); Rossouw and
Botton (2013); Tan et al. (2014); Mkhitaryan et al. (2021)), which exhibit the
characteristic linear response of harmonic oscillators, unless strong external
fields are introduced to drive them beyond the parabolic potential region. In
the opposite extreme, two- and few-level (fermionic) systems display a
paradigmatic nonlinear behavior, whereby a given excitation can block
subsequent ones. As an example, the discreteness of energy levels in
nanographenes permeates their optical response and enables nonlinear
interactions at the single-free-electron level Cox and García de Abajo (2020).
Nevertheless, fermionic excitations in systems such as atoms, molecules, and
defect states in solids possess a weak transition strength that is essentially
limited by the $f$-sum rule Noziéres (1997) and, therefore, demands the use of
low-energy electrons to yield measurable inelastic scattering signals.

Nonlinear effects open fundamental questions, such as whether an individual
electron can produce a given excitation with 100% probability, as well as the
role of the electron wave function in determining that probability. In
addition, we expect qualitatively different behavior between excitations of
bosonic and fermionic character in the nonlinear regime. Because the probe
energies required to reach a sizeable interaction strength are likely
comparable to the transition energies, recoil effects are also anticipated to
play an important role. These are relevant problems of the yet poorly explored
terrain of nonlinear and recoil phenomena taking place during the interaction
of free electrons with localized excitations.

In this Letter, we show that a free electron can excite a two-level system
with 100% probability, provided the transition symmetry and interaction
strength meet specific conditions. Based on a quantum description of free
electrons and localized excitations that rigorously incorporates nonlinear and
recoil effects, we show that the excitation probability is independent of the
electron wave function profile. Our calculations for bosonic and fermionic
systems also demonstrate that recoil effects are irrelevant unless the
electron energy is only a few times larger than the transition energy. Besides
their fundamental interest, our results suggest a way to control localized
excitations by means of free electrons, while establishing universal rules for
the maximum achievable probability depending on the symmetry of the excitation
and the electron-sample coupling strength.

Figure 1: Nonlinear effects in the excitation of two-level systems by a single
free electron. We represent the probability $P_{1}$ as a function of the
first-order (linear) probability $P_{1}^{\rm lin}$ for a lossless point-like
two-level system with and without inclusion of backscattering. The
$P_{1}=P_{1}^{\rm lin}$ probability for a bosonic excitation is also shown for
comparison. Many multiple uncorrelated electrons produce a probability of 1/2.

##  II Nonperturbative excitation probability including recoil

We consider a collimated e-beam focused down to a small lateral size at the
region of interaction with the sampled structure, such that we can ignore its
dynamics in a plane perpendicular to the beam direction $z$. We further assume
a nonlossy specimen comprising a discrete set of states $|j\rangle$ of
energies $\hbar\omega_{j}$ and initially placed in the the ground state
$|0\rangle$. The Hamiltonian of the combined electron-sample system can be
written

$\displaystyle\hat{\mathcal{H}}=$
$\displaystyle\hbar\int\\!dq\;\varepsilon_{q}\,|q\rangle\langle
q|+\hbar\sum_{j}\omega_{j}|j\rangle\langle j|$
$\displaystyle+\hbar\int\\!dq\int\\!dq^{\prime}\sum_{jj^{\prime}}G_{qj,q^{\prime}j^{\prime}}|qj\rangle\langle
q^{\prime}j^{\prime}|,$

where the electron is represented by orthonormal momentum states $|q\rangle$
of energies $\hbar\varepsilon_{q}$, whereas $G_{qj,q^{\prime}j^{\prime}}$ are
electron-sample coupling coefficients. Expanding the wave function of the
combined system as

$\displaystyle|\Psi(t)\rangle=\int\\!dq\sum_{j}{\rm e}^{-{\rm
i}(\varepsilon_{q}+\omega_{j})t}\alpha_{qj}(t)\,|qj\rangle,$

inserting it into the Schrödinder equation
$\hat{\mathcal{H}}|\Psi(t)\rangle={\rm i}\hbar|\dot{\Psi}(t)\rangle$, and
adopting the initial conditions
$\alpha_{qj}(-\infty)=\alpha_{q}^{0}\,\delta_{j0}$ [i.e., with the specimen in
the ground state $j=0$ and an incident electron wave function
$\psi^{0}(z,t)\propto\int\\!dq\,\alpha_{q}^{0}\,{\rm e}^{{\rm
i}(qz-\varepsilon_{q}t)}$], we find the post-interaction solution (see
Appendix)

$\displaystyle\alpha_{qj}(\infty)$
$\displaystyle=\alpha_{q}^{0}\,\delta_{j0}-2\pi{\rm
i}\frac{M_{q\tilde{q}_{j},j}}{v_{\tilde{q}_{j}}}\;\alpha_{\tilde{q}_{j}}^{0},$
(1)

where the coefficients $M_{qq^{\prime},j}$ are independent of the incident
electron state and satisfy to the self-consistent Lippmann–Schwinger Messiah
(1966) relation

$\displaystyle
M_{qq^{\prime},j}=G_{qj,q^{\prime}0}-\\!\int\\!dq^{\prime\prime}\sum_{j^{\prime}}\frac{G_{qj,q^{\prime\prime}j^{\prime}}M_{q^{\prime\prime}q^{\prime},j^{\prime}}}{\varepsilon_{q^{\prime\prime}q^{\prime}}+\omega_{j^{\prime}0}-{\rm
i}0^{+}}$ (2)

with $\varepsilon_{qq^{\prime}}=\varepsilon_{q}-\varepsilon_{q^{\prime}}$ and
$\omega_{jj^{\prime}}=\omega_{j}-\omega_{j^{\prime}}$. Here,
$v_{q}=d\varepsilon_{q}/dq$ is the group velocity of the $q$ electron
component, while $\tilde{q}_{j}$ is implicitly defined by
$\varepsilon_{\tilde{q}_{j}}=\varepsilon_{q}+\omega_{j0}$ with
$\tilde{q}_{j}>0$ (i.e., $\alpha_{\tilde{q}_{j}}^{0}$ only contains forward
propagating components).

We are interested in the probability $P_{j}$ for a sampled system initially
prepared in its ground state $|0\rangle$ to be left in state $|j\rangle$ after
the interaction has taken place. We thus write
$P_{j}=\int\\!dq\;\big{|}\alpha_{qj}(\infty)\big{|}^{2}$, which upon insertion
of Eq. (1), leads to a decomposition of the probability in incident-momentum
components according to (see Appendix)

$\displaystyle P_{j}=\int_{q_{\rm
min}^{j}}^{\infty}\\!dq\,\big{|}\alpha_{q}^{0}\big{|}^{2}\,P_{q,j},$ (3)

where

$\displaystyle
P_{q,j}=\frac{4\pi^{2}}{v_{q_{j}}v_{q}}\left(\left|M_{q_{j},q,j}\right|^{2}+\big{|}M_{-q_{j},q,j}\big{|}^{2}\right)$
(4)

for excited states $j\neq 0$. Here, the final electron wave vector $q_{j}>0$
is defined through $\varepsilon_{q_{j}}=\varepsilon_{q}-\omega_{j0}$, and a
minimun incident wave vector $q_{\rm min}^{j}$ is imposed by the threshold
excitation condition $\varepsilon_{q_{\rm min}^{j}}=\omega_{j0}$. The first
and second terms in Eq. (4) correspond to the contributions of forward and
backward electron scattering (i.e., final wave vectors $q_{j}$ and $-q_{j}$,
respectively). This result reveals a trivial role of the incident electron
wave function: each initial wave vector component contributes to the
excitation probability in proportion to $\big{|}\alpha_{q}^{0}\big{|}^{2}$
[see Eq. (9)], with no dependence on the phase of $\alpha_{q}^{0}$ [i.e., on
the profile of the incident wave function $\psi^{0}(z,t)$]. We remark that
this conclusion is derived from a nonperturbative formalism that rigorously
accounts for nonlinear and recoil effects.

Because the energy spread of the incident beam plays a trivial role, we limit
our discussion to monochromatic electrons of energy $\hbar\varepsilon_{q_{0}}$
with $\big{|}\alpha_{q}^{0}\big{|}^{2}=\delta(q-q_{0})$, so that the
excitation probability reduces to $P_{j}=P_{q_{0},j}$, subject to the
condition $q_{0}>q_{\rm min}^{j}$. In addition, we focus on two-level systems,
although the present formalism can be readily applied to multilevel
configurations. We thus concentrate on the excitation probability $P_{1}$ and
also consider the linear probability $P_{1}^{\rm lin}$ for reference, obtained
from Eqs. (2) and (4) by neglecting the integral term in the former (see
Appendix).

##  III Point-like interaction limit

As a first tutorial step, we obtain a closed-form solution when the
interaction is localized to just one point, so that the coupling coefficients
$G_{qj,q^{\prime}j^{\prime}}$ are independent of $q$ and $q^{\prime}$. Then,
the excitation probability reduces to (see Appendix)

$\displaystyle P_{1}=\frac{P_{1}^{\rm lin}}{\left(1+P_{1}^{\rm
lin}/2\right)^{2}},$

which presents a single maximum $P_{1}=1/2$ as a function of the linear
probability at $P_{1}^{\rm lin}=2$, as shown in Fig. 1. Only for this case, we
include backscattering in the linear probability. This result already reveals
that maximum excitation is only achieved for a very specific value of the
coupling coefficient or, alternatively, $P_{1}^{\rm lin}$. Interestingly, the
presence of two inelastic channels (forward and backward scattering) limits
the maximum probability to 50%. Indeed, if we disregard backscattering, which
should be reasonable for energetic electrons, a similar analysis leads to (see
Appendix)

$\displaystyle P_{1}=\frac{P_{1}^{\rm lin}}{\left(1+P_{1}^{\rm
lin}/4\right)^{2}},$

whose maximum value is now $P_{1}=1$, obtained at $P_{1}^{\rm lin}=4$. For
comparison, we show the $P_{1}=P_{1}^{\rm lin}$ line corresponding to a
bosonic mode, and obviously, all of these results are in mutual agreement in
the $P_{1}^{\rm lin}\ll 1$ limit. Incidentally, the average population of the
excited state in a two-level system interacting with many multiple
uncorrelated electrons is $1/2$ García de Abajo and Di Giulio (2020).

Figure 2: Excitation of two-level systems by a single free electron. We
consider the configuration illustrated in the inset of panel (a) and study the
excitation probability for different transition symmetries. (a),(b) Electron-
sample interaction coefficients $G_{10}(z)$ and $G_{q1,00}$ in the real-space
(a) and momentum-space (b) representations, respectively, for dipolar ($p_{x}$
and $p_{z}$) and quadrupolar ($d_{z^{2}}$, $d_{xz}$, and $d_{x^{2}-y^{2}}$)
excitations with different angular symmetries. (c),(d) Evolution of the
excited state occupation as a function of position along the electron
trajectory $z$ and impact parameter $R_{e}$ for two selected dipolar and
quadrupolar excitation symmetries (see labels). (e)-(i) Dependence of the
post-interaction excitation probability (at $z\rightarrow\infty$) as a
function of impact parameter $R_{e}$ and linear probability $P_{1}^{\rm lin}$
for all nonvanishing dipolar and quadrupolar excitation symmetries under the
investigated beam-sample configuration. The probability reaches 100% at the
positions indicated by the white dots in (e)-(i). We take fixed values of
$P_{1}^{\rm lin}=2.5$ and 1.5 in (c) and (d), respectively.

##  IV Nonlinear e-beam excitation without recoil

As we show below, recoil effects can be neglected if the electron energy
exceeds several times the transition energy. We can then linearize the
electron energy difference as
$\varepsilon_{qq^{\prime}}\approx(q-q^{\prime})v$, where $v$ is the electron
velocity. Considering a small sampled system, whose interaction with low-
energy electrons can be described through the Coulomb potential, we find the
associated coupling coefficients to only depend on the wave vector difference
$q-q^{\prime}$ and take the form $G_{qj,0j^{\prime}}\propto\big{(}{\rm
sign}\\{q\\}\big{)}^{\sigma}|q|^{l}K_{m}(|q|R_{e})$, where $R_{e}$ is the
beam-sample distance, $(l,m)$ are the angular momentum numbers associated with
the excitation symmetry, $\sigma$ takes values of 0 or 1, and a constant of
proportionality depending on the details of the system is taken to be absorbed
in $P_{1}^{\rm lin}$. In particular, we consider excitations of dipolar
[$p_{x}$ and $p_{z}$, corresponding to $(l,m,\sigma)=(1,1,0)$ and $(1,0,1)$,
respectively] and quadrupolar [$d_{z_{2}}$, $d_{xz}$, and $d_{x^{2}-y^{2}}$,
corresponding to $(2,0,0)$, $(2,1,1)$, and $(2,2,0)$] character, with a
geometrical configuration as shown in the inset of Fig. 2(a) (see Appendix for
details, and Fig. 2(b) for the associated momentum-space coupling
coefficients).

Under these conditions, the wave function of the system admits the form (see
Appendix)

$\displaystyle\langle z|\Psi(t)\rangle=\psi^{0}(z,t)\sum_{j}f_{j}(z)\,{\rm
e}^{-{\rm i}\omega_{j0}(z/v-t)}\,{\rm e}^{-{\rm i}\omega_{j}t}\,|j\rangle,$

where the space-dependent functions $f_{j}(z)$ evolve as

$\displaystyle\frac{d\,f_{j}(z)}{dz}$ $\displaystyle=-\frac{{\rm
i}}{v}\sum_{j^{\prime}}G_{jj^{\prime}}(z)\;{\rm e}^{{\rm
i}\omega_{jj^{\prime}}z/v}\,f_{j^{\prime}}(z),$ (5)

and we introduce real-space coupling coefficients
$G_{jj^{\prime}}(z)=\int\\!dq\;G_{qj,0j^{\prime}}\,{\rm e}^{{\rm i}qz}$ [see
Fig. 2(a)]. Finally, the excitation probability is simply given by
$P_{j}=\big{|}f_{j}(\infty)\big{|}^{2}$, while the linear limit reduces to
$P_{j}^{\rm lin}=(4\pi^{2}/v^{2})\big{|}G_{0j,\omega_{j0}/v,0}\big{|}^{2}$ for
$j\neq 0$ (see Appendix).

We numerically integrate Eq. (5) for two-level systems with the excitation
symmetries noted above to obtain the universal plots of $P_{1}$ presented in
Figs. 2(e)-2(i) as a function of the dimensionless parameters
$\omega_{10}R_{e}/v$ and $P_{1}^{\rm lin}$. Remarkably, $P_{1}$ reaches a
single maximum of 100% at a specific $(R_{e},P_{1}^{\rm lin})$ point [white
dots in Figs. 2(e)-2(i)]. The position of this maximum occurs at values of
$P_{1}^{\rm lin}$ that are in the range of those obtained in the point-
interaction limit (Fig. 1), while the impact parameter $R_{e}$ lies close to
the stationary points of $G_{q1,00}$ as a function of $qR_{e}$ for a wave
vector transfer $q=\omega_{10}/v$ determined by the nonrecoil approximation
[cf. the maxima of the curves in Fig. 2(b) and the abscissas of the white dots
in Figs. 2(e)-2(i)]. Two of the studied symmetries have this maximum at
$R_{e}=0$, accompanied by a lack of any zeros in the real-space profile of the
corresponding coupling coefficients [Fig. 2(a)], in contrast to the other
excitations under consideration. Incidentally, $P_{1}$ presents multiple
maxima as we move along $R_{e}$ for fixed $P_{1}^{\rm lin}$, the magnitudes of
which decrease with increasing impact parameter. This is the result of a
complex evolution of the position-dependent probability $|f_{1}(z)|^{2}$ along
the electron path, which exhibits oscillations before reaching an asymptotic
value of $P_{1}$ at large $z$ [see examples of this dynamics in Figs.
2(c),2(d)].

Figure 3: Recoil effects in near-edge excitation. (a) Excitation probability
$P_{1}$ for a two-level sample as a function of incident electron energy
(horizontal axis, normalized to the excitation energy $\hbar\omega_{10}$) and
linear excitation probability $P_{1}^{\rm lin}$. (b) Probability extracted
from panel (a) for fixed $P_{1}^{\rm lin}=1$ (black curve) compared with the
excitation probability for a bosonic mode with the same $0\rightarrow 1$
matrix element (green curves, comprising a decomposition in the contributions
$P_{n}$ of different Fock states $|n\rangle$, as well as the final average
population $\langle n\rangle=\sum_{n=1}^{\infty}n\,P_{n}$). Dashed curves
inticate the $\varepsilon_{0}\gg\omega_{10}$ nonrecoil limit. We consider
dipolar excitations of $p_{x}$ symmetry and a normalized impact parameter
$\omega_{10}R_{e}/v=0.2$ in all cases.

##  V Effect of electron recoil

The solution of Eqs. (2) and (4) for a two-level system produces an excitation
probability that is substantially reduced with respect to the nonrecoil limit
when the incident electron energy $\hbar\varepsilon_{0}$ approaches the
excitation threshold $\hbar\omega_{10}$. We illustrate this effect in Fig.
3(a) for $\omega_{10}R_{e}/v=0.2$ and $p_{x}$ transition symmetry over a wide
range of coupling strengths (vertical axis), but we find that this conclusion
is general upon extensive numerical inspection of different $R_{e}$ values.
The nonrecoil result is however recovered when $\varepsilon_{0}$ is just a few
times $\omega_{10}$. A similar effect of recoil is observed in the excitation
of a bosonic mode [Fig. 3(b)], although the interplay between different Fock
states $|n\rangle$ leads to a more complex evolution characterized by sharp
oscillations in both the total excitation probability and the partial
contribution coming from each $|n\rangle$ state. These oscillations are
attenuated as $\varepsilon_{0}$ increases, leading to a Poissonian
distribution García de Abajo (2013); Di Giulio et al. (2019) (see also
Appendix).

##  VI CONCLUSIONS

In conclusion, a wealth of phenomena unfolds from the interaction between free
electrons and few-level systems. In particular, we have shown that achieving
complete excitation of a single transition in a specimen by an individual free
electron is not simply a matter of increasing the interaction strength, but it
also requires a specific balance that depends on the symmetry of the excited
mode. In addition, the excitation probability is independent of the electron
wave function profile even when fully accounting for nonlinear and recoil
effects. Low-energy electrons in the $<100\,$eV range are promising to explore
these effects, as they can generate multiple excitations of a single plasmon
mode in atomically thin nanostructures García de Abajo (2013). Excitons in
two-dimensional materials Tizei et al. (2015) offer a potentially practical
candidate to study the iteration of free electrons with few-level systems,
while defect states in those materials, already explored with tunneling
microscopes Schuler et al. (2019), are robust two-level systems that could be
investigated with low-energy electrons in a reflection configuration. Free-
electron interaction with diluted atomic or molecular gases could also serve
as a platform to study the coupling strength, while incipient electron
microscopy studies on optical atomic lattices and condensates Gericke et al.
(2008) could be extended to measure inelastic scattering and explore the
physics portrayed in the present work.

## ACKNOWLEDGMENTS

This work has been supported in part by the European Research Council
(Advanced Grant 789104-eNANO), the European Commission (Horizon 2020 Grants
101017720 FET-Proactive EBEAM and 964591-SMART-electron), the Spanish MICINN
(PID2020-112625GB-I00 and Severo Ochoa CEX2019-000910-S), the Catalan CERCA
Program, and Fundaciós Cellex and Mir-Puig.

## APPENDIX

We provide a self-contained derivation of the theory used in the main text to
describe the interaction between a single collimated free electron and a
system comprising a discrete set of quantum states. The excitation
probabilities are shown to only depend on the spectral distribution of the
incident electron, but not on the phase and shape of its wave function, even
when recoil, nonlinear, and relativistic effects are accounted for in a
rigorous manner. We also present explicit expressions for the coupling
coefficients between the electron and selected excitations with well-defined
multipolar symmetry in the nonrelativistic limit. A simpler solution to the
interaction problem is further elaborated within the nonrecoil approximation.
Finally, we offer details of a numerical implementation using nonrelativistic
kinematics, but fully incorporating recoil and nonlinear effects.

##  Appendix A Free-electron interaction with a discrete-level system beyond
the linear nonrecoil regime

We study the interaction between a collimated free electron and a sample
comprising a discrete set of states $|j\rangle$ of energies $\hbar\omega_{j}$
by adopting the following assumptions:

  1. (i) —

Longitudinal motion.—The electron is tightly focused in the transverse plane
(perpendicular to the propagation direction $z$) down to a small region in
which the interaction with the sample is approximately homogeneous [i.e.,
independent of the transverse coordinates ${\bf R}=(x,y)$ across the electron
beam (e-beam)]. In addition, the transverse electron wave function remains
nearly unchanged during the interaction time, so that we can dismiss it as
well as any change in the transverse electron energy.

  2. (ii) —

Nonlossy specimen.—Inelastic decay of the excited sample states plays a
negligible role during the interaction time.

Item (i) allows us to describe the electron in terms of a basis set of
momentum states $|q\rangle$ of energies $\hbar\varepsilon_{q}$ (with
$\varepsilon_{q}=c\sqrt{(m_{\rm e}c/\hbar)^{2}+q^{2}}$ or
$\varepsilon_{q}=\hbar q^{2}/2m_{\rm e}$ within relativistic or
nonrelativistic kinematics, respectively) and wave functions $\langle
z|q\rangle={\rm e}^{{\rm i}qz}/\sqrt{2\pi}$ satisfying the orthonormalization
relation $\langle q|q^{\prime}\rangle=\delta(q-q^{\prime})$. In addition,
point (ii) permits describing the evolution of the system by solving the
Schrödinger equation with the total Hamiltonian

$\displaystyle\hat{\mathcal{H}}=\hbar\int\\!dq\;\varepsilon_{q}\,|q\rangle\langle
q|+\hbar\sum_{j}\omega_{j}|j\rangle\langle j|$
$\displaystyle+\hbar\int\\!dq\int\\!dq^{\prime}\sum_{jj^{\prime}}G_{qj,q^{\prime}j^{\prime}}|qj\rangle\langle
q^{\prime}j^{\prime}|,$ (6)

where $G_{qj,q^{\prime}j^{\prime}}=\hbar^{-1}\langle
qj|\hat{\mathcal{H}}|q^{\prime}j^{\prime}\rangle$ are electron-sample coupling
coefficients.

We proceed by writing the wave function of the combined electron-sample system
as

$\displaystyle|\Psi(t)\rangle=\int\\!dq\sum_{j}{\rm e}^{-{\rm
i}(\varepsilon_{q}+\omega_{j})t}\alpha_{qj}(t)\,|qj\rangle,$ (7)

which, upon insertion into $\hat{\mathcal{H}}|\Psi(t)\rangle={\rm
i}\hbar|\dot{\Psi}(t)\rangle$, leads to the equation of motion

$\displaystyle\dot{\alpha}_{qj}(t)$ $\displaystyle=-{\rm
i}\int\\!dq^{\prime}\sum_{j^{\prime}}G_{qj,q^{\prime}j^{\prime}}\,\alpha_{q^{\prime}j^{\prime}}(t)\,{\rm
e}^{{\rm i}(\varepsilon_{qq^{\prime}}+\omega_{jj^{\prime}})t}$ (8)

for the expansion coefficients in Eq. (7), where we use the compact notation
$\varepsilon_{qq^{\prime}}=\varepsilon_{q}-\varepsilon_{q^{\prime}}$ and
$\omega_{jj^{\prime}}=\omega_{j}-\omega_{j^{\prime}}$. Obviously, the
Hamiltonian in Eq. (6) is Hermitian, so the normalization condition
$\int\\!dq\sum_{j}|\alpha_{qj}|^{2}=1$ is maintained during time propagation.

We are interested in studying the probability

$\displaystyle P_{j}=\int\\!dq\;\big{|}\alpha_{qj}(\infty)\big{|}^{2}$ (9)

that a sample initially prepared in its ground state $|0\rangle$ is left in a
state $|j\rangle$ after the interaction has taken place. Consequently, we set
$\alpha_{qj}(-\infty)=\alpha_{q}^{0}\,\delta_{j0}$, where the coefficients
$\alpha_{q}^{0}$ define the incident electron wave function

$\displaystyle\psi^{0}(z,t)=\int\\!dq\;\alpha_{q}^{0}\;\frac{{\rm e}^{{\rm
i}(qz-\varepsilon_{q}t)}}{\sqrt{2\pi}},$ (10)

whose normalization ($\int dz\;\big{|}\psi^{0}(z,t)\big{|}^{2}=1$) imposes the
condition $\int\\!dq\;\big{|}\alpha_{q}^{0}\big{|}^{2}=1$.

At this point, we anticipate a solution of the form

$\displaystyle\alpha_{qj}(t)=\alpha_{q}^{0}\,\delta_{j0}-\int\\!dq^{\prime}\;M_{qq^{\prime},j}\;\alpha_{q^{\prime}}^{0}\;\frac{{\rm
e}^{{\rm
i}(\varepsilon_{qq^{\prime}}+\omega_{j0})t}}{\varepsilon_{qq^{\prime}}+\omega_{j0}-{\rm
i}0^{+}},$ (11)

which is suggested by iteratively integrating Eq. (8) in a perturbation-theory
approach. Indeed, inserting Eq. (11) into Eq. (8) and integrating over time,
we find the self-consistent Lippmann–Schwinger equation Messiah (1966)

$\displaystyle
M_{qq^{\prime},j}=G_{qj,q^{\prime}0}-\int\\!dq^{\prime\prime}\sum_{j^{\prime}}\frac{G_{qj,q^{\prime\prime}j^{\prime}}M_{q^{\prime\prime}q^{\prime},j^{\prime}}}{\varepsilon_{q^{\prime\prime}q^{\prime}}+\omega_{j^{\prime}0}-{\rm
i}0^{+}}$ (12)

for the coefficients in Eq. (11), which are independent of the incident
electron state. Now, applying the identity $\exp({\rm i}\theta t)/(\theta-{\rm
i}0^{+})\xrightarrow[t\to\infty]{}2\pi{\rm i}\delta(\theta)$ with
$\theta=\varepsilon_{qq^{\prime}}+\omega_{j0}$ to Eq. (11), we find the post-
interaction solution

$\displaystyle\alpha_{qj}(\infty)$
$\displaystyle=\alpha_{q}^{0}\,\delta_{j0}-2\pi{\rm
i}\int\\!dq^{\prime}\;M_{qq^{\prime},j}\;\alpha_{q^{\prime}}^{0}\;\delta(\varepsilon_{qq^{\prime}}+\omega_{j0})$
$\displaystyle=\alpha_{q}^{0}\,\delta_{j0}-2\pi{\rm
i}\frac{M_{q\tilde{q}_{j},j}}{v_{\tilde{q}_{j}}}\;\alpha_{\tilde{q}_{j}}^{0},$

where $v_{q}=d\varepsilon_{q}/dq$ is the group velocity of the $q$ electron
component, and the rightmost expression is obtained by manipulating the
$\delta$-function as
$\delta(\varepsilon_{qq^{\prime}}+\omega_{j0})=\big{[}\delta(q^{\prime}-\tilde{q}_{j})+\delta(q^{\prime}+\tilde{q}_{j})\big{]}/v_{\tilde{q}_{j}}$
with $\tilde{q}_{j}$ implicitly defined by
$\varepsilon_{\tilde{q}_{j}}=\varepsilon_{q}+\omega_{j0}$. Here, we consider
an incident electron that only contains $\tilde{q}_{j}>0$ components in
$\alpha_{\tilde{q}_{j}}^{0}$ (i.e., it moves toward increasing $z$), but the
final state can receive both $q>0$ (forward scattering) and $q<0$
(backscattering) contributions. Finally, inserting this result into Eq. (9),
using the identity $v_{\tilde{q}_{j}}d\tilde{q}_{j}=v_{q}dq$, and changing the
variable of integration as $\tilde{q}_{j}\rightarrow q$, the post-interaction
occupation probability of level $j$ reduces to

$\displaystyle P_{j}$ $\displaystyle=\int_{q_{\rm
min}^{j}}^{\infty}\\!\frac{dq}{v_{q_{j}}v_{q}}$ (13)
$\displaystyle\times\left[\left|\left(\delta_{j0}\,v_{q}-2\pi{\rm
i}\,M_{q_{j},q,j}\right)\right|^{2}+4\pi^{2}\,\big{|}M_{-q_{j},q,j}\big{|}^{2}\right]\;\big{|}\alpha_{q}^{0}\big{|}^{2},$

where $q_{j}$ depends on $q$ and is now defined by
$\varepsilon_{q_{j}}=\varepsilon_{q}-\omega_{j0}$ (i.e., $q_{j}=(m_{\rm
e}c/\hbar)\big{[}\big{(}\sqrt{1+(\hbar q/m_{\rm
e}c)^{2}}-\hbar\omega_{j0}/m_{\rm e}c^{2}\big{)}^{2}-1\big{]}^{1/2}$ or
$q_{j}=\sqrt{q^{2}-2m_{\rm e}\omega_{j0}/\hbar}$ within relativisitic or
nonrelativistic kinematics, respectively); a wave vector threshold $q_{\rm
min}^{j}$ is imposed by the minimum electron energy capable of exciting the
$j$ level in the sample (i.e., $\varepsilon_{q_{\rm min}^{j}}=\omega_{j0}$);
and the first and second terms inside the square brackets account for forward
and backward electron scattering contributions with final electron wave
vectors $q_{j}$ and $-q_{j}$, respectively. The result embodied in Eq. (13) is
general within the approximations in points (i) and (ii) above, and it shows a
trivial dependence on the incident electron wave function: the contribution of
each incident wave vector component $q$ to the excitation probability is
weighted by its strength $\big{|}\alpha_{q}^{0}\big{|}^{2}$, so the phase in
$\alpha_{q}^{0}$ (i.e., the wave function profile) does not play any role at
all. We stress that the present derivation demonstrates that this conclusion
is maintained even when rigorously accounting for nonlinear electron-sample
interactions and recoil effects.

###  A.1 Solution for monochromatic electrons

Because the energy spread of the incident electron plays a trivial role, we
can limit our discussion to monochromatic electrons of incident energy
$\hbar\varepsilon_{q_{0}}$, characterized by
$\big{|}\alpha_{q}^{0}\big{|}^{2}=\delta(q-q_{0})$. Then, the excitation
probability reduces to

$\displaystyle P_{j}=$
$\displaystyle\frac{\Theta(\epsilon_{q_{0}}-\omega_{j0})}{v_{q_{j}}v_{q_{0}}}$
(14) $\displaystyle\times\left[\left|\delta_{j0}\,v_{q_{0}}-2\pi{\rm
i}\,M_{q_{j},j}\right|^{2}+4\pi^{2}\,\big{|}M_{-q_{j},j}\big{|}^{2}\right],$

where the coefficients $M_{qj}\equiv M_{qq_{0},j}$ satisfy the self-consistent
equation

$\displaystyle
M_{qj}=G_{qj,q_{0}0}-\int\\!dq^{\prime}\sum_{j^{\prime}}\frac{G_{qj,q^{\prime}j^{\prime}}M_{q^{\prime}j^{\prime}}}{\varepsilon_{q^{\prime}q_{0}}+\omega_{j^{\prime}0}-{\rm
i}0^{+}},$ (15)

which we write directly from Eq. (12). We remark again that the first and
second terms inside the square brackets of Eq. (14) correspond to forward and
backward electron scattering contributions, with associated final electron
wave vectors determined as illustrated in the following sketch:

Obviously, backscattering is generally negligible, unless the electron energy
is close to the threshold $\hbar\omega_{j0}$. For reference, we also consider
the first-order (linear) excitation probabilities

$\displaystyle P_{j}^{\rm
lin}=\frac{4\pi^{2}}{v_{q_{0}}v_{q_{j}}}\left(\left|G_{q_{j}j,q_{0}0}\right|^{2}+\left|G_{-q_{j}j,q_{0}0}\right|^{2}\right),$
(16)

derived from Eq. (14) for $j\neq 0$ by neglecting the integral term in Eq.
(15) (i.e., setting $M_{qj}=G_{qj,q_{0}0}$).

###  A.2 Point-like excitation limit with backscattering

Interestingly, exact solutions can be obtained when the interaction is
localized to just a single point at $z=0$. Indeed, a point-like interaction is
characterized by coupling coefficients $G_{qj,q^{\prime}j^{\prime}}\equiv
G_{jj^{\prime}}$ that are independent of $q$ and $q^{\prime}$ [i.e., giving
rise to a real-space interaction coefficients of the form
$G_{jj}(z)\propto\delta(z)$ according to Eq. (26) in Sec. B below]. In this
scenario, further adopting nonrelativisitc kinematics, so that
$\varepsilon_{q}$ is quadratic in $q$, we find that Eq. (15) reduces to
$M_{j}=G_{j0}-2\pi{\rm
i}\sum_{j^{\prime}}G_{jj^{\prime}}\,M_{j^{\prime}}/v_{q_{j^{\prime}}}$ not
(a), whose solutions $M_{j}\equiv M_{qj}$ are also independent of $q$. In
particular, for a two-level system comprising the states $j=0$ and 1, the only
nonzero coupling elements are $G_{10}=G^{*}_{01}$, so combining the solution
for $M_{1}$ with Eq. (14), we find the excitation probability

$\displaystyle P_{1}=\frac{P_{1}^{\rm lin}}{\left(1+P_{1}^{\rm
lin}/2\right)^{2}}.$

We note that forward and backward scattering events contribute identically to
the linear probability under the assumption of a $q$-independent coupling,
which leads to $P_{1}^{\rm
lin}=8\pi^{2}\big{|}G_{10}\big{|}^{2}/v_{q_{0}}v_{q_{1}}$ [see Eq. (16)].
Then, a single absolute maximum $P_{1}=1/2$ is obtained for $P_{1}^{\rm
lin}=2$. Reassuringly, this approximation maintains the overall probability
$P_{0}+P_{1}=1$.

###  A.3 Point-like excitation limit without backscattering

When considering energetic electrons, a point-like scatterer is still
described by $q$-independent coupling coefficients for relatively small
changes in wave vector. However, backscattering (i.e., with $q<0$ and
$q^{\prime}>0$) involves large momentum transfers, for which
$G_{qj,q^{\prime}j^{\prime}}$ becomes negligible at large incident electron
energy. Then, we obtain $q$-independent coefficients $M_{j}\equiv
M_{qj}=G_{j0}-{\rm
i}\pi\sum_{j^{\prime}}G_{jj^{\prime}}\,M_{j^{\prime}}/v_{q_{j^{\prime}}}$ for
$q>0$ (notice the absence of a factor of 2 with respect to the result
including backscattering), while $M_{qj}=0$ for $q<0$. Finally, we find the
solution

$\displaystyle P_{1}=\frac{P_{1}^{\rm lin}}{\left(1+P_{1}^{\rm
lin}/4\right)^{2}}$

with $P_{1}^{\rm lin}=4\pi^{2}\big{|}G_{10}\big{|}^{2}/v_{q_{0}}v_{q_{1}}$
[i.e., also neglecting backscattering in Eq. (16)], which now presents a
single absolute maximum $P_{1}=1$ when $P_{1}^{\rm lin}=4$. In addition, the
condition $P_{0}+P_{1}=1$ is also verified.

##  Appendix B Numerical implementation

The coupling coefficients $G_{qj,q^{\prime}j^{\prime}}$ are generally smooth
functions of $q$ and $q^{\prime}$, so we can solve Eq. (15) by discretizing
$q$ through a set of $2N+1$ points $p_{l}=q_{0}+lh$, each of them representing
an interval $p_{l}-h/2<q<p_{l}+h/2$ labeled by $l=-N,\cdots,N$. Consecutive
points are separated by a constant spacing $h=\Delta/(2N+1)$, and together
they cover a finite range of size $\Delta$ centered at the incident electron
wave vector $q_{0}$. We can then recast Eq. (15) into the finite matrix
equation

$\displaystyle M_{j}=g_{j}-\sum_{j^{\prime}}S_{jj^{\prime}}\cdot
M_{j^{\prime}},$ (17)

where $M_{j}$ and $g_{j}$ denote vectors of $2N+1$ components $M_{j,l}\equiv
M_{q_{l}j}$ and $g_{j,l}\equiv G_{p_{l}j,q_{0}0}$, respectively, whereas
$S_{jj^{\prime}}=G_{jj^{\prime}}\cdot\Delta_{j^{\prime}}$ are
$(2N+1)\times(2N+1)$ square matrices defined in terms of the coupling matrices
$G_{jj^{\prime},ll^{\prime}}\equiv G_{p_{l}j,p_{l^{\prime}}j^{\prime}}$ and
the diagonal matrix

$\displaystyle\Delta_{j,ll^{\prime}}=\delta_{ll^{\prime}}\int_{p_{l}-h/2}^{p_{l}+h/2}\frac{dq}{\varepsilon_{qq_{0}}+\omega_{j0}-{\rm
i}0^{+}}.$

In practice, a strong electron-sample interaction is likely involving
relatively small electron velocities, and therefore, we can work in the
nonrelativistic limit and write $\varepsilon_{q}=\hbar q^{2}/2m_{\rm e}$ and
$v_{q}=\hbar q/m_{\rm e}$. The matrix elements of $\Delta_{j}$ are then given
by the closed-form expressions

$\displaystyle\Delta_{j,ll^{\prime}}=\delta_{ll^{\prime}}\;\frac{m_{\rm
e}}{\hbar|q_{j}|}\times\left\\{\begin{array}[]{ll}\log\left|\big{[}p_{l}^{2}-(|q_{j}|-h/2)^{2}\big{]}\big{/}\big{[}p_{l}^{2}-(|q_{j}|+h/2)^{2}\big{]}\right|+{\rm
i}\pi\,\theta_{jl},&\quad\quad\varepsilon_{q_{0}}>\omega_{j0},\\\ &\\\
2\left(\tan^{-1}\big{[}(p_{l}+h/2)/|q_{j}|\big{]}-\tan^{-1}\big{[}(p_{l}-h/2)/|q_{j}|\big{]}\right),&\quad\quad\varepsilon_{q_{0}}<\omega_{j0},\end{array}\right.$
(21)

where $|q_{j}|=\sqrt{|q_{0}^{2}-2m_{\rm e}\omega_{j0}/\hbar|}$,
$\theta_{jl}=1$ if either $p_{l}-h/2<|q_{j}|<p_{l}+h/2$ or
$p_{l}-h/2<-|q_{j}|<p_{l}+h/2$, and $\theta_{jl}=0$ otherwise. In addition, we
use the explicit expressions offered in Sec. D for the coupling coefficients
$G_{qj,q^{\prime}j^{\prime}}$ associated with selected excitations of well-
defined multipolar symmetries. Incidentally, $\Delta_{j,ll^{\prime}}$ may
diverge in the unlikely event that $q_{j\neq 0}$ coincides with $p_{l}\pm h/2$
for any $l$, a problem that we avoid by slightly changing $N$.

We apply this procedure to study two different types of samples in the main
text, for which the evaluation of the coupling matrices $G_{jj^{\prime}}$ and
the numerical solution of Eq. (17) are simplified by the following
considerations:

  * •

Two-level sample. For a system comprising two states $j=0$ and 1, we have
$G_{10}=G_{01}^{\dagger}$ and $G_{00}=G_{11}=0$. The solution of Eq. (17) for
the excited state vector is $M_{1}=\big{(}1-S_{10}\cdot
S_{01}\big{)}^{-1}\\!\\!\cdot g_{1}$ and requires performing one matrix
inversion and one product between square matrices, each of them involving
$\sim N^{3}$ complex-number multiplications.

  * •

Single-boson sample. Coupling to a bosonic mode of frequency $\omega_{b}$ can
also be described using the present formalism. The bosonic Fock states
$|j\rangle$ appear at equally spaced frequencies given by
$\omega_{j0}=j\omega_{b}$, whereas the coupling matrix elements only connect
consecutive states, and consequently, we have

$\displaystyle
G_{jj^{\prime}}=\left\\{\begin{array}[]{rl}0,&\quad\quad|j-j^{\prime}|\neq
1,\\\ \sqrt{j}\,G_{10},&\quad\quad j^{\prime}=j-1,\\\
\sqrt{j+1}\,G_{10}^{\dagger},&\quad\quad j^{\prime}=j+1.\end{array}\right.$
(25)

These matrices are expressed in terms of $G_{10}$, which describes
$0\rightarrow 1$ transitions. Then, Eq. (17) becomes a block-tridiagonal
system, which we solve following the standard forward-sweep, backward-
substitution method, requiring a total of $\sim(n+1)\times N^{3}$ complex-
number multiplications when the boson ladder is cut at $j=n$.

Incidentally, for these two types of samples, assuming the system to be
initially prepared in the ground state $j=0$ and considering a given final
state $j=j_{1}$, all of the terms in the scattering series obtained from Eq.
(17) by Taylor expanding the self-consistent equation $M=(1+S)^{-1}\cdot
g=g-S\cdot g+S^{2}\cdot g+\cdots$ present the same difference between the net
number of up ($j\rightarrow j+1$) and down ($j\rightarrow j-1$) jumps, so the
overall probabilities are left unchanged when $G_{10}$ is multiplied by an
arbitrary phase factor ${\rm e}^{{\rm i}\varphi}$ (i.e., $G_{jj^{\prime}}$
introduces a factor ${\rm e}^{-{\rm i}\varphi}$ or ${\rm e}^{{\rm i}\varphi}$
for up or down jumps, respectively, and the final amplitude is modified by an
overall amplitude ${\rm e}^{-{\rm i}j_{1}\varphi}$).

Using the methods discussed above, we find convergent results within the scale
of the plots in the main text using $2N+1\sim 500$ discretization points and a
wave vector range $\Delta$ of a few times
$q_{0}\big{(}1-\sqrt{1-\omega_{10}/\epsilon_{q_{0}}}\big{)}$. In addition, for
sufficiently large electron energy $\varepsilon_{q_{0}}\gg\omega_{10}$, this
procedure produces results in excellent quantitative agreement with the
solution found in the nonrecoil approximation [see Eq. (27) in Sec. C below].

##  Appendix C Nonrecoil approximation

A simpler solution is obtained if the energy spread of the incident e-beam and
the sample excitation energies are both small compared with the average
electron energy. We can then approximate
$\varepsilon_{q}-\varepsilon_{q^{\prime}}\approx(q-q^{\prime})v$, where
$v=\sum_{q}v_{q}\big{|}\alpha_{q}^{0}\big{|}^{2}$ is the average electron
velocity. In addition, we use the fact that the coefficients
$G_{qj,q^{\prime}j^{\prime}}$ only depend on the difference of wave vectors
$q-q^{\prime}$ not (b). Applying these considerations to Eq. (7), we can write

$\displaystyle\langle z|\Psi(t)\rangle=\frac{{\rm e}^{{\rm i}\bar{q}z-{\rm
i}\varepsilon_{\bar{q}}t}}{\sqrt{2\pi}}\sum_{j}\varphi_{j}(z,t)\,{\rm
e}^{-{\rm i}\omega_{j}z/v}\,|j\rangle,$

where $\bar{q}$ is the average wave vector, and we define the electron wave
function

$\displaystyle\varphi_{j}(z,t)$ $\displaystyle=\int\\!dq\;{\rm e}^{{\rm
i}(q-\bar{q}+\omega_{j}/v)(z-vt)}\alpha_{qj}(t)$

associated with the each sample state $j$. Likewise, multiplying both sides of
Eq. (8) by ${\rm e}^{{\rm i}(q-\bar{q}+\omega_{j}/v)(z-vt)}$ and integrating
over $q$, we obtain

$\displaystyle(\partial_{t}+v\partial_{z})\,\varphi_{j}(z,t)$
$\displaystyle=-{\rm i}\sum_{j^{\prime}}G_{jj^{\prime}}(z)\,{\rm e}^{{\rm
i}\omega_{jj^{\prime}}z/v}\,\varphi_{j^{\prime}}(z,t),$

where the coupling coefficient

$\displaystyle G_{jj^{\prime}}(z)=\int\\!dq\;G_{qj,0j^{\prime}}\,{\rm e}^{{\rm
i}qz}$ (26)

is expressed in the real-space representation.

These equations admit a solution of the form
$\varphi_{j}(z,t)=\phi^{0}(z-vt)\,{\rm e}^{{\rm
i}\omega_{0}(z/v-t)}\,f_{j}(z)$, consisting of an overall factor than only
depends on $z-vt$, accompanied by space-dependent functions $f_{j}(z)$ that
satify

$\displaystyle\frac{d\,f_{j}(z)}{dz}$ $\displaystyle=-\frac{{\rm
i}}{v}\sum_{j^{\prime}}G_{jj^{\prime}}(z)\;{\rm e}^{{\rm
i}\omega_{jj^{\prime}}z/v}\,f_{j^{\prime}}(z).$ (27)

Incidentally, we have explicitly indicated a factor ${\rm e}^{{\rm
i}\omega_{0}(z/v-t)}$ in $\varphi_{j}(z,t)$ to simplify the initial condition
for a sample prepared in the $|0\rangle$ state (see below). Also, we can set
$\phi^{0}(z-vt)=\int\\!dq\;\alpha_{q}^{0}\;{\rm e}^{{\rm
i}(q-\bar{q})(z-vt)}$, as suggested by the incident electron wave function in
Eq. (10), which takes the form $\psi^{0}(z,t)=\big{(}{\rm e}^{{\rm
i}\bar{q}z-{\rm i}\varepsilon_{\bar{q}}t}/\sqrt{2\pi}\big{)}\phi^{0}(z-vt)$ in
the nonrecoil approximation. Therefore, the combined electron-sample wave
function reduces to

$\displaystyle\langle z|\Psi(t)\rangle=\psi^{0}(z,t)\sum_{j}f_{j}(z)\,{\rm
e}^{-{\rm i}\omega_{j0}(z/v-t)}\,{\rm e}^{-{\rm i}\omega_{j}t}\,|j\rangle,$
(28)

which is a superposition of time-dependent states ${\rm e}^{-{\rm
i}\omega_{j}t}\,|j\rangle$, each of them accompanied by an electron wave
function $\psi^{0}(z,t)\,{\rm e}^{-{\rm i}\omega_{j0}(z/v-t)}$ that reflects
the associated change in electron momentum by $-\hbar\omega_{j0}/v$. Once
$f_{j}(z)$ is obtained by solving Eq. (27) with the initial conditions
$f_{j}(-\infty)=\delta_{j0}$ (i.e., the sample in the $|0\rangle$ state), the
post-interaction probabilities are finally given by
$P_{j}=\big{|}f_{j}(\infty)\big{|}^{2}$.

This analysis can be readily applied to any initial pure state $|j_{i}\rangle$
by simply substituting $\omega_{0}$ by $\omega_{j_{i}}$ in these expressions,
and also extended to start with a coherent superposition of sample states
$\sum_{j_{i}}a_{j_{i}}{\rm e}^{-{\rm i}\omega_{j_{i}}t}|j_{i}\rangle$ by
separately propagating each initial state $|j_{i}\rangle$ and then weighting
the resulting $f_{j}(z)$ with the coefficients $a_{j_{i}}$, so that we obtain
the wave function

$\displaystyle\langle
z|\Psi(t)\rangle=\psi^{0}(z,t)\sum_{jj_{i}}F_{jj_{i}}(z)\,a_{j_{i}}\,{\rm
e}^{-{\rm i}\omega_{jj_{i}}(z/v-t)}\,{\rm e}^{-{\rm
i}\omega_{j}t}\,|j\rangle,$

where $F_{jj_{i}}(z)$ is given by the coefficient $f_{j}(z)$ obtained from Eq.
(27) with the initial conditions $f_{j}(-\infty)=\delta_{jj_{i}}$. For
incident nonochromatic electrons, the probability of finding the sample in
state $j$ after the interaction has taken place becomes
$P_{j}=\sum_{j_{i}}\big{|}F_{jj_{i}}(\infty)\,a_{j_{i}}\big{|}^{2}$, whereas
the electron energy-loss probability reduces to $\Gamma_{\rm
EELS}(\omega)=\sum_{jj_{i}}\big{|}F_{jj_{i}}(\infty)\,a_{j_{i}}\big{|}^{2}\delta(\omega-\omega_{jj_{i}})$,
which in practice needs to be convoluted with the zero-loss peak of the
microscope.

Again, we corroborate the trivial role played by the incident electron wave
function: in the nonrecoil approximation, the final total wave function is the
incident electron wave function multiplied by factor [the $j$ sum in Eq. (28)]
that only depends on the incident electron state through the velocity $v$.
From a quantum-optics perspective Di Giulio et al. (2019), this result
reflects the fact that the electron acts on the sample as a classical source
[i.e., regarding Eq. (27) as the optical Bloch equations of the system] if its
velocity is taken to be constant (i.e., in the nonrecoil approximation), so
that starting with $j=0$, and in virtue of energy conservation, the
interaction simply causes the electron wave function to undergo rigid shifts
$\approx-\hbar\omega_{j0}/v$ in momentum when the sample is excited to a state
$j\neq 0$.

For reference, the linear probability for $j\neq 0$ now corresponds to the
solution of Eq. (27) obtained by replacing $f_{j^{\prime}}(z)$ by
$\delta_{j^{\prime}0}$ in the right-hand side, and consequently, we find
$f_{j}(\infty)=(-{\rm i}/v)\int dz\,G_{j0}(z)\;{\rm e}^{{\rm
i}\omega_{j0}z/v}=(-2\pi{\rm i}/v)\,G_{q_{j}j,q_{0}0}$ with
$q_{j}=q_{0}-\omega_{j0}/v$ [i.e., the inverse Fourier transform of Eq. (26)],
which leads to $P_{j}^{\rm
lin}=(4\pi^{2}/v^{2})\big{|}G_{q_{j}j,q_{0}0}\big{|}^{2}=(4\pi^{2}/v^{2})\big{|}G_{0j,\omega_{j0}/v,0}\big{|}^{2}$.
This result agrees with Eq. (16) when neglecting the backscattering term
$\propto\big{|}G_{-q_{j}j,q_{0}0}\big{|}^{2}$ and setting $v_{q}\approx v$.

###  C.1 Interaction with a bosonic mode

The solution to Eq. (27) becomes analytical for a sample hosting a single
boson mode of frequency $\omega_{b}$, which comprises an infinite number of
discrete Fock states $|j\rangle$ of frequencies $j\omega_{b}$ with
$j=0,1,\dots$, coupled by coefficients $G_{jj^{\prime}}(z)$ that satisfy Eq.
(25). Indeed, for such bosonic system, Eq. (27) reduces to

$\displaystyle\frac{d\,f_{j}(z)}{dz}$
$\displaystyle=\sqrt{j}\,u^{*}(z)\,f_{j-1}(z)-\sqrt{j+1}\,u(z)\,f_{j+1}(z)$

with $u(z)=({\rm i}/v)G_{10}^{*}(z){\rm e}^{-{\rm i}\omega_{b}z/v}$, so it
admits the closed-form solution Di Giulio et al. (2019) $f_{j}(z)={\rm
e}^{{\rm i}\chi(z)}{\rm
e}^{-|\beta_{0}(z)|^{2}/2}\big{[}\beta_{0}^{*}(z)\big{]}^{j}/\sqrt{j!}$, where
$\chi(z)=\int_{-\infty}^{z}dz^{\prime}\int_{-\infty}^{z^{\prime}}dz^{\prime\prime}\,{\rm
Im}\\{u^{*}(z^{\prime})u(z^{\prime\prime})\\}$ is a global phase and
$\beta_{0}(z)=\int_{-\infty}^{z}dz^{\prime}\,u(z^{\prime})$. The boson is then
evolving as a coherent state of varying amplitude $\beta_{0}(z)$, featuring a
Possonian distribution of occupation numbers $P_{j}(z)=|f_{j}(z)|^{2}={\rm
e}^{-|\beta_{0}(z)|^{2}}\big{|}\beta_{0}(z)\big{|}^{2j}/j!$ and average
population $\sum_{j=1}^{\infty}j\,P_{j}(z)=|\beta_{0}(z)|^{2}$ that leads to
the post-interaction value

$\displaystyle\langle j\rangle$
$\displaystyle=\sum_{j=0}^{\infty}j\,P_{j}=|\beta_{0}(\infty)|^{2}$
$\displaystyle=\frac{1}{v^{2}}\left|\int dz\,G_{10}(z){\rm e}^{{\rm
i}\omega_{b}z/v}\right|^{2}$
$\displaystyle=\frac{4\pi^{2}}{v^{2}}\big{|}G_{01,\omega_{b}/v,0}\big{|}^{2}=P_{1}^{\rm
lin},$

where the integral is again identified with the inverse Fourier transform
connecting $G_{10}(z)$ to $G_{q1,00}=G_{01,-q,0}$ [see Eq. (26)].

##  Appendix D Nonretarded coupling coefficients for small samples

In the nonretarded limit, the electron-sample interaction is mediated by the
Coulomb potential, so the corresponding matrix elements in Eq. (6) read

$\displaystyle G_{qj,q^{\prime}j^{\prime}}$
$\displaystyle=\frac{1}{\hbar}\langle
qj|\hat{\mathcal{H}}|q^{\prime}j^{\prime}\rangle$
$\displaystyle=-\frac{e}{2\pi\hbar}\int dz_{e}\;{\rm e}^{{\rm
i}(q^{\prime}-q)z_{e}}\int d^{3}{\bf r}\frac{\langle j|\hat{\rho}({\bf
r})|j^{\prime}\rangle}{|{\bf r}_{e}-{\bf r}|}$
$\displaystyle=-\frac{e}{\pi\hbar}\int d^{3}{\bf r}\;\langle j|\hat{\rho}({\bf
r})|j^{\prime}\rangle$ (29) $\displaystyle\quad\quad\quad\times
K_{0}(|q-q^{\prime}||{\bf R}_{e}-{\bf R}|)\;{\rm e}^{{\rm i}(q^{\prime}-q)z},$

where $\hat{\rho}({\bf r})=\sum_{i}q_{i}\,\delta({\bf r}-{\bf r}_{i})$ is the
charge density operator, expressed as a sum over electrons and ions of charges
$q_{i}$ and positions ${\bf r}_{i}$ in the specimen, and ${\bf R}_{e}$ is the
e-beam impact parameter, set to a fixed value in accordance with assumption
(i) in Sec. A. As anticipated in Sec. C, $G_{qj,q^{\prime}j^{\prime}}$ only
depends on the wave vector difference $q-q^{\prime}$. Incidentally, we neglect
coupling terms with $j=j^{\prime}$, although one could conceivably imagine a
system in which the time-dependent energy shifts produced by the presence of
the electron in the states $j$ could affect the excitation dynamics, as
recently investigated for the nonlinear interaction between low-energy free
electrons and nanographenes Cox and García de Abajo (2020).

For a small sample compared with the impact parameter $R_{e}$, the coupling
coefficients in Eq. (29) can be approximated by only retaining the first
nonvanishing contribution to the Taylor expansion of
$K_{0}(|q-q^{\prime}||{\bf R}_{e}-{\bf R}|){\rm e}^{{\rm i}(q^{\prime}-q)z}$
around ${\bf r}=0$. The zeroth-order term cancels for a neutral sample (i.e,
$\int d^{3}{\bf r}\;\hat{\rho}({\bf r})=\sum_{i}q_{i}=0)$, while the linear
term yields

$\displaystyle G_{qj,q^{\prime}j^{\prime}}=-\frac{e}{\pi\hbar}\,{\bf
p}_{jj^{\prime}}\cdot$
$\displaystyle\bigg{[}|q-q^{\prime}|K_{1}(|q-q^{\prime}|R_{e})\,\hat{\bf
R}_{e}$ (30) $\displaystyle+{\rm
i}(q^{\prime}-q)\,K_{0}(|q-q^{\prime}|R_{e})\,\hat{\bf z}\bigg{]},$

where ${\bf p}_{jj^{\prime}}=\int d^{3}{\bf r}\;\langle j|\hat{\rho}({\bf
r})|j^{\prime}\rangle\;{\bf r}$ is the transition dipole. Inserting this
result in Eq. (26) to obtain the space-dependent interaction in the nonrecoil
approximation, we find the expected dipole potential
$G_{jj^{\prime}}(z)=-(e/\hbar)\,{\bf p}_{jj^{\prime}}\cdot({\bf
R}_{e}+z\,\hat{\bf z})/(R_{e}^{2}+z^{2})^{3/2}$.

A more general treatment that is suitable for dealing with multipolar
excitations can be followed by expanding the Coulomb potential in Eq. (29)
using spherical harmonics as Jackson (1999)

$\frac{1}{|{\bf r}_{e}-{\bf
r}|}=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\frac{4\pi}{2l+1}\,(r^{l}/r_{e}^{l+1})\,Y_{lm}(\hat{\bf
r}_{e})Y_{lm}^{*}(\hat{\bf r})$

under the assumption that $r_{e}>r$ (i.e., provided the sample can be
inscribed in a sphere that is not intersected by the e-beam). Inserting this
expansion into Eq. (29) and using the analytical expression for the $z_{e}$
integral of ${\rm e}^{{\rm i}(q^{\prime}-q)z_{e}}Y_{lm}(\hat{\bf
r}_{e})/r_{e}^{l+1}$ derived in Ref. García de Abajo (1999), we obtain

$\displaystyle
G_{qj,q^{\prime}j^{\prime}}=-\frac{e}{\pi\hbar}\sum_{l=1}^{\infty}\sum_{m=-l}^{l}Q_{lm,jj^{\prime}}\frac{(-{\rm
i})^{l+m}\;{\rm e}^{{\rm i}m\varphi_{{\bf
R}_{e}}}}{\sqrt{(l-m)!(l+m)!}}(q-q^{\prime})^{l}K_{m}(|q-q^{\prime}|R_{e})\times\left\\{\begin{matrix}1,&\quad
q-q^{\prime}>0,\\\ (-1)^{m},&\quad q-q^{\prime}<0,\end{matrix}\right.$ (31)

where

$\displaystyle Q_{lm,jj^{\prime}}=\sqrt{\frac{4\pi}{2l+1}}\int d^{3}{\bf
r}\;\langle j|\hat{\rho}({\bf r})|j^{\prime}\rangle\,r^{l}Y_{lm}^{*}(\hat{\bf
r}),$

$\varphi_{{\bf R}_{e}}$ is the azimuthal angle of ${\bf R}_{e}$, and we
eliminate the $l=0$ term because it vanishes due to charge neutrality in the
sample. Incidentally, we follow the notation of Ref. Messiah (1966) for the
spherical harmonics $Y_{lm}(\hat{\bf r})$, which differs by a factor
$(-1)^{m}$ from that in Ref. Abramowitz and Stegun (1972). Detailed inspection
shows that the $l=1$ term in this expression reproduces the result in Eq.
(30), with $Q_{1m,jj^{\prime}}$ reducing to the components of the transition
dipole.

$\begin{matrix}\text{multipolar}\\\ \text{order}\end{matrix}$ |  $\begin{matrix}\text{charge density}\\\ \text{symmetry}\end{matrix}$ |  $\begin{matrix}\text{charge density}\\\ \langle j|\hat{\rho}({\bf r})|j^{\prime}\rangle\propto\end{matrix}$ |  $g_{jj^{\prime}}\times G_{jj^{\prime}}(z)$ |  $g_{jj^{\prime}}\times G_{qj,0j^{\prime}}$
---|---|---|---|---
dipole |  $p_{x}$ |  $x/r$ |  $\pi\,R_{e}/(R_{e}^{2}+z^{2})^{3/2}$ |  $|q|\;K_{1}\big{(}|q|R_{e}\big{)}$
dipole |  $p_{z}$ |  $z/r$ |  ${\rm i}\pi\,z/(R_{e}^{2}+z^{2})^{3/2}$ |  $q\;K_{0}\big{(}|q|R_{e}\big{)}$
quadrupole |  $d_{z^{2}}$ |  $3z^{2}/r^{2}$ |  $\pi(R_{e}^{2}-2z^{2})/(R_{e}^{2}+z^{2})^{5/2}$ |  $q^{2}\;K_{0}\big{(}|q|R_{e}\big{)}$
quadrupole |  $d_{xz}$ |  $xz/r^{2}$ |  $3{\rm i}\pi\,z\,R_{e}/(R_{e}^{2}+z^{2})^{5/2}$ |  $q^{2}\;K_{1}\big{(}|q|R_{e}\big{)}\;{\rm sign}\\{q\\}$
quadrupole |  $d_{x^{2}-y^{2}}$ |  $(x^{2}-y^{2})/r^{2}$ |  $3\pi\,R_{e}^{2}/(R_{e}^{2}+z^{2})^{5/2}$ |  $q^{2}\;K_{2}\big{(}|q|R_{e}\big{)}$
Table 1: Coupling coefficients for dipolar and quadrupolar
$j^{\prime}\rightarrow j$ transitions with different symmetries. The e-beam
moves along the $z$ axis and intersects the $x$ axis at a distance $R_{e}$
from the sample. A $q$\- and $z$-independent normalization constant
$g_{jj^{\prime}}$ is introduced in the two rightmost columns, encapsulating
the dependence on the radial matrix elements. We note that
$G_{jj^{\prime}}(z)$ is obtained from $G_{qj,0j^{\prime}}(z)$ by using Eq.
(26). The two remaining quadrupolar transitions $d_{xy}$ and $d_{yz}$ do not
couple to the electron with the specified trajectory.

Beyond dipoles, we also explore quadrupolar excitations characterized by
charge densities $\langle j|\hat{\rho}({\bf r})|j^{\prime}\rangle$ with a
spatial angular dependence given by $Y_{20}\propto(3z^{2}/r^{2}-1)$,
$(Y_{2,-1}-Y_{21})\propto xz/r^{2}$, and
$(Y_{2,-2}+Y_{22})\propto(x^{2}-y^{2})/r^{2}$, as found, for example, in
transitions from a hydrogenic $s$ orbital to $d_{z^{2}}$, $d_{xz}$, and
$d_{x^{2}-y^{2}}$ states. We take $\varphi_{{\bf R}_{e}}=0$ (i.e., the e-beam
crosses the $x$ axis), so that the other two possible quadrupolar excitations
($d_{xy}$ and $d_{yz}$) do not couple to the electron because of symmetry
mismatch. In Table 1, we summarize the associated coupling coefficients in
both momentum and real-space representations [$G_{qj,0j^{\prime}}$ and
$G_{jj^{\prime}}(z)$, respectively] for the nonzero dipolar and quadrupolar
transitions under the noted conditions. In this work, we use
$G_{qj,0j^{\prime}}$ in combination with Eq. (17) to produce numerical
nonlinear results including recoil effects, while numerical integration of Eq.
(27) with $G_{jj^{\prime}}(z)$ as input allows us to compute excitation
probabilities in the nonrecoil approximation. We remark that the normalization
constants $g_{jj^{\prime}}$ in Table 1 encapsulate all the factors that
accompany the $(q-q^{\prime})^{l}K_{m}(|q-q^{\prime}|R_{e})$ dependence for
each angular symmetry in Eq. (31), and in particular, they are proportional to
the multipolar transition strength $Q_{lm,jj^{\prime}}$.

In the main text, rather than specifying $g_{jj^{\prime}}$, we express our
results as a function of the first-order (linear) excitation probability
$P_{1}^{\rm lin}$ given by Eq. (16), so we take

$\displaystyle|g_{jj^{\prime}}|=2\pi\sqrt{\frac{|g_{jj^{\prime}}G_{q_{j}j,q_{0}0}|^{2}+|g_{jj^{\prime}}G_{-q_{j}j,q_{0}0}|^{2}}{v_{q_{0}}v_{q_{j}}P_{j}^{\rm
lin}}}$ (32)

with $g_{jj^{\prime}}G_{\pm q_{j}j,q_{0}0}\equiv g_{jj^{\prime}}G_{\pm
q_{j}-q_{0},j,00}$ explicitly given by the rightmost column of Table 1 with
$q=\pm q_{j}-q_{0}$. This prescription needs to be modified in the nonrecoil
approximation, as we argue in Sec. C, so we calculate $|g_{jj^{\prime}}|$
after eliminating the backscattering term
$|g_{jj^{\prime}}G_{-q_{j}j,q_{0}0}|^{2}$ in Eq. (32) to write
$|g_{jj^{\prime}}|=(2\pi/v)|g_{jj^{\prime}}G_{q_{j}j,q_{0}0}|\big{/}\big{(}P_{j}^{\rm
lin}\big{)}^{1/2}$.

## References

  * Nellist and Pennycook (1998) P. D. Nellist and S. J. Pennycook, Phys. Rev. Lett. 81, 4156 (1998).
  * Batson et al. (2002) P. E. Batson, N. Dellby, and O. L. Krivanek, Nature 418, 617 (2002).
  * Krivanek et al. (2014) O. L. Krivanek, T. C. Lovejoy, N. Dellby, T. Aoki, R. W. Carpenter, P. Rez, E. Soignard, J. Zhu, P. E. Batson, M. J. Lagos, et al., Nature 514, 209 (2014).
  * Lagos et al. (2017) M. J. Lagos, A. Trügler, U. Hohenester, and P. E. Batson, Nature 543, 529 (2017).
  * Krivanek et al. (2019) O. L. Krivanek, N. Dellby, J. A. Hachtel, J.-C. Idrobo, M. T. Hotz, B. Plotkin-Swing, N. J. Bacon, A. L. Bleloch, G. J. Corbin, M. V. Hoffman, et al., Ultramicroscopy 203, 60 (2019).
  * Hage et al. (2020) F. S. Hage, G. Radtke, D. M. Kepaptsoglou, M. Lazzeri, and Q. M. Ramasse, Science 367, 1124 (2020).
  * Hachtel et al. (2019) J. A. Hachtel, J. Huang, I. Popovs, S. Jansone-Popova, J. K. Keum, J. Jakowski, T. C. Lovejoy, N. Dellby, O. L. Krivanek, and J. C. Idrobo, Science 363, 525 (2019).
  * Mkhitaryan et al. (2021) V. Mkhitaryan, K. March, E. Tseng, X. Li, L. Scarabelli, L. M. Liz-Marzán, S.-Y. Chen, L. H. G. Tizei, O. Stéphan, J.-M. Song, et al., Nano Lett. 21, 2444 (2021).
  * Yan et al. (2021) X. Yan, C. Liu, C. A. Gadre, L. Gu, T. Aoki, T. C. Lovejoy, N. Dellby, O. L. Krivanek, D. G. Schlom, R. Wu, et al., Nature 589, 65 (2021).
  * Egerton (1996) R. F. Egerton, _Electron Energy-loss Spectroscopy in the Electron Microscope_ (Plenum Press, New York, 1996).
  * Egerton (2003) R. F. Egerton, Micron 34, 127 (2003).
  * Erni and Browning (2005) R. Erni and N. D. Browning, Ultramicroscopy 104, 176 (2005).
  * Brydson (2001) R. Brydson, _Electron Energy Loss Spectroscopy_ (BIOS Scientific Publishers, Oxford, 2001).
  * García de Abajo (2010) F. J. García de Abajo, Rev. Mod. Phys. 82, 209 (2010).
  * Hage et al. (2018) F. S. Hage, R. J. Nicholls, J. R. Yates, D. G. McCulloch, T. C. Lovejoy, N. Dellby, O. L. Krivanek, K. Refson, and Q. M. Ramasse, Sci. Adv. 4, eaar7495 (2018).
  * García de Abajo and Di Giulio (2020) F. J. García de Abajo and V. Di Giulio, ACS Photonics 8, 945 (2020).
  * 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. B. Hayun, O. Reinhardt, S. Tsesses, and I. Kaminer, Nature 582, 50 (2020).
  * Auad et al. (2022) Y. Auad, C. Hamon, M. Tencé, H. Lourenço-Martins, V. Mkhitaryan, O. Stéphan, F. J. García de Abajo, L. H. G. Tizei, and M. Kociak, Nano Lett. 22, 4149 (2022).
  * Rez et al. (2016) P. Rez, T. Aoki, K. March, D. Gur, O. L. Krivanek, N. Dellby, T. C. Lovejoy, S. G. Wolf, and H. Cohen, Nat. Commun. 7, 10945 (2016).
  * Haiber and Crozier (2018) D. M. Haiber and P. A. Crozier, ACS Nano 12, 5463 (2018).
  * Li et al. (2020) N. Li, X. Guo, X. Yang, R. Qi, T. Qiao, Y. Li, R. Shi, Y. Li, K. Liu, Z. Xu, et al., Nat. Mater. 20, 43 (2020).
  * Konečná et al. (2021) A. Konečná, J. Li, J. H. Edgar, F. J. García de Abajo, and J. A. Hachtel, Small 17, 2103404 (2021).
  * Bosman et al. (2007) M. Bosman, V. J. Keast, M. Watanabe, A. I. Maaroof, and M. B. Cortie, Nanotechnology 18, 165505 (2007).
  * Nelayah et al. (2007) J. Nelayah, M. Kociak, O. Stéphan, F. J. García de Abajo, M. Tencé, L. Henrard, D. Taverna, I. Pastoriza-Santos, L. M. Liz-Marzán, and C. Colliex, Nat. Phys. 3, 348 (2007).
  * Rossouw and Botton (2013) D. Rossouw and G. A. Botton, Phys. Rev. Lett. 110, 066801 (2013).
  * Tan et al. (2014) S. F. Tan, L. Wu, J. K. W. Yang, P. Bai, M. Bosman, and C. A. Nijhuis, Science 343, 1496 (2014).
  * Rocca (1995) M. Rocca, Surf. Sci. Rep. 22, 1 (1995).
  * Tromp (2019) R. Tromp, in _Handbook of Microscopy_ , edited by P. W. Hawkes and J. C. H. Spence (Springer, Berlin, 2019), pp. 565–604.
  * García de Abajo (2013) F. J. García de Abajo, ACS Nano 7, 11409 (2013).
  * Cox and García de Abajo (2020) J. D. Cox and F. J. García de Abajo, Nano Lett. 20, 4792 (2020).
  * Barwick et al. (2009) B. Barwick, D. J. Flannigan, and A. H. Zewail, Nature 462, 902 (2009).
  * García de Abajo et al. (2010) F. J. García de Abajo, A. Asenjo-Garcia, and M. Kociak, Nano Lett. 10, 1859 (2010).
  * Feist et al. (2015) A. Feist, K. E. Echternkamp, J. Schauss, S. V. Yalunin, S. Schäfer, and C. Ropers, Nature 521, 200 (2015).
  * 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).
  * 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, et al., Nature 600, 23/30 December (2021).
  * Konečná et al. (2019) A. Konečná, V. Di Giulio, V. Mkhitaryan, C. Ropers, and F. J. García de Abajo, ACS Photonics 7, 1290 (2019).
  * Noziéres (1997) P. Noziéres, _Theory of Interacting Fermi Systems_ (Taylor & Francis, Boca Raton, 1997).
  * Messiah (1966) A. Messiah, _Quantum Mechanics_ (North-Holland, New York, 1966).
  * Di Giulio et al. (2019) V. Di Giulio, M. Kociak, and F. J. García de Abajo, Optica 6, 1524 (2019).
  * Tizei et al. (2015) L. H. G. Tizei, Y.-C. Lin, M. Mukai, H. Sawada, A.-Y. Lu, L.-J. Li, K. Kimoto, and K. Suenaga, Phys. Rev. Lett. 114, 107601 (2015).
  * Schuler et al. (2019) B. Schuler, K. A. Cochrane, C. Kastl, E. S. Barnard, E. Wong, N. J. Borys, A. M. Schwartzberg, D. F. Ogletree, F. J. García de Abajo, and A. Weber-Bargioni, Sci. Adv. 6, eabb5988 (2019).
  * Gericke et al. (2008) T. Gericke, P. Würtz, D. Reitz, T. Langen, and H. Ott, Nat. Phys. 4, 949 (2008).
  * not (a) This result comes from the $\delta$-function contribution of $(\varepsilon_{q^{\prime}q_{0}}+\omega_{j^{\prime}0}-{\rm i}0^{+})^{-1}={\rm PV}\left\\{(\varepsilon_{q^{\prime}q_{0}}+\omega_{j^{\prime}0})^{-1}\right\\}+{\rm i}\pi\delta(\varepsilon_{q^{\prime}q_{0}}+\omega_{j^{\prime}0})$ to Eq. (15), where PV stands for the principal value, which vanishes identically.
  * not (b) Explicit expressions for $G_{qj,q^{\prime}j^{\prime}}$ in the nonretarded limit (see Sec. D) show that these coupling coefficients only depend on the difference of wave vectors $q-q^{\prime}$. This result is also maintained when including retardation in the nonrecoil approximation, since then we can write $G_{qj,q^{\prime}j^{\prime}}=\frac{1}{\hbar}\langle qj|\hat{\mathcal{H}}|q^{\prime}j^{\prime}\rangle=(e{\bf v}/\hbar c)\cdot\langle qj|\hat{\mathcal{A}}|q^{\prime}j^{\prime}\rangle=(e{\bf v}/2\pi\hbar c)\cdot\int dz_{e}\,{\rm e}^{{\rm i}(q^{\prime}-q)z_{e}}\langle j|\hat{\mathcal{A}}({\bf r}_{e})|j^{\prime}\rangle$, where $\hat{\mathcal{A}}({\bf r}_{e})$ is the vector potential operator mediating the electron-sample interaction [see F. J. García de Abajo and V. Di Giulio, ACS Photonics 8, 945 (2020)].
  * Jackson (1999) J. D. Jackson, _Classical Electrodynamics_ (Wiley, New York, 1999).
  * García de Abajo (1999) F. J. García de Abajo, Phys. Rev. B 59, 3095 (1999).
  * Abramowitz and Stegun (1972) M. Abramowitz and I. A. Stegun, _Handbook of Mathematical Functions_ (Dover, New York, 1972).