Theory of nonlinear polaritonics: $\chi^{(2)}$ scattering on a $\beta$-SiC surface

TL;DR

This paper develops a quantum-level theoretical framework for second-order nonlinear scattering in inhomogeneous dielectric environments, with specific focus on $eta$-SiC surfaces, enabling practical calculations of nonlinear optical processes like SHG and DFG.

Contribution

It introduces a comprehensive quantum description of $ ext{chi}^{(2)}$ nonlinearities applicable to complex geometries, validated with $eta$-SiC surface experiments, and compatible with finite element methods.

Findings

Calculated second harmonic generation rate matching experimental data.

Determined difference frequency generation rates for surface phonon polaritons.

Provided a practical method for integrating nonlinear scattering calculations into existing simulation tools.

Abstract

In this article we provide a practical prescription to harness the rigorous microscopic, quantum level descriptions of light-matter systems provided by Hopfield diagonalisation for quantum description of nonlinear scattering. A general frame to describe the practically important second-order optical nonlinearities which underpin sum and difference frequency generation is developed for arbitrarily inhomogeneous dielectric environments. Specific attention is then focussed on planar systems with optical nonlinearity mediated by a polar dielectric $\beta$-SiC halfspace. In this system we calculate the rate of second harmonic generation and the result compared to recent experimental measurements. Furthermore the rate of difference frequency generation of subdiffraction surface phonon polaritons on the $\beta$-SiC halfspace by two plane waves is calculated. The developed theory is easily integrated with commercial finite element solvers, opening the way for calculation of second-order nonlinear scattering coefficients in complex geometries which lack analytical linear solutions.

Equations

$${\hat{\mathcal{H}}}_0$$

$${\hat{\mathcal{K}}}_n=\int \mathrm{d}\mathbf{r}\left[{\mathbf{\alpha }}_n\cdot \hat{\mathbf{D}}+{\mathbf{\beta }}_n\cdot \hat{\mathbf{H}}+{\mathbf{\gamma }}_n\cdot \hat{\mathbf{P}}+{\mathbf{\eta }}_n\cdot \hat{\mathbf{X}}\right],$$

$$\left[{\hat{\mathcal{K}}}_n,{\hat{\mathcal{H}}}_0\right]=\hslash {\omega }_n{\hat{\mathcal{K}}}_n,$$

$$\left[{\hat{\mathcal{K}}}_m,{\hat{\mathcal{K}}}_n^{†}\right]$$

$$\hat{\mathbf{E}}$$

$$\hat{\mathbf{X}}$$

$${\omega }_{\text{T}}^2={\omega }_{\text{L}}^2-\frac{{\kappa }^2}{\rho {ϵ}_0}.$$

$$ϵ\left(\omega ,z>0\right)$$

$$ϵ\left(\omega ,z<0\right)$$

$${ϵ}_i\left(\omega \right)k_0^2=k_{\parallel }^2+k_{i,z}^2,$$

$${\mathbf{\alpha }}_{1,\mathbf{k}}^{\text{TMv}}\left({\mathbf{r}}_{\parallel },z>0\right)$$

$${\mathbf{\alpha }}_{2,\mathbf{k}}^{\text{TMv}}\left({\mathbf{r}}_{\parallel },z<0\right)$$

$${\mathbf{p}}_{i\pm }=\frac{k_{\parallel }{\mathbf{e}}_z\mp k_{i,z}{\mathbf{e}}_{\parallel }}{k_0\sqrt{{ϵ}_i}},$$

$$N_{\mathbf{k}}^{\text{TMv}}=\sqrt{\frac{1}{\hslash {ϵ}_0\omega \mathcal{V}}},$$

$$k_{\parallel }=\frac{\omega }{c}\sqrt{\frac{{ϵ}_2(\omega )}{{ϵ}_2(\omega )+1}},$$

