Vibrational surface EELS probes confined Fuchs-Kliewer modes

Recently, two reports have demonstrated the amazing possibility to probe vibrational excitations from nanoparticles with a spatial resolution much smaller than the corresponding free-space phonon wavelength using electron energy loss spectroscopy (EELS). While Lagos et al. evidenced a strong spatial and spectral modulation of the EELS signal over a nanoparticle, Krivanek et al. did not. Here, we show that discrepancies among different EELS experiments as well as their relation to optical near- and far-field optical experiments can be understood by introducing the concept of confined bright and dark Fuchs-Kliewer modes, whose density of states is probed by EELS. Such a concise formalism is the vibrational counterpart of the broadly used formalism for localized surface plasmons; it makes it straightforward to predict or interpret phenomena already known for localized surface plasmons such as environment-related energy shifts or the possibility of 3D mapping of the related surface charge densities.

that the local dielectric constant (ω) = (ω, q = 0) (where ω is the energy and (ω) is equal to its value at zero transferred momentum q) is sufficient to describe electromagnetic excitations in a finite system. In Ibach's simple geometry, ω s was such that (ω s ) = −1.
Fuchs and Kliewer demonstrated the amazing efficiency of the LCDM to describe more complicated geometries, such as slabs [11] and infinite cylinders [12]. Already in these simple The Fuchs-Kliewer work has been extended with an impressive success [10] to the description of surface plasmons (SP) in simple systems such as slabs and cylinders [12,13] (see Figure 1b). Stimulated by the development of the research on plasmons in nanoparticles systems, several simulation schemes basically relying on the LCDM (boundary element model, BEM [5,14,15] and discrete dipole approximation [16]) have been extensively used to simulate optical and EELS spectra dominated by localized SPs confined on nanoparticles. BEM simulations have been recently extended to the phonon range for STEM-EELS [2] using the MNPBEM [17] implementation. Now, beyond their unique simulation capabilities, LCDM derived theories have offered a deep understanding of localized SP physics. In particular, they made explicit the link between STEM-EELS and optical near-field spectroscopies as both are related to the elec-tromagnetic local density of states (EMLDOS) [6,7], and showed that EELS is related to the extinction cross-section for dipolar modes [18,19].
The goal of this paper is to show how the reasoning once made to explain SP confinement in nanoparticles and interpret STEM-EELS experiments can now be used to rationalize the interpretation of surface STEM-EELS vibrational experiments in nano-objects and predict new physical effects.
In the following, we will introduce the confined FK (cFK) modes as surface phonons whose properties are mostly defined by the classical confinement that they experience in particles much smaller than the free-space equivalent wavelength. In this sense, if normal phonon modes are conceptually related to bulk plasmon modes and FK modes (also known as surface phonons) to surface plasmons, cFKs are the phononic counterpart to localized SPs. For the sake of simplicity we will neglect retardation in the following. As we will show, this is justified by the relatively small sizes of phononic nanoparticles studied in the literature [1,2]. A rigorous definition of the cFK modes can then be given in the quasi-static (QS) approximation using a modal decomposition form, first introduced in the case of confined SP [4,5,7], see Annex. cFKs are then defined as a set of eigencharges σ i and eigenvalues λ i , i being the mode index. In the general case, λ i , which depends only on the geometry of the nanoparticle, has to be determined numerically, and corresponding eigenenergies can be deduced through a simple implicit relation between λ i and the energy dependent dielectric constant (see Annex). In the case of a model Drude-Lorentz dielectric constant, a general expression for the cFK eigenergies is (see Annex): cFK energies lie between the bulk LO and TO energies, as −1 < λ i < 1 [4], and we directly see that the energy of two well-known FK modes for an infinitely thin slab, describing the charge-antisymmetric and -symmetric modes (see Figure 1a), are retrieved for λ i = ±1. In addition, other simple cases can be straightforwardly deduced. λ i = 0 corresponds to the above-mentioned surface phonon [9] in a Drude-Lorentz model, λ i = −1/3 [20] corresponds to the dipolar mode of a sphere ( = −2, ).
To exemplify the interest of this approach, we start with the case of nanorods that has been widely investigated in surface plasmon physics [21], and especially by EELS [22,23].
The simplicity of the structure makes it easy to understand the intimate link between shape and modes structures, and we adapt it here to the case of a phononic material following arguments for localized SPs found in [20]. Modes in a nanorod of radius r and length L are similar to the FK modes of the infinite rod, except that the confinement restricts the available wavevectors to multiple of 1/2L. This is exemplified in Figure 1c where the discrete modes dispersion relation, simulated for a large set of nanorods lengths, overlaps the one of an infinite rod. Such modes are the cFK modes of the nanorod. The cFK modes disperse between ω T O and ω s , in analogy with the corresponding dispersion in localized SP in nanorods restricted between 0 and ω sp [20]. Similarly to the corresponding localized SP modes, each mode with eigenvalue λ i corresponds to an oscillation of the surface eigencharge, as depicted in Figure 1d. We see on the prototypical case of a nanorod that the QS approximation is much more justified for cFK than for localized plasmons for objects of same sizes: the length (top scale in Figure 1d) of a typical nanorod is much smaller than the equivalent free-space wavelength of the cFK (right scale in Figure 1d). Another difference is the pile up of low order modes for long nano-antennas close to ω T O which is obviously absent for localized surface plasmons.  [24] reveal a series of peaks. As seen on Table I, a direct comparison of their energy values with that of the cFK deduced from Eq. 1 based on the sole knowledge of the λ i , ω T O , ω LO and ∞ shows an almost perfect agreement. This validates conceptually our approach, and also allows to use a simple EELS modal decomposition (see Annex 3) for EELS simulations.
In Figure 2a, we also compare EELS to macroscopic optical quantities such as the absorption, extinction and scattering cross-sections calculated in the retarded approximation.
As in the case of EELS, the spectra do not peak at the normal modes energies ω LO and ω T O . Instead, they are dominated by the cFK modes, in analogy with the well-known case of a slab spectrum dominated by the FK modes [11] or more generally for an ensemble of nanoparticle [25]. This is particularly justified from the modal decomposition of the cross-sections, see Eq. 4 in Annex and [18]: the optical cross-sections are proportional to a spectral function peaking at the dipolar cFK modes energy. Contrary to the case of EELS, only the dipolar modes are observable (but a very slight contribution from the third order mode). The spectra obviously show a large dependence on the incoming polarization. For polarizations along the nanorod axis, the dipolar mode of the low energy branch is excited.
For a polarization perpendicular to it, the dipolar modes of the other branches, almost all arising at ω s [26], are excited, see Figure 2a. This points to the fact that EELS is sensible to both bright (i. e. optically active) and dark (i.e. not optically active) cFKs, in contrast to optical far-field techniques.
Obtaining truly dark (non-emitting/absorbing) localized SPs is difficult due to the relatively large sizes of plasmonic particles [18] with respect to the corresponding free space wavelengths. In contrast, for the cFKs where the QS approximation is justified for much larger particles sizes, almost only dipolar modes are bright. We note that the scattering cross-section is several order of magnitude smaller than the extinction one. This is basically related to the fact that, other things being equal, the ratio between scattering and extinction scales as 1/ω 3 , where ω is the energy of interest. This makes extinction and absorption cross-sections almost identical at the low energy of the phonon regime, making EELS very close to the absorption cross-section for dipolar cFK modes (see also the analytical proof in the Annex). We note that this contrasts with the case of a silver plasmonic nanorod of the same size (see Figure 3). In this case, scattering has a major contribution in the extinction cross-section.
We can now clarify the type of selection rules when exciting cFK optically or with electrons. To start with, in the QS approximation, only dipolar modes can be excited by a plane wave, and the electrical polarization of the plane wave must be aligned with the dipole direction. Away from the QS, similar symmetry arguments arise: even modes (mode 2 and 4 on Figure 1d) cannot be excited by a plane wave with electrical field in the plane containing the axis of the nanoantenna, while odd modes (1 and 3) can be excited. Tilting the beam direction with respect to the antenna axis will break the symmetry and make it possible to also detect even order modes. More generally, for optical experiments, the selection rules are completely determined by the general symmetry of the surface charge distribution with respect to the plane wave direction and polarization.
The interplay between the symmetries of the incoming electron electrical field and the surface eigencharges is different. As with optics, cFK modes are also probed by EELS, but contrary to optics, EELS is sensible to all modes even in the QS approximation. Also, the symmetry of the surface eigencharges impacts rather the spatial distribution of the EELS signal. Indeed, EELS maps ( Figure 2b) closely resemble the EMLDOS projected along the electron propagation direction z (zEMLDOS, Figure 2c), with the EMLDOS spatial and spectral distribution being essentially determined by the size, shape and symmetries of the object of interest. The resemblance between EELS and zEMLDOS is expected by analogy with the localized SP case, where also a general analytical relation between these two quantities can be determined [6]. Much as in the case of localized SPs [7], EELS as well as near-field optical techniques do not map directly the eigencharges [27]. Rather, they map the related zEMLDOS, itself related to the z-projection of the electric eigenfield in the QS limit [7,15]. An even more precise description of EELS of cFK in terms of electromagnetic quantities is given by the almost identity between EELS and the z-integrated eigenpotentials [28], see Figure 2d.
We can sum up the results exemplified on the nanorods but valid for any kind of phononic nanoobject.
First, surface EELS and optical IR absorption, extinction and scattering are probing the same physical excitations, namely cFK. The symmetry of the cFK surface eigencharges, which depends on the global shape and symmetry of the subtending particle determines the coupling strength of the cFK with the probing electrons or photons. This is in stark contrast with IR absorption or bulk EELS [2,29,30], which are probing normal modes, which depend on local (atomic) symmetries, i. e. the bulk material properties. This is also a main difference between our work, which relates surface vibrational EELS to the concept of EMLDOS, and recent theoretical works describing the link between bulk EELS to the concept of phononic density of states (pDOS ). Again, pDOS is dependent on the atomic structure symmetry while EMLDOS is dependent on the global (shape) symmetry of the nanoparticle. Also, for similar reasons, surface EELS is completely different to Raman spectroscopy which probes bulk properties of atomic oscillations, although following selection rules different to that of bulk IR absorption. Note that the LCDM can also be used to predict the bulk EELS experiments results through a term proportional to −Im(1/ (ω)), giving essentially a peak at ω LO in the Drude-Lorentz model. The intensity of the related peak maybe influenced by the screening at the surface, a phenomenon handled in the LCDM theory and known as "begrenzung" effect [2]. There are however several limits explaining the need to develop dedicated theories for bulk phonons beyong the LCDM [2,29,30], related to the interpretation of angular resolved experiments [2,29,30].
Second, EELS maps are close to that obtained with the near-field optical measurement which are related to the EMLDOS [31], and map quantities close to the cFK electric eigenfields, and more precisely the eigenpotentials, along the electron direction integrated on the electron beam path (see an analytical proof in Eq. 3 and [28]). The typical spatial extent of the EELS signal is related to that of the EMLDOS, and almost identical to that of the integrated eigenpotentials.
Third, due to the large free space wavelength of the cFK compared to typical dimensions of nano-objects, the QS approximation holds essentially true for sub-micron nanoparticles, and any nanoparticle can be described by a series of eigencharges and related λ i that only depends on the shape of the nanoparticle.
In addition, this theory works well for understanding cFK, but will obviously fail to describe long-wavelength, propagating surface phonons that may arise in the particular case of very large particle or slabs. In the case of slabs or infinite cylinders, however, alternatives rigorous retarded theories exist [11]. The differences in the predictions between a quasi-static (such as presented here) and retarded formalism weakly affect lowest energy, charge-symmetric modes that are usually dominant in slabs and cylinders. The resonances energy did not change as a function of the electron beam position whether it was impinging the objects or in vacuum close to them. The 173 meV resonance was attributed to the LO normal mode of hBN, and the other compared to IR results without further assignment. Following the reasoning of this paper, one can rationalize these results, see also Table II. The 173 meV (hBN) modes and 138 meV (SiO 2 ) are likely to be chargesymmetric (lower branch in Figure 1, λ i close to −1) FK modes. Indeed, with the help of equation 1 (see Table II), one can directly deduce that their energies are between the ω T O and ω s (and very close to ω T O = 169.5 meV in the case of hBN) but largely different from ω LO , see Table II. For symmetry reasons, the dipole strength of the charge-antisymmetric mode vanishes with the thickness of the slab [34]. It might explain why this mode was not reported in [1]. On the other hand, as summarized in Table II, Batson and Lagos [35] reported the measurement of two peaks on an h-BN flake, the first at 187 meV (below ω s ) and the second at 203 meV (above ω s ). These are likely to be charge symmetric and antisymmetric modes reciprocally -as confirmed by preliminary simulations in [35]-for a slightly thicker slab (as the symmetric mode energy is at higher energy and the symmetric mode is still weaker but now measurable). It is worth noting that in these cases, the energy of the modes depends on the geometry and symmetry of the nano-object, and we expect of course the observation of thickness dependent modes when more experimental works will be available in the literature. Finally, no modes energy spatial variation has been reported on these two sorts of slabs [1,35]. Recently, Schmidt et al. [36] showed that the plasmonic modes in thin objects with edges can be decomposed in slabs modes and edges modes independently. The slabs modes follow the infinite slabs dispersion relations, and edges the nanoantennas ones [37]. The modes of lowest energy branches have the same charge symmetry with respect to the slab or cylinder mid-plane, so that the slab and edge lowest energy modes share the same symmetry. Translated to surface phonons in SiO 2 slabs it means that we should expect two different modes of same symmetry with respect to the slab mid-plane; however, both dispersion curves are very close (see e.g Figure 1a), and for very thin objects both slabs and edge modes energy tend to a unique and same value (ω T O ), making it difficult to detect experimentally any spectral or spatial variation except an intensity decrease in vacuum.  Table III how well Eq. 1 reproduces our simulations and Lagos's ones, themselves pointed to be in very good agreement with experiments ( [2]). Our theory gives however a stronger insight into the nature of the probed modes. In Lagos et al. [2], modes are denominated through their EELS spatial distribution, with no discussion on their symmetries, which are known to be complex for cubes plasmons [38,39]. Indeed, as shown in Figure 5, the corner mode can be decomposed in dipolar, quadrupolar and octupolar contributions (see also Table III) that are degenerated in the quasistatic approximation. Because one of its components is dipolar, the corner mode is likely to be bright (ie theoretically measurable through an IR extinction experiment) although weakly scattering compared to a plasmonic cube of the same size.
Quite interestingly, the edge mode is in fact composed of a large number of cFKs of close λ i , see Table III and SI. The symmetry of all these constituting modes makes the edge mode a dark one. Concerning the face mode, the number of polygons required for convergence was to high to deduce a definite value or set of values for λ i . However, this highest energy mode has an energy very close to ω s for MgO, corresponding to λ i = 0 (see Annex). This is expected from localized SPs analogy, as high momenta modes converge systematically to this value.
We now turn to a point which has not been considered so far but may have important implications for the interpretation of the forthcoming experiments. Indeed, the effect of the substrate, known to be essential in plasmon physics, has not been discussed in the context of surface vibrational STEM-EELS experiments. It is well-known that localized SP energy and spatial distribution drastically depend on the close presence of other materials, like a substrate or an embedding matrix. In Figure 4b, we show the effect of embedding a phononic nanorod into a material of constant dielectric constant different to one. It produces an expected redshift of the excitation, yet still constrained between ω T O and ω LO . The case of a nanoparticle on a substrate is more subtle. In particular, in the case of a nanocube, it is well-known from localized SP physics that the modes will split into modes at low energy localized close to the substrate (proximal modes) and at higher energy close to the vacuum (distal modes) [38]. In [2], only the distal modes were reported, although both types of modes are actually predicted (see Figure 4). We note that the distal modes energies are very close to the mode of a free space cube, explaining the good agreement between our theory, Lagos and our simulations without substrate, and experimental results. Observation of the proximal band would however require a spectral resolution even better than actually available.
Finally, the theory presented here can be extended to understand more complicated situations. This is in analogy with the success of the theory presented for localized SPs [4][5][6][7]19], which has been extended to the 3D mapping of the EMLDOS [40] or of the surface eigencharges [8], the simulation of the cathodoluminescence signals [18,19], the interaction of surface excitations with phase-shaped incoming beams [27], or the coupling between localized SP. Also, this model can be refined by developing a retarded model or a non-local approximation extension [41].
. From this set, that can be determined numerically [4,5,17], one can deduce all eigen-quantities such as the eigenpotential or the electrical eigenfield E i ( r) at all points r, or any observable such as the EMLDOS ρ αα ( r, ω) (here α represents the projection direction): the EELS probability (simplified here to the case where the beam is outside of the object of interest) [7]: where v is the speed of the electron, z the direction of electron propagation, R ⊥ the position of the beam in the plane perpendicular to z and the extinction cross-section, which is equal to the absorption cross section in the QS limit reads [18]: where A i is a mode dependent prefactor, and the sum runs over the dipolar d cFK modes only. (ω) = ∞ (1 + then the spectral function then takes the simple form of a lorentzian peaking at the cFK mode energy ω i (solution of equation 1, this is the energy of the ith cFK in absence of dissipation), weighted by some energy independent prefactor.
EMLDOS, EELS and absorption cross-section can be straightforwardly deduced from this expression of the spectral function.
The above deductions can be extended analytically to the case where the object of interest is embedded in a medium. Similar developments (see SI of [19] or [40]) can be done in the retarded regime assuming a model dielectric function.    The case of hBN is a bit more involved, as hBN is an uniaxial anisotropic material.

Analogy between localized SP and cFK modes
Nevertheless, the FK theory can be extended to anisotropic materials for slabs [12]. The charge symmetric mode converge to the in-plane TO mode energy ω T O ⊥ and the charge antisymmetric to the out-of-plane LO mode energy ω LO || . The terminology ⊥, || is related to the anisotropy axis. Likewise, the surface phonon energy will be a combination of in-plane and out-of-plane phonons energy given by the condition √ ⊥ || = −1, with ⊥ and || the in and out of plane dielectric constant [42]. We note that an HREELS study [43] reported a value for the LO mode of a single hBN sheet around 173 meV, much as the value reported by [1].
Given the similarities pointed out in the paper between HREELS and STEM-EELS and the symmetry arguments, the reported LO mode is most likely to rather be a charge-symmetric FK mode.   , from retarded simulations with experimental dielectric constant found in [24], from retarded simulation in [2] and experimental results from [2]. Inputs for equation 1 are ω T O = 50.7 meV, ω LO = 91.3 meV, ∞ = 3.01 [24]. Energies are given in meV units. Note the apparent discrepancy for the face mode values between simulations and experiments, proven in [2] to be an effect of finite spectral resolution in the experiments.

Simulations
Dispersion relations in Figure 1a and b have been calculated using formulas from [46] and using a Drude model adapted to silver and a Drude-Lorentz adapted to MgO. All the other simulations have been carried out using the MNPBEM toolbox [17] using experimental values for the dielectric function of the MgO [24]. Figure 1d has been calculated using the quasistatic eigensolver while Figure 1c, Figure 2, Figure 4 and Figure 5 employ a retarded formulation of the Maxwell equations. Rods have been simulated using approximately 1000 polygons, cubes in vacuum with 5000 polygons and cubes on substrate with 5000 polygons as well. We simulated a 100 nm length cube with approximately 6000 polygons and calculated the corresponding eigencharges and geometrical eigenvalues λ i using the plasmonmode solver.
The radii of curvature of the cube corners in the xy plane are fixed at 3 nm. The rounding in the yz (resp. xz) direction is not precisely controlled within the MNPBEM toolbox [17] (when using the tripolygon and edgeprofile functions). However we estimated the radius of curvature in these planes to be much shorter than 3 nm. Because of the slight asymmetry of the mesh, the three dipoles (resp. quadrupole and edge dipolar) are not slightly degenerated, see λ i values on figure 5.