Probing and steering bulk and surface phonon polaritons in uniaxial materials using fast electrons: hexagonal boron nitride

We theoretically describe how fast electrons couple to polaritonic modes in uniaxial materials by analyzing the electron energy loss (EEL) spectra. We show that in the case of an uniaxial medium with hyperbolic dispersion, bulk and surface modes can be excited by a fast electron traveling through the volume or along an infinite interface between the material and vacuum. Interestingly, and in contrast to the excitations in isotropic materials, bulk modes can be excited by fast electrons traveling outside the uniaxial medium. We demonstrate our findings with the representative uniaxial material hexagonal boron nitride. We show that the excitation of bulk and surface phonon polariton modes is strongly related to the electron velocity and highly dependent on the angle between the electron beam trajectory and the optical axis of the material. Our work provides a systematic study for understanding bulk and surface polaritons excited by a fast electron beam in hyperbolic materials and sets a way to steer and control the propagation of the polaritonic waves by changing the electron velocity and its direction.


I. INTRODUCTION
Polar materials have become of high interest in the field of nanophotonics due to their ability to support phonon polaritons, quasi-particles which result from the coupling between electromagnetic waves and crystal lattice vibrations 1,2 with a characteristic wavelength lying in the mid-infrared region. These quasi-particles can enhance the electromagnetic field deep below the diffraction limit with large quality factors compared to infrared plasmons [3][4][5] , making them promising building blocks for infrared nanophotonics applications [6][7][8][9] .
One interesting 2D polar material is hexagonal boron nitride (h-BN) because of its high quality phonon polaritons and the easy preparation of the single atomic layers made by exfoliation [10][11][12][13][14][15] . Besides being widely used in heterostructures 16 , h-BN is emerging by itself as a versatile material offering novel optical and electrooptical functionalities. The crystal layer structure that constitues h-BN, mediated via van der Waals forces, produces an uniaxial optical response of the material. This implies that the dielectric function of h-BN needs to be described by a diagonal tensor ε ↔ with two principal axes 11,12 : ε x = ε y = ε ⊥ and ε z = ε . When Re(ε ) · Re(ε ⊥ ) < 0, phonon polaritons can propagate inside the material and exhibit a hyperbolic dispersion 6,17 , that is, the relationship between the different components of the polariton wavevector k(ω) = (k x , k y , k z ) traces a surface in momentum space which corresponds to hyperboloids. For h-BN, one can find two energy bands (the Reststrahlen bands) where one of the principal components of the dielectric tensor is negative. Each Reststrahlen band is defined by the energy region between the transverse and longitudinal optical phonon energy, TO and LO, respectively (TO ⊥ and LO ⊥ for the upper Reststrahlen band and TO and LO for the lower Reststrahlen band, see Fig. 1). On the other hand, when Re(ε ) · Re(ε ⊥ ) > 0, the isofrequency surfaces traced by the polariton wavevector in momentum space are ellipsoids.  1 depicts ε ⊥ and ε (see appendix A for expressions and parameters of the dielectric tensor components), which represent the in-plane and out-of-plane dielectric components of h-BN, respectively. The energy range in Fig. 1, shaded in red, corresponds to the lower Reststrahlen band (94.2-102.3 meV) where the real part of the out-of-plane permittivity is negative, leading to isofrequency surfaces in the form of two-sheet hyperboloids (inset Type I). The energy region shaded in grey corresponds to the upper Reststrahlen band (168. .1 meV), where the real part of the in-plane permittiv-arXiv:2006.05359v1 [cond-mat.mtrl-sci] 9 Jun 2020 ity is negative and the isofrequency surfaces correspond to one-sheet hyperboloids (inset Type II).
Hyperbolic phonon polaritons excitable on h-BN within the range of 90 200 meV might be a key to many novel photonic technologies relying on the nanoscale confinement of light and its manipulation. As a result, efficient design and utilization of h-BN structures require spectroscopic studies with adequate spatial resolution. This can be provided, for instance, by electron energy loss spectroscopy (EELS) using electrons as localized electromagnetic probes. Recently, instrumental improvements in EELS performed in a scanning transmission microscope (STEM-EELS) allowed to spatially map phonon polaritons 18 and hyperbolic phonon polaritons in h-BN 19,20 . The focused electron beam of an electron microscope has thus become a suitable probe to access the spectral information of low-energy excitations in technologically relevant materials, with nanoscale spatial resolution. Thus, EEL spectra in phononic materials can be of paramount importance to reveal the properties of phonon polariton excitations.
In this work we first show that a fast electron traveling through bulk h-BN can excite volume phonon polaritons inside and outside the h-BN Reststrahlen bands. Our analysis reveals that the excitation of the volume polariton modes is strongly dependent on the electron velocity and also on the orientation between the electron beam trajectory and the h-BN optical axis. We then study the formation of wake patterns in the field distribution induced by the electron beam at h-BN. Our methodology allows us to connect the excitation of these wake fields with the different electron energy loss mechanisms experienced by the fast electron in the medium: (i) excitation of phonon polaritons or (ii) Cherenkov radiation. We also discuss the emergence of asymmetric wake patterns exhibited by the induced electromagnetic field when the electron beam trajectory sustains an angle relative to the h-BN optical axis. Finally, in the last two section of the paper we show that a fast electron beam interacting with a semi-infinite h-BN interface excite Dyakonov surface phonon polartions within the h-BN upper Reststrahlen band. We further demonstrate that the probing electron traveling above the h-BN in aloof trajectories excites volume phonon polaritons (remotely activation). All these findings offer a way to steer and control the propagation of the polaritonic waves and reveal the importance of the anisotropic optical response of the material in the EELS analysis.

A. Bulk modes in h-BN
According to Maxwell's equations in momentumfrequency (k−ω) space, the dispersion relation for a wave propagating in the volume of an anisotropic material can be found from the following relationship: 21 (1) where ↔ G −1 is the inverse of the Green's tensor, k(ω) = (k x , k y , k z ) is the wavevector of the wave, k 0 = ω/c is the magnitude of the wavevector in vacuum, c is the speed of light, det[x ↔ ] stands for the determinant of a matrix, ⊗ is the tensor product and I ↔ is the identity tensor. Particularly, for an uniaxial medium, the dielectric response can be described in tensor form as ε ↔ (ω) = diag[ε ⊥ , ε ⊥ , ε ]. For this case, two solutions (modes) arise from Eq. (1), yielding the dispersion relation for ordinary waves and the dispersion relation for extraordinary waves Equation (2) represents concentric spheres in k -space for a given energy ω (with ε ⊥ > 0), while Eq. (3) represents hyperboloids or ellipsoids in the reciprocal space depending on the sign of the dielectric components ε and ε ⊥ . Thus, we conclude that the isofrequency surfaces of the polariton wavevector k(ω) in momentum space (for an uniaxial medium) represent geometrically spheres, ellipsoids or hyperboloids. For h-BN, the insets in Fig. 1 depict the isofrequency surfaces for each energy region inside and outside the Reststrahlen bands.

