Electron beams traversing spherical nanoparticles: analytic and numerical treatment

,


I. INTRODUCTION
In recent decades, electron-beam spectroscopy has emerged as a revolutionary tool for the optical characterization of materials.Swift electrons passing in close proximity or through a specimen undergo energy loss owing to energy transfer to the optical modes sustained in the material [1].From localized and propagating surface plasmons in metallic structures [2][3][4][5], to Mie resonances in dielectric resonators [6][7][8] and phonon polaritons in polar crystals [9][10][11], electron-beam spectroscopy has proven quintessential for mapping collective excitations in a broad spectral range that spans from ultraviolet to far-infrared frequencies.
With the diffraction limit ultimately being controlled by the de Broglie wavelength, highly energetic electrons are excellent probes to study the optical properties of truly nanoscale structures, with atomic spatial resolution and sub-meV energy resolution [12,13].In electron energy-loss spectroscopy (EELS), the sample is excited by a high-energy electron beam (30 − 300 keV), and the energy lost to the interaction is measured in a transmission electron microscope (TEM) setup [14,15].EELS allows, thereby, the detection of both radiative and dark modes, including longitudinal bulk plasmons (BPs) [16], breathing modes [17], or antibonding modes in nanoparticle (NP) dimers [18].Optical excitations in thick samples can be imaged in cathodoluminescence (CL) spectroscopy, performed in scanning electron microscopes (SEMs) at intermediate beam energies (1 − 50 keV) [19,20].In CL measurements, the signal collected is the result of far-field photon emission from the sample, originating from the radiative decay of the excited modes.Recent advances in instrumentation have even added temporal resolution in EEL and CL spectra, introducing the field of ultrafast electron microscopy (UEM) [21,22].
Considering the recent progress in electron spectroscopy techniques, robust analytic and computational tools are evidently required to interpret the plethora of experimental data.While first theoretical efforts were performed within the non-retarded approximation for the description of plasmons in thin films [3], the theory was gradually generalized to account for collective excitations in diverse media and geometries [23,24], also considering retardation effects [25,26].Quantum approaches [27,28] and analytic solutions including relativistic effects for simple geometries were later developed [29,30], allowing the combination of high-velocity electron beams with both common and less conventional materials, including dielectric media, polar crystals, graphene and other twodimensional (2D) materials [8,11,31,32].Most theoretical relativistic descriptions have focused on aloof electron trajectories, that do not penetrate the specimen, in contrast to experimental practices, where the electron beam is typically scanned over the entire sample area [14].However, aloof electron trajectories oftentimes fail to capture intriguing phenomena associated with bulk properties, such as BPs and bulk phonons [9,18], and other sources of electron-induced photon emission, like Cherenkov or transition radiation [29,[33][34][35].
Despite the undeniable advantage of analytic solutions in data analysis, their applicability is limited to a handful of highly symmetric geometries.Over the years, differ-ent numerical schemes have been consolidated to complement analytic approaches in simulating the electromagnetic properties of nanophotonic systems [36,37].These include the boundary element method (BEM) [38,39], the finite element method (FEM) [40,41], or the finite-difference frequency-domain (FDFD) and finitedifference time-domain (FDTD) methods [42][43][44], many of which have been employed successfully for CL and EELS simulations [38,[45][46][47].An alternative route is offered by the discontinuous Galerkin time-domain (DGTD) method, which employs the Galerkin scheme to solve Maxwell's equations in the time domain [48][49][50][51][52].This method combines the flexible space discretization of finite elements, with the memory efficiency and the ability to include nonlinearities, characteristic of timedomain methods.As a consequence, DGTD offers great versatility in simulating objects of complex geometry and nonlinear response.
In this work, we present and compare an analytic approach and the DGTD method for the study of spherical nanostructures excited by aloof and penetrating electron beams.Following the work of García de Abajo for aloof electron beams [30], we derive analytic formulas for the energy loss and photon emission probability, generalized here to account for penetrating trajectories.We then validate the implementation of electron-beam excitation of nanostructures in DGTD [49] by comparing the EEL and CL spectra produced by the two methods for a perfectly spherical plasmonic NP featuring localized surface plasmons (LSPs) and BPs.Finally, we apply the numerical method to study the optical response of a NP with surface roughness, showcasing the ability of the DGTD method to emulate scenarios aligned with realistic experimental conditions that involve imperfect structures [53][54][55].

A. Analytic approach
As a first step to examine the agreement and complementarity of analytic and numerical tools, we outline the modeling of the physical system and the assumptions made in each method.As a testbed, we consider a perfectly spherical metal NP of radius R embedded in air.The NP is characterized by a unity relative permeability (µ = 1), while its relative permittivity depends on the angular frequency ω as described by the Drude model with plasma frequency ω p and damping rate τ −1 .
In the analytic calculation, the electron beam is modeled as a single elementary point charge −e, with kinetic energy on the order of keV.Since the energy transferred to plasmon resonances is typically a few eV, we invoke the non-recoil approximation [1], and assume that the electron travels with constant velocity v, say along the z-axis.Therefore, it follows a straight trajectory r e = r 0 + vt, where r 0 = (b, ϕ 0 , z 0 = −∞) is its initial position in cylindrical coordinates.The impact parameter b indicates the distance between the electron trajectory and the center of the sphere, as shown in the schematics of Fig. 1.In passing, we mention that the angular coordinate ϕ 0 need not be specified, as it does not enter the calculation, due to the symmetry of the problem.
The fast electron drives plasmon oscillations in the metal, generating an induced electric and magnetic field, E ind and H ind respectively, eventually resulting in the emission of radiation to the environment (labeled as region II in Fig. 1).Following the basic steps of Mie theory [56], one can decompose the induced electromagnetic field into transverse electric (TE) and transverse magnetic (TM) components.In region II the induced fields take the general form of outgoing spherical waves [57] where ℓ, m are the angular momentum quantum numbers, and a II ℓm /b II ℓm are the expansion coefficients of the TE/TM components corresponding to modes of electric/magnetic multipole character (see Appendix C for the full derivation).Furthermore, in Eqs.(2) h + ℓ is the spherical Hankel function of the first kind, X ℓm are the vector spherical harmonics of Eq. (A6), k 0 = ω/c is the wave number in free space, c = 1/ √ ε 0 µ 0 the speed of light in vacuum, where ε 0 and µ 0 denote the vacuum permittivity and permeability, respectively, and Z 0 = µ 0 /ε 0 is the impedance in free space.The energy radiated in the far field can be found by integrating the Poynting flux at a spherical surface of radius r → ∞, in the normal direction r.Then the (CL) probability of collecting a photon of energy ℏω is given by where dΩ denotes the infinitesimal solid angle.By inserting Eqs.(2) into Eq.( 3), and evaluating the result in the far field (k 0 r → ∞), we find where a II ℓm and b II ℓm are given by Eqs.(D8) in Appendix D. Apart from the emission of radiation, part of the energy transferred to the optical modes of the NP dissipates non-radiatively, owing to the intrinsic losses within the material.The total energy lost can be calculated by the work done by the electron against the induced field along the entire electron trajectory.Then the (EEL) probability of the electron losing energy ℏω is given by The integral in Eq. ( 5) can be decomposed into three terms (see details in Appendix D) Here, Γ bulk is related to the bulk modes of the unbound medium, reduced by the Begrenzung term Γ Begr that accounts for the presence of a boundary [58].The Γ surf term contains the contribution from modes excited by the part of the electron trajectory lying externally to the NP, and is, thus, associated with the excitation of LSPs.
The terms entering Eq. ( 6) are given by the following formulas and Here, 2z e = 2 √ R 2 − b 2 is the length of the electron path inside the NP, and γ = 1/[1 − εβ 2 ] 1/2 with β = v/c are the Lorentz kinematic factors (γ 0 is evaluated in free space).In Eqs.(7b) and (7c) K m is the modified Bessel function of the second kind and Y m ℓ are the spherical harmonics.In addition, we have set r = √ b 2 + z 2 , and θ = arccos(z/r), while analytic expressions for coefficients a I ℓm , b I ℓm , M ℓm , N ℓm , H ± ℓm , and J ± ℓm can be found in Appendices A, C and D.
It is important to note here that the aforementioned decomposition of the EEL spectra introduces a free parameter; assuming that upon losing energy ℏω the electron transfers a transverse (with respect to the electron trajectory) momentum q to excite an optical mode, q c is the maximum transverse momentum collected.In an experiment, the momentum cutoff is determined by the half-aperture collection angle φ of the microscope spectrometer, as where m e is the electron mass.We may freely choose the value for this momentum cutoff, making sure that it aligns with the typical values for the collection angle in scanning TEM (STEM) setups, which are on the order of a few mrad.Naturally, this introduces a level of arbitrariness in the EEL spectra, as different values of q c lead to different peak intensities at the BP energy.

B. DGTD simulation
We complement our analytical work with numerical simulations of the electromagnetic problem of a NP excited by a moving Gaussian charge distribution.To this end, we employ the DGTD method, which combines a piecewise polynomial spatial interpolation on an unstructured tetrahedral mesh with a Runge-Kutta time integrator to obtain a high-order accurate explicit solver for Maxwell's equations in time domain Here, j is the total current density that encompasses both any current associated with the excitation source, as well as dispersive polarization currents.The resulting method is memory-efficient compared to traditional finite elements and especially well-suited for the calculation of wide-band spectra.
One key difference between the Mie-based theory and the DGTD simulations, that can potentially lead to deviations between the two approaches, is the implementation of the excitation source.In the numerical treatment, we model the electron beam with a Gaussian charge distribution of the form with width σ e = 5 nm.This choice essentially prevents numerical artifacts that would arise in a point-charge particle modeling, while also being compatible with the typical spot size in CL experiments [8].Thereby, we introduce a new spatial scale, which needs to be considered when it is comparable to the mesh element size and the characteristic lengths of the physical system, e.g. the radius of the NP and the impact parameter.For very large distances between the electron and the NP (i.e., b − R ≫ σ e ), the source resembles a point charge, and we thus expect an excellent agreement with analytic results.However, in the opposite scenario, the finite width of the electron beam becomes important, and the corresponding fields do not accurately match those of a point charge.

A. EEL and CL spectroscopy of perfectly spherical NPs
In what follows, we analyze the response of a metallic NP excited by a fast electron beam, as predicted by both the analytic and the numerical approach.We consider a smooth sphere of radius R = 75 nm with plasma energy ℏω p = 5 eV and damping rate τ −1 corresponding to an energy of ℏτ −1 = 50 meV.The parameters mimic typical plasmonic metals and are chosen for the purpose of illustration, while the particular values have no consequences for our general conclusions.Fig. 2 shows the scattering spectrum (panel a) under plane-wave illumination, and the CL and EEL probability (panels b and c, respectively) calculated for a low electron velocity v = 0.33c (kinetic energy ≈ 30 keV) intersecting the NP at b = 35 nm from its center.Eq. ( 4) suggests that the total photon emission probability can be decomposed into contributions of pairwise orthogonal electric and magnetic multipoles of order ℓ.Fig. 2b reveals that the CL spectrum is composed of the contributions of the first four (ℓ = 1, 2, 3, 4) electric-type modes, appearing at ap- proximately 2, 2.8, 3.1, and 3.2 eV, associated with the excitation of LSP resonances, while higher-order (ℓ > 4) multipoles contribute negligibly to the spectrum.In comparison with the scattering cross section of the NP shown in Fig. 2a (calculated using Mie theory [59]), the CL spectrum provides very similar information.This is somewhat expected, since both calculations are based on the collection of far-field radiation.We observe, nonetheless, two notably distinct features in the CL spectrum of Fig. 2b.Firstly, the electron source excites more efficiently the ℓ = 2 and 3 modes, whereas in the scattering spectrum of Fig. 2a the dipolar mode peak features the highest intensity.The relative peak intensities in the CL and EEL spectra depend strongly on the impact parameter, which determines the arrangement of the polarization charges in the material [60].Secondly, we observe a small redshift of the dipolar (ℓ = 1) mode in CL due to retardation, stemming from the fact that the speed of the electron is only a fraction of the speed of light.
In Fig. 2c we present the EEL spectrum, decomposed as described in Eq. ( 6).In the zoom-in area of the figure we observe the excitation of numerous higher-order multipoles, appearing as sharp peaks at energies up to 3.5 eV.With increasing multipole order, the wavelength of the corresponding mode reduces, and higher-order modes experience the curved surface of the NP as increasingly more flat.As a result, they accumulate at the energy corresponding to that of a surface plasmon polariton (SPP) at a planar interface, at ℏω spp = ℏω p / √ 2 ≈ 3.5 eV.Above the SPP energy, the spectrum exhibits a pronounced peak at the BP energy ℏω bp ≈ ℏω p = 5 eV, pertaining to the excitation of BPs in the volume of the NP.Since BPs are longitudinal modes, they do not couple to far-field radiation and, therefore, can be detected only in EELS.At the same energy, we observe the expected negative peak related to the Begrenzung term, reducing the BP peak in the total EEL probability [58].
Having a clear picture of the origin of all spectral features, in Figs.3a and b we compare the CL and EEL probability, respectively, of the same metal NP, as calculated employing the analytic (Mie) and the DGTD method.We test the agreement between the two calculations probing various impact parameters, that range from b = 125 nm to 10 nm, corresponding to aloof electron trajectories (violet and blue spectra), grazing (green), and penetrating (yellow and red).Overall, the DGTD method reproduces the positions of the LSP modes and the corresponding CL probabilities of the Mie calculations.In the CL spectra of Fig. 3a we find an excellent agreement, with a relative error [as given by Eq. (S.11) in the SI] of around 1% for aloof and penetrating electron trajectories, according to Table I.The highest error is acquired when the electron grazes the surface of the NP, passing exactly at b = R (middle panel in Fig. 3a).This point reflects an important limit in the capabilities of the DGTD method; in the grazing trajectory, due to the finite width of the electron beam, half of the Gaussian charge density distribution lies inside the NP, while the other half lies outside, leading to numerical inconsistencies.One may avoid this point, which is inevitably difficult to resolve, by slightly adjusting the impact parameter by half the Gaussian width.
Regarding compatibility of the EEL spectra, the two rightmost panels in Fig. 3b reveal an excellent agreement for aloof electron trajectories.We consistently find a higher error for grazing trajectories (middle panel in Fig. 3b); here, the EEL spectrum calculated with DGTD exhibits a numerical artifact at the BP energy (gray dashed line), resulting once again from the fact that just a fraction of the electron charge density distribution penetrates the NP, exciting only partially the bulk mode.A substantial deviation between the two methods is found in the EEL spectra for both grazing and penetrating electron trajectories (three leftmost panels in Fig. 3b) between the SPP and the BP energy, denoted by the gray dotted and dashed lines, respectively.The disagreement is, naturally, reflected in the large relative errors presented in Table I accordingly.At these impact parameters, the condition b − R ≫ σ e is not fulfilled, hence the field of the electron deviates considerably from that of a point charge.Moreover, the beam width σ e becomes important compared to the mesh element size, since the surface of the NP is more finely discretized than the surrounding medium.Finally, there exists an additional source of error in the evaluation of the BP contribution, that stems from the rather arbitrary choice of the transverse momentum cutoff q c in the analytic approach.In contrast, in the DGTD implementation there is a respective internal limit, associated with σ e .As a result, we consistently find higher relative errors in EELS in comparison with CL, and for grazing and penetrating electron trajectories as compared to aloof.This is evident for both low-energy electron beams as in Fig 3 , as well as higher energies that are more realistic for EEL measurements.
Admittedly, the accumulation point of high-order multipoles at ℏω spp is hard to resolve in both methods.On the one hand, the analytic Mie calculation assumes a moving point charge, which can, in principle, excite an infinite number of multipoles, resulting in a sharp highintensity peak at energy ℏω spp .On the other hand, the finite mesh size and beam width implemented in DGTD imposes a limitation to the number of multipoles that can be resolved for a given discretization, since high order multipoles associated with field variation shorter than the mesh element size at the surface cannot be captured without the use of very high order polynomials.In EELS and CL experiments, there exists an analogous limitation, associated with the finite width of the electron beam employed, as well as the geometric imperfections of the NP.The versatility of DGTD allows us not only to adjust the electron beam width according to the experimental setup, but also to mimic NPs with surface roughness, as we discuss in Section III B. In the analytic approach, the smearing of higher-order modes and the overall quenching of the sharp peak at the accumulation point can too be reproduced, once we consider the nonlocal response of the material; this can be done particularly easily within Mie theory [61][62][63].Nonlocal effects manifest as increased damping and uneven energy shifts of high-order modes, and, therefore, lead to the suppression of the individual modes, as well as the reduction of their overlap at ℏω spp [64].
The difficulty in resolving the high-order multipoles even in the analytic calculation is clear in the convergence study presented in Figs. 4 and 5. Due to multipole orthogonality, the sole contribution of increasing ℓ max is to amplify the peak at ℏω spp .This suggests monitoring the total EEL spectrum integrated over energy (area under the curve) as a proxy for convergence.For aloof trajectories, the asymptotics of the Hankel functions suggest exponential convergence, which is in agreement with Figs.4a and b; the EEL probability converges at ℓ max = 10.In contrast, Fig. 4c shows that in the case of the grazing trajectory the convergence order breaks 1.5 2.5 3.5 4.5 5.5 Energy (eV) down to a square root law.By extrapolation (magenta line) we can conclude that even for ℓ max = 63 the analytic result is still converged only up to about 15%.Fig. 4d corroborates that the area missing from the converged value corresponds to the higher-order modes piling up at the SPP energy.Finally, we note that Fig. 4a exhibits the same square root convergence before the curve flattens off.The exact value of ℓ max where this transition happens increases as the impact parameter approaches R.
For the penetrating trajectories shown in Fig. 5 we find signs of the same slow convergence around the SPP accumulation point.This is masked in our convergence plots by the fact that the BP peak diverges.The source of the divergence lies in the decomposition of the EEL probability into two competing contributions stemming from the bulk and the Begrenzung terms.As illustrated in Figs.5a and c, the two terms produce divergences of opposite sign; the negative Begrenzung term diverges linearly for increasing multipole order ℓ max , whereas the positive bulk term diverges logarithmically for increasing momentum cutoff q c [see Eqs.(7a) and (7c)].As computational resources do not allow driving ℓ max and q c to infinity, once one of the two parameters is truncated to a certain cutoff value, the other has to be adjusted accordingly.Figs.5b and d show that the different values of ℓ max and q c , respectively, affect the spectra for energies above the SPP mode.

B. EELS of NPs with surface roughness
The synthesis of metallic NPs, as the ones studied in the present work, is routinely done with colloidal chemistry, for a large variety of materials and NP shapes [65,66].However, despite being able to accurately control the NP size, assuring a smooth surface is rather challenging.Typically the structures exhibit protuberances on the surface, which can be responsible for symmetry breaking [67,68], hot spots in dimers [69,70] and picocavities [71,72], or energy shifts of the LSPs [73].Within the DGTD method, surface texture can be easily implemented on top of the perfectly spherical mesh and incorporated in the numerical calculations.Here, we follow the prescription presented in Ref. [55] to implement the desired roughness.Given an initial smooth sphere of radius R, we introduce deviations from its nominal value in both the normal and the lateral direction, which are defined by the root-mean squared (rms) roughness value and the correlation length l, respectively.The rms parameter is related to the local variations of the radius that follow a Gaussian white noise distribution, while l describes the average distance between neighboring bumps on the surface.
In Fig. 6 we explore the effect of surface roughness on a spherical NP, described by the same Drude permittivity as in the previous section, now excited by an electron beam traveling with velocity v = 0.7c (kinetic energy ≈ 200 keV) at distance b = 100 nm from its center.We probe meshes of two degrees of surface roughness on top of the NP of nominal radius R = 75 nm, namely rms = 2 nm and 4 nm, while the correlation length is fixed at l = 10 nm in both cases.Since the breaking of the spherical symmetry introduces a dependence on the electron propagation direction, and on the mesh morphology, in Fig. 6 we plot the average EEL probability, corresponding to the average values obtained for 6 different meshes.The resulting spectra for rms = 2 nm and 4 nm (dark red and blue curves, respectively) deviate notably from that of a smooth sphere (gray shaded area).Firstly, we observe an increasing redshift of the spectra with increasing degree of roughness, in agreement with experimental observations of corrugated plasmonic NPs [73].The energy shift is most evident for the dipolar mode at around 2 eV, and is the result of the area increase of the rough NP as compared to the smooth one.Evaluation of this area from the mesh parameters yields an effective radius of R eff = 76.2 nm for rms = 2 nm, and R eff = 79.6 nm for rms = 4 nm.Indeed the EEL spectra of smooth spheres of said effective radii reproduce accurately the position of the dipolar mode (see Appendix E).In addition to the redshift, the spectra of the corrugated NPs feature a large number of low-intensity peaks.These new spectral features are the result of two factors.Due to the breaking of the spherical symmetry, the prior degenerate modes associated with the same angular momentum ℓ but different m number, now exhibit a small energy difference.As a result of the lift of the degeneracy, the sharp peaks observed in the spectrum of the smooth NP are suppressed and, depending on the size of this energy difference with respect to the linewidth of the degenerate mode, they are either split or broadened.Moreover, additional spectral features may arise from hot spots, namely protuberances of large curvature that strongly enhance and confine the incident field [72].It is important to note here that, as the features become increasingly smaller, a rigorous description of the system requires the implementation of nonlocal effects in the method [74,75], which effectively introduces a cutoff in the contribution of large-wavevector components [76,77].

IV. CONCLUSIONS
We have presented an analytic and a numerical method for the study of spherical structures excited by fast electron beams.Based on Mie theory, we have derived formulas for the calculation of the EEL and CL probability that are valid for both aloof and penetrating electron beams, as is typically the practice in EEL and CL measurements.Focusing on the plasmon oscillations of a metallic NP as a testbed, we compared the analytic theory with numerical simulations performed using the DGTD method, and found excellent agreement.We discussed the applicability and limitations of each method, particularly for grazing trajectories and at energies near the surface-and bulk-plasmon resonances.Finally, we showcased the flexibility of the DGTD method by studying a NP with different degrees of surface corrugation, which can lead to resonance shifts and splittings due to the lifting of mode degeneracy.We thus believe that both methods are essential and complementary for exploring collective optical excitations in matter, and for interpreting experimental observations.
where k 0 is the wave number in free space, and the scalar functions ψ M/E satisfy the expressions In the absence of charges and currents, the electric field satisfies the vector wave equation ∇ 2 + k 2 0 E(r, ω) = 0.As a result, the scalar functions satisfy the scalar wave equation whose Green's function in spherical coordinates is given by [57] where Y m ℓ are the scalar spherical harmonics.In Eq. (A4), we have introduced the angular momentum quantum numbers ℓ and m, and r </> stands for the smaller/greater between two points r and r ′ .Eq. (A3) admits standing and outgoing spherical waves as solutions, represented by the spherical Bessel function j ℓ , and the spherical Hankel function of the first kind h + ℓ , respectively.Thus, we expand the electric field at points r < |r e | as follows where X ℓm (θ, ϕ) are the vector spherical harmonics, defined as with X 00 = 0.The expansion coefficients in Eq. (A5) are given by and with γ 0 = 1/ 1 − β 2 and β = v/c being the Lorentz kinematic factors, and K m the modified Bessel function of the second kind.In Eqs.(A7) we have set taking relativistic and retardation effects into account.
Even though in the present work we focus on metallic NPs that sustain LSP and BP modes, the expressions derived are general and can be used to study the response of any linear and isotropic material subjected to EEL and CL spectroscopy.We consider the spherical NP of radius R investigated in Section III, with relative permittivity ε and relative permeability µ = 1.The NP is excited by an electron beam that passes through the bulk of the material, entering at point (b, ϕ 0 , −z e ) and exiting at (b, ϕ 0 , z e ) in cylindrical coordinates with respect to the NP center, as illustrated in Fig. 7a.We assume that the recoil of the fast electron in a single scattering event is negligible due to its high velocity, and we therefore work within the non-recoil approximation [58].
To compute the total electromagnetic field arising from the interaction of a fast electron with a spherical NP, we follow the basic steps of Mie theory [59].According to it, the fields inside and surrounding the NP are expanded in terms of spherical waves of unknown expansion coefficients, which are then found by matching the fields at the surface of the NP according to appropriate boundary conditions.We first express the direct field of the electron E 0 in terms of spherical waves.Clearly the electromagnetic fields generated by the electron while moving within the NP are different to the ones generated when it travels in the surrounding medium.The electric field E II 0 corresponding to the external trajectory (for R < r < |r e |) is given by Eq. (A5), with modified expansion coefficients and with , and θ = arccos(z/r).In Eq. (C1b) we have set that holds for any type of spherical Bessel function f ℓ (k n r) evaluated in any medium n (the prime here and on any other Bessel function denotes the derivative of the function with respect to the argument).In particular, in Eq. (C1b) we have F ± ℓm (k n z) = H ± ℓm (k 0 z) and f ℓ = h + ℓ .The electric field generated by the fast electron while traveling inside the NP (for |r e | < r < R) acquires the form of outgoing spherical waves, as where k = √ εk 0 is the wave number inside the NP, and the expansion coefficients are given by b 0,I ℓm = − B .Similar expressions hold for the magnetic field H, which can be obtained from Faraday's law H = ∇ × E/(iωµ 0 ).Therefore, the total fields inside (I) and surrounding the NP (II), close to the boundary, take the form where Z = µ 0 /(εε 0 ) is the impedance inside the NP.At the interface r = R, the total fields must satisfy the boundary conditions that require the continuity of the tangential components of the E and H field [57], leading to the following set of equations Using the orthonormality relations of the vector spherical harmonics [57], we obtain the unknown expansion coefficients We now have all ingredients required to evaluate the fields of Eqs.(C5), and calculate the EEL an CL probability.

Appendix D: EEL and CL probability EEL Probability
The fast electron loses part of its kinetic energy due to energy transfer to the optical modes sustained in the structure; here, the plasmon oscillations sustained in metals.The probability of the electron losing energy ℏω is related to the work done by the electron moving against the induced field along the electron trajectory [30], and is given by Eq. ( 5), or reshaped as dz Re e −iωz/v ẑ • E ind (r e , ω) .

(D1)
The induced field can be found as the total field minus the electron field E air 0 in the absence of the structure, i.e.E ind = E tot − E air 0 .Eventually, in the two regions the induced field is E II ind = E II B , and and Eq.(D1) yields where we have split the integral running along the electron trajectory, into the path of the electron within and outside the NP.We further decompose the EEL probability into three terms as follows In the decomposition presented in Eq. (D3), the Γ bulk term corresponds to the energy lost to the excitation of bulk modes of the unbound medium modified and reduced by the Begrenzung term Γ Begr , that accounts for the presence of a boundary.The remaining terms are grouped together and referred to as the "surface term" Γ surf , to indicate that they only contain contributions to the EEL probability pertaining to the excitation of LSPs (see Fig. 2c).We stress, however, that said name serves classification purposes here, and it might not be appropriate in the study of non-plasmonic materials that support modes of different characteristics.For high-index dielectric nanospheres, for instance, that host Mie resonances in their volume, the description "surface" is not appropriate.
To derive the bulk contribution to EELS, we make use of Eq. (B2), which determines the field generated by the fast electron traveling distance 2z e within an infinite medium of relative permittivity ε and relative permeability µ = 1.The corresponding EEL probability Γ bulk (ω) can be found via the momentum-resolved EEL probability P EEL (q, ω) Then the EEL probability Γ EEL (ω) is obtained by leading to Eq. (7a).The isolation of the bulk term inevitably leads to the introduction of the free parameter q c ; considering that the electron transfers transverse (with respect to the electron trajectory) momentum q upon losing energy ℏω, we impose a cutoff q c to the maximum transverse momentum collected, in order to ensure finiteness of the EEL probability [58].Its value can be determined by the half-aperture collection angle of the microscope spectrometer, according to Eq. ( 8).
Inserting the fields of Eqs.(C5) in Eq. (D3) we obtain where we have further decomposed the surface and Begrenzung terms into contributions from electric (E) and magnetic-type (M) multipoles, that take the form where the notation a/b 0,I ℓm | air indicates evaluation of the terms given by Eq. (C4) in air (ε = 1, k = k 0 ).This correction ensures a zero probability in the absence of the structure.
For an aloof electron trajectory, we set z e = 0 and find the expressions presented in Refs.[1,30].The final result agrees with the expression for Γ CL presented in Ref. [79].
As a final comment, we mention that the EEL and CL probabilities in Eqs.(D5), (D7) and (D9) are given in units of seconds (s), whereas in all relevant figures in the main text they are presented in units of inverse energy (eV) −1 .To perform the unit conversion one can simply divide by the reduced Planck constant ℏ. color).By evaluating the area from the mesh parameters we obtain an effective radius of R eff = 76.2 nm for rms = 2 nm, and R eff = 79.6 nm for rms = 4 nm.In order to confirm our argument, in Fig. 8 we juxtapose the EEL probability calculated analytically for a smooth NP of radius equal to the calculated R eff in each case.Indeed, the redshift of the lower-energy peak (electric dipolar mode) is very well reproduced by smooth NPs of the same area.Due to the symmetry breaking, higher-order modes exhibit a more complex behavior; additionally to the shift, they also feature energy splits.

FIG. 1 .
FIG. 1. Schematic illustration of the geometry under study; a metallic NP of radius R = 75 nm, and permittivity described by the Drude model of Eq. (1), is excited by an electron beam passing with velocity v at impact parameter b with respect to its center.The color dots correspond to the selection of impact parameters examined in our study.Labels I and II indicate the region inside the NP and the surrounding medium (air), respectively.

FIG. 2 .
FIG. 2. (a)Scattering cross section, normalized to the geometrical cross section πR 2 , of the metallic NP of Fig.1under plane-wave excitation (gray shaded area).The spectrum is decomposed into the contributions from the first three multipoles (blue, red and green curves, for ℓ = 1, 2, 3 respectively).In panels (b) and (c) we consider an electron traveling with velocity v = 0.33c (kinetic energy ≈ 30 keV) passing through the NP of Fig.1at impact parameter b = 35 nm.(b) CL probability (yellow shaded area) and multipole decomposition (blue, red, green, and violet curves).(c) EEL probability (yellow shaded area) decomposed into the surface, bulk, and Begrenzung contributions (green, pink and blue curves), with the inset focusing on the energy window 1 − 4 eV.The bulk term is given by Eq. (7a) for qc = 0.71 nm −1 and the surface and Begrenzung terms are obtained by Eqs.(7b) and (7c) for multipole order cutoff ℓmax = 63.

2 FIG. 3 .
FIG. 3. (a)CL and (b) EEL probability in the interaction between the metallic NP of Fig.1and an electron traveling with velocity v = 0.33c (kinetic energy ≈ 30 keV) considering 5 different impact parameters, as denoted in the labels and illustrated in the insets of each panel.The shaded areas correspond to results of the analytic Mie-type calculation, whereas the solid lines show the corresponding DGTD simulations.The vertical gray dotted and dashed lines in the three leftmost panels of (b) trace the energy of the SPP (ℏωspp ≈ 3.5 eV) and the BP mode (ℏω bp ≈ 5 eV), respectively.In the two leftmost panels of (b) we use qc = 0.71 nm −1 , while ℓmax = 63 in all panels.

FIG. 4 .
FIG. 4. (a, c) Area under the curve of the EEL spectra of Fig. 3b, (a) for the aloof electron trajectory at b = 125 nm, (c) for the grazing trajectory at b = 75 nm, showcasing the convergence of the EEL probability for increasing values of the multipole cutoff ℓmax, plotted versus 1/ √ ℓmax.In both panels, the magenta line is fitted to the data points marked as black bullets, while the red crosses represent data points excluded from the fitting.The vertical gray lines serve as guides to the eye for the position of ℓmax = 1, 2, . . .9. (b, d) EEL probability calculated for impact parameters (b) b = 125 nm, and (d) b = 75 nm, and for selected values of ℓmax, as denoted in the labels.

2 FIG. 5 .
FIG.5.(a, c) Area under the curve of the EEL spectrum of the leftmost panel in Fig.3b(b = 10 nm), showcasing the divergence of the EEL probability for increasing values of (a) the multipole cutoff ℓmax, and (c) the transverse momentum cutoff qc plotted versus ln(qc).The data follow a linear divergence trend with respect to their corresponding horizontal axes.In both panels, the magenta line is fitted to the data points marked as black bullets, while the red crosses represent data points excluded from the fitting.(b, d) EEL probability calculated at b = 10 nm for selected values of (b) ℓmax, and (d) qc.In panels (a, b) we scan over ℓmax, while keeping a fixed value qc = 0.7 nm −1 , whereas in (c, d) we scan over qc for a fixed value ℓmax = 63, as denoted in the labels.In panel (c) the vertical lines correspond to the selected qc values presented in (d), following the same color coding.

FIG. 6 .
FIG.6.Average EEL probability in the interaction between a spherical NP featuring surface roughness and an electron beam passing with velocity v = 0.7c (kinetic energy ≈ 200 keV) at distance b = 100 nm.The solid lines correspond to NPs, whose shapes deviate from that of a perfectly smooth sphere of radius R = 75 nm (gray shaded area) by root-mean square roughness values rms = 2 nm (dark red curve) and rms = 4 nm (dark blue curve).The average EEL probability corresponds to the average values obtained with DGTD for 6 different rough meshes characterized by the same rms.

FIG. 7 .
FIG. 7. (a) Schematic illustration of the electron beam traveling with velocity v = vẑ, and traversing the spherical NP of radius R. In the illustration, we assume that the electron travels on the x-z plane, entering at point (x = b, z = −ze) and exiting at (x = b, z = ze).(b) Schematic illustration of an electron beam of velocity v = vẑ traversing the bulk of an unbound medium of optical parameters ε and µ = 1.

2 .
CL probability Upon coupling to the propagating electron, the radiative excited modes decay by emitting radiation to the far field.The probability Γ CL of emitting a photon of energy ℏω is obtained via the Poynting flux of Eq. (3).By inserting the induced field E II ind = E II B , and letting r → ∞, (D9) As in the EELS calculation, we remove the contribution of the direct electron field by modifying coefficients a II ℓm and b II ℓm according to Eqs. (D8).

FIG. 8 .
FIG. 8. Average EEL probability in the interaction between a spherical NP featuring surface roughness and an electron beam passing with velocity v = 0.7c (kinetic energy ≈ 200 keV) at distance b = 100 nm.The solid lines correspond to NPs, whose shapes deviate from that of a perfectly smooth sphere of radius R = 75 nm by (a) rms = 2 nm (dark red curve) and (b) rms = 4 nm (dark blue curve).The pink (panel a) and cyan areas (panel b) show the analytically derived EEL probability for a smooth sphere of effective radius (a) R eff = 76.2 nm, and (b) R eff = 79.6 nm.

TABLE I .
Relative error between the analytic and the DGTD calculations of ΓCL and ΓEEL for varying b.