$${\mathbf{\alpha }}_{1,{\mathbf{k}}_{\parallel }}^{\mathrm{S}}\left({\mathbf{r}}_{\parallel },z>0\right)$$

$${\mathbf{\alpha }}_{2,{\mathbf{k}}_{\parallel }}^{\mathrm{S}}\left({\mathbf{r}}_{\parallel },z<0\right)$$

$$N_{{\mathbf{k}}_{\parallel }}^{\text{S}}$$

$$=\frac{k_0}{k_{\parallel }}\sqrt{\frac{{ϵ}_2\left(\omega \right)}{\hslash {ϵ}_0\omega {\mathcal{V}}_{\mathrm{eff}}\left(\omega \right)}},$$

$${\mathbf{\alpha }}_{2,{\mathbf{k}}_{\parallel }}^{\mathrm{S}}\left({\mathbf{r}}_{\parallel },z<0\right)=\frac{e^{i{\kappa }_{2,z}z}{e}^{i{\mathbf{k}}_{\parallel }\cdot {\mathbf{r}}_{\parallel }}}{\sqrt{\hslash {ϵ}_0\omega {\mathcal{V}}_{\mathrm{eff}}\left(\omega \right)}}\left({\mathbf{e}}_{\parallel }-\frac{k_{\parallel }}{{\kappa }_{2,z}}{\mathbf{e}}_z\right),$$

$${\mathcal{H}}_{\mathrm{NL}}^{(N)}$$

$${\chi }^{(2)}\left({\omega }_1,{\omega }_2,{\omega }_3\right)$$

$$+G_2\left[D\left({\omega }_1\right)D\left({\omega }_2\right)+D\left({\omega }_1\right)D\left({\omega }_3\right)+D\left({\omega }_2\right)D\left({\omega }_3\right)\right]$$

$$+G_{3}D\left({\omega }_1\right)D\left({\omega }_2\right)D\left({\omega }_3\right),$$

$$D\left(\omega \right)=\frac{{\omega }_{\text{TO}}^2}{{\omega }_{\text{TO}}^2-{\omega }^2},$$

$${\mathcal{H}}_{\text{NL}}^{(2)}=\int \mathrm{d}\mathbf{r}\sum _{\mathbf{n},{\mathbf{j}}_{\mathbf{n}}}{ϵ}_0{\chi }^{(2)}\left({\omega }_{n_1},{\omega }_{n_2},{\omega }_{n_3}\right){\hat{\mathrm{E}}}_{n_1}^{j_{n_1}}{\hat{\mathrm{E}}}_{n_2}^{j_{n_2}}{\hat{\mathrm{E}}}_{n_3}^{j_{n_3}},$$

$${\hat{\mathrm{E}}}_n^j=\hslash {\omega }_n\left[{\bar{\alpha }}_n^j{\hat{\mathcal{K}}}_n+{\alpha }_n^j{\hat{\mathcal{K}}}_n^{†}\right],$$

$$G_1$$

$$G_2$$

$$G_3$$

$${\Xi }^{n_1,n_2,n_3}=$$

$$\times {\mathcal{T}}^{n_1,n_2,n_3}\delta \left({\mathbf{k}}_{\parallel }^{n_1}+{\mathbf{k}}_{\parallel }^{n_2}-{\mathbf{k}}_{\parallel }^{n_3}\right)$$

$$\times \left[\mathcal{P}\left(\frac{1}{i\Delta k_z^{n_1,n_2,n_3}}\right)-\pi \delta \left(\Delta k_z^{n_1,n_2,n_3}\right)\right],$$

$$\Delta k_z^{n_1,n_2,n_3}=k_{2,z}\left({\mathbf{k}}_{\parallel }^{n_1},k_{1,z}^{n_1}\right)+k_{2,z}\left({\mathbf{k}}_{\parallel }^{n_2},k_{1,z}^{n_2}\right)+k_{2,z}^{n_3},$$

