Theory of Optical Nonlocality in Polar Dielectrics

TL;DR

This paper develops a nonlocal theoretical framework for phonon polaritons in polar dielectrics, revealing unique phenomena in nanostructures that differ from plasmonic systems, with implications for mid-infrared nanophotonics.

Contribution

It introduces a nonlocal theory for phonon polaritons in polar dielectrics, addressing limitations of local models and uncovering new nonlocal effects in nanostructures.

Findings

Nonlocal effects significantly alter phonon polariton behavior in nanostructures.

Polar dielectrics exhibit negative dispersion of phononic modes, unlike plasmonic systems.

The theory explains phenomena in nanospheres and thin films that local models cannot capture.

Abstract

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.

Equations

$$\mathbf{D}\left(\mathbf{r},\omega \right)={ϵ}_0{ϵ}_{\mathrm{LRA}}\left(\mathbf{r},\omega \right)\mathbf{E}\left(\mathbf{r},\omega \right),$$

$${ϵ}_{\mathrm{L}}\left(\omega ,k\right)={ϵ}_{\infty }\frac{{\omega }_{\mathrm{L}}^2-\omega \left(\omega +i\gamma \right)-{\beta }_{\mathrm{L}}^2{k}^2}{{\omega }_{\mathrm{T}}^2-\omega \left(\omega +i\gamma \right)-{\beta }_{\mathrm{L}}^2{k}^2}.$$

$${ϵ}_{\mathrm{T}}\left(\omega ,k\right)={ϵ}_{\infty }\frac{{\omega }_{\mathrm{L}}^2-\omega \left(\omega +i\gamma \right)-{\beta }_{\mathrm{T}}^2{k}^2}{{\omega }_{\mathrm{T}}^2-\omega \left(\omega +i\gamma \right)-{\beta }_{\mathrm{T}}^2{k}^2}.$$

$$\omega =\sqrt{{\omega }_{\mathrm{P}}^2+{\beta }^2{k}^2},$$

$$l=\frac{1}{\mathrm{Im}\left[k\right]},$$

$${\beta }_{\mathrm{L}}k=\sqrt{{\omega }_{\mathrm{L}}^2-{\omega }^2}-\frac{i\omega \gamma }{2\sqrt{{\omega }_{\mathrm{L}}^2-{\omega }^2}}+\mathcal{O}\left({\gamma }^2\right).$$

$$l_{\mathrm{ph}}=\frac{2{\beta }_{\mathrm{L}}\sqrt{{\omega }_{\mathrm{L}}^2-{\omega }^2}}{\omega \gamma }.$$

$$l_{\mathrm{pl}}=\frac{\beta }{\sqrt{{\omega }_{\mathrm{P}}^2-{\omega }^2}}.$$

$${\gamma }_{\mathrm{NL}}=\gamma +\frac{\mathrm{A}{\beta }_{\mathrm{L}}}{\mathrm{R}},$$

$$r_{\mathrm{NL}}=\frac{\mathrm{A}{\beta }_{\mathrm{L}}}{\gamma }\approx 6.1\mathrm{nm},$$

$${\omega }_n=\sqrt{{\omega }_{\mathrm{L}}^2-{\left(\frac{n\pi {\beta }_{\mathrm{L}}}{d}\right)}^2},$$

$$\left[{\omega }_{\mathrm{T}}^2-\omega \left(\omega +i\gamma \right)+{\beta }_{\mathrm{L}}^2\nabla \left(\nabla \cdot \right)-{\beta }_{\mathrm{T}}^2\nabla \times \nabla \times \right]\mathbf{X}=\frac{\mu }{\rho }\mathbf{E},$$

$$\nabla \cdot \bar{\mathbf{\tau }}=\nabla \cdot \left[{\beta }_{\mathrm{T}}^2\left(\nabla \mathbf{X}+{\left(\nabla \mathbf{X}\right)}^{\mathrm{T}}\right)+\left({\beta }_{\mathrm{L}}^2-2{\beta }_{\mathrm{T}}^2\right)\nabla \cdot \mathbf{X}\bar{\mathrm{I}}\right],$$

$$\mathbf{P}=\mu \mathbf{X}+{ϵ}_0\left({ϵ}_{\infty }-1\right)\mathbf{E},$$

$${\mathbf{P}}_{\mathrm{L}}=\left[{ϵ}_0\left({ϵ}_{\infty }-1\right)+\frac{{\mu }^2/\rho }{{\omega }_{\mathrm{T}}^2-\omega \left(\omega +i\gamma \right)-{\beta }_{\mathrm{L}}^2{k}^2}\right]{\mathbf{E}}_{\mathrm{L}}.$$

$$\mathbf{D}={ϵ}_0ϵ\left(\omega ,k\right)\mathbf{E}={ϵ}_0\mathbf{E}+\mathbf{P},$$

$${ϵ}_{\mathrm{L}}\left(\omega ,k\right)={ϵ}_{\infty }+\frac{{\mu }^2/\rho {ϵ}_0}{{\omega }_{\mathrm{T}}^2-\omega \left(\omega +i\gamma \right)-{\beta }_{\mathrm{L}}^2{k}^2},$$

$${\omega }_{\mathrm{L}}^2={\omega }_{\mathrm{T}}^2+\frac{{\mu }^2}{\rho {ϵ}_0{ϵ}_{\infty }}.$$

$${\int }_{\delta \Omega }\left(\mathbf{E}\times \mathbf{H}\right)\cdot \hat{\mathbf{n}}\mathrm{dS}=-{\int }_{\Omega }\left[{ϵ}_0\mathbf{E}\cdot \dot{\mathbf{E}}+\mathbf{E}\cdot \dot{\mathbf{P}}+{\mu }_0\mathbf{H}\cdot \dot{\mathbf{H}}\right]\mathrm{dV},$$

