Layer-Resolved Absorption of Light in Arbitrarily Anisotropic Heterostructures

We present a generalized formalism to describe the optical energy flow and spatially resolved absorption in arbitrarily anisotropic layered structures. The algorithm is capable of treating any number of layers of arbitrarily anisotropic, birefringent, and absorbing media and is implemented in an open access computer program. We derive explicit expressions for the transmitted and absorbed power at any point in the multilayer structure, using the electric field distribution from a 4 $\times$ 4 transfer matrix formalism. As a test ground, we study three nanophotonic device structures featuring unique layer-resolved absorption characteristics, making use of (i) in-plane hyperbolic phonon polaritons, (ii) layer-selective, cavity-enhanced exciton absorption in transition metal dichalcogenide monolayers, and (iii) intersubband-cavity polaritons in quantum wells. Covering such a broad spectral range from the far-infrared to the visible, the case studies demonstrate the generality and wide applicability of our approach.


I. INTRODUCTION
The absorption of light in thin layers of strongly anisotropic materials has received enormous attention since the dawn of two-dimensional (2D) materials and their heterostructures 1 . Tremendous progress has been reported using 2D materials for numerous nanophotonic applications such as hybrid graphene-based photodetectors 2,3 , optoelectronic 4,5 and photovoltaic 6,7 devices employing transition metal dichalcogenide (TMDC) monolayers, enhanced lightmatter interaction using photonic integration with optical cavities 8,9 , approaches towards TMDC-based nanolasers 10,11 , hyperlensing 12,13 based on hyperbolic polaritons 14 , bio-sensing 15 , and thermoelectric applications using black phosphorus (BP) monolayers 16,17 . In light of these thriving developments and the great potential entailed in nanophotonic technology, a robust and consistent theoretical framework for the description of light-matter interaction in layered heterostructures of anisotropic materials is of central importance.
In order to understand, analyze and predict the optical response of multilayer structures, the transfer matrix formalism has proven to be of great utility [18][19][20] . In isotropic layered media, a 2 × 2 transfer matrix fully describes any light-matter interaction, and with knowledge of the local electric and magnetic fields, the optical power flow can be described by the Pointing vector S 21,22 . However, what already proves intricate in isotropic multilayers, becomes even more cumbersome when the materials are uniaxial or even biaxial requiring a 4 × 4 transfer matrix formalism [23][24][25] , as it is the case for many state-of-theart nanophotonic materials like hexagonal boron nitride (hBN) 26 , molybdenum trioxide (MoO 3 ) 27 , or BP 28 . In consequence, to the best of our knowledge, previous approaches aiming at the analytical computation of light absorption in anisotropic multilayers are restricted to special cases [29][30][31][32] , whereas a fully generalized formalism applicable to any number of layers of media with arbi-trary permittivity has not been proposed so far.
In this work, we derive explicit expressions for the layer-resolved transmittance and absorption in stratified heterostructures of arbitrarily anisotropic, birefringent, and absorbing media, using the electric field distribution provided by our previous transfer matrix formalism 25,33 . Our algorithm is numerically stable, yields continuous solutions, and is implemented in an open access computer program 34,35 , enabling a robust and consistent framework that is capable of treating light of any polarization impinging at any incident angle onto any number of arbitrarily anisotropic, birefringent, and absorbing layers. To demonstrate the capabilities of our algorithm, we present and discuss simulation results for three nanophotonic device structures, featuring several phenomena such as azimuth-dependent hyperbolic phonon polaritons in a MoO 3 / aluminum nitride (AlN) / silicon carbide (SiC) heterostructure, layer-selective exciton absorption of molybdenum disulfide (MoS 2 ) monolayers in a Fabry-Pérot cavity, and strong light-matter coupling between a cavity mode and an epsilon-near-zero mode in a doped gallium nitride (GaN) multi-quantum well system. Section II summarizes the transfer matrix framework that is used to calculate the electric field distribution and the momenta of the eigenmodes in an anisotropic multilayer system. Based on this theory, section III introduces the calculation of the layer-resolved transmittance and absorption. In section IV, the simulation results are presented.

II. TRANFER MATRIX FRAMEWORK
The 4 × 4 transfer matrix formalism comprising the calculation and sorting of the eigenmodes and the treatment of singularities (Section II A), the calculation of reflection and transmission coefficients (Section II B) and of the electric fields (Section II C) is based on our previous work 25 , and therefore is here only briefly summarized in order to provide the necessary framework for the calculation of the layer-resolved absorption (Section III).