B. Electron energy loss probability
Fast electron beams can couple to bulk polaritonic modes sustained in anisotropic media. We can observe this by analyzing the energy losses experienced by the electron when traveling in such media. Electron energy losses, ∆E EELS , can be calculated within classical electrodynamics as the work performed by the induced electromagnetic field, E ind (r; t), on the probing electron [22][23][24][25][26] where the integration is performed along the electron beam trajectory r e (t), e is the elementary charge, and E ind (r; t) is evaluated along r e (t). Notice that we approximate the electron beam as a classical point charge. The high-energy currents used in typical EELS experiments justify this approximation [27][28][29] . If we Fourier transform E ind (r; t) → E ind (r; ω) in Eq. (4), the electron energy losses can be written as where one identifies the electron energy loss (EEL) probability per unit path, Γ(ω), as withv the unit vector in the same direction as the electron velocity v and t e is the time for the electron to travel a distance dL. Hence, to calculate Γ(ω) one needs to know the induced electric field E ind (r; ω). We derive below the expressions of the total electric field E tot (r; ω) and Γ(ω) for an electron beam trajectory parallel to the optical axis of h-BN, as depicted in Fig. 2. It follows from Maxwell's equations that the field produced by the fast electron plus the induced electric field, namely the total electromagnetic field (E tot (r; ω)) is given by with ε 0 the vacuum permittivity and ρ(k; ω) = −2eπδ(ω − k · v) the charge density of the probing electron. The integration in Eq. (7) extends over the whole reciprocal space and the delta function introduced by the charge density assures conservation of energy and momentum. Indeed, one finds that in the non-relativistic limit the energy that the electron with initial velocity v transfers to the medium upon loosing momentum k is with p = mv the initial momentum of the fast electron. By neglecting recoil of the incident electron, from Eq. (8) one arrives to the so-called non-recoil approximation where ω = k · v. Note that the z -component of the wavevector is fixed by k z = ω/v when the electron travels in the z -direction.
To calculate the bulk loss probability Γ bulk (ω) experienced by the fast electron in the anisotropic medium we substitute Eq. (7) into Eq. (6). Notice that a fast electron traveling in vacuum looses no energy, this allows to use E tot (r; ω) instead of E ind (r; ω) in Eq. (6). Using the cylindrical symmetry of the field produced by the fast electron one finds that where is the probability for the electron to transfer a transverse momentum k ⊥ (to the electron trajectory) upon lossing energy ω. We will refer to this quantity as the momentum-resolved loss probability. In Eq. (10) , and φ is the angle between k ⊥ and the k x -axis, with k c ⊥ the maximum perpendicular momentum of the electrons selected by the collection aperture of the EELS spectrometer.
Particularly, when the electron beam trajectory points out in the same direction as the h-BN optical axis (v = vẑ), expressions for E tot (r; ω) and thus for Γ bulk (ω) can be found in a closed form (see appendix B for the analytical formula of the Green's tensor in uniaxial anisotropic media): where g(r; is written in cylindrical coordinates r = (R, z) = (x, y, z), R = x 2 + y 2 , K 0 (x), K 1 (x) are the zero and first order modified Bessel functions of the second kind, sgn stands for the sign function and γ ⊥ = 1/ 1 − v 2 ε ⊥ /c 2 is the Lorentz factor. The momentum-resolved loss probability becomes and, substituting Eq. (13) into Eq. (9), one finds that 4 The non-retarded versions of P bulk (k ⊥ ; ω) and Γ bulk (ω) can be obtained by setting k 0 equal to zero in Eqs. (13) and (14). The spectrum of the momentum-resolved loss probability and the EEL probability provide valuable information which reveals the properties of the modes of the anisotropic material. We thus explore in the following the connection between the dispersion relation of the h-BN excitations in the upper Reststrahlen band with these two quantities.

C. Upper Reststrahlen band
In the following we address the electron energy losses in h-BN and the connection of these losses with the isofrequency surfaces of the material. We first show in Fig. 3a the isofrequency curve of a h-BN phonon polariton for an energy in the upper Reststrahlen band (red curve). We chose 195 meV as a representative value of this band. When a fast electron beam is used to probe these excitations in the medium, the velocity of the electron determines the momentum transfer, as k · v = ω (Eq. 8).
If the electron is traveling along the z -direction, then k z = ω/v (blue horizontal line in Fig. 3a). Following Eq. (3), this also sets the value of the k ⊥ momentum component (k 2 The intersections between k z = ω/v and the isofrequency curves in the upper Reststrahlen band establish a relationship between the energy ω of the hyperbolic phonon polariton and its perpendicular momentum component k ⊥ . In the left panel of Fig. 3b we plot this relationship (blue dashed line) and the momentum-resolved loss probability P bulk (k ⊥ ; ω) (light blue-yellow contour plot) for v = 0.1c. We note that the highest values of P bulk (k ⊥ ; ω) coincide with the blue dashed line and its asymptotic behavior approaches LO ⊥ for large k ⊥ . This demonstrates that electron energy losses in the upper band are due to phonon polariton excitations. We confirm this by integrating P bulk (k ⊥ ; ω) over k ⊥ up to a cutoff momentum k c ⊥ , which yields the EEL probability Γ bulk (ω) (right panel of Fig. 3b). A clear peak can be observed at the longitudinal optical phonon. This energy loss peak is slightly asymmetric with a broader tail inside the Reststrahlen band compared to that outside the band. Importantly, at energies above LO ⊥ no losses are found. This can be understood with the help of the isofrecuency curves in Fig. 3a. For instance, at energy 205 meV (black dashed line, above the upper Reststrahlen band) the ellipse does not intersect the blue horizontal line and therefore there is no excitation above the upper band. For energies below TO ⊥ , the ellipses may intersect or not the blue horizontal line of k z depending on the particular energy. For instance, at an energy of 160 meV (green dashed line, below the upper Reststrahlen band in Fig. 3a) the ellipse does not cut k z = ω/(0.1c) and therefore there is no excitation induced in that case. However, for lower energies the isofrequency surfaces can cut the k z line, and therefore an anisotropic dielectric mode can be excited (tail below 170 meV in Fig. 3b). We learn from this analysis that the excitation of the phonon polariton modes close to the upper Reststrahlen band is highly dependent on the topology (hyperbolic or elliptic) of the isofrequency surfaces.
The dependency of phonon polaritons excitation on the isofrequency surface allows to control the polaritonic modes as we discuss now in Fig. 3c, where we show the real part of the z -component of the total electric field at ω = 195 meV (representing the energy within the hyperbolic dispersion regime in the upper Reststrahlen band), induced by a fast electron with velocity v = 0.1c. A schematic representation of such electron beam trajectory is displayed in the inset of Fig. 3c. We observe two important features: the formation of a wake pattern and an oscillatory behavior of the field in the zdirection. This spatial periodicity is connected with the parallel momentum component ( k z = ω/v) transferred by the electron since the observed wavelength along the z -axis is λ z = 2π/k z . This implies that the wavelength λ z decreases with increasing energy of the phonon polariton. Furthermore, the direction of the wake pattern is governed by the polariton phase velocity (v p parallel to k(ω), black arrow). The outward direction (relative to the electron beam trajectory) of the wavefronts is determined by the sign of the radial component of v p relative to the radial component of the energy flow (given by the Pointing vector S = E × H parallel to the group velocity v g = ∇ k ω 30,32-35 , magenta arrow). We recognize in Fig. 3a that the group and the phase velocities are nearly perpendicular, and their projection onto the radial axis are parallel, leading to a wave propagating away from the electron beam trajectory (positive phase and positive group velocity with respect to the energy propagation direction). It is worth noting that the projection of the group and the phase velocities onto the beam trajectory direction (z -direction) leads to positive positive phase and negative group velocities relative to S z (Fig. 3a).
As pointed out, for each energy ω, the velocity of the fast electron determines (primarily) the polariton wavevector parallel to the beam trajectory, k z , and consequently the perpendicular wavevector k ⊥ (according to Eq. (3)). To emphasize the velocity dependency, we perform the same analysis ( Fig. 3d-f) as in Fig. 3a-c but increasing the electron velocity to v = 0.5c. In Fig. 3d a zoom into the isofrequency curve of Fig. 3a is presented, together with the value k z (horizontal blue line) determined by the electron velocity v = 0.5c. The increase of the electron velocity leads to the excitation of 195 meV polaritons with reduced momentum (determined by the intersection of the blue horizontal line and the red isofrequency curve). By calculating the momentum-resolved loss probability P bulk (k ⊥ ; ω) (left panel of Fig. 3e) and the EEL probability Γ bulk (ω) (right panel of Fig 3e) we . The black arrows represent the polariton wavevector k(ω), θ k is the angle between k(ω) and the kz-axis, the magenta arrows represent the group velocity vg and the orange arrows the Poynting vector S. Panel (d) shows a zoom into (a). In (d) the horizontal blue line represents kz = ω/(0.5c) . The contour plot (left panel) in (b) shows the momentum-resolved loss probability P bulk (k ⊥ ; ω) normalized to the maximum value (3 a.u.) for v = 0.1c. The right panel in (b) shows the energy loss probability Γ bulk (ω) obtained by integrating P bulk (k ⊥ ; ω) over k ⊥ up to k c ⊥ = 0. find the same behavior as in Fig. 3b for v = 0.1c, except for a one order of magnitude reduction in both k ⊥ and the value of the loss probability. The differences in the properties of the phonon polaritons launched by the fast electron at both electron velocities are distinguishable in Fig. 3f, where we show the real part of the z -component of the total electric field induced by the fast electron with v = 0.5c at energy 195 meV. The spatial period λ z of the polariton is longer compared to that in Fig. 3c as a result of the increase in the electron velocity (smaller k z transferred). Interestingly, the di-rection of the wake field is quite similar to that of panel 3c. This behavior is a specific feature of hyperbolic polaritons since the intersection of the blue line both for v = 0.1c and v = 0.5c occur at the asymptote of the hyperbola which results in polariton wavevectors that have very similar propagation direction but different absolute values.

D. Lower Reststrahlen band
In Fig. 4a we show the isofrequency curve of h-BN phonon polaritons for an energy in the lower Reststrahlen band (red line). Note that the hyperbolas are rotated by 90 • as compared to the upper Reststrahlen band (see Fig. 3a and 3d). However, the momentum k z transferred by the fast electron to the phonon polaritons is still given by the crossing of the hyperbolas with the horizontal blue line (representing k z = ω/v for v = 0.1c in Fig. 4a). From Eq. (3) we obtain the polariton perpendicular momentum k ⊥ , which is shown in Fig. 4b as a function of energy ω (dashed blue curve). We also plot the momentum-resolved loss probability P bulk (k ⊥ ; ω) for energies within the lower Reststrahlen band. Notice that the highest values of P bulk (k ⊥ ; ω) (red and yellow colors in the contour plot) coincide perfectly with the blue dashed curve, demonstrating that the electron energy losses in the lower band are also governed by polariton excitations. However, in contrast to the upper band, we find that the dashed blue curve has a negative slope, dω/dk ⊥ < 0, indicating that the group and the phase velocities are antiparallel (have opposite sign) along the radial direction. We will show below with Fig. 4c that the phase velocity in the radial direction is indeed antiparallel (negative) relative to the Poynting vector (energy flow) while the group velocity in the radial direction is parallel (positive), which is a consequence of the phase and group velocity vectors being perpendicular to each other and rotated by 90 • degrees compared to the upper Reststrahlen band.
To obtain spectroscopic information on the excitations in the lower Reststrahlen band, we calculate the EEL probability Γ bulk (ω) by integration of P bulk (k ⊥ ; ω) in momentum space (right panel in Fig. 4(b)). Contrary to the upper Reststrahlen band, we observe a uniform and relatively small loss probability between TO and LO without the appearance of a sharp peak around LO . We explain this finding by: 1) the large cutoff momenta ( k c ⊥ ) imposed by the aperture of the microscope detector and 2) the relationship between the energy and the transverse momentum of the polaritons in the lower band (see Fig. 4b left panel). Indeed, we observe in Fig. 4b that the asymptotic behavior of the blue dashed line tends to TO for large k ⊥ . This shows that low energy hyperbolic phonon polaritons, close to TO , largely contribute to the energy losses for large k c ⊥ values. Contrary to the upper band, where the high momenta contribution to the electron energy losses comes from polaritons with high energy, close to LO ⊥ (Fig. 3b, left panel). We address the reader to appendix C where we show Γ bulk (ω) in the lower Reststrahlen band for different cutoff values.
The excitation of phonon polaritons (within the lower Reststrahlen band) by the probing electron can be observed in Fig. 4c, where we show the real part of the z -component of the total electric field induced at energy ω = 100 meV. Analogously to the upper band, the oscillatory behavior of the field distribution along the z -direction is governed by the transferred momentum k z . Interestingly, the wake pattern is reversed compared to that for the upper Reststrahlen band (Figs. 3c and 3f), i.e., the wavefronts are propagating toward the electron beam [35][36][37] . By plotting the group and phase velocity vectors onto the field plots (purple and black arrows, respectively; also plotted in Fig. 4a), we can clearly recognize that the projections of both vectors onto the radial axis (perpendicular to the electron beam trajectory) are antiparallel. This leads to a negative phase and positive group velocity relative to the Poynting vector direction (which points always away from the electron beam to preserve causality) along the radial axis. The negative phase velocity in the radial direction is a direct result of the phase velocity vector being nearly pependicular to the Poynting vector, both being rotated by 90 • as compared to the upper Reststrahlen band (where both phase and group velocities are positive relative to energy propagation in the radial direction, see Figs. 3c and 3f).
When the velocity of the electron is increased up to 50 % the speed of light, the k z component of the wavevector parallel to the beam trajectory is reduced. In this case, the matching between the red hyperbola and the horizontal blue line is prevented as observed in Fig. 4d. This mismatch of energy and momentum forbids the excitation of hyperbolic phonon polaritons. However, the blue line intersects the elliptical isofrequency surface of anisotropic bulk phonon polaritons (dielectric) above and below the lower Reststrahlen band (black and green dashed curves calculated for 105 and 92.5 meV, respectively). The matching of energy and momentum at the intersections of the elliptical isofrequency surfaces leads to the excitation of the dielectric modes, as demonstrated by calculating the momentum-resolved loss probability P bulk (k ⊥ ; ω) (left panel of Fig. 4e). This loss probability is determined by the relationship between the energy of the elliptical polartions and the perpendicular momentum component (dashed blue lines, showing ω(k ⊥ ) of the elliptical polaritons). The integration of P bulk (k ⊥ ; ω) in the reciprocal space subsequently yields small energy loss probabilities outside the Reststrahlen band, whereas inside the Reststrahlen band the loss probability is negligible due to absence of polariton excitations.
In Fig. 4f we show the total electric field induced by the electron beam for energies inside (marked 2) and outside (marked 1 and 3) the lower Reststrahlen band. We can observe the formation of wake patterns only for those energies where the dielectric modes are excited (marked as 1 and 3). Importantly, the wake wavefronts propagate outward the beam trajectory as a consequence of the group (v g ) and phase velocities (v p ) being parallel (positive) relative to the Poyting vector in the radial direction (Fig. 4d). We can also notice that the projection of these velocity vectors onto the z -direction is positive. This demonstrates that the radial and z projections of v p and v g for elliptical polaritons are positive, contrary to the hyperbolic regime (Reststrahlen bands) where one of the components is negative (Figs. 3c,f and 4c). x + k 2 y . The horizontal blue line represents the momentum kz = ω/(0.1c) transferred by the fast electron to the polaritons when it travels along the zdirection with velocity 0.1c. The blue line is evaluated at energy 100 meV. The black arrows represent the polariton wavevector k(ω), θ k is the angle between k(ω) and the kz-axis, the magenta arrows represent the group velocity vg and the orange arrows the Poynting vector S. Panel (d) shows a zoom into (a). In (d) the horizontal blue line represents kz = ω/(0.5c) evaluated at energy 92.5 meV. The contour plot (left panel) in (b) shows the momentum-resolved loss probability P bulk (k ⊥ ; ω) normalized to the maximum value (3 a.u) for v = 0.1c. The right panel in (b) shows the energy loss probability Γ bulk (ω) obtained by integrating P bulk (k ⊥ ; ω) over k ⊥ up to k c ⊥ = 0.

