Theory of phonon sidebands in the absorption spectra of moiré exciton–polaritons

Excitons in twisted bilayers of transition metal dichalcogenides have strongly modified dispersion relations due to the formation of periodic moiré potentials. The strong coupling to a light field in an optical cavity leads to the appearance of moiré polaritons. In this paper, we derive a theoretical model for the linear absorption spectrum of the coupled moiré polariton–phonon system based on the time-convolutionless (TCL) approach. Results obtained by numerically solving the TCL equation are compared to those obtained in the Markovian limit and from a perturbative treatment of non-Markovian corrections. A key quantity for the interpretation of the findings is the generalized phonon spectral density. We discuss the phonon impact on the spectrum for realistic moiré exciton dispersions by varying twist angle and temperature. Key features introduced by the coupling to phonons are broadenings and energy shifts of the upper and lower polariton peak and the appearance of phonon sidebands between them. We analyze these features with respect to the role of Markovian and non-Markovian effects and find that they strongly depend on the twist angle. We can distinguish between the regimes of large, small, and intermediate twist angles. In the latter phonon effects are particularly pronounced due to dominating phonon transitions into regions which are characterized by van Hove singularities in the density of states.

For larger twist angles and, thus, shorter distances to the neighboring minima the interaction between excitons in different minima increases leading to delocalized exciton states and a more 2D-like behavior [6,23,24].In addition to the modification of the electronic states, the moiré potential also impacts the crystal lattice, leading to the formation of moiré phonon modes [12,14] and a modified exciton-phonon coupling [13,27,28].
When integrating the twisted heterostructure in a photonic cavity formed, for example, by two distributed Bragg reflectors (DBRs), as sketched in Fig. 1 (b), the moiré excitons interact with photons of the cavity leading to the formation of moiré excitonpolaritons [29][30][31].Due to their low effective mass they are interesting candidates for Bose-Einstein condensation [32,33].
In this contribution we investigate the phonon impact on the optical properties of moiré exciton-polaritons for varying twist angles and temperatures.After the introduction of the model (Sec.2.1) and the optical absorption spectrum (Sec.2.2), in Sec.2.3 a correlation expansion up to second order Born approximation [34,35] together with a time-convolutionless (TCL) approach [35,36] are used to derive a closed set of coupled differential equations for the dynamics of the polariton polarization with timedependent coefficients.This approach has been successfully applied to excitons and polaritons in single TMDC layers [35,36], and to the calculation of nonlinear optical signals [37].Here, it will be extended to polaritons in twisted bilayers.In Sec.2.4 the parameters chosen for our calculations are specified.
In Sec.3.1 the TCL equations are solved numerically for a realistic moiré exciton dispersion.This allows us to discuss the optical spectra for varying twist angles and temperatures.In Sec.3.2 the numerically calculated results are compared with a perturbative solution which helps us in obtaining a better understanding of the features in the spectra and in particular of the significance of (non-)Markovian processes.In Sec.3.3 we then analyze the twist angle dependence of the key features of the spectra, in particular the convergence of the perturbative approach as well as the peak positions and widths.Interestingly, the phonon coupling leads to a clear separation of regimes for small and large twist angles, separated by a transition region where the phonon influence is particularly pronounced.Finally, in Sec. 4 we summarize our results with some concluding remarks.

Hamiltonian
We consider excitons in a twisted TMDC bilayer that is embedded in a photonic cavity as sketched in Fig. 1 (b).Due to the twist between the monolayers by an angle ϑ a two-dimensional periodic structure forms, characterized by two moiré lattice vectors R 1 and R 2 , as schematically plotted in Fig. 1 (a).The periodic structure leads to a periodic moiré potential for the excitons, which can be written as where G j is the jth reciprocal moiré lattice vector and ṼG j is the corresponding Fourier component of the moiré potential.In the following we will assume equal lattice constants of the monolayers.Then, for all twist angles except for multiples of 60 • the periodicity of the moiré potential is larger than the periodicity of the monolayers, leading to a 1st moiré Brillouin zone (MBZ) B m which is smaller than the 1st Brillouin zone of the monolayers.This periodic moiré potential leads to the formation of moiré excitons [4-6, 8, 23], described by the creation (annihilation) operators X † n,k (X n,k ) with band index n, center of mass wave vector k with k ∈ B m , and dispersion relation ℏω n,k .The dispersion strongly depends on the twist angle ϑ between the two layers; it is flat for small ϑ [23,24] and becomes more 2D-like with increasing ϑ due to the decreasing period length of the moiré potential in real space.Examples of exciton dispersions for different values of ϑ can be found in Fig. A1 in the Appendix A.
The derivation of the moiré exciton operators and their dispersion relation starting from the electron-hole picture for a homogeneous bilayer without moiré potential is given in Appendix A. We basically follow the derivation in Refs.[6,23].In a first step, operators for excitons of the homogeneous bilayer are introduced by solving a Wannier equation incorporating the electron-hole interaction but without a moiré potential.In a second step, the moiré excitons are introduced by diagonalizing the Hamiltonian of the homogeneous excitons extended by the moiré potential.This involves a back-folding of the homogeneous exciton bands into the 1st MBZ.Note, that in this paper we will use the expressions homogeneous exciton and moiré exciton to distinguish between the excitons of the 2D bilayer without and with the impact of the moiré potential, respectively.
The moiré excitons interact with the phonon modes of the heterostructure.The phonons are described by the creation (annihilation) operators b † j,Q (b j,Q ) for a phonon in the branch j with wave vector Q and dispersion relation ℏΩ j,Q .We neglect the influence of the moiré lattice on the phonons, therefore the phonon wave vector Q is not restricted to the 1st MBZ.The phonons couple to the moiré excitons via the coupling constant ℏG (n,n ′ ) j,k,Q .In Appendix B this coupling constant is derived starting from the electron/hole-phonon coupling Hamiltonian and inserting the definition of the moiré exciton operators from Appendix A.
Embedding the moiré heterostructure in a photonic cavity finally leads to a coupling of the moiré excitons to cavity photons.The photons are described by the creation (annihilation) operators a † K (a K ) for a photon in the considered cavity mode with inplane wave vector K.We restrict the photons to the subspace of a single mode in the out-of-plane direction with discrete wave vector K z and energy ℏω 0 = ℏcK z , c being the velocity of light in the medium, and a continuum of in-plane modes with wave vectors K, leading to the dispersion relation The in-plane photon wave vectors are not restricted to the 1st MBZ, but for each K one can find a reciprocal lattice vector G m , such that In general, moiré excitons with wave vector k couple to all photon modes K that differ only by a reciprocal lattice vector G m .However, assuming, that the 1st MBZ is sufficiently large, the photons from a neighboring MBZ with K = k + G m and G m ̸ = 0 are effectively decoupled from the excitons in the 1st MBZ due to the steep photon dispersion.Therefore, the coupling between exciton and photon modes that differ by a reciprocal lattice vector G m ̸ = 0 can usually be neglected and it is sufficient to take into account photon modes with k ∈ B m .The exciton-photon coupling is then described by the coupling strength λ n,k .In systems where this is not fulfilled or where the cavity field is resonant to excitons close to the edge of the 1st MBZ, a coupling of multiple photon modes to the exciton via photon-umklapp processes may occur, as discussed in Ref. [33].Here, however, we will not consider this case.
The total Hamiltonian for the exciton-photon-phonon system then reads where The exciton-photon part, i.e., the first three sums in Eq. 3, can be diagonalized by introducing the polariton operators with the normalization condition Λ β * k;Λ,n β k,Λ,n ′ = δ n,n ′ , the number of moiré exciton modes N , and the quantum number for the polariton branch Λ with Λ = 0, 1, . . ., N .Exciton and photon operators are obtained from the polariton operators by the inverse transformation For small excitation densities the exciton operators can be treated as bosonic operators with X n,k , X † n ′ ,k ′ ≈ δ n,n ′ δ k,k ′ .This translates to a bosonic commutation relation for the polariton operators The Hamiltonian in this basis then reads with the polariton dispersion relation ℏϵ Λ,k , that can be calculated from the exciton and photon dispersion relations.The polariton-phonon coupling reads with G chosen such that k + Q + G ∈ B m .As seen in Eq. 7, the phonons lead to transitions within one (Λ = Λ ′ ) or between different polariton branches Λ and Λ ′ .This could also involve transitions to momentum-dark polariton branches [38] which, however, will not be considered here.

