In the Reststrahlen region, between the transverse and longitudinal phonon frequencies, polar dielectric materials respond metallically to light, and the resulting strong light–matter interactions can lead to the formation of hybrid quasiparticles termed surface phonon polaritons. Recent works have demonstrated that when an optical system contains nanoscale polar elements, these excitations can acquire a longitudinal field component as a result of the material dispersion of the lattice, leading to the formation of secondary quasiparticles termed longitudinal-transverse polaritons. In this work, we build on previous macroscopic electromagnetic theories, developing a full second-quantized theory of longitudinal-transverse polaritons. Beginning from the Hamiltonian of the light–matter system, we treat distortion to the lattice, introducing an elastic free energy. We then diagonalize the Hamiltonian, demonstrating that the equations of motion for the polariton are equivalent to those of macroscopic electromagnetism and quantize the nonlocal operators. Finally, we demonstrate how to reconstruct the electromagnetic fields in terms of the polariton states and explore polariton induced enhancements of the Purcell factor. These results demonstrate how nonlocality can narrow, enhance, and spectrally tune near-field emission with applications in mid-infrared sensing.

## I. INTRODUCTION

Surface phonon polaritons (SPhPs) are hybrid light–matter excitations, formed when a photon interacts with the optical phonon modes of a polar lattice. Like plasmons in the visible spectral region, they allow for miniaturization of photonic resonators below the diffraction limit,^{1} with applications in mid-infrared sensing,^{2} nonlinear optics,^{3–5} and the design of thermal emitters.^{6–9}

In a simple model of polar optics, the lattice is considered to be non-dispersive and is described by a frequency-dependent dielectric function parameterized by the optical phonon frequencies at the zone-center. This model works extremely well in regimes where the dispersion of the optic phonon branches can be neglected; however, when the system approaches the nanoscale, this is no longer necessarily the case. Systems where material dispersion is important are termed optically nonlocal, meaning that the applied field at one point can affect the response at another. This is a consequence of energy transport in the matter, for example, through bulk plasma waves in nano-plasmonic systems.^{10} This transfer of energy from the photon field leads to enhanced damping and spectral shifts and ultimately places an upper-limit on how tightly the field can be confined.^{11–14}

Nonlocal effects in polar systems are fundamentally different from those in plasmonic systems because of the opposite dispersion of bulk optical phonon and plasma waves. In the Reststrahlen region where the dielectric function is negative and SPhPs exist, polar crystals also support propagative Longitudinal Optical (LO) phonon excitations, allowing for resonant coupling between discrete LO phonon modes and SPhPs. These nonlocal effects in polar systems were first demonstrated in a recent work, which showed strong coupling between localized SPhPs in a 4H–SiC nanopillar and zone-folded optical phonons.^{15} The resulting excitations are hybrid modes termed longitudinal-transverse polaritons (LTPs) and have a transverse electric field from their photonic component and a longitudinal one from the phonon field. This unique property has led to LTPs being proposed as a platform for mid-infrared optoelectronics as a result of the possibility of exciting them through longitudinal electrical currents while still outcoupling to transverse free-space radiation in the far-field.^{16}

A recent series of publications have studied LTPs in more general systems, demonstrating them to be a general feature of polar resonators at the nanoscale,^{17–19} and a similar phenomenology has also recently been observed in a Yukawa fluid.^{20} These works follow the approach of nonlocal plasmonics, starting from the macroscopic Maxwell equations, introducing new macroscopic fields to describe phonon modes in the lattice and matching fields at material boundaries considering the flow of energy in the system. They have been utilized to explain the emergence of anomalous modes in complex crystal hybrid structures and macroscopic systems comprised of hundreds of nanoscopic polar layers^{17,21} and have proved able to calculate the electromagnetic response.

Polariton systems are typically described in second-quantization formalisms,^{22–25} which are attractive because they allow for a transparent understanding of energy distribution between the different light and matter excitations from which the polariton is composed and for calculation of nonlinear polariton–polariton^{26–28} or electron–polariton scattering.^{29} Although recently a second-quantization theory of LTPs was presented, this was based on a demonstration of equivalence between Maxwell’s equations and a model Hamiltonian in a planar waveguide^{30} and is not suitable for extension to describe the complex modal geometries typically utilized in SPhP optics.

