  Phonon polaritons (PhPs) in freestanding and supported multilayers (MuLs) of
hexagonal boron nitride (hBN) are systematically studied using a macroscopic
optical-phonon model. The PhP properties such as confinement, group velocity,
propagation quality factor (PQF) and wavelength scaling are studied. Owing to
the nonlocal high-frequency screening, there is an upper frequency limit making
the two-dimensional (2D) PhPs have a frequency band, and also a maximum PQF
occurs near the centre frequency. The substrate's dielectric response should be
included to accurately calculate the PhP properties. While the simple
electrostatic approximation (ESA) is a proper treatment for PhP frequencies
$\omega$ above $\omega_0$ (e.g. $\omega>1.03\omega_0$ for the 30-layers;
$\omega_0$ is the $\Gamma$ point optical phonon frequency), it fails to
describe the PhP properties near $\omega_0$ and the effect of retardation
should be included for an accurate description. The PhP wavelength versus the
layer thickness near $\omega_0$ deviates significantly from a linear scaling
law given by the ESA due to strong phonon-photon coupling. The calculated PhP
dispersion, PQF and scaling are compared with experimental data of a number of
spectroscopic studies and good agreement is obtained. While the frequency of
incident light should be near the centre frequency to maximize the PQF, the PhP
wavelength, confinement and propagation length can be engineered by varying the
MuL thickness and its dielectric environment.