Absorption spectrum
We are interested in the absorption spectrum of light impinging on the cavity under an angle θ, i.e., light with an in-plane wave vector k with k = (ω/c 0 ) sin(θ), c 0 being the vacuum speed of light.The light couples through the top mirror to the cavity photons with the same in-plane wave vector.The absorption spectrum α k (ω) is then obtained from the response of the cavity mode to a delta-pulse like excitation at time t = 0, which effectively sets an initial value of the cavity field amplitude ⟨a k ⟩ (t = 0).We will measure the photon energy ω relative to the cavity mode frequency ω 0 at k = 0.The spectrum is then given by the real part of the Fourier transform of the field amplitude, where we have used the inverse transformation from the polariton to the photon operator according to Eq. 5.In this paper we will concentrate on the spectra α 0 (ω) for normal incidence, i.e., for k = 0.

Equations of motion
As we have seen in Eq. 9, the quantity of interest for the calculation of the absorption spectrum is the expectation value of the polariton operator P Λ,k .Therefore, we now derive the equation of motion for the polariton polarization PΛ,k = ⟨P Λ,k ⟩ e iω 0 t , conveniently defined in a frame rotating with the frequency of the light field at k = 0. Applying Heisenberg's equation of motion to the polariton operator leads to with the moiré polariton dispersion relative to the photon mode ∆ Λ,k = ϵ Λ,k −ω 0 and the phonon-assisted polariton polarizations (PAPPs) in the same rotating frame, defined as where G is chosen such that k − Q − G ∈ B m and where we have used the relation g j,k+Q+G,−Q .When calculating the equation of motion for S(Λ,Λ ′ ,±) j,k,Q we find higher-order phononassisted expectation values, which finally results in an infinite hierarchy of equations of motion.This hierarchy can be truncated using a correlation expansion [34,35], assuming that correlations among an increasing number of operators become less important.Due to the assumption of a small number of photons and excitons and since we are interested in the linear response spectrum, we neglect all terms with more than one polariton operator.We also restrict our calculations to the 2nd order Born approximation [35], i.e., we neglect all terms that are higher than 2nd order in the polariton-phonon coupling.In combination with the small polariton density this implies that all expectation values of b j,Q , b † j,−Q and all their combinations are constant (except for self oscillations) [35].Assuming that the phonons are initially in a thermal state with temperature T , the expectation value b with . Equations 10 and 12 now form a closed set of equations, that have to be solved for every set of k, Q, j, Λ and Λ ′ .
By formally integrating Eq. 12 we obtain with the initial condition S(Λ,Λ ′ ,±) j,k,Q (0) = 0, in agreement with the assumption of an ultrashort excitation of the light mode at t = 0.The PAPP at time t thus depends on the polarization PΛ ′′ ,k (t − τ ) at all earlier times t − τ with 0 ≤ τ ≤ t, therefore including memory effects in the dynamics.Within the 2nd order Born approximation this expectation value can be approximated by PΛ ′′ ,k (t−τ ) ≈ PΛ ′′ ,k (t) exp [i∆ Λ ′′ ,k τ ], since all other terms lead to contributions that are higher than 2nd order in the polaritonphonon coupling.In the theory of open quantum systems this is called the timeconvolutionless (TCL) approach [39] and it has been found to provide a good description of exciton and polariton spectra in TMDCs [35,36].
Inserting the previous result into Eq. 10 leads to the closed equation for the with time-dependent coefficients where we have introduced the generalized phonon spectral density (gPSD) ρ Λ,Λ ′ (k, Ω), defined as describing phonon-induced transitions within or between the polariton branches at given moiré exciton wave vectors k.Within this approximation the complete phonon influence is dynamically included in the time-dependent coefficients Γ Λ,Λ ′ (k, t) [35,36].
Although memory effects in the polarization have been eliminated in the TCL approach, the coefficients Γ Λ,Λ ′ (k, t) still carry information on the time of the excitation process in their time dependence.This information can be eliminated by additionally performing a Markov approximation, which consists in replacing the coefficients by their long-time limit Γ Λ,Λ ′ (k, t → ∞).Formally, this is achieved by replacing the time integral in Eq. 14 with the help of the Dirac identity according to In this limit we obtain The real part describes the dephasing due to real (i.e., energy conserving) transitions (Ω = 0), while the imaginary part gives rise to energy shifts associated with virtual (i.e., energy non-conserving) transitions.We can therefore interpret the variable Ω in the gPSD [Eq.15] as the energy mismatch of the respective phonon-induced transition.We will come back to this point later when interpreting the phonon impact on our calculated spectra.As discussed in Appendix C, the Markov limit according to Eqs. 16a and 16b can be taken as the starting point for a perturbative solution of Eq. 13.In Sec.3.2 we will compare the numerical results with those obtained from such a perturbative solution, which will provide us with a deeper understanding of the various contributions of the spectrum.