A. Matrix Formalism
The incident medium is taken to be non-absorptive with isotropic (relative) permittivity ε inc , while all other media can feature absorption and fully anisotropic (relative) permittivity tensorsε. Each permittivity tensor ε i of medium i with principle relative permittivities in the crystal frame ε x , ε y , and ε z can be rotated into the lab frame using a three-dimensional coordinate rotation matrix 24 (Eq. 2 in Ref. 25 ). In the following, media with a diagonal permittivity tensor in the lab frame are referred to as non-birefringent, while media with a permittivity tensor that features non-zero off-diagonal elements is called birefringent. Furthermore, all media are assumed to have an isotropic magnetic permeability µ.
The coordinate system in the lab frame is defined such that the multilayer interfaces are parallel to the xy-plane, while the z-direction points from the incident medium towards the substrate and has its origin at the first interface between incident medium and layer i = 1. The layers are indexed from i = 1 to i = N , and the thickness of each layer is d i . Furthermore, i = 0 refers to the incident medium and i = N +1 to the substrate. The plane of incidence is the x-z-plane, yielding the following wavevector k i in layer i: where ω is the circular frequency of the incident light, c is the speed of light in vacuum, ξ = √ ε inc sin(θ) is the inplane x-component of the wavevector which is conserved throughout the entire multilayer system, θ is the incident angle, and q i is the dimensionless z-component of the wavevector in layer i. In any medium, the propagation of an electromagnetic wave is described by exactly four eigenmodes j = 1, 2, 3, 4 with different z-components q ij of the wavevector. These four q ij can be obtained for each medium i individually by solving the eigenvalue problem of a characteristic matrix ∆ (Eq. 11 in Ref. 25 ), as has been derived originally by Berreman 23 . However, for media with highly dispersive permittivities, the four obtained eigenvalues q ij and their related eigenmodes can switch their order as a function of frequency, and thus have to be identified in an unambiguous manner. Following Li et al. 36 , the modes are separated into forward and backward propagating waves according to the sign of q ij (Eq. 12 in Ref. 25 ). We assign the forward propagating (transmitted) waves to be described by q i1 and q i2 , and the backward propagating (reflected) waves by q i3 and q i4 . Furthermore, each pair is sorted by the polarization of the corresponding mode, utilizing the electric fields given by the eigenvectors Ψ ij (Eq. 13 in Ref. 25 ). In non-birefringent media, the two modes are separated into p-polarized (q i1 and q i3 ) and s-polarized (q i2 and q i4 ) waves by analyzing the xcomponent of their electric fields. For birefringent media, on the other hand, the sorting is realized by analyzing the x-component of the Poynting vector S ij = E ij × H ij , and the modes are separated into ordinary (q i1 and q i3 ) and extraordinary (q i2 and q i4 ) waves 25 .
In the case of non-birefringent media, the four solutions q ij become degenerate, leading to singularities in the formalisms of previous works 23,24,37 . To resolve this problem, we follow the solution presented by Xu et al. 38 . Using the appropriately sorted q ij , obtained as described above, the eigenvectors γ ij of the four eigenmodes in each layer i are: with the values of γ ijk given by Xu et al. 38 (Eq. 20 in Ref. 25 ), and k = 1, 2, 3 being the x, y, and z components of γ ij . Furthermore, γ ij has to be normalized: We note that this normalization is essential to ensure a correct calculation of the cross-polarization components of the transfer matrix. The normalized electric field eigenvectors γ ij , being free from singularities, replace the eigenvectors Ψ ij for all further calculations in the formalism. At each interface, the boundary conditions for electric and magnetic fields allow to connect the fields of the two adjacent layers i − 1 and i. Formulated for all four modes simultaneously, the boundary conditions are: where A i is a 4 × 4 matrix calculated from the eigenvectors γ ijk 38 (Eq. 22 in Ref. 25 ), and E i is a dimensionless 4-component electric field vector containing the amplitudes of the resulting electric fields of all four modes. In the following, we refer to E i as the amplitude vector, and its components are sorted as follows: where ⇒ (⇐) stands for the forward propagating, transmitted (backward propagating, reflected) modes, and p, s refers to the p-and s-polarized modes in non-birefringent media, while o, e indicates the ordinary and extraordinary modes in birefringent media. By multiplying A i−1 −1 on both sides of Eq. 4, we find the implicit definition of the interface matrix L i , which projects the amplitude vector in medium i onto the amplitude vector in medium i − 1: For the transition between two birefringent or between two non-birefringent media, the projection of a wave of one particular polarization in layer i only yields a finite amplitude in layer i − 1 of the mode of same polarization, i.e. s/e ↔ s/e, and p/o ↔ p/o. For the transition between a birefringent and a non-birefringent medium, on the other hand, the interface matrix L i projects a mode of one particular polarization in layer i onto both polarization states in layer i − 1. This cross-polarized projection occurs because in birefringent media, the in-plane directions of the ordinary and extraordinary eigenmodes are rotated ( = nπ, n ∈ N 0 ) with respect to the directions of the p-and s-polarized eigenmodes in the non-birefringent medium.
The propagation of all four eigenmodes through layer i is described by the propagation matrix P i 24 : where the rotation of polarization in birefringent media arises due to a phase difference that is accumulated during propagation through the medium, because of different propagation speeds of the ordinary and extraordinary modes (q i1 = q i2 and q i3 = q i4 ). The transfer matrix T i of a single layer i is defined as: and the full transfer matrix Γ of all N layers is where A 0 −1 (A N +1 ) ensures the correct mode projection between the multilayer system and the incident medium (substrate).

