Probing Quantum Optical Excitations with Fast Electrons

# Probing Quantum Optical Excitations with Fast Electrons

Valerio Di Giulio ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain    Mathieu Kociak Laboratoire de Physique des Solides, CNRS, Université Paris Sud XI, F 91405 Orsay, France    F. Javier García de Abajo ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain ICREA-Institució Catalana de Recerca i Estudis Avançats, Passeig Lluís Companys 23, 08010 Barcelona, Spain
July 2, 2019
###### Abstract

Probing optical excitations with nanometer resolution is important for understanding their dynamics and interactions down to the atomic scale. Electron microscopes currently offer the unparalleled ability of rendering spatially-resolved electron spectra with combined meV and sub-nm resolution, while the use of ultrafast optical pulses enables fs temporal resolution and exposure of the electrons to ultraintense confined optical fields. Here, we theoretically investigate fundamental aspects of the interaction of fast electrons with localized optical modes that are made possible by these advances. We use a quantum-optics description of the optical field to predict that the resulting electron spectra strongly depend on the statistics of the sample excitations (bosonic or fermionic) and their population (Fock, coherent, or thermal), whose autocorrelation functions are directly retrieved from the ratios of electron gain intensities. We further explore feasible experimental scenarios to probe the quantum characteristics of the sampled excitations and their populations.

## I Introduction

Electron energy-loss spectroscopy (EELS) performed in electron microscopes is a fertile source of information on the dielectric properties of materials down to the nanometer scale Howie (2003); Mkhoyan et al. (2007); Krivanek et al. (2014). This technique is widely used to identify chemical species with atomic resolution through their characteristic high-energy core losses Muller et al. (2008); Zhu et al. (2012). Additionally, low-loss EELS provides insight into the spatial and spectral distributions of plasmons in metallic nanostructures García de Abajo (2010); Rossouw and Botton (2013); Kociak and Stephan (2014); Anton Hörl and Hohenester (2015); Guzzinati et al. (2017), and more recently also of phonons in polaritonic materials Krivanek et al. (2014); Lagos et al. (2017) thanks to remarkable advances in instrument resolution. In a parallel effort, the ultrafast dynamics of nanostructured materials and their influence on optical near-fields can be studied by synchronizing the arrivals of fs light and electron pulses at the sample Barwick et al. (2009); Ryabov and Baum (2016); Kozák et al. (2017). Indeed, although photons and electrons interact extremely weakly in free space due to the lack of energy-momentum matching, the evanescent field components produced upon light scattering by material structures breaks the mismatch, giving rise to efficient light-electron interaction, and effectively producing exchanges of multiple quanta between the electron and the optical field, accompanied by a complex sub-fs dynamics Barwick et al. (2009); García de Abajo et al. (2010); Park et al. (2010); Kirchner et al. (2014); Piazza et al. (2015); Feist et al. (2015); Lummen et al. (2016); Echternkamp et al. (2016); Feist et al. (2017); Priebe et al. (2017); Pomarico et al. (2018); Vanacore et al. (2018); Cai et al. (2018); Morimoto and Baum (2018). Based on this principle, photon-induced near-field electron microscopy (PINEM) is performed by analyzing the resulting multiple gain and loss features in the electron spectra.

PINEM experiments have so far relied on coherent light sources (i.e., lasers), for which the measured spectra are well reproduced by assuming sample bosonic excitations that are coherently populated with a large number of quanta. The probability of each electron spectral peak associated with a net exchange of quanta is then simply given by the squared Bessel function , where a single parameter captures the strength of the electron-light interaction, mediated by the optical electric field component along the direction of electron propagation for an electron velocity and photon frequency Vanacore et al. (2018). For nanometer-sized samples (e.g., nm) illuminated at optical frequencies (eV), a field amplitude V/m renders . Eventually, even the zero loss peak (ZLP, corresponding to ) is fully depleted for . The underlying physics is thus described in terms of a classical optical field interacting with the electron through sample-mediated harmonic evanescent fields. However, we expect new physics to arise when departing from this regime by considering anharmonic states of the illuminated sample, such as those associated with fermionic excitations Tizei and Kociak (2013); Meuret et al. (2015); Bourrellier et al. (2016) or when the external light source is not in a coherent state such as that furnished by a laser.