Assumptions and parameters
In moiré crystals the twist angle can be seen as a control parameter that has an impact on the moiré exciton band structure, and therefore the dispersion of the polariton branches [6,15,24,29], the exciton wave functions [24,29], the phonon dispersion [12,40] and the exciton-phonon coupling [13], i.e., all parameters that contribute to the linear absorption spectrum of Eq. 9.
While until now all derivations were general, for the following calculations we consider parameters motivated by moiré excitons in an AA (R-type [41]) stacked MoSe 2 /WSe 2 hetero-bilayer.For the moiré exciton dispersion we take into account a single homogeneous exciton band with an effective mass of M = 0.84 m 0 [24], m 0 being the free electron mass, describing a quadratic dispersion of the homogeneous excitons.Note that a constant energy offset will later be chosen such that the moiré exciton at k = 0 is in resonance with the cavity mode at K = k = 0.A situation where light-matter coupling is resonant at another wave vector K ̸ = 0, as for example realized in Ref. [42] for a TMDC monolayer coupled to an optical cavity, can also lead to interesting effects due to varying proportions of the exciton and light contributions to the polaritons, but will not be discussed in the present work.
The homogeneous excitons experience the periodic moiré potential of Eq. 1.We approximate the potential by taking into account the reciprocal lattice vectors G j to the first neighbors of the 1st MBZ with j = 1, . . ., 6.The numbering is such that G j and G j+1 enclose an angle of 60 • .The Fourier components of the potential are chosen as ṼG j = Ṽ exp(iψ) for j = 1, 3, 5, ṼG j = Ṽ exp(−iψ) for j = 2, 4, 6, and ṼG j = 0 otherwise, with Ṽ = 11.8 meV and phase angle ψ = 79.5 • , [18,24].By solving the eigenvalue equation (A.7) for every k in the 1st MBZ numerically, we obtain the moiré exciton dispersion relation ℏω 1,k .
Figure 2 shows the curvature of the moiré exciton dispersion ℏω 1,k at the γ point in the direction towards the κ (blue solid line) and the m (red dashed line) point and the energies at the κ (green solid line) and m point (violet dashed line) relative to the γ point energy as a function of the twist angle in the range 1 • ≤ ϑ ≤ 9 • .The top axis denotes the absolute value of the real space moiré lattice constant, i.e., the real space moiré periodicity for a given twist angle.The full exciton dispersion relations for three examplary values of the twist angle ϑ can be found in Fig. A1 in Appendix A.
Our first observation is that the curvature is always the same for both depicted high-symmetry directions, showing that the dispersion is isotropic around the γ point of the MBZ.At small twist angles the curvatures at the γ point and the energies at the edge of the 1st MBZ vanish, which reflects the flat moiré exciton dispersion in this limit.When increasing the twist angle, both the curvatures and the energies at the MBZ boundaries increase.At large twist angles the curvatures saturate, i.e., in this region the moiré exciton dispersion in a certain region around the γ point can be approximated in each direction by a parabola with twist-angle independent curvature reflecting the dispersion relation of the homogeneous excitons.The energies at the κ and m points increase with growing angle due to the increasing size of the MBZ.According to the general properties of band structures, the slope of the exciton dispersion at the MBZ boundary vanishes in the direction perpendicular to the boundary, which leads to an energy maximum at the κ point and a saddle point associated with a van Hove singularity in the density of states at the m point.We set the cavity mode energy ℏω 0 to be in resonance with the moiré exciton band at k = 0.The moiré exciton-photon coupling constant is set to λ 1,k = 10 meV, such that the polariton splitting at k = 0 is on the same order of magnitude as in Refs.[29,30].
A sketch of the exciton (dashed red), photon (dashed blue) and polariton (solid red and blue) dispersion in the 1st and 2nd MBZ is shown in Fig. 3.We refer to the region in the MBZ where the exciton-light coupling strongly modifies the dispersion as the polariton region (red area) and to the regions where the dispersion is only weakly affected as exciton (gray area) or photon region (yellow area), respectively.For clarity these areas are only marked in the 2nd MBZ at G = G 1 (right).
For the phonons we assume a single acoustic branch with index j = 0 and linear dispersion Ω 0,Q = c s |Q| that is unaffected by the moiré potential and isotropic in all in-plane directions with sound velocity c s , i.e., the Debye model.The phonon wave vector Q is thus not restricted to the 1st MBZ, but extends over the full Brillouin zone of a single layer.
As derived in Appendix B, the moiré exciton-phonon coupling constant can be written as where g is the coupling constant of electrons/holes, µ e/h = m e/h /(m e + m h ) is the ratio of the electron (hole) mass m e (m h ) to the total exciton mass, F (Q) is the form factor of the homogeneous exciton, resulting from the electron-hole attraction (see Appendix A), and the form factor f (n,n ′ ) k,Q originates from the moiré potential confinement.For the electron/hole-phonon coupling constants in Eq. 17 we assume the common deformation potential coupling [35] g with the mass density of a single layer ρ, normalization volume V , and deformation potential for electrons/holes in the long-wavelength limit D e/h [43], corresponding to the first order perturbation theory of the scattering potential in the phonon wave number.The phonon parameters of the MoSe 2 /WSe 2 bilayer are approximated by the ones of MoSe 2 with density ρ = 4.26 × 10 −7 g/cm 2 , speed of sound c s = 4.1 nm/ps and average deformation potentials D (1) e = 2.4 eV and D (1) h = −1.98 eV [35].The exciton form factor F (Q) is obtained from the 1s exciton wave function of the homogeneous exciton.It satisfies the condition F (0) = 1 and it is localized roughly on the scale of the inverse of the exciton Bohr radius.Here, we approximate it by a Gaussian choosing with the localization length a = 2 nm.The moiré form factor f in Eq. 17 gives rise to an additional twist angle dependence of the moiré exciton-phonon coupling.It also satisfies the condition f k,0 = 1.We assume that the moiré form factor will not strongly affect the general shape of the moiré exciton phonon coupling constant.Therefore we do not expect new features in the resulting spectra and set f k,Q = 1 to focus on the twist angle dependence of the moiré exciton dispersion.
The moiré exciton-phonon coupling leads to phonon induced transitions between polariton states either within the same branch or connecting different branches, including both normal and phonon umklapp processes, as sketched schematically in Fig. 3.
Initially, no excitons and photons are present in the system, while the phonons are in a thermal state with temperature T .We assume that a short laser pulse, applied perpendicular to the cavity surface brings the photons into a coherent state at k = 0 with a small coherent amplitude justifying the limit of small polariton densities.Since Eq. 13 does not couple different k, the assumption of a perpendicular excitation guarantees that only polariton polarizations with k = 0 are excited in the system at all times.We additionally introduce a constant damping rate for each polariton polarization by replacing ℏ∆ Λ,k → ℏ∆ Λ,k − iγ 0 with γ 0 = 0.4 meV to take other dephasing mechanism like photon loss of the cavity or direct photon emission from the exciton to the continuum of unconfined photon modes into account, which results in non-vanishing spectral line widths.

