Probing Quantum Optical Excitations with Fast Electrons

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 [1][2][3]. This technique is widely used to identify chemical species with atomic resolution through their characteristic high-energy core losses [4,5]. Additionally, low-loss EELS provides insight into the spatial and spectral distributions of plasmons in metallic nanostructures [6][7][8][9][10], and more recently also of phonons in polaritonic materials [3,11] 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 [12][13][14]. 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 [12,[15][16][17][18][19][20][21][22][23][24][25][26][27]. 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 * Corresponding author: javier.garciadeabajo@nanophotonics.es associated with a net exchange of quanta is then simply given by the squared Bessel function J 2 (2|β|), where a single parameter β = (e/ ω) dz E z (z)e −iωz/v captures the strength of the electron-light interaction, mediated by the optical electric field component E z (z) along the direction of electron propagation z for an electron velocity v and photon frequency ω [25]. For nanometer-sized samples (e.g., ∆z ∼ 100 nm) illuminated at optical frequencies ( ω ∼ 1 eV), a field amplitude E ∼ 10 7 V/m renders |β| ∼ e∆zE/ ω ∼ 1. Eventually, even the zero loss peak (ZLP, corresponding to = 0) is fully depleted for |β| ≈ 1.2. 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 [28][29][30] or when the external light source is not in a coherent state such as that furnished by a laser. As an interesting avenue in this direction, electron-photon entanglement has been recently predicted to influence the interaction with an electron when the sample is previously excited by a trailing electron [31].
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 from 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 func-tions 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 ω 0 interacting with a focused beam electron of momentum and kinetic energy tightly peaked around k 0 and E 0 , 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 ψ 0 (r, t) = e i(k0·r−E0t/ ) φ 0 (r − vt), where φ 0 is a slowly varying function of the moving-frame position r − vt. Further adopting the nonrecoil approximation ( ω 0 E 0 ) and neglecting inelastic boson losses (>ps lifetimes) during the interaction time (in the fs range), we can write the system HamiltonianĤ =Ĥ 0 +Ĥ 1 as the sum of unperturbed and interaction parts (see Appendix A for a self-contained derivation from the Dirac equation) whereâ † andâ are creation and annihilation operators of the boson mode satisfying the commutation relations of a quantum harmonic oscillator, v = ( k 0 /m e )/(1 + E 0 /m e c 2 ) is the electron velocity, andĤ 1 describes the minimal coupling between the electron and the electromagnetic vector potentialÂ associated with the sample excitations. This formalism is similar to previous works [15,32], with the important difference that we incorporate a free-electron boson term and the vector potential now becomes an operator, where E 0 (r) is the single-mode electric field amplitude. Upon inspection, taking v along z, we find the wave function of the sample-electron system to admit the solution where f n represents the amplitude of the boson Fock state |n combined with a change ω 0 in electron energy. Inserting Eq. (1) into the Schrödinger equationĤ|ψ = i ∂|ψ /∂t, we find that it is indeed a solution, provided the amplitudes satisfy the equation where u(z) = (e/ ω 0 )E 0z (z)e −iω0z/v . Interestingly, this expression guarantees that n+ is conserved along the interaction (i.e., the number of excitations in the electronboson system is preserved), which in turn allows us to treat each initial population p n of the boson state |n 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 f n (−∞) = δ 0 √ p n , with the electron prepared in the incident state = 0. After propagation according to Eq. (2), the transmitted EELS probability reduces to Γ(ω) = is the probability for the electron to change its energy by ω 0 . Noticing that Eq. (2) is formally equivalent to the Schrödinger equation for a classically driven quantum oscillator (see Appendix B), we find the analytical solution [33,34] where the sum is limited to the range max{0, − } ≤ n ≤ n, χ is a global phase that is irrelevant in this study, and is an electron-boson coupling coefficient (like the PINEM β coefficient, but with the electric field normalized to one quantum); this result, which is in excellent agreement with direct numerical integration of Eq. (2), shows that the interaction depends exclusively on the initial mode population p n and the parameter β 0 defined by Eq. (4).

III. WEAK COUPLING LIMIT
The first-order perturbation solution of Eq. (2) readily leads to loss and gain probabilities (see Appendix C) P −1 = (n + 1)|β 0 | 2 and P 1 =n|β 0 | 2 , respectively, wherē n = ∞ n=1 np n is the average mode population and β 0 is given by Eq. (4). We then find that both loss and gain peaks increase in strength with |β 0 | 2 , but their difference is independent ofn. In fact the electron-boson interaction strength is determined byn|β 0 | 2 , which allows us to separate the regimes of weak and strong coupling depending on |β 0 | 2 andn, as shown in Fig. 1. Additionally, we observe that the ratio of gains to losses approaches 1 in then 1 limit, which is consistent with the behavior of the weak-coupling ration/(n + 1).
Continuing the perturbation series, we find that the leading order contribution to the > 0 gain peak is directly proportional to the th -order autocorrelation function at zero delay g ( ) = [35]; we find the powerful result P /P 1 = g ( ) /( !) 2 (see Appendix C), which allows us to extract the autocorrelation functions from ratios of measured gain intensities.
Incidentally, for a fermionic sample mode (e.g., a twolevel system) under external cw illumination, we have f n>1 = 0, p 1 =n, and p 0 = 1 −n, which, to first-order perturbation in the electron-fermion interaction, readily lead to P −1 = (1 −n)|β 0 | 2 and P 1 =n|β 0 | 2 . As the light intensity increases, the steady-state populations of a lossy fermion approach p 0 = p 1 = 1/2 (see Appendix E), leading to loss and gain peaks of equal strength.

IV. INTERACTION WITH A DIPOLAR EXCITATION
For completeness, we discuss a mode described by a transition dipole p at the origin, so that the single-mode electric field reduces to E 0 = pk 2 0 + p · ∇∇ e ik0r /r, where k 0 = ω 0 /c, r = √ b 2 + z 2 , and b 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 [6]. Inserting this field into Eq. (4), we find the coupling pa- where the modified Bessel functions K m are evaluated at ζ = ω 0 b/vγ, we introduce the relativistic Lorentz factor γ = 1/ 1 − v 2 /c 2 , and the electron beam is taken to move along z and cross the x axis. Reassuringly, the resulting lowest-order loss probability for the ground state sample |β 0 | 2 coincides with the integrated EELS probability for a dipolar scatterer (see Appendix F), thus corroborating that the expression assumed for E 0 has the correct single-dipole-mode normalization. The values of β 0 considered in this study are in the range estimated via this model for the dipolar response of metallic (plasmonic) and dielectric (Mie resonances) optical cavities interacting with few to 100s keV electrons (see Appendices G and H).

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 occupationn: (1) a Fock state described by p n = δ n,n , as one would obtain by coupling to quantum emitters (see below); (2) a coherent state characterized by a Poissonian distribution p n = e −nnn /n! [33], 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 T commensurate with ω 0 /k B . The latter is characterized by Examples of these distributions are plotted in Fig. 2(a-c) forn =2, 10, and 50. We further present in Fig. 2(do) the evolution of the transmitted electron spectra for each of the distributions as a function of the coupling parameter √n |β 0 | (vertical scales). The spectra become more asymmetric for smallern because the number of gains cannot substantially exceedn (see also Fig. 1), as observed in recent experiments [36], while the number of losses increases indefinitely with |β 0 |.
In contrast, in then 1 limit, the electron spectra become symmetric with respect to = 0 [ Fig. 2(g,k,o)] when | | n. We can then approximate the square roots of Eq. (2) by √ n + , leading to the same equation as for PINEM [16,25,32,37], whose solution reads (see Appendix D) where the prefactor √ p n+ denotes the initial population of the cavity mode before interaction. Inserting this into Eq. (3) and taking then 1 limit, we find the analytical expressions (see Appendix D) depending on the statistics of the mode population, where β = √n β 0 . For Fock and coherent distributions, this expression coincides with the well-known PINEM probability [16,25,32] for an optical field amplitude E z = √n E 0z ; the resulting spectra [ Fig. 2(g,k)] present the predicted [15] and subsequently measured [19] quantumbilliard oscillatory structure as a function of both and field strength, as previously studied for model multilevel atoms [38]. In contrast, a thermal boson distribution leads to a monotonic decrease with increasing and |E z | [ Fig. 2(n,o)]. In all cases, we find an average | | ∼ |β|.

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, similar to those already probed by EELS [39]) hosting a spectrally isolated mode (frequency ω 0 , inelastic damping rate κ) fed by a number N of 3-level QEs, as illustrated in Fig. 3(a). A similar scheme applies to 4-level systems, which are also extensively used in experimental realizations of QEs coupled to optical cavities [40]. The emitters are initialized in their excited state at t = 0 (e.g., by optical pumping with a π pulse), from which they decay to an intermediate state by coupling to the cavity at a rate g (same for all QEs for simplicity) with a resonant transition frequency ω 0 . 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 p m n for a cavity Fock state |n with m remaining excited emitters follows the equation of motion dp m n /dt = g n(m + 1)p m+1 n−1 − (n + 1)mp m n + κ (n + 1)p m n+1 − np m n , which we solve numerically to obtain the time-dependent distribution p n = N m=0 p m n . Examples of the evolution of the resulting average mode populationn and second-order autocorrelation g (2) are shown in Fig. 3(b,c). The latter starts at g (2) = 2(1 − 1/N ) at t = 0 and in the absence of damping evolves toward 1−1/N at long times (see Appendix I), as expected for an assembly of N single-photon emitters. For finite cavity damping,n reaches a maximum < N , from which it exhibits an exponential decay, while g (2) eventually jumps to large values whenn becomes very small. For a sufficiently large number of QEs, we find a substantial average population while g (2) varies from nearly 2 down to a quantum regimes characterized by < 1. This evolution strongly affects the resulting electron spectra as a function of the time t 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 1/g and 1/κ. The spectra initially resemble those of the thermal distributions of Fig. 2, as expected from the g (2) ≈ 2 values, and gradually become FIG. 3: Interaction with an optical cavity coupled to pumped quantum emitters (QEs). (a) We consider a bosonic optical cavity (e.g., a Mie resonator) sustaining a single mode (frequency ω0, inelastic decay rate κ) and infiltrated with N three-level QEs. Optical pumping prepares the emitters in their upper energy state at time t = 0, from which they decay to an intermediate level by resonant coupling to the cavity mode at a rate g. (b,c) Temporal evolution of the average cavity mode population n (b) and second-order autocorrelation function at zero delay g (2) (c) for N = 2, 10, and 50 with κ = 0 (solid curves) and κ/N g = 0.1 (dashed curves). (d-i) Evolution of the populations pn (d-f) and the electron spectra (g-i) as a function of the delay time for N = 50 and different cavity decay rates: κ/N g = 0 (d,g), 0.01 (e,h) and 0.1 (f,i). We assume an electron-mode coupling |β0| = 0.7. similar to those of Fock states; for finite cavity damping, they undergo an attenuation similar to Fig. 2 To illustrate the feasibility of the assumed parameters, we consider a Mie resonance in a silicon sphere ( = 12, radius a) for which we estimate (see Appendix H) a sharp resonance of quality factor ω 0 /κ ∼ 10 4 for a size parameter ω 0 a/c = 2.6775, which causes an increase in the decay rate of QEs embedded in it (Purcell effect [41]) by an enhancement factor g/g 0 ∼ 100, so for natural QE decay rates g 0 ∼ GHz and optical frequencies ω 0 ∼ 100 THz, we have κ/g ∼ 0.1, while the coupling parameter for this mode is |β 0 | ∼ 0.03 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, for which we have elaborated a model based on realistic optical cavities (e.g., Mie resonators) infiltrated with gain atoms or molecules (e.g., rhodamine B). We predict that the autocorrelation functions of the population distributions are directly retrievable 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. 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 [15,16,25,32,37]. Here, we present a selfcontained 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 [37]. We represent the optical field by the scalar and vector potentials ϕ and A, and describe the interaction with the electron through the Dirac equation [42,43] where p = −i ∇ is the momentum operator and are Dirac matrices defined in terms of the 2 × 2 identity and Pauli matrices I 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 [43]) of well-defined relativistic momentum k and energy E k = c m 2 e c 2 + 2 k 2 : where are the plane-wave spinor amplitudes,ŝ is a twocomponent spin unit vector, V is the normalization volume, A k = (E k + m e c 2 )/2E k , and B k = c/ 2E k (E k + m e c 2 ). These plane waves are solutions of the noninteracting Dirac equation so when inserting Eq. (A2) into Eq. (A1), we can replace m e c 2 β + c α · p by E k inside the k sum.
We focus in this work on incident electrons prepared in a state involving wave vectors k tightly packed around a central value k 0 , and further adopt the nonrecoil approximation by assuming that the interaction with the optical field effectively produces k components also satisfying |k − k 0 | k 0 . Under these conditions, we approximate the energy in Eq. (A4) by the first-order Taylor expansion We further notice that k inside the k sum can be substituted by −i∇ outside it. Putting these elements together, we can recast Eq. (A1) as We now approximate Ψ k ≈ Ψ k0 in Eq. (A2) because only a small error is made when replacing k by k 0 in the plane wave spinor amplitude [Eq. (A3)], so we have where ψ(r, t) = V −1/2 k ψ k e ik·r−iE k t/ is a scalar electron wave function. Inserting Eq. (A6) into Eq. (A5), multiplying both sides of the resulting equation by Ψ † k0 from the left, and noticing the relations Ψ † k Ψ k = 1 and Ψ † k αΨ k = ( c/E k ) k satisfied by the spinor plane waves, Eq. (A5) leads to the scalar equation which we can recast using the electron velocity vector v = ( k 0 /m e )/(1 + E 0 /m e c 2 ) and the ensuing relations k 0 = m e vγ and E 0 = m e c 2 γ with γ = 1/ 1 − v 2 /c 2 as This expression, which results from the assumption of the nonrecoil approximation, constitutes an effective Schrödinger equation for the scalar amplitude ψ(r, t) of the electron spinor [see Eq. (A6)], but we remark that it fully incorporates relativistic kinematics (i.e., E 0 and k 0 are the relativistic electron energy and momentum). Incidentally, the nonrecoil approximation allows us to assume that the spinor Ψ k0 is conserved during the interaction [see Eq. (A6)], and therefore, the electron spinŝ does not change.
We now extend Eq. (A7) to account for nonclassical fields by working in a gauge in which ϕ = 0 and substituting the vector potential by the operator [44] where the sum is extended over electromagnetic boson modes j of creation and annihilation operatorŝ a † j andâ j , frequency ω j , and associated electric field E j (r). The wave function of the entire electron-field system |ψ(r, t) = {nj } ψ {nj } (r, t)|{n j } is thus expanded to describe a distinct scalar electron wave function ψ {nj } (r, t) for each of the possible number states |{n j } of the boson ensemble, so that we finally write the Schrödinger equation where the first term inĤ 0 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 j = 0.

