\section*{Results and Discussion}
In this Letter we will consider the zincblende polytype $\beta$-SiC, whose dielectric function is negative in the Reststrahlen band between the transverse and longitudinal optic phonon frequencies $\omega_{\mathrm{TO}} = 797.5/$cm and $\omega_{\mathrm{LO}}= 977/$cm. We consider the high-frequency dielectric constant $\epsilon_{\infty} = 6.49$ and transverse phonon damping rate $\gamma = 4$/cm\cite{Paarmann16}. In this region $\beta$-SiC supports surface phonon polaritons, allowing sub-diffraction light localisation by transient storage of energy in the kinetic and potential energy of the ionic lattice. 
Surface phonon polaritons are a promising platform for mid-infrared nonlinear optics due to their ultra-low mode volumes, high quality factors, and the intrinsic anharmonicity of the host polar dielectric. As we lack measurements of the required nonlinear parameters for $\beta$-SiC, in this Letter we utilise the {\it ab initio} results of Vanderbilt {\it et al.}, who studied second and third order anharmonic phonon-phonon scattering in various materials with a diamond crystal structure. They observed a certain universality in coefficients governing the mechanical anharmonicity due to phonon-phonon scattering on rescaling by the bond length and spring constant \cite{Vanderbilt86}. Our recent results demonstrate that such values lead to reasonable estimates, at least for second order nonlinearities, also for $\beta$-SiC  \cite{Gubbin17b}, and we expect the same should hold for third order nonlinearities. Intriguingly, the predicted third-order nonlinear anharmonic scattering coefficients are attractive, contrasting the repulsive material interactions in microcavity polariton systems. As shown in the Supplemental Information the $\chi^{(3)}$ nonlinear susceptibility is composed of a sum of resonant terms, with poles of order 1 to 4 close to the optical photon frequency $\omega_{\mathrm{TO}}$. Given that all modes participating in the processes studied in this Letter lie in the Reststrahlen band, in the neighbourhood of the highest fourth order pole of $\chi^{(3)}$, we can ignore lower order contributions to the optical nonlinearity and consider solely the pure mechanical contribution due to phonon-phonon scattering. Furthermore in $\beta$-SiC only the $xxxx, xyxy, xxyy$ components of the nonlinear tensor are non-zero when accounting for screening \cite{Vanderbilt86}.\\ 
In nonlinear optical experiments the available pump fluence is a key parameter determining which phenomena are within experimental reach. We therefore use parameters representative of the free-electron laser previously utilized to probe $\chi^{(2)}$ nonlinearities in polar dielectrics by the group of A. Paarmann \cite{Passler17}, specifically a pulsed excitation of duration $1$pS, energy $2\mu$J, and an elliptical excitation spot with semi-major and semi-minor axes $500\mathrm{\mu}$m and $150\mathrm{\mu}$m respectively corresponding to a pulse fluence of $0.84$mJ/cm$^2$.

\subsection*{Kerr Self-Interaction}
\begin{figure}
	\includegraphics[width=0.35\textwidth]{Figure1}
	\caption{\label{fig:Fig1} a. Sketch of self-interaction, as a frequency shift dependent on excitation number $\mathcal{N}$ from bare frequency $\omega_{\mathbf{p}}$. b. Reflectance of the bilayer SiC/ vacuum interface illustrated in the inset. c. Reflectance of the vacuum/ $200$nm SiC/ vacuum trilayer shown in the inset. On plots b and c the solid line indicates the poles of the reflectance coefficient, signifying the polariton dispersion. Both colour maps are plotted on a logarithmic scale. d. The frequency shift calculated for the bilayer mode and the film modes (solid lines indicating the lower branch and dashed lines the upper) from Eq.~\ref{eq:SPMre} for pump parameters representative of a free electron laser \cite{Passler17}.}
