Not yet reviewed by Pith; the record is open.

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.