Appendix B: Analytical solution of Eq. (2) in the main text
Equation (2) is obtained from Eqs. (A8)-(A10) by assuming a single mode j = 0 and expressing the wave function in terms of amplitudes f n through Eq. (1). Importantly, Eq. (2) is formally equivalent to the Schrödinger equation for the coupling of a quantum harmonic oscillator to a classical perturbation g(t).
In order to show this in a tutorial way, we review some textbook physics and consider simpler forms of the free and interaction HamiltoniansĤ 0 = ω 0â †â and H 1 = g * (t)â † +g(t)â, and write the Schrödinger equation (i ∂/∂t)|ψ S = Ĥ 0 +Ĥ 1 |ψ S with the wave function |ψ S = n α n e −inω0t |n expanded in terms of number states |n . The time-dependent amplitudes α n satisfy the equation which leads to Eq. (2) by performing the substitutions t → z/v, g(t)e −iω0t → −i v u(z), and α n → f n . Considering that +n is conserved during propagation according to Eq. (2), we need to separately apply this equivalence to each + n channel, the value of which is determined by the incident electron state (here assumed to be = 0) and each number state of the boson mode. In order to obtain an analytical solution, it is convenient to move to the interaction picture, where the wave function |ψ I = e iĤ0t/ |ψ S satisfies (i ∂/∂t)|ψ I = H 1I |ψ I andĤ 1I = e iĤ0t/ Ĥ 1 e −iĤ0t/ = g * (t)e iω0tâ † + g(t)e −iω0tâ . In terms of number states, we have |ψ I (t) = n α n (t)|n .
The sought-after analytical solution is then obtained as [34,45] |ψ is the time-evolution operator, is an effective coupling coefficient, and we have introduced a generally ignored time-dependent phase Incidentally, this solution can be directly verified by taking the time derivative ofŜ, for which it is useful to apply the relation [ex,ŷ] = Cex together with the Stone-von Neumann theorem ex +ŷ = exeŷe −C/2 , which are valid when the commutator [x,ŷ] = C is a c-number. The latter theorem also allows us to writê S = e iχ e β * 0â † e −β0â e −|β0| 2 /2 . By Taylor expanding the exponential operators in this expression, we finally obtain the analytical result [34,45] n|Ŝ|n 0 =e iχ n 0 !n! e −|β0| 2 /2 (−β 0 ) n0−n (B2) where the sum is restricted to the range max{0, n−n 0 } ≤ n ≤ n. Equation (B2) is then applied in the main text to n 0 = + n (i.e., the conserved initial sum of and n). Finally, we extend to integration region from t 0 → −∞ to t → ∞ and then Eq. (B1) reduces to Eq. (4) in the main text.
Appendix C: 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 f n are small during electronboson interaction. When preparing the incident electron in = 0 and the boson in state |n , the nonvanishing elements of the perturbation series f n = ∞ s=0 f n,s satisfy the equations where u(z) = (e/ ω 0 )E 0z (z)e −iω0z/v . After interaction, the first-order (s = 1) amplitudes reduce to f n+1,1 where β 0 = dz u(z) [see Eq. (4) in the main text] and n = ∞ n=0 n p n is the average population corresponding to the boson occupation distribution p n . These expressions allow us to identify the weak-coupling limit condition as √n |β 0 | 1. In the above series expansion, the lowest-order contribution to the > 0 gain peak corresponds to the amplitude f n− , , which satisfies the equation By iteratively solving this concatenated series of equations, we find the post-interaction solution where the rightmost expression is derived upon examination of the symmetry of the -dimensional integrand upon permutation of its arguments [46], which allows us to push the upper integration limits to ∞ by creating ! copies of it. The intensity of the > 1 gain peak thus becomes P = |β 0 | 2 n(n − 1) · · · (n − + 1) /( !) 2 , where denotes the average over the mode population. This leads to the powerful result where g ( ) = n(n − 1) · · · (n − + 1) /n is the th -order correlation function at zero delay, which can then be directly inferred from a ratio of peak intensities measured in the transmitted electron spectrum. We present some illustrative examples in Fig. 4. For coherent states (i.e., a Poissonian distribution), we have g ( ) = 1 (with > 0) [35], leading to gain peak intensity ratios P /P 1 = (1/ !) 2 . In contrast, for a thermal distribution one has g ( ) = ! [35], which produces more intense gain peaks with P /P 1 = 1/ ! We stress that these results are valid only for weak interactions, as we are assuming that |β| = √n |β 0 | 1. Incidentally, if we also assumen 1 then P is given by analytical expressions for coherent and thermal mode populations [see Eq. (f) in the main text and Sec. D below], for which, using the small argument approximation ≈ (z/2) / ! of both J (z) and I (z) with ≥ 0 [47], we find P ≈ |β| 2 /( !) 2 and P ≈ |β| 2 / !, respectively, therefore directly recovering the above results for the P /P 1 ratio. Additionally, we note that Fock states lead to the same ratio as coherent states in then 1 limit because they are characterized by g ( ) =n(n − 1) . . . (n − + 1)/n ≈ 1 (with > 0).
Appendix D: Electron-boson interaction in then 1 limit

