Theory of Optical Nonlocality in Polar Dielectrics

Sub-wavelength confinement of mid-infrared light can be achieved exploiting the metal-like optical response of polar dielectric crystals in their Reststrahlen spectral region, where they support evanescent modes termed surface phonon polaritons. In the past few years the investigation of phonon polaritons localised in nanoresonators and layered heterostructures has enjoyed remarkable success, highlighting them as a promising platform for mid-infrared nanophotonic applications. Here we prove that the standard local dielectric description of phonon polaritons in nanometric objects fails due to the nonlocal nature of the phonon response and we develop the corresponding nonlocal theory. Application of our general theory to both dielectric nanospheres and thin films demonstrates that polar dielectrics exhibit a rich nonlocal phenomenology, qualitatively different from the one of plasmonic systems, due to the negative dispersion of phononic optical modes.


Introduction
Recent progress in nanophotonics was accelerated by the development of light concentrators able to localise electromagnetic energy on length scales substantially smaller than the diffraction limit [1]. Self-sustained electromagnetic oscillations on deep sub-wavelength scales are made possible by storing part of the electromagnetic energy into kinetic energy of charged particles [2]. At near-infrared and shorter wavelengths the use of electron gases in metals and doped semiconductors led to the spectacular development of plasmonics [3][4][5]. Plasmonic concentrators rely on hybridisation between photons and the free electron gas of the supporting metal, meaning noble metal plasmonics is inefficient in the mid-infrared spectral region where the photon energy is entirely out of resonance with the plasma frequency. A promising alternative in this region are concentrators based on surface phonon polaritons (SPhPs), which achieve light confinement in the mid-infrared region with comparatively limited losses by storing electromagnetic energy as vibrations of the crystal lattice [2,[6][7][8]. Localization of SPhPs in user defined resonators [9,10] allows for extremely small mode volumes [11,12] in highly tuneable [13][14][15][16], low-loss modes, with applications in sensing [17], nonlinear optics [18][19][20], and the creation of nanophotonic circuitry [21]. Typical electromagnetic theories describing light concentrators are parameterised by frequency dependant dielectric functions imposing spectrally and spatially local relationships between electric field and polarization. Such an approximation, which assumes that the length-scales of photon fields and of the microscopic process giving rise to the dielectric response are well separated, fails in plasmonic systems whose size approaches 10nm [22], or which contain geometrical singularities resulting in field hotspots of the same dimension [23]. This leads local response theories to * Corresponding author: s.de-liberato@soton.ac.uk predict erroneous modal frequencies and a significant overestimation of the achievable field enhancement [24][25][26][27]. These effects largely arise from the presence of longitudinal plasma waves in the electron gas, excited at the edge of the scatterer by the transverse photon field. Nonlocal plasmonic theories correct these errors by considering the microscopic dynamics of the electron gas, typically by coupling a hydrodynamic Drude model to Maxwell's equations [28][29][30][31][32]. Physically these models account for pressure in the electron gas, resulting in a smearing of electromagnetic hotspots, plasmon spill-in induced frequency shifts, and the emergence of confined longitudinal modes [33][34][35]. Polar dielectrics support longitudinal optical (LO) phonons, analogous to longitudinal plasma waves in noble metals. In nano-objects these can also be expected to hybridise with the photon field and perturb the system response. As SPhP concentrators improve it is thus necessary to understand the limits of a description based on a local dielectric response. Hybridisation of LO phonons and SPhPs has been suggested as a possible path toward electrical injection of SPhPs, and it was recently demonstrated by using elongated silicon carbide (SiC) polytypes [36], where a weak phonon mode [37] naturally crosses the SPhP dispersion. In this paper we develop a general theory of the nonlocal optical response of polar dielectrics, coupling the ionic equations of motion of the crystal lattice to Maxwell's equations, taking into account the optical phonons' dispersion. We firstly apply such a theory to calculate the extinction spectrum of 3C-SiC nanospheres. We find that confined LO phonon modes, redshifted from the local Fröhlich resonance, appear in the extinction cross section. We also observe size-dependent nonlocal broadening, not present in the hydrodynamic plasmonic model. Not observed is the strong size-dependant frequency shift seen in plasmonic nonlocality, a consequence of the narrow width of the Reststrahlen band. Secondly, we study the nonlocal response of the epsilon-near-zero (ENZ) resonance of a thin polar dielectric film [16,38]. The ENZ mode of these structures is quasi-resonant with the LO phonon frequency, making it extremely sensitive to changes in the longitudinal mode frequency.  1: Illustration of the different impact of nonlocality in plasmonic (a) and phononic (b) systems. On each plot the shading indicates the region where evanescent optical modes are supported. The solid curves shows (on a different wavevector scale) the dispersion of the bulk longitudinal modes and the dashed line shows the dispersion of the surface polariton at a planar interface. In the metallic case the energy of the plasma wave increases with wavevector, leading to the existence of modes with real wavevector for frequencies above the plasma frequency ωP and with imaginary wavevector below. Surface plasmons are supported for frequencies ω < ωP, and therefore cannot couple with propagative plasma waves. In the phononic case the opposite is true. Longitudinal phonons with real wavevectors exist for frequencies below the LO phonon frequency ωL and with imaginary wavevector above. SPhPs supported in the region ω < ωL are thus resonantly coupled with propagative longitudinal waves. The insets in each plot sketch the shape of the longitudinal electric field in the region where surface polaritons are supported.
Our theory demonstrates that the impact of nonlocality on the optical response of phonon polaritons is qualitatively different from the case of nonlocal plasmonics. This can be intuitively understood from Fig. 1a and Fig. 1b, which illustrate the effect of modal dispersion in both plasmonic and phononic systems. Due to the positive dispersion of plasma frequency in noble metals, it is not possible to have resonant coupling between longitudinal propagative modes and the lower-lying surface modes. Optical phonons instead have negative dispersions, meaning a set of propagative longitudinal modes exists in the same frequency region as the SPhP mode, yielding hybrids with mixed longitudinal-transverse nature. Given the by now relatively well-known phenomenology of nonlocal plasmonics, in this paper we will use it as a comparison to highlight the unique features of nonlocality in polar dielectrics

Continuum Dielectric Model
In a generic linear material the electric field E (r, t) and the displacement field D (r, t) are linked by the con-stitutive relation The usually applied local response approximation (LRA) consists in neglecting the nonlocal nature of the dielectric response, leading to the expression in which ǫ LRA (r, ω) is the frequency-and positiondependant local dielectric function. Although such a formalism can be used to study SPhPs in inhomogeneous systems [39], it remains accurate only while we can neglect nonlocal effects due to the phonon dispersion. In order to move beyond the local theory it is necessary to select a specific model to describe the dynamics of the lattice. We treat our polar crystal as an isotropic lattice, described in the continuum limit by [40][41][42][43][44][45][46].
in which X is the relative ionic displacement, ω T is the transverse optical (TO) phonon frequency and ρ and µ are the effective mass and charge densities. Note that here and in the following we suppress space (wavevector) and time (frequency) dependencies where possible. For simplicity in this paper we will consider a single loss rate γ, independent of both frequency and polarization. The spatially dispersive terms in Eq. (3) are calculable as the divergence of a tensorτ in whichĪ is the identity tensor. In a model of acoustic phononsτ is simply the mechanical stress tensor, parameterised by the Lamé coefficients which describe the acoustic phonon velocities. As here we consider optical phonons the coefficients β L , β T are instead phenomenological velocities which describe LO and TO phonon dispersions respectively [47]. Note that this formalism only accounts for terms up to quadratic order in the phonon dispersion, potentially leading to deviation far from either phonon frequency. At this point we introduce the polarisation field through the equation in which ǫ ∞ is the high-frequency dielectric constant. As a consequence of Eq. (3), dispersive LO phonons exist at all the (complex) frequencies satisfying where the LO frequency ω L is defined by the relation Transverse fields in the lattice are analogously described by a spatially dispersive dielectric function calculable from Eqs. 3 and 5 The key point to note from Eqs. 6 and 8 is that an increase in the wavevector k results in a red shift in the phonon frequencies. This differs from plasmonic systems which exhibit a blue shift in the plasma frequency with increasing wavevector where ω P is the zone-centre value and β is a characteristic velocity governing the dispersion of the plasma. Surface plasmons, mediated by the negative dielectric function of the metal, are only supported for ω SP < ω P , a region where Eq. (9) only has evanescent solutions (Fig. 1a).
Longitudinal plasma waves then decay with typical skin depth of the order 1Å. Phonon modes in polar dielectrics exhibit instead a red shift with increasing wavevector (Fig. 1b) meaning LO modes with real wavevector coexist spectrally with SPhP modes.
To illustrate the striking impact of this difference on the nonlocal optical features we show a comparison of extinction spectra calculated utilising nonlocal Mie theory for a 10nm diameter gold (Fig. 2a) and 3C-SiC (Fig. 2b) spheres. The theory underlying this calculation is presented in a later section and the local parameters for 3C-SiC have been found in the literature [37]. For the nonlocal velocities we utilise β T = 9.15 × 10 6 cm/s, β L = 15.39 × 10 6 cm/s, approximated from studies of phonon dispersion in 3C-SiC [48,49]. For the LO phonon damping rate we extrapolate the low temperature results of Debarnardi et al. to room temperature, assuming a 2 phonon decay channel [50]. The nonlocal parameters for gold are taken from the work of Cirací [24]. In the metallic case the plasmon resonance blue shifts slightly as a result of electron spill-in and additional peaks appear for frequencies larger than ω P , where Eq. (9) has real wavevector solutions. Note that the additional cusp above the dipolar resonance in Fig. 2a corresponds to a higher order mode. The peaks above the plasma frequency correspond to confined longitudinal plasma modes in the metallic sphere, excited by the photon field at the particle boundary. These results are analogous to those previously presented for metallic nanowires of similar dimension [22]. In the polar dielectric case instead, only a minor red shift of the SPhP is present but additional discrete peaks are observed below the LO phonon frequency where Eq. (6) has real wavevector solutions. These correspond to confined LO phonon modes in the bulk of the nanoparticle.

Additional Boundary Conditions
In the LRA mode amplitudes at the interface between two piecewise continuous media are determined by application of the Maxwell boundary conditions (MBC) to electric and magnetic fields E, and H. As a nonlocal formalism entails consideration of additional fields in each layer additional boundary conditions (ABC) are necessary to uniquely determine all the mode amplitudes. In order to find the boundary conditions we use the ionic equation of motion Eq. (3) to derive a hybrid Poynting vector describing the total energy transported in both the electromagnetic and mechanical fields [51] where the tensorτ is defined by Eq. (4). In order to ensure that interfaces act as neither energy sources or sinks, we impose continuity of the normal component of S. This requires, in addition to the standard MBC, continuity of X and of the normal components of the tensorτ .
At an interface between nonlocal and local media the boundary conditions are under-specified. This is a problem well known in the case of exciton polaritons, where multiple solutions have been proposed, derived either from microscopic considerations [52,53] or by application of symmetry conditions to the fields [28,54]. The most commonly utilised ABC are the Pekar-Ridley [52,55] and Fuchs-Kliewer [28] which, assuming for definiteness an interface lying in the xy plane read X = 0 and ∂ z X {x,y} = X z = 0 respectively. Both of these choices satisfy ρτẊ ·ẑ = 0, ensuring conservation of energy across the interface when the mechanical excitations cannot propagate (e.g., at a dielectric-vacuum interface).
In the isotropic model we are considering the additional high-momentum transverse mode coupled by the nonlocality is closely resonant with the bare TO phonon, and thus it has a negligible associated electric field. ABC then predominantly mix the low-momentum transverse photon-like mode and the dispersive LO excitation. The condition on X {x,y} is then trivially satisfied without affecting the optical response. This means that calculations of observables such as reflectance utilising the Fuchs-Kliewer and Pekar-Ridley boundary conditions are practically indistinguishable. Fuchs-Kliewer type ABC are derived by symmetry arguments, and they are equivalent to assuming specular reflection from the interface, which is the relevant symmetry in our case. We thus utilise this condition in the following instead than the Pekar-Ridley condition typically applied in plasmonic systems [29].

Scattering from Spherical Resonators
In this Section we exploit an extended Mie theory, detailed in Appendix A, to calculate the scattering and extinction cross sections of nanoscopic polar dielectric spheres. In order to provide a better understanding of the results, a simplified analytic model valid in the quasistatic limit is also formulated [56]. The approach we take is conceptually similar to previous nonlocal extensions of Mie theory in cylindrical polar coordinates [57,58] which have been utilised to describe nonlocal effects in noble metallic systems [22]. In Fig. 3a we plot the extinction cross section σ ext , for 3C-SiC spheres of 5, 10 and 20nm radii. In the local case deep-subwavelength spheres exhibit a single dipolar Fröhlich resonance where ǫ LRA (ω) = −2, which for 3C-SiC is at ≈ 934/cm. In the nonlocal extinction spectra this resonance remains, however in smaller particles additional closely spaced peaks appear for frequencies below that of the zone-centre LO phonon. These correspond to confined longitudinal modes propagating into the bulk of the nanosphere. We can see that as the sphere radius decreases the Fröhlich resonance splits into two distinct peaks. This is more visible in In Fig. 3b where we consider smaller radii and zoom in on the peak region. This phenomenon, absent in the plasmonic case, is due to the existence of resonant longitudinal modes which can therefore become strongly coupled to the Fröhlich resonance. Note that the results for the 1nm radius are only tentative, as our continuum model in Eq. (3) can not be considered accurate when for sizes of the order of the lattice spacing. When strong coupling is not achieved we recover instead a behaviour redolent of dark state hybridisation, in which a broad bright state hybridises with narrow dark states. This results in the typical Fano scattering spectra. Here the Fröhlich resonance which radiates in the far field plays the part of the bright state, while the longitudinal modes are dark. Fano interference acts to reduce the bright states far-field scattering, while increasing the absorption as energy is funnelled into non-radiative dark states [59]. This is visible in the comparison of scattering and absorption efficiencies for a 2nm radius sphere in In order to gain a better understanding of these results it is informative to consider a simplified semi-analytical model. We thus treat the scattering problem in the quasistatic limit, assuming the impinging radiation field does not vary appreciably over the sphere diameter. This is a good approximation for the parameters considered in this paper, with radii R normalised over the resonant wavelength of the order of 0.01. As detailed in Appendix B in this limit the polarisability of a sphere in vacuum has the form whereǫ is the local dielectric function renormalised by nonlocal effects through the nonlocal parameter with j 1 (x) and j ′ 1 (x) a spherical Bessel function of the first kind and its derivative, and ξ the longitudinal mode wavevector defined by the equation Within this formalism we can quickly evaluate the extinction cross section For the parameters presented in Fig. 3a the extinction cross section given in Eq. (15) is indistinguishable from the full Mie theory.

Enhanced Nonlocal Broadening
In a hydrodynamic model of nonlocal plasmonics the electron dynamics are driven by pressure, resulting in convective charge flow. This results in a shift in the real mode frequency of the plasmon mode when fields are sufficiently confined for such flow to play a role. Another nonlocal effect often observed in the optics of small plasmonic particles is a size-dependent damping, initially described by Kreibig [60]. This can either be introduced phenomenologically as part of a hydrodynamic nonlocal treatment of the structure or treated more thoroughly as part of a generalised nonlocal optical response theory by accounting for entropy driven diffusive electron flow [25] or by accounting for the viscosity of the electron gas [61]. This results in an increasing linewidth as any induced charge will over time diffuse away, diminishing the plasmon mode. Any broadening is size-dependent, requiring electrons to diffuse over a comparable length to the field confinement length over the bare mode lifetime. The reason size-dependent damping does not emerge from a purely hydrodynamic model in the absence of diffusion is because, as previously explained, the longitudinal plasma waves excited in the electron gas obey Eq. (9), and they are thus evanescent in nature. It is easy to show that for an evanescent mode at a planar interface, the normal component of the time-average mechanical Poynting vector in Eq. (10) is zero. In the polar system studied in this paper the LO modes excited by the SPhPs are propagative in nature, and therefore they transport energy away from that interface. This allows them to act as an additional loss channel, demonstrating hydrodynamic Kreibig damping in the absence of diffusive effects, which remains important even for relatively large particles. In Fig. 4a we plot the quality factor of the fundamental Fröhlich resonance for 3C-SiC spheres with radii between 10 and 1000nm, obtained from the extinction cross section in Eq. (15). In this region the particle is sufficiently large to support many closely spaced longitudinal modes, whose finite linewidths mean they form a continuum and do not appear as discrete peaks in the extinction cross section. The results show that for particles with radii less than 1µm, coupling to these non-resonant modes leads to a broadening relative to the local result which increases as the particle radius decreases. We can form a model analogous to that of Kreibig to fit the data, assuming a nonlocal damping rate of the form The dimensionless constant A describes the magnitude of the size-dependent damping. Fitting this formula to the numerical data we find A = 0.03. We can utilise this to calculate a typical lengthscale over which Kreibig damping becomes dominant which is comparable to that for diffusive damping in plasmonic systems where the characteristic velocity and bare damping rate are both enhanced by around 2 orders of magnitude [25].

Frequency Shifts and Purcell Factor
In plasmonic systems the dominant nonlocal effect is typically a blue shifting mode frequency [22], manifested in the quasistatic model as an increase in the effective local dielectric functionǫ L (ω). In the case of polar dielectrics the Fröhlich resonance also undergoes a nonlocal shift, toward the red in this case, but of much smaller amplitude. This can be verified from Fig. 5a, in which we compare plasmonic and phononic shifts as a function of the nonlocal parameter δ NL . Such an effect is due to the narrow nature of the Reststrahlen band, which leads to a faster frequency dispersion of the local dielectric function. The system can thus accommodate a change in the dielectric function due to a certain value of the nonlocal parameter δ NL by a much smaller shift in frequency.
We also calculated, using the quasistatic approximation the total emission rate of a dipole 5nm away from the surface of the sphere with averaged orientation as a function of radius [62,63]. The results, normalised to the local case, are shown in Fig. 5b, where for comparison we also show the result for a gold sphere. Note that quasistatic calculation only considers decay into the dipolar mode of the sphere, ignoring quenching due to higher order resonances. In both the metallic and polar dielectric systems a drop in the decay rate is observed, due to smearing of field hotspots due to nonlocality. In the dielectric case though, we see a number of extra resonances with enhanced decay rate appearing, due to the existence of quasi-resonant, hybridised LO phonon modes.

Epsilon Near Zero Modes in Thin SiC Films
We pass now to investigate nonlocal effects in a ultrathin 3C-SiC film in vacuum, which will allow us to gain an intuitive understanding of nonlocal physics in both thin films and semiconducting heterostructures. Optically thick films support degenerate SPhP modes at each interface. In a thin film of thickness d these hybridise yielding two distinct SPhP branches, interpreted as symmetric and antisymmetric superpositions of the modes on each interface. The symmetric mode is blue-shifted from the bulk mode and in the d → 0 limit exists at the asymptotic zone-centre longitudinal optical phonon frequency. As the dielectric function of the film vanishes here this mode is often referred to as an ENZ mode [16,38].
A sufficiently thin film also acts as a Fabry-Perot for longitudinal phonons, supporting a series of discrete modes with quantised out-of-plane wavevector, satisfying k z = nπ d , with n a positive integer. For wavevectors in the film plane k ≪ nβLπ d the LO phonon modes have frequencies a result analogous to that observed for ENZ modes in thin conductive films, where additional resonances obeying a similar equation to were observed above the film plasma frequency [64]. In a sufficiently thin sample the fundamental LO resonance ω 1 can differ appreciably from the zone-center frequency ω L , leading to a shift in the upper edge of the Reststrahlen band. As the ENZ mode is quasi-resonant with the LO phonon, it can be expected to be extremely sensitive to such a shift.
To study this we calculate the frequency of the ENZ mode at constant in-plane wavevector k c = 1200cm −1 for a 3C-SiC film as a function of inverse film thickness in the local and nonlocal cases, the accompanying theory is presented in Appendix C. In order to access the ENZ mode which lies outside the light line we probe utilising an Otto prism coupling configuration, schematically shown in the inset of Fig. 6a, considering the film of thickness d to occupy a region −d < z < 0, vacuum to occupy the regions z < −d and 0 < z < h and a high refractive index (n = 2.4) KRS5 prism to occupy the region z > h. Effects of the illumination setup are removed from the result by normalising the fields to the power transmitted by the prism. Results are shown in Fig. 6a. In the local case the result is well known, the mode frequency tends to the zone-centre ω L as d → 0 [38]. The nonlocal result initially tracks the local one before red-shifting for films less that 10nm thick, following the effective edge of the Reststrahlen band ω 1 . One of the most interesting features of ENZ modes is the large normal field enhancement achievable thanks to the high index contrast between the film and the cladding. In Fig. 6b we show the magnitude of the normal component of the electric field at the film centre, normalised over incoupled power, as a function of the inverse 3C-SiC film thickness using both a local and nonlocal theory. In the local case, as expected a linear increase in the field strength is observed. In the nonlocal case, an increase in the enhancement is observed. This can initially come as a surprise, given that nonlocal effects normally play an important role in reducing the field enhancement by smearing the charge distribution. Such an effect is another consequence of the propagative nature of longitudinal waves in polar dielectric, which couples the fundamental longitudinal resonance of the layer to the incoming radiation. This can be verified by noticing that the transverse, photonic component of the normal electric field instead decreases above 1/d ≈ 0.2nm −1 .

Discussion and Conclusions
We studied the emergence of nonlocal effects in polar dielectric systems, demonstrating they become important for technologically relevant nanostructures. In particular we demonstrated that the negative dispersion of optical phonons, and the existence of a finite transverse frequency bounding the Reststrahlen band from below, lead to nonlocal phenomenology remarkably different from the well-known case of plasmonic nonlocality. The application of our theory to dielectric nanospheres and thin films demonstrated the appearance of resonant longitudinal modes which can couple with the object's fundamental resonances, leading to strong coupling or the appearance of Fano resonances. The propagative na-ture of the longitudinal modes also leads to a hydrodynamic Kreibig damping and it could increase the maximal field enhancement of ENZ modes. It remains an open question whether at least some of the nonlocal phenomenology due to the propagating longitudinal modes can be captured by using local models as in the plasmonic case [65], thus allowing for simple integration into existing commercial numerical software. We hope the general theory we developed will be of immediate utility to the growing community investigating SPhPs localised at the nano-scale, also providing an alternative route toward the generation of longitudinaltransverse SPhPs which have recently been proposed as a mechanism to efficiently excite SPhPs through electrical currents [36].
In this Appendix we give some details of the vectorial Mie theory used to calculate the optical response of a dielectric sphere. Although the case of an homogeneous sphere in vacuum can be treated using simplified Pekar-Ridley boundary conditions, in the general case the mechanical ABC render it necessary to construct a set of spherical vector harmonics which also yield spherical vector harmonics under application of the divergence and curl operators. A suitable choice are where Y plm (θ, φ) are scalar spherical harmonics with p = {e, o} determining the parity in φ [66] Y elm = P lm (cos θ) cos mφ, Y olm = P lm (cos θ) sin mφ, (A4) and P lm (cos θ) are the associated Legendre polynomials. Any field F can the be expanded as in which the expansion coefficients can be found exploiting the orthogonality of the spherical vector harmonics and analogously for V plm (r) [66]. Modes proportional to Φ plm are always transverse while those proportional to Y plm , Ψ plm can be decomposed into transverse and longitudinal components satisfying respectively. We can then expand the impinging plane wave as where E 0 is the field intensity and k the wavenumber outside the sphere, j l (ρ) is a spherical Bessel function of the first kind and ρ = kr. We can expand the scattered field from the particle as where we enforced transversality through Eq. (A7), a l , b l are the Mie scattering coefficients and h (1) l (ρ) is a spherical Hankel function of the first type, chosen to ensure an outgoing spherical wave in the limit ρ → ∞. For a sphere of radius R, the transverse electric field inside the particle can be written as where ρ T = ǫ (ω)k 0 r and c l , d l are unknown coefficients. The LO phonon electric field inside the particle is given by where we used the longitudinal condition Eq. (A8) and ρ L = ξr, where ξ is the longitudinal wavevector defined in Eq. (14). The three sufficient boundary conditions are continuity of the tangential components of H and E, and of the normal component of ǫ ∞ E. Assuming the particle and background are non-magnetic this leads to the equation where the longitudinal wavevector ξ is defined in Eq. (14).
We can now solve this equation to calculate the charge density in a polar dielectric sphere of radius R embedded in a non-resonant dielectric with dielectric constant ǫ b . We start by writing the general solution of Eq. (B4) in the form where Y lm (θ, φ) are spherical harmonics, j l (k L r) the spherical Bessel functions and B l are coefficients to be determined. Outside the sphere q(r > R) = 0 and the electrostatic potential thus satisfies the Laplace equation ∇ 2 φ(r > R) = 0. Considering the sphere to be illuminated by a constant field E = E 0ẑ we can determine that the electric potential in the far field is given by φ (r → ∞) = −E 0 r cos θ. We can therefore expand the potential outside the sphere as where α is the sphere polarizability. The electric potential inside the sphere can be instead determined by applying Eq. (B3) and Eq. (B5) Finally we can calculate the potential field corresponding to the ionic displacement using Eq. (B1) The MBC and the continuity of the normal component of the polarization then yield the set of equations where the dash indicates a total derivative with respect to the argument. We can thus solve for the polarisability α, which takes the form in Eq. (11) once we introduce the LRA dielectric function where γ is the damping rate. This particle polarisability relates to the extinction cross section as and furthermore we can relate it to the decay rate of an emitter orientated perpendicular to the particle surface at short distance from the particle [63] in which γ 0 is the free-space emission rate and z 0 the distance of the emitter from the center of the sphere. This quantity is just the Purcell factor, related to the mode volume V eff , the free-space wavelength λ 0 , and quality factor Q through In this Appendix the electric field modes utilised in the calculations for the ENZ mode are detailed. We probe the ENZ mode utilising a prism-coupling configuration. A phonon-active medium (region 2) occupies −d < z < 0. Non-active cladding regions 1 z < −d and 3 h > z > 0 are assumed to be positive dielectrics. Region 4 z > h contains a high-index prism. Fields can be decomposed into into longitudinal and transverse components. In the non-active regions electric fields are given by in which k x is the in-plane wavevector, ω the mode frequency, B i are the unknown coefficients to be determined by application of the appropriate boundary conditions, and the out-of-plane wavevector in each region is given by where ǫ i is the layer's non-dispersive dielectric constant. As we assume illumination from medium 4 we set B a 4 = 1. In the phonon active layer characterised by the nonlocal, dispersive dielectric function ǫ (ω, k) given by Eq. (8) the spatially dispersive Helmholtz equation has two solutions for a fixed value of ω. This means there are two transverse branches of form weighted by undetermined coefficients B {a,b} ± , with wavevectors k ± = [k x , 0, iα ± ]. The polarisation field in the active layer can be calculated from Eq. (5) and Eq. (3) as The electric field of the longitudinal phonon field can be written as iξ k xẑ e ikxx+ξ(z+d) e iωt (C11) where we defined the out-of-plane wavevector from the dispersion relation in Eq. (6) and A {a,b} 2 are further undetermined coefficients. As the longitudinal mode has no associated displacement field in the absence of free charge we can calculate the polarisation field as Finally we apply the boundary conditions at the layer boundaries. Starting with the MBC we find that continuity of the in-plane magnetic field yields while that of the in-plane electric field yields