\end{figure}
The Kerr effect manifests as an intensity dependant frequency shift due to the excitations self-interaction, interpretable as an intensity dependant shift in the effective refractive index \cite{DeLeon14b}. The nonlinear Hamiltonian in Eq.~\ref{eq:HamNL}, specialised onto the case of self-interaction ($\mathbf{p}_j = \mathbf{p}$, $\forall j$) and linearised, leads to renormalised modal frequencies $\omega_{\mathbf{p}}' = \omega_{\mathbf{p}}\left(1 + \Delta_{\mathbf{p}}\right)$ as sketched in Fig.~\ref{fig:Fig1}a. The fractional shift is derived in the Supplementary Information and given by
\begin{equation}
\label{eq:SPMre}
	\Delta_{\mathbf{p}} =  \frac{9}{2} \frac{\hbar \omega_{\mathbf{p}} \mathcal{N}_{\mathbf{p}}}{\epsilon_0 \mathrm{V}_{\mathbf{p}}^2}  \sum_{ijkl}  \int \mathrm{d}\mathbf{r} \;  \chi_{ijkl}^{(3)} \left(\mathbf{r}, \omega_{\mathbf{p}},\omega_{\mathbf{p}},\omega_{\mathbf{p}},\omega_{\mathbf{p}}\right)   \bar{\alpha}_{\mathbf{p}}^{i}\left(\mathbf{r}\right) \bar{\alpha}_{\mathbf{p}}^{j}\left(\mathbf{r}\right) \alpha_{\mathbf{p}}^k\left(\mathbf{r}\right) \alpha_{\mathbf{p}}^l\left(\mathbf{r}\right),
\end{equation}
with $\mathcal{N}_{\mathbf{p}}$ the population in the mode under consideration.
In order to assess whether such effects are large enough to be measured we initially consider two 1D inhomogeneous systems, consisting of either a planar $\beta$-SiC halfspace (inset of Fig.~\ref{fig:Fig1}b) or a thin, suspended, $\beta$-SiC film, such as those usually utilised for fabrication of planar photonic crystals \cite{Calusine14} (inset of  Fig.~\ref{fig:Fig1}c). The halfspace system supports a single evanescent surface phonon polariton while the film supports both symmetric and anti-symmetric superpositions of the surface phonons at each interface. The modes of the system, found considering the poles of the reflection coefficients for the stack, are shown in Figs.~\ref{fig:Fig1}b and \ref{fig:Fig1}c for the bilayer system and a suspended SiC film of thickness $200$nm. \\
Utilising Eq.~\ref{eq:SPMre} we may now calculate the expected fractional frequency shifts, shown in Fig.~\ref{fig:Fig1}d. Larger shifts are observed for thinner films, understandable as a result of increased in-plane confinement for the upper polariton branch and a depression in frequency toward the pole in material $\chi^{(3)}$ at $\omega_{\mathrm{TO}}$ for the lower branch. The decrease in shift for the upper branch at low in-plane wavevector is a result of the mode moving close to the lightline. These modes can be excited via prism coupling \cite{Passler17}. It was assumed that excitation occurred near critical coupling, where approximately $80\%$ of the input beam couples to the surface mode. \\
\begin{figure}
\includegraphics[width=0.35\textwidth]{Figure2}
\caption{\label{fig:Fig2} a. Real resonant frequencies of the lowest order transverse (m=1, red circles) and longitudinal (m=0, blue circles) modes of the single, cylindrical resonator shown in the inset as a function of diameter. The cylinder height is held at $800$nm and the corners are smoothed with radius $2$nm. b. Fractional frequencies shifts for the same cylinder calculated by Eq.~\ref{eq:SPMre}}	
\end{figure}
In order to further explore the nonlinear shift we now consider the modes of sub-wavelength resonators, which support optical resonances with ultra-small photonic mode volumes and correspondingly high energy localisation \cite{Chen14, Gubbin17a}. Here we simply study an isolated cylindrical resonator whose diameter is varied from $800$nm to $50$nm but our results should provide an accurate approximation to those expected for experimentally practical resonator arrays \cite{Gubbin16a}. The modes of such a resonator vary azimuthally as $e^{i m \phi}$ and they are easily probed utilising   COMSOL Multiphysics, allowing for evaluation of the integrated quantities entering Eq.~\ref{eq:SPMre}.\\
The modal frequencies for the lowest order longitudinal $(m=0)$ and transverse $(m=1)$ modes are shown in Fig.~\ref{fig:Fig2}a. The transverse mode is weakly affected by the geometrical change as the cylinder is quasistatic for all sizes studied, while the longitudinal mode shifts to the transverse optic phonon frequency as the diameter decreases. In Fig.~\ref{fig:Fig2}b we show the expected frequency shift. The shift for the transverse mode is far larger than for the longitudinal mode despite the longitudinal mode lying closer to the pole in the $\chi^{(3)}$ at $\omega_{\mathrm{TO}}$ throughout. This is a result of increased optical confinement and modal overlap. The peak in the shift of the dispersionless transverse mode near $600$nm is a result of competition between an increased optical confinement with decreasing diameter and a decrease in the incoupling as a result of falling optical cross section.