Here, we discuss the interaction of electron beams with individual optical modes and predict nontrivial characteristics of this interaction when the modes are excited through external illumination depending on the mode nature and population statistics. Specifically, the electron spectra resulting form the interaction with bosonic and fermionic excitations exhibit a radically different dependence on the external light intensity. Additionally, the electron spectra for bosonic modes depend dramatically on the photon statistics, giving rise to a varied phenomenology of asymmetric gain and loss peaks at low fluences and distinct intensity dependences under strong pumping. Interestingly, the autocorrelation functions can be directly retrieved from ratios of measured electron gain intensities. We further propose a feasible experimental realization of these ideas based on a sample consisting of an optical cavity that is fed by optically pumped three-level quantum emitters (QEs).

## Ii Interaction between a beam electron and a bosonic excitation

We consider a sample characterized by a single boson mode of frequency interacting with a focused beam electron of momentum and kinetic energy tightly peaked around and , respectively. Assuming the sample to have an extension along the beam direction sufficiently small as to preserve the transversal beam profile in the sample region, we can write the incident electron wave function as , where is a slowly varying function of the moving-frame position . Further adopting the nonrecoil approximation () and neglecting inelastic boson losses (ps lifetimes) during the interaction time (in the fs range), we can write the system Hamiltonian as the sum of unperturbed and interaction parts (see Appendix A for a self-contained derivation from the Dirac equation)

 ^H0 =ℏω0^a†^a+E0−ℏv⋅(i∇+k0), ^H1 =(ev/c)⋅^A,

where and are creation and annihilation operators of the boson mode, is the electron velocity, and describes the minimal coupling between the electron and the electromagnetic vector potential associated with the sample excitations. This formalism is similar to previous works García de Abajo et al. (2010, 2016), with the important difference that we incorporate a free-electron boson term and the vector potential now becomes an operator,

 ^A=(−ic/ω0)[→E0(r)^a−→E∗0(r)^a†],

where is the single-mode electric field amplitude. Upon inspection, taking along , we find the wave function of the sample-electron system to admit the solution

 |ψ(r,t)⟩=ψ0(r,t)∞∑ℓ=−∞∞∑n=0eiω0[ℓ(z/v−t)−nt]fnℓ(z)|n⟩, (1)

where represents the amplitude of the boson Fock state combined with a change in electron energy. Inserting Eq. (1) into the Schrödinger equation , we find that it is indeed a solution, provided the amplitudes satisfy the equation

 dfnℓdz=eℏω0[√nu∗fn−1ℓ+1−√n+1ufn+1ℓ−1], (2)

where . Interestingly, this expression guarantees that is conserved along the interaction (i.e., the number of excitations in the electron-boson system is preserved), which in turn allows us to treat each initial population of the boson state as an independent channel. However, if we are just interested in the transmitted electron spectrum, we can dismiss the relative phases of these channels and initialize the amplitudes as , with the electron prepared in the incident state . After propagation according to Eq. (2), the transmitted EELS probability reduces to , where

 Pℓ=∞∑n=max{ℓ,0}∣∣fnℓ(∞)∣∣2 (3)

is the probability for the electron to change its energy by . In general, we obtain by direct numerical integration of Eq. (2) because analytical solutions are only possible when the mode is initially unpopulated Carruthers and Nieto (1965); García de Abajo (2013) or approximated to always be in coherent states Kfir (2019).

## Iii Weak coupling limit

The first-order perturbation solution of Eq. (2) readily leads to loss and gain probabilities (see Appendix B) and , respectively, where is the average mode population and

 β0=eℏω0∫dzE0ze−iω0z/v (4)

is an electron-boson coupling coefficient (like the PINEM coefficient, but with the electric field normalized to one quantum). We then find that both loss and gain peaks increase in strength with , but their difference is independent of . In fact the electron-boson interaction strength is determined by , which allows us to separate the regimes of weak and strong coupling depending on and , as shown in Fig. 1. Additionally, we observe that the ratio of gains to losses approaches 1 in the limit, which is consistent with the behavior of the weak-coupling ratio .