Linear absorption spectra
We start the results part by discussing linear absorption spectra α 0 (ω) of moiré exciton polaritons obtained from a full numerical solution of the TCL equation ( 13), and analyzing the impact of phonons seen in these spectra as a function of various parameters.As seen in Eqs. 13, 14 and 15, the spectrum depends on the gPSD ρ Λ,Λ ′ (k, Ω), which includes the moiré polariton-phonon coupling constants g (Λ,Λ ′ ) j,k,Q , the phonon dispersion relations Ω j,Q and the relative polariton dispersion ∆ Λ,k = ϵ Λ,k − ω 0 with respect to the photon mode frequency ω 0 .Considering all assumptions from the previous section, we notice that ∆ Λ,k is effectively the only quantity that explicitly depends on the twist angle ϑ, as seen in Fig. 2.
In Fig. 4, the top row (a, d, g) shows the linear absorption spectra for the temperatures T = 4, 70, 200 K (from bright to dark colors) and twist angles ϑ = 5 • , 3.5 • , and 1 • (from left to right).The middle row (b, e, h) depicts the corresponding gPSDs ρ Λ,Λ (k, Ω) at k = 0 for the lower (Λ = 0, dashed red) and the upper (Λ = 1, solid blue) polariton as a function of the energy mismatch ℏΩ.Note that for our set of parameters, in particular due to the resonant exciton-photon coupling, ρ 0,Λ (k, Ω) and ρ 1,Λ (k, Ω) are real and identical (not shown), therefore it is sufficient to show and discuss the diagonal elements ρ Λ,Λ .Finally, the bottom row (c, f, i) illustrates the polariton dispersion relation of the respective structure with exemplary phonon transitions (arrows).The dispersion relations immediately reveal the qualitative differences between the three chosen angles: While for ϑ = 5 • there is a large overlap in the energies of the lower (LP) and upper (UP) polariton branches, for ϑ = 1 • the LP dispersion is essentially flat in the exciton region, resulting in a pronounced energy gap between LP und UP.The intermediate angle ϑ = 3.5 • marks the transition between these two limiting cases, where the top of the LP branch is approximately at the same energy as the bottom of the UP branch.
We start with the discussion of the spectra at ϑ = 5 • [Fig.4 (a)].Here, in the spectra we find two main peaks, directly reflecting absorption by the UP (right peak) and the LP (left peak).The peak position of the LP slightly depends on temperature, while its width is mainly unaffected.In contrast, the peak position of the UP remains essentially constant, but it gets significantly broader with increasing temperature.Furthermore, we find a continuous background extending from the upper peak to lower energies, highlighted in the inset.All these effects are induced by phonon-assisted intra-and inter-band transitions and can be explained using the gPSD in Fig. 4 (b) in combination with the polariton dispersion in (c).
The gPSD for the UP, ρ 11 [solid blue line in Fig. 4 (b)], is a broad distribution that abruptly starts at ℏΩ ≈ −10 meV, slowly decays until a sharp cut-off at ℏΩ ≈ 26 meV, and exhibits two distinct peaks at ℏΩ ≈ 15 meV and 20 meV.The left peak vanishes at low temperatures, which will be discussed later.The gPSD for the LP, ρ 00 (dashed red line), has the same shape as the one of the UP but it is shifted by approximately 20 meV = 2λ 1,k to higher energies (not shown entirely).At first glance, in comparison to the spectra in (a), the energetic ordering of the gPSD contributions seems counterintuitive.However, it can be understood by having a look at the phonon processes entering the gPSD.
As can be seen from the delta-distribution in its definition in Eq. 15, ρ Λ,Λ ′ (k, Ω) describes phonon-induced transitions from an initial polariton state with energy ℏϵ Λ ′ ,k to all possible final states with energies ℏϵ Λ ′′ ,k−Q−G under the emission or absorption of a phonon with energy ℏΩ j,±Q and Ω describes the deviation from energy conservation, i.e., the energy mismatch of the process, in agreement with the finding in Eq. 16a that Ω = 0 corresponds to real, i.e., energy-conserving transitions, which survive in the Markov limit.Based on this interpretation, exemplary phonon-induced intra-and interband transitions starting from k = 0 are sketched in Fig. 4 (c).The phonon emission or absorption is sketched by dotted arrows with a slope determined by the phonon dispersion relation, i.e., by the sound velocity.Dotted arrows with a downward slope represent phonon emission processes and those with upward slope phonon absorption.Arrows that are 'reflected' at the boundary of the MBZ are umklapp processes which will have an impact on the gPSD in particular for small Brillouin zones, where the involved phonon wave vector Q exceeds the 1st MBZ.The energy mismatch ℏΩ is indicated by a vertical, solid or dashed arrow.A vanishing energy mismatch, i.e., ℏΩ = 0, corresponds to energy-conserving phonon transitions, while ℏΩ ̸ = 0 describes off-resonant transitions which are in general suppressed due to the 1/Ω scaling in Eq. 14.The blue (red) arrows correspond to transitions starting at the γ point of the UP (LP), where the small horizontal offset of the arrow is chosen for a better visibility.
With this in mind we can now explain the general shape of the gPSD in Fig. 4 (b), starting with the gPSD for the UP (solid blue lines).The abrupt onset at ℏΩ ≈ −10 meV and end at ℏΩ ≈ 26 meV can be traced back to the finite energy width of the LP band.The lower cut-off directly corresponds to transitions starting from the γ point of the UP and ending in the exciton region at the energy 0 meV of the LP.The high energy cut-off stems from processes which end at the edge of the 1st MBZ (κ point) in the LP branch.The energy spread of the moiré exciton dispersion thus determines the width of the gPSD.Strictly speaking, transitions to the polariton region of the LP are also possible, and they are in fact present in the region between −20 meV and −10 meV.However, due to the large curvature of the LP branch and the resulting very small density of states in the polariton region these contributions have a negligible weight in the gPSD.Consequently, the general shape of the gPSD is dominated by transitions to the exciton part.
Energy-conserving transitions, i.e., transitions with Ω = 0, are possible from the bottom of the UP to the LP branch, both by normal and by umklapp processes, as indicated by the dotted lines on the left side, starting from the bottom of the UP branch.The non-vanishing value of the gPSD of the UP at ℏΩ = 0 gives rise to a strong phonon impact on the UP resonance in the spectrum [Fig.4 (a)], resulting in its pronounced broadening with increasing temperature.
We now come back to the peaks in the gPSD seen in Fig. 4 (b) at ℏΩ ≈ 15 meV and 20 meV.They describe transitions to the flat part of the LP branch close to the edge of the 1st MBZ, namely to the m point of the LP, while absorbing or emitting a phonon.This leads to van Hove singularities in the gPSD, but at this point the processes are strongly energy non-conserving and therefore they do not have a significant impact on the spectrum in Fig. 4 (a).From the asymmetry in the relative peak heights in (b) for small temperatures we can conclude, that the higher energetic peak, i.e., the one characterized by ℏΩ = 20 meV in (c), mainly belongs to phonon emission and the lower one labeled by ℏΩ = 15 meV to phonon absorption processes.In the limit T → 0 K phonon absorption processes are not possible anymore, while phonon creation remains.
When comparing the gPSD of the UP, ρ 11 (solid blue line), with the one of the LP, ρ 00 (dashed red line), we directly see that they mainly differ by their energetic positions.Because transitions to the polariton part are again negligible due to the small density of states, the gPSD essentially starts at an energy of 10 meV with transitions to the bottom of the exciton part of the dispersion and it ends at an energy of 46 meV with transitions to the κ point (not shown).Similar to the case of the UP gPSD, peaks reflecting the van Hove singularities appear, now at energies ℏΩ ≈ 35 meV and 40 meV and therefore outside the plotted energy range.
While the gPSD of the UP is non-vanishing around ℏΩ = 0, i.e., energy-conserving phonon transitions are possible, the gPSD of the LP is shifted to higher energies, such that only energy non-conserving transitions appear in the spectrum.These energy nonconserving transitions give rise to the energy renormalization [see Eq. 16b], i.e., to the temperature-dependent shift of the LP peak in the spectrum in (a).This explains why the LP peak exhibits mainly an energy shift while the UP peak shows a prominent broadening.Energy non-conserving transitions also lead to the continuous background at energies below the UP peak in the spectrum, representing a phonon sideband.Its origin and its specific shape will be discussed later.
The spectra for a twist angle of ϑ = 3.5 • are shown in Fig. 4 (d).Overall, we find the same features as before, but the UP peak is broader than in the spectra for ϑ = 5 • and some additional structures emerge in the phonon sideband region, as highlighted in the inset.When investigating the corresponding gPSD in (e) we again find two pronounced peaks for the UP (solid blue lines) and LP (dashed red lines), reflecting transitions to the van Hove singularities at the m points, which are now shifted to lower energies due to the reduction of the size of the MBZ and the decreased curvature of the moiré exciton dispersion.This can also be seen in the polariton dispersion in Fig. 4 (f), where the energy of the LP at the edge of the MBZ is much smaller compared to the band structure for ϑ = 5 • and it is now close to the energy of the UP at the γ point.Most of the gPSD peaks are energy non-conserving, such that they have a relatively small impact on the spectra.They lead to the small structures in the continuous phonon sideband below the UP peak and contribute to the shift of the polariton peak positions.However, the peak in the gPSD of the UP at ℏΩ ≈ 0 is energy-conserving and therefore leads to the strong broadening of the UP peak in the spectra in Fig. 4 (d).Note that this peak in the gPSD emerges from a phonon emission process [highest dashed blue arrow in (f)], explaining the significant spectral broadening of the UP peak even at low temperatures.
For almost entirely flat bands at ϑ = 1 • the spectra are shown in Fig. 4 (g).We find that the UP and LP peaks are not significantly broadened and the continuous phonon background is replaced by one peak for T = 4 K, two for T = 70 K, and three for T = 200 K in the middle between the polariton peaks.In the gPSD in (h) we do not find energy-conserving transitions anymore, i.e., we have vanishing values at ℏΩ = 0, only a double peak structure centered at approximately −10 meV for the UP (blue) and around +10 meV for the LP (red).The transitions leading to these peaks are sketched in the corresponding polariton dispersion in (i).Due to the opening of a polariton band gap, energy-conserving transitions to the exciton part are not possible anymore, neither for the LP nor for the UP.This explains the small width of the UP peak in (g) compared to the previous twist angles and the dominance of the energy shifts for both peaks.
The peak structure in the phonon sideband region of the spectrum between the two polariton peaks can also be traced back to the gPSD.It can be understood best on the basis of the perturbative solution described in the following section.