$$\displaystyle\int_{\delta\Omega}\left(\mathbf{E}\times\mathbf{H}\right)\cdot\hat{\mathbf{n}}\,\mathrm{dS}=-\int_{\Omega}\biggr{[}\epsilon_{0}\epsilon_{\infty}\mathbf{E}\cdot\dot{\mathbf{E}}+\mu_{0}\mathbf{H}\cdot\dot{\mathbf{H}}$$

$$\displaystyle\quad\quad\quad+\rho\left(\omega_{\mathrm{T}}^{2}-\omega\left(\omega+i\gamma\right)\right)\mathbf{X}\cdot\dot{\mathbf{X}}+\rho\left(\nabla\cdot\bar{\boldsymbol{\tau}}\right)\cdot\dot{\mathbf{X}}\biggr{]}\mathrm{dV}.$$

$${\int }_{\Omega }\mathrm{dV}\rho \left[\nabla \cdot {\bar{\mathbf{\tau }}}_i\right]{\dot{\mathrm{X}}}_i$$

$$={\int }_{\delta \Omega }\mathrm{dS}\rho \left[{\dot{\mathrm{X}}}_i{\bar{\mathbf{\tau }}}_i\right]\cdot \hat{\mathbf{n}}-{\int }_{\Omega }\mathrm{dV}\rho {\bar{\mathbf{\tau }}}_i\cdot \nabla {\dot{\mathrm{X}}}_i,$$

$$\mathbf{S}=\mathbf{E}\times \mathbf{H}+\rho \bar{\mathbf{\tau }}\dot{\mathbf{X}},$$

$$\alpha =4\pi {\mathrm{R}}^3\frac{\widetilde{ϵ}\left(\omega \right)-1}{\widetilde{ϵ}\left(\omega \right)+2},$$

$$\widetilde{ϵ}\left(\omega \right)=\frac{{ϵ}_{\mathrm{LRA}}\left(\omega \right)}{1+{\delta }_{\mathrm{NL}}},$$

$${\delta }_{\mathrm{NL}}=\frac{1}{\xi \mathrm{R}}\frac{{ϵ}_{\mathrm{LRA}}\left(\omega \right)-{ϵ}_{\infty }}{{ϵ}_{\infty }}\frac{j_1\left(\xi \mathrm{R}\right)}{j_1^{\prime }\left(\xi \mathrm{R}\right)},$$

$${\xi }^2=\left[{\omega }_{\mathrm{L}}^2-\omega \left(\omega +i\gamma \right)\right]/{\beta }_{\mathrm{L}}^2.$$

$${\sigma }_{\mathrm{ext}}=\frac{1}{\pi {\mathrm{R}}^2}\frac{\omega }{c}\left[\frac{{\omega }^3}{6\pi c^3}|\alpha {|}^2+\mathrm{Im}\left(\alpha \right)\right].$$

$${\mathbf{Y}}_{plm}\left(\theta ,\varphi \right)$$

$${\mathbf{\Psi }}_{plm}\left(\theta ,\varphi \right)$$

$${\mathbf{\Phi }}_{plm}\left(\theta ,\varphi \right)$$

Full text

I Introduction

Nanophotonics is predicated on the ability to concentrate and control light on length scales substantially below the diffraction limit Ballarini2019 . The necessary deep sub-wavelength self-sustaining electromagnetic oscillations are possible through cycling electromagnetic energy into the kinetic energy of charged particles Khurgin2017 as in the ever-growing field of plasmonics Schuller2010 ; Brongersma2015 . Plasmonic concentrators, reliant on hybridisation between photons and a free electron gas, are inefficient for operating frequencies out of resonance with the gases plasma frequency which is typically in the near ultra-violet. A promising alternative in the mid-infrared are surface phonon polariton (SPhP) light concentrators, which permit deep sub-wavelength light confinement by storing electromagnetic energy as vibrations of the crystal lattice Greffet2002 ; Hillenbrand2002 ; Caldwell2015a . Localization of SPhPs in user defined resonators Caldwell2013 allows for extremely small mode volumes Ellis2016 ; Gubbin2017 in highly tuneable low-loss modes Spann2016 ; Gubbin2016 ; Dunkelberger2018 ; Passler2018 , with applications in sensing Berte2018 , nonlinear optics Gubbin2017b ; Razdolski2018 , and the creation of nanophotonic circuitry Li2016 ; Li2018 ; Chaudhary2019 .

Typical electromagnetic theories are parameterised by frequency dependant dielectric functions in which a spatially local relationship between electric field and polarization is implicit. This approximation is known to fail in plasmonic systems of nanometric length scale Raza2011 ; Fernandez-Dominguez2012 , when longitudinal plasma waves induced in the electron gas by the transverse photon field at the particle boundary become significant. Ignoring this nonlocal effect leads to erroneous predictions of modal frequencies and significant overestimation of the achievable field confinement Ciraci2012 ; Mortensen2014 ; Wubs2017 ; Maack2017 .

Polar dielectrics support longitudinal optical (LO) phonons, analogous to longitudinal plasma waves in noble metals. As SPhP concentrators grow smaller the dielectric local theory used to model them is also expected to break down, as LO phonons hybridise with the photon field and perturb the system response. Far from being a nuisance, hybridisation of LO phonons and SPhPs has been suggested as a possible path toward electrical injection of SPhPs, and was recently demonstrated by using elongated silicon carbide (SiC) polytypes Gubbin2019 , where a weak phonon mode Paarmann2016 naturally crosses the SPhP dispersion.