In then
1 limit, the bulk of the electron-boson interaction involves high n'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 (|β 0 | 1), which is still compatible with a high total effective interactionn|β 0 | 2 ∼ 1 for sufficiently largen. We can then approximate both √ n and √ n + 1 by √ n + in Eq. (2) of the main text; for each value of n + , which is conserved during propagation of the electron amplitudes f n along z, the resulting equation coincides with Eq. (3) of Ref. 32 for the PINEM interaction with an optical field E z = E 0z √ n + , and therefore, we take from that reference the solution f n (∞) = √ p n+ e i arg{−β0} J (2 √ n + |β 0 |), where β 0 is the electron-mode interaction parameter defined by Eq. (4) in the main text and p n+l is the population distribution of the mode Fock state |n + before interaction with the electron. Because n | |, we can further approximate n + ≈ n and write the probability associated with a net number of quanta exchanges [Eq. (3) in the main text] as For a boson prepared in the Fock state |n , the interaction probabilities reduce to P = J 2 (2|β|), where β = √n β 0 and obviously the average mode population isn = n. Likewise, for a coherent state the population distribution approaches a Gaussian p n ≈ e −(n−n) 2 /2n / √ 2πn in then 1 limit [48], the width of which (∼ √n ) becomes increasingly small compared with the average population n as this one increases; we can thus approximate n ≈n inside the Bessel function of Eq. (D1), which leads to P ≈ J 2 (2|β|), 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 p n = (1 − e −θ ) e −nθ with average population n = 1/(e θ − 1), where θ = ω 0 /k B T and T is the mode temperature. We approach then 1 limit at high temperatures, for which θ 1, and consequently, θ ≈ 1/n and p n ≈ e −n/n /n. Inserting these expressions into Eq.
(D1) and approximating the sum as an integral with the change of variable n/n = x 2 , we find where again β = √n β 0 , and the rightmost analytical equality (Eq. (6.633-2) of Ref. [49]) allows us to express the result in terms of the modified Bessel function I .
Appendix E: Steady-state population of optically pumped bosonic and fermionic modes