Comparison with the Markovian limit and a perturbative solution
In the previous section the spectra have been obtained from a numerical solution of Eq. 13.There, we have seen that the phonon impact on the moiré exciton-polaritons depends non-trivially on the twist angle.On the basis of this solution, together with the corresponding gPSD, we were able to clearly interpret important features like the energy shifts and the broadenings of the polariton peaks.However, other features, in particular the detailed structure of the phonon sidebands turned out to be not so obvious.
The basic equation ( 13) of our model is obtained in the TCL scheme; it is time-local but it is still non-Markovian.In Sec.2.3 we have derived the Markovian limit and we have seen that -as usual -this limit gives rise to energy shifts and broadenings of the spectral lines.In the following we will use the Markovian limit as the starting point for a perturbative solution of Eq. 13.We will then compare the exact numerical results with various orders of the perturbative solution, allowing us to obtain a better understanding of the role of non-Markovian effects and, in particular, of the sideband structure in the spectra.
Starting point for the perturbative solution is the separation of the time-dependent coefficients Γ Λ,Λ ′ (k, t) from Eq. 14 into a time-independent Markovian part given by Eq. 2.3 and a time-dependent non-Markovian correction δΓ Λ,Λ ′ (k, t) according to This leads to the equation of motion where includes the bare polariton part and the Markovian part of the phonon interaction.
Neglecting the non-Markovian part, Eq.21 is a linear differential equation with timeindependent coefficients in an (N + 1)−dimensional space, which for sufficiently small values of N can be solved analytically by a suitable rotation in that space, leading to Lorentzian peaks in the rotated basis.Thus we can deduce, that the Markovian part Γ Markov Λ,Λ ′ leads to a polaronic energy shift due to its imaginary part, a dephasing due to its real part, and a mixing of the polariton branches due to its off-diagonal nature.The non-Markovian part δΓ Λ,Λ (k, t), which will be treated as a perturbation, then dynamically includes corrections to the Markovian line shifts and line broadenings and it gives rise to phonon sidebands [35].In the following we will give a brief overview of the main steps in the derivation of this perturbative approach.Details are given in Appendix C.
First, a change of basis is performed, which diagonalizes the Markovian part and allows us to obtain an analytical solution of this part.This analytical solution corresponds to the 0th order in the perturbation.Afterwards we perform a Laplace transform [see C. As already mentioned, in the 0th order, i.e., in the Markovian limit, the excitonphonon coupling leads to polaron-shifted polariton peaks broadened by the dephasing due to energy-conserving, phonon-induced transitions.The energy renormalization mainly depends on the area under the function ρ Λ,Λ ′ (0, Ω)/Ω while the dephasing is determined by ρ Λ,Λ ′ (0, 0).The non-Markovian part then enters in the nth order with n ≥ 1 and involves a convolution of ρ Λ,Λ ′ (0, Ω)/Ω with the solution in the (n − 1)th order and a correction term involving ρ Λ,Λ ′ (0, 0).
In Fig. 5 we compare the spectra obtained from the numerical solution of Eq. 13 (dashed black lines) with the ones obtained by the perturbative solution according to C.3 in increasing perturbation orders (solid blue, orange and red lines).The temperatures T = 4 K (top row) and T = 200 K (bottom row) and twist angles of ϑ = 5 • (a, d), ϑ = 3.5 • (b, e), and ϑ = 1 • (c, f) are the same as in Fig. 4.
We find that at low temperature (T = 4 K, top row) both polariton peaks are already very well reproduced by the 0th order, i.e., by the Markovian limit (blue solid line) for all twist angles.In particular, for a large angle of ϑ = 5 • [Fig.5 (a)] all curves are essentially indistinguishable showing that non-Markovian effects are completely negligible here.For ϑ = 3.5 • [Fig.5 (b)] the Markovian result slightly underestimates the height of the UP peak and overestimates the spectrum between the two peaks, but already the 1st order corrections lead to a perfect agreement with the numerical result.Higher order corrections turn out to be negligible here.For ϑ = 1 • [Fig.5 (c)] the peaks are again well reproduced in the Markovian limit, but the full solution exhibits a weak additional peak between the two polariton peaks.This is a phonon sideband and it reflects the absorption of a photon under the simultaneous emission of a phonon.It can be understood from the corresponding gPSD, shown in Fig. 4 (h) (light curves).We see that the gPSD of the UP (blue) exhibits a peak at ℏΩ ≈ −8 meV and the gPSD of the LP (red) a peak at ℏΩ ≈ 12 meV.Since the 1st order correction is given by a convolution of the gPSD with the 0th order spectrum, these peaks lead to a phonon sideband of the UP about 8 meV below the main UP peak and a phonon sideband of the LP about 12 meV above the main LP peak.Since the main peaks are separated by 20 meV, the two sidebands overlap resulting in the single additional peak in Fig. 5 (c) at ℏω ≈ 2 meV.Again, the 1st order correction already leads to a perfect agreement with the numerical solution, showing that higher order corrections are negligible.The situation changes at the higher temperature of T = 200 K (bottom row), where we find clear deviations of the Markovian result (blue lines) from the results of the full numerical solution (black dashed lines).At a large twist angle of ϑ = 5 • [Fig.5 (d)] the UP and the LP peak, in particular their asymmetry, are well reproduced already in the Markovian limit.Between these peaks a weak but clearly visible broad phonon sideband appears, which reflects the broad gPSD in Fig. 4 (b) (dark curves) with the rather sharp onset at −10 meV for the UP and 10 meV for the LP.Again, the sidebands of LP and UP overlap; they are characterized by an onset at ℏω ≈ 0 and they are completely described by the 1st order correction with vanishing higher order contributions.
At an intermediate twist angle of ϑ = 3.5 • [Fig.5 (e)] we find that the Markovian approximation strongly overestimates the broadening of the UP peak.In fact, the gPSD in Fig. 4 (e) (dark curves) is characterized by a rather large value of ρ 1,1 (0, 0) leading to a strong dephasing.This is reduced by non-Markovian effects, which leads to an increase of the spectrum at the UP peak and a decrease at the LP peak.Furthermore, a sideband with an onset at ℏω ≈ 0 and a peak at ℏω ≈ 6 meV reflecting the shape of the gPSD shows up.In this case of a quite efficient exciton-phonon coupling, even the 3rd perturbation order is not sufficient to perfectly reproduce the numerical solution.
Finally, at a small twist angle of ϑ = 1 • [Fig.5 (f)] the positions and the widths of the main peaks are well reproduced by the 0th order, i.e., the Markovian approximation.However, we find a redistribution of oscillator strength to the sideband region leading to a reduced height of the main peaks.The sideband now exhibits a three peak structure, which is qualitatively well described in 1st order and quantitatively reproduced in 3rd order perturbation theory.At first sight this seems surprising since the gPSD in Fig. 4 (h) (dark curves) exhibits a two-peak structure both for the UP and the LP, centered around ∓10 meV.However, we notice that the main peaks in the spectrum now exhibit pronounced shifts by about ±1.5 meV for the UP and LP, respectively.Therefore, the sidebands corresponding to the LP and the UP do not completely overlap anymore completely.Instead, the higher energy peak of the LP sideband overlaps with the lower energy peak of the UP, resulting in the three-peak structure with an enhanced middle peak.

Twist angle dependence
In the previous sections we have studied the phonon influence on the optical absorption spectra of moiré exciton polaritons for several twist angles and temperatures.We have seen that the phonon coupling essentially modifies the spectra with respect to three key features: (i) it leads to a temperature dependent broadening of the polariton peaks; (ii) it gives rise to temperature dependent shifts of the peaks; and (iii) it leads to the appearance of phonon sidebands in the spectra, which manifest themselves either in new peaks or shoulders in the region between the main polariton peaks.We were able to clearly interpret these features on the basis of our TCL approach which revealed that the broadening and the shift of the peaks are in most cases well described by the Markovian limit, while the phonon sidebands as well as corrections to the peak heights and widths are consequences of the non-Markovian contributions to the equations of motion.In this section we take a more general view on the twist angle dependence of the phonon influence.First, we analyze the deviations between the perturbative and the numerical solution and then we examine the twist angle dependence of width and position of the LP and UP peaks in the spectrum.
In Fig. 6 we plot the normalized deviation of the nth perturbation order compared to the numerical solution, calculated by as a function of the twist angle for several temperatures and perturbation orders.We can clearly distinguish three regimes of twist angles: In the regimes of large (ϑ ≳ 4.5 • ) or small (ϑ ≲ 2.5 • ) twist angles the deviation of a given order n is essentially independent of the twist angle while in the intermediate regime (2.5 • ≲ ϑ ≲ 4.5 • ) pronounced variations as a function of the twist angle appear.
In the regime of large twist angles (ϑ ≳ 4.5 • ) already the 0th order (solid blue), i.e., the Markovian limit, is in general a good approximation.Only for the highest temperature there are some deviations in the 0th order, which are essentially removed by the 1st order (dashed orange).The reason is that in this regime the exciton band width is larger than the polariton band gap, such that there are always energy-conserving transition between the UP and the LP branch.These reflect the main phonon influence and are typically well described in the Markovian limit.Only at higher temperatures a phonon sideband builds up, which is then well described in the 1st order.
At small twist angles (ϑ ≲ 2.5 • ) the deviations of the Markovian limit from the full result are larger, in particular at higher temperatures, which can be traced back to pronounced phonon sidebands in this regime.However, the perturbative approach converges rather quickly and in 3rd order (dotted red) the deviations from the numerically obtained result become negligibly small for all temperatures.
In the intermediate range of twist angles (2.5 • ≲ ϑ ≲ 4.5 • ) pronounced deviations of the Markovian limit occur and also the convergence of the perturbative solution is less good.In particular, at higher temperatures even the 3rd order still exhibits deviations of up to almost 30%.In this regime the impact of phonons on the spectrum is particularly strong because now energy-conserving transitions to the van Hove singularities are relevant.Due to these efficient scattering channels also the non-Markovian terms are rather large resulting in the poor convergence of the perturbative approach.Finally, we analyze the phonon impact on the width and position of the UP and LP peaks.Without coupling to phonons, these are two Lorentzian peaks at ±10 meV with a width of about 0.8 meV, determined by the background dephasing 2ℏγ 0 .Figure 7 shows the phonon impact on the peak width in (a) and the peak position in (b) of the UP (blue) and LP (red) peak as a function of the twist angle, again for different temperatures (with bright colors for low and dark colors for high temperatures).
We can clearly distinguish the same three regimes of twist angles as in Fig. 6: In the regimes of large (ϑ ≳ 4.5 • ) or small (ϑ ≲ 2.5 • ) twist angles the peak widths and peak positions are essentially independent of the twist angle while in the intermediate regime (2.5 • ≲ ϑ ≲ 4.5 • ) pronounced variations of both quantities as functions of the twist angle are seen.
For small twist angles (ϑ ≲ 2.5 • ) we find a negligible phonon influence on the peak width.The peak width is independent of both twist angle and of temperature; it is essentially given by the constant background dephasing rate 2ℏγ 0 .Due to the flat exciton band and the resulting energy gap in the polariton dispersion relation no real, i.e., energy-conserving, phonon-induced transitions are possible.In contrast, the peak positions exhibit pronounced, temperature dependent shifts leading to an increasing separation between LP and UP peak with increasing temperature.This can be understood from the gPSD in Fig. 4 (h).The UP and the LP contribution are almost mirror symmetric, the UP part being located entirely at negative energies and the LP part at positive energies.Therefore, the Markovian energy shifts according to Eq. 16b are of similar size but with opposite sign, as is indeed seen in Fig. 7 (b).
In contrast, for large twist angles (ϑ ≳ 4.5 • ) we find a strongly temperature dependent width of the UP peak (a) and energy shift of the LP peak (b), while the impact on the width of the LP and the shift of the UP peak is negligible.Here, energyconserving transitions from the bottom of the UP branch to the LP branch are always possible while from the bottom LP branch no such transitions to regions with large density of states are possible.On the other hand, as can be seen in Fig. 4 (b), the gPSD of the LP is non-vanishing only for positive energies leading to the observed energy shift, while the gPSD of the UP is almost symmetric around Ω = 0, which leads to efficient cancellations in Eq. 16b.
At intermediate twist angles (2.5 • ≲ ϑ ≲ 4.5 • ) we find a significantly increased energy shift and broadening of the UP peak, while the LP peak is only weakly affected.This is caused by the pronounced impact of phonon transitions to the boundary of the 1st MBZ, giving rise to van Hove singularities at or close to Ω = 0 in the gPSD.At the same time, the gPSD is still very asymmetric with respect to Ω = 0, leading to a strong energy shift.In the recently introduced language of twisted 2D materials one might therefore call this regime the magic angle for phonon scattering, which becomes particularly prominent here.In Appendix D we present further quantities, namely the gPSD at Ω = 0 and the phonon sideband strength, that also show the transition from small to large twist angles around ϑ = 3.5 • .

Conclusion
In this paper we have derived a numerical model for the linear absorption spectrum of excitons in a TMDC moiré lattice strongly coupled to an optical cavity, where we found that the impact of phonons depends on several quantities, e.g., the moiré exciton dispersion, moiré exciton wave functions and the exciton-phonon coupling.These quantities can directly or indirectly be adjusted via the twist angle of the heterostructure.In the present paper we focused on the phonon impact associated with the varying moiré exciton dispersion.The phonon influence can be expressed in terms of a generalized phonon spectral density, which turned out to be a useful concept for the understanding and interpretation of the spectra.In a comparison with a perturbative solution we were able to identify Markovian processes as the main origin of the energy renormalizations of the polariton peaks and their broadenings.Non-Markovian processes lead to corrections of the Markovian shifts and broadenings.Most importantly, however, they give rise to phonon sidebands in the spectral region between the polariton peaks.
We found fundamentally different optical spectra for small, intermediate and large twist angles with a strongly increased phonon impact for intermediate values around 3.5 • which one might call magic angle for phonon scattering.For different designs a different value for the magic angle is possible.On the one hand, spectra at small twist angles are governed by the presence of a polariton band gap, suppressing phonon transitions and thus leading to sharp polariton peaks.Between the UP and the LP peak characteristic peaked phonon sidebands appear, reflecting phonon assisted transitions to the flat, excitonic part of the polariton dispersion.On the other hand, spectra for large twist angles are dominated by energy-conserving phonon transitions leading to a broadening of the UP peak.Here, the phonon sidebands are broad structures reflecting the wide range of possible phonon transitions.Around the magic angle for phonon scattering phonon processes are most pronounced.This is caused by energy-conserving transitions to the edge of the 1st MBZ, i.e., into van Hove singularities, which then lead to increased phonon scattering resulting in a pronounced peak broadening of the UP.Due to the efficient coupling to phonons the convergence of the perturbative approach is reduced.
Our results show that the rich interplay between moiré exciton-polaritons, which for small twist angles might be considered as qubit arrays for quantum technological applications [25], and phonons is a key ingredient to fully understand the detected optical response.Using the twist angle in the heterostructure as a control parameter the phonon impact might be promoted or suppressed, depending on the desired functionality of the pursued device.
with the electron (hole) dispersion treated in the effective mass approximation ϵ 2m h ), the energy of the band gap E gap , the effective mass of the electron (hole) m e (m h ) and the interaction potential V Q between electrons and holes.
We start by introducing the homogeneous exciton creation operator Y † l,κ = K Φ l (K − µ e κ)c † K d † κ−K with µ e/h = m e/h /M and total mass M = m e + m h , that creates an eigenstate of the Hamiltonian without the moiré potential, i.e., H free + H e−h , with center of mass wave vector κ.The corresponding annihilation operator is Y l,κ .Note that we have introduced the term homogeneous exciton to distinguish it from the moiré exciton that includes the moiré potential.The one-exciton state Y † l,κ |0⟩ satisfies the eigenvalue equation with ground state (electron-hole vacuum) |0⟩, eigenenergies E l,κ , and quantum number of the relative motion l.This leads to a Wannier equation [45] for the coefficients Φ l , i.e., the exciton wave function in reciprocal space, according to with reduced mass µ = m e m h /M and the energies of the homogeneous excitons These homogeneous excitons additionally experience the moiré potential H m .Assuming only pair excitations of electrons and holes and a small exciton density, the electronic transitions resulting from H m can be approximated by [46] With these approximations the Hamiltonian can be written in terms of homogeneous exciton operators G is a Fourier component of the moiré potential for excitons, resulting from the moiré potentials of electrons and holes and modified by the exciton form factor F l,l ′ (Q).As seen in Eq.A.4, this moiré potential in general couples exciton states that differ by a reciprocal moiré lattice vector G j both within the same exciton band (l = l ′ ) and in different bands (l ̸ = l ′ ).In TMDCs the binding energy of the 1s exciton is typically very large [47,48], such that we can restrict ourselves to the 1s exciton with l = l ′ = 1s and we neglect interband transitions that are induced by the moiré potential, thus assuming a single form factor F (Q).In this case the Hamiltonian can be further simplified to where we dropped all band indices l, l ′ .This interaction is of the same form as, for example, derived in Ref. [6] and describes an exciton in a periodic moiré potential with reciprocal lattice vectors G j .Note that Eq.A.5 can also directly be interpreted as a homogeneous exciton coupled to a moiré potential for excitons [see Eq. 1] As seen in Eq.A.5, the moiré potential only couples excitons that differ by a reciprocal lattice vector of the moiré lattice G j .Each wave vector κ can be uniquely written as a sum of a contribution within the 1st MBZ, k ∈ B m and a reciprocal The first term agrees with the imaginary part of the Markovian limit of Γ Λ,Λ ′ [see Eq. 16b].Using Eqs 2.3 and 20, we thus obtain for the non-Markovian part Based on this separation, we obtain the equation of motion ( 21) with the bare polariton and the Markovian part included in the matrix . This diagonalization is guaranteed to work as long as A Λ,Λ ′ (k) has pairwise distinct eigenvalues [49].Since in general A Λ,Λ ′ is not Hermitian or symmetric, the eigenvalues will be complex.Due to the structure of the equation, we can identify these eigenvalues D Λ (k) with the peak positions (real parts) and line widths (imaginary parts) of the polariton, such that the diagonalization is possible as long as the polaritons have distinct energies.
Via the Laplace transformation the differential equation from C.1 can be converted into a Fredholm equation [44] xk and the rotated gPSD ρ Note that the Markovian part of the system is completely included in f k and the non-Markovian part in the operator K k .Assuming that K k → λK k is a small perturbation, C.3 can be solved via the power series [44] which is an iterative solution since Φ .Both Eqs. 13 and C.3 describe the polarization of the polaritons with wave vector k including the influence of phonons, taking the polariton and phonon dispersion and the coupling strength to the phonons into account.These quantities are all hidden in the gPSD ρ Λ,Λ ′ (k, Ω) and the rotated gPSD ρ ′ Λ,Λ ′ (k, Ω), respectively.The optical absorption spectrum (Eq.9) is then given by a superposition of all Φ We arrive at an approximate solution by only taking a finite number of Φ The first term in C.8 includes the Markovian part of the phonon interaction and leads to peaks at s = −iD Λ (k), which depend on the temperature.The second term includes the non-Markovian part that also gives a contribution to the polariton peaks via the terms ∼ 1 s+iD Λ (k) and also an additional term that depends on the shape of the rotated gPSD dΩ . By inserting f Λ ′ ,k (s + iΩ) from C.4 into the latter term, one can see that it leads to a strong contribution at s = −i [Ω + D Λ ′ (k)] weighted by the gPSD at Ω, resulting in phonon sidebands.
Assuming that the gPSD has localized maxima at Ω j these would narrow to additional peaks at s = −i [Ω j + D Λ ′ (k)] after this first iteration step.Notably, if Ω j ≈ −D Λ ′ (k) this leads to additional peaks at s ≈ 0. GPSD contributions at Ω = 0 would lead to increased corrections at one of the polariton peaks.These features are discussed in Sec.3.2.