E. Induced wake patterns and Cherenkov radiation
We have shown in sections IIC (Figs. 3c,f) and IID (Figs. 4c,f,) that the field distributions produced by a fast electron traveling through h-BN can exhibit wake patterns. The excitation of these patterns (for energies inside and outside the Reststrahlen bands) is connected to the different mechanisms of energy losses experienced by the fast electron in the h-BN. In the following we discuss this connection.
First, it is worth noting that the excitation of the wake fields inside the Reststrahlen bands occurs for energies where electron losses appear (compare Fig. 4c with the image in Fig. 4f labeled as 2). As we pointed out, the electron energy losses within the Reststrahlen bands correspond to the excitation of hyperbolic phonon polaritons. This implies that the wake fields are associated to the excitation of coherent-charge density fluctuations [37][38][39][40][41][42] in the h-BN, namely, the phonon polaritons.
In contrast to the wake fields inside the Reststrahlen bands, the emergence of the wake patterns outside the bands (see Fig. 4, images labeled as 1 and 3) occurs due to a different physical process to that of the excitation of hyperbolic phonon polaritons. Outside the Reststrahlen bands the h-BN dielectric function is purely dielectric and thus the electron energy losses correspond to the radiation emitted by the electron when it passes through the medium with velocity larger than the speed of light (in the h-BN). This mechanism is known as Vavilov-Cherenkov radiation [43][44][45][46][47][48][49] . We have confirmed that the losses in this energy range are present even in the absence of damping in the material (not shown), confirming that the losses are due to Cherenkov radiation in this case. This only happens for velocities which fulfill being consistent with the condition for excitation of Cherenkov radiation 50 . Finally, one can also note that the excitation of the wake fields in the lower Reststrahlen band depends on the electron velocity (compare Figs. 4c and 4f label 2). Indeed, for energies in the lower band one can deduce from Eq. (11) that the wake patterns appear under the following condition: where only the real part of the dielectric function is considered. Interestingly, one can observe that the velocity of the fast electron fulfills different conditions for different energy ranges (compare Eqs. (15) and (16)). This difference is a direct consequence of the distinct physical processes in the excitation of the wake fields. The different nature of the excitation of the wake fields outside and inside the Reststrahlen bands is also reflected in the angle θ w = 90 • − θ k that the wake patterns sustain with respect to the electron beam trajectory. An analysis of this angle and its relationship with Eqs. (15) and (16) is developed at the end of appendix D.