In this article we develop a macroscopic nonlocal theory describing the optical response of polar crystals. Initially we use it to study 3C-SiC nanospheres discovering rich nonlocal phenomenology. We then study epsilon-near-zero (ENZ) resonances in AlN thin films and build on these results to model crystal hybrids comprised of aluminium nitride (AlN) / gallium nitride (GaN) atomic-scale superlattices. Our method provides an agile, numerically lightweight toolset to simulate and understand the nonlocal response of nanoscale polar structures, without the drawbacks of numerically taxing molecular dynamics Chalopin2012 or density functional calculations Paudel2009 ; Ratchford2019 . By comparison to experimental data we demonstrate our approaches validity for a variety of systems down to the nanometre scale. Our theory allows us to explain recently reported discrepancies between predictions of the dielectric theory and the experimental data Ratchford2019 , providing the first experimental evidence of nonlocal effects in polar dielectrics and demonstrating its relevance as a design tool for SPhP devices. For the sake of clarity most of the technical details of calculations have been placed into the Appendices.

II Nonlocal dielectric response

The optical response of polar dielectrics can be described in the local response approximation (LRA) by assuming the electric $\mathbf{E}\left(\mathbf{r},t\right)$ and displacement $\mathbf{D}\left(\mathbf{r},t\right)$ fields are spatially local

[TABLE]

in which $\epsilon_{\mathrm{LRA}}\left(\mathbf{r},\omega\right)$ is the frequency and position dependant local dielectric function. This means that the material response at $\mathbf{r}$ depends solely on the driving field at $\mathbf{r}$. Although such a formalism can be used to study SPhPs in inhomogeneous systems Gubbin2016b , it remains accurate only while dispersion of the material properties can be neglected. To move beyond this local theory we model the polar crystal as an isotropic continuum coupled to Maxwell’s equations. Such an approach, detailed in Appendices A and B, leads to a theory describing a dispersive LO phonon branch existing at all the frequencies $\omega$ and (complex) wavevectors $\mathbf{k}$ satisfying $\epsilon_{\mathrm{L}}\left(\omega,k\right)=0$ Li1998 , where $\epsilon_{\mathrm{L}}$ is the longitudinal dielectric function

[TABLE]

Transverse fields in the lattice are instead described by the transverse dielectric function

[TABLE]

and satisfy dispersion relation $\epsilon_{\mathrm{T}}\left(\omega,k\right)\omega^{2}=c^{2}k^{2}$. In Eqs. (2)-(3) $\beta_{\mathrm{L}}$ and $\beta_{\mathrm{T}}$ are phenomenological velocities describing the LO and TO phonon dispersions respectively, both of which are treated in a quadratic approximation TralleroGiner1992 . Both velocities may be extracted from ab initio calculations or experimental measurements of the polar lattice phonon dispersion, by fitting a quadratic expansion to the low-wavevector region.

Such an approach is analogous to established theories in nonlocal plasmonics, where a hydrodynamic description of the electron gas is typically employed Fuchs1971 ; Ciraci2013 ; Schnitzer2016 , leading to smearing of electromagnetic hotspots, plasmon spill-in induced frequency shifts, and the emergence of confined longitudinal modes GarciadeAbajo2008 ; McMahon2009 . Notwithstanding the technical similarities, we find the impact of nonlocality on the optical response of phonon polaritons qualitatively differs from nonlocal plasmonics. This can be intuitively understood by considering Fig. 1a, which illustrates longitudinal mode dispersion in plasmonic and phononic systems. In a metal the bulk plasma frequency blue shifts with increasing wavevector

[TABLE]

where $\beta$ is a characteristic velocity governing the dispersion of the plasma wave. This equation has real wavevector solutions for $\omega>\omega_{\mathrm{P}}$, illustrated by grey shading in Fig. 1a. Surface plasmons, mediated by the negative dielectric function of the metal, are supported in the region $\omega_{\mathrm{SP}}<\omega_{\mathrm{P}}$, illustrated by the green shading. Here Eq. (4) only has evanescent solutions and longitudinal plasma waves decay with typical skin depth of order $1\;\mathrm{\AA}$.

Optical phonons exhibit a red shift with increasing wavevector, as clear from Eq. 2. This means propagative longitudinal phonons coexist with the region where the dielectric function is negative, illustrated by the hatched region in Fig. 1a. Local theories imply delta-like screening charges at the particle surface. The charge profile is smoothed in a nonlocal theory by longitudinal modes excited at the particle edge which push charge away from the surface Ciraci2013 . In metallic particles this electronic spill-in causes a nonlocal change in electron density, leading to blue-shifted mode frequencies for $\omega<\omega_{\mathrm{P}}$ Raza2013 . Additional peaks appear for $\omega>\omega_{\mathrm{P}}$, where Eq. (4) has real wavevector solutions and the longitudinal waves can propagate into the sphere. Peaks correspond to quantised plasma waves in the sphere Raza2011 ; Rojas1988 . In the polar case we observe equivalent effective ionic charge spill-in, with the difference that the equation determining the longitudinal modes $\epsilon_{\mathrm{L}}(\omega,k)=0$ has real wavevector solutions for $\omega<\omega_{\mathrm{L}}$. The modes responsible for the spill-in are thus discrete propagative modes occupying the same spectral region as SPhPs.

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 10 nm diameter gold and 3C-SiC spheres in Fig. 1b and Fig. 1c respectively. The theory underlying this calculation is presented in the Appendices.

II.1 Nonlocal Length Scales

Nonlocal systems can be assessed by a simple figure of merit, the skindepth of nonlocal excitations

[TABLE]

where the wavevector is assumed perpendicular to the surface of the metal or polar crystal. To understand this we consider the LO phonon dispersion, given by the zeros of the longitudinal dielectric function in Eq. 2. The quality factor of optical phonon resonances is usually large enough to justify a lowest order expansion in the loss rate around $\gamma=0$: $\gamma\ll\omega,\omega_{\mathrm{L}}$. We obtain

[TABLE]