In this work, we develop a full quantum theory of LTPs. Starting from the Hamiltonian of a spectrally dispersive and spatially inhomogeneous polar dielectric, we introduce a new free energy to account for the elastic energy of optical phonons. We then diagonalize this Hamiltonian, demonstrating that the resulting equations of motion for the Hopfield fields are equivalent to the nonlocal Maxwell equations.^{17} Operators for the physical fields are reconstructed from the Hopfield fields, the boundary conditions satisfied by those fields are determined, and the quantization conditions for the polaritonic modes are derived. Finally, we utilize our theory to demonstrate the nonlocal emission from dipoles embedded near a polar nanolayer.

## II. THEORY

Our starting point is the Hamiltonian density in a local non-magnetic polar dielectric, including spectral dispersion. We consider the material to be piecewise homogeneous such that we may employ Fourier analysis in each region. The system’s quantum Hamiltonian is given by

where $D\u0302(H\u0302)$ is the electric displacement (magnetic) field operator, $X\u0302$ is the ionic displacement operator, $Q\u0302$ is the momentum operator, *ω*_{L} is the longitudinal optical phonon frequency at the zone-center, *ϵ*_{∞} is the high-frequency dielectric constant, and *ρ* is the crystal mass density.^{25} The oscillator strength *κ* characterizes the width of the Reststrahlen that relates to the phonon frequencies through

where in the latter equality, we utilized the Lydanne–Sachs–Teller relation, linking the crystal’s transverse and longitudinal optic phonon frequencies (*ω*_{T} and *ω*_{L}) using the static dielectric function of the lattice *ϵ*_{st}.^{31} The Hamiltonian [Eq. (1)] accounts for the kinetic energy contribution of ions oscillating in place through the term proportional to the lattice momentum P. It does not account for lattice distortion or for energy transported in finite-wavevector phonon waves. As propagating phonons are analogous to propagating elastic waves, we account for this additional energy by considering an additional contribution to the Hamiltonian of the lattice,

where $C\u0304ijkl$ is an effective elasticity tensor satisfying the symmetry conditions^{18}

and $S\u0302ij$ are scalar stress operators, defined as

Note that this is a first order approximation to the elastic free energy, which will result in a quadratic phonon dispersion relation. This approximation is at the same level as that utilized in classical theories of polar nonlocality.^{17} The full nonlocal Hamiltonian is a sum of the two Hamiltonian densities,

### A. Equations of motion

The goal of this paper is to diagonalize Eq. (6), writing it in terms of a series of bosonic ladder operators, which describe the polaritonic eigenmodes of the system. To that end, following previous approaches,^{25} we introduce the general polaritonic operator as a linear superposition of the free fields,

where Greek symbols are Hopfield fields describing the weighting of the fields comprising the eigenmode. If $K\u0302$ is an eigenmode of $H\u0302$, it satisfies the equation of motion,

where *ω* is the polariton frequency. The lengthy quantization procedure, carried out in Appendix A, yields equations of motion for the Hopfield fields of the polariton,

where $\tau \u0304$ is the stress tensor of the polar lattice, related to the effective elasticity tensor through

and we define the composite field

This field can be interpreted as the sum of the displacement field associated with the polariton’s transverse photonic component, proportional to ** α**, and the longitudinal electric field associated with the polariton’s LO phonon component, given by the latter term. The result is a new field

**, which acts as the full electric field of the polariton.**

*θ*Note that in these equations of motion, we can immediately identify the Hopfield fields as analogs of the classical fields. The magnetic field analog is ** β**, while the electric field is

**. The other coefficients refer to the ionic position and momentum. It will be useful to the reader to keep this interpretation in mind through the rest of this article. We now show that the equations of motion lead to the dispersion relations predicted by classical approaches.**

*θ*^{17}For simplicity, we focus on isotropic crystal structures, for which the components of $C\u0304ijkl$ are given by

where the dots indicate other permutations of indices through Eq. (4). These components relate to the longitudinal and transverse phonon velocities in the quadratic dispersion limit *β*_{L}, *β*_{T} introduced phenomenologically in other works^{17,30} through

Finally, through Eq. (10), the components of the effective stress tensor can be derived,

where $Sij$ are scalar stresses of the Hopfield field *ζ*,

One class of solutions of Eqs. (9a)–(9d) is transverse fields, for which the Hopfield field for the classical electric field satisfies ∇ · ** θ** = 0. In this case, we can solve the equations of motion directly, recovering the nonlocal wave equation for transverse fields,

where *ϵ*_{T} is the transverse dielectric function of the lattice, given by

and *β*_{T} can be interpreted as the transverse optical phonon velocity in the quadratic limit.

The equations of motion also support longitudinal solutions for which the classical electric Hopfield field satisfies ∇ × ** θ** = 0. This results in the dispersion relation