$${\mathcal{T}}^{n_1,n_2,n_3}$$

Full text

I Introduction

At the quantum level light-matter interactions are conveniently described in a microscopic Hopfield picture, taking matter degrees of freedom as bosonic fields coupled to the photons. The system Hamiltonian is then diagonalised utilising coupled light-matter quasiparticles, termed polaritons Hopfield58 . This prescription can be followed for quantization of fields in homogeneous, linear dielectrics in the absence of Huttner91 and including Huttner92 losses. Hopfield approaches are particularly advantageous as frames to describe nonlinear processes, which can be understood as scattering cross-sections for the polaritons Carusotto13 ; Ciuti00 ; DeLiberato09 ; DeLiberato13 ; Barachati15 ; Savenko11 ; Kavokin10 ; Kavokin12 ; Liew13 .

Exploiting recent advances in photonics, demonstrated archetypically by the plasmons extant at metal-dielectric interfaces formed by hybridisation of photons with coherent oscillations of the electron gas, energy can be confined in deep subdiffraction volumes. This allows for downscaling of photonic devices to the nanoscale Gramotnev10 ; Kim14 and enables measurable nonlinear scattering at diminished pump threshold. Plasmonic nonlinear devices often exploit intrinsic metal nonlinearities, arising from free carriers or band to band transitions and first demonstrated for harmonic generation mediated by Kretschmann-coupled surface plasmons of a sub-wavelength silver film Tsang96 . Metal-mediated plasmon generation by four-wave mixing Palomba08 and frequency difference generation Constant16 have also been demonstrated, however the centrosymmetry of plasmonic metals precludes second order nonlinear effects in the bulk and limits studies to symmetry-breaking interfaces Butet15 . Instead, theoretical studies in nonlinear plasmonics focus toward third-order nonlinear processes such as self-phase modulation DeLeon14 or other Kerr nonlinearities Huang09 .

Interfaces between polar and non-polar dielectrics also permit deep subdiffraction localisation by hybridising photons with coherent oscillations of the polar crystal atomic lattice, in surface phonon polaritons. Supported in the Reststrahlen region, between transverse and longitudinal optic phonon frequencies which ordinarily occupy the midinfrared spectral region, these modes absent themselves from the Ohmic losses inherent in plasmonic systems and render achievable quality factors an order of magnitude larger. Moreover, interplay between extreme energy localisation in user-defined deep sub-wavelength resonators Caldwell13 ; Gubbin16b and bulk second-order nonlinearity in polar crystals make phonon polaritons an excellent testbed for midinfrared nonlinear optics Paarmann16 ; Razdolski16 .

These inhomogeneous photonic systems require polaritonic descriptions which also account for real-space variation in the dielectric environment, a necessity recently demonstrated in description of the healing length in polaritonic condensates Elistratov16 . Following initial efforts Suttorp04 ; Todorov14 , a full polaritonic description was recently developed for non-magnetic, arbitrarily inhomogeneous, linearly polarisable media Gubbin16c , opening the path for quantum descriptions of nonlinear processes in sub-wavelength devices. In this Paper we embark on development of a quantum, microscopic description of second-order nonlinear scattering in inhomogeneous phonon-polariton systems. In Sec. II, after having sketched the main points of the theory of Hopfield diagonalisation in linear inhomogeneous media previously reported Gubbin16c , we specialise our discussion to the $\beta$-SiC/vacuum interface which has been the subject of contemporary study Paarmann16 , deriving its relevant polaritonic modes. In Sec. III we develop the theory of polaritonic $\chi^{(2)}$ scattering, exploiting data regarding plane-wave scattering from recent second harmonic generation experiments Paarmann16 to both validate our theory and determining the phenomenological couplings that parameterise it. Exploiting those coefficients, that are in good agreement with data from ab initio simulations present in the literature Vanderbilt86 , we will then be able to investigate nonlinear processes involving subdiffraction phonon polariton modes.