Fermionic mode
A two-level system (states j = 0,1 of energies ε j ) coupled to a monochromatic light field E(t) = E 0 e −iωt + E * 0 e iωt constitutes a textbook example of light-matter interactions, commonly described through the optical Bloch equations [50]. The density matrix of the system satisfies the equation of motion where the HamiltonianĤ = j ε j |j j| + g(t) σ † +σ incorporates the interaction with the transition dipole p through the coupling energy g(t) = −E(t)p (we assume E 0 p to be real and the field aligned with the dipole), and we defineσ = |0 1| andσ † = |1 0|. Additionally, we account for inelastic 1 → 0 transition losses at a rate κ through a Lindbladian term in Eq. (E1). Writing the density matrix in the interaction picture aŝ ρ = j,j ρ jj e i(ε j −εj )t |j j |, Eq. (E1) reduces to the Bloch equationṡ n = (−2/ )Im{ρ 10 ge −iω0t } − κn, for the average populationn = ρ 11 ≡ p 1 and the coherence ρ 10 ; the other two elements of the density matrix are given by ρ 00 ≡ p 1 = 1 −n and ρ 01 = ρ * 10 . At resonance (ω = ω 0 ), adopting the rotating-wave approximation (RWA), we have ge ±iω0t ≈ −E 0 p, leading to the steady-state solution (ṅ =ρ 10 = 0) n = 1 2 which depends on the ratio of the light intensity I = (c/2π)|E 0 | 2 to the saturation intensity of the system I s = c( κ) 2 /16πp 2 . 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 f n vanishing unless n = 0 or 1. In the weak electron-fermion coupling regime, this results in loss and gain probabilities where |β 0 | 2 accounts for the electron-mode coupling strength [see Eq. (4)] andn is given by Eq. (E2).