F. Asymmetric wake patterns induced by tilting the electron beam trajectory
As we pointed out in the last sections, the excitation of hyperbolic phonon polaritons can be controlled by the velocity of the fast electrons. In the following we study how steering of phonon polaritons can be controlled via the angle α between the electron beam trajectory and the h-BN optical axis.
When the electron travels at an angle α relative to the h-BN optical axis (illustrated in Fig. 5), the condition for the conservation of energy and momentum given by by the electron to the phonon polaritons (along the beam trajectory given byv) is still given by kv = ω/v. The polariton wavevector can be obtained from the intersection between the blue plane k · v = ω and the isofrequency surfaces (red hyperboloids in Figs. 6a and 6d). Interestingly, we observe that the intersections are not cylindrically symmetric with respect to the z -axis (Figs. 6a,c). This implies that the polaritonic wave will propagate asymmetrically with respect to the electron beam trajectory. Indeed, depending on the direction of propagation, the intersection between the blue planes and the red hyperboloids in Figs. 6a,d will occur at wavevectors k (1) and k (2) whose z -component can be the same (symmetrical case) or different (asymmetrical case). To better understand the different asymmetries in the propagation of the polaritonic wave we refer the reader to appendix D, where we show the intersection of the blue plane and the red hyperboloids (Figs. 6a,d) for selected directions of the wavevector. Notice that the symmetric case is similar to the one we discussed in sections IIC and IID. Therefore we will focus here on the analysis of the polariton propagation direction which shows the largest asymmetry, that is, the k y k z -plane.
We show in Fig. 6a the plane k · v = ω for v = 0.1c (blue surface) and the isofrequency hyperboloid (red surface) for a representative energy in the upper Reststrahlen band ( ω = 180 meV). In Fig. 6b, we plot the projection of the intersection between the blue plane and the red hyperboloid in the k y k z -plane. The blue dashed line represents the electron beam trajectory and the black dashed line the k z -axis. One can notice that the matching between the blue solid line and the red hyperbola (Fig. 6b) occurs at wavevectors k (1) and k (2) whose z -component is different. Thus, the projections onto the z -axis of the phase velocities v (1) p and v (2) p (parallel to k (1) and k (2) , black arrows) are also different. Due to the hyperbolic shape of the isofrequency curve the z -component of the group velocities v (1) g (parallel to the Poynting vector S (1) , right orange arrow in Fig. 6b) and v (2) g (parallel to the Poynting vector S (2) , left orange arrow in Fig. 6b) are also asymmetric. This difference (asymmetry) in the components of the two phase and group velocities leads to a highly asymmetric propaga- tion of the polaritonic wave with respect to the electron beam trajectory.
The dependency of the polaritonic waves on the angle α can be observed in Fig. 6c, where we plot the real part of the z -component of the total electric field produced by the fast electron at energy ω = 180 meV and v = 0.1c when α = 20 • and α = 45 • . Similar to the parallel trajectory (sections IIC and IID), one can notice the formation of wake patterns and the spatial periodicity of the field. This periodicity is determined by the momentum transferred along the beam trajectory ( kv = ω/v) since the corresponding wavelength is λv = 2π/kv. Thus, the higher the energy of the polariton (the velocity of the electron) is, the smaller (bigger) λv will be. The wake patterns formed by the field distribution are clearly asymmetric with respect to the beam trajectory. We observe (Fig. 6c) that the wake fields exhibit largest asymmetry as α is increased from 25 • (Fig.  6c, left panel) to 45 • (Fig. 6c, right panel). This is a direct consequence on how the electron transfers different momentum components, k y and k z , to the polaritonic excitation (see Figs. 6b). One can notice from Fig.  6b that as α is increased, k (1) z ≈ 0 and k (2) tends to the asymptote of the red hyperbola. Therefore, for large angles α the polaritonic wave will propagate relative to the beam trajectory with a phase velocity close to zero on one side of the beam trajectory and with a constant phase velocity on the other side of the beam trajectory.
These findings explain the absence of the wavefronts in Fig. 6c for α = 45 • at the left side of the electron beam. It is worth noting that Fig. 6c corresponds to the propagation of the polaritonic wave in the yz-plane. However, for other propagation directions, the field distributions will be different.
In Figs. 6d-f we show the same analysis (electron beam trajectory tilted an angle α with respect to the h-BN optical axis) for a representative energy within the lower Reststrahlen band (100 meV). Importantly, for this case the projections onto the y-axis of the phase velocities v (1) p and v (2) p are antiparallel (negative) to the y-component of the Poynting vectors S (1) and S (2) . This yields an asymmetric wave propagating with negative phase velocity (Fig. 6f).
Additionally, the electron velocity v allows to control the momentum transfer by the fast electron to the phonon polaritons (Eq. (8)). Indeed, one can obtain the relationship between v and the excitation of the asymmetric wake patterns by analyzing the wake angles θ (Figs. 6c,f). We refer the reader to appendix D where we derived this relationship. For completeness, we show in appendix E the electron energy losses experienced by a fast electron traveling through h-BN in tilted trajectories with respect to the h-BN optical axis.
We have found that the excitation of the polaritonic wave is highly dependent on the orientation of the electron beam trajectory with respect to the h-BN crystallographic arrangement. Thus, while the speed of the electron serves as a means to excite the polaritonic wave or not, the orientation of the electron beam trajectory can serve to control the direction of the polaritonic excitation.