II Linear Regime

II.1 General theory

Following a prior treatment Gubbin16c , we describe a non-magnetic, inhomogeneous light-matter system, characterised by a single dipolar resonance with Hamiltonian

[TABLE]

where $\hat{\mathbf{D}}$ $(\hat{\mathbf{H}})$ is the electric displacement (magnetic) field operator and $\hat{\mathbf{X}}$ $(\hat{\mathbf{P}})$ is the material displacement (momentum) field operator. The system is fully parameterised by spatially dependant density $\rho$, longitudinal frequency $\omega_{\mathrm{L}}$, and macroscopic dipole moment $\kappa$. In the spirit of the Hopfield prescription, the Hamiltonian in Eq. 1 can be diagonalised in terms of linear superpositions of light and matter conjugate operators through the polariton operator

[TABLE]

where the Greek symbols are the space-dependent Hopfield coefficients Gubbin16c and $n$ indexes the different positive-frequency normal modes. Such operators obey the Heisenberg equation

[TABLE]

which are equivalent to Maxwell equations over the Hopfield coefficients in an inhomogeneous dielectric. Together with the requirement that the polaritonic operators obey bosonic commutation relations

[TABLE]

this allows us to solve the system and to recover the properly normalised light and matter fields as linear superpositions of polaritonic modes weighted by the wavefunctions $\bm{\alpha}_{n}$

[TABLE]

where the transverse resonance of the matter is

[TABLE]

II.2 Application to surface phonon polaritons

The theory sketched above can be naturally applied to the case of an halfspace filled with a polar dielectric, sketched in Figure 1a, identifying $\omega_{\text{L}}$ and $\omega_{\text{T}}$ with longitudinal and transverse optical phonon frequencies $\omega_{\text{LO}}$ and $\omega_{\text{TO}}$. Developing Eq. 3 this would lead to the solve Maxwell equations with the inhomogeneous dielectric function

[TABLE]

where we labelled the mediums $1$ and $2$ and include off-resonance effects through the phenomenological high frequency dielectric constant $\epsilon_{\infty}$. These effects could be formally included by considering an additional matter resonance with longitudinal frequency far above the Reststrahlen band. In this geometry the eigenproblem in Eq. 3 has been solved Gubbin16c , giving rise to seven classes of solutions. Two (TMv and TEv) describe photons incident upon the surface coming from the vacuum and their respective reflected and transmitted waves, with either Transverse Magnetic or Transverse Electric polarizations. Four (TMl, TEl, TMu, TEu) describe excitations incident upon the surface from the dielectric side and their respective reflected and transmitted waves, indexed by both the polarization and by the polaritonic branch, either lower or upper, they belong to. An example of each group is shown in Figure 1a. The last one (S) is instead a surface bound evanescent solution describing the branch of surface phonon polaritons. The mode index $n$ will thus consist of the solution class and the relevant wavevector: in-plane for the evanescent solutions and both in- and out-of-plane in the medium of origin for the propagative ones. The solution of Eq. 3 leads to the Helmholtz equation, linking for each solution the in-plane wavevector $k_{\parallel}$ and the, generally complex, out-of-plane components $k_{i,z}$ in each halfspace $i$ by

[TABLE]

where $k_{0}=\tfrac{\omega}{c}$ is the vacuum wavelength. As an example in the case of propagative, TM polarised radiation incident from $z>0$ (TMv) with in- and out-of-plane wavevectors $\mathbf{k}_{\parallel}$ and $k_{1,z}$, the solution of Eq. 3 is

[TABLE]

where $N_{\mathbf{k}}^{\text{TMv}}$ is a normalization constant to be determined by the modal orthonormality condition in Eq. 4, $t_{12}^{\text{TM}}\;(r_{12}^{\text{TM}})$ are the Fresnel transmission (reflection) coefficients for a TM polarised beam propagating from medium 1 to medium 2 and the unit vectors for TM polarisation are given by