Appendix D. Twist-angle dependencies of additional quantities
In Sec.3.3 we investigated the twist-angle dependence of the UP and LP peak position and width and discussed the general shape of these curves.We found a twist angle region with a significantly increased energy shift and broadening of the UP, which we identified as the transition region.We traced these effects back to energy-conserving phonon-assisted transitions to the van Hove singularities.In this section we present the twist-angle dependencies of further quantities that show the same transition between small and large twist angles.
For reference, Fig. D1(a) shows the peak width of the UP and LP (same as Fig. 7(a)).In addition, Fig. D1(b) depicts the gPSD value ρ Λ,Λ (0, Ω) at Ω = 0, and (c) the absorption spectrum integrated from −5 meV to +5 meV, i.e., across the isolated phonon sidebands, for different temperatures.The gPSD at Ω = 0 (b) directly reflects energy-conserving phonon transitions and is the main contribution of the gPSD to the broadening of the polariton peaks.Therefore, we find a similar trend as in (a) with a slightly better resolution of the double peak structure in the transition region.These peaks can be directly traced back to phonon-assisted transitions to van Hove singularities in the polariton dispersion, where one peak corresponds to phonon absorption and the other to phonon emission processes.According to Eqs. (C.1) and (C.3) the gPSD contributes to the spectrum in multiple ways, namely, its value at Ω = 0, the area under the curve ρ Λ,Λ ′ (0, Ω)/Ω and a convolution with this function (see Appendix C).This directly explains the slight differences between Figs.D1(a) and (b).
For small twist angles we find the emergence of phonon sidebands in the spectra in the energy range around ℏω = 0 (see Fig. 4(c)).Otherwise, for large twist angles, the sidebands continuously range from the UP peak to lower energies (see Fig. 4(a)).This transition in the character of phonon sidebands illustrates that the spectral behavior is dominated by energy-conserving transitions at large twist angles and by the existence of a polariton band gap at small twist angles.Figure D1(c) shows the integrated spectrum around the phonon sidebands in the center between UP and LP peak.We choose the energy range from −5 meV to +5 meV, i.e., such that this quantity contains all spectral contributions close to ℏω = 0, but not the main polariton peaks.While at large twist angles the only contribution is given by a fraction of the continuous phonon sideband ranging from the UP peak to lower energies resulting in the plateau at smaller values, at small twist angles it includes the whole phonon sideband leading to the plateau at larger values (see Fig. D1(c)).We again find that the transition between these two flat regions happens around ϑ ≈ 3.5 • .Contrary to (a) and (b), where we find sharp features (peaks) exactly at that twist angle, the transition in (c) is smoother.This supports our findings, that the energy-conserving transitions mainly contribute to the UP peak, while the phonon sidebands at ℏω ≈ 0 result from the formation of the polariton band gap.