III. EXCITATION OF DYAKONOV SURFACE PHONON POLARITONS IN H-BN BY A LOCALIZED BEAM OF FAST ELECTRONS
We next study the EELS signal when the electron beam is traveling above an h-BN semi-infinite surface. The interface between vacuum and h-BN lies on the yzplane, as depicted in Fig. 7, with the y-axis in the direction of ε ⊥ and the z -axis in the direction of the h-BN optical axis. The electron travels in vacuum at a distance x 0 from the surface (we will refer to this distance as the impact parameter) with velocity v parallel to the optical axis of the h-BN. A schematic representation of the probing electron-surface system is shown in Fig. 7.
Surfaces of uniaxial materials with optical axis parallel to the surface support a specific kind of surface waves, the so-called Dyakonov waves 52,53 . When either ε ⊥ or ε is negative (as in the case of the Reststrahlen bands in h-BN), surface polaritons called Dyankonov surface polaritons 54 can propagate along the surface. Recently, Dyakonov surface phonon polaritons have been observed by scattering-type scanning near-field optical microscopy (s-SNOM) at the edges of h-BN flakes 13,55 as well as by STEM-EELS 19,56 . In the latter experiments, the probing electrons were passing outside the flake edge in a perpendicular trajectory. However, the excitation and detection of Dyakonov surface phonon polaritons with an electron beam parallel to an extended surface has not been described yet.
In the following, we first describe the Dyakonov surface phonon polariton modes that exist at the interface between h-BN and vacuum. We then show that a localized beam of fast electron can couple to these polaritons and we analyse the corresponding EEL spectra and their polariton wake patterns. Importantly, we find that surface Dyakonov phonon polaritons are excited only in the upper Reststrahlen band. Therefore, our analysis and calculations are restricted to this energy range.

A. Surfaces modes in h-BN
According to Dyakonov's theory 52 , the interface described in Fig. 7 supports electromagnetic waves that propagate along it and their associated electromagnetic fields decay exponentially perpendicular to the interface 52,53,57,58 . These surface waves can be expressed as a linear superposition of the four following modes propagating along the interface: (i) a transverse electric (TE) mode, (ii) a transverse magnetic (TM) mode [the corresponding fields decay into the vacuum, upper half space in Fig. 7 labeled I], (iii) an ordinary mode, and (iv) an extraordinary mode [the corresponding fields decay exponentially into h-BN, lower half space in Fig.  7 labeled II]. Following this scheme the electric field in each media can be written as, where harmonic dependency in time has been assumed, and A TE , A TM , A o , A e are the amplitudes of each mode.
The wavevector of each mode is given by where κ I , κ o II , κ e II > 0 and k y , k z ∈ C need to fulfill the following conditions Applying boundary conditions imposed by Maxwell's equations at the interface between vacuum and h-BN, one obtains the following relationship 52,58,59 It is worth noting that Dyakonov's original work 52 was derived for positive ε ⊥ and ε . However Eq. (20) is still valid when ε ⊥ < 0 and ε > 0 57,59 . Since negative values in the real part of the dielectric components support the excitation of polaritonic states, Dyakonov surface waves sustained in h-BN in the mid-infrared region are thus called Dyakonov surface phonon polaritons.
In Fig. 8 we plot the isofrequency contour (red curve) of the h-BN surface polariton for an energy within the upper Reststrahlen band (193 meV), obtained from Eqs. (19a)-(19c) and (20). For comparison, we show a cut (k y k z -plane) of the isofrequency surface of the hyperbolic volume polariton (black dashed line obtained from Eq. (3)). We find that the isofrequency curve of the surface polariton is a hyperbola, similar to that of the volume polariton particularly for small momenta. At large momenta, on the other hand, the opening angle of the isofrequency contour of the surface polariton, θ s , is smaller than that of the volume polariton θ v , demonstrating that the dispersion of Dyakonov phonon polaritons is different to the one obtained for the bulk hyperbolic phonon polaritons.