[TABLE]

where the subscript $i$ indicates the medium of propagation and the index $\pm$ the direction. Here $(\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z})$ form a right-handed orthonormal basis and $\mathbf{e}_{\parallel}$ is the in-plane unit vector parallel to $\mathbf{k}_{\parallel}$. The TE polarised solutions can be found by simply replacing the ${\mathbf{p}}_{i\pm}$ with the plane unit vectors perpendicular to $\mathbf{k}_{\parallel}$ and the Fresnel coefficients with their TE polarised analogues. For real incident wavevectors the normalisation constant can be calculated as

[TABLE]

where $\mathcal{V}$ is the volume of quantization. The evanescent, surface bound modes (S) shown in Figure 2 with in-plane wavevector $\mathbf{k}_{\parallel}$ and dispersion obeying

[TABLE]

have instead the form

[TABLE]

with the normalisation

[TABLE]

where $\mathcal{A}$ is the quantization area and $\nu$ is the ratio of the group and phase velocities in the dielectric halfspace Gubbin16c . We can therefore write the quantized field in the dielectric halfspace in the simple form

[TABLE]

with $\kappa_{2,z}=ik_{2,z}$ real and $\mathcal{V}_{\mathrm{eff}}$ can be interpreted as a dispersive effective mode volume for the evanescent mode Archambault10 .

III Results and discussion

III.1 Nonlinear Quantum Theory

III.1.1 General Theory

The quadratic treatment implicit in Eq. 1 amounts to the lowest order term in the Taylor expansion of the full Hamiltonian. The most general $N^{\text{th}}$ order neglected term in the case of a local interaction will consist of the product of light or matter operators

[TABLE]

where the spatial coordinates are indexed by $j$s and the $t$s index the different light and matter operators, in our case $\hat{\bm{\mathcal{O}}}_{t}\in\left[\hat{\mathbf{D}},\hat{\mathbf{H}},\hat{\mathbf{X}},\hat{\mathbf{P}}\right]$. Using the expression of the field operators in Eq. 5 and their equivalents for the conjugate fields, we can rewrite Eq. 16 in the form of a $N^{\text{th}}$ order polaritonic scattering term

[TABLE]

where the scattering tensor $\Xi$ can be found by carrying out integrals over the relevant spatial variables of the usual nonlinear tensor $\Phi$.

In the following we will aim to describe $\chi^{(2)}$ ($N=3$) effects in inhomogeneous $\beta$-SiC. As the interaction is mediated by electric field coupling to charges position, there are only three nonzero field combinations in Eq. 16: the purely mechanical nonlinearity $\mathbf{\hat{X}\hat{X}\hat{X}}$, the electrical one $\mathbf{\hat{X}\hat{X}\hat{E}}$, and the Raman scattering $\mathbf{\hat{X}\hat{E}\hat{E}}$. Due to the crystal lattice symmetry, in each case the nonlinear tensor is characterised by the single parameter coupling three orthogonal field components Roman06 , leaving us with three independent coupling coefficients to be fixed.

Inserting Eq. 5 inside of Eq. 16 and introducing the second order susceptibility

[TABLE]

with

[TABLE]

the interaction Hamiltonian can be put in the usual nonlinear optics form

[TABLE]

where from Eq. 5

[TABLE]

is the Cartesian component $j$ of the electric field in the polariton mode $n$, and the sum is over all the triplets of polariton modes with different Cartesian coordinates. Introducing the Born effective charge per primitive cell $Z$ Lee03 , primitive cell reduced mass $M$, and primitive cell volume $v$, we can relate the three macroscopic dimensionless coupling parameters in Eq. 18 to the microscopic parameters usually derived in lattice dynamics simulations as

[TABLE]