Appendix F: Interaction with a dipolar mode and ensuing EELS probability
We consider a dipolar mode of frequency ω 0 characterized by a transition electric dipole moment p placed at the origin. As an ansatz, we write the single-mode electric field as the one produced by this dipole, E 0 = [k 2 0 p + (p · ∇)∇]e ik0r /r, where k 0 = ω 0 /c. The interaction parameter β 0 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 z derivatives in the expression for E 0 can be replaced by iω 0 /v. Finally, assuming a distance b from the electron beam to the dipole and using the integral dz e ik0r−iω0z/v /r = 2K 0 (ζ), where r = √ b 2 + z 2 , ζ = ω 0 b/vγ, and γ = 1/ 1 − v 2 /c 2 (see Eq. (3.914-4) in Ref. [49], which we use here under the assumption that k 0 has an infinitesimal positive imaginary part), we readily find the expression which is used in the main text.
In the weak-coupling regime (see Sec. E), the integral of the EELS probability over the mode spectral peak when the mode is initially depleted (n = 0) 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 p x = p y = p z = p; this leads to where f (ζ) = K 2 1 (ζ) + K 2 0 (ζ)/γ 2 . In order to corroborate the correctness of the normalization of E 0 in the above ansatz, we compare P isotropic −1 with the result derived from the classical EELS probability for an isotropic dipolar particle [6], where α(ω) is the polarizability. Linear response theory allows us to write the latter as [52]  Plasmons in metallic nanoparticles constitute excellent candidates to explore the interaction between electrons and optical cavities. Here, we estimate the coupling parameter |β 0 | [see leading factor in Eq. (F1)] for two types of metallic nanoparticles in which the aspect ratio allows one to tune their frequency, also affecting the values of |β 0 |; in particular, we consider prolate ellipsoids and spherical shells made of silver (permittivity (ω) ≈ b −ω 2 p /ω(ω+i/τ ) with b ≈ 4, ω p = 9.17 eV, and /τ = 21 meV [53]). Following similar methods as those of Ref. [54], a prolate ellipsoid of volume V is found to exhibit a normal-to-the-symmetry-axis resonance frequency , where 0 = 1 − 1/L, L = (r 2 /2)∆ −3 π/2 − arctan(1/∆) − ∆/r 2 is the depolarization factor for a height-to-diameter aspect ratio r > 1, and ∆ = √ r 2 − 1. For a spherical metal shell of small thickness t compared with the radius a, filled with a core dielectric c , and treated in the t a limit, we find ω 0 ≈ ω p / √ c + 2 2t/a and p ≈ 3 ω 0 a 3 /2( c + 2). Numerical inspection of these two types particles yields optimum coupling values that can reach |β 0 | ∼ 0.1 with plasmon energies in the 1 eV region, particle diameters of ∼ 20 nm, and electron energies ∼ 10 keV when considering either disk-like prolate ellipsoids (aspect ratio r ∼ 5) or thin shells (t/a ∼ 0.2) filled with silica ( c = 2). The values adopted in the main paper are commensurate with these parameters.