Continuing the perturbation series, we find that the leading order contribution to the gain peak is directly proportional to the -order autocorrelation function at zero delay Loudon (2000); we find the powerful result (see Appendix B), which allows us to extract the autocorrelation functions from ratios of measured gain intensities.

Incidentally, for a fermionic sample mode (e.g., a two-level system) under external cw illumination, we have , , and , which, to first-order perturbation in the electron-fermion interaction, readily lead to and . As the light intensity increases, the steady-state populations of a lossy fermion approach (see Appendix D), leading to loss and gain peaks of equal strength.

## Iv Interaction with a dipolar excitation

In what follows we focus on a mode described by a transition dipole at the origin, so that the single-mode electric field reduces to , where , , and is the electron impact-parameter. This is an excellent approximation for plasmonic particles and Mie resonators in a spectral region dominated by a dipolar excitation García de Abajo (2010). Inserting this field into Eq. (4), we find the coupling parameter , where the modified Bessel functions are evaluated at , we introduce the relativistic Lorentz factor , and the electron beam is taken to move along and cross the axis. Reassuringly, the resulting lowest-order loss probability for the ground state sample coincides with the integrated EELS probability for a dipolar scatterer (see Appendix E), thus corroborating that the expression assumed for has the correct single-dipole-mode normalization. In what follows, we take a fixed electron energy (10 keV), dipole orientation (), and impact parameter (), and further use as a strength parameter, for which we use feasible values estimated for metallic (plasmonic) or dielectric (Mie resonances) optical cavities (see Appendices F and G).

## V Dependence on boson population statistics

The boson mode can be populated in different ways, and in particular when pumped by external illumination, the photon statistics is directly imprinted on it. Here, we study three representative distributions with average occupation : (1) a Fock state described by , as one would obtain by coupling to quantum emitters (see below); (2) a coherent state characterized by a Poissonian distribution Glauber (1963), as provided by external laser illumination; and (3) a thermal distribution, as produced by illuminating with chaotic light or by heating the sample mode at a temperature commensurate with . The latter is characterized by , corresponding to a Bose-Einstein average occupation , where . Examples of these distributions are plotted in Fig. 2(a-c) for 2, 10, and 50. We further present in Fig. 2(d-o) the evolution of the transmitted electron spectra for each of the distributions as a function of the coupling parameter (vertical scales). The spectra become more asymmetric for smaller because the number of gains cannot substantially exceed (see also Fig. 1), as observed in recent experiments Das et al. (2019), while the number of losses increases indefinitely with .

In contrast, in the limit, the electron spectra become symmetric with respect to [Fig. 2(j-l)] when . We can then approximate the square roots of Eq. (2) by , leading to the same equation as for PINEM Park et al. (2010); Park and Zewail (2012); García de Abajo et al. (2016); Vanacore et al. (2018), whose solution reads (see Appendix C)

 fnℓ(∞)=√pn+ℓeiℓarg{−β0}Jℓ(2√n+ℓ|β0|),

where the prefactor denotes the initial population of the cavity mode before interaction. Inserting this into Eq. (3) and taking the limit, we find the analytical expressions (see Appendix C)

 (5)

depending on the statistics of the mode population, where . For Fock and coherent distributions, this expression coincides with the well-known PINEM probability Park et al. (2010); García de Abajo et al. (2016); Vanacore et al. (2018) for an optical field amplitude ; the resulting spectra [Fig. 2(j,k)] present the predicted García de Abajo et al. (2010) and subsequently measured Feist et al. (2015) quantum-billiard oscillatory structure as a function of both and field strength, as previously studied for model multilevel atoms Eberly et al. (1977). In contrast, a thermal boson distribution leads to a monotonic decrease with increasing and [Fig. 2(l)]. In all cases, we find an average . The limit [Eq. (5)] is nearly reached under the conditions of Fig. 2 for (cf. two rightmost panel columns).

## Vi Interaction with an optical cavity populated through pumped QEs