where *ϵ*_{L} is the longitudinal dielectric function of the lattice, given by

where *β*_{L} is the transverse optical phonon velocity. These transverse and longitudinal dispersion relations are exactly those predicted from the classical theory, validating the approach taken in this section. The LTP excitations which are the focus of this work are hybrid modes with both transverse and longitudinal components, mixed by the boundary conditions on the Hopfield fields ** α, β, γ, ζ**, which must enforce the continuity of the quantum fields as predicted in standard theories of polar nonlocality.

^{17}These boundary conditions are derived in Sec. II C.

### B. Field construction

where the index *n* can be either discrete or continuous depending on the mode and system under study with associated eigenfrequencies *ω*_{n} and wavevectors *k*_{n}. Then, an arbitrary field operator $Y\u0302$ can be written as a sum over the polaritonic operators as

in which the expansion coefficients are given by

The electric displacement field, magnetic field, ionic displacement, electric field operator $E\u0302$, and polarization field operator $P\u0302$ can be written as

where in the final relation, we used the relation between the material polarization density operator $P\u0302=\kappa X\u0302$ and the following relation derived from Eqs. (9a)–(9b) by elimination of ** β**:

Again, note the physical significance of the polariton Hopfield fields. The displacement field is proportional to the electric Hopfield field ** θ** multiplied by the system’s transverse dielectric function, and the magnetic field operator is proportional to the magnetic Hopfield field

**.**

*β*### C. Boundary conditions

We have derived the equations of motion for the nonlocal polariton system [Eqs. (9a)–(9d)]. This equation set describes the relation between the Hopfield fields in a homogeneous system. To solve a physical problem, we also need to derive boundary conditions, describing how polaritonic Hopfield fields behave at material interfaces. As has been discussed at length in previous works, the standard boundary conditions utilized in electromagnetic theory are insufficient to describe nonlocal systems due to the introduction of additional fields describing the lattice distortion.^{17} The appropriate new boundary conditions can be derived starting from the classical Poynting relation, which links the rate of change of electromagnetic energy density in a volume Ω with the energy flux passing through enclosing surface dΩ with outgoing surface normal unit vector $n\u20d7$, given by

where we utilized the standard constitutive relation for a nonlocal dielectric to eliminate the displacement field.^{17} In this equation, terms proportional to a field multiplied by its time derivative describe the rate of change of energy densities in Ω. The term on the right proportional to *κ* can be expressed as the sum of such a density and a transport term utilizing the following relation, derived by elimination of ** γ, β** from Eq. 9(d):

to obtain

The dots represent terms proportional to the remaining combinations of operators ($K\u0302n\u2020K\u0302n,K\u0302nK\u0302n\u2020,K\u0302nK\u0302n$). For Eq. (34) to be valid, terms proportional to each permutation of operators on each side must balance. Rewriting the volume integral of the second term, we find

where index *i* refers to the Cartesian components of *ζ*_{n}. The remaining volume integral describes energy stored in elastic distortions of the lattice, while the surface integral describes energy transported in finite-wavevector phonon modes. Bringing the surface integral over to the left, write a composite Poynting-like vector for the Hopfield fields,

where we collected the components proportional to $K\u0302n\u2020K\u0302n\u2020$ and expanded the electromagnetic Poynting vector utilizing Eqs. 32(b) and 32(e) as

whereas for the right-hand side, we collect terms proportional to $K\u0302n\u2020K\u0302n\u2020$. Continuity of *σ*_{n} requires that tangential components of *θ*_{n}, *σ*_{n} are continuous across material interfaces. As can be seen from inspection of Eqs. (32a) and (32b), these are the standard Maxwell boundary conditions on the electromagnetic fields **E, H**. The second term fixes the additional boundary conditions for the nonlocal problem, requiring that *ζ*_{n} and the normal component of $\tau \u0304n$ are continuous. The former from Eq. (32c) enforces continuity of the lattice ionic displacement, while the latter constrains the parallel and perpendicular components of the stress at the interface, both in agreement with classical theories of polar nonlocality.^{17} Physically, the additional nonlocal boundary conditions are required to describe the additional phonon fields in the dielectric and to weight them with respect to the photonic fields. They determine the degree of nonlocal mixing.

### D. LTP quantization

Now, we have all the tools necessary to find the Hopfield fields for a given inhomogeneous system. There is still one step remaining as for operators $K\u0302n$ to describe fundamental bosonic excitations of the nonlocal system, they also need to be quantized as bosonic fields,^{25}

where the sign function yields $sgn\omega n=\omega n/|\omega n|$. Utilizing the definition of the polaritonic operators in terms of the quantum fields, we derive the result