Appendix H: Interaction strength of quantum emitters and beam electrons with dielectric optical cavities
High-index dielectric nanostructures can trap light with small radiative leakage. For example, a Si sphere of radius a modeled with a permittivity = 12 exhibits a narrow resonance at a size parameter ρ 0 = ω 0 a/c ≈ 2.6775 with a quality factor (frequency divided by width) Q = ω 0 /κ ∼ 10 4 . Analytical EELS calculations based on a previously published formula [6] predict a peak-integrated excitation probability ∼ 10 −3 (i.e., |β 0 | ∼ 0.03) for grazingly passing electrons of 100-200 keV kinetic energy (see illustrative calculation in Fig.  5). Besides this relatively high probability, the large Q value of the Mie resonance under consideration produces a high Purcell enhancement in quantum emitters (QEs) when they are embedded inside the structure, implying nearly perfect QE-cavity coupling and negligible radiative losses. Indeed, following a quantum optics formalism for dispersionless and lossless dielectrics [55], we can express the electromagnetic Green tensor as G(r, r , ω) = j f j (r) ⊗ f * j (r )/(ω 2 − ω 2 j + i0 + ), where the sum extends over photon modes j in the presence of the cavity, and the mode functions are orthonormalized according to d 3 r f j (r)·f * j (r) (r) = δ jj ; we now argue that the modes contributing to the Mie resonance ω 0 should have similar spatial profiles inside the cavity, so we approximate G(r, r , ω) ≈ f 0 (r) ⊗ f * 0 (r )/(ω(ω + iκ) − ω 2 0 ) by phenomenologically introducing the resonance width κ; the Purcell enhancement factor EF is then proportional to the local-density of optical states (LDOS) normalized to the vacuum value [56], which can be calculated from G as EF ≈ (−6πc 3 /ω)Im{n · G(r, r, ω) ·n}, with r = r corresponding to the position of the QE; for resonant coupling ω = ω 0 , arguing from the normalization condition that |f 0 | 2 ∼ 1/ V , where V = 4πa 3 /3 is the cavity volume, we find EF ∼ 9Q/2 ρ 3 0 ∼ 250 for the cavity under consideration. We now envision QEs with a natural decay rate g 0 ∼ GHz, whose coupling rate increases to g = EF × g 0 ∼ 10 2 GHz; for an optical frequency in the ω 0 ∼ 100 THz range, the cavity damping rate is κ ∼ ω 0 /Q ∼ 10GHz, thus leading to small values of κ/g ∼ 0.1 similar to those assumed in the main text.  We consider an optical cavity hosting a bosonic excitation mode of frequency ω 0 coupled to N three-level QEs that are prepared in their excited state at time t = 0 [see Fig. 3(a) in the main text]. The QEs can decay to their intermediate states by coupling to the cavity at a rate g, which we assume to be the same for all QEs. The transition frequency in this decay is taken to coincide with the cavity mode ω 0 . Once in their intermediate states, the QEs experience internal decay at a fast rate compared with g, so that they no longer interact with the cavity. The present analysis also applies to 4-level atoms such as those used in lasers with their intermediate-states transition matching ω 0 . We further incorporate an inelastic decay of the cavity mode at a rate κ, associated with either intrinsic losses in the materials or radiative emission.
We characterize the QEs-cavity system in terms of the time-dependent probabilities p m n (t) for the cavity to be in the Fock state |n while m emitters are still in their excited states. The equation of motion for these probabilities can be readily written as dp m n dt =g n(m + 1)p m+1 n−1 − (n + 1)mp m n + κ (n + 1)p m n+1 − np m n (I1) with the initial condition p m n (0) = δ n0 δ mN . The p m+1 n−1 term in the right-hand side of Eq. (I1) contributes to increase the p m n probability due to the decay of an emitter [i.e., when going from m + 1 to m, so it must be proportional to the number of excited QEs before the decay, m + 1] that transfers its energy to the cavity [i.e., from n − 1 to n, so it is proportional to n]; this type of process also subtracts probability from p m n , as described by the negative term proportional to g in Eq. (I1); similar reasoning leads to the remaining two terms proportional to κ to account for inelastic cavity losses. Upon inspection, we find that the total probability nm p m n = 1 is conserved. Additionally, the total number of excitations in the system cannot exceed N , so p m n = 0 if n ≥ N or m ≥ N .
We are interested in the cavity occupation distribution p n = m p m n , the average populationn = n n p n , and how they influence the coupling to a passing elec-tron. As a rough way of characterizing the population distribution, we consider the instantaneous second-order correlation g (2) = (n 2 −n)/n 2 , wheren 2 = n n 2 p n . At small times right after pumping the QEs to their excited states, g (2) can be obtained analytically by Taylor expanding p m n ; we findn = N gt + O(t 2 ) andn 2 −n = 2N (N − 1)(gt) 2 + O(t 3 ), from which g (2) = 2(1 − 1/N ) at t = 0 + . By numerically solving Eq. (I1), we observe that g (2) is generally a decreasing function of gt (see Fig. 2 in the main text). In particular, for κ = 0 the number of excitations in the system must be conserved (i.e., n + m = N ), and eventually all of them are transferred to the cavity, leading to the asymptotic solution p m n (∞) = δ nN δ m0 , which in turn produces g (2) = 1−1/N . Incidentally, a large increase in g (2) is eventually observed whenn drops to very small values in lossy cavities.
We present additional results for cavities with N = 100 QEs in Fig. 6.