Figure 1 .
Figure 1.(a) Sketch of two twisted hexagonal lattices with the resulting moiré lattice vectors R 1,2 .(b) Side view of the moiré bilayer in an optical cavity formed by two DBRs.

Figure 2 .Figure 3 .
Figure 2. Curvature of the lowest moiré exciton band at k = 0 in the direction of the κ (solid blue line) and m point (dashed red line) of the 1st MBZ and exciton energy at the κ (solid green line) and m point (dashed violet line) relative to the energy at the γ point as a function of the twist angle ϑ (bottom scale) or, equivalently, the moiré lattice constant 2π/|G 1 | (top scale).The inset shows a sketch of the 1st MBZ with the definition of the κ and m points.
2], which results in a Fredholm equation (C.3) [44], where the non-Markovian part contributes via a convolution operator K k [see C.5].By treating the non-Markovian part as a small perturbation we obtain an iterative solution of the Fredholm equation [see Eqs.C.6 and C.7].For more details we refer to Appendix C.

Figure 6 .
Figure 6.Normalized deviation of the perturbative method from the numerical solution (see Eq. 23) as a function of the twist angle ϑ for several temperatures and perturbation orders.

Figure 7 .
Figure 7. Twist angle scan of the width (a) and position (b) of the UP and LP peaks for different temperatures.

Figure A1 .
Figure A1.Moiré exciton dispersion relations for twist angles of 1 • , 3.5 • , and 5 • (from left to right).The lowest bands (red lines) have been used for the calculation of the polaritons.

Figure D1 .
Figure D1.Twist angle scan of the peak width (a) (identical to Fig. 7 (a)) and ρ Λ,Λ (0, Ω) at Ω = 0 of the UP and LP peaks for different temperatures.(c) Twist angle scan of the phonon sideband-integrated spectrum for different temperatures.The twist angle ϑ = 3.5 • is marked with a vertical dashed line.