To proceed, it is useful to introduce a further layer of abstraction. Although an LTP is comprised of distinct longitudinal and transverse parts, this is not captured by the Hopfield fields |**Ψ**_{n}⟩ = |*α*_{n}, *β*_{n}, *γ*_{n}, *ζ*_{n}⟩, which describe the fully coupled fields of the polariton. We can obtain a more physically meaningful result by segregating the fields into longitudinal and transverse components whose electric Hopfield fields satisfy $\u2207\u22c5\theta nT=0$ and $\u2207\xd7\theta nL=0$, respectively. We can then quantize these components separately. For transverse fields utilizing the transverse dielectric function defined in Eq. (26), we find

Note that utilizing

we can re-write this as in terms of both the electric and magnetic Hopfield fields as

which in the local limit *k* → 0 is the transverse field energy density in a dispersive dielectric.^{32}

For longitudinal fields, the Hopfield field $\beta nT$, which gives the magnetic field contribution to the polariton, goes to zero, so we can write

where we assumed an isotropic lattice. Note that for a pure longitudinal excitation, whose frequency satisfies $\u03f5Lk,\omega =0$, this can be simplified to

where we defined the Fröhlich coupling constant utilizing the Lydanne–Sachs–Teller relation in Eq. (2),

It is important to make clear that in these quantization relations, the frequencies are those of the hybrid LTP, determined by application of the full boundary conditions not the individual longitudinal or transverse fields.^{16} Finally, having quantized these excitations, we introduce additional polariton-like Hopfield fields describing the longitudinal (transverse) composition of the *n*th LTP mode L_{i,n} (T_{j,n}), which satisfy the normalization

where the index *i* rules over all the longitudinal components and *j* over all transverse components of the polariton. Then, if we expand the index of the Hopfield fields to consider both the LTP branch *n* and the contributing mode *ij* as $\alpha i,nL,\alpha j,nT$ satisfy Eq. (41), the full polariton field given by

is also quantized.

## III. PURCELL ENHANCED EMISSION NEAR AN EPSILON-NEAR-ZERO WAVEGUIDE

Now, we apply our theory to a physical system. We focus on the system of an isotropic polar dielectric film of thickness *d* sandwiched between positive dielectric cladding layers. For simplicity, we consider the case where only nonlocal longitudinal excitations are present in the system, a good approximation for systems driven near the LO phonon frequency where the TO phonon dispersion does not extend. The system under study is shown in Fig. 1(a). The film considered to occupy −*d* < *z* < 0 has a dispersive dielectric function *ϵ*, while the cladding regions defined as *z* < −*d*, *z* > 0 have frequency independent dielectric constant *ϵ*_{C}.

In the thick film limit, this kind of symmetric planar waveguide supports two spectrally degenerate SPhP modes, with the standard dispersion for the evanescent waves on a planar bilayer. When the film thickness *d* is less than the skin depth of the modes, the SPhPs on each interface hybridize into symmetric and antisymmetric superpositions, lifting the spectral degeneracy. As a result of repulsion between like charges on each side of the film, the symmetric mode blue shifts toward the LO phonon frequency, which marks the upper edge of the Reststrahlen region, in which evanescent modes are supported. In the thin film limit, this mode sits close to the LO phonon frequency where the polar film’s local dielectric function goes to zero. For this reason, it is termed epsilon-near-zero excitation. These modes are especially interesting in nanophotonics as they allow for strong enhancement of the out-of-plane field in the nanolayer with applications in sensing and mid-infrared nano-waveguiding.^{33,34} These systems are especially interesting in optical nonlocality, where classical studies have demonstrated hybridization between them and localized longitudinal phonons in the nanolayer.^{17,30}

### A. Modes of the system

We construct the eigenmode of the epsilon-near-zero mode to satisfy the Maxwell boundary condition on the magnetic Hopfield field ** β** as the LO field has no associated

**field. The eigenmode is a sum of SPhPs localized at each film/cladding interface. We denote the wavevector in the plane of the waveguide as $k=kxx\u20d7+kyy\u20d7$, and for brevity, write the 3D wavevector in the form $k,kz$, where**

*β**k*

_{z}is the out-of-plane component, and the position $r,z$. The electric Hopfield field

**associated with the SPhP localized at**

*θ**z*= −

*d*is given by

while that located at *z* = 0 has

In the above equations, *α* (*α*_{C}) is the out-of-plane wavevector in the film (cladding),