In the propagative phonon regime $\omega<\omega_{\mathrm{L}}$ the square-root yields a real quantity and the imaginary wavevector is given by the second term, leading to length scale

[TABLE]

In the plasmonic case, using Eq. 4, we obtain the equivalent of Eq. (6), but with imaginary and real parts inverted. For $\omega<\omega_{\mathrm{P}}$ we find

[TABLE]

Evaluating at the Fröhlich frequency in each material we find $\mathrm{l_{ph}}\approx 12$ nm for 3C-SiC and $\mathrm{l_{pl}}\approx 1\;\mathrm{\AA}$ for gold.

III Scattering from SiC nanospheres

In Fig. 2a we plot the extinction cross section $\sigma_{\mathrm{ext}}$, for 3C-SiC spheres of $5,10$ and $20$ nm radius, calculated using an extended Mie theory detailed in Appendix C. In the local case deep-subwavelength spheres exhibit a dipolar Fröhlich resonance where $\epsilon_{{\mathrm{LRA}}}\left(\omega_{\mathrm{F}}\right)=-2$, at $\omega_{\mathrm{F}}\approx 934\;\mathrm{cm}^{-1}$ for 3C-SiC. In nonlocal extinction spectra this resonance remains, however in smaller particles additional closely spaced peaks appear for $\omega<\omega_{\mathrm{L}}$. These are the result of the screening charge induced at the particle surface. The induced charges repulse, producing density waves which spread from the interface into the nanoparticle bulk. As in this spectral region $\epsilon_{\mathrm{L}}\left(\omega,k\right)=0$ has real wavevector solutions the excited longitudinal waves are propagative and have a discrete spectrum. These peaks become prominent when the particle radius reduces below $\approx 10$ nm, in agreement with the prediction of Eq. 7.

As the sphere radius decreases we observe level repulsion and the onset of strong coupling between the Fröhlich resonance and the longitudinal modes. In a bulk crystal such coupling is forbidden due to the orthogonal nature of longitudinal and transverse modes, which are only coupled by the system’s interfaces. Accordingly, the effect becomes more visible in Fig. 2b where we consider smaller radii and focus on the peak region. Note that this phenomenon is absent in the plasmonic case, being due to the existence of discrete propagative longitudinal modes which can become resonantly coupled to the Fröhlich resonance. The coupling is proportional to the SPhP line-shape and is therefore enhanced for longitudinal modes near to resonance with the SPhP Giannini2011 .

When the SPhP is not on resonance with the LO modes coupling decreases and we instead recover behaviour typical of dark state hybridisation, in which a broad bright state hybridises with narrower dark states resulting in a Fano scattering spectraGiannini2011 . Here the Fröhlich resonance, radiative in the far field, plays the part of the bright state while longitudinal modes are dark. This behaviour is analogous to the weak phonon coupling previously observed in larger 4H-SiC particles, where the weak phonon frequency was determined by the lattice periodicity rather than particle morphology Gubbin2019 ; Ellis2016 . Fano interference acts to reduce the bright states far-field scattering, while increasing absorption as energy is funnelled into non-radiative dark states Gennaro2014 ; Simoncelli2018 . This is visible in the comparison of scattering and absorption efficiencies for a $4$ nm radius sphere in Fig. 2c, where both cross sections are normalised to unity in order that they may be presented on the same vertical scale. On this plot vertical lines indicate the pure longitudinal resonances, calculated assuming the longitudinal electric field vanishes at the particle boundary. As shown in the inset, the longitudinal modes at $905$ and $950\;\mathrm{cm}^{-1}$ suppress scattering while enhancing absorption, giving a Fano line profile.

III.1 Frequency Shifts

In plasmonic systems the dominant nonlocal effect is typically a blue shift in mode frequency Raza2011 , manifested in a quasi-static model as an increase in the effective local dielectric function. In the case of polar dielectrics the Fröhlich resonance also undergoes a nonlocal shift, in this case to the red, but of much smaller amplitude. This can be verified from Fig. 3a, in which we compare the absolute plasmonic and phononic frequency shifts from their local values $\omega_{\mathrm{F}}$ as a function of the absolute value of the dimensionless parameter quantifying nonlocal effects $\delta_{\mathrm{NL}}$, precisely defined in Appendix C. The reduced shift in the polar dielectric case 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 accommodate a change in the dielectric function due to a certain value of the nonlocal parameter $\delta_{\mathrm{NL}}$ by a smaller shift in frequency.

III.2 Enhanced Nonlocal Broadening

Another nonlocal effect often observed in the optics of small plasmonic particles is a size-dependent damping, initially described by Kreibig Kreibig1969 . This can be introduced phenomenologically in a hydrodynamic nonlocal theory as a diffusive term Mortensen2014 ; Ciraci2017 . Any broadening is size-dependent, requiring electrons to diffuse during the bare mode lifetime over a distance comparable to the field confinement length. The reason why size-dependent damping does not emerge from a pure hydrodynamic model in the absence of diffusion is because the longitudinal plasma waves excited in the electron gas obey Eq. (4), and are thus evanescent. In the polar system studied in this paper the LO modes excited by the SPhPs are instead propagative, and can thus transport energy away from the interface. This allows them to act as an additional loss channel, leading to hydrodynamic Kreibig damping, which remains important even for relatively large particles.

In order to calculate the size-dependent broadening we use a semi-analytical approximation of the full nonlocal Mie theory, described in Appendix C, which yields the frequency dependant polarisability of the sphere. SPhP mode-frequency and linewidth are then extracted by a Lorentzian fit. In Fig. 3b we plot the resulting quality factor of the fundamental Fröhlich resonance for 3C-SiC spheres with radii between $10$ and $1000$ nm. In this region the particle is sufficiently large to support many closely spaced longitudinal modes, whose finite linewidth 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\;\mu$m, coupling to these non-resonant modes leads to a broadening relative to the local result which grows with decreasing particle radius. We can form a model analogous to that of Kreibig to fit the data, assuming nonlocal damping rate of form

