  Phonon polariton modes in layered anisotropic heterostructures are a key
building block for modern nanophotonic technologies. The light-matter
interaction for evanescent excitation of such a multilayer system can be
theoretically described by a transfer matrix formalism. This method allows to
compute the imaginary part of the p-polarized reflection coefficient
Im$(r_{pp})$, which is typically used to analyze the polariton dispersion of
the multilayer structure, but lacks the possibility to access the
layer-resolved polaritonic response. We present an approach to compute the
layer-resolved polariton resonance intensity in aribtrarily anisotropic layered
heterostructures, based on calculating the Poynting vector extracted from a
transfer matrix formalism. Our approach is independent of the experimental
excitation conditions, and fulfills an empirical conservation law. As a test
ground, we study two state-of-the-art nanophotonic multilayer systems, covering
strong coupling and tunable hyperbolic surface phonon polaritons in twisted
\MoO~double layers. Providing a new level of insight into the polaritonic
response, our method holds great potential for understanding, optimizing and
predicting new forms of polariton heterostructures in the future.
