Probing Quantum Optical Excitations with Fast Electrons

TL;DR

This paper theoretically investigates how fast electrons interact with localized optical modes, revealing that electron spectra depend on the quantum statistics and populations of sample excitations, enabling quantum characterization at nanometer scales.

Contribution

It introduces a quantum-optics framework to analyze electron spectra interactions with optical modes, linking spectral features to quantum statistics and populations of excitations.

Findings

Electron spectra depend on the quantum statistics of excitations.

Ratios of electron gain intensities reveal autocorrelation functions.

Feasible experiments can probe quantum properties of excitations.

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.

Figures

Equations

$${\hat{\mathcal{H}}}_0$$

$${\hat{\mathcal{H}}}_1$$

$$\hat{\mathbf{A}}=(-\mathrm{i}c/{\omega }_0)\left[{\vec{\mathcal{E}}}_0(\mathbf{r})\hat{a}-{\vec{\mathcal{E}}}_0^{\ast }(\mathbf{r}){\hat{a}}^{†}\right],$$

$$|\psi (\mathbf{r},t)⟩={\psi }_0(\mathbf{r},t)\sum _{\ell =-\infty }^{\infty }\sum _{n=0}^{\infty }{\mathrm{e}}^{\mathrm{i}{\omega }_0[\ell (z/v-t)-nt]}{f}_{\ell }^n(z)|n⟩,$$

$$\frac{d{f}_{\ell }^n}{dz}=\sqrt{n}{u}^{\ast }{f}_{\ell +1}^{n-1}-\sqrt{n+1}u{f}_{\ell -1}^{n+1},$$

$$P_{\ell }=\sum _{n=\max \{0,-\ell \}}^{\infty }{\left|f_{\ell }^n(\infty )\right|}^2$$

$$f_{\ell }^n(\infty )=$$

$$\times \sum _{n^{\prime }}\frac{{\left(-|{\beta }_0{|}^2\right)}^{n^{\prime }}}{n^{\prime }!(\ell +n^{\prime })!(n-n^{\prime })!},$$

$${\beta }_0=\frac{e}{\hslash {\omega }_0}\int dz{\mathcal{E}}_{0z}{\mathrm{e}}^{-\mathrm{i}{\omega }_{0}z/v}$$

$$f_{\ell }^n(\infty )=\sqrt{p_{n+\ell }}{\mathrm{e}}^{\mathrm{i}\ell \arg \{-{\beta }_0\}}{J}_{\ell }\left(2\sqrt{n+\ell }|{\beta }_0|\right),$$

$$\displaystyle P_{\ell}=\left\{\begin{array}[]{ll}J_{\ell}^{2}(2|\beta|),&\quad\quad\text{Fock, coherent}\\ {\rm e}^{-2|\beta|^{2}}I_{\ell}(2|\beta|^{2}),&\quad\quad\text{thermal}\end{array}\right.$$

$$\left[m_{\mathrm{e}}{c}^2\beta +c\vec{\alpha }\cdot (\mathbf{p}+\frac{e}{c}\mathbf{A})-e\phi \right]\Psi =\mathrm{i}\hslash \frac{\partial \Psi }{\partial t},$$

$$\displaystyle\beta=\left[\begin{array}[]{cc}\mathcal{I}&0\\ 0&-\mathcal{I}\end{array}\right],\quad\quad\vec{\alpha}=\left[\begin{array}[]{cc}0&\vec{\sigma}\\ \vec{\sigma}&0\end{array}\right]$$

$$\Psi =V^{-1/2}\sum _{\mathbf{k}}{\psi }_{\mathbf{k}}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{r}-\mathrm{i}{E}_{k}t/\hslash }{\Psi }_{\mathbf{k}},$$

$$\displaystyle\Psi_{\bf k}=\left[\begin{array}[]{c}A_{k}\;\hat{\bf s}\\ B_{k}\;\vec{\sigma}\cdot{\bf k}\;\hat{\bf s}\end{array}\right]$$

$$\left(m_{\mathrm{e}}{c}^2\beta +c\vec{\alpha }\cdot \mathbf{p}\right){\Psi }_{\mathbf{k}}=E_k{\Psi }_{\mathbf{k}},$$

$$\left[E_0-({\hslash }^2{c}^2/E_0){\mathbf{k}}_0\cdot (\mathrm{i}\nabla +{\mathbf{k}}_0)+e\vec{\alpha }\cdot \mathbf{A}-e\phi \right]\Psi =\mathrm{i}\hslash \frac{\partial \Psi }{\partial t}.$$

$$\Psi \approx \psi (\mathbf{r},t){\Psi }_{{\mathbf{k}}_0},$$

$$\displaystyle\big{[}$$

$$\displaystyle+(\hbar ec/E_{0})\,{\bf k}_{0}\cdot{\bf A}-e\varphi\big{]}\psi({\bf r},t)={\rm i}\hbar\frac{\partial\psi({\bf r},t)}{\partial t},$$

$$\left[E_0-\hslash \mathbf{v}\cdot (\mathrm{i}\nabla +{\mathbf{k}}_0)+(e\mathbf{v}/c)\cdot \mathbf{A}-e\phi \right]\psi (\mathbf{r},t)$$

$$=\mathrm{i}\hslash \frac{\partial \psi (\mathbf{r},t)}{\partial t}.$$

$$\hat{\mathbf{A}}=\sum _j(-\mathrm{i}c/{\omega }_j)\left[{\vec{\mathcal{E}}}_j(\mathbf{r}){\hat{a}}_j-{\vec{\mathcal{E}}}_j^{\ast }(\mathbf{r}){\hat{a}}_j^{†}\right],$$

$$({\hat{\mathcal{H}}}_0+{\hat{\mathcal{H}}}_1)|\psi (\mathbf{r},t)⟩=\mathrm{i}\hslash \frac{\partial |\psi (\mathbf{r},t)⟩}{\partial t}$$

$${\hat{\mathcal{H}}}_0$$

$${\hat{\mathcal{H}}}_1$$

$$\mathrm{i}\hslash \frac{d{\alpha }_n}{dt}=\sqrt{n}{g}^{\ast }{\mathrm{e}}^{\mathrm{i}{\omega }_{0}t}{\alpha }_{n-1}+\sqrt{n+1}g{\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}{\alpha }_{n+1},$$

$$|{\psi }_{\mathrm{I}}(t)⟩=\sum _n{\alpha }_n(t)|n⟩.$$