[TABLE]

where $\gamma$ is the local damping rate and $\mathrm{R}$ is the nanosphere radius.

The dimensionless constant $\mathrm{A}$ describes the magnitude of the size-dependent damping. Fitting this formula to the numerical data we find $\mathrm{A}=0.03$. We can utilise this to calculate a typical length scale for which damping via LO phonon emission becomes dominant

[TABLE]

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 Mortensen2014 .

IV Nonlocality in planar nitride systems

IV.1 Nonlocality in AlN ENZ modes

Now we investigate nonlocal effects in freestanding ultrathin AlN films. Details of the approach used are provided in Appendix E. Optically thick films support degenerate SPhP modes at each interface. When the thickness $d$ decreases, eventually these hybridise yielding two distinct SPhP branches which are symmetric and antisymmetric superpositions of the modes on each interface. The symmetric mode is blue-shifted from the mode of the thick film and in the $d\to 0$ limit exists at the asymptotic zone-centre LO phonon frequency. As the dielectric function of the film vanishes at this frequency the mode is often referred to as an epsilon-near-zero (ENZ) mode Campione2015 ; Passler2018 .

Just like the SiC spheres studied in the previous section, a thin film also acts as a Fabry-Perot cavity for longitudinal phonons, supporting a series of discrete modes with quantised out-of-plane wavevector, satisfying $\eta=\frac{n\pi}{d}$, with $n$ a positive integer. For wavevectors in the film plane $k_{x}\ll\eta$ the LO phonon modes have frequencies

[TABLE]

a result analogous to that observed for ENZ modes in thin conductive films, where additional resonances obeying a similar equation were observed above the film plasma frequency deCeglia2018 . In a sufficiently thin sample the fundamental LO resonance $\omega_{1}$ can differ appreciably from the zone-center frequency $\omega_{\mathrm{L}}$, leading to a shift in the upper edge of the Reststrahlen band. This is illustrated for an AlN film by the dots in Fig. 4, which show the $n=1$ solution of Eq. 11 as a function of inverse film thickness. As the ENZ mode lies close to the LO phonon frequency, 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_{\parallel}c=1200\;\mathrm{cm}^{-1}$ for an AlN film as a function of inverse film thickness in the local and nonlocal cases. In order to access the ENZ mode, which lies outside the light line, we utilise an Otto prism coupling configuration, schematically shown at the top of Fig. 4, considering the AlN film of thickness $d$, surrounded by vacuum with a high refractive index ($n=2.4$) KRS5 prism a distance $h$ above. Results shown in Fig. 4 demonstrate that nonlocal effects shift the ENZ resonance from its $\omega_{\mathrm{L}}$ frequency, having it to follow instead the effective edge of the Reststrahlen band $\omega_{1}$ for films less that $10$ nm thick. In the following we will see that similar shifts due to the ENZ confinement can be observed also by embedding the AlN between other materials, at frequencies for which the Reststrahlen bands do not overlap.

IV.2 Nonlocality in AlN/GaN atomic-scale superlattices

Recently the fabrication and optical characterisation of atomic-scale superlattices or crystal hybrids, comprised of alternating nanoscopic layers of GaN and AlN grown on a semi-insulating SiC substrate were reported Ratchford2019 . Such a work demonstrated the failure of a standard local response theory to describe the observed optical properties of the superlattice, requiring instead density functional perturbation theory to correctly reproduce its features Paudel2009 .

Using a nonlocal extension of the scattering matrix approach we simulated the two heterostructures, labelled A and B, characterised in Ref. Ratchford2019 . Heterostructure A consists of 50 alternating layers of GaN ($2.2$ nm)/ AlN ($2.2$ nm) on a $50$ nm AlN buffer layer. Heterostructure B consists of 500 alternating layers of GaN ($1.4$ nm)/ AlN ($1.2$ nm). Both heterostructures are grown on a deep SiC substrate. In Fig. 5 we show experimental data, courtesy of the authors, with blue-solid lines and we compare them to our optimised nonlocal model (red dashed lines) and a separately optimised local theory (green dot-dashed lines). More details on the fitting procedure are presented in Appendix F.

Not only does the nonlocal theory provide a better fit, but we can recognise two sets of features which are correctly reproduced by the nonlocal theory while not captured at all by the local one. The first is photonic hybridisation with longitudinal modes, leading to visible peaks in heterostructure A around $835\;\mathrm{cm}^{-1}$ and $865\;\mathrm{cm}^{-1}$. These are analogous to the nonlocality-induced additional peaks visible in Fig. 2. The second is the strong redshift of the resonance around $855\;\mathrm{cm}^{-1}$ for heterostructure B, which we already discussed in Fig. 4 for the freestanding film. As previously noted, our nonlocal theory treats phonon dispersions as quadratic. This approximation reproduces well the phonon dispersion near the zone-centre phonon frequencies, but it diverges from the physical dispersion closer to the Brillouin zone edge Karch1994 . One consequence of the quadratic approximation is the failure of our model to correctly reproduce the dip observed at $795\;\mathrm{cm}^{-1}$ in Fig. 5b. The modes confined in the layers of the narrower structure have in fact larger wavevectors through Eq. (11), thus probing a region in which the quadratic dispersion underestimate the longitudinal phonon frequency, pushing the longitudinal modes at lower frequencies. The eventual breakdown of the quadratic approximation is also why we constrain the plot region in Fig. 5 to around the SiC Reststrahlen region. It is a priori possible to overcome such a problem by considering wavevector-dependent, piecewise constant effective phonon velocities, able to reproduce the phonon dispersion over larger portions of the Brillouin zone, but the advantages of such an approach have to be weighted against the increase of required fitting parameters.

