Phonon polaritons (PHPs) in freestanding and supported multilayers of hexagonal boron nitride 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 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 center frequency. The substrate’s dielectric response should be included to accurately calculate the PHP properties. While the simple electrostatic approximation (ESA) is valid for PHPs with frequencies [Formula: see text] above [Formula: see text] (e.g., [Formula: see text] for the 30-layers; [Formula: see text] is the [Formula: see text] point optical-phonon frequency), it fails to describe the PHP properties near [Formula: see text] and the effect of retardation should be included for a proper description. The PHP wavelength vs layer thickness near [Formula: see text] deviates significantly from a linear scaling law given by the ESA due to strong coupling of photons and longitudinal optical phonons. The calculated PHP dispersion and scaling are compared with experimental data of a number of spectroscopic studies and found to be in good agreement for most of the results. While the frequency of incident light should be near the center frequency to maximize the PQF, the PHP wavelength, confinement, and propagation length can be engineered by varying the multilayer thickness and its dielectric environment.