As a feasible system to explore the above ideas, we consider an optical cavity (e.g., a Mie resonator) hosting a spectrally isolated mode (frequency , inelastic damping rate ) fed by a number of 3-level QEs, as illustrated in Fig. 3(a) (a similar scheme applies to 4-level QEs). The emitters are initialized in their excited state at (e.g., by optical pumping with a pulse), from which they decay to an intermediate state by coupling to the cavity at a rate (same for all QEs for simplicity) with a resonant transition frequency . We further assume fast internal decay from the intermediate to the ground state, so that each QE interacts with the cavity only once. The combined probability for a cavity Fock state with remaining excited emitters follows the equation of motion , which we solve numerically to obtain the time-dependent distribution . Examples of the evolution of the resulting average mode population and second-order autocorrelation are shown in Fig. 3(b-d). The latter starts at at and in the absence of damping evolves toward at long times (see Appendix H), as expected for an assembly of single-photon emitters. For finite cavity damping, reaches a maximum , from which it exhibits an exponential decay, while eventually jumps to large values when becomes very small. For a sufficiently large number of QEs, we find a substantial average population while varies from nearly 2 down to a quantum regimes characterized by . This evolution strongly affects the resulting electron spectra as a function of the time at which the electron-boson interaction occurs after pumping the QEs, shown in Fig. 3(g-i) under the assumption that the interaction time is small compared with both and . The spectra initially resemble those of the thermal distributions of Fig. 2, as expected from the values, and gradually become similar to those of Fock states; for finite cavity damping, they undergo an attenuation similar to Fig. 2(g) as decreases [right part of Figs. 3(h,i)]. To illustrate the feasibility of the assumed parameters, we consider a Mie resonance in a silicon sphere (, radius ) for which we estimate (see Appendix G) a sharp resonance of quality factor for a size parameter , which causes an increase in the decay rate of QEs embedded in it (Purcell effect Purcell (1946)) by an enhancement factor , so for natural QE decay rates  GHz and optical frequencies THz, we have , while the coupling parameter for this mode is with 100-200 keV electrons.

## Vii Conclusion

In summary, we have shown that the interaction of electron beams with near optical fields depends on both the quantum nature of the sample excitations (fermionic vs bosonic) and the statistics of their populations. For bosonic modes, the spectral distribution of losses and gains varies dramatically when comparing Fock, coherent, and thermal distributions. Our simulations reveal that these regimes can be explored by populating an optical cavity through optically pumped quantum emitters. We remark that the autocorrelation functions of the population distributions are directly retrieved from the peak intensities in the electron spectra, thus opening an unexplored avenue for the study of the ultrafast plasmon, hot-electron, and phonon dynamics in optically pumped nanostructures using time-resolved electron microscope spectroscopy.

## Appendix A Interaction between fast electrons and optical modes: A relativistic quantum-optics approach

The interaction of a beam electron with the evanescent field surrounding an illuminated nanostructure has been extensively studied under the assumption of a classical light field García de Abajo et al. (2010); Park et al. (2010); Park and Zewail (2012); García de Abajo et al. (2016); Vanacore et al. (2018). Here, we present a self-contained derivation of the theory needed to describe the interaction with nonclassical evanescent fields. But first, we start by revisiting the interaction with a classical field in a fully relativistic approach Park and Zewail (2012). We represent the optical field by the scalar and vector potentials and , and describe the interaction with the electron through the Dirac equation Messiah (1966); Sakurai (1994)

 [mec2β+c→α⋅(p+ecA)−eφ]Ψ=iℏ∂Ψ∂t, (6)

where is the momentum operator and

 β=[I00−I],→α=[0→σ→σ0]

are Dirac matrices defined in terms of the identity and Pauli matrices and . We now expand the electron 4-component spinor in terms of positive-energy plane waves (i.e., we neglect electron-positron pair creation processes Sakurai (1994)) of well-defined relativistic momentum and energy :

 Ψ=V−1/2∑kψkeik⋅r−iEkt/ℏΨk, (7)

where

 Ψk=[Ak^sBk→σ⋅k^s] (8)