\subsection*{Parametric Scattering}
\begin{figure} 
\includegraphics[width=0.8\textwidth]{Figure3}
\caption{\label{fig:Fig3} a. An illustration of parametric scattering between 3 modes with different values of the azimuthal  quantum number $m$. Illustrated are the $z-$component of the electric field for these modes in a cylinder of height $800$nm and of diameter $850$nm. b. Plot of the two sides of the threshold condition in Eq.~\ref{eq:Threshold} as a function of the pulse energy. The stimulated in-scattering is show by a solid line and the loss rate by a dashed line. They cross at the threshold power, indicated by the vertical dot-dashed line. The shading indicates the region experimentally accessible with the laser power reported by the group of A. Paarmann\cite{Passler17}.}
\end{figure}
Having shown that the $\chi^{(3)}$ nonlinearity in SiC is expected to be large enough to lead to experimentally measurable effects, we now investigate the feasibility of a polariton parametric oscillator based on phonon-polaritons. Parametric oscillators exploit optical nonlinearities to convert photons from a pump beam scattering into photons in signal and idler channels. This requires the three modes involved to respect phase matching conditions that can put strong constraints on the feasibility of practical devices. In microcavity polariton systems phase matching can be achieved at a magic angle, relying on the kink in the polariton dispersion \cite{Baumberg00}, or utilising multiple cavities \cite{Diederichs06}. In this regard localised phonon polaritons are a much more flexible system in which to study parametric oscillation due to the morphologically dependant nature of their resonances.\\
We study the general process where two excitations in the central pump mode with frequency $\omega_{C}$, containing a coherent beam with $\mathcal{N}_C\gg1$ excitations,
scatter to one in the upper (idler) and one in the lower (signal) empty modes with frequencies $\omega_{L}$ and $\omega_{U}$. Using the Fermi golden rule the spontaneous emission rate can be put in the form
\begin{align}
	\mathcal{W} =  \frac{3^4 \pi \hbar^2 \omega_{C}^2 \mathcal{N}_C^2}{2^3\epsilon_0^2 \mathrm{V}_{C}^2} \int & \mathrm{d} \omega \;  \omega \left(2 \omega_{C} - \omega \right) \rho_{L}\left(\omega \right)  \rho_{U}\left(2 \omega_{C} - \omega \right) 
	\lvert \tilde{\chi}^{(3)} \left(\omega_{C},\omega_{C},\omega,2 \omega_{C} - \omega\right)   \rvert^2,
	\label{eq:FGRp}