where $\alpha_{TO}$ is the Raman polarisability per primitive cell, $\mu^{(2)}$ is the second order dipole moment per primitive cell, and $\phi^{(3)}$ the third order lattice potential Roman06 ; Flytzanis72 . Notice that we have ignored the contribution due to the high-frequency static second order susceptibility $\chi_{\infty}^{(2)}$, that would be negligible in the resonant regime studied.

III.1.2 Second Harmonic Generation at a $\beta$-SiC/vacuum interface

We will now apply the theory developed to the case experimentally investigated by Paarmann et al Paarmann16 : second harmonic generation at a $\beta$-SiC/vacuum interface. We have thus to calculate the scattering tensor $\Xi$ by performing the integral in Eq. 20 using the TM polarised plane-wave eigenmodes Eq. II.2 or their TE polarised analogues into the nonlinear Hamiltonian and perform the space integral. As explained in Sec. II B, each modal index $n$ stands for the solution class and the relevant two- or three-dimensional wavevector in the medium of origin. For simplicity we consider the case illustrated in Figure 1a where the harmonic pump photons are azimuthally coplanar (all the $\mathbf{k}_{\parallel}$ are aligned). Remembering that the nonlinear third order tensor in zincblende crystals is characterised by a single value, coupling three orthogonally polarised fields Vanderbilt86 , we can calculate the scattering coefficient of two TEv harmonic pumps (TE polarised incident from $z>0$, indexed by $1$ and $2$) generating a TMu second harmonic (TM polarised upper polariton mode, indexed by $3$)

[TABLE]

where $\mathcal{P}$ indicates the principal value and the choice of a TMu second harmonic is fixed by the requirements of having a single rays escaping toward $z>0$ and a resonance frequency larger than $\omega_{\text{LO}}$. The mismatch in the out-of-plane wavevectors in Eq. 25 is given by

[TABLE]

and the Fresnel tensor has the form

[TABLE]

with

[TABLE]

and the non-indexed out-of-plane wavevectors and the frequencies are calculated through the Helmholtz equation in Eq. 8. Utilising Fermi’s golden rule and summing over final states with different out-of-plane wavevectors we can now calculate the scattering rate in the TM upper polariton mode of frequency $\omega$

[TABLE]

where $N_{j}$ is the photon number in the pump beam $j$ and

[TABLE]

with $\mathrm{v_{g}}(\omega)$ the group velocity. Note that in Eq. 29 the term proportional to $\delta\left(\Delta k_{z}^{n_{1},n_{2},n_{3}}\right)$ disappears because for the second harmonic generation we wish to consider $\Delta k_{z}^{n_{1},n_{2},n_{3}}=0$ cannot be satisfied together with the conservation of energy enforced by $\delta(\omega_{n_{1}}+\omega_{n_{2}}-\omega)$. Specialising Eq. 29 to the case of second harmonic generation with harmonic pump frequency $\omega_{0}$, and defining the energy density for each pump $\mathcal{F}_{j}=\hbar\omega_{0}N_{j}/\mathcal{V}$, we finally obtain the emission rate per unit of surface

[TABLE]

In order to quantitatively fit the experimental results Paarmann16 we now need to introduce the effect of losses into Eq. 31. While we could have used from the beginning the lossy version of the theory Gubbin16c , this would have been needlessly cumbersome given the low losses at the relevant frequencies in SiC, that mainly only result in a broadening of the resonances. We can thus simply generalise Eq. 31 to the case of weak loss by considering a material damping $\gamma$ in the dielectric function in Eq. II.2

[TABLE]

and keeping into account the presence of the complex poles in the dispersive functions $D\left(\omega\right)$ comprising $\chi^{(2)}$ in Eq. 18

[TABLE]

