# Atomiclike Selection Rules in Free Electron Scattering

Simon Garrigou, Hugo Lourenço-Martins

► To cite this version:

Simon Garrigou, Hugo Lourenço-Martins. Atomiclike Selection Rules in Free Electron Scattering. Physical Review Letters, 2025, 134 (25), pp.256902. 10.1103/12y2-99jx. hal-05374151

HAL Id: hal-05374151

https://hal.science/hal-05374151v1

Submitted on 20 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Atomic-like selection rules in free electron scattering

Simon Garrigou¹ and Hugo Lourenço-Martins¹, *

¹CEMES-CNRS, Université de Toulouse, CNRS, Toulouse, France

(Dated: November 19, 2024)

Phase-shaped electron energy-loss spectroscopy (PSEELS) measures the scattering probability of structured free electron beams by a target. Over the last decade, it was shown that this scheme can be employed to emulate polarized optical spectroscopies with electrons, and therefore to transpose macroscopic optical concepts - such as dichroism - down to the deep sub-wavelength scale. In this work, we theoretically demonstrate that PSEELS can, in fact, go way further than mimicking optics and enables to map a plethora of so far inaccessible nano-optical quantities such as the electric quadrupolar momentum.

Spontaneous emission refers to the process in which a quantum emitter (QE, i.e. a few-levels bound electron system) transits between two states $|e\rangle \to |g\rangle$ by emitting a photon. This free-space decay rate $\Gamma_0$ can be enhanced by weakly coupling the emitter to any type of dielectric environment [1]. Indeed, this medium increases the local density $\rho$ of electromagnetic states (EMLDOS) in which the emitter can decay into, an effect intuitively captured in the celebrated Purcell formula [2, 3], see Fig. 1(a):

$$ \Gamma = \frac{2\omega}{3\hbar\varepsilon_0} |d|^2 \rho_d(\mathbf{r}_0, \omega) \quad (1) $$

with $\mathbf{r}_0 \in \mathbb{R}^3$ the position of the emitter, $\hbar$ the reduced Planck constant, $\omega$ the angular frequency of emitted radiation, $\varepsilon_0$ the vacuum dielectric permittivity, $d = \langle e|\dot{d}|g \rangle$ the dipole moment of the transition and $\rho_d$ the EMLDOS along the dipole moment.

This whole phenomenology is not restricted to bound electron systems such as quantum dots [4], atoms [5] or molecules [6], but can also be observed with free electron states in a scheme called electron energy-loss spectroscopy (EELS) in a transmission electron microscope (TEM). EELS consists in analysing the weak energy-losses undergone by swift electrons ($v \sim 80\%$ of the speed of light) during their inelastic interaction with a nano-object, thus revealing its sub-ångström dynamics, as recently demonstrated on single atoms [7]. Remarkably, a connection between the z-Fourier transform $\tilde{\rho}$ of the EMLDOS and the observable measured in EELS - the probability $\Gamma^{\text{EELS}}$ for the electron to lose an energy $\hbar\omega$ - was established [8]:

$$ \Gamma^{\text{EELS}} = \frac{2\pi e^2}{\hbar\omega} \tilde{\rho}_z(\mathbf{R}_0, q_z, \omega) \quad (2) $$

where $q_z = k_f - k_i = \omega/v$ is the transferred momentum, $k_i, k_f$ the initial and final electron momenta, $z$ conventionally represents the propagation axis of the electron beam, and $\mathbf{R}_0$ the impact parameter of the beam in its transverse plane. The similarity between equations (1) and (2) suggested a deep analogy between bound and free electron spontaneous scattering, albeit without explicitly involving any transition dipole moment - a key ingredient of the Purcell effect (1).

This difficulty was lifted over the last decade by a series of theoretical works [9–11] which proposed to shape [12, 13] and post-select [14, 15] the electron beam wavefront. Indeed, swift electron beams being highly paraxial, this approach simply decorates the initial EELS planewave scattering problem $|k_i\rangle \to |k_f\rangle$ by the symmetry of the transverse wavefunctions $|k_i\rangle \otimes |\Psi_{i,\perp}\rangle \to |k_f\rangle \otimes |\Psi_{f,\perp}\rangle$, see Fig. 1(b). This additional degree of freedom introduces new selection rules to the transition probabilities (2). In particular, it was shown that the probability of transitioning from an Hermite-Laguerre-Gauss (HLG) state $|\Psi\rangle$ of order $i = 1$ - representing topologically singular beams carrying one quantum of linear or orbital angular momenta [16, 17] - to an unshaped Gaussian wavefront $|G\rangle$ can be written [18–21]:

$$ \Gamma^{\text{pEELS}} = \frac{2\pi q_z^2}{\hbar\omega} |d_{\Psi}|^2 \rho_d(\mathbf{R}_0, q_z, \omega) \quad (3) $$

where $d_{\Psi} = \langle G|\dot{d}|\Psi \rangle$. This scheme - referred to as polarized EELS (pEELS) - unlocks the measurement of the EMLDOS along any polarization axis in EELS [22–24] and completes the analogy between free electron scattering (3) and the Purcell phenomenology (1).

In this letter, we demonstrate that EELS can be pushed beyond the dipole order (3) by considering the scattering of HLG states of arbitrary initial $i$ and final $j$ orders with $|i - j| > 1$, thus surpassing the Purcell analogy. Our formalism puts the emphasis on the quantum numbers (i.e. linear and angular momentum) conservation, therefore making selection rules explicit at all orders and providing an intuitive picture of the scattering mechanisms. In the spirit of atomic spectroscopy, this simplicity enables us to draw look-up tables providing any experimentalist with the nano-optical quantity probed by EELS for any integer $i$ and $j$. Eventually, we close the letter with a series of three examples, first checking that both EELS and pEELS are recovered for $|i - j| = 0$ and $|i - j| = 1$ respectively, then demonstrating that the transition $|i - j| = 2$ measures the quadrupolar component of nano-optical fields - thus showing that transitions satisfying $|i - j| > 1$ access quantities unreachable with most of all-optical techniques.
<figure><img src="image_2.png" /><figcaption>Fig. 1 Bound and free electrons spontaneous scattering. a. Purcell effect: a non-trivial dielectric environment of EMLDOS ρ - represented here by a cavity - enhances the spontaneous decay rate of a quantum emitter (QE). b. PSEELS: the exact same configuration as in a. but replacing the bound electron states by free electron orbitals i.e. planewaves |k⟩ decorated by structured wavefronts |Ψ⊥⟩. Ṽ<sub>proj</sub> the projected potential defined in (4)</figcaption></figure>

We begin this letter by presenting our model: the exact details of the derivation being given in the supplementary material, we will hereafter only focus on the main steps and approximations. We start by considering a fast electron beam of longitudinal (z)-momentum $k_i$ and arbitrary wavefront $\Psi_{i,\perp}$. Within the paraxial approximation, the electron wavefunction can be separated between longitudinal and transverse states $|\psi\rangle = |k\rangle \otimes |\Psi_{\perp}\rangle$. The first step of the derivation is to compute the interaction probability between the electron beam and a target represented by a set of quantum states $\{|n\rangle\}_{n \in \mathbb{N}}$. Here, we make no assumption on the nature of the target which can be any type of sample (e.g. phononic, photonic or plasmonic material, molecule, van der Waals material, etc), we will always describe it as a polariton field [25–27], and thus generically denotes them as *photonic states*. We moreover impose a set of common approximations valid for swift electron scattering [28]: Quasi-static ($c \to +\infty$), single scattering (weak interaction), non-recoil (large electron velocity $v$) and local potential (narrow beam waist) approximations. Within this framework, one can show that the electron-target interaction potential reduces to a projected potential:

$$ \hat{V}_{\text{proj}}(q_z, \omega) = \int dz \langle z | \hat{V} | z \rangle e^{iq_z z} \quad (4) $$

where the longitudinal degree of freedom of the electron has been absorbed in the definition of an effective projected potential $\hat{V}_{\text{proj}}$ [29–32], translating the delocalization of interaction along the electron path, see Fig. 1(b). The interaction probability is then given by a direct ap-

plication of the Fermi’s golden rule:

$$ \begin{aligned} \Gamma^{\text{PSEELS}}(\omega) &= \sum_{n}^{\tilde{n}} \left| \frac{1}{\hbar v} \langle \Psi_{f,\perp} | \langle n | \hat{V}_{\text{proj}} | 0 \rangle | \Psi_{i,\perp} \rangle \right|^2 \\ &= \sum_{n}^{\tilde{n}} |g_{0,n}|^2 \end{aligned} \quad (5) $$

where the energy conservation is encompassed in the new definition of the sum $\sum_{n}^{\tilde{n}} = \sum_{n} \delta(\omega_0 - \omega_n - \omega)$. One can see that, when preparing $|\Psi_{i,\perp}\rangle$ and post-selecting $|\Psi_{f,\perp}\rangle$ transverse electron states (see Fig. 1(d)) - a scheme denoted as phase-shaped EELS (PSEELS) - one selects a specific transition occurring between discrete electronic states and for all accessible photonic states $|n\rangle$, represented by the coupling constant $g_{0,n}$. The total loss-probability (5) is then given by the sum over all the possible final photonic states. For pedagogy, in the rest of the paper, we will drop this sum and focus on one specific transition in the target $g_{0,n}$. From (5), the discrete electron transverse states essentially behave as atomic orbitals, imposing specific selection rules by their symmetry. However, while atomic orbitals are three dimensional objects with shapes fixed by the nature of the atom, free electron orbitals (i.e $|\Psi_{i,\perp}\rangle$ and $|\Psi_{f,\perp}\rangle$) are two dimensional and their shape can be chosen arbitrarily through beam shaping and post-selection.

In order to explicit the connection between free electron orbitals and atomic orbitals, as well as the corresponding selection rules, we will now perform a multipolar development of the projected potential. Indeed, since the electron beam waist (Angstrom) is small compared to the electrostatic potential variations scale (nanometer), one can consider the projected potential and its derivative as constant over the electron wavefunction, and represented by the function $\hat{V}_{\text{proj}}(\mathbf{R}_0, q_z, \omega)$, where $\mathbf{R}_0$ is the impact parameter of the electron in the transverse plane. In this limit, one can perform a multipolar development equivalent to the Power-Zienau-Woolley form of the minimal-coupling [1]:

$$ \hat{V}_{\text{proj}}(q_z, \omega) = \underbrace{\hat{\rho}\hat{\Phi}(q_z, \omega)}_{\hat{V}_0} - \underbrace{\hat{d} \cdot \hat{E}(q_z, \omega)}_{\hat{V}_1} - \underbrace{\overleftrightarrow{Q} \nabla \cdot \hat{E}(q_z, \omega)}_{\hat{V}_2} + \dots \quad (6) $$

where the term involving the $k^{\text{th}}$ derivative of the potential will be noted $\hat{V}_k$. As we will confirm further, the term $\hat{d}$ in $\hat{V}_1$ is the dipolar operator at the origin of (3). Now, each term of $\hat{V}_{\text{proj}}$ is composed of two operators, acting either on the electron transverse degree of freedom or the photonic degrees of freedom. Plugged back into (5), the latter gives expectation values of the type $\langle n | \hat{A} | 0 \rangle = A_{n0}$, with $\hat{A}$ a field operator, e.g $\hat{E}$ or $\hat{\Phi}$. In later developments, we will focus on a single $n$ photonic mode of the system and omit the $n0$ subscript for clarity.

The last step of our derivation consists in giving an alternative and more intuitive form to the electron $i^{\text{th}}$
momentum operators $\hat{\rho}, \hat{\mathbf{d}}, \hat{\vec{Q}}$, .... This is done by describing the transverse electron field with a ladder operator formalism typically encountered in the two dimensional quantum harmonic oscillator. Introducing the annihilation $\hat{a}_u$ and creation $\hat{a}_u^\dagger$ operators along direction $u$. With $u \in \{x, y\}$ or $u \in \{\circ, \mathring{\circ}\}$, linear and circular basis respectively. The electronic position operators can be re-written as [33]:

$$ \hat{x} = \frac{w_0}{2} (\hat{a}_x + \hat{a}_x^\dagger) = \frac{w_0}{2\sqrt{2}} (\hat{a}_\circ + \hat{a}_\circ + \hat{a}_\circ^\dagger + \hat{a}_\circ^\dagger) \quad (7) $$

$$ \hat{y} = \frac{w_0}{2} (\hat{a}_y + \hat{a}_y^\dagger) = \frac{i w_0}{2\sqrt{2}} (\hat{a}_\circ - \hat{a}_\circ - \hat{a}_\circ^\dagger + \hat{a}_\circ^\dagger) \quad (8) $$

With the circular operators $\hat{a}_\circ = \frac{1}{\sqrt{2}}(\hat{a}_x - i\hat{a}_y)$ and $\hat{a}_\circ = \frac{1}{\sqrt{2}}(\hat{a}_x + i\hat{a}_y)$. Hence, the dipolar term of the interaction potential in linear basis now simply reads:

$$ -\hat{\mathbf{d}}.\mathbf{E} = \frac{ew_0}{2} [(\hat{a}_x + \hat{a}_x^\dagger)E_x + (\hat{a}_y + \hat{a}_y^\dagger)E_y] \quad (9) $$

Note that here, since the photonic degrees of freedom have been absorbed in $\mathbf{E}$, the previous expression corresponds to a pure transverse electron operator. Every term (quadrupolar, octupolar...) of (6) is now expressed as a sum/product of these ladder operators. Remarkably, one can show that each term $\hat{V}_k$ of the interaction potential (6) only involves products $\hat{a}^{(k)}$ of k ladder operators, annihilation and/or creation along directions $u \in \{x, y, \circ, \mathring{\circ}\}$:

$$ \langle n | \hat{V}_{\text{proj}} | 0 \rangle = \sum_k \langle n | \hat{V}_k [\hat{a}^{(k)}] | 0 \rangle \quad (10) $$

The integer $k$ will be called the order of the transition in the following. For instance, the term (9) represents a transition of order 1 where ladder operators only appear to the power one.

Now, why is this new form of the electron operators in terms of ladder operator advantageous? As mentioned in the introduction, the vast majority of electron beam wavefronts considered in both the experimental and theoretical literature are Hermite- (HG) and Laguerre-Gauss (LG) states. The key property of these sets of states is that they also form two eigenbasis of the 2D quantum harmonic oscillator. Thus, as demonstrated in the pioneering work of G. Nienhuis and L. Allen [34], all the properties of the HG and LG light beams can be expressed in term of ladder operators, in complete analogy with the quantum harmonic oscillator. Here, we reuse this technique in our model to describe not light but electron HLG beams, the method being essentially the same. In particular, every electron transverse state can be built by the successive application of ladder operators on a simple Gaussian state $|G\rangle$:

$$ |HG_{n_x, n_y}\rangle = \frac{1}{\sqrt{n_x! n_y!}} (\hat{a}_x^\dagger)^{n_x} (\hat{a}_y^\dagger)^{n_y} |G\rangle \quad (11) $$

$$ |LG_{n_\circ, n_\circ}\rangle = \frac{1}{\sqrt{n_\circ! n_\circ!}} (\hat{a}_\circ^\dagger)^{n_\circ} (\hat{a}_\circ^\dagger)^{n_\circ} |G\rangle \quad (12) $$

In this text, we use the notation $|LG_{n_\circ, n_\circ}\rangle$ which deviates from common $|LG_p^l\rangle$, the correspondence being given by $p = \min(n_\circ, n_\circ)$ and $l = n_\circ - n_\circ$. In the exact same way as we did for the interaction potential, we define the order $i$ of a pure HG or LG beam profile as the number of creation operators required to build it from the Gaussian state. Any transverse state can be transformed into another transverse state by application of ladder operators. Then, the order $|i-j|$ of the transition between two transverse states $i$ and $j$ is defined as the number of operators needed to transit from $i$ to $j$. A first order transition corresponds to $|HG_{1,0}\rangle \to |G\rangle$ or $|LG_{1,0}\rangle \to |LG_{1,1}\rangle$ for instance.

At this stage, the strength of this approach starts to be visible: both the multipolar interaction potential and the transverse electron states are now expressed through the same ladder operators: computing transition probability now only relies on simple Hilbert space algebra. But the deep interest of this model is revealed when calculating any transition between two e.g. HG states of order $i$ and $j$, one gets to leading order:

$$ g_{0,i \to j} = \sum_k \mathcal{M}^{(k)} [\hat{V}_k] \delta_{k, |i-j|} \quad (13) $$