\end{align}
where $\mathcal{N}_C$ is the number of excitations in the pumped mode, $\rho_{\mathbf{p}}\left(\omega\right)$ represents the density of final states of mode $\mathbf{p}$, and $\tilde{\chi}^{(3)}$ is the modal third-order nonlinear coefficient as defined in the Supplementary Information. The product of two output modes densities implies non-negligible emission  only if the parametric scattering condition $\omega_{L} + \omega_{U} = 2 \omega_{C}$ is satisfied within the modes linewidths.\\
Although our theory is applicable to any resonator geometry and can easily be used for practical devices such as periodic arrays of resonators on a same material substrate where the analogy with microcavity polaritonics is most apparent \cite{Gubbin16a}, here we apply it to the illustrative system of an isolated cylindrical resonator. 
The modes of a cylinder vary azimuthally as $e^{i m \phi}$, in this case we consider two transverse $m=1$ modes scattering to one longitudinal $m=0$ mode and one $m=2$ mode as illustrated in Fig.~\ref{fig:Fig3}a. Such a combination has no azimuthal variation, except that arising from the $\chi^{(3)}$ tensor, ensuring a non-vanishing matrix element. The resonant frequencies are calculated utilising a pole-finding algorithm and the COMSOL multiphysics finite element solver and for a diameter of $850$nm we obtain $~911.5/\mathrm{cm}, \; 916.5/\mathrm{cm}, \; 921.5/\mathrm{cm}$ for the $m=1,2,3$ modes respectively with Q-factors $123,\; 133,\; 154$, fulfilling the parametric scattering condition. Following the procedure outlined in the Supplementary Information we can then write
\begin{align}
	\mathcal{W} &= \frac{3^4 \hbar^2 \omega_{C}^2 \mathcal{N}_C^2}{2^3 \epsilon_0^2 \mathrm{V}_{C}^2} \lvert \tilde{\chi}^{(3)} \left(\omega_{C},\omega_{C},\omega_{L}, \omega_{U}\right) \rvert^2  \frac{\gamma_U \gamma_L}{\left(\gamma_U + \gamma_L\right)^2}\left[ \omega_{L} \mathcal{Q}_{U} + \mathcal{Q}_{L} \omega_{U}\right],
\end{align}
where the $\mathcal{Q}_{\mathbf{p}}$ are the quality factors of the respective output modes. We can utilise this scattering rate, in conjunction with the known pump rate  and loss rates calculated from simulations of the linear electromagnetics of the resonator to form a set of rate equations linking the population of the pumped mode $\mathcal{N}_C$ to the ones of the output modes $\mathcal{N}_L$ and $\mathcal{N}_U$. Solving these equations in the undepleted regime under the assumption that the loss rates for each mode are approximately equal and given by $\gamma$, we can derive the threshold condition
\begin{equation}
	2 \mathcal{W} \mathcal{N}_{C,\mathrm{th}}^2  = \gamma,
	\label{eq:Threshold}
\end{equation}
from which we can calculate the required threshold energy per pulse utilising the free-electron lasers characteristic pulse length and spot size extrapolated to continuous wave pumping as $139\mathrm{nJ}$. We plot the two sides of this equation in Fig.~\ref{fig:Fig3}b as a function of the pulse energy, indicating the threshold by the vertical dot-dashed line and shading the region below the single pulse energy of $6 \mathrm{\mu J}$ accessible using a mid-infrared free electron laser  \cite{Passler17}.  
The threshold lies deeply inside the accessible region and this hints tantalisingly that such processes are experimentally attainable in such sub-diffraction resonator systems. Note that here we considered only one particular set of modes fulfilling phase matching conditions, chosen for sake of simplicity. Conservation relations could in fact also be achieved in different geometries, for example by coupling the modes of discrete surface phonon resonators to propagative surface phonon polaritons on a planar interface \cite{Gubbin16a}.\\

\section*{Conclusion}
In conclusion, we have presented a theory of four-wave-mixing in arbitrary phonon-polariton systems. We applied this theory to the measurement of  Kerr self-interaction \cite{DeLeon14b} for both propagative and localised  modes, showing that the relevant shifts are within experimental reach. We further studied parametric scattering between discrete resonator modes, again showing parametric oscillator regime is achievable for reasonable experimental parameters. The theory considered above can be easily applied to realistic experimental geometries where phase matching can be achieved by hybridising the modes of discrete resonators with a propagating phonon polariton \cite{Gubbin16a}. Localised phonon polaritons, thanks to their small mode volumes, long lifetimes, and large nonlinearities, present themselves as an ideal platform for mid-infrared nonlinear polaritonics. We hope the present work will stimulate further investigation in this fascinating domain, which we expect to permit translation of many of the groundbreaking results obtained with microcavity polaritons to a novel and still challenging frequency domain.
 
\section*{Acknowledgement}
The authors thank Joshua D. Caldwell for useful feedback. We acknowledge support from EPSRC grant EP/M003183/1. S.D.L. is Royal Society Research Fellow.
 
