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.

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 layers17,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–polariton26–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 waveguide30 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.

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

(1)

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.25 The oscillator strength κ characterizes the width of the Reststrahlen that relates to the phonon frequencies through

(2)

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,

(3)

where C̄ijkl is an effective elasticity tensor satisfying the symmetry conditions18 

(4)

and Ŝij are scalar stress operators, defined as

(5)

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,

(6)

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,

(7)

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,

(8)

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,

(9a)
(9b)
(9c)
(9d)

where τ̄ is the stress tensor of the polar lattice, related to the effective elasticity tensor through

(10)

and we define the composite field

(11)

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̄ijkl are given by

(12)
(13)
(14)

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 works17,30 through

(15)
(16)
(17)

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

(18)
(19)
(20)
(21)
(22)
(23)

where Sij are scalar stresses of the Hopfield field ζ,

(24)

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,

(25)

where ϵT is the transverse dielectric function of the lattice, given by

(26)

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

(27)

where ϵL is the longitudinal dielectric function of the lattice, given by

(28)

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.

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

(29)

where the index n can be either discrete or continuous depending on the mode and system under study with associated eigenfrequencies ωn and wavevectors kn. Then, an arbitrary field operator Ŷ can be written as a sum over the polaritonic operators as

(30)

in which the expansion coefficients are given by

(31)

The electric displacement field, magnetic field, ionic displacement, electric field operator Ê, and polarization field operator P̂ can be written as

(32a)
(32b)
(32c)
(32d)
(32e)

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 β:

(33)

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.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, given by

(34)

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):

(35)

to obtain

(36)

The dots represent terms proportional to the remaining combinations of operators (K̂nK̂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

(37)

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,

(38)

where we collected the components proportional to K̂nK̂n and expanded the electromagnetic Poynting vector utilizing Eqs. 32(b) and 32(e) as

(39)

whereas for the right-hand side, we collect terms proportional to K̂nK̂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.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.

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

(40)

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

(41)

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

(42)

Note that utilizing

(43)

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

(44)

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

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

(45)

where we assumed an isotropic lattice. Note that for a pure longitudinal excitation, whose frequency satisfies ϵLk,ω=0, this can be simplified to

(46)

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

(47)

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 nth LTP mode Li,n (Tj,n), which satisfy the normalization

(48)

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

(49)

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.

FIG. 1.

(a) A sketch of the thin film–emitter system. A negative dielectric film with transverse dielectric function ϵ2 < 0 sits between −d < z < 0. Symmetric positive dielectric cladding with dielectric constant ϵ1 occupies z < −d and z > 0. The film supports a SPhP excitation (red) and a discrete spectrum of localized LO phonon modes (blue), which are hybridized through the film boundaries with coupling frequency Ω. We consider an emitter embedded in the upper, a positive dielectric, and halfspace, separated from the film by a height h. (b) The inverse of the polariton dispersion G as defined in  Appendix D for a 10 nm 3C–SiC film between two high-index (nC = 2.6) cladding regions. The dark regions illustrate the modal frequencies. The horizontal dashed green lines indicate the LO phonon branches, and the black dashed line indicates the epsilon-near-zero dispersion. Polariton frequencies are shown by the blue solid lines.

FIG. 1.

(a) A sketch of the thin film–emitter system. A negative dielectric film with transverse dielectric function ϵ2 < 0 sits between −d < z < 0. Symmetric positive dielectric cladding with dielectric constant ϵ1 occupies z < −d and z > 0. The film supports a SPhP excitation (red) and a discrete spectrum of localized LO phonon modes (blue), which are hybridized through the film boundaries with coupling frequency Ω. We consider an emitter embedded in the upper, a positive dielectric, and halfspace, separated from the film by a height h. (b) The inverse of the polariton dispersion G as defined in  Appendix D for a 10 nm 3C–SiC film between two high-index (nC = 2.6) cladding regions. The dark regions illustrate the modal frequencies. The horizontal dashed green lines indicate the LO phonon branches, and the black dashed line indicates the epsilon-near-zero dispersion. Polariton frequencies are shown by the blue solid lines.

Close modal

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 kz 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

(50)

while that located at z = 0 has

(51)

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

(52)
(53)

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:

(54)

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

(55)

where n is an integer defining the mode number, ξ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

(56)

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

(57)

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

(58)

which can be manipulated into the form30 

(59)

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 frequency30 

(60)

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 Tk, Lk. Then, the total electric Hopfield field of the LTP is given by

(61)

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.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 wavevector34 undergoing avoided crossings40 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

(62)

where h.c. refers to the Hermitian conjugate, k is the wavevector in the waveguide plane, Tk is the transverse weighting derived in  Appendix D, Ck is the quantization constant derived in  Appendix B, uk 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

(63)

in which the transition matrix element is given by

(64)

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 nk 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

(65)

Here, the nk 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 = E2E1 is given by

(66)

Evaluating the inner integral, we obtain

(67)

where we defined unitless projections of the dipole moment operator d12,∥ and d12,z using d12,k̂=d12d12,cosθ and defined d12,z=d12d12ẑ, d12,=d12d12ẑ, and the effective length

(68)

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

(69)

we recover the Purcell factor

(70)

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 

FIG. 2.

Calculated Purcell enhancement for 10 nm [(a)–(c)], 5 nm [(d)–(f)], and 2 nm [(g)–(i)] 3C–SiC films for a dipole located 50 nm [(a), (d), and (h)], 25 nm [(b), (e), and (i)], and 10 nm [(c), (f), and (j)] above the film. The green curves indicate the local enhancement, while the red curves show the nonlocal. Enhancement for dipoles polarized parallel and perpendicular to the interface is indicated by the solid and dashed lines, respectively.

FIG. 2.

Calculated Purcell enhancement for 10 nm [(a)–(c)], 5 nm [(d)–(f)], and 2 nm [(g)–(i)] 3C–SiC films for a dipole located 50 nm [(a), (d), and (h)], 25 nm [(b), (e), and (i)], and 10 nm [(c), (f), and (j)] above the film. The green curves indicate the local enhancement, while the red curves show the nonlocal. Enhancement for dipoles polarized parallel and perpendicular to the interface is indicated by the solid and dashed lines, respectively.

Close modal

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

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 by25 

(A1a)
(A1b)

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)]:

(A2a)
(A2b)
(A2c)

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

(A3a)
(A3b)
(A3c)

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

(A4a)
(A4b)
(A4c)
(A4d)

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

(A5)

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

(A6)

and eliminating ζ from Eq. (A4a),

(A7)

Finally, on inversion, this leads to

(A8)

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

(A9)

We can extend this theory to the nonlocal case, considering also commutation between the polaritonic operator K̂ 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,

(A10)

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

(A11)

we commute

(A12)

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

(A13)

and the identity

(A14)

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

(A15)

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

(A16)

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

(A17)

where we noted the major symmetry of the stiffness tensor C̄ijkl=C̄klij and the minor symmetries of the stiffness tensor C̄ijkl=C̄jikl=C̄ijlk and, in the final step, utilized the tensor divergence

(A18)

This modifies the fourth polaritonic equation of motion to

(A19)

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

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

(B1)
(B2)

while the total transverse field is θkT=θl,k+θu,k. The magnetic Hopfield field β=iω×θ can be written as

(B3)
(B4)

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×θ at the film interfaces, recovering the dispersion relation for odd-parity modes in the planar waveguide,

(B5)

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,

(B6)

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

(B7)

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

(B8)

and if we define the quantization length

(B9)

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

(B10)
(B11)

where we defined

(B12)

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

(B13)

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

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

(C1)

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 = /d. We hypothesize an electric Hopfield field of the form

(C2)

where Bk,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

(C3)

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

(C4)

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

(C5)

where we defined the effective mode length as

(C6)

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

(C7)

where

(C8)

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.

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 θ×z and an additional nonlocal boundary condition. As discussed in the main body, the appropriate choice of additional condition is on ϵθz. The boundary condition on θ×z can be written as

(D1)

where Lk, Tk are expansion coefficients weighting the longitudinal and transverse components of the LTP. Note that in the transverse limit L → 0, this gives

(D2)

which is the standard dispersion relation for the ENZ mode in a symmetric trilayer waveguide.34 The boundary condition on ϵθz gives

(D3)

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

(D4)

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 Lk, Tk longitudinal-transverse expansion coefficients. This can be done using either equation, for example,

(D5)

and the normalization condition

(D6)

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

(D7)

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

(D8)

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

1.
J. D.
Caldwell
 et al., “
Low-loss, extreme subdiffraction photon confinement via silicon carbide localized surface phonon polariton resonators
,”
Nano Lett.
13
,
3690
3697
(
2013
).
2.
R.
Berte
 et al., “