where $\mathcal{M}^{(k)}$ corresponds to the transition of order $i$ involving only the multipolar term $\hat{V}_k$. The analytical expressions of the amplitudes $\mathcal{M}^{(k)}$ can be quite complex - especially for increasing order $k$ - and are exactly derived in the supplementary material. However, one can already notice that, remarkably, the interaction probability will be non-zero if and only if the order of the transition in the transverse electron state matches the order of a multipolar term of the projected potential, expressed through the condition $\delta_{k, |i-j|}$, with $\delta_{\mu, \nu}$ the Kronecker delta. The latter is reminiscent of the emergence of selection rules and is responsible for a pairing between the order of the transition $|i-j|$ and a specific multipolar component of the field, as shown on tables I and II. We visualize here that different experiments sharing the same $|i-j|$ probe in fact the same quantity, the deciding factor being the nature of the transition more than that of the initial and final state.

This pairing is not accidental but instead a consequence of the Noether theorem. For instance, consider a transition between two HG states with $\Delta n_x = 1$ and $\Delta n_y = 0$. Then, the corresponding amplitude will scale as $g_{0,n_x n_y \to n'_x n'_y}^{HG} \propto E_x$. Keeping in mind that $n_x$ represents a quantum of linear momentum in the direction $x$, the selection rule expresses the fact that if the electron interacts with a field oriented along the $x$-direction, it must undergo a transverse momentum kick in the same direction.

In order to further explore the mechanics behind these conservation/selection rules, we will exemplify three transitions in the rest of this letter. We start with the simplest
<table><caption>**Table I** |$\Psi_{f,\perp}$〉 → |$\Psi_{i,\perp}$〉 PSEELS transition amplitudes for Hermite-Gauss beam shapes. The amplitudes are given only as the probed nanophotonic quantity, for brevity we remove their dependency ($\mathbf{R}_0, q_z, \omega$). Beam profiles are plotted in amplitude of the wavefunction, red for positive, blue for negative. ($n_x, n_y$) index of the wavefront is displayed in the bottom left corner.</caption><thead><tr><th>$|\Psi_{i,\perp}\rangle$</th><th>$(0,0)$</th><th>$(1,0)$</th><th>$(1,1)$</th></tr></thead><tbody><tr><th>$(0,0)$</th><td>$\Phi$</td><td>$E_x$</td><td>$\partial_x E_y + \partial_y E_x$</td></tr><tr><th>$(1,0)$</th><td>$E_x$</td><td>$\Phi$</td><td>$E_y$</td></tr><tr><th>$(1,1)$</th><td>$\partial_x E_y + \partial_y E_x$</td><td>$E_y$</td><td>$\Phi$</td></tr><tr><th>$(0,1)$</th><td>$E_y$</td><td>$\partial_x E_y + \partial_y E_x$</td><td>$E_x$</td></tr><tr><th>$(2,0)$</th><td>$\partial_x E_x$</td><td>$E_x$</td><td>$\partial_x E_y + \partial_y E_x$</td></tr></tbody></table>

<table><caption>**Table II** |$\Psi_{f,\perp}$〉 → |$\Psi_{i,\perp}$〉 PSEELS transition amplitudes for Laguerre-Gauss beam shapes. The amplitudes are simplified to the probed nanophotonic quantity, for brevity we remove their dependency ($\mathbf{R}_0, q_z, \omega$). Beam profiles are plotted in the squared norm of the wavefunction, the orbital angular momentum of the mode is displayed in the upper left corner and ($n_\odot, n_\phi$) index of the wavefront is displayed in the bottom left corner.</caption><thead><tr><th>$|\Psi_{i,\perp}\rangle$</th><th>0<br>$(0,0)$</th><th>+1<br>$(1,0)$</th><th>0<br>$(1,1)$</th></tr></thead><tbody><tr><th>0<br>$(0,0)$</th><td>$\Phi$</td><td>$E_\odot$</td><td>$\partial_\odot E_\odot + \partial_\odot E_\phi$</td></tr><tr><th>+1<br>$(1,0)$</th><td>$E_\odot$</td><td>$\Phi$</td><td>$E_\odot$</td></tr><tr><th>0<br>$(1,1)$</th><td>$\partial_\odot E_\phi + \partial_\odot E_\phi$</td><td>$E_\odot$</td><td>$\Phi$</td></tr><tr><th>-1<br>$(0,1)$</th><td>$E_\odot$</td><td>$\partial_\odot E_\phi$</td><td>$E_\odot$</td></tr><tr><th>+2<br>$(2,0)$</th><td>$\partial_\odot E_\phi$</td><td>$E_\odot$</td><td>$\partial_\odot E_\phi$</td></tr></tbody></table>

transition: $|G\rangle \to |G\rangle$. The initial and final beams being the same, it corresponds to a $|i-j|=0^{\text{th}}$ order transition. Thus, only the first term of the multipolar development contributes i.e.:

$$g_{0,00\to00}^{HG} = \frac{e}{\hbar v} \Phi \underbrace{\langle G | \mathbf{1} | G \rangle}_{\delta_{0,0} \delta_{0,0}} \quad (14)$$

where we omitted the dependence over $\mathbf{R}_0, q_z$ and $\omega$ for brevity. From this equation, we immediately get: $\Gamma^{\text{PSEELS}} = |g_{0,00\to00}^{HG}|^2 = \frac{2\pi e^2}{\hbar\omega} \tilde{\rho}_z(\mathbf{R}_0, q_z, \omega)$. We retrieve thus the conventional EELS probability (2). This is

expected since in an EELS experiment, the initial beam is roughly gaussian and with no post selection the outgoing beam will be a superposition of transverse shapes. The dominant term in probability will come from the Gaussian part of the outgoing beam. More generally, in any transition without transverse momentum exchange $\Delta n_x = \Delta n_y = 0$, the transverse electromagnetic field must have no influence and the probability must only depend on $E_z$. These EELS-like transitions are highlighted by yellow cells in tables I and II.

We now move to the next order an examine $|i-j|=1^{\text{st}}$ order transition, for instance to a transition $|LG_{10}\rangle \to |G\rangle$, being the interaction of a left vortex beam with the target followed by the post-selection of a Gaussian wavefront. We readily obtain:

$$g_{0,10\to00}^{LG} = \frac{ew_0}{2\hbar v} E_\odot \underbrace{\langle G | \hat{a}_\odot | LG_{10} \rangle}_{\delta_{1-1,0} \delta_{0,0}} \quad (15)$$

with $E_\odot$ the left circularly polarized electric field component. This transition corresponds to an electron state losing a quantum $\Delta n_\odot = -1$ of orbital angular momentum (OAM) and an optical field gaining a quantum of spin angular momentum (polarization) associated to the corresponding field amplitude $E_\odot$. This effect - usually referred to as spin-orbit coupling - again illustrates the Noether theorem and was already observed in electron beams, e.g. [35]. More generally, for arbitrary impinging HLG state, we retrieve the pEELS probability (3) highlighted by blue cells in tables I and II. Thus, one can see that at order 1, only the dipolar term (9) contributes. A quantum field theory based treatment of this problem identified $\langle \hat{d} \rangle$ to the polarization of the effective photon of the interaction - a quantity usually called optical polarization analogue (OPA, [18–21]). This shows that the electron beam then effectively behaves as a nanoscale source of polarized white light, thus enabling electron microscopes to push most of the polarization-selective optical experiment down to the deep sub-wavelength regime.

So far, we have successfully recovered all the important results of EELS and pEELS of the literature, we shall now go further. Indeed, our model enables us to explore higher arbitrary orders without any further computational effort. We illustrate this on the next $|i-j|=2^{\text{nd}}$ order, in particular the transition $|HG_{11}\rangle \to |G\rangle$. We immediately get:

$$g_{0,11\to00}^{HG} = \frac{ew_0^2}{8\hbar v} (\partial_x E_y + \partial_y E_x) \quad (16)$$

The right-hand side of the equation corresponds to the definition the local quadrupolar symmetry of a near-field shown on figure Fig.2(b). Thus enabling to measure this quantity at a precise position on the sample and then to map it by scanning the electron beam over the sample as shown in 2(c,d).
<figure><img src="image_3.png" /><figcaption>Fig. 2 (a) 3D representation of four silver nanorods in-plane. (b) 3D representation of the electric field lines in the central gap, xy-plane top view with field direction arrows on the field lines. Colors of the field line representing the orientation angle $\theta_{xy}$ of the electric field projected on the xy-plane. (c,d) Spatially resolved maps across the gap: (c) PSEELS probability of transfer between a Hermite-Gauss $|HG_{1,1}\rangle$ and a Gaussian transverse state $|G\rangle$ for a 2nm wide electron beam; (d) $|\partial_y E_x + \partial_x E_y|^2$, predicted to be proportional to the quantity measured in (c). Simulations done under MNPBEM library [36].</figcaption></figure>

Thinking back in terms of conservation law: during this quadrupolar transition the electron exchanges two quanta of linear momentum with the field $\Delta n_x = -1$ and $\Delta n_y = -1$. However, the electron still loses one quantum of energy $\hbar\omega$ i.e. one photon. Thus, one can see that for order > 1 transitions - highlighted by green cells in tables I and II - the standard interpretation in terms of OPA employed in pEELS [18–21] fails: it is not possible to identify all the exchanged quanta to the sole polarization of an effective photon. This shows that for $k > 1$, PSEELS goes further than mimicking optics and enters a regime unreachable with standard all-optical methods, becoming able to map so far inaccessible nano-optical quantities.

Our formalism gives a simultaneous and intuitive treatment of all the possible transitions between Hermite-Gauss or Laguerre-Gauss states in PSEELS. It enables the calculation of the analytical transition amplitude between any wavefronts as a sum of all the possible elementary HLG processes (13) by expressing the wavefronts on the HG or LG basis. Conversely, it enables us to infer - without the need of heavy numerical computations - which

beams to use to access any nano-optical quantities in PSEELS. Any optical quantity $C$ which can be expressed as a linear combination of the potential and its derivative $C(\{\phi, E, \partial_u E, \dots, \partial_u^k E, \dots\})$ is given by the exact same linear combination of HLG wavefronts associated with the corresponding transitions.

In all of our developments the target degrees of freedom have been traced over $\sum_n \langle n | \cdot | 0 \rangle$. It should be kept in mind that EELS actually produces a sum of entangled target-electron states [37, 38] e.g. for an initial electron-target state of the form $|HG_{n_x, n_y}, 0\rangle$, one gets:

$$|\psi_f\rangle = \sum_{n_x', n_y', n} g_{0, n_x n_y \to n_x' n_y', n}^{HG} |HG_{n_x', n_y', n}\rangle \quad (17)$$

Thus, the selection of a specific transition $n_x n_y \to n_x' n_y'$ not only gives access to specific information in EELS encoded in $g_0$, but also post-selects the corresponding quantum state of the target. Applying the latest developments in electron-heralded photon state generation [39, 40], the post-selection of specific free electron transitions (e.g. $|i-j|=2$) could enable the generation of quantum states of light with any desired symmetry (e.g. quadrupolar). Eventually, common schemes in quantum optics harness the desexcitation of atoms with cavities to generate specific quantum state of light [41], albeit with parameters set by the intrinsic properties of the atom. Employing free electron orbitals instead of atomic orbitals would enable to bring flexibility to this method: here the nature of the "effective atom" is set by the electron beam shaper and sorter.

This project has been funded in part by the French National Agency for Research under the program of future investment QUENOT (ANR-20-CE30-0033). The authors thank Mathieu Kociak and Jesse Groenen for illuminating discussions on the topic and for the help and support during the redaction of the manuscript.

## References

* hugo.lourenco-martins@cnrs.fr

[1] Lukas Novotny and Bert Hecht. *Principles of Nano-Optics*. Cambridge University Press, Cambridge, reissue edition edition, June 2011.

[2] E M Purcell, H C Torrey, and R V Pound. Resonance Absorption by Nuclear Magnetic Moments in a Solid. *Phys. Rev.*, 69(1-2):37–38, January 1946.

[3] Ulrich Hohenester. *Nano and Quantum Optics: An Introduction to Basic Principles and Theory*. Graduate Texts in Physics. Springer International Publishing, 2020.

[4] A. Dousse, L. Lanco, J. Suffczyński, E. Semenova, A. Miller, A. Lemaitre, I. Sagnes, C. Roblin, J. Bloch, and P. Senellart. Controlled Light-Matter Coupling for a Single Quantum Dot Embedded in a Pillar Microcavity Using
Far-Field Optical Lithography. *Physical Review Letters*, 101(26):267404, December 2008.

[5] J. Gallego, W. Alt, T. Macha, M. Martinez-Dorantes, D. Pandey, and D. Meschede. Strong Purcell Effect on a Neutral Atom Trapped in an Open Fiber Cavity. *Physical Review Letters*, 121(17):173603, October 2018.

[6] Anna Rosławska, Tomáš Neuman, Benjamin Doppagne, Andrei G. Borisov, Michelangelo Romeo, Fabrice Scheurer, Javier Aizpurua, and Guillaume Schull. Mapping Lamb, Stark, and Purcell Effects at a Chromophore-Picocavity Junction with Hyper-Resolved Fluorescence Microscopy. *Physical Review X*, 12(1):011012, January 2022.

[7] F. S. Hage, G. Radtke, D. M. Kepaptsoglou, M. Lazzeri, and Q. M. Ramasse. Single-atom vibrational spectroscopy in the scanning transmission electron microscope. *Science*, 367(6482):1124–1127, March 2020.

[8] F. J. García de Abajo and M. Kociak. Probing the Photonic Local Density of States with Electron Energy Loss Spectroscopy. *Physical Review Letters*, 100(10):106804, 2008.

[9] a. Asenjo-Garcia and F. J. García De Abajo. Dichroism in the interaction between vortex electron beams, plasmons, and molecules. *Physical Review Letters*, 113(6):1–5, 2014.

[10] Daniel Ugarte and Caterina Ducati. Controlling Multipolar Surface Plasmon Excitation through the Azimuthal Phase Structure of Electron Vortex Beams. 205418:1–9, 2015.

[11] Ruben Van Boxem, Bart Partoens, and Jo Verbeeck. Inelastic electron-vortex-beam scattering. *Physical Review A*, 91(3):032703, March 2015.

[12] Jo Verbeeck, Armand Béché, Knut Müller-Caspary, Giulio Guzzinati, Minh Anh Luong, and Martien Den Hertog. Demonstration of a 2 × 2 programmable phase plate for electrons. *Ultramicroscopy*, 190:58–65, July 2018.

[13] Chu-Ping Yu, Francisco Vega Ibañez, Armand Béché, and Johan Verbeeck. Quantum wavefront shaping with a 48-element programmable phase plate for electrons. *SciPost Physics*, 15(6):223, December 2023.

[14] Vincenzo Grillo, Amir H. Tavabi, Federico Venturi, Hugo Larocque, Roberto Balboni, Gian Carlo Gazzadi, Stefano Frabboni, Peng-Han Lu, Erfan Mafakheri, Frédéric Bouchard, Rafal E. Dunin-Borkowski, Robert W. Boyd, Martin P. J. Lavery, Miles J. Padgett, and Ebrahim Karimi. Measuring the orbital angular momentum spectrum of an electron beam. *Nature Communications*, 8(1):15536, May 2017.

[15] Amir H. Tavabi, Paolo Rosi, Enzo Rotunno, Alberto Roncaglia, Luca Belsito, Stefano Frabboni, Giulio Pozzi, Gian Carlo Gazzadi, Peng-Han Lu, Robert Nijland, Moumita Ghosh, Peter Tiemeijer, Ebrahim Karimi, Rafal E. Dunin-Borkowski, and Vincenzo Grillo. Experimental Demonstration of an Electrostatic Orbital Angular Momentum Sorter for Electron Beams. *Physical Review Letters*, 126(9):094802, March 2021.

[16] K. Y. Bliokh, I. P. Ivanov, G. Guzzinati, L. Clark, R. Van Boxem, A. Béché, R. Juchtmans, M. A. Alonso, P. Schattschneider, F. Nori, and J. Verbeeck. Theory and applications of free-electron vortex states. *Physics Reports*, 690:1–70, May 2017.

[17] Yijie Shen, Xuejiao Wang, Zhenwei Xie, Changjun Min, Xing Fu, Qiang Liu, Mali Gong, and Xiaocong Yuan. Optical vortices 30 years on: OAM manipulation from topological charge to multiple singularities. *Light: Science & Applications*, 8(1):90, October 2019.

[18] Hugo Lourenço-Martins, Davy Gérard, and Mathieu Koziak. Optical polarization analogue in free electron beams. *Nature Physics*, 17(5):598–603, May 2021.

[19] Marc R. Bourgeois, Austin G. Nixon, Matthieu Chalifour, Elliot K. Beutler, and David J. Masiello. Polarization-Resolved Electron Energy Gain Nanospectroscopy With Phase-Structured Electron Beams. *Nano Letters*, page acs.nanolett.2c02375, August 2022.

[20] Marc R. Bourgeois, Austin G. Nixon, Matthieu Chalifour, and David J. Masiello. Optical Polarization Analogs in Inelastic Free Electron Scattering, May 2023.