B. Electron energy loss probability
The excitation of Dyakonov surface phonon polaritons by fast electron beams can be revealed by electron energy loss spectra. In the following we describe the strategy to obtain the momentum-resolved loss probability, P surf (k y ; ω), and the EEL probability, Γ surf (ω), when the probing electron travels above the h-BN surface (see Fig.  7).
To calculate Γ surf (ω), following Eq. (6), one needs to obtain the induced electric field, E ind (r e ; ω), along the electron beam trajectory. To that end we obtain E ind (r; ω) by solving Maxwell's equations in the presence of vacuum-h-BN interface, assuming that the electron travels in vacuum with constant velocity v and impact parameter x 0 (Fig. 7). Considering the boundary conditions at the interfaces (x = 0), one finds the induced electric field in vacuum (region I in Fig. 7): with b I , d I , g I being the coefficients involving the dielectric functions at both sides of the interface andρ = −2πeδ(ω − k z v)e −κIx0 /ε 0 . We refer to appendices F and G for a detailed derivation of the total and induced electric fields at each half space (vacuum and h-BN). By Fourier transforming E ind I (x, k y , k z ; ω) → E ind I (r; ω) in Eq. (21) and inserting its value into Eq. (6), we find that Γ surf (ω) can be written as where is the probability that the electron transfers a transverse momentum k y (y-component of the momentum) upon loosing energy ω. Notice that the z -component of the wavevector in P surf (k y ; ω) is fixed by k z = ω/v, implying that the electron still transfers a parallel momentum equal to ω/v. The integration of Eq. (22) is performed up to the cutoff value k c y , which is determined by the aperture of the EELS detector.
As we discussed in section II, the spectrum of the momentum-resolved loss probability and the EEL probability provides information on the properties of the excited modes in the anisotropic medium. We thus show in the following the relationship between these two quantities and the excitation of Dyakonov surface phonon polaritons.

C. Excitation of surface phonon polaritons
As pointed out above, the parallel momentum k z transferred by the fast electron to the phonon polaritons is determined by the relation k z = ω/v. Similarly to the bulk analysis of section II, this relationship represents a horizontal line in the k y k z representation of Fig. 8. Thus, the transferred momentum can be determined by the crossing between this horizontal line (k z = ω/v) and the isofrequency hyperbolas obtained from Eqs. (19a)-(19c) and (20). From the latter equations one can obtain the relationship between the y-component of the polariton wavevector, k y , and ω, which is shown in the left panel of Fig. 9(a) (dashed blue curve). We also plot the momentum-resolved loss probability P surf (k y ; ω) for energies around the upper Reststrahlen band. The probing electron is traveling above the h-BN surface with an impact parameter of 10 nm and a velocity v = 0.1c. Some similarities between P surf (k y ; ω) and P bulk (k ⊥ ; ω) ( Fig.  3b, left panel) become apparent. For instance, the highest values of P surf (k y ; ω) (red and yellow colors in Fig.  9a) coincide perfectly with the blue dashed curve, demonstrating that the electron energy losses in the upper band are caused mainly due to the excitation of hyperbolic phonon polaritons. However, by comparing Figs. 3b and 9a we recognize that the asymptotic behavior (at large momenta) of P surf (k y ; ω) occurs at a lower energy compared to the asymptotic behavior of P bulk (k ⊥ ; ω). While P bulk (k ⊥ ; ω) tends to the LO ⊥ phonon energy, P surf (k y ; ω) tends to the surface optical (SO ⊥ ) phonon energy given by the condition ε ⊥ (ω SO ⊥ ) = −1 (derived from Eqs. (19a)-(19c) and (20) for large momenta). Importantly, the latter is a fingerprint of the excitation of surface polariton modes. In our case (electron traveling in vacuum above the h-BN surface) these surface modes correspond to Dyakonov surface phonon polaritons. We confirm this by integrating P surf (k y ; ω) over k y up to a cutoff momentum k c y , which yields the EEL probability Γ surf (ω) (right panel of Fig. 9b). A clear peak can be observed at the SO ⊥ phonon energy. This loss peak is slightly asymmetric with a broader tail for lower energies in the Reststrahlen band compared to that for larger energies in the band. Notice that for energies above SO ⊥ the loss spectrum displays a shoulder arising from background losses present in the entire upper band at small momentum (red blurred area for small momentum in the left panel of Fig. 9a).
The excitation of Dyakonov surface phonon polaritons (within the upper Reststrahlen band) by the probing electron can be observed in real space in Fig. 9b, where we show the real part of the z -component of the induced electric field at energies 193 meV (marked as 1) and 198 meV (marked as 2). The top panels correspond to the evaluation of Re(E ind z (r; ω)) in the yz-plane (inplane at the interface), and the bottom panels to the evaluation in the xz-plane (lateral view, containing the electron trajectory). One can recognize from the in-plane views (Figs. 9b marked as 1) the formation of wake patterns and the oscillatory behavior of the induced field in the z -direction. Similarly to the field distribution shown in Fig. 3c, the spatial periodicity along the zdirection is connected with the parallel wavevector component k z = ω/v, since λ z = 2π/k z . Moreover, the wake wavefronts show interesting propagation patterns both in the transverse direction from the beam trajectory as well as into the h-BN.
In the top panel of Fig. 9b (image labeled as 1), the wavefronts along the y-direction propagate with positive phase and group velocities relative to the Poynting vector. Indeed, we find that the dashed blue curve in Fig.  9a has a positive slope (dω/dk y > 0), indicating that the projections onto the y-axis of the group and phase velocities are parallel (positive). We also notice that Dyakonov surface phonon polaritons are confined to the interface with penetration of the field into the h-BN interface (Fig.  9b, bottom image labeled as 1). For energies larger than that of the SO ⊥ phonon, Dyakonov surface phonon polaritons are not excited (Fig. 9b, image labeled as 2). Thus, the induced field distributions for those energies correspond to the reflection of the electron electromagnetic field at the h-BN surface (Fig. 9b, top image labeled as 2). We can also notice that the field penetrates into the h-BN (bottom panel 2 of Fig. 9b), which is connected with the presence of the red blurred region corresponding to the losses appearing for lower momenta in Fig. 9a  (left panel).
When the velocity of the probing electron is increased up to 50 % the speed of light, the momentum parallel to the beam trajectory, k z , is reduced and so does the k y component of the Dyakonov surface phonon polariton. By calculating the momentum-resolved loss probability P surf (k y ; ω) (left panel of Fig. 9c) and the EEL probability Γ surf (ω) one obtains a similar behavior as in Fig.  9a for v = 0.1c, except for a one order of magnitude reduction of both k y and the value of the loss probability.
The differences in the properties of the Dyakonov surface phonon polaritons launched by the fast electron beam can be observed in Fig. 9d, where we show the real part of the z -component of the induced electric field for v = 0.5c at energies 193 meV and 200 meV. Notice that the spatial periodicity λ z of the polariton is longer in this situation compared to that in Fig. 9b as a result of the increased electron velocity. Also, the penetration of the field into the h-BN medium is larger compared to that in Fig. 9b. This increase in the penetration depth can be attributed to the increase of the background losses present in the entire upper band (blurred red are in the left panel of Fig. 9c). For completeness and similar to the analysis presented above, we study in appendix H the excitation of phonon polaritons in the lower Reststrahlen band by a fast electron traveling in aloof trajectory. In the appendix we show the momentum-resolved loss probability (P surf (k y ; ω)), the EEL probability (Γ surf (ω)) and the wake patterns for this energy range.