The calculated parameter for the AlN longitudinal velocity $\beta_{\mathrm{L}}^{\mathrm{AlN}}=5.1\times 10^{5}\;\mathrm{cm\;s}^{-1}$ and damping $\gamma^{\mathrm{AlN}}=10.3\;\mathrm{cm}^{-1}$ leads through Eq. (7) to a predicted skin-depth of $\mathrm{l}_{\text{ph}}=1.5\;\mathrm{nm}$ at the $n=1$ resonant frequency $855\;\mathrm{cm}^{-1}$. As for dielectric nanospheres, also in this case the use of Eq. (7) provides the right order of magnitude for the length scale at which nonlocal phenomena become important.

V Discussion and Conclusions

We studied the emergence of nonlocal effects in isotropic polar dielectric systems, demonstrating their importance for technologically relevant nanostructures. 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 a rich and unique nonlocal phenomenology. In particular opening the possibility to resonantly couple photonic excitations to propagative longitudinal modes makes it possible to engineer hybrid longitudinal-transverse SPhPs, recently highlighted as an important stepping stone toward the realisation of electrically pumped SPhP devices Gubbin2019 , in generic nanoscale objects.

By modelling recently published experimental data, demonstrating that some of their unexplained features result directly from nonlocal effects, showing their immediate importance to the design of dielectric nanodevices, and providing a first experimental verification that our theory provides quantitatively correct results for systems down to few atomic lattices in size. This can mean that the physics at that level is already acceptably well described by a continuous theory or, alternatively, that the average effect of the microscopic degrees of freedom can be described by a renormalization of the effective nonlocal system parameters. Future experimental and numerical investigations will be required to gain a better understanding of the boundary of applicability of our approach. When used in conjunction with experimental data or microscopic calculations Chalopin2012 ; Paudel2009 , required to fix some of the model parameters, our theory thus allows for agile design and simulation of nanoscopic polar devices. Although in this paper all materials are treated as isotropic our model could be extended to consider anisotropic systems. This would entail numerical calculation of the modal dispersion, replacing the ionic equation of motion with it’s anisotropic analogue. We were able to implement our nonlocal theory in a number of numerical and analytical techniques adapted to different geometries (Mie theory, quasi-static, scattering matrix). 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 Luo2013 , thus allowing for simple integration into existing commercial numerical software.

Acknowledgements

S.D.L. is a Royal Society Research Fellow. The authors acknowledge support from the Innovation Fund of the EPSRC Programme EP/M009122/1 and from the Royal Society grant RGF\EA\181001.

Appendix A Continuum dielectric model

We treat our polar crystal as an isotropic lattice, described in the continuum limit by Trallero-Giner1992 ; Li1998

[TABLE]

in which $\mathbf{X}$ is the relative ionic displacement, $\omega_{\mathrm{T}}$ is the transverse optical (TO) phonon frequency and $\rho$ and $\mu$ are the effective mass and charge densities. We suppress space (wavevector) and time (frequency) dependencies where possible. For simplicity we consider a single loss rate $\gamma$, independent of both frequency and polarization. The spatially dispersive terms in Eq. (12) are calculable as the divergence of a tensor $\bar{\boldsymbol{\tau}}$

[TABLE]

in which $\bar{\mathrm{I}}$ is the identity tensor. In a model of acoustic phonons $\bar{\boldsymbol{\tau}}$ is the mechanical stress tensor, parameterised by the Lamé coefficients which describe the acoustic phonon velocities. As here we consider optical phonons the coefficients $\beta_{\mathrm{L}},\beta_{\mathrm{T}}$ are instead phenomenological velocities which describe LO and TO phonon dispersions respectively TralleroGiner1992 . 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. We can then introduce the total polarisation $\mathbf{P}$, which is a sum of the of the ionic polarisation $\mu\mathbf{X}$ and the high-frequency non-resonant background

[TABLE]

in which $\epsilon_{\infty}$ is the high-frequency dielectric constant.

The longitudinal and transverse dielectric functions, given in Eq. (2) and Eq. (3) of the main text, can then be derived. For the longitudinal mode $\nabla\times\mathbf{X}=0$, and Fourier transforming Eq. (12) we can write

[TABLE]

Using the standard constitutive relation

[TABLE]

we find

[TABLE]

which simplifies to that given in Eq. (2) defining the LO phonon frequency

[TABLE]

A similar relation for the transverse dielectric function $\epsilon_{\mathrm{T}}\left(\omega,k\right)$ in Eq. (3) can be derived considering instead $\nabla\cdot\mathbf{X}=0$ in Eq. (12).

Appendix B 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 $\mathbf{E}$, and $\mathbf{H}$. As a nonlocal formalism entails consideration of additional fields in each layer, additional boundary conditions (ABC) are necessary to uniquely determine the model. In order to find the correct boundary conditions we use the ionic equation of motion Eq. (12) to derive a hybrid Poynting vector describing the total energy transported in both the electromagnetic and mechanical fields. We start from the Poynting’s theorem

[TABLE]

which relates the change in energy enclosed in volume $\mathrm{\Omega}$ with the flux through it’s surface $\delta\Omega$. Substituting Eq. (12) into Eq. (19) via the constitutive relation Eq. (14), an approach proposed by Loudon Loudon1970 , we can write

[TABLE]

The final term in this equation describes mechanical transport. It can be recast as

[TABLE]

where the latter term describes a nonlocal adjustment of the energy density and the former is an energy flux. Combining the flux term with the local Poynting flux from Eq. (19) we derive the hybrid Poynting vector

[TABLE]