are the plane-wave spinor amplitudes, is a two-component spin unit vector, is the normalization volume, , and . These plane waves are solutions of the noninteracting Dirac equation

 (mec2β+c→α⋅p)Ψk=EkΨk, (9)

so when inserting Eq. (7) into Eq. (6), we can replace by inside the sum.

We focus in this work on incident electrons prepared in a state involving wave vectors tightly packed around a central value , and further adopt the nonrecoil approximation by assuming that the interaction with the optical field effectively produces components also satisfying . Under these conditions, we approximate the energy in Eq. (9) by the first-order Taylor expansion , where . We further notice that inside the sum can be substituted by outside it. Putting these elements together, we can recast Eq. (6) as

 [E0−(ℏ2c2/E0)k0⋅(i∇+k0)+e→α⋅A−eφ]Ψ=iℏ∂Ψ∂t. (10)

We now approximate in Eq. (7) because only a small error is made when replacing by in the plane wave spinor amplitude [Eq. (8)], so we have

 Ψ≈ψ(r,t)Ψk0, (11)

where is a scalar electron wave function. Inserting Eq. (11) into Eq. (10), multiplying both sides of the resulting equation by from the left, and noticing the relations and satisfied by the spinor plane waves, Eq. (10) leads to the scalar equation

 [ E0−(ℏ2c2/E0)k0⋅(i∇+k0) +(ℏec/E0)k0⋅A−eφ]ψ(r,t)=iℏ∂ψ(r,t)∂t,

which we can recast using the electron velocity vector and the ensuing relations and with as

 [E0−ℏv⋅(i∇+k0)+(ev/c)⋅A−eφ]ψ(r,t) =iℏ∂ψ(r,t)∂t. (12)

This expression, which results from the assumption of the nonrecoil approximation, constitutes an effective Schrödinger equation for the scalar amplitude of the electron spinor [see Eq. (11)], but we remark that it fully incorporates relativistic kinematics (i.e., and are the relativistic electron energy and momentum). Incidentally, the nonrecoil approximation allows us to assume that the spinor is conserved during the interaction [see Eq. (11)], and therefore, the electron spin does not change.

We now extend Eq. (12) to account for nonclassical fields by working in a gauge in which and substituting the vector potential by the operator Milonni (1994)

 ^A=∑j(−ic/ωj)[→Ej(r)^aj−→E∗j(r)^a†j],

where the sum is extended over electromagnetic boson modes of creation and annihilation operators and , frequency , and associated electric field . The wave function of the entire electron-field system is thus expanded to describe a distinct scalar electron wave function for each of the possible number states of the boson ensemble, so that we finally write the Schrödinger equation

 (^H0+^H1)|ψ(r,t)⟩=iℏ∂|ψ(r,t)⟩∂t

with

 ^H0 =∑jℏωj^a†j^aj+E0−ℏv⋅(i∇+k0), ^H1 =(ev/c)⋅^A,

where the first term in is introduced to account for the noninteracting optical modes. We adopt these expressions in the main text assuming a single dominant boson mode corresponding to the index value .

## Appendix B Electron-boson interaction in the weak-coupling limit

A perturbative solution can be produced for Eq. (2) of the main text in the weak-coupling limit, provided the variations of all amplitudes are small during electron-boson interaction. When preparing the incident electron in and the boson in state , the nonvanishing elements of the perturbation series satisfy the equations

 dfn+1,1−1dz=√n+1u∗, dfn−1,11dz=−√nu, dfn+2,2−2dz=√n+2u∗fn+1,1−1, dfn,20dz=√nu∗fn−1,11−√n+1ufn+1,1−1, dfn−2,22dz=−√n−1ufn−1,11, …,

where . After interaction, the first-order () amplitudes reduce to and , which upon insertion into Eq. (3) lead to

 P−1=|β0|2(1+¯n), (13a) P1=|β0|2¯n, (13b)

where [see Eq. (4) in the main text] and is the average population corresponding to the boson occupation distribution . These expressions allow us to identify the weak-coupling limit condition as .

In the above series expansion, the lowest-order contribution to the gain peak corresponds to the amplitude , which satisfies the equation

 dfn−ℓ,ℓℓdz=−√n−ℓ+1ufn−ℓ+1,ℓ−1ℓ−1.