[21] Austin G. Nixon, Matthieu Chalifour, Marc R. Bourgeois, Michael Sanchez, and David J. Masiello. Inelastic scattering of transversely structured free electrons from nanophotonic targets: Theory and computation. *Physical Review A*, 109(4):043502, April 2024.

[22] Michael Schüler and Jamal Berakdar. Energy-loss spectroscopy of C60 fullerenes with twisted electrons: Influence of orbital-angular-momentum transfer on plasmon generation. *Physical Review A*, 94(5):052710, November 2016.

[23] Matteo Zanfrognini, Enzo Rotunno, Stefano Frabboni, Alicia Sit, Ebrahim Karimi, Ulrich Hohenester, and Vincenzo Grillo. Orbital Angular Momentum and Energy Loss Characterization of Plasmonic Excitations in Metallic Nanostructures in TEM. *ACS Photonics*, 6(3):620–627, March 2019.

[24] Andrea Konečná, Mikołaj K. Schmidt, Rainer Hillenbrand, and Javier Aizpurua. Probing the electromagnetic response of dielectric antennas by vortex electron beams. *Physical Review Research*, 5(2):023192, June 2023.

[25] D. N. Basov, M. M. Fogler, and F. J. García de Abajo. Polaritons in van der Waals materials. *Science*, 354(6309):aag1992, October 2016.

[26] Nicholas Rivera and Ido Kaminer. Light-matter interactions with photonic quasiparticles. *Nature Reviews Physics*, 2(10):538–561, October 2020.

[27] Ofer Kfir, Valerio Di Giulio, F. Javier García de Abajo, and Claus Ropers. Optical coherence transfer mediated by free electrons. *Science Advances*, 7(18):eabf6380, April 2021.

[28] F. J. García de Abajo. Optical excitations in electron microscopy. *Reviews of Modern Physics*, 82(1):209–275, February 2010.

[29] Scott D Findlay, Mark P Oxley, and Leslie J Allen. Modeling Atomic-Resolution Scanning Transmission Electron Microscopy Images. *Microscopy and Microanalysis*, 14(1):48–59, February 2008.

[30] C. Dwyer. Multislice theory of fast electron scattering incorporating atomic inner-shell ionization. *Ultramicroscopy*, 104(2):141–151, September 2005.

[31] W. Coene and D. Van Dyck. Inelastic scattering of high-energy electrons in real space. *Ultramicroscopy*, 33(4):261–267, October 1990.

[32] A. Lubk and J. Rusz. Jacob's ladder of approximations to paraxial dynamic electron scattering. *Physical Review B*, 92(23):235114, December 2015.

[33] Bernard Diu and Franck Laloë Claude Cohen-Tannoudji. *Cohen-Tannoudji, Diu and Laloë - Quantum Mechanics (Vol. I, II and III, 2nd Ed.).* December 2019.

[34] G. Nienhuis and L. Allen. Paraxial wave optics and harmonic oscillators. *Physical Review A*, 48(1):656–665, July 1993.

[35] G. M. Vanacore, G. Berruto, I. Madan, E. Pomarico,
P. Biagioni, R. J. Lamb, D. McGrouther, O. Reinhardt, I. Kaminer, B. Barwick, H. Larocque, V. Grillo, E. Karimi, F. J. García de Abajo, and F. Carbone. Ultrafast generation and control of an electron vortex beam via chiral plasmonic near fields. *Nature Materials*, 18(6):573–579, June 2019.

[36] Ulrich Hohenester and Andreas Trügler. MNPBEM – A Matlab toolbox for the simulation of plasmonic nanoparticles. *Computer Physics Communications*, 183(2):370–381, February 2012.

[37] Ofer Kfir. Entanglements of Electrons and Cavity Photons in the Strong-Coupling Regime. *Physical Review Letters*, 123(10):103602, September 2019.

[38] Andrea Konečná, Fadil Iyikanat, and F. Javier García de Abajo. Entangling free electrons and optical excitations. *Science Advances*, 8(47):eabo7853, November 2022.

[39] Armin Feist, Guanhao Huang, Germaine Arend, Yujia Yang, Jan-Wilke Henke, Arslan Sajid Raja, F. Jasmin Kappert, Rui Ning Wang, Hugo Lourenço-Martins, Zheru Qiu, Junqiu Liu, Ofer Kfir, Tobias J. Kippenberg, and Claus Ropers. Cavity-mediated electron-photon pairs. *Science*, 377(6607):777–780, August 2022.

[40] Germaine Arend, Guanhao Huang, Armin Feist, Yujia Yang, Jan-Wilke Henke, Zheru Qiu, Hao Jeng, Arslan Sajid Raja, Rudolf Haindl, Rui Ning Wang, Tobias J. Kippenberg, and Claus Ropers. Electrons herald nonclassical light. *September 2024*.

[41] G. S. Agarwal and K. Tara. Nonclassical properties of states generated by the excitations on a coherent state. *Physical Review A*, 43(1):492–497, January 1991.

[42] Axel Lubk. Chapter Two - Paraxial Quantum Mechanics. In Peter W. Hawkes, editor, *Advances in Imaging and Electron Physics*, volume 206, pages 15–58. Elsevier, January 2018.

[43] Arthur Losquin and Mathieu Kociak. Link between Cathodoluminescence and Electron Energy Loss Spectroscopy and the Radiative and Full Electromagnetic Local Density of States. *ACS Photonics*, 2(11):1619–1627, November 2015.

[44] Henrik Bruus and Karsten Flensberg. *Many-Body Quantum Theory in Condensed Matter Physics: An Introduction*. Oxford University Press Oxford, September 2004.

[45] Franck Aguilar, Hugo Lourenço-Martins, Damián Montero, Xiaoyan Li, Mathieu Kociak, and Alfredo Campos. Selective Probing of Longitudinal and Transverse Plasmon Modes with Electron Phase-Matching. *The Journal of Physical Chemistry C*, 127(45):22252–22264, November 2023.

[46] Giulio Guzzinati, Armand Béché, Hugo Lourenço-Martins, Jérôme Martin, Mathieu Kociak, and Jo Verbeeck. Probing the symmetry of the potential of localized surface plasmon resonances with phase-shaped electron beams. *Nature Communications*, 8(1):14999, April 2017.

[47] C. Fabre and N. Treps. Modes and states in quantum optics. *Reviews of Modern Physics*, 92(3):035005, September 2020.

[48] E. Karimi, R. W. Boyd, P. de la Hoz, H. de Guise, J. Řeháček, Z. Hradil, A. Aiello, G. Leuchs, and L. L. Sánchez-Soto. Radial quantum number of Laguerre-Gauss modes. *Physical Review A*, 89(6):063813, June 2014.

[49] M. W. Beijersbergen, L. Allen, H. E. L. O. van der Veen, and J. P. Woerdman. Astigmatic laser mode converters and transfer of orbital angular momentum. *Optics Communications*, 96(1):123–132, February 1993.

[50] L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman. Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes. *Physical Review A*, 45(11):8185–8189, June 1992.

[51] E. Abramochkin and V. Volostnikov. Beam transformations and nontransformed beams. *Optics Communications*, 83(1):123–135, May 1991.
# Contents

<table><tr><td>Acknowledgments</td><td>5</td></tr><tr><td>Phase shaped electron energy loss spectroscopy</td><td>9</td></tr><tr><td>EELS approximations</td><td>9</td></tr><tr><td>Electron description</td><td>9</td></tr><tr><td>Interaction probability for general electron beam</td><td>10</td></tr><tr><td>Interaction Hamiltonian: Multipolar expansion</td><td>12</td></tr><tr><td>Projected potential: intuitive interpretation</td><td>12</td></tr><tr><td>Mutipolar expansion</td><td>12</td></tr><tr><td>Multipolar EELS</td><td>12</td></tr><tr><td>Gaussian order EELS</td><td>12</td></tr><tr><td>Field operators</td><td>14</td></tr><tr><td>Transverse electron wavefunction description</td><td>15</td></tr><tr><td>Harmonic oscillator reminders</td><td>15</td></tr><tr><td>Ladder operator to build wavefunctions</td><td>16</td></tr><tr><td>Cartesian : Hermite-Gauss</td><td>16</td></tr><tr><td>Circular : Laguerre-Gauss</td><td>16</td></tr><tr><td>Basis interconnectivity</td><td>18</td></tr><tr><td>General transverse wavefunction expression</td><td>19</td></tr><tr><td>Order of a pure transverse state</td><td>19</td></tr><tr><td>Transition order between pure states</td><td>19</td></tr><tr><td>Multipolar interaction potential under ladder operator formalism</td><td>19</td></tr><tr><td>Interaction potential with operator algebra</td><td>19</td></tr><tr><td>Dipolar term</td><td>20</td></tr><tr><td>Dipolar potential : Hermite-Gauss</td><td>20</td></tr><tr><td>Dipolar potential : Laguerre-Gauss</td><td>20</td></tr><tr><td>Transition dipole moments</td><td>21</td></tr><tr><td>Quadrupolar term</td><td>21</td></tr><tr><td>Quadrupolar potential : Hermite-Gauss</td><td>22</td></tr><tr><td>Quadrupolar potential : Laguerre-Gauss</td><td>22</td></tr><tr><td>Higher orders</td><td>22</td></tr><tr><td>Shaped EELS interaction computation scheme</td><td>23</td></tr><tr><td>Transition amplitudes and probability</td><td>23</td></tr><tr><td>Computation example</td><td>23</td></tr><tr><td>Electrostatic term</td><td>24</td></tr><tr><td>Dipolar Hermite transitions</td><td>25</td></tr><tr><td>Dipolar Laguerre transitions</td><td>25</td></tr><tr><td>Quadrupolar Hermite transitions</td><td>25</td></tr><tr><td>Quadrupolar Laguerre transitions</td><td>26</td></tr><tr><td>Transition amplitudes to Gaussian final state</td><td>26</td></tr><tr><td>Higher k orders of potential impact on lower |i - j| order transitions amplitudes</td><td>26</td></tr><tr><td>Multi mode EELS</td><td>27</td></tr><tr><td>Identifying physical quantities measured under these techniques</td><td>27</td></tr><tr><td>Mapping</td><td>27</td></tr><tr><td>Another point of view : maps and spectra</td><td>27</td></tr></table>
Phase shaped electron energy loss spectroscopy

Accounting for the transverse profile changes of the electron beam through an energy loss interaction enriches the EELS experiment. The more complex the probe, the finer details we get. The goal of this section is to express properly the EELS interaction probability taking into account the transverse phase profile of the electron and its change during interaction.

EELS approximations

In the EELS scheme, me consider a fast electron beam passing nearby a target (atom, molecule, dielectric or plasmonic structure...). The target presents a nano-localized optical field represented by a set of quantum states {|n⟩}n∈N. It can be excited by exchanging a virtual photon with the fast electron. As a consequence we describe the target as a polariton field [25–27], and generically denotes its states as *photonic states*. To describe the system we apply a set of common approximations.

**Quasi-static:** By comparing the typical extension of the system’s optical field *D* to the wavelength of the interaction photon λ = 2πc ω , we can neglect the retardation effect, this is the **Quasi-static** approximation :

$$
\frac{\omega D}{c} \ll 1 \tag{18}
$$

With *c*, the celerity of light and *ω* the angular frequency corresponding to the probed energy.

**Single scattering:** EELS scattering amplitudes are weak, it can be considered that the electron interacts once at maximum with the sample during its full propagation, this is the first Born approximation

$$
\hat{U} \approx 1 - \frac{i}{\hbar} \int \hat{V}(t) dt \qquad (19)
$$

With $\hat{U}$ the system (electron + target) time evolution operator and $\hat{V}$ the Coulomb potential responsible for the interaction. $\hbar$ is the reduced Planck constant.

**Non-recoil:** In a transmission electron microscope (TEM) where these experiments are done, the typical speeds of the electron v/c ~ 0.5 – 0.8 enable us to neglect the lateral momentum transfers, this is the **non-recoil** approximation

$$
q_z = k_f - k_i \approx \frac{\omega}{v} \tag{20}
$$

With $k_i$ and $k_f$ the z-component of the wavevector for respectively the initial and final state of the electron and $q_z$ the transferred momentum.

**Paraxial approximation** : TEM experiments involve small angles, justifying the use of the paraxial approximation [42] enabling to decouple the transverse and longitudinal part of the electron wavefunction

$$
\psi(\mathbf{r}) \approx \frac{1}{\sqrt{L}} e^{ikz} \Psi_{\perp}(\mathbf{R}) \qquad (21)
$$

With *L* the quantization length for the longitudinal wavefunction. In Dirac notation

$$
|\psi_i\rangle = |\mathbf{k}_i^z\rangle \otimes |\Psi_{i,\perp}\rangle \qquad (22)
$$

In this study we consider a **non diffracting** electron beam, the transverse wavefunction is unchanged during propagation.

Electron description

Following these approximations, the electron beam is described by a wavefunction separated in a transverse and longitudinal part as described by (22). In Dirac notation it is the outer product of a transverse state and a longitudinal state: under the different approximations these two subspaces do not mix during free propagation.

$$
|\psi_i\rangle = |k_i\rangle \otimes |\Psi_{i,\perp}\rangle \qquad (23)
$$
## Interaction probability for general electron beam

We consider the system composed of a target or photonic cavity and a free electron. EELS can be described as a quantum process of transition between two states of the system. With $|i\rangle = |\psi_i, 0\rangle$ the initial state and $|f\rangle = |\psi_f, n\rangle$ a final excited state. Respectively composed of the electron state $|\psi_i\rangle$ or $|\psi_f\rangle$ and the photonic states of the target $|0\rangle$ or $|n\rangle$. We note the respective energies of these states: $\hbar\epsilon_i$, $\hbar\epsilon_f$, $\hbar\omega_0$, and $\hbar\omega_n$. From the Fermi's Golden rule, the EELS interaction rate reads

$$ \frac{d\Gamma^{\text{EELS}}}{dt} = \frac{2\pi}{\hbar} \sum_f \left|\langle f | \hat{V} | i \rangle\right|^2 \delta(\hbar\epsilon_i - \hbar\epsilon_f + \hbar\omega_0 - \hbar\omega_n) \quad (24) $$

Considering the propagation of the electron over a distance L during a time T, the interaction probability reads

$$ \Gamma^{\text{EELS}} = \frac{2\pi T}{\hbar^2} \sum_f \left|\langle f | \hat{V} | i \rangle\right|^2 \delta(\epsilon_i - \epsilon_f + \omega_0 - \omega_n) \quad (25) $$

With T the electron interaction time on the target.

$$ \Gamma^{\text{EELS}} = \frac{2\pi T}{\hbar^2} \sum_f \int d\omega \left|\langle f | \hat{V} | i \rangle\right|^2 \delta(\epsilon_i - \epsilon_f - \omega) \delta(\omega_0 - \omega_n - \omega) \quad (26) $$

Initially the photonic state is empty of any excitation and the electron wavefunction is that of a fast electron in a microscope, eventually with a shaped transverse profile. The final state is obtained after a virtual photon emitted by the electron has populated a photonic state of the target. Hence during the process, the electron loses a quantum of energy, its transverse profile can change and a photonic state $|n\rangle$ for the target is populated by one excitation. The interaction leads to a superposition of final states corresponding to the different interaction paths. The interaction probability for each final state can be retrieved by sorting the outgoing electron beam and selecting the right state. This is done via an electron spectrometer for energy resolution (fixed $\omega$) and post-selecting the transverse state with a state sorter for the electron beam.

$$ \Gamma^{\text{EELS}} = \frac{2\pi T}{\hbar^2} \int d\omega \sum_{\psi_f, n} \left|\langle \psi_f, n | \hat{V} | \psi_i, 0 \rangle\right|^2 \delta(\epsilon_i - \epsilon_f - \omega) \delta(\omega_0 - \omega_n - \omega) \quad (27) $$

We evaluate the interaction for a specific energy. In consequence of the paraxial approximation we separate the transverse and longitudinal parts of the electron wavefunction. To simplify notations we adopt the convention $|a\rangle \otimes |b\rangle = |a\rangle|b\rangle = |a, b\rangle$

$$ \Gamma^{\text{EELS}}(\omega) = \frac{2\pi T}{\hbar^2} \sum_{k_f, \Psi_{f,\perp}, n} \left|\langle k_f | \langle \Psi_{f,\perp}, n | \hat{V} | \Psi_{i,\perp}, 0 \rangle | k_i \rangle\right|^2 \delta(\epsilon_i - \epsilon_f - \omega) \delta(\omega_0 - \omega_n - \omega) \quad (28) $$

The longitudinal states for the free electron are expressed in a plane-wave basis. The sum can hence be made into an integral

$$ \sum_{k_f} = \frac{L}{2\pi} \int dk_f \quad (29) $$

Using the completeness relation $1 = \int dz|z\rangle\langle z|$ and the longitudinal state position representation $\langle z|k_f\rangle = \frac{1}{\sqrt{L}}e^{ik_f z}$