The total symmetric field electric field coefficient in the waveguide is the linear superposition *θ*^{T} = *θ*_{u,k} + *θ*_{l,k}. This collective excitation is quantized in Appendix B and can be written in the following form:

where $AkT$ is the quantization factor, calculated by independent quantization of the transverse field, and $ukz$ are unit vectors describing the out-of-plane mode functions of the coupled excitations.

As the cladding layers are phonon-inactive, the thin film also acts as a closed cavity for LO phonons, supporting modes with frequencies

where *n* is an integer defining the mode number, *ξ*_{n} = *nπ*/*d* is the quantized out-of-plane wavevector of the phonon, and **k** is the wavevector in the plane of the film. If we consider only symmetric modes in the waveguide *n* ∈ odd, the electric Hopfield field ** θ** of the phonons can be written as

where $Bk,nL$ is a constant determined from the quantization Eq. (45) derived in Appendix C and given by

where S is the quantization surface and $Lk,nL$ is a quantization length for the mode also defined in Appendix C, which accounts for the frequency shift of the mode from the zone-center.

Both sets of excitations are hybridized through the Maxwell boundary condition on the tangential component of ** θ** and the additional boundary conditions discussed in Sec. II C. In keeping with classical studies of polar nonlocality, we utilize the additional boundary condition on the normal component of

**, which can be recast as a constraint on the normal component of**

*ζ**ϵ*

_{∞}

**. Application of these boundary conditions yields the dispersion relation**

*θ*which can be manipulated into the form^{30}

by expanding the hyperbolic functions into Mittag–Leffler series and exploiting the analytic dispersion of the SPhP $\omega kT$ in the local limit. As has previously been noted, this is the dispersion relation arising from a Hamiltonian where a single photonic mode is coupled to a bath of LO phonons with coupling frequency^{30}

As discussed in Appendix D, the full field of the polariton can be constructed from the separately quantized longitudinal and transverse fields of the LTP constituents by introduction of the weighting coefficients T_{k}, L_{k}. Then, the total electric Hopfield field of the LTP is given by

where the transverse and longitudinal mode functions $\theta kT,\theta kL$ are defined and quantized individually in Appendixes B and C, respectively, and the index *n* runs over the discrete longitudinal phonon branches.

### B. Nonlocal Purcell enhancement

It is well known that emitters placed near to interfaces supporting plasmon polaritons can exhibit enhanced spontaneous emission due to the strong local field enhancement.^{35} The enhancement can be described utilizing the Purcell factor, which is the ratio of the emission rate of a dipole near to a photonic structure to that of the same dipole embedded in a homogeneous medium. It describes the enhancement in spontaneous emission induced by the resonator and is an especially interesting quantity in systems containing nanoscale features. In these systems, predictions of large Purcell enhancements obtained with local optical theories can be modified when nonlocality is taken into account. Studies of plasmonic systems have demonstrated that the effect of nonlocality is to strongly diminish the achievable Purcell factor as a consequence of nonlocal quenching and transfer of energy out of the electromagnetic field into the kinetic energy of free electrons,^{36,37} possibly altering the emission regime.^{38,39} Nonlocality is expected to have a qualitatively different effect in polar systems, mostly because of the propagative nature of LO phonon modes within the Reststrahlen region. This allows energy transferred from an emitter into the system’s matter degrees of freedom to be coherently recycled as part of the LTP mode. Additionally, as the Purcell factor can be particularly strong in regions where the polariton group velocity is low, LTP systems could allow for enhanced spontaneous emission around their localized LO phonon frequencies, leading to large Purcell factors, which can be tuned across the Reststrahlen region.

To demonstrate this, we apply our theory to the calculation of the Purcell enhancement for the system studied in Sec. III A. An example spectrum, calculated from the dispersion relation found in Appendix D, is shown in Fig. 1 for a 10 nm 3C–SiC film embedded between two high-index cladding regions. The epsilon-near-zero SPhP dispersion red-shifts with the increasing wavevector^{34} undergoing avoided crossings^{40} with each coupled LO phonon branch.^{30} As each polariton branch heads toward its asymptotic frequency, its group velocity tends to zero.

In general, we can write the quantized transverse electric field operator of the LTP as

where h.c. refers to the Hermitian conjugate, **k** is the wavevector in the waveguide plane, T_{k} is the transverse weighting derived in Appendix D, C_{k} is the quantization constant derived in Appendix B, **u**_{k} is the out-of-plane mode function derived in Appendix D, and *ω*_{k} is the LTP frequency derived from the nonlocal boundary conditions on the Hopfield equations in Appendix D.

We consider an emitter embedded in the upper positive dielectric halfspace, separated from the nonlocal film by a distance *h*. To derive meaningful results, we need to include broadening for the LTP modes. Although this could have been achieved extending the nonlocal theory in Sec. II through coupling the coherent system to a continuum of bath modes and Fano diagonalizing the polariton system,^{25,41} this is beyond the scope of this paper and is, nevertheless, not crucial due to the large quality factors of SPhPs and LTPs. Here, we instead include losses phenomenologically through a Lorentzian density of states for the transition $X$. The emission rate is given by

in which the transition matrix element is given by

where $d\u0302$ is the dipole moment; states are labeled |*i*, *j*⟩, where *i* labels the level of the dipole and *j* labels the LTP system; and *n*_{k} is the LTP population. Using Eq. (62), defining $d\u030212=\u27e81|d\u0302|2\u27e9$, and using the definition of the electric field [Eq. (62)], we can write

Here, the *n*_{k} term accounts for stimulated emission and the unity term accounts for spontaneous emission. Ignoring the stimulated emission, we find that the spontaneous emission rate for a dipole with transition frequency *ω*_{0} = *E*_{2} − *E*_{1} is given by

Evaluating the inner integral, we obtain

where we defined unitless projections of the dipole moment operator *d*_{12,∥} and *d*_{12,z} using $d12,\u2225\u22c5k\u0302=d12d12,\u2225\u2061cos\u2061\theta $ and defined $d12,z=d12d12\u22c5z\u0302$, $d12,\u2225=d12\u2212d12z\u0302$, and the effective length

Dividing by the spontaneous emission rate of a dipole embedded in an infinite medium,

we recover the Purcell factor

To demonstrate this result, we calculate the Purcell factor for the thin-nonlocal film system for a few dipole-film separations. We have previously demonstrated that decreasing the thickness of a nonlocal film leads to a spectral spreading of the longitudinal mode frequencies away from the zone-center longitudinal phonon, in conjunction with an increase in the LTP coupling strength as a result of the enhanced electromagnetic fields at the film interfaces.^{16} The results are shown in Fig. 2, where rows demarcate the different film thicknesses (10, 5, and 2 nm) and columns are for film–emitter separations (50, 25, and 10 nm). Emission rates for dipoles polarized perpendicular to the film are shown by the dashed lines, while those parallel to the film are indicated by the solid lines. In the local case (red curves), the Purcell enhancement exhibits a single peak for each combination of film thickness and separation; this represents that the frequency where the photonic modes overlap with the emitter is strongest. The result is a peak because the field confinement increases with the increasing in-plane wavevector (lower frequency) until at a large enough in-plane wavevector, the emitter no longer sits in the near-field of the mode. In the nonlocal case (green curves), the Purcell enhancement breaks into distinct narrow peaks; these result from emission into the different polaritonic branches of the system. For the 10 nm film, there are many branches as the LO modes of the nanolayer are spectrally close. As the emitter–film separation decreases, the number of peaks decreases and nonlocal emission outpaces the local result. This is because emitters very close to the nanolayer can emit into the low group-velocity flatband polariton regions where the density of final states is large. These flatband regions correspond to the peaks for 10 nm separation. Similar physics is observed for the 5 nm film in the middle row. There are less peaks in the nonlocal emission spectrum because the LO modes of the layer are more spaced out. For the 2 nm film, the LO modes are very far apart in frequency, and only one mode remains on the plot. As the emitter moves closer to the film, the nonlocal peak red-shifts toward the flatband polariton region. This is a result of an enhanced density of states at the flatband polariton region; however, there the fields become increasingly evanescent, only coupling to nearby emitters. This, however, happens more slowly as the LTP modes remain dispersive to a larger in-plane wavevector for thinner layers.^{19}

This plot shows that although the LTP is less photonic in nature ($|Tk\u2225|2<1$) than the local epsilon-near-zero (ENZ) excitation, it is able to enhance the Purcell factor. This is because near to the avoided crossing regions of the dispersion and at large in-plane wavevectors, the LTP group velocity drops, enhancing the density of LTP states and leading to a peak in emission. In the regions between polariton branches, emission is suppressed, leading to multi-peak Purcell spectra. In thinner films where coupling between LO phonons and SPhPs is enhanced, this results in the narrowband emission enhancement.

## IV. CONCLUSION

In this work, we have presented a full quantum theory of longitudinal-transverse polaritons, directly derived from consideration of the free fields in an inhomogeneous nonlocal medium. We derived the equations of motion for the Hopfield fields of the polariton, demonstrating them to be equivalent to the nonlocal macroscopic Maxwell equations shown previously. This point is important because it means that these equations can be solved by standard numerical means, allowing our theory to be integrated with the finite element or finite difference time domain solvers typically used to study SPhP resonators.^{25} We also derived equations which allow for the quantization of LTP modes. Finally, we derived the nonlocal Purcell enhancement for a simple one-dimensional waveguide, demonstrating enhanced emission for dipoles close to the nonlocal layer. This is an example of how the tunability of the LTP system could be exploited for enhanced light–matter interaction. The quantum theory developed in this paper is important for the study of novel and emergent phenomena in LTP systems. In a polar system, nonlinear electron–electron scattering can be mediated by the charged ions comprising the lattice, leading to emission of LO phonons. If these LO phonons are part of an LTP, this could permit for direct electrical excitation or detection of SPhP modes, potentially underpinning a novel generation of mid-infrared optoelectronic devices.^{42} A quantum theory such as that presented here is necessary to calculate this emission rate. Similarly, it can also be applied to descriptions of nonlinear LTP–LTP scattering mediated by anharmonic interactions between the underlying optical phonons, something which could allow for the formation of an LTP condensate analogous to the exciton–polariton condensates observed in the previous decade in the visible spectral region.^{43}

## ACKNOWLEDGMENTS

S.D.L. was supported by a Royal Society Research Fellowship and the Philip Leverhulme Prize. The authors acknowledge support from the Royal Society under Grant No. RGF∖EA∖181001 and the Leverhulme Trust under Grant No. RPG-2019-174.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: NONLOCAL QUANTIZATION

In this section, we derive the nonlocal equations of motion for the Hopfield fields, given in Eqs. (9a)–(9d). To achieve this, we require the standard commutation relations between the canonical variables of the local Hamiltonian [Eq. (1)], which are given by^{25}

where the differential operator ∇′ operates on **r**′. These equations allow us to calculate the commutator of the polaritonic operator [Eq. (7)] with the local Hamiltonian [Eq. (1)]. As the commutation relation [Eq. (8)] is lengthy, we proceed term by term, calculating the commutator of each polariton component with the Hamiltonian. Utilizing the commutator between the electromagnetic fields [Eq. (A1a)], we can find the following contributions to the polariton equation of motion [Eq. (8)]:

and furthermore, through Eq. (A1b), we can find

If we now collect terms proportional to the remaining operators to form Eq. (8), we find four local equations of motion,

These equations can be solved by the introduction of the novel electric Hopfield field ** θ**,

and eliminating ** ζ** from Eq. (A4a),

Finally, on inversion, this leads to

where we recognized the local dielectric function of the polar lattice in the lossless limit,

We can extend this theory to the nonlocal case, considering also commutation between the polaritonic operator $K\u0302$ and the nonlocal component of the Hamiltonian [Eq. (6)]. To do this, we need to derive additional commutation relationships between the ionic displacement field and its spatial derivatives, which enter into the nonlocal component of the Hamiltonian. In the general case, we can the commutator between momentum and the first derivative of the ionic position,

where we noted that the derivative operator $\u2202m\u2032$ operates on **r**′ rather than **r** and, therefore, commutes with functions of **r**. We also need to calculate commutation relations involving the stresses $S\u0302ij$. Writing the nonlocal Hamiltonian density entering Eq. (6) fully,

we commute

in which we utilized the commutation relation [Eq. (A10)]

and the identity

Note that in Eq. (A12), our underlying assumption of piecewise homogeneity prevented any material parameters contributing to the right-hand side.

Unfortunately, the derivatives in Eq. (A12) are on the field operators X rather than the Hopfield field *ζ*, which precludes us from calculating the equation of motion. We can overcome this problem by collecting terms proportional to *ζ*_{i} to find the integrated quantity

where we assumed that the Hopfield field and its derivatives vanish at the boundaries of the system. Similar results are derivable for the other terms in Eq. (A12). Finally, if we define the elements of the material stress tensor ** τ** as

we can re-write the result of Eq. (A12),

where we noted the major symmetry of the stiffness tensor $C\u0304ijkl=C\u0304klij$ and the minor symmetries of the stiffness tensor $C\u0304ijkl=C\u0304jikl=C\u0304ijlk$ and, in the final step, utilized the tensor divergence

This modifies the fourth polaritonic equation of motion to

which is that utilized in the main body of this paper.

### APPENDIX B: QUANTIZATION OF THE ENZ

The epsilon-near-zero (ENZ) mode is a linear superposition of SPhPs at each interface of the thin film studied in the main body of this paper. The electric Hopfield field of the two SPhPs comprising the ENZ is defined in the main body of this paper as

while the total transverse field is $\theta kT=\theta l,k+\theta u,k$. The magnetic Hopfield field $\beta =i\omega \u2207\xd7\theta $ can be written as

which as mentioned in the main body is continuous at material boundaries, while the tangential electric Hopfield field ** θ** is not. In the local case, we would apply continuity of $z\u20d7\xd7\theta $ at the film interfaces, recovering the dispersion relation for odd-parity modes in the planar waveguide,

As we are interested in LTP modes, we instead derive the quantization condition in the general case. We start by integrating the fields in the barriers,

where we expanded in powers of *d* and retained terms to linear order. Here, S is the in-plane quantization surface. Finally, we integrate over the film

where again terms up to order *d* were retained. In the lossless cladding, $I\alpha C=0,R\alpha C=\alpha C$ and similarly for the film, so

and if we define the quantization length

we can write, for example, the quantized electric field coefficients

where we defined

Note that in the main body of this paper, we denote the quantized transverse field by

where $uz$ is the out-of-plane mode function of the excitation.

### APPENDIX C: NONLOCAL QUANTIZATION

We consider a single longitudinal phonon branch in the central layer. The electric potential associated with LO phonons localized in the thin film −*d* < *z* < 0 is given by

where **k** is the in-plane wavevector and *ξ* is the out-of-plane wavevector. This relates to the electric field through **E** = −∇*ϕ*. If the cladding layers act as hard barriers for the phonon, the out-of-plane wavevector of the mode will be quantized, as for photons in a Fabry–Pèrot cavity, characterized by integer index *n* and out-of-plane wavevector *ξ*_{n} = *nπ*/*d*. We hypothesize an electric Hopfield field of the form