B. Reflection and Transmission Coefficients
The full transfer matrix Γ projects the amplitude vector in the substrate E + N +1 onto the amplitude vector in the incident medium E − 0 : where E − i−1 and E + i denote the fields on both sides of the interface between layer i − 1 and i, respectively. Following the equations presented by Yeh 24 , the transmission (t) and reflection (r) coefficients for incident p-or s-polarization can be calculated in terms of the matrix elements of Γ as follows: where the subscripts refer to the incoming and outgoing polarization state, respectively. The transmission coefficients describe the transmitted electric field amplitude into p-and s-polarized states in the case of a nonbirefringent substrate, and into the ordinary and extraordinary eigenstates in the case of a birefringent substrate. We note that the indices of Γ differ from the equations reported by Yeh 24 in order to account for a different order of the eigenmodes in the amplitude vector (Eq. 5).

C. Electric Field Distribution
Employing the interface and propagation matrices L i and P i , respectively, an amplitude vector can be projected to any z-point in the multilayer system (please note the published erratum that corrects the calculation of the electric field distribution in our original work 33 ). However, due to the rotation of polarization in birefringent media and thus the mixing of polarization states, in general, the cases of incident p-and s-polarization have to be treated seperately. As a starting point, the transmission coefficients can be utilized to formulate the amplitude vector E + N +1 for either p-or s-polarized incident light in the substrate at the interface with layer N as follows: where the reflected (⇐) components are set to zero, since no light source is assumed to be on the substrate side of the multilayer system. Furthermore, in order to obtain the electric field amplitudes as a function of z, the propagation through layer i is calculated by means of the propagation matrix P i : with 0 < z < d i being the relative z-position in layer i. Starting from E + N +1 , the interface matrices L i and propagation matrices P i subsequently propagate the amplitude vector towards the incident medium. In the reverse direction, the inverse propagation matrix P N +1 −1 allows to calculate the E-fields in the substrate. As a result, the four mode amplitudes are obtained as a function of z within each layer.
In order to obtain the three-components E x , E y , and E z of the electric field for each of the four modes j, the four mode amplitudes have to be multiplied with their respective eigenmode vector γ ij (Eq. 3). This yields for the electric fields E of the four modes for each layer i, as a function of z, and for either p-or s-polarized incident light: where we have omitted the index i and the z dependence for the sake of readability. The full electric field E i (z) for either p-or s-polarized incident light in layer i at point z is given by the sum of all four electric field vectors: The in-plane components of the sum, E x and E y , are continuous throughout the entire multilayer structure, as it is required by Maxwells boundary conditions.