$$ \Gamma^{\text{EELS}}(\omega) = \frac{LT}{\hbar^2} \sum_{n, \Psi_{f,\perp}} \int dk_f \left| \iint dzdz' \langle k_f | z \rangle \langle \Psi_{f,\perp}, n | \langle z | \hat{V} | z \rangle | \Psi_{i,\perp}, 0 \rangle \langle z' | k_i \rangle \right|^2 \delta(\epsilon_i - \epsilon_f - \omega) \delta(\omega_0 - \omega_n - \omega) \quad (30) $$

$$ \Gamma^{\text{EELS}}(\omega) = \frac{T}{\hbar^2 L} \sum_{n, \Psi_{f,\perp}} \int dk_f \left| \langle \Psi_{f,\perp}, n | \iint dzdz' e^{-i(k_f z - k_i z')} \langle z | \hat{V} | z \rangle | \Psi_{i,\perp}, 0 \rangle \right|^2 \delta(\epsilon_i - \epsilon_f - \omega) \delta(\omega_0 - \omega_n - \omega) \quad (31) $$

In the narrow electron beam limit, we consider that the potential is local. $\langle z | \hat{V} | z' \rangle = \hat{V}(z, z') = \delta(z - z') \hat{V}(z)$. Hence
$$ \Gamma^{\text{EELS}}(\omega) = \frac{1}{\hbar^2 v} \sum_{n, \Psi_{f,\perp}} \int dk_f \left| \langle \Psi_{f,\perp}, n | \int dz e^{-i(k_f - k_i)z} \langle z | \hat{V} | z \rangle | \Psi_{i,\perp}, 0 \rangle \right|^2 \delta(\epsilon_i - \epsilon_f - \omega) \delta(\omega_0 - \omega_n - \omega) \quad (32) $$

With $v = \frac{L}{T}$

$$ \Gamma^{\text{EELS}}(\omega) = \frac{1}{\hbar^2 v} \sum_{n, \Psi_{f,\perp}} \int dk_f \iint dz dz' e^{-i(k_f - k_i)(z-z')} \langle \Psi_{i,\perp}, 0 | \langle z' | \hat{V}^\dagger | z' \rangle | \Psi_{f,\perp}, n \rangle \\ \times \langle \Psi_{f,\perp}, n | \langle z | \hat{V} | z \rangle | \Psi_{i,\perp}, 0 \rangle \delta(\epsilon_i - \epsilon_f - \omega) \delta(\omega_0 - \omega_n - \omega) \quad (33) $$

With $q_z = k_f - k_i$ and $\epsilon_i - \epsilon_f = -q_z v$, we obtain

$$ \int dk_f e^{-i(k_f - k_i)(z - z')} \delta(\epsilon_i - \epsilon_f - \omega) = \frac{1}{v} e^{i \frac{\omega}{v} (z - z')} \quad (34) $$

The expression simplifies to

$$ \Gamma^{\text{EELS}}(\omega) = \frac{1}{\hbar^2 v^2} \sum_{n, \Psi_{f,\perp}} \left| \langle \Psi_{f,\perp}, n | \int dz e^{i \frac{\omega}{v} z} \langle z | \hat{V} | z \rangle | \Psi_{i,\perp}, 0 \rangle \right|^2 \delta(\omega_0 - \omega_n - \omega) \quad (35) $$

Where we identify the projected potential [29–32], z-Fourier transform of the interaction potential.

$$ \hat{V}_{\text{proj}}(q_z, \omega) = \int dz \langle z | \hat{V} | z \rangle e^{iq_z z} \quad (36) $$

$$ \Gamma^{\text{EELS}}(\omega) = \frac{1}{\hbar^2 v^2} \sum_{n, \Psi_{f,\perp}} \left| \langle \Psi_{f,\perp} | \langle n | \hat{V}_{\text{proj}} | 0 \rangle | \Psi_{i,\perp} \rangle \right|^2 \delta(\omega_0 - \omega_n - \omega) \quad (37) $$

The final electron beam is in a superposition of transverse profiles. Post-selecting the electron transverse state to a single state leads to

$$ \Gamma^{\text{PSEELS}}(\omega) = \frac{1}{\hbar^2 v^2} \sum_n \left| \langle \Psi_{f,\perp} | \langle n | \hat{V}_{\text{proj}} | 0 \rangle | \Psi_{i,\perp} \rangle \right|^2 \delta(\omega_0 - \omega_n - \omega) \quad (38) $$

We now call this type of experiments Phase Shaped EELS (PSEELS). This formula contains the contribution of every photonic states at energy $\hbar\omega$. The non-zero contributions are those for which the transition moment for the photonic state matches the transition moment of the free electron transverse profile. This leads to a huge selectivity upon proper choice of initial and final states for the electron. For brevity, we name $\omega_{0n} = \omega_0 - \omega_n$ and encompass the energy conservation in a new sum definition

$$ \sum_n \delta(\omega_0 - \omega_n - \omega) = \sum_{\substack{n \\ \omega = \omega_{0n}}} = \tilde{\sum}_n \quad (39) $$

$$ \Gamma^{\text{PSEELS}}(\omega) = \tilde{\sum}_n \left| \frac{1}{\hbar v} \langle \Psi_{f,\perp} | \langle n | \hat{V}_{\text{proj}} | 0 \rangle | \Psi_{i,\perp} \rangle \right|^2 \quad (40) $$

$$ \Gamma^{\text{PSEELS}}(\omega) = \tilde{\sum}_n \left| \frac{1}{\hbar v} \langle \Psi_{f,\perp} | \langle n | \hat{V}_{\text{proj}} | 0 \rangle | \Psi_{i,\perp} \rangle \right|^2 = \tilde{\sum}_n |g_{0,n}|^2 \quad (41) $$

With $g_{0,n}$ the coupling constant. The next steps to study phase shaped interactions is to describe the transverse profiles of the electron beam and to express the projected potential. We will see in the following that the transverse profiles can be expressed through a ladder operator formalism and the projected potential can be made into a multipolar development.
# Interaction Hamiltonian: Multipolar expansion

## Projected potential: intuitive interpretation

The EELS measurement is done after the electron propagation on a two dimensional detector, effectively summing the contributions along the z coordinate modulated by the electron phase. The interaction potential summed this ways leads to the projected potential, z-Fourier transform of the interaction potential, which emerges naturally from the EELS computation. It contains all the physical information of the interaction and accounts for the delocalized nature of the interaction.

## Multipolar expansion

The electron beam waist $w_0$ (picometer) is small compared to the electrostatic potential variations scale (nanometer). Therefore the potential seen by the electron $\langle n|\hat{V}_{\text{proj}}|0\rangle = \hat{V}_{\text{proj}}(\mathbf{R}_0, q_z, \omega)$ is considered constant over the waist of the electron beam, the same applies for its derivatives and thus the projected potential can be expressed via a multipolar development (equivalent to the Power-Zienau-Woolley form of the minimal-coupling) [1].

$$ \hat{V}_{\text{proj}}(q_z, \omega) = \underbrace{\hat{\rho}\hat{\Phi}(q_z, \omega)}_{\hat{V}_0} - \underbrace{\hat{d} \cdot \hat{\mathbf{E}}(q_z, \omega)}_{\hat{V}_1} - \underbrace{\hat{\vec{Q}} \cdot \nabla \cdot \hat{\mathbf{E}}(q_z, \omega)}_{\hat{V}_2} + \dots + \text{Magnetic terms} \quad (42) $$

Each term is composed of an operator acting on the electron transverse degree of freedom and an operator acting on the photonic degrees of freedom. If specified to a plasmonic structure, the photonic operator will lead to the electrostatic potential and its derivatives.

## Multipolar EELS

The PSEELS coupling constant as seen in (41) can be expressed for each photonic state $|n\rangle$ with the help of (42) as a multipolar EELS coupling constant. Each contribution will give informations on the physics of the system.

$$ \begin{aligned} g_{0,n}(q_z, \omega) &= \frac{1}{\hbar v} \langle \Psi_{f,\perp} | \langle n | \hat{\rho} \hat{\Phi} - \hat{d} \cdot \hat{\mathbf{E}} + \dots | 0 \rangle | \Psi_{i,\perp} \rangle \\ &= g_{0,n}^{(0)} + g_{0,n}^{(1)} + g_{0,n}^{(2)} + \dots \end{aligned} \quad (43) $$

As a consequence, The PSEELS probability (40) is also expressed as a multipolar probability

$$ \begin{aligned} \Gamma^{\text{PSEELS}}(q_z, \omega) &= \sum_n^{\sim} \left| \frac{1}{\hbar v} \langle \Psi_{f,\perp} | \langle n | \hat{\rho} \hat{\Phi} - \hat{d} \cdot \hat{\mathbf{E}} + \dots | 0 \rangle | \Psi_{i,\perp} \rangle \right|^2 \\ &\approx \Gamma^{(0)} + \Gamma^{(1)} + \Gamma^{(2)} + \dots \end{aligned} \quad (44) $$

Given the narrow beam condition needed for the multipolar development, the cross terms are either zero or orders of magnitude below the previous PSEELS order. We can neglect them to obtain this elegant order by order EELS expression.

## Gaussian order EELS

We show here that the interaction due to the lowest order of the potential multipolar development is equivalent to classical EELS. We consider an initial and final Gaussian electron beam transverse profile, roughly equivalent to non shaped EELS. In a simplified vision, considering only the lowest order of the interaction

$$ \Gamma^{(0)}(q_z, \omega) = \sum_n^{\sim} \left| \frac{1}{\hbar v} \langle \Psi_{f,\perp} | \langle n | \hat{\rho} \hat{\Phi} | 0 \rangle | \Psi_{i,\perp} \rangle \right|^2 \quad (45) $$
The $\hat{\rho}$ operator acts on the transverse degree of freedom of the electron and the $\hat{\phi}$ acts on the photonic state only. They can be separated as follows

$$ \Gamma^{(0)}(q_z, \omega) = \sum_n^c \left| \frac{1}{\hbar v} \langle \Psi_{f,\perp} | \hat{\rho} | \Psi_{i,\perp} \rangle \langle n | \hat{\Phi} | 0 \rangle \right|^2 \quad (46) $$

We choose the initial and final transverse states as Gaussian states $|G\rangle$ which wavefunction has a waist $w_0$ and is centered at $\mathbf{R}_0$.

$$ \Gamma^{(0)}(q_z, \omega) = \sum_n^c \left| \frac{e}{\hbar v} \int d\mathbf{R} \Psi_{f,\perp}^*(\mathbf{R}) \Psi_{i,\perp}(\mathbf{R}) \langle n | \hat{\Phi}(\mathbf{R}) | 0 \rangle \right|^2 \quad (47) $$

$$ \Psi_{f,\perp}(\mathbf{R}) = \Psi_{i,\perp}(\mathbf{R}) = \langle \mathbf{R} | G \rangle = G(\mathbf{R}) = \frac{1}{w_0} \sqrt{\frac{2}{\pi}} e^{-(\mathbf{R}-\mathbf{R}_0)^2/w_0^2} \quad (48) $$

Under the narrow beam condition, where the typical length of variation for the potential is great compared to the beam waist, the transverse profile squared norm tends to a dirac delta function.

$$ \left|G(\mathbf{R})\right|^2 \underset{w_0 \to 0}{\longrightarrow} \delta(\mathbf{R} - \mathbf{R}_0) \quad (49) $$

Hence the electrostatic potential is probed at the beam position.

$$ \Gamma^{(0)}(q_z, \omega) = \sum_n^c \left| \frac{e}{\hbar v} \langle n | \hat{\Phi}(\mathbf{R}_0) | 0 \rangle \right|^2 \quad (50) $$

This expression is in fact classical EELS, the following computations are aimed at retrieving the link to the local density of states (LDOS) as in [43].

$$ \Gamma^{(0)}(q_z, \omega) = \frac{e^2}{\hbar^2 v^2} \sum_n \langle 0 | \hat{\Phi}^{\dagger}(\mathbf{R}_0) | n \rangle \langle n | \hat{\Phi}(\mathbf{R}_0) | 0 \rangle \quad (51) $$

We "un-project" the potential

$$ \Gamma^{(0)}(q_z, \omega) = \frac{e^2}{\hbar^2 v^2} \iint dzdz' e^{iq_z(z'-z)} \sum_n \langle 0 | \hat{\Phi}^{\dagger}(\mathbf{R}_0, z) | n \rangle \langle n | \hat{\Phi}(\mathbf{R}_0, z') | 0 \rangle \quad (52) $$

The electrostatic potential is hermitian, hence $\Phi^\dagger = \Phi$. Using the formula for the electrostatic propagator in the framework of Kubo's linear response theory [44]

$$ \mathrm{Im}\{W(r, r', \omega)\} = -\frac{\pi}{\hbar} \sum_n \langle 0 | \phi(\mathbf{r}) | n \rangle \langle n | \phi(\mathbf{r}') | 0 \rangle \quad (53) $$

$$ \Gamma^{(0)}(q_z, \omega) = \frac{e^2}{\hbar \pi v^2} \iint dz dz' e^{iq_z(z'-z)} \mathrm{Im}\{-W(\mathbf{R}_0, \mathbf{R}_0, z, z', \omega)\} \quad (54) $$

$$ q_z = \frac{\omega}{v} $$

$$ \Gamma^{(0)}(q_z, \omega) = \frac{e^2}{\hbar \pi \omega^2} \iint dz dz' q_z^2 e^{iq_z(z'-z)} \mathrm{Im}\{-W(\mathbf{R}_0, \mathbf{R}_0, z, z', \omega)\} \quad (55) $$

By double integration by part along $z$ and $z'$.
$$ \Gamma^{(0)}(q_z, \omega) = \frac{e^2}{\hbar\pi\omega^2} \iint dzdz' e^{iq_z(z'-z)} \mathrm{Im}\{-\partial_z\partial_{z'}W(\mathbf{R}_0, \mathbf{R}_0, z, z', \omega)\} \quad (56) $$

Leading to the identification of the Green dyad

$$ G_{zz} = \frac{1}{4\pi\omega^2} \partial_z \partial_{z'} W(\mathbf{r}, \mathbf{r}') \quad (57) $$