Sub-nanometer thin oxide film sensing with localized surface phonon polaritons
,”
ACS Photonics
5
,
2807
2815
(
2018
).
3.
C. R.
Gubbin
and
S.
De Liberato
, “
Theory of four-wave-mixing in phonon polaritons
,”
ACS Photonics
5
,
284
288
(
2017
).
4.
I.
Razdolski
 et al., “
Second harmonic generation from strongly coupled localized and propagating phonon-polariton modes
,”
Phys. Rev. B
98
,
125425
(
2018
).
5.
S.
Kitade
,
A.
Yamada
,
I.
Morichika
,
K.
Yabana
, and
S.
Ashihara
, “
Nonlinear shift in phonon-polariton dispersion on a SiC surface
,”
ACS Photonics
8
,
152
157
(
2021
).
6.
J.-J.
Greffet
 et al., “
Coherent emission of light by thermal sources
,”
Nature
416
,
61
64
(
2002
).
7.
J. A.
Schuller
,
T.
Taubner
, and
M. L.
Brongersma
, “
Optical antenna thermal emitters
,”
Nat. Photonics
3
,
658
661
(
2009
).
8.
C.
Arnold
 et al., “
Coherent thermal infrared emission by two-dimensional silicon carbide gratings
,”
Phys. Rev. B
86
,
035316
(
2012
).
9.
G.
Lu
 et al., “
Engineering the spectral and spatial dispersion of thermal emission via polariton–phonon strong coupling
,”
Nano Lett.
21
,
1831
1838
(
2021
).
10.
C.
Ciracì
,
J. B.
Pendry
, and
D. R.
Smith
, “
Hydrodynamic model for plasmonics: A macroscopic approach to a microscopic problem
,”
ChemPhysChem
14
,
1109
1116
(
2013
).
11.
C.
Ciracì
 et al., “
Probing the ultimate limits of plasmonic enhancement
,”
Science
337
,
1072
1074
(
2012
).
12.
A. I.
Fernández-Domínguez
,
A.
Wiener
,
F. J.
García-Vidal
,
S. A.
Maier
, and
J. B.
Pendry
, “
Transformation-optics description of nonlocal effects in plasmonic nanostructures
,”
Phys. Rev. Lett.
108
,
106802
(
2012
).
13.
N. A.
Mortensen
,
S.
Raza
,
M.
Wubs
,
T.
Søndergaard
, and
S. I.
Bozhevolnyi
, “
A generalized non-local optical response theory for plasmonic nanostructures
,”
Nat. Commun.
5
,
3809
(
2014
).
14.
Y.
Luo
,
A. I.
Fernandez-Dominguez
,
A.
Wiener
,
S. A.
Maier
, and
J. B.
Pendry
, “
Surface plasmons and nonlocality: A simple model
,”
Phys. Rev. Lett.
111
,
093901
(
2013
).
15.
C. R.
Gubbin
 et al., “
Hybrid longitudinal-transverse phonon polaritons
,”
Nat. Commun.
10
,
1682
(
2019
).
16.
C. R.
Gubbin
,
S.
De Liberato
, and
T.
Folland
, “
Perspective: Phonon polaritons for infrared optoelectronics
,”
J. Appl. Phys.
(in press) (
2022
).
17.
C. R.
Gubbin
and
S.
De Liberato
, “
Optical nonlocality in polar dielectrics
,”
Phys. Rev. X
10
,
021027
(
2020
).
18.
C. R.
Gubbin
and
S.
De Liberato
, “
Nonlocal scattering matrix description of anisotropic polar heterostructures
,”
Phys. Rev. B
102
,
235301
(
2020
).
19.
C. R.
Gubbin
and
S.
De Liberato
, “
Impact of phonon nonlocality on nanogap and nanolayer polar resonators
,”
Phys. Rev. B
102
,
201302
(
2020
).
20.
E. V.
Yakovlev
 et al., “
Direct experimental evidence of longitudinal and transverse mode hybridization and anticrossing in simple model fluids
,”
J. Phys. Chem. Lett.
11
,
1370
1376
(
2020
).
21.
D. C.
Ratchford
 et al., “