By iteratively solving this concatenated series of equations, we find the post-interaction solution

 fn−ℓ,ℓℓ(∞)=(−1)ℓ√n(n−1)…(n−ℓ+1) ×∫∞−∞dz1∫z1−∞dz2⋯∫zℓ−1−∞dzℓu(z1)u(z2)⋯u(zℓ) =(−β0)ℓℓ!,

where the rightmost expression is derived upon examination of the symmetry of the -dimensional integrand upon permutation of its arguments Abrikosov et al. (1965), which allows us to push the upper integration limits to by creating copies of it. The intensity of the gain peak thus becomes , where denotes the average over the mode population. This leads to the powerful result

 PℓPℓ1=g(ℓ)(ℓ!)2,

where is the -order correlation function at zero delay, which can then be directly inferred from a ratio of peak intensities measured in the transmitted electron spectrum.

For coherent states (i.e., a Poissonian distribution), we have (with ) Loudon (2000), leading to gain peak intensity ratios . In contrast, for a thermal distribution one has Loudon (2000), which produces more intense gain peaks with We stress that these results are valid only for weak interactions, as we are assuming that .

Incidentally, if we also assume then is given by analytical expressions for coherent and thermal mode populations [see Eq. (5) in the main text and Sec. C below], for which, using the small argument approximation of both and with Abramowitz and Stegun (1972), we find and , respectively, therefore directly recovering the above results for the ratio, which are valid for gain peaks with arbitrary . Additionally, we note that Fock states lead to the same ratio as coherent states in the limit because they are characterized by (with ).

## Appendix C Electron-boson interaction in the ¯n≫1 limit

In the limit, the bulk of the electron-boson interaction involves high ’s, which we consider to be much larger than the net number of quanta exchanges . This condition is satisfied if the interaction-strength parameter is small (), which is still compatible with a high total effective interaction for sufficiently large . We can then approximate both and by in Eq. (2) of the main text; for each value of , which is conserved during propagation of the electron amplitudes along , the resulting equation coincides with Eq. (3) of Ref. García de Abajo et al., 2016 for the PINEM interaction with an optical field , and therefore, we take from that reference the solution , where is the electron-mode interaction parameter defined by Eq. (4) in the main text and is the population distribution of the mode Fock state before interaction with the electron. Because , we can further approximate and write the probability associated with a net number of quanta exchanges [Eq. (3) in the main text] as

 Pℓ≈∞∑n=0pnJ2ℓ(2√n|β0|). (14)

For a boson prepared in the Fock state , the interaction probabilities reduce to , where and obviously the average mode population is . Likewise, for a coherent state the population distribution approaches a Gaussian in the limit Feller (1968), the width of which () becomes increasingly small compared with the average population as this one increases; we can thus approximate inside the Bessel function of Eq. (14), which leads to , that is, the same result as for the Fock state. We thus conclude that in the large average population limit both Fock and coherent states of the boson mode produce the same types of electron spectra as observed in PINEM experiments.

The situation is however different for a chaotic thermal distribution with average population , where and is the mode temperature. We approach the limit at high temperatures, for which , and consequently, and . Inserting these expressions into Eq. (14) and approximating the sum as an integral with the change of variable , we find

 Pℓ≈∫∞0xdxe−x2J2ℓ(2x|β|)=e−2|β|2Iℓ(2|β|2),

where again , and the rightmost analytical equality (Eq. (6.633-2) of Ref. Gradshteyn and Ryzhik (2007)) allows us to express the result in terms of the modified Bessel function .

## Appendix D Steady-state population of optically pumped bosonic and fermionic modes

### d.1 Fermionic mode

A two-level system (states ,1 of energies ) coupled to a monochromatic light field constitutes a textbook example of light-matter interactions, commonly described through the optical Bloch equations Scully and Zubairy (1997). The density matrix of the system satisfies the equation of motion

 d^ρdt=iℏ[^ρ,^H]+γ2(2^σ^ρ^σ†−^σ†^σ^ρ−^ρ^σ†^σ), (15)