IV. REMOTE EXCITATION OF BULK PHONON POLARITONS
We have shown in Figs. 9b and 9d that the electric field penetrates into the bulk of the h-BN semi-infinite surface, which is surprising, as one does not expect the excitation of volume modes in isotropic materials for electron beam trajectories outside the material. By comparing the angles of the wake patterns, we demonstrate that indeed volume modes are excited in h-BN with external beam trajectories.
We first calculated the angle θ w of the wake wavefronts produced by the fast electron traveling through bulk h- BN with v = 0.5c at ω = 193 meV (Figs. 10a and 10c), obtaining a value of θ w = 32.35 • . We compare θ w with the angles of the wake wavefronts produced by the fast electron traveling in an aloof trajectory 10 nm above the h-BN surface (Fig. 10b and 10d). From this comparison we find that: (i) the angle θ ws = 24.67 • of the wake pattern at the h-BN surface (Fig. 10b) is different from θ w , and (ii) the angle of the wake pattern excited inside the h-BN is the same as θ w (Fig. 10d). This implies that volume modes are excited by the fast electron traveling in trajectories outside the anisotropic medium.
Importantly, these findings open the possibility of remotely exciting volume phonon polaritons. In contrast to isotropic materials, where an aloof electron beam only couples to surface modes, for anisotropic materials the energy and momentum matching between the electron and the polaritons allows for launching of bulk excitations.