III. LAYER-RESOLVED TRANSMITTANCE AND ABSORPTION
The total reflectance R of the multilayer system for a given ingoing and outgoing polarization a and b, respectively, can be readily calculated from the corresponding reflection coefficient (Eqs. [11][12][13][14]: The transmittance T , which is the transmitted power into the substrate, on the other hand, in general is not given by the electric field intensity T = |t| 2 (except for the special case if the substrate is vacuum, ε = 1). Instead, T can be calculated from the time-averaged Poynting vector S 39-41 , which describes the direction and magnitude of the energy flux of an electromagnetic wave at any point z in the structure: The full electric field E (for incident polarization a) as a function of z in each layer i was calculated in the previous section (Eq. 18), and the full magnetic field H is obtained as follows using Maxwells equations: where k ij are the wavevectors in layer i of the four modes j, see Eq. 1. Because E and H are known in each layer i and as a function of z from the transfer matrix formalism, the Poynting vector can be evaluated likewise yielding S i (z), which will be used in the following to calculate the transmittance and absorption at any point z in the multilayer system. It is important to note that while E and H can be calulcated for each of the four modes j individually, this mode separation in general -specifically, in the case of birefringent media -is not possible for the Poynting vector S. In other words, in birefringent media, the sum of the Poynting vectors of the four modes is not equal to the Poynting vector calculated from the total fields E (Eq. 18) and H (Eq. 21). The difference arises because in birefringent media, E ⊥ H. Therefore, the cross products E × H between different modes j are no longer zero. For the correct calculation of the Poynting vector in birefringent media, it is thus necessary to calculate the cross product of the total fields E and H, as shown in Eq. 20. Interestingly, this means that the energy flux in birefringent media cannot be split up into the ordinary and extraordinary eigenmodes, but has to be considered as a single quantity. In the following, we therefore discuss the transmittance and absorption for s-or p-polarized incident light without differentiating between the eigenmodes anymore.
An exception is the incident medium, which is set to be isotropic. Here, the Poynting vector can be calculated for each mode individually, and for the purpose of normalizing the transmitted power, we calculate the Poynting vector of the incident light S inc (in layer i = 0 at position z = 0) for either p-or s-polarization as follows: In a stratified multilayer system, the transmitted energy is given by the z-component of the Poynting vector. Thus, we note that alternatively to Eq. 19, the reflectance R can be calculated from the Poynting vector: where the minus sign accounts for the negative zdirection of the reflected light. S b refl,z is the z-component of the Poynting vector of the reflected light of polarization b = p, s (in layer i = 0 at position z = 0) given by: As discussed above, for anisotropic media, a separation of the energy flow into the different eigenmodes of polarization b is not generally possible. Therefore, all following equations calculate the total transmittance or absorption for the respective incident polarization a.
The transmittance T into the substrate i = N + 1 at the interface with layer N for incident light of polarization a is given by: where D = N i=1 d i is the thickness of the multilayer system. The full z-dependence of the transmittance can be evaluated by using the z-component of the Poynting vector S i,z (z) at a certain z-position in layer i: With this, the absorption A of the entire multilayer system, that is up to the last interface between layer N and the substrate, is given by and the z-resolved absorption A i (z) in each layer i is where R a = R ap + R as is the total reflectance. Note that Eq. 28 describes the total absorption starting from z = 0 at the first interface up to the specified position z in layer i. The layer-i-resolved absorption, on the other hand, is given by: where d 1.
is the thickness of all layers through which the incident light has propagated before reaching the layer i.
Before we study three nanophotonic device structures in the following section, we calculate the layer-resolved absorption for two simple test structures in Fig. 1. The first is a typical TMDC heterostructure, comprising monolayers of tungsten diselenide (WSe 2 ), MoS 2 , and tungsten disulfide (WS 2 ) sandwiched between hBN monolayers (sketched in Fig. 1a), where each TMDC monolayer features individual exciton absorption lines. In Fig. 1b, the reflection and transmittance of the entire structure is plotted, and Fig. 1c shows the absorption spectra of each TMDC monolayer. The total absorption (gray line in Fig. 1c) obtained from the reflectance and transmittance spectra (Eq. 27) exhibits three indistinguishable absorption features, whereas the layerresolved absorption calculations unravel the absorption spectrum, allowing to identify the contribution of each TMDC monolayer.
In Fig. 1d-f, we show the absorption in a polar dielectric heterostructure of SiC, AlN, and GaN thin films on a Si substrate probed at infrared (IR) frequencies. Polar dielectric crystals feature an IR-active transverse optical (TO) phonon mode, where light is predominantly absorbed, whereas the longitudinal optical (LO) phonon mode is not IR-active and thus featureless in a bulk crystal. Thin films, on the other hand, support the so called Berreman mode in proximity to the LO frequency, leading to a strong absorption feature at oblique incidence [43][44][45][46] . Thus, the reflectance and transmittance spectra (Fig. 1e) of the polar dielectric heterostructure are of complicated shape, exhibiting six different features arising from the three different polar crystal thin films. The layer-resolved absorption calculations split these features into three spectra with two absorption peaks each (Fig. 1f), allowing to identify the respective polar crystal thin film that leads to the absorption at its respective TO and LO frequencies.

IV. SIMULATIONS OF NANOPHOTONIC DEVICES
The transfer matrix formalism and the calculation of the layer-resolved absorption and transmittance presented in the previous section can be applied for any wavelength and any number of layers, consisting of birefringent or non-birefringent media described by an arbitrary permittivity tensorε i . As case studies, in this section we describe three selected nanophotonic device structures. The first example discusses the hyperbolic The absorption in the MoO 3 , AlN, and SiC layers as a function of ω and Φ and for fixed θ and d i is shown in Fig. 2d-f. The reflectance of the entire system is plotted in Fig. 2c. As required by energy conservation, the sum of the absorbed power in the three polar crystals equals the attenuated power visible as absorption dips in the reflectance. However, while the reflectance only yields the total absorption, the layer-resolved calculations allow to identify the exact position of a power drain in a multilayer system.
In particular, the MoO 3 /AlN/SiC heterostructure features several sharp absorption lines at 660, 800, 920, and 980 cm −1 that are mostly independent of Φ, and one prominent absorption line that strongly varies with Φ, indicating that the latter depends on in-plane anisotropy (ε x = ε y ) while the former do not. In the multilayer sample, only the MoO 3 layer exhibits in-plane anisotropy, while AlN and SiC are c-cut uniaxial crystals (principle relative permittivities ε x , ε y , and ε z are shown in Fig. 2b). Indeed, the Φ-dependent feature is mostly absorbed in the MoO 3 layer. This feature is the hyperbolic phonon polariton (hPhP) supported by the air/MoO 3 interface arising in the in-plane reststrahlen bands of MoO 3 . Its tunability in frequency arises from the large in-plane anisotropy of MoO 3 , which leads to a Φ-dependent effective permittivity sensed by the hPhP upon azimuthal rotation. Notably, below the SiC TO frequencies (∼ 800 cm −1 ), a significant part of the hPhP leaks into the SiC substrate, while above ω SiC TO the mode is mostly confined to the MoO 3 layer. The high confinement above ω SiC TO occurs because of the negative permittivity in the SiC reststrahlen band, while below ω SiC TO , the mode can penentrate the substrate. Interestingly, this mode penetration happens across the AlN layer, where only a small part of the mode is absorbed. Since AlN features its reststrahlen bands across the entire frequency range of the hPhP supported by the MoO 3 , and thus evanescently attenuates all modes, the hPhP appears to tunnel through the AlN layer to be absorbed by the SiC substrate.
The multilayer system presented here only scratches the surface of various possible material configurations and compositions that can potentially be employed for tailoring surface and interface polariton resonances 20,52-57 .
In particular the emerging field of volume-confined hyperbolic polaritons 14,20,58,59 , enabled by the anisotropic permittivity of the supporting media, holds great potential for future nanophotonic applications, such as subdiffraction imaging and hyperlensing 12,13 . Providing the full layer-resolved ab-sorption, our algorithm paves the way to predict and study hyperbolic polariton modes in any anisotropic stratified heterostructure.

B. Layer-selective absorption of MoS2 excitons in a
Fabry-Pérot cavity TMDC monolayers such as MoS 2 feature strong exciton resonances at visible frequencies 42,47 . By inserting these monolayers into van der Waals heterostructures forming a Fabry-Pérot cavity, the light-matter interaction enabled by the TMDC exciton can be strongly enhanced 60 . Here, we embed two MoS 2 monolayers into a d c = 1.9 − 2.4 µm thick hBN cavity with a SiO 2 backreflector on a Si substrate 61 , as sketched in Fig. 3a. Employing the presented formalism, the absorption A i in the two MoS 2 monolayers as a function of photon energy and cavity thickness d c can be calculated.
In the photon energy range of 1.7 − 2.2 eV, MoS 2 features two excitons (A and B) at ∼ 1.9 and ∼ 2.1 eV. These excitons are apparent as resonance peaks in the isotropic relative permittivity plotted in Fig. 3b and re- sulting in peaks in an absorption spectrum. However, the Fabry-Pérot cavity creates a static modulation of the electric field enhancement with peaks and nodes as a function of z position, and thus the resulting absorption in a MoS 2 monolayer not only depends on the photon energy, but also sensitively depends on the z position of the MoS 2 monolayer in the cavity. Taking advantage of this field modulation, we place one MoS 2 monolayer (layer 4) at the center of the cavity where the cavity modes alternate between node and peak with maximal amplitude, and the other MoS 2 monolayer (layer 2) in close proximity where the cavity features a node when there is a peak in the center, and vice versa.
The cavity modes yield the periodic modulation in photon energy and cavity thickness d c that can be seen in the reflectance (transmittance) maps shown in Fig. 3f and g (l and m), for incoming s-and p-polarized light, respectively. The different modulation contrast for R s and R p (T s and T p ) arises from the large incident angle of θ = 70 • , which was chosen to optimize the absorption A s in the MoS 2 monolayers for s-polarized incident light. For smaller incident angles, the differences for sand p-polarization decrease, but with a reduction in the absorption A s .
In Fig. 3h,i (j,k) the absorption for s-and p-polarized incident light in the first MoS 2 monolayer, A s,p 2 (second MoS 2 monolayer, A s,p 4 ), is shown. Due to the choice of the z positions of the two MoS 2 monolayers, each film is sensitive to only every second cavity mode, where layer 2 (first MoS 2 monolayer) absorbs those modes that are not absorbed by layer 4 (second MoS 2 monolayer). Additionally to the absorption modulation imposed by the cavity, the A and B excitons of MoS 2 yield two absorption features at their respective energies, marked by dotted vertical lines in Fig. 3h-k. For the optimized case of s-polarized incident light, the MoS 2 monolayers reach cavity-enhanced absorption values of up to 20% at both exciton energies. At a cavity thickness of d c = 2.15 µm for s-polarization, layer 2 only absorbs at the energy of exciton A, while in layer 4, absorption only occurs at the energy of exciton B. This layer-selective absorption is further illustrated in the absorption spectra (solid lines) for a fixed cavity thickness of d c = 2.15 µm shown in Fig. 3d and e for s-and p-polarized light, respectively. For both polarizations, a high contrast between the two layers at each exciton absorption line is achieved, yielding efficient layer-selectivity.
Finally, we compare these results obtained for an isotropic permittivity model 42 with the same calculations performed for an anisotropic model of MoS 2 47 . In Fig. 3c, the in-plane (ε x,y ) and out-of-plane (ε z ) permittivity values taken from Funke et al. 47 are plotted. While ε x,y are qualitatively the same as the values from Jung et al. 42 (Fig. 3b), ε z differs strongly, taking the almost constant value of ε z = 1 + 0 i. Even though the difference in ε z is substantial, the absorption spectra shown in Fig. 3d,e (dashed lines) are qualitatively identical to the spectra calculated from an isotropic permittivity model (solid lines). This confirms that spectroscopic measurements of TMDC monolayers are mostly insensitive to their out-of-plane permittivity 42 , with the exception of cases where ε z features a zero-crossing, the so-called epsilon-near-zero (ENZ) frequency, giving rise to drastic optical responses such as enhanced higherharmonic generation 45,62,63 . However, this is not the case for MoS 2 , and therefore the spectra are almost identical.
In this example, the layer-resolved absorption calculations from our algorithm provide the essential information for simulating the layer-selective exciton absorption and optimizing the system parameters. In the thriving field of 2D nanophotonics featuring TMDC van der Waals heterostructures, where structures are optimized for maximal light harvesting 6,7 , optoelectronic devices 4,5 or nanolasers 10,11 , such layer-resolved absorption calcu-lations promise to be of essential importance. Due to the generality of the presented algorithm, the light-matter interaction in any 2D heterostructure can be readily investigated, highlighting the broad applicability of our approach.

C. Strong Coupling in a Multi Quantum Well -Cavity System
Doped semiconductor quantum wells (QWs) support transitions between consecutive quantum-confined electronic levels, called intraband or inter-subband (ISB) transitions. Contrary to interband transitions, ISB transitions do not only depend on the bandgap properties of the semiconductor, but also on the electronic confinement inside the QWs, and thus offer a great frequency tunability by changing the width, but also the doping level inside the QW. They play a major role in semiconductor optic devices operating in the IR where semiconducting materials with a suitable bandgap are lacking, and are the building block of quantum well IR photodetectors 64 and quantum cascade lasers 65 . They also offer a practical platform to study the optical properties of dense confined electron gases 66 , which notably lead to the demonstration of the strong 67,68 and ultra-strong light-matter coupling regimes 69,70 . It is remarkable that such fundamental electrodynamical phenomena are directly observable on semiconductor devices [71][72][73] . One peculiar aspect of these ISB transitions is that they only couple to the component of the electric field along the confinement direction of the QW structure. Hence, the optical properties of a doped QW can be described by an effective permittivity tensor with different in-plane (ε x,y ) and out-of-plane (ε z ) values, which has been realized by several permittivity models [74][75][76][77] . The description of light propagation in stratified anisotropic media containing such QWs requires a complex formalism, such as the transfer matrix method the here presented algorithm builds on 25 .
We focus here on an existent experimental configuration to further demonstrate the potential of our formalism for the case of strong light-matter coupling between a cavity mode and a collective intersubband excitation in a multi-quantum well (MQW) superlattice. The system is composed of a GaN cavity formed by a 500 nm thick, Si-doped GaN slab and a Au mirror, as sketched on the top of Fig. 4. The doped GaN layer is modeled using the Drude model, and acts as a low-index mirror. The two mirrors are separated by a 1 µm thick GaN spacer, forming an empty cavity. The system sustains a guided transverse magnetic (TM) mode, where the electric field is confined mostly between the two mirrors and its outof-plane component is maximal near the Au mirror. This guided mode can be probed in a reectance experiment, as discussed in the following. In order to probe large internal angles of incidence experimentally, the sample has to be prepared in a prism shape, for example by cleaving the incident GaN layer facets. The calculated p-polarized reflectance R p is shown in Fig. 4a, evidencing the dispersion relation of the cavity mode with varying incidence angle inside the GaN substrate. The layer-resolved absorption spectra as a function of the angle of incidence for this system are reported in Fig. 4b,c and reveal that absorption occurs mostly in the doped GaN mirror. Notably, for an internal incidence angle of 48 • at the guided mode frequency, all the light is dissipated in the mirrors, leading to a minimum of zero reflectance in Fig. 4a.
We now turn to the system's response when the cavity is partially filled with a MQW structure (Fig. 4dg). The superlattice is composed of twenty repetitions of a 3 nm thick GaN QW, Si-doped with a concentration of 2 × 10 13 cm −2 and twenty-one loss-less 10 nm thick Al 0.26 Ga 0.74 N barriers. Since the guided mode is a TM mode, it naturally provides a component of the electric field along the z-direction for non-zero angles of incidence, which fulfills the ISB transition selection rule. In order to maximize the coupling between the ISB transition in the MQW and the cavity mode, the superlattice is placed where the z-component of the electric field is the largest, that is just below the Au mirror, as shown at the bottom of Fig. 4. The QW dielectric tensor is modeled using a semi-classical approach 75 . We selected the QW dimension and doping level in a way that it sustains a strong, collective electronic excitation known as an intersubband plasmon 66 . The components of the real part of the dielectric permittivity tensor of all the layers are presented in Fig. 4h,i. Note that while the ISB transition in the quantum well has a resonance frequency of ∼ 1600 cm −1 , the Coulomb interaction between the QW electrons results in an ENZ mode at ∼ 2030 cm −1 , marked by the dotted line in Fig. 4i. The ENZ mode dominates the optical response of the QWs.
We show in Fig. 4d the calculated reflectance R p of the cavity containing the MQW. The white dotted line shows the ENZ frequency. A clear anti-crossing can be seen between the cavity mode and the ENZ resonance, which is characteristic of the strong light-matter coupling regime, resulting in two polariton branches. The minimal separation between the two branches amounts to a vacuum Rabi splitting 2Ω R = 200 cm −1 . The layer-resolved absorption spectra as a function of the angle of incidence are shown in Fig. 4e-g for the GaN mirror, the MQW superlattice, and the Au mirror, respectively. When the cavity mode is far detuned from the ENZ frequency, the absorption occurs mostly in the two mirrors, and especially in the GaN mirror, as for the empty cavity case. The situation changes dramatically when the cavity mode is tuned near the ENZ frequency. The absorption then mostly occurs in the MQW. It is however important to note that the absorption is maximal at the frequencies of the two polaritons and not at the ENZ frequency, as it would occur in the case of a weak coupling between the cavity mode and the QW resonance. For an internal incidence angle of 48 • , the maximum absorption inside the MQW stack is now 0.5 at each of the frequencies of the two polaritons.
The algorithm presented here thus allows to directly calculate the light absorption in the complex case of a doped MQW superlattice strongly coupled to a cavity mode. In addition to the known features, such as the avoided crossing in the angle-dependent reflectance spectrum, we can directly calculate the absorption inside the MQW active region, which can be usefully linked to the detected photocurrent in the perspective of using such structures in photodetector devices 67,71,73 .

V. DISCUSSION
We have presented three nanophotonic devices based on anisotropic multilayer structures made from metals, polar dielectrics, and TMDC monolayers, covering incident wavelengths from the far-IR up to the visible. Our formalism allows to calculate the transmittance and absorption in any layer, giving unprecedented insight into the physics of light propagation in anisotropic, and even birefringent, stratified systems.
In recent years, the field of nanophotonics has intensely investigated the optical response of 2D heterostructures. A particularly thriving subject has been polaritonic excitations, which can be supported by a broad variety of systems including slabs of metals, doped semiconductors, polar dielectrics, 2D materials such as TMDC monolayers and their stratified heterostructures 17,54,78 . Key features of polaritons for nanophotonic technologies are their high spatial confinement and field enhancement, which are driven by the particular design of the multilayer materials, stacking order, and layer thicknesses. Potential applications of such systems range from sensing 79,80 and solar cells 81 , over optoelectronic devices 3 and beam manipulation via metamaterials 82 , to waveguiding 56,83 , and ultrafast optical components 45,84 . However, due to the lack of a general formalism, the optical response of these polaritonic multilayer systems is often either approximated by effective, isotropic permittivity models [85][86][87] , or described by specifically derived formulas [29][30][31][32] . Our generalized formalism allows for a precise, layer-resolved study that includes any isotropic, anisotropic or even birefringent response of any number of layers, and thus holds great potential for the prediction and analysis of polariton modes in stratified heterostructures. This is especially relevant for systems where one or more materials feature anisotropic permittivity with spectral regions of hyperbolicity where the principle real permittivities have opposite signs, such as hBN or MoO 3 . Recently, these hyperbolic materials have attracted increasing interest 14,20,58,59 due to the existence of hyperbolic polaritons, featuring novel properties for nanophotonic applications such as subdiffraction imaging and hyperlensing 12,13 . Because of the strong anisotropy of these systems, an effective permittivity approach is not purposeful. Here, our formalism provides the essential theoretical framework that is necessary to model, predict and analyze the optical response of such hyperbolic heterostructures, as we have discussed exemplarily for a MoO 3 /AlN/SiC system (Fig. 2). Providing the full layer-resolved information about the field distribution and power flow of the excited polaritons, our method allows to readily and concisely model and study hyperbolic polariton modes in any anisotropic stratified heterostructure.
Furthermore, the layer-resolved absorption formalism provides a description for designing optoelectronic devices such as detectors, for which the photoresponse is linked to the light absorption solely in the active region of the device. In the case of the MQW system, optimizing the overall light absorption by minimizing both the reflectance and transmittance of the system would not be sufficient to model the performances of a photodetector in the strong light-matter coupling regime 71,73 . Unfortunately, these are the only quantities that can be probed in a reflectance experiment. Calculating the layer-resolved absorption in the structure allows to optimize the cavity and MQW design, aiming at minimizing the light absorption in the cavity mirrors while maximizing the absorption in the MQW. This is well exemplified in Fig  4 e,f where we can see that the light is preferably dissipated in either the doped GaN mirror or the MQW superlattice depending on the detuning between the cavity mode and the QW resonance. This behavior cannot be deduced only from the reflectance measurement simulated in Fig. 4d. Fitting experimental reflectance data using our formalism 25 would allow to retrieve the amount of light dissipated inside the active region only from the experimentally observable quantities, and to estimate figures of merit such as the quantum efficiency of the device. The present method thus provides a convenient way to model and optimize complex, optically anisotropic heterostructures for optoelectronic devices.

VI. CONCLUSION
In this work, we have derived explicit expressions for the calculation of the layer-resolved transmittance and absorption of light propagating in arbitrarily anisotropic, birefringent, and absorbing multilayer media. The algorithm relies on the electric field distribution computed from a 4 × 4 transfer matrix formalism 25,33 , yielding a robust and consistent framework for light-matter interaction in stratified systems of arbitrary permittivity, which is implemented in an open access computer program 34,35 . As case studies, we applied the algorithm to simulations of three nanophotonic device structures featuring hyperbolic phonon polaritons in a polar dielectric heterostructure, MoS 2 excitons in a Fabry-Pérot cavity, and ENZ resonances in a cavity-coupled multiquantum-well, where we observed azimuth-dependent hyperbolic polariton tunneling, layer-selective exciton absorption, and strong coupling between ENZ and cavity modes. Allowing for a detailed analysis of the layerresolved electric field distribution, transmittance, and absorption of light in any multilayer system, our algorithm holds great potential for the prediction of nanophotonic light-matter interactions in arbitrarily anisotropic stratified heterostructures.

VII. ACKNOWLEDGMENTS
We thank M. Wolf, S. Wasserroth and R. Ernstorfer (FHI Berlin) for careful reading of the manuscript and M. Wolf and the Max Planck Society for supporting this work.