Polaritonic quantization in nonlocal polar materials

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,¹ with applications in mid-infrared sensing,² 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.¹⁰ 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.¹⁵ 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.¹⁶ 

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.²⁰ 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.²⁹ 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³⁰ 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.¹⁷ 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.

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̂(Ĥ)$ is the electric displacement (magnetic) field operator, $X̂$ is the ionic displacement operator, $Q̂$ 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.²⁵ 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.³¹ 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̄ijkl$ is an effective elasticity tensor satisfying the symmetry conditions¹⁸ 

and $Ŝij$ 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.¹⁷ The full nonlocal Hamiltonian is a sum of the two Hamiltonian densities,

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,²⁵ 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̂$ is an eigenmode of $Ĥ$⁠, 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 $τ̄$ 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.¹⁷ For simplicity, we focus on isotropic crystal structures, for which the components of $C̄ijkl$ 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.¹⁷ These boundary conditions are derived in Sec. II C.

We can solve the equations of motion [Eqs. (9a)–(9d)] in terms of a set of eigenmodes,

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 $Ŷ$ 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 $Ê$⁠, and polarization field operator $P̂$ can be written as

where in the final relation, we used the relation between the material polarization density operator $P̂=κX̂$ 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 β.

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.¹⁷ 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⃗$⁠, given by

where we utilized the standard constitutive relation for a nonlocal dielectric to eliminate the displacement field.¹⁷ 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̂n†K̂n,K̂nK̂n†,K̂nK̂n$⁠). 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̂n†K̂n†$ 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̂n†K̂n†$⁠. 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 $τ̄n$ 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.¹⁷ 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.

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̂n$ to describe fundamental bosonic excitations of the nonlocal system, they also need to be quantized as bosonic fields,²⁵ 

where the sign function yields $sgnωn=ωn/|ω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 $∇⋅θnT=0$ and $∇×θ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.³² 

For longitudinal fields, the Hopfield field $β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 $ϵLk,ω=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.¹⁶ Finally, having quantized these excitations, we introduce additional polariton-like Hopfield fields describing the longitudinal (transverse) composition of the nth 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 $αi,nL,αj,nT$ satisfy Eq. (41), the full polariton field given by

is also quantized.

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

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⃗+kyy⃗$⁠, 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³⁰ 

by expanding the hyperbolic functions into Mittag–Leffler series and exploiting the analytic dispersion of the SPhP $ω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³⁰ 

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 $θkT,θkL$ are defined and quantized individually in  Appendixes B and  C, respectively, and the index n runs over the discrete longitudinal phonon branches.

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.³⁵ 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³⁴ undergoing avoided crossings⁴⁰ with each coupled LO phonon branch.³⁰ 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̂$ 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̂12=⟨1|d̂|2⟩$⁠, 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 ω₀ = E₂ − E₁ 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,∥⋅k̂=d12d12,∥⁡cos⁡θ$ and defined $d12,z=d12d12⋅ẑ$⁠, $d12,∥=d12−d12ẑ$⁠, 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.¹⁶ 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.¹⁹ 

This plot shows that although the LTP is less photonic in nature (⁠$|Tk∥|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.

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.²⁵ 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.⁴² 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.⁴³ 

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.

The authors have no conflicts to disclose.

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

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²⁵ 

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 θ,

The equation of motion for θ can be derived combining Eqs. (A4c) and (A4d) to find

and eliminating ζ from Eq. (A4a),

Finally, on inversion, this leads to