Using the lossy version of Eq. 31 we can fit to the experimental data Paarmann16 , where a free electron laser is used to measure the second harmonics generated in reflection at a $\beta$-SiC surface. Notice that form the the form of the normal modes in Eq. II.2, and as clearly shown in Figure 1a, the emitted TMu modes contain both a part propagating inside the dielectric and one part radiating in vacuum. As the experiment collects and measures only the latter part, in order to correctly fit the experiments we need to multiply the emission rate obtained from Eq. 31 for the Fresnel power transmittance coefficient for the upward travelling outgoing beam $\lvert t_{21}^{\mathrm{TM}}\rvert^{2}$. In Figure 1b we show the theoretical result from Eq. 31 (blue solid line) compared to the experimental data, kindly made available by Paarmann et al Paarmann16 (red dots). Our theory correctly replicates the position of the two main resonances seen in the experimental data without adjustable parameters. The peak near $\omega_{\mathrm{LO}}$ results from a resonance in the pump Fresnel tensor, while the one at $\omega_{\mathrm{TO}}$ is a result of strong resonance of the second order susceptibility, plot for a second harmonic process in Figure 1c, due to the pole in Eq. 33.

In order to quantitatively reproduce also the peak intensities we used the coupling coefficients $G_{1},\;G_{2},\;G_{3}$ introduced in in Eq. 18 as fitting parameters, yielding a good quantitative agreement and thus fixing the ratios of the different coupling constants. Combining these results with measurements of the Faust-Henry coefficient $G_{1}/\chi_{\infty}^{(2)}$ which describes the Raman anharmonicity Yugami86 and of the unclampled-ion, strain free, electro-optic susceptibility Tang91

[TABLE]

it is therefore possible to calculate the absolute magnitude of the coupling parameters $G_{1},\;G_{2},\;G_{3}$ which fully define the second-order nonlinear susceptibility in the neighbourhood of the lattice optical phonon resonances. The third-order lattice potential $\phi^{(3)}$ is thus easily calculable from the fit by Eq. 24 and it could be also independently calculated by ab initio methods. While such ab initio simulations for $\beta$-SiC are not present in the literature, and they are beyond the scope of the present paper, an estimate can be derived by the frozen phonon calculations of Vanderbilt et al. for the third-order lattice potentials of monoatomic crystals in diamond structure Vanderbilt86 . In such a work the authors note that the zone-center nonlinear coefficients rescaled on the bond length $a$ and spring constant $g$

[TABLE]

present a remarkable universality for materials as different as C, Si, and Ge, for which $\tilde{\phi}^{(3)}\approx-4$. Extrapolating this universality to be at least partially valid for $\beta$-SiC, that shares both crystal structure and atomic constituents with them, we obtain a qualitative agreement with our result $\tilde{\phi}^{(3)}=-1.35$.

III.2 Difference Generation of SPhPs

Having fixed material parameters through reproducing experimental results we are now in a position to make predictions about the rates for novel non-linear processes. The excitation of subdiffraction surface modes by free space radiation was first proposed by Novotny et al., utilising a four-wave mixing procedure to excite the plasmons of a gold halfspace Palomba08 . A four-wave mixing procedure is necessary to conserve energy on excitation of modes in the visible spectral region utilising common sources emitting in the UV-visible spectral region. Subdiffraction modes in the midinfrared spectral region however can be generated by three-body scattering of pump beams operating in the visible spectral region. This was recently demonstrated to excite the low-energy plasmons of monolayer graphene from 5 to 50THz Constant16 . Utilising the quantization for the evanescent S modes Eq. 13 we can describe to the process illustrated in Figure2a, describing difference generation of surface phonon polaritons. At the quantum level the process can be described by decay of a photon from the pump into a surface mode accompanied by stimulated emission of a photon into the probe. Considering for definiteness the matrix elements for such a process involving TE polarised pump (1) and probe (2) beams with a surface evanescent signal (3) we find

[TABLE]

where the Fresnel tensor now reads

[TABLE]