$$|{\psi }_{\mathrm{I}}(t)⟩=\hat{S}(t,t_0)|{\psi }_{\mathrm{I}}(t_0)⟩,$$

$$\hat{S}={\mathrm{e}}^{\mathrm{i}\chi }{\mathrm{e}}^{{\beta }_0^{\ast }{\hat{a}}^{†}-{\beta }_0\hat{a}}$$

$${\beta }_0(t,t_0)=\frac{\mathrm{i}}{\hslash }{\int }_{t_0}^{t}d{t}^{\prime }g(t^{\prime }){\mathrm{e}}^{-\mathrm{i}{\omega }_0{t}^{\prime }}$$

$$\chi (t,t_0)=\frac{-1}{\hslash }{\int }_{t_0}^{t}d{t}^{\prime }\mathrm{Re}\left\{{\beta }_0(t^{\prime },t_0)g^{\ast }(t^{\prime }){\mathrm{e}}^{\mathrm{i}{\omega }_0{t}^{\prime }}\right\}.$$

$$⟨n|\hat{S}|n_0⟩=$$

Full text

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 $\ell$ quanta is then simply given by the squared Bessel function $J_{\ell}^{2}(2|\beta|)$, where a single parameter $\beta=(e/\hbar\omega)\int dz\,\mathcal{E}_{z}(z){\rm e}^{-{\rm i}\omega z/v}$ captures the strength of the electron-light interaction, mediated by the optical electric field component $\mathcal{E}_{z}(z)$ along the direction of electron propagation $z$ for an electron velocity $v$ and photon frequency $\omega$ Vanacore et al. (2018). For nanometer-sized samples (e.g., $\Delta z\sim 100\,$nm) illuminated at optical frequencies ($\hbar\omega\sim 1\,$eV), a field amplitude $\mathcal{E}\sim 10^{7}\,$V/m renders $|\beta|\sim e\Delta z\mathcal{E}/\hbar\omega\sim 1$. Eventually, even the zero loss peak (ZLP, corresponding to $\ell=0$) is fully depleted for $|\beta|\approx 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 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. 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 Kfir (2019).

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 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 $\omega_{0}$ interacting with a focused beam electron of momentum and kinetic energy tightly peaked around $\hbar{\bf 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 $\psi_{0}({\bf r},t)={\rm e}^{{\rm i}({\bf k}_{0}\cdot{\bf r}-E_{0}t/\hbar)}\phi_{0}({\bf r}-{\bf v}t)$, where $\phi_{0}$ is a slowly varying function of the moving-frame position ${\bf r}-{\bf v}t$. Further adopting the nonrecoil approximation ($\hbar\omega_{0}\ll E_{0}$) and neglecting inelastic boson losses ($>$ps lifetimes) during the interaction time (in the fs range), we can write the system Hamiltonian $\hat{\mathcal{H}}=\hat{\mathcal{H}}_{0}+\hat{\mathcal{H}}_{1}$ as the sum of unperturbed and interaction parts (see Appendix A for a self-contained derivation from the Dirac equation)

[TABLE]

where $\hat{a}^{\dagger}$ and $\hat{a}$ are creation and annihilation operators of the boson mode satisfying the commutation relations of a quantum harmonic oscillator, ${\bf v}=(\hbar{\bf k}_{0}/m_{\rm e})/(1+E_{0}/m_{\rm e}c^{2})$ is the electron velocity, and $\hat{\mathcal{H}}_{1}$ describes the minimal coupling between the electron and the electromagnetic vector potential $\hat{\bf A}$ 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,

[TABLE]

where $\vec{\mathcal{E}}_{0}({\bf r})$ is the single-mode electric field amplitude. Upon inspection, taking ${\bf v}$ along $z$, we find the wave function of the sample-electron system to admit the solution

[TABLE]

where $f_{\ell}^{n}$ represents the amplitude of the boson Fock state $|n\rangle$ combined with a change $\ell\hbar\omega_{0}$ in electron energy. Inserting Eq. (1) into the Schrödinger equation $\hat{\mathcal{H}}|\psi\rangle={\rm i}\hbar\partial|\psi\rangle/\partial t$, we find that it is indeed a solution, provided the amplitudes satisfy the equation

[TABLE]

where $u(z)=(e/\hbar\omega_{0})\mathcal{E}_{0z}(z){\rm e}^{-{\rm i}\omega_{0}z/v}$. Interestingly, this expression guarantees that $n+\ell$ 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 $p_{n}$ of the boson state $|n\rangle$ 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_{\ell}^{n}(-\infty)=\delta_{\ell 0}\sqrt{p_{n}}$, with the electron prepared in the incident state $\ell=0$. After propagation according to Eq. (2), the transmitted EELS probability reduces to $\Gamma(\omega)=\sum_{\ell=-\infty}^{\infty}P_{\ell}\,\delta(\omega+\ell\omega_{0})$, where

[TABLE]

is the probability for the electron to change its energy by $\ell\hbar\omega_{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 Glauber (1963); Lucas et al. (1970)

[TABLE]

where the sum is limited to the range ${\rm max}\{0,-\ell\}\leq n^{\prime}\leq n$, $\chi$ is a global phase that is irrelevant in this study, and

[TABLE]

is an electron-boson coupling coefficient (like the PINEM $\beta$ 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 $\beta_{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}=(\bar{n}+1)|\beta_{0}|^{2}$ and $P_{1}=\bar{n}|\beta_{0}|^{2}$, respectively, where $\bar{n}=\sum_{n=1}^{\infty}np_{n}$ is the average mode population and $\beta_{0}$ is given by Eq. (4). We then find that both loss and gain peaks increase in strength with $|\beta_{0}|^{2}$, but their difference is independent of $\bar{n}$. In fact the electron-boson interaction strength is determined by $\bar{n}|\beta_{0}|^{2}$, which allows us to separate the regimes of weak and strong coupling depending on $|\beta_{0}|^{2}$ and $\bar{n}$, as shown in Fig. 1. Additionally, we observe that the ratio of gains to losses approaches 1 in the $\bar{n}\gg 1$ limit, which is consistent with the behavior of the weak-coupling ratio $\bar{n}/(\bar{n}+1)$.

Continuing the perturbation series, we find that the leading order contribution to the $\ell>0$ gain peak is directly proportional to the $\ell^{\rm th}$-order autocorrelation function at zero delay ${\rm g}^{(\ell)}=\sum_{n=\ell}^{\infty}n(n-1)\dots(n-\ell+1)p_{n}/\bar{n}^{2}$ Loudon (2000); we find the powerful result $P_{\ell}/P_{1}^{\ell}={\rm g}^{(\ell)}/(\ell!)^{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 two-level system) under external cw illumination, we have $f_{\ell}^{n>1}=0$, $p_{1}=\bar{n}$, and $p_{0}=1-\bar{n}$, which, to first-order perturbation in the electron-fermion interaction, readily lead to $P_{-1}=(1-\bar{n})|\beta_{0}|^{2}$and $P_{1}=\bar{n}|\beta_{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 ${\bf p}$ at the origin, so that the single-mode electric field reduces to $\vec{\mathcal{E}}_{0}=\left({\bf p}k_{0}^{2}+{\bf p}\cdot\nabla\nabla\right){\rm e}^{{\rm i}k_{0}r}/r$, where $k_{0}=\omega_{0}/c$, $r=\sqrt{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 García de Abajo (2010). Inserting this field into Eq. (4), we find the coupling parameter $\beta_{0}=(-2e\omega_{0}/\hbar v^{2}\gamma)\left[{\rm i}p_{x}K_{1}(\zeta)+(p_{z}/\gamma)K_{0}(\zeta)\right]$, where the modified Bessel functions $K_{m}$ are evaluated at $\zeta=\omega_{0}b/v\gamma$, we introduce the relativistic Lorentz factor $\gamma=1/\sqrt{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 $|\beta_{0}|^{2}$ coincides with the integrated EELS probability for a dipolar scatterer (see Appendix F), thus corroborating that the expression assumed for $\vec{\mathcal{E}}_{0}$ has the correct single-dipole-mode normalization. The values of $\beta_{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 occupation $\bar{n}$: (1) a Fock state described by $p_{n}=\delta_{n,\bar{n}}$, as one would obtain by coupling to quantum emitters (see below); (2) a coherent state characterized by a Poissonian distribution $p_{n}={\rm e}^{-\bar{n}}\bar{n}^{n}/n!$ 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 $T$ commensurate with $\hbar\omega_{0}/{k_{\rm B}}$. The latter is characterized by $p_{n}=\left(1-{\rm e}^{-\theta}\right)\,{\rm e}^{-n\theta}$, corresponding to a Bose-Einstein average occupation $\bar{n}=1/\left({\rm e}^{\theta}-1\right)$, where $\theta=\hbar\omega_{0}/{k_{\rm B}}T$. Examples of these distributions are plotted in Fig. 2(a-c) for $\bar{n}=$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 $\sqrt{\bar{n}}|\beta_{0}|$ (vertical scales). The spectra become more asymmetric for smaller $\bar{n}$ because the number of gains cannot substantially exceed $\bar{n}$ (see also Fig. 1), as observed in recent experiments Das et al. (2019), while the number of losses increases indefinitely with $|\beta_{0}|$.

In contrast, in the $\bar{n}\gg 1$ limit, the electron spectra become symmetric with respect to $\ell=0$ [Fig. 2(g,k,o)] when $|\ell|\ll\bar{n}$. We can then approximate the square roots of Eq. (2) by $\sqrt{n+\ell}$, 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 D)

[TABLE]

where the prefactor $\sqrt{p_{n+\ell}}$ denotes the initial population of the cavity mode before interaction. Inserting this into Eq. (3) and taking the $\bar{n}\gg 1$ limit, we find the analytical expressions (see Appendix D)

[TABLE]

depending on the statistics of the mode population, where $\beta=\sqrt{\bar{n}}\beta_{0}$. 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 $\mathcal{E}_{z}=\sqrt{\bar{n}}\mathcal{E}_{0z}$; the resulting spectra [Fig. 2(g,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 $\ell$ 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 $\ell$ and $|\mathcal{E}_{z}|$ [Fig. 2(n,o)]. In all cases, we find an average $\langle|\ell|\rangle\sim|\beta|$. The $\bar{n}\gg 1$ limit [Eq. (7)] is nearly reached under the conditions of Fig. 2 for $\bar{n}=50$ (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, similar to those already probed by EELS Hyun et al. (2008)) hosting a spectrally isolated mode (frequency $\omega_{0}$, inelastic damping rate $\kappa$) 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 Rejaei and Khavasi (2015). The emitters are initialized in their excited state at $t=0$ (e.g., by optical pumping with a $\pi$ 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 $\omega_{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_{n}^{m}$ for a cavity Fock state $|n\rangle$ with $m$ remaining excited emitters follows the equation of motion $dp_{n}^{m}/dt=g\left[n(m+1)p_{n-1}^{m+1}-(n+1)mp_{n}^{m}\right]+\kappa\left[(n+1)p_{n+1}^{m}-np_{n}^{m}\right]$, which we solve numerically to obtain the time-dependent distribution $p_{n}=\sum_{m=0}^{N}p_{n}^{m}$. Examples of the evolution of the resulting average mode population $\bar{n}$ and second-order autocorrelation ${\rm g}^{(2)}$ are shown in Fig. 3(b,c). The latter starts at ${\rm 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, $\bar{n}$ reaches a maximum $<N$, from which it exhibits an exponential decay, while ${\rm g}^{(2)}$ eventually jumps to large values when $\bar{n}$ becomes very small. For a sufficiently large number of QEs, we find a substantial average population while ${\rm 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/\kappa$. The spectra initially resemble those of the thermal distributions of Fig. 2, as expected from the ${\rm g}^{(2)}\approx 2$ values, and gradually become similar to those of Fock states; for finite cavity damping, they undergo an attenuation similar to Fig. 2(g) as $\sqrt{\bar{n}}|\beta_{0}|$ 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 ($\epsilon=12$, radius $a$) for which we estimate (see Appendix H) a sharp resonance of quality factor $\omega_{0}/\kappa\sim 10^{4}$ for a size parameter $\omega_{0}a/c=2.6775$, which causes an increase in the decay rate of QEs embedded in it (Purcell effect Purcell (1946)) by an enhancement factor $g/g_{0}\sim 100$, so for natural QE decay rates $g_{0}\sim$ GHz and optical frequencies $\omega_{0}\sim 100\,$THz, we have $\kappa/g\sim 0.1$, while the coupling parameter for this mode is $|\beta_{0}|\sim 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.

APPENDIX

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 $\varphi$ and ${\bf A}$, and describe the interaction with the electron through the Dirac equation Messiah (1966); Sakurai (1994)

[TABLE]

where ${\bf p}=-{\rm i}\hbar\nabla$ is the momentum operator and

[TABLE]

are Dirac matrices defined in terms of the $2\times 2$ identity and Pauli matrices $\mathcal{I}$ and $\vec{\sigma}$. 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 $\hbar{\bf k}$ and energy $E_{k}=c\sqrt{m_{\rm e}^{2}c^{2}+\hbar^{2}k^{2}}$:

[TABLE]

where

[TABLE]

are the plane-wave spinor amplitudes, $\hat{\bf s}$ is a two-component spin unit vector, $V$ is the normalization volume, $A_{k}=\sqrt{(E_{k}+m_{\rm e}c^{2})/2E_{k}}$, and $B_{k}=\hbar c/\sqrt{2E_{k}(E_{k}+m_{\rm e}c^{2})}$. These plane waves are solutions of the noninteracting Dirac equation

[TABLE]

so when inserting Eq. (14) into Eq. (8), we can replace $m_{\rm e}c^{2}\beta+c\vec{\alpha}\cdot{\bf p}$ by $E_{k}$ inside the ${\bf k}$ sum.

We focus in this work on incident electrons prepared in a state involving wave vectors ${\bf k}$ tightly packed around a central value ${\bf k}_{0}$, and further adopt the nonrecoil approximation by assuming that the interaction with the optical field effectively produces ${\bf k}$ components also satisfying $|{\bf k}-{\bf k}_{0}|\ll k_{0}$. Under these conditions, we approximate the energy in Eq. (18) by the first-order Taylor expansion $E_{k}\approx E_{0}+(\hbar^{2}c^{2}/E_{0}){\bf k}_{0}\cdot({\bf k}-{\bf k}_{0})$, where $E_{0}=E_{k_{0}}$. We further notice that ${\bf k}$ inside the ${\bf k}$ sum can be substituted by $-{\rm i}\nabla$ outside it. Putting these elements together, we can recast Eq. (8) as

[TABLE]

We now approximate $\Psi_{\bf k}\approx\Psi_{{\bf k}_{0}}$ in Eq. (14) because only a small error is made when replacing ${\bf k}$ by ${\bf k}_{0}$ in the plane wave spinor amplitude [Eq. (17)], so we have

[TABLE]

where $\psi({\bf r},t)=V^{-1/2}\sum_{\bf k}\psi_{\bf k}\;{\rm e}^{{\rm i}{\bf k}\cdot{\bf r}-{\rm i}E_{k}t/\hbar}$ is a scalar electron wave function. Inserting Eq. (20) into Eq. (19), multiplying both sides of the resulting equation by $\Psi_{{\bf k}_{0}}^{\dagger}$ from the left, and noticing the relations $\Psi_{\bf k}^{\dagger}\Psi_{\bf k}=1$ and $\Psi_{\bf k}^{\dagger}\vec{\alpha}\Psi_{\bf k}=(\hbar c/E_{k})\,{\bf k}$ satisfied by the spinor plane waves, Eq. (19) leads to the scalar equation

[TABLE]

which we can recast using the electron velocity vector ${\bf v}=(\hbar{\bf k}_{0}/m_{\rm e})/(1+E_{0}/m_{\rm e}c^{2})$ and the ensuing relations $\hbar{\bf k}_{0}=m_{\rm e}{\bf v}\gamma$ and $E_{0}=m_{\rm e}c^{2}\gamma$ with $\gamma=1/\sqrt{1-v^{2}/c^{2}}$ as

[TABLE]

This expression, which results from the assumption of the nonrecoil approximation, constitutes an effective Schrödinger equation for the scalar amplitude $\psi({\bf r},t)$ of the electron spinor [see Eq. (20)], but we remark that it fully incorporates relativistic kinematics (i.e., $E_{0}$ and $\hbar{\bf k}_{0}$ are the relativistic electron energy and momentum). Incidentally, the nonrecoil approximation allows us to assume that the spinor $\Psi_{{\bf k}_{0}}$ is conserved during the interaction [see Eq. (20)], and therefore, the electron spin $\hat{\bf s}$ does not change.

We now extend Eq. (21) to account for nonclassical fields by working in a gauge in which $\varphi=0$ and substituting the vector potential by the operator Milonni (1994)

[TABLE]

where the sum is extended over electromagnetic boson modes $j$ of creation and annihilation operators $\hat{a}^{\dagger}_{j}$ and $\hat{a}_{j}$, frequency $\omega_{j}$, and associated electric field $\vec{\mathcal{E}}_{j}({\bf r})$. The wave function of the entire electron-field system $|\psi({\bf r},t)\rangle=\sum_{\{n_{j}\}}\psi_{\{n_{j}\}}({\bf r},t)|\{n_{j}\}\rangle$ is thus expanded to describe a distinct scalar electron wave function $\psi_{\{n_{j}\}}({\bf r},t)$ for each of the possible number states $|\{n_{j}\}\rangle$ of the boson ensemble, so that we finally write the Schrödinger equation

[TABLE]

with

[TABLE]

where the first term in $\hat{\mathcal{H}}_{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. (22)-(24) by assuming a single mode $j=0$ and expressing the wave function in terms of amplitudes $f_{\ell}^{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 $\hat{\mathcal{H}}_{0}=\hbar\omega_{0}\hat{a}^{\dagger}\hat{a}$ and $\hat{\mathcal{H}}_{1}=g^{*}(t)\hat{a}^{\dagger}+g(t)\hat{a}$, and write the Schrödinger equation $({\rm i}\hbar\partial/\partial t)|\psi_{\rm S}\rangle=\left(\hat{\mathcal{H}}_{0}+\hat{\mathcal{H}}_{1}\right)|\psi_{\rm S}\rangle$ with the wave function $|\psi_{\rm S}\rangle=\sum_{n}\alpha_{n}{\rm e}^{-{\rm i}n\omega_{0}t}|n\rangle$ expanded in terms of number states $|n\rangle$. The time-dependent amplitudes $\alpha_{n}$ satisfy the equation

[TABLE]

which leads to Eq. (2) by performing the substitutions $t\rightarrow z/v$, $g(t){\rm e}^{-{\rm i}\omega_{0}t}\rightarrow-{\rm i}\hbar v\,u(z)$, and $\alpha_{n}\rightarrow f_{\ell}^{n}$. Considering that $\ell+n$ is conserved during propagation according to Eq. (2), we need to separately apply this equivalence to each $\ell+n$ channel, the value of which is determined by the incident electron state (here assumed to be $\ell=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 $|\psi_{\rm I}\rangle={\rm e}^{{\rm i}\hat{\mathcal{H}}_{0}t/\hbar}|\psi_{\rm S}\rangle$ satisfies $({\rm i}\hbar\partial/\partial t)|\psi_{\rm I}\rangle=\hat{\mathcal{H}}_{1{\rm I}}|\psi_{\rm I}\rangle$ and $\hat{\mathcal{H}}_{1{\rm I}}={\rm e}^{{\rm i}\hat{\mathcal{H}}_{0}t/\hbar}\hat{\mathcal{H}}_{1}{\rm e}^{-{\rm i}\hat{\mathcal{H}}_{0}t/\hbar}=g^{*}(t){\rm e}^{{\rm i}\omega_{0}t}\hat{a}^{\dagger}+g(t){\rm e}^{-{\rm i}\omega_{0}t}\hat{a}$. In terms of number states, we have

[TABLE]

The sought-after analytical solution is then obtained as Carruthers and Nieto (1965); Lucas et al. (1970)

[TABLE]

where

[TABLE]

is the time-evolution operator,

[TABLE]

is an effective coupling coefficient, and we have introduced a generally ignored time-dependent phase

[TABLE]

Incidentally, this solution can be directly verified by taking the time derivative of $\hat{S}$, for which it is useful to apply the relation $[{\rm e}^{\hat{x}},\hat{y}]=C{\rm e}^{\hat{x}}$ together with the Stone-von Neumann theorem ${\rm e}^{\hat{x}+\hat{y}}={\rm e}^{\hat{x}}{\rm e}^{\hat{y}}{\rm e}^{-C/2}$, which are valid when the commutator $[\hat{x},\hat{y}]=C$ is a c-number. The latter theorem also allows us to write $\hat{S}={\rm e}^{{\rm i}\chi}{\rm e}^{\beta_{0}^{*}\hat{a}^{\dagger}}{\rm e}^{-\beta_{0}\hat{a}}{\rm e}^{-|\beta_{0}|^{2}/2}$. By Taylor expanding the exponential operators in this expression, we finally obtain the analytical result Carruthers and Nieto (1965); Lucas et al. (1970)

[TABLE]

where the sum is restricted to the range ${\rm max}\{0,n-n_{0}\}\leq n^{\prime}\leq n$. Equation (26) is then applied in the main text to $n_{0}=\ell+n$ (i.e., the conserved initial sum of $\ell$ and $n$). Finally, we extend to integration region from $t_{0}\rightarrow-\infty$ to $t\rightarrow\infty$ and then Eq. (25) 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_{\ell}^{n}$ are small during electron-boson interaction. When preparing the incident electron in $\ell=0$ and the boson in state $|n\rangle$, the nonvanishing elements of the perturbation series $f_{\ell}^{n}=\sum_{s=0}^{\infty}f_{\ell}^{n,s}$ satisfy the equations

[TABLE]

where $u(z)=(e/\hbar\omega_{0})\mathcal{E}_{0z}(z){\rm e}^{-{\rm i}\omega_{0}z/v}$. After interaction, the first-order ($s=1$) amplitudes reduce to $f_{-1}^{n+1,1}(\infty)=\sqrt{n+1}\,\beta_{0}$ and $f_{1}^{n-1,1}(\infty)=-\sqrt{n}\,\beta_{0}$, which upon insertion into Eq. (3) lead to

[TABLE]

where $\beta_{0}=\int dz\,u(z)$ [see Eq. (4) in the main text] and $\bar{n}=\sum_{n=0}^{\infty}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 $\sqrt{\bar{n}}|\beta_{0}|\ll 1$.

In the above series expansion, the lowest-order contribution to the $\ell>0$ gain peak corresponds to the amplitude $f_{\ell}^{n-\ell,\ell}$, which satisfies the equation

[TABLE]

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

[TABLE]

where the rightmost expression is derived upon examination of the symmetry of the $\ell$-dimensional integrand upon permutation of its $\ell$ arguments Abrikosov et al. (1965), which allows us to push the upper integration limits to $\infty$ by creating $\ell!$ copies of it. The intensity of the $\ell>1$ gain peak thus becomes $P_{\ell}=|\beta_{0}|^{2\ell}\left\langle n(n-1)\cdots(n-\ell+1)\right\rangle/(\ell!)^{2}$, where $\langle\rangle$ denotes the average over the mode population. This leads to the powerful result

[TABLE]

where ${\rm g}^{(\ell)}=\left\langle n(n-1)\cdots(n-\ell+1)\right\rangle/\bar{n}^{\ell}$ is the $\ell^{\rm 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 ${\rm g}^{(\ell)}=1$ (with $\ell>0$) Loudon (2000), leading to gain peak intensity ratios $P_{\ell}/P_{1}^{\ell}=(1/\ell!)^{2}$. In contrast, for a thermal distribution one has ${\rm g}^{(\ell)}=\ell!$ Loudon (2000), which produces more intense gain peaks with $P_{\ell}/P_{1}^{\ell}=1/\ell!$ We stress that these results are valid only for weak interactions, as we are assuming that $|\beta|=\sqrt{\bar{n}}|\beta_{0}|\ll 1$.

Incidentally, if we also assume $\bar{n}\gg 1$ then $P_{\ell}$ 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 $\approx(z/2)^{\ell}/\ell!$ of both $J_{\ell}(z)$ and $I_{\ell}(z)$ with $\ell\geq 0$ Abramowitz and Stegun (1972), we find $P_{\ell}\approx|\beta|^{2\ell}/(\ell!)^{2}$ and $P_{\ell}\approx|\beta|^{2\ell}/\ell!$, respectively, therefore directly recovering the above results for the $P_{\ell}/P_{1}^{\ell}$ ratio. Additionally, we note that Fock states lead to the same ratio as coherent states in the $\bar{n}\gg 1$ limit because they are characterized by ${\rm g}^{(\ell)}=\bar{n}(\bar{n}-1)\dots(\bar{n}-\ell+1)/\bar{n}^{\ell}\approx 1$ (with $\ell>0$).

Appendix D Electron-boson interaction in the $\bar{n}\gg 1$ limit

In the $\bar{n}\gg 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 $\ell$. This condition is satisfied if the interaction-strength parameter is small ($|\beta_{0}|\ll 1$), which is still compatible with a high total effective interaction $\bar{n}|\beta_{0}|^{2}\sim 1$ for sufficiently large $\bar{n}$. We can then approximate both $\sqrt{n}$ and $\sqrt{n+1}$ by $\sqrt{n+\ell}$ in Eq. (2) of the main text; for each value of $n+\ell$, which is conserved during propagation of the electron amplitudes $f_{\ell}^{n}$ along $z$, the resulting equation coincides with Eq. (3) of Ref. García de Abajo et al., 2016 for the PINEM interaction with an optical field $\mathcal{E}_{z}=\mathcal{E}_{0z}\sqrt{n+\ell}$, and therefore, we take from that reference the solution $f_{\ell}^{n}(\infty)=\sqrt{p_{n+\ell}}{\rm e}^{{\rm i}\ell\arg\{-\beta_{0}\}}J_{\ell}(2\sqrt{n+\ell}|\beta_{0}|)$, where $\beta_{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+\ell\rangle$ before interaction with the electron. Because $n\gg|\ell|$, we can further approximate $n+\ell\approx n$ and write the probability associated with a net number $\ell$ of quanta exchanges [Eq. (3) in the main text] as

[TABLE]

For a boson prepared in the Fock state $|n\rangle$, the interaction probabilities reduce to $P_{\ell}=J_{\ell}^{2}(2|\beta|)$, where $\beta=\sqrt{\bar{n}}\beta_{0}$ and obviously the average mode population is $\bar{n}=n$. Likewise, for a coherent state the population distribution approaches a Gaussian $p_{n}\approx{\rm e}^{-(n-\bar{n})^{2}/2\bar{n}}/\sqrt{2\pi\bar{n}}$ in the $\bar{n}\gg 1$ limit Feller (1968), the width of which ($\sim\sqrt{\bar{n}}$) becomes increasingly small compared with the average population $\bar{n}$ as this one increases; we can thus approximate $n\approx\bar{n}$ inside the Bessel function of Eq. (28), which leads to $P_{\ell}\approx J_{\ell}^{2}(2|\beta|)$, 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-{\rm e}^{-\theta})\,{\rm e}^{-n\theta}$ with average population $\bar{n}=1/({\rm e}^{\theta}-1)$, where $\theta=\hbar\omega_{0}/{k_{\rm B}}T$ and $T$ is the mode temperature. We approach the $\bar{n}\gg 1$ limit at high temperatures, for which $\theta\ll 1$, and consequently, $\theta\approx 1/\bar{n}$ and $p_{n}\approx{\rm e}^{-n/\bar{n}}/\bar{n}$. Inserting these expressions into Eq. (28) and approximating the sum as an integral with the change of variable $n/{\bar{n}}=x^{2}$, we find

[TABLE]

where again $\beta=\sqrt{\bar{n}}\beta_{0}$, 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 $I_{\ell}$.

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

E.1 Fermionic mode

A two-level system (states $j=0$,1 of energies $\hbar\varepsilon_{j}$) coupled to a monochromatic light field ${\bf E}(t)={\bf E}_{0}{\rm e}^{-{\rm i}\omega t}+{\bf E}_{0}^{*}{\rm e}^{{\rm i}\omega t}$ 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

[TABLE]

where the Hamiltonian $\hat{\mathcal{H}}=\sum_{j}\hbar\varepsilon_{j}|j\rangle\langle j|+g(t)\left(\hat{\sigma}^{\dagger}+\hat{\sigma}\right)$ 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 $\hat{\sigma}=|0\rangle\langle 1|$ and $\hat{\sigma}^{\dagger}=|1\rangle\langle 0|$. Additionally, we account for inelastic $1\rightarrow 0$ transition losses at a rate $\kappa$ through a Lindbladian term in Eq. (29). Writing the density matrix in the interaction picture as $\hat{\rho}=\sum_{j,j^{\prime}}\rho_{jj^{\prime}}{\rm e}^{{\rm i}(\varepsilon_{j^{\prime}}-\varepsilon_{j})t}|j\rangle\langle j^{\prime}|$, Eq. (29) reduces to the Bloch equations

[TABLE]

for the average population $\bar{n}=\rho_{11}\equiv p_{1}$ and the coherence $\rho_{10}$; the other two elements of the density matrix are given by $\rho_{00}\equiv p_{1}=1-\bar{n}$ and $\rho_{01}=\rho_{10}^{*}$. At resonance ($\omega=\omega_{0}$), adopting the rotating-wave approximation (RWA), we have $g{\rm e}^{\pm{\rm i}\omega_{0}t}\approx-E_{0}p$, leading to the steady-state solution ($\dot{\bar{n}}=\dot{\rho}_{10}=0$)

[TABLE]

which depends on the ratio of the light intensity $I=(c/2\pi)|E_{0}|^{2}$ to the saturation intensity of the system $I_{\rm s}=c(\hbar\kappa)^{2}/16\pi 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 $\hat{a}$ by the anticommutating fermionic operator $\hat{\sigma}$. This prescription leads to Eq. (2) with $f_{\ell}^{n}$ vanishing unless $n=0$ or 1. In the weak electron-fermion coupling regime, this results in loss and gain probabilities

[TABLE]

where $|\beta_{0}|^{2}$ accounts for the electron-mode coupling strength [see Eq. (4)] and $\bar{n}$ is given by Eq. (30).

E.2 Bosonic mode

A bosonic mode is described by an expression identical to Eq. (29) with $\hat{\sigma}$ and $\hat{\sigma}^{\dagger}$ replaced by the annihilation and creation operators $a$ and $a^{\dagger}$, which now satisfy the commutation relation $[\hat{a},\hat{a}^{\dagger}]=1$. Additionally, $j$ must be replaced by a boson occupation index $n=0,1,\dots$ that labels Fock states $|n\rangle$, while $\varepsilon_{j}$ must be substituted by the ladder frequencies $(n+1/2)\omega_{0}$. The boson density matrix admits a rigorous analytical solution Carruthers and Nieto (1965); García de Abajo (2013): $\hat{\rho}=|\xi(t)\rangle\langle\xi(t)|$, where $|\xi\rangle={\rm e}^{-|\xi|^{2}/2}\sum_{n=0}^{\infty}\left(\xi^{n}/\sqrt{n!}\right){\rm e}^{-{\rm i}n\omega_{0}t}|n\rangle$ is a coherent state of amplitude $\xi$ (in the interaction picture) satisfying $a|\xi\rangle={\rm e}^{-{\rm i}\omega_{0}t}\xi|\xi\rangle$ Glauber (1963) and describing a Poisonian distribution with occupation probabilities $p_{n}=\langle n|\hat{\rho}|n\rangle={\rm e}^{-|\xi|^{2}}|\xi|^{2n}/n!$ and average population $\bar{n}=|\xi|^{2}$. By inserting this expression of $\hat{\rho}$ into the equation of motion, one finds the solution $\xi(t)=(-{\rm i}/\hbar)\int_{-\infty}^{t}dt^{\prime}g(t^{\prime}){\rm e}^{{\rm i}\omega_{0}t^{\prime}-(t-t^{\prime})\kappa/2}$ for the mode amplitude driven by and external classical perturbation $g(t)$. In particular, for the monochromatic light field considered above, this integral leads to \xi(t)=(E_{0}p/\hbar)\big{[}{\rm e}^{{\rm i}(\omega_{0}-\omega)t}/(\omega_{0}-\omega-{\rm i}\kappa/2)+{\rm e}^{{\rm i}(\omega_{0}+\omega)t}/(\omega_{0}+\omega-{\rm i}\kappa/2)\big{]}. Assuming again resonant illumination ($\omega=\omega_{0}$) and neglecting the off-resonance term $\propto{\rm e}^{2{\rm i}\omega_{0}t}$, the average population reduces to $\bar{n}=|\xi|^{2}\approx(2E_{0}p/\hbar\kappa)^{2}=I/2I_{\rm s}$. The interaction with an electron in the weak coupling regime is then given by Eqs. (27) with the mode occupation written as $\bar{n}=I/2I_{\rm s}$.

Appendix F Interaction with a dipolar mode and ensuing EELS probability

We consider a dipolar mode of frequency $\omega_{0}$ characterized by a transition electric dipole moment ${\bf p}$ placed at the origin. As an ansatz, we write the single-mode electric field as the one produced by this dipole,

[TABLE]

where $k_{0}=\omega_{0}/c$. The interaction parameter $\beta_{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 $\vec{\mathcal{E}}_{0}$ can be replaced by ${\rm i}\omega_{0}/v$. Finally, assuming a distance $b$ from the electron beam to the dipole and using the integral $\int dz\,{\rm e}^{{\rm i}k_{0}r-{\rm i}\omega_{0}z/v}/r=2K_{0}(\zeta)$, where $r=\sqrt{b^{2}+z^{2}}$, $\zeta=\omega_{0}b/v\gamma$, and $\gamma=1/\sqrt{1-v^{2}/c^{2}}$ (see Eq. (3.914-4) in Ref. Gradshteyn and Ryzhik (2007), which we use here under the assumption that $k_{0}$ has an infinitesimal positive imaginary part), we readily find the expression

[TABLE]

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 ($\bar{n}=0$) reduces to $P_{-1}=|\beta_{0}|^{2}=(2e\omega_{0}/\hbar v^{2}\gamma)^{2}\left[|p_{x}|^{2}K_{1}^{2}(\zeta)+(|p_{z}|/\gamma)^{2}K_{0}^{2}(\zeta)\right]$. 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

[TABLE]

where $f(\zeta)=K_{1}^{2}(\zeta)+K_{0}^{2}(\zeta)/\gamma^{2}$.

In order to corroborate the correctness of the normalization of $\vec{\mathcal{E}}_{0}$ in the above ansatz, we compare $P_{-1}^{\rm isotropic}$ with the result derived from the classical EELS probability for an isotropic dipolar particle García de Abajo (2010), $\Gamma_{\rm EELS,dip}(\omega)=(1/\hbar\pi)(2e\omega/v^{2}\gamma)^{2}f(\omega b/v\gamma)\,{\rm Im}\{\alpha(\omega)\}$, where $\alpha(\omega)$ is the polarizability. Linear response theory allows us to write the latter as Pines and Nozières (1966) $\alpha(\omega)=(|p|^{2}/\hbar)\left[1/(\omega_{0}-\omega-{\rm i}0^{+})+1/(\omega_{0}+\omega+{\rm i}0^{+})\right]$ in terms of the mode dipole $p$ and frequency $\omega_{0}$, which upon insertion into the spectral integral $P_{-1}=\int_{0}^{\infty}d\omega\;\Gamma_{\rm EELS,dip}(\omega)$ reproduces Eq. (31), therefore confirming the ansatz.

Appendix G 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 $|\beta_{0}|$ [see leading factor in Eq. (31)] for two types of metallic nanoparticles in which the aspect ratio allows one to tune their frequency, also affecting the values of $|\beta_{0}|$; in particular, we consider prolate ellipsoids and spherical shells made of silver (permittivity $\epsilon(\omega)\approx\epsilon_{\rm b}-\omega_{\rm p}^{2}/\omega(\omega+{\rm i}/\tau)$ with $\epsilon_{\rm b}\approx 4$, $\hbar\omega_{\rm p}=9.17\,$eV, and $\hbar/\tau=21\,$meV Johnson and Christy (1972)). Following similar methods as those of Ref. Dias and García de Abajo (2019), a prolate ellipsoid of volume $V$ is found to exhibit a normal-to-the-symmetry-axis resonance frequency $\omega_{0}=\omega_{\rm p}/\sqrt{\epsilon_{\rm b}-\epsilon_{0}}$ and an effective transition dipole $p\approx(1-\epsilon_{0})\sqrt{\hbar\omega_{0}V/8\pi(\epsilon_{\rm b}-\epsilon_{0})}$, where $\epsilon_{0}=1-1/L$, $L=(r^{2}/2)\Delta^{-3}\left[\pi/2-\arctan(1/\Delta)-\Delta/r^{2}\right]$ is the depolarization factor for a height-to-diameter aspect ratio $r>1$, and $\Delta=\sqrt{r^{2}-1}$. For a spherical metal shell of small thickness $t$ compared with the radius $a$, filled with a core dielectric $\epsilon_{\rm c}$, and treated in the $t\ll a$ limit, we find $\omega_{0}\approx\left(\omega_{\rm p}/\sqrt{\epsilon_{\rm c}+2}\right)\sqrt{2t/a}$ and $p\approx\sqrt{3\hbar\omega_{0}a^{3}/2(\epsilon_{\rm c}+2)}$. Numerical inspection of these two types particles yields optimum coupling values that can reach $|\beta_{0}|\sim 0.1$ with plasmon energies in the 1 eV region, particle diameters of $\sim 20\,$nm, and electron energies $\sim 10\,$keV when considering either disk-like prolate ellipsoids (aspect ratio $r\sim 5$) or thin shells ($t/a\sim 0.2$) filled with silica ($\epsilon_{\rm 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 $\epsilon=12$ exhibits a narrow resonance at a size parameter $\rho_{0}=\omega_{0}a/c\approx 2.6775$ with a quality factor (frequency divided by width) $Q=\omega_{0}/\kappa\sim 10^{4}$. Analytical EELS calculations based on a previously published formula García de Abajo (2010) predict a peak-integrated excitation probability $\sim 10^{-3}$ (i.e., $|\beta_{0}|\sim 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 Glauber and Lewenstein (1991), we can express the electromagnetic Green tensor as $G({\bf r},{\bf r}^{\prime},\omega)=\sum_{j}{\bf f}_{j}({\bf r})\otimes{\bf f}_{j}^{*}({\bf r}^{\prime})/(\omega^{2}-\omega_{j}^{2}+{\rm i}0^{+})$, where the sum extends over photon modes $j$ in the presence of the cavity, and the mode functions are orthonormalized according to $\int d^{3}{\bf r}\;{\bf f}_{j}({\bf r})\cdot{\bf f}_{j^{\prime}}^{*}({\bf r})\epsilon({\bf r})=\delta_{jj^{\prime}}$; we now argue that the modes contributing to the Mie resonance $\omega_{0}$ should have similar spatial profiles inside the cavity, so we approximate $G({\bf r},{\bf r}^{\prime},\omega)\approx{\bf f}_{0}({\bf r})\otimes{\bf f}_{0}^{*}({\bf r}^{\prime})/(\omega(\omega+{\rm i}\kappa)-\omega_{0}^{2})$ by phenomenologically introducing the resonance width $\kappa$; the Purcell enhancement factor EF is then proportional to the local-density of optical states (LDOS) normalized to the vacuum value García de Abajo and Kociak (2008), which can be calculated from $G$ as EF$\,\approx(-6\pi c^{3}/\omega){\rm Im}\{\hat{\bf n}\cdot G({\bf r},{\bf r},\omega)\cdot\hat{\bf n}\}$, with ${\bf r}={\bf r}^{\prime}$ corresponding to the position of the QE; for resonant coupling $\omega=\omega_{0}$, arguing from the normalization condition that $|{\bf f}_{0}|^{2}\sim 1/\epsilon V$, where $V=4\pi a^{3}/3$ is the cavity volume, we find EF$\,\sim 9Q/2\epsilon\rho_{0}^{3}\sim 250$ for the cavity under consideration. We now envision QEs with a natural decay rate $g_{0}\sim\,$GHz, whose coupling rate increases to $g={\rm EF}\times g_{0}\sim 10^{2}$GHz; for an optical frequency in the $\omega_{0}\sim 100\,$THz range, the cavity damping rate is $\kappa\sim\omega_{0}/Q\sim 10$GHz, thus leading to small values of $\kappa/g\sim 0.1$ similar to those assumed in the main text.

Appendix I Population of a cavity coupled to an ensemble of quantum emitters

We consider an optical cavity hosting a bosonic excitation mode of frequency $\omega_{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 $\omega_{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 $\omega_{0}$. We further incorporate an inelastic decay of the cavity mode at a rate $\kappa$, 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_{n}^{m}(t)$ for the cavity to be in the Fock state $|n\rangle$ while $m$ emitters are still in their excited states. The equation of motion for these probabilities can be readily written as

[TABLE]

with the initial condition $p_{n}^{m}(0)=\delta_{n0}\delta_{mN}$. The $p_{n-1}^{m+1}$ term in the right-hand side of Eq. (32) contributes to increase the $p_{n}^{m}$ 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_{n}^{m}$, as described by the negative term proportional to $g$ in Eq. (32); similar reasoning leads to the remaining two terms proportional to $\kappa$ to account for inelastic cavity losses. Upon inspection, we find that the total probability $\sum_{nm}p_{n}^{m}=1$ is conserved. Additionally, the total number of excitations in the system cannot exceed $N$, so $p_{n}^{m}=0$ if $n\geq N$ or $m\geq N$.

We are interested in the cavity occupation distribution $p_{n}=\sum_{m}p_{n}^{m}$, the average population $\bar{n}=\sum_{n}n\,p_{n}$, and how they influence the coupling to a passing electron. As a rough way of characterizing the population distribution, we consider the instantaneous second-order correlation ${\rm g}^{(2)}=(\bar{n^{2}}-\bar{n})/\bar{n}^{2}$, where $\bar{n^{2}}=\sum_{n}n^{2}p_{n}$. At small times right after pumping the QEs to their excited states, ${\rm g}^{(2)}$ can be obtained analytically by Taylor expanding $p_{n}^{m}$; we find $\bar{n}=Ngt+O(t^{2})$ and $\bar{n^{2}}-\bar{n}=2N(N-1)(gt)^{2}+O(t^{3})$, from which ${\rm g}^{(2)}=2(1-1/N)$ at $t=0^{+}$. By numerically solving Eq. (32), we observe that ${\rm g}^{(2)}$ is generally a decreasing function of $gt$ (see Fig. 2 in the main text). In particular, for $\kappa=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_{n}^{m}(\infty)=\delta_{nN}\delta_{m0}$, which in turn produces ${\rm g}^{(2)}=1-1/N$. Incidentally, a large increase in ${\rm g}^{(2)}$ is eventually observed when $\bar{n}$ drops to very small values in lossy cavities.

We present additional results for cavities with $N=100$ QEs in Fig. 6.

Acknowledgements.

We thank Vahagn Mkhitaryan, Matteo Lostaglio, and Sophie Meuret for helpful and enjoyable discussions. This work has been supported in part by the Spanish MINECO (MAT2017-88492-R and SEV2015- 0522), ERC (Advanced Grant 789104-eNANO), the Catalan CERCA Program, and Fundació Privada Cellex. V.D.G. acknowledges support from the EU through a Marie Skłodowska-Curie grant (COFUND-DP, H2020-MSCA-COFUND-2014, GA n 665884).