V. SUMMARY
We have thoroughly analyzed the excitation of optical phonon polaritons in hexagonal boron nitride by focused electron beams for two relevant situations: when the electron travels through the h-BN bulk and when it travels in vacuum above a semi-infinite h-BN surface. For the bulk situation, we have observed that the electron couples to volume phonon polaritons. We demonstrated that the excitation of these polaritonic modes is strongly dependent on the electron velocity and on the angle between the optical axis of h-BN and the trajectory of the electron beam. Furthermore, we have shown that Dyakonov surface phonon polaritons can be excited by a fast electron traveling above the h-BN surface. Interestingly, aloof electron beams are capable of exciting volume polaritons in the h-BN.
By a detailed mode analysis, we showed that the electron beam transfers a specific momentum to the modes. This momentum transfer determines the properties of the excited phonon polaritons, and thus controls their phase and group velocities, as well as their propagation direction. Importantly, we found that the propagation of the polaritonic waves is highly asymmetric with respect to the electron beam trajectory when the trajectory sustains an angle relative to the h-BN optical axis.
Our findings may offer a way to steer and control the propagation of the polaritonic waves excited in hyperbolic materials. Although we studied the specific material h-BN, our findings can be generalized and can serve as an guide for the correct interpretation of the different excited modes and loss channels encountered in EELS experiments of uniaxial materials. The two components of the h-BN dielectric function can be described by a Drude-Lorentz model as 11 with ω LO , ω TO the phonon LO, TO energies, respectively, ε ∞ is the high-frequency dielectric permittivity and γ is the damping constant. The values used for each constant are presented in Table I. In this work we use the following Fourier transform convention where F(r; t) is a smooth vector field representing the electric or magnetic fields and V stands for the volume in the euclidean space R 3 . Thus, the Green's tensor satisfying the wave equation 60-63 can be expressed in k − ω space as follows From Eq. (B3) one can deduce that the inverse of the Green's tensor for an uniaxial medium can be decomposed in the form ε y and ε = ε z . This tensor decomposition allows for finding the following closed expression for ↔ G(k; ω) 21 :   Fig. 11 we show the EEL probability (Γ bulk (ω), given by Eq. (14)) in the vicinity of the lower Reststrahlen band for different cutoff values k c ⊥ : (a) 1 × 10 −2Å −1 , (b) 1 × 10 −3Å −1 , (c) 1 × 10 −4Å −1 and (d) For the calculation of Γ bulk (ω) we consider v = 0.1c.
One can observe that for small cutoff momentum the EEL probability of the LO phonon energy is better defined. Whereas for large cutoff momenta the sharp peak in Fig. 11d broadens. However, cutoff values of 1 × 10 −4Å −4 or ×10 −4Å −5 are not experimentally feasible. When the electron beam trajectory makes an angle α relative to the h-BN optical axis, the propagation of the phonon polaritons (excited by the fast electron) is highly asymmetric with respect to the beam trajectory. We analyze these asymmetries in the following.
The propagation of the polaritonic wave is governed by its phase velocity and thus, by the polariton wavevector k(ω) = (k x , k y , k z ) which fulfills Eq. (3). When the hyperbolic phonon polaritons are excited by an electron beam, the components of k(ω) have also to fulfill Eq. (8), that is, the components of k(ω) can be obtained from the following two expressions where we assume that the electron velocity is v = v(0, sin α, cos α). Moreover, if we decompose k(ω) in cylindrical coordinates as k(ω) = (q cos φ, q sin φ, k z ), with φ the azimuthal angle of the symmetry axis, and substitute it into Eqs. (D1a) and (D1b), we obtain to the following system of equations for q, φ and k z . Notice that the variable q corresponds to k ⊥ for trajectories parallel to the h-BN optical axis. However, for the oblique trajectory q = ( k x , k y ) is no longer orthogonal to the beam trajectory and thus we avoid referring to it as the transverse momentum. One can deduce from Eqs. (D2a) and (D2b) that the solutions have cylindrical symmetry (symmetric with respect to the z -axis) when α = 0 • . For cases where α = 0 • , this symmetry is broken and the solutions depend on the azimuthal angle φ. We explore this dependency below. In Fig. 12 we show the intersection between the h-BN isofrequency hyperboloids (red surfaces, Figs. 12a,c) and the plane k · v = ω determined by the electron beam trajectory (blue surfaces, Figs. 12a,c). Notice that the direction of the electron beam trajectory is orthogonal to the blue plane k · v = ω. We analyze an electron beam with velocity v = 0.1c and a trajectory angle of α = 20 • . Finally we chose two representative energies, one in the upper Reststrahlen band at 180 meV (Fig. 12a) and the other one in the lower Reststrahlen band at 100 meV (Fig. 12c). The grey 2D plots in Figs. 12b and 12d show the intersection between the red hyperboloid and the blue plane along four different directions determined by the azimuthal angle φ: 0 • , 60 • , 90 • and 150 • . In the 2D projections the blue dashed lines depict the beam trajectory, as viewed from the direction determined by φ. The polariton wavevector along each particular direction can be obtained from the intersection between the blue lines and the red hyperbolas. Importantly, one can recognize from the 2D projections that: 1. The intersection between the blue line and the red hyperbola is asymmetric with respect to the z -axis for φ ∈ (0 • , 180 • ), as we observe in Figs. 12(b) and 12(d) for φ = 60 • , 90 • , 150 • .
2. The direction of largest asymmetry occurs at φ = 90 • (k y k z -plane) and the direction of symmetric propagation occurs at φ = 0 • (k x k z -plane).
3. The intersections between the blue lines and the red hyperbolas are also asymmetric (or symmetric) with respect to the electron beam trajectory (blue dashed line).
To better understand the asymmetries in the propagation of the polaritonic waves, we focus on the direction of largest asymmetry: φ = 90 • (equivalently the k y k zplane). From Eqs. (D1a) and (D1b) one can obtain the following two solutions for the polariton wavevector in the k y k z -plane From Eqs. (D3a) and (D3b) one can recognize that k z , showing the asymmetry in the propagation of the polaritonic wave. Moreover, the angles θ (1) k and θ (2) k defined by k (1) , k (2) vectors with respect to the electron beam trajectory (see Figs. 6b and 6e) satisfy the following relations tan(θ When α = 0 • , one can deduce from Eqs. (D5a) and (D5b) that tan θ tan θ Therefore θ (2) k = θ k for this particular case of symmetric propagation. Notice that θ k is also preserved in any other azimuthal direction.
We can observe from Eqs. (D3a) and (D3b) that k (1) , k (2) depend on the electron velocity v. This dependency provides information on the condition that the electron velocity needs to satisfy for the electron beam to excite the polaritonic waves. Indeed, by imposing real value solutions to Eqs. (D5a) and (D5b), one obtains the following condition on v: This last relationship results in the following inequality when α = 0 • , which coincides exactly with the first inequality in Eq. (16) obtained in the main text. As we discuss in section IIE, Eq. (D8) reveals the condition on the electron velocity for exciting phonon polaritons or emitting Cherenkov radiation.
We show now that we can recover the properties of the excited wave in an isotropic dielectric medium from the previous expressions. Assuming that the medium has dielectric function equal to ε ⊥ = ε = ε > 0, the condition (D8) results in the canonical relation for Cherenkov radiation: v > c/ √ ε. Moreover, the two wavevector solutions k (1) , k (2) given by Eqs. (D3a) and (D3b) result in with It is worthwhile noting that M is an orthogonal matrix. This implies that the angles θ (1) k and θ (2) k are always equal. Thus, the propagation of the wake patterns excited in an isotropic dielectric media is always cylindrically symmetric with respect to the electron beam trajectory.
Appendix E: Momentum-resolved loss and EEL probabilities for electron trajectories oblique to the optical axis of h-BN As we show in appendix D, the cylindrical symmetry in the propagation of the phonon polariton wave is broken when the electron beam trajectory is not parallel to the h-BN optical axis. This break in symmetry means that the momentum-resolved loss probability P bulk (q; ω) is no longer constant along the azimuthal direction but it depends on the angle φ 64 . Notice also that the momentum q = ( k x , k y ) is no longer perpendicular to the beam trajectory r e (t) = vt(0, sin α, cos α). In fact, the two orthogonal directions to r e (t) are: (i) the xdirection and (ii) the direction set by the unit vector n α = (0, cos α, − sin α). Thus, the two transverse components (to the beam trajectory) of the polariton wavevector are k x and k α = k ·n α = k y cos α − k z sin α. (E1) Furthermore, the components of the polariton wavevector k(ω) excited by the fast electron beam need to satisfy Eq. (D1b). By solving Eqs. (D1b) and (E1) one finds that k y and k z can be written in terms of k α as Following Eq. (10), we can define the probability for the fast electron to transfer a transverse momentum ( k x , k α ) upon loosing energy ω as where ↔ G * = ↔ G(k x , k * y , k * z ) and k * y , k * z are given by Eqs. (E2a) and (E2b), respectively. On the other hand, the electron energy loss probability, Γ bulk α (ω), can be obtained by integrating P bulk α (k x , k α ; ω) over the momentum coordinates (Eq. (9) where the last equality follows by expressing q in cylindrical coordinates. Notice that the integration over the magnitude of q is performed up to the cutoff value q c .
In Fig. 13 we show the momentum-resolved loss probability P bulk α (k x , k α ; ω) and the EEL probability Γ bulk α (ω) for representative energies inside the Reststrahlen bands when v = 0.1c and two different trajectory angles α: 20 • and 45 • . One can observe in the figure that the EEL features are similar but the momentum-resolved loss probability shows asymmetries for different energies in the Reststrahlen bands.
By solving the linear system of equations set by the boundary conditions (Eq. (F4)), one finds that each coefficient in Eqs. (F2a)-(F2f) can be expressed as A II =ρ a II , B I =ρ b I , C II =ρ c II D I =ρ d I , F II =ρf II , G I =ρg I , withρ = −2πeδ(ω −k z v)e −κIx0 /ε 0 . Thus, we obtain that the induced electric fields in vacuum (labeled as I) and h-BN (labeled as II) are given by (Eqs. (F2a)-(F2f)) E ind II (x, k y , k z ; ω) = (a II , c II , 0)ρ e κ o II x (G1b) Substituting Eq. (G1a) into Eq. (6), one obtains that the EEL probability Γ surf (ω) can be written as Γ surf (ω) = e π ω Re E ind I (r e ; ω) ·ẑ e −iωte (G2) = k c y 0 dk y P surf (k y ; ω), with k c y the maximum momentum of the electrons that can pass through the collection aperture of the detector in the y-direction, and P surf (k y ; ω) = − e 2 π 2 ε 0 ωv Re g I e −2κIx0 where k z = ω/v is the momentum transferred by the electron to the polaritons along the beam trajectory.
Appendix H: Electron energy loss probability for energies around the lower Reststrahlen band for electron trajectories above the surface of h-BN In the left panel of Fig. 14a we show the momentumresolved loss probability, P surf (k y ; ω), for an electron traveling above an h-BN surface for energies around the lower Reststrahlen band. The probing electron travels above the surface at an impact parameter of 10 nm and v = 0.1c. The blue dashed line corresponds to the bulk phonon polariton dispersion (Eq. (3)). We can recognize some similarities between P surf (k y ; ω) and P bulk (k ⊥ ; ω) (compare the left panels of Figs. 4b and 14a). For instance, the maximum values of P surf (k y ; ω) are close to the bulk dispersion (blue dashed line). Interestingly, this bulk dispersion corresponds to the envelope curve of P surf (k y ; ω) implying that electron energy losses in the lower band are mainly due to bulk hyperbolic phonon polariton excitations. To obtain spectroscopic information on the excitations in the lower band, we calculate the EEL probability Γ surf (ω) by integrating P surf (k y ; ω) over k y up to a cutoff k c y (right panel in Fig. 14a). Similarly to Γ bulk (ω) (Fig. 4b, right panel), Γ surf (ω) (right panel in Fig. 14a) exhibits a uniform loss probability between TO and LO which depends on the selected cutoff momenta k c y . In Fig. 14b we show the real part of the z -component of the induced electric field for the same electron velocity and impact parameter as in Fig. 14a, for two different energies marked as 1 and 2 in panel a. We can recognize the excitation of the wake fields in the h-BN surface for those energy losses (compare the top panels labeled as 1 and 2 in Fig 14b). The bulk nature of the excited modes is revealed in the bottom panels of Fig. 14b, where we show z -component of the real part of the induced electric field Re(E ind z (r; ω)) in the xz-plane. In this lateral view of the field distribution one notice the excitation and propagation of the field into the bulk from the h-BN surface.
Panels 14c and 14d show P surf (k y ; ω), Γ surf (ω) and the induced field distribution when the electron velocity is 0.5c. It is worth noting that the blue dashed line superimposed on P surf (k y ; ω) (Fig. 14c, left panel) corresponds to another branch of Dyakonov's dispersion relation given by Eqs. (19a)-(19c) and (20).