where the Hamiltonian incorporates the interaction with the transition dipole through the coupling energy (we assume to be real and the field aligned with the dipole), and we define and . Additionally, we account for inelastic transition losses at a rate through a Lindbladian term in Eq. (15). Writing the density matrix in the interaction picture as , Eq. (15) reduces to the Bloch equations

 ˙¯n=(−2/ℏ)Im{ρ10ge−iω0t}−κ¯n, ˙ρ10=(−i/ℏ)(1−2¯n)geiω0t−κρ10/2

for the average population and the coherence ; the other two elements of the density matrix are given by and . At resonance (), adopting the rotating-wave approximation (RWA), we have , leading to the steady-state solution ()

 ¯n=1211+Is/I, (16)

which depends on the ratio of the light intensity to the saturation intensity of the system .

The interaction with a fast electron can be described following exactly the same formalism discussed in the main text for a bosonic mode, but replacing the commutating bosonic operator by the anticommutating fermionic operator . This prescription leads to Eq. (2) with vanishing unless or 1. In the weak electron-fermion coupling regime, this results in loss and gain probabilities

 P−1=|β0|2p0=|β0|2(1−¯n), P1=|β0|2p1=|β0|2¯n,

where accounts for the electron-mode coupling strength [see Eq. (4)] and is given by Eq. (16).

### d.2 Bosonic mode

A bosonic mode is described by an expression identical to Eq. (15) with and replaced by the annihilation and creation operators and , which now satisfy the commutation relation . Additionally, must be replaced by a boson occupation index that labels Fock states , while must be substituted by the ladder frequencies . The boson density matrix admits a rigorous analytical solution Carruthers and Nieto (1965); García de Abajo (2013): , where is a coherent state of amplitude (in the interaction picture) satisfying Glauber (1963) and describing a Poisonian distribution with occupation probabilities and average population . By inserting this expression of into the equation of motion, one finds the solution for the mode amplitude driven by and external classical perturbation . In particular, for the monochromatic light field considered above, this integral leads to . Assuming again resonant illumination () and neglecting the off-resonance term , the average population reduces to . The interaction with an electron in the weak coupling regime is then given by Eqs. (13) with the mode occupation written as .

## Appendix E Interaction with a dipolar mode and ensuing EELS probability

We consider a dipolar mode of frequency characterized by a transition electric dipole moment placed at the origin. As an ansatz, we write the single-mode electric field as the one produced by this dipole,

 →E0=[k20p+(p⋅∇)∇]eik0r/r,

where . The interaction parameter can be now calculated upon insertion of this field into Eq. (4) in the main text. Integrating by parts, we find that each of the derivatives in the expression for can be replaced by . Finally, assuming a distance from the electron beam to the dipole and using the integral , where , , and (see Eq. (3.914-4) in Ref. Gradshteyn and Ryzhik (2007), which we use here under the assumption that has an infinitesimal positive imaginary part), we readily find the expression

 β0=(−2eω0/ℏv2γ)[ipxK1(ζ)+(pz/γ)K0(ζ)],

which is used in the main text.

In the weak-coupling regime (see Sec. D), the integral of the EELS probability over the mode spectral peak when the mode is initially depleted () reduces to . For an isotropic particle characterized by a triply-degenerate mode of electric dipoles along the Cartesian directions, the probability is given by the sum over the three polarization directions, which amounts to setting ; this leads to

 Pisotropic−1=|β0|2=|2eω0p/ℏv2γ|2f(ω0b/vγ), (17)

where .

In order to corroborate the correctness of the normalization of in the above ansatz, we compare with the result derived from the classical EELS probability for an isotropic dipolar particle García de Abajo (2010), , where is the polarizability. Linear response theory allows us to write the latter as Pines and Nozières (1966) in terms of the mode dipole and frequency , which upon insertion into the spectral integral reproduces Eq. (17), therefore confirming the ansatz.

## Appendix F Electron-boson coupling strength for plasmonic cavities

Plasmons in metallic nanoparticles constitute excellent candidates to explore the interaction between electrons and optical cavities. Here, we estimate the coupling parameter