where B_{k,n} is a normalization constant to be determined from the quantization conditions. As the phonon is entirely localized within region 2, it can be quantized considering only the energy in the film, yielding

Now, we can utilize that for real *ξ*, valid for frequencies *ω* < *ω*_{L}, and for odd values of *n*,

This result allows us to identify a quantization factor similarly to that for the pure longitudinal modes on the LO phonon dispersion relation,

where we defined the effective mode length as

Here, the first term collects the wavevector dependent components from the integral, and the latter accounts for the modal drift away from the zone center *ω*_{L}.

Finally, we define a quantization constant suitable for all non-quantized (continuous) values of *ξ*, given by

where

We derived this quantity by carrying out the same calculation but not assuming that the integral over the film simplifies as in Eq. (C4). This is necessary to calculate the longitudinal and transverse Hopfield fields in Appendix D.

### APPENDIX D: NONLOCAL DISPERSION

To uniquely determine the LTP dispersion relation, we need to consider the longitudinal and transverse Hopfield fields together and apply the boundary conditions on them, as discussed in Sec. II C. In Appendix B, we constructed the ENZ field, so as its magnetic field coefficient ** β** is continuous. As the LO modes constructed in Appendix C are curl-free, they have no associated

**. We, therefore, need to construct the fields to satisfy the second Maxwell boundary condition on the parallel electric field coefficient $\theta \xd7z\u20d7$ and an additional nonlocal boundary condition. As discussed in the main body, the appropriate choice of additional condition is on $\u03f5\u221e\theta \u22c5z\u20d7$. The boundary condition on $\theta \xd7z\u20d7$ can be written as**

*β*where L_{k}, T_{k} are expansion coefficients weighting the longitudinal and transverse components of the LTP. Note that in the transverse limit L → 0, this gives

which is the standard dispersion relation for the ENZ mode in a symmetric trilayer waveguide.^{34} The boundary condition on $\u03f5\u221e\theta \u22c5z\u20d7$ gives

Dividing through and defining *ξ*′ = *iξ*, we find the standard dispersion for symmetric modes in a trilayer nonlocal waveguide,^{16}

which is the form given in the main body. This equation uniquely determines the relationship between the modes frequency and the wavevector. To find the transverse electric field, we also need to determine the L_{k}, T_{k} longitudinal-transverse expansion coefficients. This can be done using either equation, for example,

and the normalization condition

derived in the main body under the assumption that the longitudinal and transverse fields are quantized. This gives

which yields the expansion coefficients for the transverse field of the LTP, for example,

where *θ*_{l,k}, *θ*_{u,k} are the quantized transverse fields given by Eqs. (B10)]–(B13)] in Appendix B.

## REFERENCES

^{2}) scattering on a β-SiC surface