$$ \Gamma^{(0)}(q_z, \omega) = \frac{4e^2}{\hbar} \iint dzdz' e^{iq_z(z'-z)} \mathrm{Im}\{-G_{zz}(\mathbf{R}_0, \mathbf{R}_0, z, z', \omega)\} \quad (58) $$

Identifying the local density of states (LDOS) [43]

$$ \rho_{\mu}(\mathbf{r}, \mathbf{r}, \omega) = -\frac{2\omega}{\pi} \mathrm{Im} \left\{ \mu \cdot \overleftrightarrow{G}(\mathbf{r}, \mathbf{r}, \omega) \cdot \mu \right\} \quad (59) $$

$$ \Gamma^{(0)}(\mathbf{R}_0, q_z, \omega) = \frac{2\pi e^2}{\hbar\omega} \iint dzdz' e^{iq_z(z'-z)} \rho_z(\mathbf{R}_0, \mathbf{R}_0, z, z', \omega) \quad (60) $$

We retrieve the standard EELS probability being proportional to the LDOS.

$$ \Gamma^{(0)}(q_z, \omega) = \frac{2\pi e^2}{\hbar\omega} \rho_z(\mathbf{R}_0, \mathbf{R}_0, q_z, -q_z, \omega) \quad (61) $$

A transition from a gaussian state to a gaussian state is similar to the standard non shaped EELS experiment. The $\Gamma^{(0)}$ contribution is the leading order. Contributions from higher orders can be non-zero as well but will be shown later to be orders of magnitude smaller. Indeed, standard EELS experiment consists of a gaussian incoming beam and no selection for the outgoing beam. The dominant term in probability will come from the Gaussian part of the outgoing beam.

## Field operators

The field operators acting on the nano-optical field states give the expectation value $\langle n|\hat{A}|0\rangle$, (with $\hat{A}$ a field operator e.g $\hat{E}$ or $\hat{\Phi}$) proportional to the classical fields. For simpler notations we write

$$ \langle n|\hat{A}|0\rangle = A_{n0} \quad (62) $$

For instance, if the target is specified to be a plasmonic structure, then these average values of the photonic operators will lead to the classical electrostatic potential and its derivatives. Indeed, while equation (41) involves a sum over all the possible quantum states of the target $|n\rangle$, it can be straightforwardly replaced by a sum over all the possible modes $m$ in the case of a classical system:

$$ \sum_{\text{states } n} |g_{0,n}|^2 \xrightarrow[\text{classical field}]{\quad} \sum_{\text{modes } m} \beta_m |g_{0,m}|^2 \quad (63) $$

where $g_{0,m}$ denotes the scattering amplitude per mode [45] and is a $\beta_m$ a scaling factor which represents the number of states per mode. In the widespread case of classical plasmonic field [9, 10, 23, 46], $\beta_m = 4\hbar/\pi \mathrm{Im}\{-f_m\}$ with $f_m$ the modal spectral function [45]. Although the replacement rule (63) appears as a handy proxy ensuring a simple translation of our formalism between the quantum and classical regimes, it shall be kept in mind that modes and states are two completely different quantum optical concepts [47] which should not be confused.

In later developments, we will focus on a single $n$ photonic mode of the system and omit the $n0$ subscript for clarity. PSEELS becomes

$$ \Gamma^{\text{PSEELS}}(q_z, \omega) = \sum_n \left| \frac{1}{\hbar v} \langle \Psi_{f,\perp} | \hat{\rho} \Phi_{n0}(\mathbf{R}_0, q_z, \omega) - \hat{d} \cdot \mathbf{E}_{n0}(\mathbf{R}_0, q_z, \omega) + \dots | \Psi_{i,\perp} \rangle \right|^2 \quad (64) \\ \cong \Gamma^{(0)} + \Gamma^{(1)} + \Gamma^{(2)} + \dots $$
Transverse electron wavefunction description

We study transitions of the transverse electron wave-function between chosen states $|\Psi_{i,\perp}\rangle$ and $|\Psi_{f,\perp}\rangle$. Being a two dimensional problem, we can decompose every transverse wavefunction on the basis of states related to the two dimensional harmonic oscillator, inspired by [34]. This brings in the ladder operator formalism from the standard textbook problem to express arbitrary wavefunction shape as a combination of creation and annihilation operators acting on a fundamental state. In a first approach, we pre- and post-select on the pure harmonic oscillator states to later generalise to any state. We here make a brief reminder on this topic.

Harmonic oscillator reminders

Considering a massive particle in a two dimensional space subjected to a two dimensional harmonic potential

$$ \hat{V} = \frac{1}{2} m \omega^2 (\hat{x}^2 + \hat{y}^2) \qquad (65) $$

The hamiltonian of the system is

$$ \hat{H} = \hat{H}_{xy} = \hat{H}_x + \hat{H}_y = \frac{\hat{p}_x^2 + \hat{p}_y^2}{2m} + \frac{1}{2} m \omega^2 (\hat{x}^2 + \hat{y}^2) \qquad (66) $$

With $\hat{H}_x$ and $\hat{H}_y$ the 1D harmonic oscillator hamiltonian. To prepare for the electron description we use the notation

$$ \frac{\sqrt{2}}{w_0} = \sqrt{\frac{m\omega}{\hbar}} \qquad (67) $$

Where $w_0$ corresponds to the waist of the electron beam in the PSEELS study. The eigenvalues $|\chi\rangle$ of the hamiltonian are found using the time independent Schrödinger equation

$$ \hat{H}|\chi\rangle = E|\chi\rangle \qquad (68) $$

The x and y coordinates remain independent, thus the energy spectrum is the sum of two one dimension harmonic oscillator.

$$ E = \left(n_x + \frac{1}{2}\right) \hbar\omega + \left(n_y + \frac{1}{2}\right) \hbar\omega = (n_x + n_y + 1) \hbar\omega \qquad (69) $$

The eigenvalues of the hamiltonian are degenerate. Two different sets of operator are used to properly find and index the eigenstates: the "cartesian" $\{\hat{H}_{xy}, \hat{H}_x, \hat{H}_y\}$ or "angular momentum" $\{\hat{H}_{xy}, \hat{L}_z\}$ where $\hat{L}_z$ is the orbital angular momentum operator. Consequence of this fact is the existence of two complete sets of orthogonal eigenstates describing the system, the Hermite-Gauss $|HG_{n_x, n_y}\rangle$ for the "cartesian" and the circular $|LG_{n_\odot, n_\odot}\rangle$ for the angular momentum. The eigenvalue of $\hat{H}_{xy}$ is the number of quanta of vibration in the system: $n_x + n_y = n_\odot + n_\odot$; The eigenvalue of $\hat{L}_z$ is $n_\odot - n_\odot$, balance of left and right circular vibration quanta. The canonical way to study this system is to define the ladder operators in the Cartesian basis, creating or annihilating one quantum of vibration in the corresponding direction. In this system, acting on the x direction has no impact on the y direction.

$$ \hat{a}_x = \left( \frac{1}{w_0} \hat{x} + i \frac{w_0 \hat{p}_x}{2\hbar} \right) \qquad (70) $$

$$ \hat{a}_y = \left( \frac{1}{w_0} \hat{y} + i \frac{w_0 \hat{p}_y}{2\hbar} \right) \qquad (71) $$

The ladder operators obey the following commutation relation

$$ [\hat{a}_i, \hat{a}_j^\dagger] = \delta_{i,j} \mathbf{1} \qquad (72) $$

The system's hamiltonian can be expressed in term of these operators as :

$$ \hat{H} = \hbar\omega(1 + \hat{N}_x + \hat{N}_y) \qquad (73) $$

With the number operator

$$ \hat{N}_i = \hat{a}_i^\dagger \hat{a}_i \qquad (74) $$
## Ladder operator to build wavefunctions

Eigenstates of the system can be built from the fundamental eigenstate $|0\rangle = |G\rangle$ which is the 2D gaussian state [33]. Two sets of eigenstates can be used to describe properly the system, the Hermite-Gauss states for cartesian coordinates and the circular states for the "angular momentum" coordinates, the later can be mapped to Laguerre-Gauss states. Eigenstates in a given basis can be indexed by the number of each ladder operators needed to generate it starting from the fundamental gaussian mode

$$ \langle \mathbf{r} | G \rangle = G(x, y) = \frac{1}{w_0} \sqrt{\frac{2}{\pi}} e^{-(x^2+y^2)/w_0^2} \quad (75) $$

Which is normalized so that $\langle G|G \rangle = 1$.

## Cartesian : Hermite-Gauss

The first set of operator sharing a basis of eigenvectors is {$\hat{H}_{xy}, \hat{H}_x, \hat{H}_y$}. We name it the cartesian coordinate basis. Eigenstates are the Hermite-Gauss states $|HG_{n_x, n_y}\rangle$, the cartesian ladder operators act on these states following

$$ \begin{aligned} \hat{a}_x |HG_{n_x, n_y}\rangle &= \sqrt{n_x} |HG_{n_x-1, n_y}\rangle \\ \hat{a}_x^\dagger |HG_{n_x, n_y}\rangle &= \sqrt{n_x+1} |HG_{n_x+1, n_y}\rangle \\ \hat{a}_y |HG_{n_x, n_y}\rangle &= \sqrt{n_y} |HG_{n_x, n_y-1}\rangle \\ \hat{a}_y^\dagger |HG_{n_x, n_y}\rangle &= \sqrt{n_y+1} |HG_{n_x, n_y+1}\rangle \end{aligned} \quad (76) $$

Every HG state hence can be expressed from the Gaussian fundamental state

$$ |HG_{n_x, n_y}\rangle = \frac{1}{\sqrt{n_x! n_y!}} (\hat{a}_x^\dagger)^{n_x} (\hat{a}_y^\dagger)^{n_y} |G\rangle \quad (77) $$

Expressed as a wavefunction

$$ \langle \mathbf{r} | HG_{n_x, n_y} \rangle = HG_{n_x, n_y} (x, y) = \frac{H_{n_x} \left( \frac{x\sqrt{2}}{w_0} \right) H_{n_y} \left( \frac{y\sqrt{2}}{w_0} \right)}{w_0 \sqrt{2^{n_x+n_y-1} \pi n_x! n_y!}} e^{-(x^2+y^2)/w_0^2} \quad (78) $$

With $H_{n_x}(x)$ the Hermite polynomial of degree $n_x$. We here define the degree of the transverse electron wavefunction as the number of creation operators used to build it from the Gaussian state $n_x + n_y$. It is also the sum of the degrees of the Hermite polynomials describing it. Given this formula we can express every modes of this basis and we have a method to relate them to the fundamental gaussian state.

## Circular : Laguerre-Gauss

The second set of operator sharing a basis of eigenvectors is {$\hat{H}_{xy}, \hat{L}_z$}. The eigenstates have a well defined angular momentum. They can be mapped to the Laguerre Gauss modes. The orbital angular momentum operator reads:

$$ \hat{L} = \hat{R} \times \hat{P} \quad (79) $$

With $\hat{R}$ and $\hat{P}$ the position and momentum operators respectively. Projected on the z-axis:

$$ \hat{L}_z = \hat{x}\hat{p}_y - \hat{y}\hat{p}_x \quad (80) $$

Also expressed in the position representation ((r, φ) coordinates) as

$$ \langle \mathbf{r} | \hat{L}_z = \frac{\hbar}{i} \frac{\partial}{\partial \phi} \langle \mathbf{r} | \quad (81) $$
The z-orbital angular momentum operator can be expressed using cartesian ladder operators

$$ \hat{L}_z = i\hbar(\hat{a}_x\hat{a}_y^\dagger - \hat{a}_y^\dagger\hat{a}_x) \qquad (82) $$

We introduce $a_{\odot}^{\dagger}$ and $a_{\ominus}^{\dagger}$, that can be seen as "circular creation operators" [33]. Each one creates a quantum of circular vibration, that is, carrying one unit of angular momentum. The conjugate operator annihilates the corresponding quatum. We take the reversed convention of [33] as we define $\hat{a}_{\odot}^{\dagger}$ as the operator adding a positive quantum of orbital angular momentum to the system, e.g $\hat{L}_z \hat{a}_{\odot}^{\dagger}|G\rangle = +\hbar \hat{a}_{\odot}^{\dagger}|G\rangle$

$$ \begin{aligned} \hat{a}_{\odot}^{\dagger} &= \frac{1}{\sqrt{2}}(\hat{a}_{x}^{\dagger} + i\hat{a}_{y}^{\dagger}) \\ \hat{a}_{\ominus} &= \frac{1}{\sqrt{2}}(\hat{a}_{x} - i\hat{a}_{y}) \end{aligned} \qquad (83) $$

$$ \begin{aligned} \hat{a}_{\odot}^{\dagger} &= \frac{1}{\sqrt{2}}(\hat{a}_{x}^{\dagger} - i\hat{a}_{y}^{\dagger}) \\ \hat{a}_{\ominus} &= \frac{1}{\sqrt{2}}(\hat{a}_{x} + i\hat{a}_{y}) \end{aligned} \qquad (84) $$

With these we can rewrite the operators

$$ \hat{H}_{xy} = (\hat{a}_{\odot}^{\dagger}\hat{a}_{\odot} + \hat{a}_{\odot}^{\dagger}\hat{a}_{\ominus} + 1)\hbar\omega = (\hat{N}_{\odot} + \hat{N}_{\ominus} + 1)\hbar\omega \qquad (85) $$

$$ \hat{L}_z = (\hat{a}_{\odot}^{\dagger}\hat{a}_{\odot} - \hat{a}_{\odot}^{\dagger}\hat{a}_{\ominus})\hbar = (\hat{N}_{\odot} - \hat{N}_{\ominus})\hbar \qquad (86) $$

$\hat{H}_{xy}$ appears as an undistinguished count of the circular quanta in the system whereas $\hat{L}_z$ measures the balance of left and right quanta in the system. Under the present convention, a left rotating object carries a positive angular momentum. Acting on the Gaussian state of the system with these ladder operators gives new states with well-defined angular momenta as are the Laguerre-Gauss modes. We name them $|LG\rangle$, however the position of the indices brings a nuance.

$$ |LG_{n_{\odot}, n_{\ominus}}\rangle = \frac{(\hat{a}_{\odot}^{\dagger})^{n_{\odot}} (\hat{a}_{\odot}^{\dagger})^{n_{\ominus}}}{\sqrt{n_{\odot}! n_{\ominus}!}} |G\rangle \qquad (87) $$

Mapping to Laguerre-Gauss states is particularly interesting : they form a basis commonly used in the description of optical vortex beams and have a lot of useful representation tools. It can be shown that [33, 48, 49]

$$ |LG_{n_{\odot}, n_{\ominus}}\rangle = (-1)^p |LG_p^l\rangle \qquad (88) $$

With $l = (n_{\odot} - n_{\ominus}) \in \mathbb{Z}$, the orbital angular momentum quantum number and $p = \frac{1}{2}(n_{\odot} + n_{\ominus} - |n_{\odot} - n_{\ominus}|) = \min(n_{\odot}, n_{\ominus}) \in \mathbb{N}$ the radial quantum number.

The position representations of Laguerre-Gauss states are

$$ \langle r | LG_p^l \rangle = LG_p^l(r, \varphi) = \sqrt{\frac{2p!}{\pi(p+|l|)!}} \frac{1}{w_0} \left(\frac{\sqrt{2}r}{w_0}\right)^{|l|} L_p^{|l|} \left(\frac{2r^2}{w_0^2}\right) e^{-r^2/w_0^2} e^{il\varphi} \qquad (89) $$

The wavefunction is expressed in term of the generalized Laguerre polynomials $L_p^{|l|}(r)$. Note the subtle difference between these two family of states. The $|LG_{n_{\odot}, n_{\ominus}}\rangle$ straightforwardly expresses the wavefunction in term of number of left or right quanta and the Laguerre-Gauss is easier to represent since more commonly used, even if the $p$ quantum number is harder to interpret. These two sets of functions are identical to a $(-1)^p$ factor, this nuance is not exposed in the main text for brevity. The following table represents the correspondence of wavefunctions $|LG\rangle$ indexed by $n_{\odot}$ and $n_{\ominus}$ and $|LG\rangle$ indexed by $p$ and $l$ and represented as $(p, l)$ in the table
<table><thead><tr><th>n<sub>&#8455;</sub></th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>n<sub>&#8455;</sub></th></tr></thead><tbody><tr><th>0</th><td>(0,0)</td><td>(0-1)</td><td>(0,-2)</td><td>(0,-3)</td><td>(0,-4)</td><td></td></tr><tr><th>1</th><td>(0,1)</td><td>(1,0)</td><td>(1,-1)</td><td>(1,-2)</td><td></td><td></td></tr><tr><th>2</th><td>(0,2)</td><td>(1,1)</td><td>(2,0)</td><td>(2,-1)</td><td></td><td></td></tr><tr><th>3</th><td>(0,3)</td><td>(1,2)</td><td>(2,1)</td><td>(3,0)</td><td></td><td></td></tr><tr><th>4</th><td>(0,4)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>n<sub>&#8455;</sub></th><td></td><td></td><td></td><td></td><td></td><td>p = min(n<sub>&#8455;</sub>, n<sub>&#8455;</sub>)<br/>l = n<sub>&#8455;</sub> - n<sub>&#8455;</sub></td></tr></tbody></table>

Where $(p,l)$ functions have a global minus sign compared to their $n_{\phi}, n_{\phi}$ equivalent. Creating a left and a right quanta of excitation on the gaussian state leads to raising the $p$ values, hence adding a radial node to the wavefunction without changing the orbital angular momentum ("dark ring" in the intensity profile). The ladder operator rules are

$$ 
\begin{aligned}
\hat{a}_{\phi} |LG_{n_{\phi}, n_{\phi}}\rangle &= \sqrt{n_{\phi}} |LG_{n_{\phi}-1, n_{\phi}}\rangle \\
\hat{a}_{\phi} |LG_{n_{\phi}, n_{\phi}}\rangle &= \sqrt{n_{\phi}} |LG_{n_{\phi}, n_{\phi}-1}\rangle \\
\hat{a}_{\phi}^{\dagger} |LG_{n_{\phi}, n_{\phi}}\rangle &= \sqrt{n_{\phi}+1} |LG_{n_{\phi}+1, n_{\phi}}\rangle \\
\hat{a}_{\phi}^{\dagger} |LG_{n_{\phi}, n_{\phi}}\rangle &= \sqrt{n_{\phi}+1} |LG_{n_{\phi}, n_{\phi}+1}\rangle
\end{aligned}
\qquad (90) $$

## Basis interconnectivity

There is a direct mapping of LG to HG, discussed for instance in these articles [34, 49–51]. We adapt here the notations to our formalism. Defining $N = n_x + n_y = n_{\phi} + n_{\phi}$, the order of the transverse profile, and the usual $p = \min(n_{\phi}, n_{\phi})$ for LG

$$ |LG_{n_{\phi}, n_{\phi}}\rangle = (-1)^p \sum_{k=0}^{N} i^k b(n_{\phi}, n_{\phi}, k) |HG_{N-k, k}\rangle \qquad (91) $$

With

$$ b(n_{\phi}, n_{\phi}, k) = \left( \frac{(N-k)!k!}{2^N n_{\phi}! n_{\phi}!} \right)^{1/2} \frac{1}{k!} \frac{d^k}{dt^k} [(1-t)^{n_{\phi}} (1+t)^{n_{\phi}}]_{t=0} \qquad (92) $$

This relations shows that an order $N$ Laguerre-Gauss beam is expressed exclusively as a sum of order $N$ Hermite-Gauss beams. Meaning, in some sense, that there is a conserved quantity in the decomposition and every order is well separated. Interestingly there exists nearly the same relation linking HG along x,y axis and x',y' axis: tilted by 45°. It is Identical as for the LG but in the real space (no imaginary term):

$$ HG_{n,m} \left( \frac{x+y}{\sqrt{2}}, \frac{x-y}{\sqrt{2}} \right) = \sum_{k=0}^{N} b(n,m,k) HG_{N-k,k}(x,y) \qquad (93) $$

We can rewrite the HG to LG relation under the previously defined operator algebra for 2D functions.

$$ \left( \frac{(\hat{a}_{\phi}^{\dagger})^{n_{\phi}} (\hat{a}_{\phi}^{\dagger})^{n_{\phi}}}{\sqrt{n_{\phi}! n_{\phi}!}} \right) |G\rangle = (-1)^p \left[ \sum_{k=0}^{N} i^k b(n_{\phi}, n_{\phi}, k) \frac{(\hat{a}_x^{\dagger})^{(N-k)} (\hat{a}_y^{\dagger})^k}{\sqrt{(N-k)! k!}} \right] |G\rangle \qquad (94) $$

Which simplifies by expanding the $b$ term as:

$$ (\hat{a}_{\phi}^{\dagger})^{n_{\phi}} (\hat{a}_{\phi}^{\dagger})^{n_{\phi}} |G\rangle = \frac{(-1)^p}{\sqrt{2^N}} \left[ \sum_{k=0}^{N} \frac{i^k}{k!} \frac{d^k}{dt^k} [(1-t)^{n_{\phi}} (1+t)^{n_{\phi}}]_{t=0} (\hat{a}_x^{\dagger})^{(N-k)} (\hat{a}_y^{\dagger})^k \right] |G\rangle \qquad (95) $$
Hence we have the operator relation linking a product of creation operators in the LG basis to a product of creation operators in the HG basis.

$$ (\hat{a}_{\odot}^{\dagger})^{n_{\odot}} (\hat{a}_{\odot}^{\dagger})^{n_{\odot}} = \frac{(-1)^{min(n_{\odot}, n_{\odot})}}{\sqrt{2^N}} \left[ \sum_{k=0}^{N} \frac{i^k}{k!} \frac{d^k}{dt^k} [(1-t)^{n_{\odot}} (1+t)^{n_{\odot}}]_{t=0} (\hat{a}_x^{\dagger})^{(N-k)} (\hat{a}_y^{\dagger})^k \right] \quad (96) $$

At this point, electron transverse wavefunctions can be expressed on two different basis : The Hermite Gauss $|HG_{n_x,n_y}\rangle$ and the Laguerre-Gauss $|LG_{n_{\odot},n_{\odot}}\rangle$. Each element of these basis can be expressed as ladder operators acting on a gaussian state. The decomposition of Laguerre Gauss modes on Hermite-Gauss basis is a useful tool to express PSEELS experiments where the transition occurs between two states from different basis, for instance $|\Psi_{i,\perp}\rangle = |LG_{1,2}\rangle$, $|\Psi_{f,\perp}\rangle = |HG_{1,1}\rangle$.

## General transverse wavefunction expression

Any transverse electron beam can be expressed as a combination of $|HG_{n_x,n_y}\rangle$ or $|LG_{n_{\odot},n_{\odot}}\rangle$.

$$ \begin{aligned} |\Psi_{\perp}\rangle &= \sum_{n_x, n_y} c_{n_x, n_y} |HG_{n_x, n_y}\rangle \\ |\Psi_{\perp}\rangle &= \sum_{n_{\odot}, n_{\odot}} c_{n_{\odot}, n_{\odot}} |LG_{n_{\odot}, n_{\odot}}\rangle \end{aligned} \quad (97) $$

With

$$ \begin{aligned} c_{n_x, n_y} &= \langle HG_{n_x, n_y} | \Psi_{\perp} \rangle \\ c_{n_{\odot}, n_{\odot}} &= \langle LG_{n_{\odot}, n_{\odot}} | \Psi_{\perp} \rangle \end{aligned} \quad (98) $$

Which in turn can be expressed through ladder operators (77,87). In the following we will only study pure states (Laguerre- or Hermite- Gauss), generalization is then straightforward.

## Order of a pure transverse state

The order $i$ of a pure HG or LG transverse state is defined as the number of creation operators required to build it from the gaussian state. We have $i = n_x + n_y$ or $i = n_{\odot} + n_{\odot}$.

## Transition order between pure states

Any pure transverse state can be transformed into another pure transverse state by application of a product of ladder operators. Then, the order $|i-j|$ of the transition between two transverse states $i$ and $j$ is defined as the number of operators needed to transit from $i$ to $j$. A first order transition corresponds to $|HG_{1,0}\rangle \rightarrow |G\rangle$ or $|LG_{1,0}\rangle \rightarrow |LG_{1,1}\rangle$ for instance.

Rigorously, the order of the transition $|HG_{n_x,n_y}\rangle \rightarrow |HG_{n_{x'},n_{y'}}\rangle$ is defined as

$$ |i-j| = |n_{x'} - n_x| + |n_{y'} - n_y| = |\Delta n_x| + |\Delta n_y| \quad (99) $$

This subtlety in the definition of the transition order has an impact for instance on transition $|HG_{1,0}\rangle \rightarrow |HG_{0,1}\rangle$ where the simple difference $(i-j) = 0$ but the transition order $|i-j| = 2$. The definition is easily generalized to LG states.

## Multipolar interaction potential under ladder operator formalism

### Interaction potential with operator algebra

Any transverse electron beam profile can be expressed via a ladder operator combination, the same is done for the projected potential multipolar development (42) at every orders. Given the expressions for the ladder operators in the
two coordinate systems (70),(71), (83) and (84), the position operators can be re-expressed as

$$ \hat{x} = \frac{w_0}{2} (\hat{a}_x + \hat{a}_x^\dagger) = \frac{w_0}{2\sqrt{2}} (\hat{a}_\odot + \hat{a}_\odot + \hat{a}_\odot^\dagger + \hat{a}_\odot^\dagger) \quad (100) $$

$$ \hat{y} = \frac{w_0}{2} (\hat{a}_y + \hat{a}_y^\dagger) = \frac{i w_0}{2\sqrt{2}} (\hat{a}_\odot - \hat{a}_\odot - \hat{a}_\odot^\dagger + \hat{a}_\odot^\dagger) \quad (101) $$

As a consequence, each term $\hat{V}_k$ of the multipolar development can be expressed as a combination of ladder operators. Each term only involves products of k ladder operators $\hat{a}^{(k)} \equiv (\hat{a}_u \times \dots \times \hat{a}_{u'})_k$ factors

$$ \langle n | \hat{V}_{\text{proj}} | 0 \rangle = \sum_k \langle n | \hat{V}_k [\hat{a}^{(k)}] | 0 \rangle = \sum_k \hat{V}_{n0,k} [\hat{a}^{(k)}] \quad (102) $$

The integer k will be called the order of the transition in the following.

## Dipolar term

The dipolar term in the potential is expressed as

$$ \hat{V}_1 = -\hat{d} \cdot \mathbf{E}_{n0}(\mathbf{R}_0, q_z, \omega) = e (\hat{x}E_{n0,x}(\mathbf{R}_0, q_z, \omega) + \hat{y}E_{n0,y}(\mathbf{R}_0, q_z, \omega)) \quad (103) $$

From now on, the $(\mathbf{R}_0, q_z, \omega)$ dependence of the field quantities and the n0 subscript on the components of the electric field will be implicit for readability purpose. We express it in both basis of the harmonic oscillator states.

## Dipolar potential : Hermite-Gauss

$$ \hat{V}_1 = \frac{ew_0}{2} [(\hat{a}_x + \hat{a}_x^\dagger)E_x + (\hat{a}_y + \hat{a}_y^\dagger)E_y] \quad (104) $$

## Dipolar potential : Laguerre-Gauss

$$ \hat{V}_1 = \frac{ew_0}{2} [(\hat{a}_\odot + \hat{a}_\odot^\dagger)E_\odot + (\hat{a}_\odot - \hat{a}_\odot^\dagger)E_\odot] \quad (105) $$

Where we name

$$ E_\odot = \frac{1}{\sqrt{2}}(E_x + iE_y) \quad (106) $$

and

$$ E_\odot = \frac{1}{\sqrt{2}}(E_x - iE_y) \quad (107) $$

The left and right circular polarization components of the electric field respectively.
## Transition dipole moments

One can define linear and circular transition dipole moment as in [18]

$$ \begin{cases} \hat{d}_x = -\frac{ew_0}{2}(\hat{a}_x + \hat{a}_x^\dagger) \\ \hat{d}_y = -\frac{ew_0}{2}(\hat{a}_y + \hat{a}_y^\dagger) \end{cases} \qquad (108) $$

$$ \begin{cases} \hat{d}_{\phi} = -\frac{ew_0}{2}(\hat{a}_{\phi} + \hat{a}_{\phi}^{\dagger}) \\ \hat{d}_{\phi} = -\frac{ew_0}{2}(\hat{a}_{\phi} + \hat{a}_{\phi}^{\dagger}) \end{cases} \qquad (109) $$

Note that the annihilation of a left quantum of angular momentum induces the same transition as the creation of a right quantum of angular momentum. Conservation of angular momentum is here explicit : the creation of a left quantum of angular momentum in the electron transverse wavefunction induces a right circular field in the target, Resulting in zero total change of angular momentum.

## Quadrupolar term

From [1], expressing the Quadrupolar term of the Hamiltonian

$$ \hat{V}_2 = -[\hat{Q} \cdot \hat{\nabla}] \cdot \mathbf{E} \qquad (110) $$

Where the quadrupolar tensor is expressed as

$$ \hat{Q} = \frac{1}{2} q \hat{\mathbf{r}} \otimes \hat{\mathbf{r}} \qquad (111) $$

With $q = -e$ for electrons

$$ \hat{Q}_{\nu\mu} = -\frac{e}{2} \hat{r}_{\nu} \hat{r}_{\mu} \qquad (112) $$

Which can be made in a 2D representation of the Quadrupolar tensor

$$ \hat{Q} = -\frac{e}{2} \begin{bmatrix} \hat{x}^2 & \hat{x}\hat{y} \\ \hat{y}\hat{x} & \hat{y}^2 \end{bmatrix} \qquad (113) $$

Considering only one electron we get

$$ \hat{V}_2 = \frac{1}{2} e \left[ \hat{x}^2 \partial_x E_x + \hat{y}^2 \partial_y E_y + \hat{x}\hat{y}(\partial_x E_y + \partial_y E_x) \right] \qquad (114) $$

We remind that we consider implicit the $(\mathbf{R}_0, q_z, \omega)$ dependency for the electric field operators. **Nota bene** : The quadrupolar tensor can be found expressed under different ways,

$$ Q_{\nu\mu} = -\frac{e}{2}(r_{\nu}r_{\mu} - \mathbf{r}^2\delta_{\nu\mu}) \qquad (115) $$

Which is not equivalent to (112) since we are working in the xy-plane and $\partial_x E_x + \partial_y E_y = -iq_z E_z$ in vacuum (Maxwell Gauss).
## Quadrupolar potential : Hermite-Gauss

Injecting (100) and (101) in (114), we obtain

$$ \hat{V}_2 = \frac{ew_0^2}{8} [\partial_x E_x(\hat{a}_x + \hat{a}_x^\dagger)^2 + \partial_y E_y(\hat{a}_y + \hat{a}_y^\dagger)^2 + (\partial_x E_y + \partial_y E_x)(\hat{a}_x + \hat{a}_x^\dagger)(\hat{a}_y + \hat{a}_y^\dagger)] \quad (116) $$

Which can be developed into

$$ \hat{V}_2 = \frac{ew_0^2}{8} [\partial_x E_x(\hat{a}_x^2 + \hat{a}_x^{\dagger 2} + \{\hat{a}_x, \hat{a}_x^{\dagger}\}) + \partial_y E_y(\hat{a}_y^2 + \hat{a}_y^{\dagger 2} + \{\hat{a}_y, \hat{a}_y^{\dagger}\}) \\ + (\partial_x E_y + \partial_y E_x)(\hat{a}_x + \hat{a}_x^{\dagger})(\hat{a}_y + \hat{a}_y^{\dagger})] \quad (117) $$

With {$\hat{a}_x, \hat{a}_x^\dagger$} being the anti-commutator. We note that {$\hat{a}_x, \hat{a}_x^\dagger$} = 1 + 2\hat{N}_x

$$ \hat{V}_2 = \frac{ew_0^2}{8} \left[ \partial_x E_x (\hat{a}_x^2 + \hat{a}_x^{\dagger 2} + 2\hat{N}_x + 1) + \partial_y E_y (\hat{a}_y^2 + \hat{a}_y^{\dagger 2} + 2\hat{N}_y + 1) \right. \\ \left. \qquad + (\partial_x E_y + \partial_y E_x)(\hat{a}_x \hat{a}_y + \hat{a}_x \hat{a}_y^\dagger + \hat{a}_x^\dagger \hat{a}_y + \hat{a}_x^\dagger \hat{a}_y^\dagger) \right] \quad (118) $$

## Quadrupolar potential : Laguerre-Gauss

Expressing the quadrupolar term with circular ladder operators by using (83) and (84):

$$ \hat{V}_2 = \frac{ew_0^2}{8} \left[ (\hat{a}_{\odot} + \hat{a}_{\odot}^{\dagger})^2 \partial_{\odot} E_{\odot} + (\hat{a}_{\odot} + \hat{a}_{\odot}^{\dagger})^2 \partial_{\odot} E_{\odot} + \frac{1}{2} \{\hat{a}_{\odot} + \hat{a}_{\odot}^{\dagger}, \hat{a}_{\odot} + \hat{a}_{\odot}^{\dagger}\} (\partial_{\odot} E_{\odot} + \partial_{\odot} E_{\odot}) \right] \quad (119) $$

Or expanded

$$ \hat{V}_2 = \frac{ew_0^2}{8} \left[ (\hat{a}_{\odot}^2 + \hat{a}_{\odot}^{\dagger 2} + 2\hat{a}_{\odot}\hat{a}_{\odot}^{\dagger}) \partial_{\odot} E_{\odot} + (\hat{a}_{\odot}^2 + \hat{a}_{\odot}^{\dagger 2} + 2\hat{a}_{\odot}\hat{a}_{\odot}^{\dagger}) \partial_{\odot} E_{\odot} \right. \\ \left. \qquad + (\hat{a}_{\odot}\hat{a}_{\odot} + \hat{a}_{\odot}^{\dagger}\hat{a}_{\odot}^{\dagger} + 1 + \hat{N}_{\odot} + \hat{N}_{\odot})(\partial_{\odot} E_{\odot} + \partial_{\odot} E_{\odot}) \right] \quad (120) $$

Where we used the Wirtinger derivatives defined here as

$$ \begin{aligned} \partial_{\odot} &= \frac{\partial}{\partial \odot} = \frac{1}{\sqrt{2}} \left( \frac{\partial}{\partial x} + i \frac{\partial}{\partial y} \right) \\ \partial_{\odot} &= \frac{\partial}{\partial \odot} = \frac{1}{\sqrt{2}} \left( \frac{\partial}{\partial x} - i \frac{\partial}{\partial y} \right) \end{aligned} \quad (121) $$

These derivative acting on $E_{\odot}$ and $E_{\odot}$ in their cartesian decomposition give

$$ \begin{aligned} \partial_{\odot} E_{\odot} &= \frac{1}{2}(\partial_x + i\partial_y)(E_x + iE_y) = \frac{1}{2}(\partial_x E_x - \partial_y E_y + i\partial_x E_y + i\partial_y E_x) \\ \partial_{\odot} E_{\odot} &= \frac{1}{2}(\partial_x - i\partial_y)(E_x - iE_y) = \frac{1}{2}(\partial_x E_x - \partial_y E_y - i\partial_x E_y - i\partial_y E_x) \end{aligned} \quad (122) $$

$$ \partial_{\odot} E_{\odot} + \partial_{\odot} E_{\odot} = \frac{1}{2}(\partial_x - i\partial_y)(E_x + iE_y) + \frac{1}{2}(\partial_x + i\partial_y)(E_x - iE_y) = \partial_x E_x + \partial_y E_y \quad (123) $$

## Higher orders

We don't treat explicitly higher orders of the multipolar development in this text. They are quite easily accessible through application of the above formalism. They may reveal interesting properties to be probed. Considering the dipolar moment to be the first order of the development, then for any $k > 1$

$$ \hat{V}_k = -\frac{q}{k!} \hat{\mathbf{r}} \cdot [\hat{\mathbf{r}} \cdot \nabla]^{k-1} \mathbf{E} \quad (124) $$

Development of higher orders under the ladder operator formalism get cumbersome but one could think of using a formal calculation algorithm.
## Shaped EELS interaction computation scheme

We now have the expression for the PSEELS probability (64), the decomposition in two different basis of any transverse state for the electron beam through a combination of ladder operators (77,87) and finally the expression of every order of the projected potential with ladder operators. We have all the tools to compute the transition probability between any initial and final electron transverse state. We here give as an example the computations for orders (0),(1),(2).

## Transition amplitudes and probability

To simplify the notations, we will consider the amplitude $g_{0,i \to j}$ associated to a transition between order i and j beams. The corresponding EELS probabilities are the squared modulus of these amplitudes. To the leading order, one gets:

$$g_{0,i \to j} = \sum_k \mathcal{M}^{(k)}[\hat{V}_k] \delta_{k,|i-j|} \qquad (125)$$

where $\mathcal{M}^{(k)}$ corresponds to the transition of order $k$ involving only the multipolar term $\hat{V}_k$. i and j are the initial and final transverse electron profile orders as defined in section (3.5). $\delta_{k,|i-j|}$ is the Kronecker delta, defined as

$$\delta_{a,b} = \begin{cases} 1 & \text{if } a = b \\ 0 & \text{if } a \neq b \end{cases} \qquad (126)$$

There are non zero contributions to $g_{0,i \to j}$ from higher order terms in the potential $k > |i-j|$. We will show in section (5.2) that they are negligible under this approximation scheme and that (125) an excellent approximation.

We add a superscript to the $\mathcal{M}$ amplitudes to specify which basis is considered (HG or LG) and a subscript to indicate the indices of the initial and final transverse profiles. We choose to remove the $[\hat{V}_k]$ dependence as it is redundant with the $(k)$ superscript indicating the order of the $\mathcal{M}$ term. As an example, $\mathcal{M}_{20,10}^{(1)HG}$ is the amplitude of the transition $|HG_{20}\rangle \to |HG_{10}\rangle$ given by the dipolar term of the potential. The general amplitudes read:

$$\mathcal{M}_{n_x n_y, n'_x n'_y}^{(k)HG} = \frac{1}{\hbar v} \langle HG_{n_x, n'_x} | \hat{V}_k | HG_{n_y, n'_y} \rangle \qquad (127)$$

$$\mathcal{M}_{n_\phi n'_\phi, n'_\phi n_\phi}^{(k)LG} = \frac{1}{\hbar v} \langle LG_{n'_\phi, n_\phi} | \hat{V}_k | LG_{n_\phi, n'_\phi} \rangle \qquad (128)$$

## Computation example

To ease the comprehension, we compute step by step an example: Considering the transition $|HG_{20}\rangle \to |HG_{10}\rangle$. It is a order $|i-j| = |(2+0)-(1+0)| = 1$ process. The transition amplitude $g_{0,2 \to 1}$ is given to leading order by $\mathcal{M}_{20,10}^{(1)HG}$. We compute the amplitude:

$$\mathcal{M}_{20,10}^{(1)HG} = \frac{1}{\hbar v} \langle HG_{1,0} | \hat{V}_1 | HG_{2,0} \rangle \qquad (129)$$

The procedure to compute explicitly this transition amplitude is now direct by using (104), applying the ladder operators on the states (76) to reach the point where no more ladder operators are present and only remain brackets of HG or LG states, being orthogonal sets they are then trivial to compute.

$$\mathcal{M}_{20,10}^{(1)HG} = \frac{ew_0}{2\hbar v} \langle HG_{1,0} | [(\hat{a}_x + \hat{a}_x^\dagger)E_x + (\hat{a}_y + \hat{a}_y^\dagger)E_y] | HG_{2,0} \rangle \qquad (130)$$
And

$$ \langle HG_{1,0} | \hat{a}_y | HG_{2,0} \rangle = \langle HG_{1,1} | HG_{2,0} \rangle = 0 = \langle HG_{1,0} | \hat{a}_y^\dagger | HG_{2,0} \rangle \quad (131) $$

$$ \langle HG_{1,0} | (\hat{a}_x + \hat{a}_x^\dagger) | HG_{2,0} \rangle = \sqrt{2} \langle HG_{1,0} | HG_{1,0} \rangle + \sqrt{3} \langle HG_{1,0} | HG_{3,0} \rangle = \sqrt{2} \quad (132) $$

Leading to

$$ \mathcal{M}_{20,10}^{(1)HG} = \frac{ew_0}{\sqrt{2}\hbar v} E_x \quad (133) $$

And so

$$ \Gamma_{20 \to 10}^{HG}(\mathbf{R}_0, \omega, q_z) = \left| \frac{ew_0}{\sqrt{2}\hbar v} E_{n_0, x}(\mathbf{R}_0, \omega, q_z) \right|^2 \quad (134) $$

Where we brought back all the implicit dependencies. The intensity of the transition from a $|HG_{20}\rangle$ to a $|HG_{10}\rangle$ is proportional to the squared absolute value of the x component of the transition electric field. In fact any transition probability will be proportional to the squared absolute value of a linear combination of electric field components or derivatives. At order $k = 1$, the electron behaves as a small dipole measuring the electric near-field at the nano scale. This shows that the electron beam effectively behaves as a nanoscale source of polarized white light, enabling a prolongation of the optics experiment results to the deep sub-wavelength scale. This property named Optical Polarization Analog is well discussed in recent literature [18, 20].

Following this analytical scheme, the contributions to the different orders of $g_0$ can be computed for any initial and final states in both basis. In the following we list the general results. The selection rules are made explicit by the Kronecker delta found in the resulting amplitudes. The interested reader wanting to use the computation formalism is advised to refer to the present section rather than the general expressions in a first approach.

## Electrostatic term

$$ \hat{V}_0 = \hat{\rho} \cdot \Phi = -e \hat{n} \Phi \quad (135) $$

With $\hat{n}$ the electron number operator, equivalent to identity when considering a single electron as we do. We compute the transition amplitude between any two Hermite-Gauss transverse states

$$ \mathcal{M}_{n_x n_y, n'_x n'_y}^{(0)HG} = \frac{-e}{\hbar v} \langle HG_{n'_x, n'_y} | HG_{n_x, n_y} \rangle \Phi \quad (136) $$

Resulting in

$$ \mathcal{M}_{n_x n_y, n'_x n'_y}^{(0)HG} = \frac{-e}{\hbar v} \delta_{n_x, n'_x} \delta_{n_y, n'_y} \Phi \quad (137) $$

With $\delta_{n_x, n'_x}$ the Kronecker delta. And for Laguerre-Gauss transverse states

$$ \mathcal{M}_{n_\odot n_\odot, n'_\odot n'_\odot}^{(0)LG} = \frac{-e}{\hbar v} \delta_{n_\odot, n'_\odot} \delta_{n_\odot, n'_\odot} \Phi \quad (138) $$

As found before, the 0th order of PSEELS is equivalent to classical EELS experiments. It corresponds to no change in the transverse wavefunction $\Delta n_x = \Delta n_y = \Delta n_\odot = \Delta n'_\odot = 0$.
## Dipolar Hermite transitions

Studying now order $k=1$ transitions carried by the dipolar term of the projected potential.

$$ \hat{V}_1 = \frac{ew_0}{2} [(\hat{a}_x + \hat{a}_x^\dagger)E_x + (\hat{a}_y + \hat{a}_y^\dagger)E_y] \qquad (139) $$

The general amplitude reads:

$$ \mathcal{M}_{n_x n_y, n'_x n'_y}^{(1)HG} = \frac{ew_0}{2\hbar v} [ (\sqrt{n_x}\delta_{n_{x-1}, n'_x}\delta_{n_y, n'_y} + \sqrt{n_x+1}\delta_{n_{x+1}, n'_x}\delta_{n_y, n'_y})E_x \\ + (\sqrt{n_y}\delta_{n_x, n'_x}\delta_{n_{y-1}, n'_y} + \sqrt{n_y+1}\delta_{n_x, n'_x}\delta_{n_{y+1}, n'_y})E_y ] \qquad (140) $$

We get that the only non zero amplitudes, hence the only transitions carrying a linear dipolar moment are those where the transverse wavefunction index changes by one unit exactly: $\Delta n_x = \pm 1$ or $\Delta n_y = \pm 1$. The measured field has the direction of the ladder operator implied in the transition, this is linear momentum conservation.

## Dipolar Laguerre transitions

For Laguerre-Gauss type transverse wavefunctions, the computations are similar.

$$ \hat{V}_1 = \frac{ew_0}{2} [(\hat{a}_\odot + \hat{a}_\odot^\dagger)E_\odot + (\hat{a}_\phi + \hat{a}_\phi^\dagger)E_\phi] \qquad (141) $$

$$ \mathcal{M}_{n_\odot n_\phi, n'_\odot n'_\phi}^{(1)LG} = \frac{ew_0}{2\hbar v} \left[ E_\odot \left( \sqrt{n_\odot} \delta_{n_\odot-1, n'_\odot} \delta_{n_\phi, n'_\phi} + \sqrt{n_\phi+1} \delta_{n_\phi, n'_\phi} \delta_{n_\odot+1, n'_\phi} \right) \right. \\ \left. + E_\phi \left( \sqrt{n_\phi+1} \delta_{n_\phi+1, n'_\phi} \delta_{n_\phi, n'_\phi} + \sqrt{n_\phi} \delta_{n_\phi, n'_\phi} \delta_{n_\phi-1, n'_\phi} \right) \right] \qquad (142) $$

Again, the non zero amplitudes are those for which $\Delta n_\phi = \pm 1$ or $\Delta n_\phi = \pm 1$. Interestingly, for the interaction amplitude to be proportional to $E_\phi$, the left circular component of the electric field, the transverse profile of the electron needs to either lose one quantum of left circular momentum or gain one quantum of right circular momentum : it is the conservation of angular momentum.

## Quadrupolar Hermite transitions

We get to order $k=2$ transition carried by the quadrupolar potential expressed in Hermite-Gauss basis:

$$ \hat{V}_2 = \frac{ew_0^2}{8} \left[ \partial_x E_x (\hat{a}_x^2 + \hat{a}_x^{\dagger 2} + 2\hat{N}_x + 1) + \partial_y E_y (\hat{a}_y^2 + \hat{a}_y^{\dagger 2} + 2\hat{N}_y + 1) \right. \\ \left. + (\partial_x E_y + \partial_y E_x)(\hat{a}_x \hat{a}_y + \hat{a}_x \hat{a}_y^\dagger + \hat{a}_x^\dagger \hat{a}_y + \hat{a}_x^\dagger \hat{a}_y^\dagger) \right] \qquad (143) $$

The optical nearfield quantity probed does not change when probing $|HG_{n_x n_y}\rangle \to |HG_{n'_x n'_y}\rangle$ or $|HG_{n_x n'_y}\rangle \to |HG_{n'_x n_y}\rangle$. The general amplitude reads:

$$ \mathcal{M}_{n_x n_y, n'_x n'_y}^{(2)HG} = \frac{ew_0^2}{8\hbar v} \left[ \begin{aligned} & \partial_x E_x \delta_{n_y, n'_y} \left( \sqrt{n_x(n_x-1)} \delta_{n_{x-2}, n'_x} + \sqrt{(n_x+1)(n_x+2)} \delta_{n_{x+2}, n'_x} + (2n_x+1) \delta_{n_x, n'_x} \right) \\ & \partial_y E_y \delta_{n_x, n'_x} \left( \sqrt{n_y(n_y-1)} \delta_{n_{y-2}, n'_y} + \sqrt{(n_y+1)(n_y+2)} \delta_{n_{y+2}, n'_y} + (2n_y+1) \delta_{n_y, n'_y} \right) \\ & + (\partial_x E_y + \partial_y E_x) \left( \sqrt{n_x n_y} \delta_{n_{x-1}, n'_x} \delta_{n_{y-1}, n'_y} + \sqrt{n_x(n_y+1)} \delta_{n_{x-1}, n'_x} \delta_{n_{y+1}, n'_y} \right. \\ & \left. \qquad + \sqrt{(n_x+1)n_y} \delta_{n_{x+1}, n'_x} \delta_{n_{y-1}, n'_y} + \sqrt{(n_x+1)(n_y+1)} \delta_{n_{x+1}, n'_x} \delta_{n_{y+1}, n'_y} \right) \end{aligned} \right] \qquad (144) $$

This term of interaction can give direct access to nearfield quantities such as $\partial_x E_x$ or $\partial_x E_y + \partial_y E_x$ corresponding to quadrupolar signature in the nearfield. Due to the mixed terms $a_x a_y^\dagger$, we can measure quadrupole systems while keeping first order shaping at the price of complex post-selection ($|HG_{10}\rangle \to |HG_{01}\rangle$).
## Quadrupolar Laguerre transitions

Also order $k=2$ transition but relying on quadrupolar potential expressed in Laguerre-Gauss basis

$$ \hat{V}_2 = \frac{ew_0^2}{8} \left[ (\hat{a}_{\odot}^2 + \hat{a}_{\odot}^{\dagger 2} + 2\hat{a}_{\odot}\hat{a}_{\odot}^{\dagger})\partial_{\odot}E_{\odot} + (\hat{a}_{\odot}^2 + \hat{a}_{\odot}^{\dagger 2} + 2\hat{a}_{\odot}\hat{a}_{\odot}^{\dagger})\partial_{\odot}E_{\odot} \right. \\ \left. + (\hat{a}_{\odot}\hat{a}_{\odot} + \hat{a}_{\odot}^{\dagger}\hat{a}_{\odot}^{\dagger} + 1 + \hat{N}_{\odot} + \hat{N}_{\odot})(\partial_{\odot}E_{\odot} + \partial_{\odot}E_{\odot}) \right] \quad (145) $$

We get to complex expressions when keeping them general

$$ \begin{aligned} \mathcal{M}_{n_{\odot}n_{\odot}, n'_{\odot}n'_{\odot}}^{(2)LG} = \frac{ew_0^2}{8\hbar v} \Biggl[ & \partial_{\odot}E_{\odot} \left( \sqrt{n_{\odot}(n_{\odot}-1)}\delta_{n_{\odot}-2, n'_{\odot}}\delta_{n_{\odot}, n'_{\odot}} + \sqrt{(n_{\odot}+1)(n_{\odot}+2)}\delta_{n_{\odot}, n'_{\odot}}\delta_{n_{\odot}+2, n'_{\odot}} \right. \\ & \left. \qquad + 2\sqrt{n_{\odot}(n_{\odot}+1)}\delta_{n_{\odot}-1, n'_{\odot}}\delta_{n_{\odot}+1, n'_{\odot}} \right) \\ & + \partial_{\odot}E_{\odot} \left( \sqrt{n_{\odot}(n_{\odot}-1)}\delta_{n_{\odot}, n'_{\odot}}\delta_{n_{\odot}-2, n'_{\odot}} + \sqrt{(n_{\odot}+1)(n_{\odot}+2)}\delta_{n_{\odot}+2, n'_{\odot}}\delta_{n_{\odot}, n'_{\odot}} \right. \\ & \left. \qquad + 2\sqrt{(n_{\odot}+1)n_{\odot}}\delta_{n_{\odot}+1, n'_{\odot}}\delta_{n_{\odot}-1, n'_{\odot}} \right) \\ & + (\partial_{\odot}E_{\odot} + \partial_{\odot}E_{\odot}) \left( \sqrt{n_{\odot}n_{\odot}}\delta_{n_{\odot}-1, n'_{\odot}}\delta_{n_{\odot}-1, n'_{\odot}} + \sqrt{(n_{\odot}+1)(n_{\odot}+1)}\delta_{n_{\odot}+1, n'_{\odot}}\delta_{n_{\odot}+1, n'_{\odot}} \right. \\ & \left. \qquad + (1+n_{\odot}+n_{\odot})\delta_{n_{\odot}, n'_{\odot}}\delta_{n_{\odot}, n'_{\odot}} \right) \Biggr] \end{aligned} \quad (146) $$

The Quadrupolar term of the potential is used to access three quantities: $\partial_{\odot}E_{\odot}$; $\partial_{\odot}E_{\odot}$; $\partial_{\odot}E_{\odot} + \partial_{\odot}E_{\odot}$.

## Transition amplitudes to Gaussian final state

As an example, we specify the previous general transition amplitudes for the case of a gaussian final state $|G\rangle$. The transition amplitudes can be reformulated as

$$ \mathcal{M}_{n_x n_y, 00}^{(1)HG} = \frac{ew_0}{2\hbar v} \left[ \delta_{n_x, 1} \delta_{n_y, 0} E_x + \delta_{n_x, 0} \delta_{n_y, 1} E_y \right] \quad (147) $$

$$ \mathcal{M}_{n_{\odot} n_{\odot}, 00}^{(1)LG} = \frac{ew_0}{2\hbar v} \left[ \delta_{n_{\odot}, 1} \delta_{n_{\odot}, 0} E_{\odot} + \delta_{n_{\odot}, 0} \delta_{n_{\odot}, 1} E_{\odot} \right] \quad (148) $$

$$ \mathcal{M}_{n_x n_y, 00}^{(2)HG} = \frac{ew_0^2}{8\hbar v} \left[ \sqrt{2}\partial_x E_x \delta_{n_x, 2} \delta_{n_y, 0} + \sqrt{2}\partial_y E_y \delta_{n_x, 0} \delta_{n_y, 2} + (\partial_x E_y + \partial_y E_x) \delta_{n_x, 1} \delta_{n_y, 1} + (\partial_x E_x + \partial_y E_y) \delta_{n_y, 0} \delta_{n_y, 0} \right] \quad (149) $$

$$ \mathcal{M}_{n_{\odot} n_{\odot}, 00}^{(2)LG} = \frac{ew_0^2}{8\hbar v} \left[ \sqrt{2}\partial_{\odot}E_{\odot}\delta_{n_{\odot}, 0}\delta_{n_{\odot}, 2} + \sqrt{2}\partial_{\odot}E_{\odot}\delta_{n_{\odot}, 2}\delta_{n_{\odot}, 0} + (\partial_{\odot}E_{\odot} + \partial_{\odot}E_{\odot})(\delta_{n_{\odot}, 1}\delta_{n_{\odot}, 1} + \delta_{n_{\odot}, 0}\delta_{n_{\odot}, 0}) \right] \quad (150) $$

## Higher k orders of potential impact on lower $|i-j|$ order transitions amplitudes

Some transitions have contributions from multiple orders of the potential, once the correct $|i-j|=k$ number of creation/anhiliation operator leads to a non-zero transition amplitude, adding a $\hat{a}\hat{a}^{\dagger}$ gives a non zero contribution but corresponds to a higher order $k+2$ of potential implied in the transition. Each order comes with a multiplication by $w_0$ and a derivative of the electric field operator. Since the typical variation length of the electric field is much larger than $w_0$, this contribution negligible in the transition amplitude. For instance the $|G\rangle \rightarrow |G\rangle$ gets a $\mathcal{M}_{00,00}^{(2)HG} = \frac{ew_0^2}{8\hbar v}(\partial_x E_x + \partial_y E_y)$ contribution from the quadrupolar term of the potential. One can check that in the narrow beam limit

$$ |\mathcal{M}_{00,00}^{(2)HG}| \ll |\mathcal{M}_{00,00}^{(0)HG}| \quad (151) $$

This shows that (125) is an excellent approximation and can be considered exact.
## Multi mode EELS

Reminding the PSEELS rate (41), we have

$$ \Gamma^{\text{PSEELS}}(\omega) = \sum_{n}^{\tilde{n}} \left| \frac{1}{\hbar v} \langle \Psi_{f,\perp} | \langle n | \hat{V}_{\text{proj}} | 0 \rangle | \Psi_{i,\perp} \rangle \right|^2 \\ = \sum_{n}^{\tilde{n}} |g_{0,n}|^2 \qquad (152) $$

The measured PSEELS interaction probability contains contributions from the different photonic states n that can be populated through the interaction. This is an 'incoherent' sum, the states are probed independently from each other, the total PSEELS interaction can be computed mode by mode, making our single mode computations before perfectly correct.

## Identifying physical quantities measured under these techniques

Proper pre and post selection of the electron transverse wavefunction can be used to choose precisely which quantity of the Near-field to probe. For instance :

$$ g_{0,11 \to 00}^{HG} \propto \partial_x E_y + \partial_y E_x \qquad (153) $$

$$ g_{0,20 \to 00}^{HG} \propto \partial_x E_x \qquad (154) $$

We draw tables of the probed quantity for each transition for HG beams and for LG beams, presented in the main text.

## Mapping

Retrieving the dependencies, the transition amplitude $g_{0,i \to j}(R_0, q_z, \omega)$ is proportional to nearfield quantities at a given position $R_0$, position of the electron beam. Scanning the electron beam over the sample will give a map of the considered quantity.

## Another point of view : maps and spectra

A map is measured at a precise energy. We can also fix the position of the electron beam and measure a spectrum for a choosen transition. The resulting spectrum gives the symmetry-matching modes excited in the sample as done experimentally in [46]. The PSEELS setup can be seen both as a tool to map certain nearfield quantities or as a tool to precisely know the symmetry of the populated mode in the sample, e.g exciting on demand quadrupolar transitions in the sample.

## References

* hugo.lourenco-martins@cnrs.fr

[1] Lukas Novotny and Bert Hecht. *Principles of Nano-Optics*. Cambridge University Press, Cambridge, reissue edition edition, June 2011.

[2] E M Purcell, H C Torrey, and R V Pound. Resonance Absorption by Nuclear Magnetic Moments in a Solid. *Phys. Rev.*, 69(1-2):37-38, January 1946.

[3] Ulrich Hohenester. *Nano and Quantum Optics: An Introduction to Basic Principles and Theory*. Graduate Texts in Physics. Springer International Publishing, 2020.

[4] A. Dousse, L. Lanco, J. Suffczyński, E. Semenova, A. Miard, A. Lemaître, I. Sagnes, C. Roblin, J. Bloch, and P. Senellart. Controlled Light-Matter Coupling for a Single Quantum Dot Embedded in a Pillar Microcavity Using Far-Field Optical Lithography. *Physical Review Letters*, 101(26):267404, December 2008.
[5] J. Gallego, W. Alt, T. Macha, M. Martinez-Dorantes, D. Pandey, and D. Meschede. Strong Purcell Effect on a Neutral Atom Trapped in an Open Fiber Cavity. *Physical Review Letters*, 121(17):173603, October 2018.

[6] Anna Roslawska, Tomáš Neuman, Benjamin Doppagne, Andrei G. Borisov, Michelangelo Romeo, Fabrice Scheurer, Javier Aizpurua, and Guillaume Schull. Mapping Lamb, Stark, and Purcell Effects at a Chromophore-Picocavity Junction with Hyper-Resolved Fluorescence Microscopy. *Physical Review X*, 12(1):011012, January 2022.

[7] F. S. Hage, G. Radtke, D. M. Kepaptsoglou, M. Lazzeri, and Q. M. Ramasse. Single-atom vibrational spectroscopy in the scanning transmission electron microscope. *Science*, 367(6482):1124-1127, March 2020.

[8] F. J. García de Abajo and M. Kociak. Probing the Photonic Local Density of States with Electron Energy Loss Spectroscopy. *Physical Review Letters*, 100(10):106804, 2008.

[9] a. Asenjo-Garcia and F. J. García De Abajo. Dichroism in the interaction between vortex electron beams, plasmons, and molecules. *Physical Review Letters*, 113(6):1-5, 2014.

[10] Daniel Ugarte and Caterina Ducati. Controlling Multipolar Surface Plasmon Excitation through the Azimuthal Phase Structure of Electron Vortex Beams. 205418:1-9, 2015.

[11] Ruben Van Boxem, Bart Partoens, and Jo Verbeeck. Inelastic electron-vortex-beam scattering. *Physical Review A*, 91(3):032703, March 2015.

[12] Jo Verbeeck, Armand Béché, Knut Müller-Caspary, Giulio Guzzinati, Minh Anh Luong, and Martien Den Hertog. Demonstration of a 2 × 2 programmable phase plate for electrons. *Ultramicroscopy*, 190:58-65, July 2018.

[13] Chu-Ping Yu, Francisco Vega Ibañez, Armand Béché, and Johan Verbeeck. Quantum wavefront shaping with a 48-element programmable phase plate for electrons. *SciPost Physics*, 15(6):223, December 2023.

[14] Vincenzo Grillo, Amir H. Tavabi, Federico Venturi, Hugo Larocque, Roberto Balboni, Gian Carlo Gazzadi, Stefano Frabboni, Peng-Han Lu, Erfan Mafakheri, Frédéric Bouchard, Rafal E. Dunin-Borkowski, Robert W. Boyd, Martin P. J. Lavery, Miles J. Padgett, and Ebrahim Karimi. Measuring the orbital angular momentum spectrum of an electron beam. *Nature Communications*, 8(1):15536, May 2017.

[15] Amir H. Tavabi, Paolo Rosi, Enzo Rotunno, Alberto Roncaglia, Luca Belsito, Stefano Frabboni, Giulio Pozzi, Gian Carlo Gazzadi, Peng-Han Lu, Robert Nijland, Moumita Ghosh, Peter Tiemeijer, Ebrahim Karimi, Rafal E. Dunin-Borkowski, and Vincenzo Grillo. Experimental Demonstration of an Electrostatic Orbital Angular Momentum Sorter for Electron Beams. *Physical Review Letters*, 126(9):094802, March 2021.

[16] K. Y. Bliokh, I. P. Ivanov, G. Guzzinati, L. Clark, R. Van Boxem, A. Béché, R. Juchtmans, M. A. Alonso, P. Schattschneider, F. Nori, and J. Verbeeck. Theory and applications of free-electron vortex states. *Physics Reports*, 690:1-70, May 2017.

[17] Yijie Shen, Xuejiao Wang, Zhenwei Xie, Changjun Min, Xing Fu, Qiang Liu, Mali Gong, and Xiaocong Yuan. Optical vortices 30 years on: OAM manipulation from topological charge to multiple singularities. *Light: Science & Applications*, 8(1):90, October 2019.

[18] Hugo Lourenço-Martins, Davy Gérard, and Mathieu Kociak. Optical polarization analogue in free electron beams. *Nature Physics*, 17(5):598-603, May 2021.

[19] Marc R. Bourgeois, Austin G. Nixon, Matthieu Chalifour, Elliot K. Beutler, and David J. Masiello. Polarization-Resolved Electron Energy Gain Nanospectroscopy With Phase-Structured Electron Beams. *Nano Letters*, page acs.nanolett.2c02375, August 2022.

[20] Marc R. Bourgeois, Austin G. Nixon, Matthieu Chalifour, and David J. Masiello. Optical Polarization Analogs in Inelastic Free Electron Scattering, May 2023.

[21] Austin G. Nixon, Matthieu Chalifour, Marc R. Bourgeois, Michael Sanchez, and David J. Masiello. Inelastic scattering of transversely structured free electrons from nanophotonic targets: Theory and computation. *Physical Review A*, 109(4):043502, April 2024.

[22] Michael Schüler and Jamal Berakdar. Energy-loss spectroscopy of C60 fullerenes with twisted electrons: Influence of orbital-angular-momentum transfer on plasmon generation. *Physical Review A*, 94(5):052710, November 2016.

[23] Matteo Zanfrognini, Enzo Rotunno, Stefano Frabboni, Alicia Sit, Ebrahim Karimi, Ulrich Hohenester, and Vincenzo Grillo. Orbital Angular Momentum and Energy Loss Characterization of Plasmonic Excitations in Metallic Nanostructures in TEM. *ACS Photonics*, 6(3):620-627, March 2019.

[24] Andrea Konečná, Mikołaj K. Schmidt, Rainer Hillenbrand, and Javier Aizpurua. Probing the electromagnetic response of dielectric antennas by vortex electron beams. *Physical Review Research*, 5(2):023192, June 2023.

[25] D. N. Basov, M. M. Fogler, and F. J. García de Abajo. Polaritons in van der Waals materials. *Science*, 354(6309):aag1992, October 2016.

[26] Nicholas Rivera and Ido Kaminer. Light-matter interactions with photonic quasiparticles. *Nature Reviews Physics*, 2(10):538-561, October 2020.

[27] Ofer Kfir, Valerio Di Giulio, F. Javier García de Abajo, and Claus Ropers. Optical coherence transfer mediated by free electrons. *Science Advances*, 7(18):eabf6380, April 2021.

[28] F. J. García de Abajo. Optical excitations in electron microscopy. *Reviews of Modern Physics*, 82(1):209-275, February 2010.

[29] Scott D Findlay, Mark P Oxley, and Leslie J Allen. Modeling Atomic-Resolution Scanning Transmission Electron Microscopy Images. *Microscopy and Microanalysis*, 14(1):48-59, February 2008.

[30] C. Dwyer. Multislice theory of fast electron scattering incorporating atomic inner-shell ionization. *Ultramicroscopy*, 104(2):141-151, September 2005.

[31] W. Coene and D. Van Dyck. Inelastic scattering of high-energy electrons in real space. *Ultramicroscopy*, 33(4):261-267, October 1990.
[32] A. Lubk and J. Rusz. Jacob's ladder of approximations to paraxial dynamic electron scattering. *Physical Review B*, 92(23):235114, December 2015.

[33] Bernard Diu and Franck Laloë, Claude Cohen-Tannoudji. *Cohen-Tannoudji, Diu and Laloë - Quantum Mechanics (Vol. I, II and III, 2nd Ed.).* December 2019.

[34] G. Nienhuis and L. Allen. Paraxial wave optics and harmonic oscillators. *Physical Review A*, 48(1):656-665, July 1993.

[35] G. M. Vanacore, G. Berruto, I. Madan, E. Pomarico, P. Biagioni, R. J. Lamb, D. McGrouther, O. Reinhardt, I. Kaminer, B. Barwick, H. Larocque, V. Grillo, E. Karimi, F. J. García de Abajo, and F. Carbone. Ultrafast generation and control of an electron vortex beam via chiral plasmonic near fields. *Nature Materials*, 18(6):573-579, June 2019.

[36] Ulrich Hohenester and Andreas Trügler. MNPBEM – A Matlab toolbox for the simulation of plasmonic nanoparticles. *Computer Physics Communications*, 183(2):370-381, February 2012.

[37] Ofer Kfir. Entanglements of Electrons and Cavity Photons in the Strong-Coupling Regime. *Physical Review Letters*, 123(10):103602, September 2019.

[38] Andrea Konečná, Fadil Iyikanat, and F. Javier García de Abajo. Entangling free electrons and optical excitations. *Science Advances*, 8(47):eabo7853, November 2022.

[39] Armin Feist, Guanhao Huang, Germaine Arend, Yujia Yang, Jan-Wilke Henke, Arslan Sajid Raja, F. Jasmin Kappert, Rui Ning Wang, Hugo Lourenço-Martins, Zheru Qiu, Junqiu Liu, Ofer Kfir, Tobias J. Kippenberg, and Claus Ropers. Cavity-mediated electron-photon pairs. *Science*, 377(6607):777-780, August 2022.

[40] Germaine Arend, Guanhao Huang, Armin Feist, Yujia Yang, Jan-Wilke Henke, Zheru Qiu, Hao Jeng, Arslan Sajid Raja, Rudolf Haindl, Rui Ning Wang, Tobias J. Kippenberg, and Claus Ropers. Electrons herald non-classical light, September 2024.

[41] G. S. Agarwal and K. Tara. Nonclassical properties of states generated by the excitations on a coherent state. *Physical Review A*, 43(1):492-497, January 1991.

[42] Axel Lubk. Chapter Two - Paraxial Quantum Mechanics. In Peter W. Hawkes, editor, *Advances in Imaging and Electron Physics*, volume 206, pages 15-58. Elsevier, January 2018.

[43] Arthur Losquin and Mathieu Kociak. Link between Cathodoluminescence and Electron Energy Loss Spectroscopy and the Radiative and Full Electromagnetic Local Density of States. *ACS Photonics*, 2(11):1619-1627, November 2015.

[44] Henrik Bruus and Karsten Flensberg. *Many-Body Quantum Theory in Condensed Matter Physics: An Introduction*. Oxford University Press Oxford, September 2004.

[45] Franck Aguilar, Hugo Lourenço-Martins, Damián Montero, Xiaoyan Li, Mathieu Kociak, and Alfredo Campos. Selective Probing of Longitudinal and Transverse Plasmon Modes with Electron Phase-Matching. *The Journal of Physical Chemistry C*, 127(45):22252-22264, November 2023.

[46] Giulio Guzzinati, Armand Béché, Hugo Lourenço-Martins, Jérôme Martin, Mathieu Kociak, and Jo Verbeeck. Probing the symmetry of the potential of localized surface plasmon resonances with phase-shaped electron beams. *Nature Communications*, 8(1):14999, April 2017.

[47] C. Fabre and N. Treps. Modes and states in quantum optics. *Reviews of Modern Physics*, 92(3):035005, September 2020.

[48] E. Karimi, R. W. Boyd, P. de la Hoz, H. de Guise, J. Réháček, Z. Hradil, A. Aiello, G. Leuchs, and L. L. Sánchez-Soto. Radial quantum number of Laguerre-Gauss modes. *Physical Review A*, 89(6):063813, June 2014.

[49] M. W. Beijersbergen, L. Allen, H. E. L. O. van der Veen, and J. P. Woerdman. Astigmatic laser mode converters and transfer of orbital angular momentum. *Optics Communications*, 96(1):123-132, February 1993.

[50] L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman. Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes. *Physical Review A*, 45(11):8185-8189, June 1992.

[51] E. Abramochkin and V. Volostnikov. Beam transformations and nontransformed beams. *Optics Communications*, 83(1):123-135, May 1991.