Controlling the infrared dielectric function through atomic-scale heterostructures
,”
ACS Nano
13
,
6730
6741
(
2019
).
22.
J. J.
Hopfield
, “
Theory of the contribution of excitons to the complex dielectric constant of crystals
,”
Phys. Rev.
112
,
1555
1567
(
1958
).
23.
F.
Alpeggiani
and
L. C.
Andreani
, “
Quantum theory of surface plasmon polaritons: Planar and spherical geometries
,”
Plasmonics
9
,
965
978
(
2014
).
24.
A.
Archambault
,
F. m. c.
Marquier
,
J.-J.
Greffet
, and
C.
Arnold
, “
Quantum theory of spontaneous and stimulated emission of surface plasmons
,”
Phys. Rev. B
82
,
035411
(
2010
).
25.
C. R.
Gubbin
,
S. A.
Maier
, and
S.
De Liberato
, “
Real-space Hopfield diagonalization of inhomogeneous dispersive media
,”
Phys. Rev. B
94
,
205301
(
2016
).
26.
C. R.
Gubbin
and
S.
De Liberato
, “
Theory of nonlinear polaritonics: χ(2) scattering on a β-SiC surface
,”
ACS Photonics
4
,
1381
1388
(
2017
).
27.
I.
Carusotto
and
C.
Ciuti
, “
Quantum fluids of light
,”
Rev. Mod. Phys.
85
,
299
366
(
2013
).
28.
L. B.
Tan
 et al., “
Interacting polaron-polaritons
,”
Phys. Rev. X
10
,
021011
(
2020
).
29.
D. K.
Efimkin
,
E. K.
Laird
,
J.
Levinsen
,
M. M.
Parish
, and
A. H.
MacDonald
, “
Electron-exciton interactions in the exciton-polaron problem
,”
Phys. Rev. B
103
,
075417
(
2021
).
30.
C. R.
Gubbin
, and
S.
De Liberato
, “
Quantum theory of longitudinal-transverse polaritons in nonlocal thin films
,”
Phys. Rev. Appl.
(in press) (
2022
).
31.
R. H.
Lyddane
,
R. G.
Sachs
, and
E.
Teller
, “
On the polar vibrations of alkali halides
,”
Phys. Rev.
59
,
673
676
(
1941
).
32.
R.
Ruppin
, “
Electromagnetic energy density in a dispersive and absorptive material
,”
Phys. Lett. A
299
,
309
312
(
2002
).
33.
N. C.
Passler
 et al., “
Strong coupling of epsilon-near-zero phonon polaritons in polar dielectric heterostructures
,”
Nano Lett.
18
,
4285
4292
(
2018
).
34.
S.
Campione
,
I.
Brener
, and
F.
Marquier
, “
Theory of epsilon-near-zero modes in ultrathin films
,”
Phys. Rev. B
91
,
121408
(
2015
).
35.
W. L.
Barnes
, “
Fluorescence near interfaces: The role of photonic mode density
,”
J. Mod. Opt.
45
,
661
699
(
1998
).
36.
C.
Tserkezis
,
N.
Stefanou
,
M.
Wubs
, and
N. A.
Mortensen
, “
Molecular fluorescence enhancement in plasmonic environments: Exploring the role of nonlocal effects
,”
Nanoscale
8
,
17532
17541
(
2016
).
37.
C.
Tserkezis
,
N. A.
Mortensen
, and
M.
Wubs
, “
How nonlocal damping reduces plasmon-enhanced fluorescence in ultranarrow gaps
,”
Phys. Rev. B
96
,
085413
(
2017
).
38.
R.
Jurga
,
S.
D’Agostino
,
F.
Della Sala
, and
C.
Ciracì
, “
Plasmonic nonlocal response effects on dipole decay dynamics in the weak- and strong-coupling regimes
,”
J. Phys. Chem. C
121
,
22361
22368
(
2017
).
39.
C.
Tserkezis
,
M.
Wubs
, and
N. A.
Mortensen
, “
Robustness of the Rabi splitting under nonlocal corrections in plexcitonics
,”
ACS Photonics
5
,
133
142
(
2018
).
40.
C. R.
Gubbin
,
F.
Martini
,
A.
Politi
,
S. A.
Maier
, and
S.
De Liberato
, “
Strong and coherent coupling between localized and propagating phonon polaritons
,”
Phys. Rev. Lett.
116
,
246402
(
2016
).
41.
S.
Rajabali
 et al., “
Polaritonic nonlocality in light-matter interaction
,”
Nat. Photonics
15
,
690
(
2021
).
42.
C. R.
Gubbin
and
S.
De Liberato
, “
Electrical generation of surface phonon polaritons
” (unpublished) (
2022
).
43.
K. S.
Daskalakis
,
S. A.
Maier
,
R.
Murray
, and
S.
Kéna-Cohen
, “
Nonlinear interactions in an organic polariton condensate
,”
Nat. Mater.
13
,
271
278
(
2014
).