and the relevant second order non-linear susceptibility of the $\beta$-SiC halfspace as a function of the signal frequency is shown in Figure 2b. This is substantially less than the nonlinear susceptibility for second harmonic generation illustrated in Figure 1c, because only the signal beam is close to the $\omega_{\text{TO}}$ resonance. The surface mode scattering rate can be now determined by assuming Lorentzian broadening for a surface mode with in-plane wavevector $\mathbf{k}_{\parallel}$

[TABLE]

where the surface mode real and imaginary frequencies $\mathrm{Re}\left[\omega^{\mathrm{S}}_{k_{\parallel}}\right]$ and $\mathrm{Im}\left[\omega^{\mathrm{S}}_{k_{\parallel}}\right]$ are solutions of Eq. 12 with the dissipative dielectric function in Eq. 32. Fixing the pump’s frequency $\omega_{1}=17000$cm*-1* and its incidence angle to $30^{\circ}$ (thus effectively fixing the in-plane wavevector) we can now vary the probe frequency and inclination, leading to resonant emission at frequency $\omega$ in the surface mode described by the emission rate per unit surface

[TABLE]

where $\mathcal{L}_{\mathrm{eff}}=\mathcal{V}_{\mathrm{eff}}/\mathcal{A}$ is the effective confinement length along the $z$ direction.

III.3 Second Harmonic Generation in Subdiffraction Resonators

From Eq. 31 and Eq. 39 we can see that the intensity of second-order nonlinear response can be modified engineering two main parameters. First the $\chi^{(2)}$ of the material facilitating the nonlinear interaction which, for phononic nonlinearities, peaks at the transverse optical phonon frequency. Second the overlap between the participating modes, leading to the effective confinement factor in Eq. 39 when one or more of the modes are evanescent or localised. Both of these factors can be improved by exploiting the morphologically dependant resonances of subdiffraction $\beta$-SiC resonators which have recently been investigated experimentally Caldwell13 ; Gubbin16a and numerically Gubbin16b ; Chen14 . These modes provide spectral tuneablility of the strong field enhancement, allowing a greater overlap of the confined mode with the resonance of the material $\chi^{(2)}$. In addition by confining the mode in three dimensions it is possible to achieve ultrasmall mode volumes up to $10^{6}$ times smaller than the wavelength limited volume in SiC, which provide strong enhancements to the local photonic density of states. Recently second harmonic generation exploiting the resonances of subdiffraction 4H and 6H-SiC resonators was demonstrated experimentally Razdolski16 . In this work the interplay between field confinement and the dispersion of the $\chi^{(2)}$ was directly observed in the form of resonant enhancement of second harmonic generation in the region of a band-folded lattice phonon branch. Our model can be readily extended to describe these systems utilising numerical simulations of the quantized system eigenmodes in the linear regime Gubbin16b . The overlap integrals in the nonlinear Hamiltonian Eq. 16 can then be readily calculated.

IV Conclusion

In conclusion we have presented a comprehensive quantum theory of $\chi^{(2)}$ scattering in inhomogeneous light-matter systems, applying it to model recent nonlinear experiments on $\beta$-SiC surfaces. The quantitative agreement between the resonant frequencies predicted by our theory and the experimental data previously reported Paarmann16 validated our theory, while a fitting of the relative peak heights allowed us to fix the phenomenological scattering parameters. With a fully parameterised theory we were then able to make quantitative prediction of novel nonlinear effects. Phonon polaritons, either on flat surfaces, or localised in sub-wavelength resonators, have the potential to become an important platform for midinfrared photonics and quantum polaritonics. The present paper is a first, important step in gaining a quantitative understanding of their scattering properties at the quantum level, necessary to open the door to future investigations on the possibility to realise quantum fluids of light Carusotto13 made of phonon polaritons, in a sort of midinfrared analogous of the very successful investigations of microcavity exciton-polaritons in the near-infrared region.

V Acknowledgements

The authors thank Alex Paarmann for provision of data used to fit our scattering coefficients. We acknowledge support from EPSRC grant EP/M003183/1. S.D.L. is Royal Society Research Fellow.