where the tensor $\bar{\boldsymbol{\tau}}$ is defined by Eq. (13). In order to ensure that interfaces act as neither energy sources or sinks, we impose continuity of the normal component of $\mathbf{S}$. This requires, in addition to the standard MBC, continuity of $\mathbf{X}$ and of the normal components of the tensor $\bar{\boldsymbol{\tau}}$.

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 Pekar1959 , by application of symmetry conditions to the fields Fuchs1971 or by considering explicitly the nonlocal surface potential near the boundary Ruppin1984 . The most commonly utilised ABC are the Pekar-Ridley Pekar1959 and Fuchs-Kliewer Fuchs1971 which, assuming for definiteness an interface lying in the $xy$ plane read $\mathbf{X}=0$ and $\partial_{z}\mathrm{X}_{\{x,y\}}=\mathrm{X}_{z}=0$ respectively. Both of these choices satisfy $\rho\bar{\boldsymbol{\tau}}\dot{\mathbf{X}}\cdot\hat{\mathbf{z}}=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 $\mathrm{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. We then use Fuchs-Kliewer type ABC as they are derived by symmetry arguments, and they are equivalent to assuming specular reflection from the interface, which is the relevant symmetry in our case.

Appendix C Scattering from dielectric spheres

In order to include nonlocal effects in the treatment of the optical response of small dielectric spheres we used an extended Mie theory, conceptually similar to previous nonlocal extensions of Mie theory utilised to describe nonlocal effects in noble metallic systems Raza2011 and described fully in the Supplementary information. In order to test our calculations, and to provide a better understanding of the results, we also used a semi-analytic model valid in the quasi-static limit, that is assuming the impinging radiation field does not vary appreciably over the sphere diameter. Raza2013 . This is a good approximation for the parameters considered in this paper, with radii $\mathrm{R}$ normalised over the resonant wavelength of the order of $0.01$. For the parameters presented in Fig. 2b the extinction cross section given in Eq. (27) is indistinguishable from the full Mie theory.

In the quasi-static limit the polarisability of a sphere in vacuum has the form

[TABLE]

where

[TABLE]

is the local dielectric function renormalised by nonlocal effects through the nonlocal parameter

[TABLE]

with $j_{1}(x)$ and $j_{1}^{\prime}(x)$ a spherical Bessel function of the first kind and its derivative, and $\xi$ is the longitudinal mode wavevector defined by the equation

[TABLE]

Within this formalism we can quickly evaluate the extinction cross section

[TABLE]

Further details on both the extended Mie theory and on the semi-analytic quasi-static theory can be found in Appendix D.

In the 3C-SiC case we use transverse optic phonon frequency $\omega_{\mathrm{T}}=796.1\;\mathrm{cm}^{-1}$, longitudinal optic phonon frequency $\omega_{\mathrm{L}}=973\;\mathrm{cm}^{-1}$ Caldwell2015a , high-frequency dielectric constant $\epsilon_{\infty}=6.52$, damping rate $\gamma=4\;\mathrm{cm}^{-1}$, and nonlocal velocities $\beta_{\mathrm{T}}^{SiC}=9.15\times 10^{5}\;\mathrm{cm\;s}^{-1},\;\beta_{\mathrm{L}}^{SiC}=15.39\times 10^{5}\;\mathrm{cm\;s}^{-1}$ approximated from studies of phonon dispersion Karch1994 . 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 Debernardi1999 .

In the gold case we use plasma frequency $\omega_{\mathrm{p}}=72670\;\mathrm{cm}^{-1}$, damping rate $\gamma=580\;\mathrm{cm}^{-1}$ Derkachova2015 and nonlocal velocity $\beta_{\mathrm{p}}=1.08\times 10^{8}\;\mathrm{cm\;s}^{-1}$ Luo2013 ; Ciraci2012 .

Appendix D Nonlocal Mie Theory

In Appendix C we derived the extinction cross section in the quasistatic limit for a nanosphere supporting longitudinal nonlocal modes. To consider transverse nonlocality it is necessary to develop a full nonlocal Mie theory. When considering the mechanical boundary conditions it is 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 Barrera1985

[TABLE]

where $\mathrm{Y}_{plm}(\theta,\phi)$ are scalar spherical harmonics with $p=\{e,o\}$ denoting even or odd parity in $\phi$

[TABLE]

and $\mathrm{P}_{lm}\left(\cos\theta\right)$ are the associated Legendre polynomials. Single-valued vector fields vanishing at infinity ($\mathbf{F}$) can be written in this basis

[TABLE]

where expansion coefficients can be found exploiting the orthogonality of the spherical vector harmonics. For example

[TABLE]

This choice of spherical harmonics allows for segregation into longitudinal and transverse modes. Modes proportional to $\boldsymbol{\Phi}_{plm}$ are always transverse while those proportional to $\mathbf{Y}_{plm},\boldsymbol{\Psi}_{plm}$ can be decomposed into transverse and longitudinal components satisfying respectively the equations

[TABLE]

This theory can be applied to any layered spherical resonator supporting longitudinal and transverse nonlocality. As an example we calculate the fields in a single layer nanosphere of radius $R$, supporting a longitudinal mode. An impinging plane wave can be expanded in the region $r>R$ as

[TABLE]

where $\mathrm{E}_{0}$ is the field intensity and $k$ the wavenumber outside the sphere, $j_{l}(\rho)$ is a spherical Bessel function of the first kind and $\rho=kr$. The scattered field in region $r>R$ can be similarly written

[TABLE]

where we enforced transversality through Eq. (35), $a_{l},b_{l}$ are the Mie scattering coefficients and $h_{l}^{(1)}(\rho)$ is a spherical Hankel function of the first type, chosen to ensure an outgoing spherical wave in the limit $\rho\to\infty$.

The transverse field in the region $r<R$ can be written

[TABLE]

where $\rho_{T}=\sqrt{\epsilon\left(\omega\right)}k_{0}r$ and $c_{l},d_{l}$ are unknown coefficients. The longitudinal phonon electric field inside the particle is given by

[TABLE]

where we used Eq. (36) and $\rho_{L}=\xi r$, where $\xi$ is the longitudinal wavevector defined in Eq. (26).

The unknown coefficients $a_{l},b_{l},c_{l},d_{l},g_{l}$ can be calculated by application of the appropriate boundary conditions. In this case these are the standard Maxwell boundary conditions on the components of $\mathbf{E}$ and $\mathbf{H}$ parallel to the sphere surface, in addition to the boundary condition on the radial component of the ionic displacement $\mathbf{X}$. In conjunction with continuity of the radial displacement field $\mathbf{D}$ this boundary condition can be recast as one on the radial component of $\epsilon_{\infty}\mathbf{E}$. This leads to extinction and scattering cross sections

[TABLE]

Appendix E Nonlocal model of a layered dielectric

The application of the nonlocal theory to a layered structure can be performed by calculating a complete set of modes for each layer and then applying MBC and ABC to match their amplitudes at the boundaries. In the case of the four layers considered for the AlN-vacuum ENZ resonance (prism–vacuum–AlN–vacuum) this leads to $12$ independent modes: two transverse counter-propagating modes for each non-active layer, and four transverse and two longitudinal modes for the AlN slab. The AlN slab (region 2) occupies $-d<z<0$, vacuum occupies regions 1 $z<-d$ and 3 $h>z>0$, and region 4 $z>h$ contains a high-index prism. The photon fields in the multilayer stack are given by

[TABLE]

in which $k_{x}$ is the in-plane wavevector, $\omega$ the mode frequency, $\mathrm{B}_{i}$ are unknown coefficients to be determined by the boundary conditions and $\epsilon_{i}$ is the $i$th layers dispersionless dielectric function. The unit vector is given by

[TABLE]

and the out-of-plane transverse wavevector in layer $i=1,3,4$ is

[TABLE]

Assuming illumination from medium $4$ we can set $\mathrm{B}_{4}^{\mathrm{a}}=1$ and recognise $\mathrm{B}_{4}^{\mathrm{b}}$ as the reflectance.

In the nonlocal AlN film two transverse and one longitudinal mode are supported. The electric field profile for TO phonons and photons can be written

[TABLE]

where $\alpha_{\pm}$ represents the two solutions to the nonlocal dispersion relation $\epsilon_{\mathrm{T}}\left(\omega,k\right)\omega^{2}=c^{2}k^{2}$. The longitudinal mode field can be written

[TABLE]

where $\mathrm{A}^{\{\mathrm{a,b}\}}$ are further undetermined coefficients and the out-of-plane wavevector

[TABLE]

is the solution to $\epsilon_{L}\left(\omega,k\right)=0$.

From these electric fields the fields entering the boundary conditions are calculable. The magnetic field is calculable by Maxwell-Faraday equation. As additional boundary conditions we constrain the in- and out-of-plane components of the ionic displacement $\mathbf{X}$. The AlN film is clad in vacuum where $\mathbf{X}=0$. As the nonlocal model eliminates charge surface charging the normal component of the displacement field

[TABLE]

is continuous at the AlN boundary. The additional boundary condition on the perpendicular component of $\mathbf{X}=0$ can therefore be recast to one on $\epsilon_{\infty}\mathrm{E}_{z}$. The final boundary condition is $\mathrm{X}_{\parallel}=0$. Applying these conditions yields a $10\times 10$ matrix problem for the unknown coefficients $\mathrm{A,B}$.

Appendix F Fitting procedure

Transfer matrix algorithms, the tool of choice to model optical properties of planar semiconducting heterostructures, are generally unstable in the limit of thick absorbing layers. As in a nonlocal model the phonon mode out-of-plane wavevectors are very large this problem emerges even when considering optically thin films. This instability forced us to utilise a scattering matrix approach. In a scattering matrix model we start by calculating the matrix of the $p^{\text{th}}$ layer $\mathrm{s}^{(p)}$ which relates the up- and down-propagating incoming modes $\mathbf{u}^{(p)}$ and $\mathbf{d}^{(p+1)}$ to the outgoing ones $\mathbf{u}^{(p+1)}$ and $\mathbf{d}^{(p)}$. The scattering matrix for the full $P$-layer stack,

[TABLE]

can then be obtained by calculating the Redheffer star product of the single-layer matrices. We then fitted the experimental data using both local and the nonlocal theories. In both cases we allowed the layers’ thickness to vary around the nominal value in order to take into account atomic intercalation or screening Chalopin2012 at the interface, which reduces the effective thickness which can be described by bulk crystal structure. Consistently with such a picture we found optimal thickness reductions roughly equal to one lattice constant per interface.

The longitudinal phonon velocity in the non-resonant GaN layer was fixed at $\beta_{\mathrm{L}}^{GaN}=6.5\times 10^{5}\;\mathrm{cm\;s}^{-1}$, approximated from numerical phonon dispersion curves Bungaro2000 . In order to reduce the fitting parameters we neglected TO phonon velocities, which have a very limited impact on the calculated spectra, thus imposing $\beta_{\mathrm{T}}^{AlN}=\beta_{\mathrm{T}}^{GaN}=0$ also in the nonlocal theory. Our fitting procedure used two unknown parameters for the local theory and three in the nonlocal one: the shift in layer thickness $\delta$, the damping rate in the AlN $\gamma^{AlN}$, and the longitudinal velocity in the AlN layer $\beta_{\mathrm{L}}^{AlN}$, the latter being unique to the nonlocal theory. Optimal parameters are found to be $\beta_{\mathrm{L}}^{AlN}=5.1\times 10^{5}\;\mathrm{cm\;s}^{-1}$, $\delta=0.94\;\mathrm{nm}$ and $\gamma^{AlN}=10.3\;\mathrm{cm}^{-1}$.