# Fileset

[PhysRevX.14.031025.pdf](https://mdr.nims.go.jp/filesets/bb3327da-2fc7-40c6-b9bd-321152081503/download)

## Creator

[Alexander Steinhoff](https://orcid.org/0000-0002-8134-4344), [Edith Wietek](https://orcid.org/0000-0001-6681-9260), Matthias Florian, [Tommy Schulz](https://orcid.org/0009-0001-9965-3495), [Takashi Taniguchi](https://orcid.org/0000-0002-1467-3105), [Kenji Watanabe](https://orcid.org/0000-0003-3701-8119), [Shen Zhao](https://orcid.org/0000-0002-7342-7887), [Alexander Högele](https://orcid.org/0000-0002-0178-9117), Frank Jahnke, [Alexey Chernikov](https://orcid.org/0000-0002-9213-2777)

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Exciton-Exciton Interactions in Van der Waals Heterobilayers](https://mdr.nims.go.jp/datasets/766185b0-9246-4cb0-b305-bdc8f9ebbaaf)

## Fulltext

Exciton-Exciton Interactions in Van der Waals HeterobilayersExciton-Exciton Interactions in Van der Waals HeterobilayersAlexander Steinhoff ,1,2,* Edith Wietek ,3 Matthias Florian,4 Tommy Schulz ,1,2 Takashi Taniguchi,5Kenji Watanabe ,6 Shen Zhao ,7 Alexander Högele ,7,8 Frank Jahnke,1,2 and Alexey Chernikov 31Institut für Theoretische Physik, Universität Bremen, 28334 Bremen, Germany2Bremen Center for Computational Materials Science, Universität Bremen, 28334 Bremen, Germany3Institute of Applied Physics and Würzburg-Dresden Cluster of Excellence ct.qmat,Technische Universität Dresden, 01062 Dresden, Germany4Department of Electrical Engineering and Computer Science, University of Michigan,Ann Arbor, Michigan 48109, USA5Research Center for Materials Nanoarchitectonics, National Institute for Materials Science,1-1 Namiki, Tsukuba 305-0044, Japan6Research Center for Electronic and Optical Materials, National Institute for Materials Science,1-1 Namiki, Tsukuba 305-0044, Japan7Fakultät für Physik, Munich Quantum Center, and Center for NanoScience (CeNS),Ludwig-Maximilians-Universität München, 80539 München, Germany8Munich Center for Quantum Science and Technology (MCQST), 80799 München, Germany(Received 17 November 2023; accepted 13 June 2024; published 14 August 2024)Exciton-exciton interactions are key to understanding nonlinear optical and transport phenomena in vander Waals heterobilayers, which emerged as versatile platforms to study correlated electronic states. Wepresent a combined theory-experiment study of excitonic many-body effects based on first-principle bandstructures and Coulomb interaction matrix elements. Key to our approach is the explicit treatment of thefermionic substructure of excitons and dynamical screening effects for density-induced energy renorm-alization and dissipation. We demonstrate that dipolar blueshifts are almost perfectly compensated bymany-body effects, mainly by screening-induced self-energy corrections. Moreover, we identify acrossover between attractive and repulsive behavior at elevated exciton densities. Theoretical findingsare supported by experimental studies of spectrally narrow, mobile interlayer excitons in atomicallyreconstructed, h-BN-encapsulated MoSe2=WSe2 heterobilayers. Both theory and experiment show energyrenormalization on a scale of a few meV even for high injection densities in the vicinity of the Motttransition. Our results revise the established picture of dipolar repulsion dominating exciton-excitoninteractions in van der Waals heterostructures and open up opportunities for their external design.DOI: 10.1103/PhysRevX.14.031025 Subject Areas: Condensed Matter Physics,Materials Science, OpticsI. INTRODUCTIONVertically stacked van der Waals heterobilayers withtype-II band alignment host layer-separated, Coulomb-correlated electron-hole pairs forming interlayer excitons(ILXs) with binding energies of more than 100 meV andlifetimes that are often drastically increased in comparisonto excitons within a single layer [1–7]. For a wide range ofelectron-hole densities below the excitonic Mott transition[8,9], ILXs constitute an interacting quantum gas withrenormalized spectral properties. Exciton-exciton inter-actions are of major importance in the context of under-standingon-site repulsion energies and emergingmany-bodystates inmoiré systems [10–14], excitonic transport [15–20],and nonlinear optical response [21–23]. Importantly, theylead to excitation-dependent shifts of the exciton resonancetoward higher energies, in close analogy to the physicsof coupled quantum well systems [24–29]. These energyshifts are broadly used as key observables to assess thestrength of the exciton-exciton coupling and determineexciton densities [30,31].Conventionally, the interaction is described in the semi-classical picture of long-range repulsion between twodipoles representing spatially indirect excitons [32,33].Corresponding theoretical descriptions range from basiccapacitor models [26,29] to correlation-corrected formulas[27,28], microscopic theories beyond the effective bosonic*Contact author: asteinhoff@itp.uni-bremen.dePublished by the American Physical Society under the terms ofthe Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution tothe author(s) and the published article’s title, journal citation,and DOI.PHYSICAL REVIEW X 14, 031025 (2024)2160-3308=24=14(3)=031025(28) 031025-1 Published by the American Physical Societyhttps://orcid.org/0000-0002-8134-4344https://orcid.org/0000-0001-6681-9260https://orcid.org/0009-0001-9965-3495https://orcid.org/0000-0003-3701-8119https://orcid.org/0000-0002-7342-7887https://orcid.org/0000-0002-0178-9117https://orcid.org/0000-0002-9213-2777https://ror.org/04ers2y35https://ror.org/04ers2y35https://ror.org/042aqky30https://ror.org/00jmfr291https://ror.org/026v1ze26https://ror.org/026v1ze26https://ror.org/05591te55https://crossmark.crossref.org/dialog/?doi=10.1103/PhysRevX.14.031025&domain=pdf&date_stamp=2024-08-14https://doi.org/10.1103/PhysRevX.14.031025https://doi.org/10.1103/PhysRevX.14.031025https://doi.org/10.1103/PhysRevX.14.031025https://doi.org/10.1103/PhysRevX.14.031025https://creativecommons.org/licenses/by/4.0/https://creativecommons.org/licenses/by/4.0/picture [34,35], and adaptations of classical approaches todescribe localized moiré excitons [36–38]. Nevertheless,the dominating picture of exciton-exciton interactionsbroadly used in the field of van der Waals heterostructuresremains based on the dipole-dipole repulsion. Onlyrecently, predictions of additional contributions emerged,associated with exchange interaction and Pauli blockadefrom the underlying fermionic substructure [18,39]. Thisstate of affairs strongly suggests to test and revise thedipolar description of the exciton-exciton interaction in vander Waals heterostructures, motivating the development ofa comprehensive theoretical approach and a well-definedexperimental setting.Here, we present a combined theory-experiment study ofdensity-dependent exciton energy renormalizations basedon a systematic many-body description of interactingexcitons as composite particles. We determine both theresulting energy shifts and scattering-induced spectralbroadening. Supported by spectroscopic measurementsof mobile excitons in atomically reconstructed, h-BN-encapsulated MoSe2=WSe2 heterobilayers, we find thatthe classical dipolar description drastically overestimatesthe interaction strength. We show that the dipole-dipolerepulsion is largely compensated by competing many-bodyeffects within the dense exciton gas, including correctionterms due to the fermionic substructure and screening.Several contributions have large absolute values yet oppo-site signs. As a net result, exciton energies shift by only afew meV even at exciton densities close to ionizationthreshold above 1012 cm−2. Moreover, we demonstratethat the interplay of different renormalization effects canresult in an effectively attractive interaction with a density-dependent crossover to the repulsive regime.To obtain quantitative results, we consider a represen-tative case of an H-stacked MoSe2=WSe2 heterobilayerencapsulated in hexagonal boron nitride (h-BN) and buildon a theory for the dense exciton gas based on two-particleGreen functions [34,40]. Key to our method developmentis the combination of this established many-body theorywith material-realistic band structure and Coulomb matrixelement calculations, augmented by frequency-dependentexcitonic screening effects. The starting point of ourcalculations is the Bethe-Salpeter equation in the absenceof photoexcited carriers:ðεeQ−k þ εhk − Eν;QÞΦν;Qðe; h;kÞ−Xq;h0;e0Ve;h;h0;e0Q−k;k;kþq;Q−k−qΦν;Qðe0; h0;kþ qÞ ¼ 0; ð1Þwhich we solve to obtain exciton wave functionsΦν;Qðe; h;kÞ and energies Eν;Q. Here, exciton eigenstatesare classified in terms of a quantum number ν for therelative electron-hole motion and the total exciton momen-tum Q. First-principle band structures εak and dielectricallyscreened Coulomb matrix elements Va;b;b0;a0k;k0;k0þq;k−q are usedin the electron-hole picture, expanding beyond the macro-scopic Coulomb interaction model and parabolic bandsemployed in the literature [18]. In the exciton representa-tion, we formulate a self-consistency equation for renor-malized exciton energies asẼν;Q ¼ Eν;Q þ ReΣðν;Q; Ẽν;Q=ℏÞ; ð2Þwith the frequency-dependent exciton self-energyΣðν;Q;ωÞ ¼ ΣHðν;QÞ þ ΣFðν;QÞ þ ΣPBðν;QÞþ ΣMW; retðν;Q;ωÞ: ð3ÞExplicit expressions as well as a detailed derivationof the various self-energy contributions are given inAppendix A. Energy renormalizations are accompaniedby excitation-induced broadening determined by Γν;Q ¼−ImΣðν;Q; Ẽν;QÞ, where 2Γν;Q corresponds to the full widthat half maximum (FWHM) of the homogeneous excitonlinewidth. The self-energy consists of Hartree (H), Fock (F),Pauli blocking (PB), and Montroll-Ward (MW) terms.The different contributions to exciton-exciton interactiondescribedby the exciton self-energy are schematically shownin Fig. 1. A part of the Hartree term represents classicalelectrostatic interaction between excitons [Fig. 1(a)]. It is ofdipole-dipole type and contains two repulsive and twoattractive contributions. Since like charges reside withinthe same layer in type-II heterobilayers, the repulsive termsdominate over the attractive interlayer terms. Therefore, thenet effect is a pronounced blueshift of exciton energies. Onthe other hand, exchange interaction between the fermionicconstituents of excitons [included in both Hartree and Fockterms, Fig. 1(b)] leads to an attractive interaction as in aFermi gas. Excitons as composite bosons also experience arepulsive bosonic exchange interaction that is sketched inFig. 1(c), represented by the Fock term in Eq. (3). This can beunderstood as the consequence of boson statistics favoringoccupation of the same states that, in turn, leads to strongerCoulomb repulsion.Pauli blocking of the fermionic phase space due to theexciton substructure is illustrated in Fig. 1(d). It acts as asource of the increased energy of the exciton transition,well known from the semiconductor Bloch equations[39,41,42]. Finally, the Montroll-Ward self-energy containsall noninstantaneous contributions of the GW self-energy,describing frequency-dependent screening of bosonicexchange interaction by excitations within the dense excitongas. As discussed further below, this is a particularlyimportant contribution to density-dependent exciton energyrenormalizations that is explicitly treated beyond a morephenomenological description of static excitonic screeningin the literature [18]. The result is a decrease of the excitonALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-2energies, conceptually similar to the excitation-induced bandgap renormalization [41–44].II. EXPERIMENTAL SETUPThe experiments are performed on a MoSe2=WSe2heterobilayer encapsulated in high-qualityh-BN.Weemployfabrication conditions of near-sixty-degree stacking favoringlarge-scale atomic reconstruction [45], thus suppressingformation of moiré patterns to achieve a flat potentiallandscape. The Hhh registry is identified by its characteristicphotoluminescence (PL) signatures, optical selection rules,and associated g factors [20]. The PL spectra feature low-energy triplet- and high-energy singlet-spin configurationswith narrow linewidths of the exciton resonances down to afew meV confirming the sample quality. As localizationof ILXs is suppressed in Hhh reconstructed samples, theexcitons are mobile, exhibiting phonon-limited free dif-fusion at lowest temperatures. This is particularly important,since the residual linewidth is likely to be determined bysmall inhomogeneous broadening in addition to phonon-induced scattering. Experimental demonstration of freeexciton diffusion at 5 K in this type of h-BN-encapsulatedmonolayers [46] and heterobilayers [20] shows that theremaining potential fluctuations occur on comparativelylarge spatial scales. We, thus, stress the absence of evidencefor exciton localization in the studied sample to comparethe obtained results with theoretical calculations based onthe assumption of spatially flat potentials. Furthermore,signatures of exciton-exciton annihilation and repulsionvalidate that the excitons can scatter and interact effi-ciently with each other [20].The measurements are performed in an optical micros-copy cryostat at the temperatures of 5 and 70 K. The sampleis excited by a pulsed, 140 fs Ti:sapphire laser with arepetition rate of 80 MHz and its excitation wavelengthtuned resonantly into the A∶1s state of MoSe2. The excita-tion energy densities ranged from0.2 to 20 μJ=cm2 per pulse.Assuming an absorbance of 11%, as estimated from reflec-tance measurements, we determine the correspondingpeak exciton densities of 6 × 1010 and 6 × 1012 cm−2 (seeAppendix B for details). The PL is spectrally dispersed by aspectrometer and detected by a streak camera and a charge-coupled device for time-resolved and time-integrated mea-surements, respectively.III. RESULTSTo quantify many-body renormalizations due tothe exciton-exciton interaction, we solve Eq. (2) self-consistently for various exciton temperatures and densities.The static renormalizations contained in the self-energy (3)are explicitly given byΣHðν;QÞ ¼Xν0;Q0�ṼðDÞ;ν;ν0;ν0;νQ;Q0;Q0;Q − ṼðXÞ;ν;ν0;ν0;νQ;Q0;Q0;Q�NXν0;Q0 ; ð4ÞΣFðν;QÞ ¼Xν0;Q0�ṼðDÞ;ν;ν0;ν;ν0Q;Q0;Q;Q0 − ṼðXÞ;ν;ν0;ν;ν0Q;Q0;Q;Q0�NXν0;Q0 ; ð5Þ(a) (b) (c) (d) (e)FIG. 1. Contributions to exciton-exciton interaction. (a) Dipole-dipole interaction (dip-dip) between excitons in states ν and ν0. Theprocess is composed of elementary Coulomb interaction between the electron of ν and the electron of ν0, proportional to the matrixelement Vee, and between the electron of ν and the hole of ν0 (Vhh) as well as corresponding electron-hole terms (Veh). The resulting shiftof the exciton energy Eν is positive. (b) Fermionic exchange interaction (ferm X) between electrons or holes that are part of two differentexcitons, resulting in a decrease of the exciton energy. (c) Bosonic exchange interaction (bos X) of the whole excitons, resulting in anincrease of the exciton energy. (d) Phase-space filling (PSF) due to fermionic constituents of exciton ν0 as experienced by exciton ν,which yields an increase of Eν. (e) Screened bosonic exchange interaction (screen) from the exciton Montroll-Ward (MW) self-energy.The interaction process is screened by a momentum- and frequency-dependent excitonic polarization ΠqðωÞ corresponding to scatteringprocesses of surrounding excitons between different internal states, which results in a lowering of exciton energy.EXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-3ΣPBðν;QÞ ¼Xν0;Q0�ṼPB;ν;ν0;ν0;νQ;Q0;Q0;Q þ ṼPB;ν;ν0;ν;ν0Q;Q0;Q;Q0�NXν0;Q0 ; ð6Þwith effective exciton-exciton interaction matrix elementsṼ as defined in Appendix A and exciton populations NXν;Q.Dielectric screening by the h-BN environment is modeledby an effectively isotropic dielectric constant εh-BN ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi4.95 × 2.86p≈ 3.76 [47]. We consider the interlayerexcitons populating the K and K0 valleys of the twolowest conduction and two highest valence bands, respec-tively. The focus on K-valley states is motivated by theirrepresentative nature for interlayer excitons and by theabsence of signatures from momentum-indirect excitons inthe experiment [20]. For the bright spin-singlet states,we obtain a binding energy of 131 meV in good agree-ment with typical experimental values [4]. Finally, welimit the exciton density range to be below the predictedexciton Mott density nMott ≈ 2 × 1012 cm−2 for thismaterial system [20].Numerical results of our calculations are shown in Fig. 2.In Fig. 2(a), we analyze the cumulative effect of allcontributions to exciton-exciton interaction that have beenintroduced in Fig. 1 for the 1s exciton with vanishing totalmomentum Q ¼ 0 and bright spin-singlet configuration.The individual contributions are presented in Fig. 2(b).The repulsive dipole-dipole interaction given by the excitonHartree self-energy without fermionic correction terms[first term in Eq. (4)] yields a linear blueshift of up to30 meV at the highest density close to 2 × 1012 cm−2. Theblueshift is reduced by the combination of fermionicexchange contributions, which are part of Hartree andFock self-energies [second terms in Eqs. (4) and (5)],respectively, and Pauli blocking. The latter originates fromthe filling of phase space by the fermionic constituents ofexcitons, as described by the self-energy in Eq. (6).Fermionic exchange and Pauli blocking, although on theorder of 50 meV each at a density of 1012 cm−2, compen-sate each other to a large extent.Bosonic exchange of excitons as described by the Fockself-energy without fermionic terms [first term in Eq. (5)]yields only a weak additional blueshift on the few-meVscale. The relative weakness compared to the directexciton-exciton interaction can be partly understoodfrom the dependence of the different self-energies onexciton populations. While the Hartree-like renormaliza-tion of a certain exciton state is approximately propor-tional to the total exciton density, ΣH;ðDÞðν;QÞ ≈ṼðDÞQ¼0Pν0;Q0 NXν0;Q0 , the Fock-like term has a strong depend-ence on the momentum and spin distribution of exci-tons, ΣF;ðDÞðν;QÞ≈Pν0;Q0 ṼðDÞjQ−Q0jδseðνÞ;seðν0ÞδshðνÞ;shðν0ÞNXν0;Q0 .dip-dipdip-dipbosscreenferm XfermbosscreenFIG. 2. Exciton energy renormalization induced by exciton-exciton interaction. (a) Cumulative density-dependent renormalization ofthe zero-momentum bright 1s-exciton (spin-singlet) energy at a temperature T ¼ 100 K, subsequently adding dipole-dipole interaction(dip-dip), fermionic exchange interaction and phase-space filling (þ ferm Xþ PSF), bosonic exchange interaction (þ bos X), andscreened bosonic exchange (þ screen). The latter represents the result of the full calculation. (b) Individual contributions to the excitonenergy renormalization corresponding to the cumulative presentation in (a) at an electron-hole pair density of 1012 cm−2. The redshiftdue to fermionic exchange and the blueshift due to phase-space filling compensate to a large extent. (c) Calculated temperature anddensity dependence of energy renormalization for the zero-momentum bright 1s exciton. The result for T ¼ 5 K has been obtained fromextrapolating the high-temperature data. (d) Calculated temperature and density dependence of the zero-momentum bright 1s-excitonbroadening.ALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-4Thus, bosonic exchange is sensitive only to a fraction of thetotal exciton density, which splits into intra- and intervalleyexcitons with like and unlike electron ðseÞ and holespins ðshÞ.Finally, dynamical screening of the Coulomb interactiondue to the presence of a polarizable exciton gas provides asubstantial contribution to the energy renormalization. It isrepresented by the Montroll-Ward self-energy (A71),which contains the noninstantaneous contributions to theexciton GW self-energy, leading to a redshift on the orderof tens of meV. Because of a sublinear increase of theredshift with the density, the cumulative energy renormal-ization is nonmonotonic. Although this self-energy isdetermined by frequency-dependent screening caused bythe polarization of the dense exciton gas, we can under-stand its attractive character in the static limit (A80). Here,the self-energy splits into a contribution that leads to staticscreening of exciton exchange (without fermionic correc-tions) and a Coulomb hole (CH) contribution. The formeryields a reduction of the weak momentum-dependentblueshift discussed above, while the latter is essentially arigid redshift of the exciton dispersion. Since the CH shift isthe dominant contribution, the Montroll-Ward self-energyis more sensitive to the total exciton density than to thedistribution of excitons over the states.To summarize these results in Figs. 2(a) and 2(b), themain competition takes place between the blueshift due todipolar interaction and the redshift induced by excitonicscreening, which are both sensitive to the total excitondensity. This leads to the conclusion that the details of theexciton band structure are likely to be secondary forthe total exciton shift. All in all, we find a remarkablystrong compensation of different renormalization effects,leaving a net exciton blueshift of only a few meV even athigh densities, close to the Mott transition. This density-dependent behavior is consistently found at both low andhigh temperatures, as shown in Fig. 2(c) for spin-brightexcitons. For all temperatures, we find a few-meV redshiftat small and intermediate densities which turns into a few-meV blueshift at high densities. This corresponds to adensity-dependent crossover from an effectively attractiveto repulsive exciton-exciton interaction. With decreasingtemperature, the attractive character of interaction becomesslightly weaker yet also closer to the thermal energies of theexcitons. Moreover, we confirm that similar behavior isexpected for “gray” (spin-triplet) interlayer excitons, pre-sented in Fig. 6 in Appendix A. Since the triplet excitonstate is about 20 meV below the bright one, it is morestrongly populated, which leads to an increase of fermionicand bosonic exchange effects among exciton triplets.As a result, the bright exciton interaction is slightly morerepulsive.Importantly, the noninstantaneous nature of theMontroll-Ward self-energy also results in a lifetime broad-ening of excitonic resonances according to Eq. (A79).Temperature- and density-dependent results are shown inFig. 2(d). Overall, we find a sublinear increase of broad-ening with the exciton density that is more pronounced athigher temperatures. At T ¼ 5 K, we find a FWHMbroadening of about 8 meV at highest densities, slightlylarger than the typical linewidths of the exciton resonancesin these systems, likely dominated by residual inhomoge-neities. At elevated temperatures and densities, however,exciton-exciton scattering causes an increase of excitonlinewidths by tens of meV, which is comparable to phonon-induced broadening [48–52]. Exciton-exciton scatteringcan, thus, provide a substantial contribution to dissipationand dephasing.We also note that exciton-exciton annihilation [53–55]should not strongly contribute to the renormalization of theexciton energies. The associated changes of the excitonlifetimes are in the range of tens to hundreds of picosecondsfor both as-exfoliated and h-BN-encapsulated samples[56]. The determined exciton-exciton annihilation coeffi-cient is 5 × 10−3 cm2=s in the studied MoSe2=WSe2heterobilayer, corresponding to recombination rates below0.1 ps−1 across the density range. The correspondingenergy values are, thus, below 0.1 meV, justifying ourtheoretical approach that does not include this effect.Theoretical results are confirmed by experimental obser-vations presented in Fig. 3. As outlined in Sec. II, we ensurethe following conditions for the experiment-theory com-parison: suppression of long-range disorder [57], absenceof moiré-like potentials from large-area atomic recon-struction [45], and the presence of mobile excitons evenat low excitation densities [20]. All three conditions areimportant to avoid density-induced filling of localized tailstates and demonstrate that the excitons can diffusesufficiently far to interact with each other. The latter isparticularly important and is the main reason for using thesame sample as studied in Ref. [20] to extract density-dependent energy shifts. In the studied case, the excitonsare shown to be mobile as well as subject to exciton-excitonrepulsion and exciton-exciton annihilation, confirmingthe above.A selection of time-integrated PL spectra for differentexciton peak densities at a temperature of 5 K is shown inFig. 3(a). The narrow linewidths allow for the distinctassignment of the PL peaks of charged (IX−T ) and neutral(IX0T) interlayer triplets of Hhh registry [20,45]. A smallblueshift of only a few meV is observed for both excitonspecies up to highest densities in the range of the Motttransition, consistent with the literature for samples withsufficiently narrow linewidths [58]. Simultaneously, therelative PL weight shifts from IX−T to IX0T with increasingdensity. We note that, for sufficiently large spectral broad-ening, this shift of the center of gravity would lead to theappearance of a much larger blueshift. To illustrate that, thespectra are convoluted with a Gaussian (full width at halfmaximum of 15 meV), so that the PL of charged andEXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-5neutral excitons merges into one peak. By concealing themultipeak structure of the emission, the density-dependentblueshift is overestimated in contrast to weak shiftsobtained for the individual peaks.In addition to analyzing time-integrated spectra, we usean alternative method of investigating the effect of densityon ILX via time-resolved PL. Taking advantage of theexciton decay, we relate the time axis to the relative changein exciton density, gaining access to quasi-instantaneousmeasurements of the energy shift at a given density. Theexciton lifetimes on the order of hundreds of picosecondsare substantially longer than relevant timescales of thetheoretically calculated phenomena. A representative time-resolved PL evolution at a nominal temperature of 5 K andpeak electron-hole pair density of 1012 cm−2 is shown inFig. 3(b). At time t ¼ 0, the injected electron-hole pair5 K70 K5 K70 K108 109 1010 1011 1012 1013e-h pair density (cm- -2)FWHM (meV)24680Energy (eV)1.3861.3881.3901.3921.3941.4121.4145 K3 meVhHhMoSe2WSe2Experiment Theory1.35 1.40 1.4570 K5210.2Energy (eV) e-h pair density (1011 cm-2) e-h pair density (1011 cm-2)energy shift (meV)energy shift (meV)(d)PL  (norm, offset)(e) (f)(c)1.3951.3901.3850.10.20.30.40.50Energy (eV)Time (ns)int PLPLnorm(b)0.5 10011.38 1.40 1.42PL (norm, offset)(a)5210.20.10.020.010.002expconvneh(1012cm-2)5 KEnergy (eV)01-1-220 4 128 0 4 1285 K 5 K01-1-22neh(1012cm-2)FIG. 3. Experimental energy shifts of interlayer excitons. (a) Time integrated spectra of K-K0 interlayer excitons in an Hhhreconstructed MoSe2=WSe2 heterobilayer for different exciton peak densities after pulsed excitation resonant with the MoSe2 intralayerresonance at the lattice temperature of 5 K. For comparison, spectra convoluted with a Gaussian (FWHM of 15 meV) are presented tosimulate additional inhomogeneous broadening. (b) Streak camera image of the spectrally and time-resolved PL at a peak density of1012 cm−2. The spectra are normalized at each time step, and the color scale is logarithmic. White dashed lines are guides to the eye,following the shift of the peak energy with time. Corresponding PL transient is presented in the right. (c) (Upper) Density-dependentenergy shift of the interlayer exciton PL peaks at two different temperatures, extracted from time-resolved data. The resonances arelabeled according to the previously identified excitons in the Hhh-reconstructed MoSe2=WSe2 heterobilayer: charged and neutralinterlayer triplets (IX0T and IX−T ) and the neutral interlayer singlet (IX0S). Dashed lines are guides to the eye. (Lower) Correspondingdensity-dependent linewidths of charged and neutral interlayer triplets at T ¼ 5 K. The error bars for the extracted peak energies arebelow 1 meV for each individual data point, and there is a spread of the measured values from the two different analysis schemes andmeasurement series. (d) Selected PL spectra of neutral interlayer triplets and singlets at 70 K. (e) Extracted energy shifts of neutralinterlayer triplets from experiment at 5 K. (f) Corresponding results from the calculations.ALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-6density is set equal to the estimated peak density.Consequently, at later times, the exciton density decreasesproportionally to the decaying PL intensity, allowing us toextract the time-dependent peak position effectively as afunction of density. This is reflected in the shift of the peakpositions toward lower energies with time, as indicated bythe dashed lines in Fig. 3(b). We evaluate the energy peakposition for time steps of 50 ps width and assign thecorresponding densities to values given by the injectedmaximum density and the relative change of PL intensity.The results of both time-integrated and time-resolvedanalysis are summarized in Fig. 3(c), upper, for allmeasured triplet and singlet interlayer exciton states. Thedensity range covers the linear regime up to the estimatedMott threshold of several 1012 cm−2 [9] and beyond. Wenote, however, that the main focus for the experiment-theory comparison is the low-to-intermediate regime ofinteracting, bound excitons below the Mott transition. Theexperiments are performed at a temperature of 5 K, to takeadvantage of the narrow linewidths and maximum PLyield,and at 70 K, where the influence of the interlayer trion inPL is negligible. Representative PL spectra of neutralinterlayer triplets and singlets at 70 K are presented inFig. 3(d). The maximum blueshift is found to be less orequal to 3 meV for all studied exciton species in goodagreement with theory with weak indications of a redshift.A direct comparison between experiment and theory ispresented in Figs. 3(e) and 3(f) for the case of the neutralinterlayer triplets at 5 K. While there are quantitativedifferences between the curves, the overall behavior isvery similar. Most importantly, the experiment does notshow large shifts that would be otherwise expected from thesimplified dipolar capacitor model.The energy shifts are accompanied by a small spectralbroadening at higher densities, similar to theoretically pre-dicted values, as illustrated in Fig. 3(c), lower. We also notethat the above observations are typical for samples withsufficiently narrow linewidths, as shown in SupplementalMaterial [59]. Also, the obtained data are subject toexperimental error resulting in the spread of the measuredvalues. Nevertheless, the experiments consistently revealsmall energy shifts in the meV range, independent of thepresence of charged excitons, observed at both low andelevated temperatures.IV. CONCLUSIONSIn conclusion, we have demonstrated that the establishedpicture of dipolar repulsion of interlayer excitons isinsufficient and severely overestimates energy shifts result-ing from exciton-exciton interaction. Among relevantcontributions, we identify fermionic and bosonic exchangeeffects as well as the phase-space filling and excitonicscreening that lead to renormalizations of comparablemagnitude but opposite signs. The main competition takesplace between the repulsive dipolar interaction and theattractive screening-induced self-energy correction, whichboth essentially depend on the total exciton density. The netresult is an energy shift of only a few meV at excitondensities up to theMott transition.Moreover,wedemonstrateconditions for the emergence of weakly attractive excitonpotentials, predicted to occur at low and intermediatedensities. Interestingly, elevated temperatures yield astronger tendency toward the attractive regime. In addition,the energy shifts are accompanied by scattering-inducedspectral broadening of exciton states with a linear temper-ature dependence. The theoretical predictions are confirmedin experiment by taking advantage of well-defined, atomi-cally reconstructed MoSe2=WSe2 heterobilayers void oflocalization with spectrally narrow resonances and mobileexcitons. The magnitude of the measured energy shiftsagrees with the results of the theoretical calculations withinthe experimental uncertainty. This holds both for low andelevated temperatures and is independent of the presence orabsence of residual doping.The developed understanding of the excitonic inter-actions has wide implications for the interpretation ofoptical and nonlinear transport phenomena as well as theoverall phase diagram of interlayer excitons. It challengesthe common notion of dominant dipole-dipole repulsionand offers a more nuanced approach to understandingexciton-exciton interactions in van der Waals heterostruc-tures. Rendering the extraction of interlayer exciton den-sities from the energy shifts via capacitor model less usefulthan previously assumed, it highlights the importance ofspectral analysis for samples featuring multiple, closelyspaced resonances. Interestingly, one can expect thatindividual contributions to the exciton energy renormali-zation could be tunable using distinct sample geometries,dielectric environments, and external fields. Making use ofthis tunability, it should even be possible to design andswitch between effectively repulsive and attractive inter-action regimes. The implications range from the possiblerealization of local compression and excitonic droplets, topotentially favorable conditions for the formation of macro-scopic many-body states such as superfluids and conden-sates. Natural extensions of the presented approach andresults toward interlayer excitons in moiré and moiré-likesuperlattices would have major consequences for the on-site interaction terms determining the correlations. Overall,exciton-exciton interactions beyond the dipolar regimeshould have a substantial impact on the behavior of denseexcitonic quantum gases on a fundamental level and behighly relevant for exploiting nonlinearities in photonic andoptoelectronic applications.ACKNOWLEDGMENTSWe acknowledge financial support from the DeutscheForschungsgemeinschaft via SPP2244 (Project-IDNo. 443405595), Emmy Noether Initiative (CH 1672/1,Project-ID No. 287022282), SFB 1277 (project B05,EXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-7Project-ID No. 314695032), the Würzburg-DresdenCluster of Excellence on Complexity and Topology inQuantum Matter (ct.qmat) (EXC 2147, Project-IDNo. 390858490), and the Munich Center for QuantumScience and Technology (MCQST) (EXC 2111, Project-IDNo. 390814868), as well as resources for computationaltime at the HLRN (Göttingen/Berlin). M. F. and S. Z.acknowledge support by the Alexander von Humboldtfoundation. K.W. and T. T. acknowledge support fromthe JSPS KAKENHI (Grants No. 20H00354,No. 21H05233, and No. 23H02052) and World PremierInternational Research Center Initiative (WPI), MEXT,Japan. A. H. acknowledges funding by the EuropeanResearch Council (ERC) (Grant AgreementNo. 772195). The authors also thank Mikhail M. Glazovand Sivan Refaely-Abramson for valuable discussions.APPENDIX A: THEORYTo obtain a material-realistic description of renormali-zation effects induced by exciton-exciton interactions, wecombine first-principle band structures and Coulomb inter-action matrix elements with a many-body theory for thedense exciton gas based on nonequilibrium Green func-tions. We essentially follow the derivation given by May,Boldt, and Henneberger [34,40], extending it with respectto (i) the generality of matrix elements and (ii) the inclusionof frequency-dependent screening effects.1. Density functional theory calculations, spin-orbitcoupling, and Coulomb matrix elementsDensity functional theory (DFT) calculations for afreestanding Hhh MoSe2=WSe2 heterobilayer are carriedout using QUANTUM ESPRESSOV.6.6 [60,61]. We apply thegeneralized gradient approximation by Perdew, Burke, andErnzerhof [62] and use an optimized norm-conservingVanderbilt pseudopotential [63] at a plane-wave cutoffof 80 Ry. Uniform meshes (including the Γ point) with18 × 18 × 1 k points are combined with a Fermi-Diracsmearing of 5 mRy. Using a fixed lattice constant ofa ¼ 3.29 Å [64] and a fixed cell height of 35 Å, forcesare minimized below 5 × 10−3 eV=Å. The D3 Grimmemethod [65] is used to include van der Waals corrections.We use RESPACK [66] to construct a lattice HamiltonianH0ðkÞ in a 22-dimensional localized basis of Wannierorbitals (dz2 , dxz, dyz, dx2−y2 , and dxy for Mo and W,respectively, and px, py, and pz for Se) from the DFTresults. We also calculate the dielectric function as well asbare and screened Coulomb matrix elements in the local-ized basis. Spin-orbit interaction is included using an on-site L · S-coupling Hamiltonian, which is added to thenonrelativistic Wannier Hamiltonian:HðkÞ ¼ I2 ⊗ H0ðkÞ þHSOC: ðA1ÞHere, I2 is the 2 × 2 identity matrix in the Hilbert spacespanned by eigenstates j↑i and j↓i of the spin z component(perpendicular to the monolayer). We assume that theCoulomb matrix in Wannier representation is spin inde-pendent and that spin-up and spin-down states are notmixed. Diagonalization of HðkÞ yields the band structureελk and the Bloch states jψλki ¼Pα cλα;kjk; αi, where thecoefficients cλα;k describe the momentum-dependent con-tribution of the orbital α to the Bloch band λ. The Blochsums jk;αi are connected to the localized basis viajk; αi ¼ ð1= ffiffiffiffiNp ÞPR eik·RjR; αi with the number of unitcells N and lattice vectors R. The SOC Hamiltonian isgiven byHSOC ¼ 1ℏ2L̃ · S ¼ 12ℏL̃ · σ ðA2Þwith the Pauli matrices σ ¼ ðσx; σy; σzÞ and the angularmomentum operator provided in Ref. [20].Starting from the density-density-like bare Coulombinteraction matrix elements in the Wannier basis,UαβðqÞ ¼XReiq·RUαββαðRÞ ¼XReiq·Rh0; αjhR; βjUðr; r0ÞjR; βij0; αi ðA3Þand the corresponding (statically) screened matrix elements VαβðqÞ, we obtain an analytic description of Coulombinteraction in freestanding transition metal dichalcogenide (TMD) heterobilayers that can be augmented by screening froma dielectric environment [8,67,68]. The parametrization of Coulomb matrix elements for a freestanding Hhh MoSe2=WSe2heterobilayer is provided in Ref. [20]. Environmental screening can be taken into account according to the Wannier functioncontinuum electrostatics approach [69] that combines a macroscopic electrostatic model for the screening by the dielectricenvironment with a localized description of Coulomb interaction. The macroscopic dielectric function of a freestandingbilayer embedded in a vertical heterostructure is obtained by solving Poisson’s equation [70].We finally compute screened Coulomb matrix elements in the Bloch state representation by a unitary transformationusing the coefficients cλα;k:Vλ;ν;ν0;λ0k1;k2;k3;k4¼ 1AXα;βðcλα;k1Þ�ðcνβ;k2Þ�Vαβk1−k4cν0β;k3cλ0α;k4; ðA4ÞALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-8where k4 ¼ k1 þ k2 − k3 þG due to momentum con-servation and A denotes the crystal area.Inspired by the discussion in Ref. [71], we assume thatBloch states are approximately spin diagonal. We assign adefinite spin to each band according to the dominantcontribution given by the coefficients cλα;k. Furthermore,we make use of the fact that Coulomb interaction isspin conserving, so that we can set Coulomb matrixelements Vλ;ν;ν0;λ0k1k2k3k4to zero if λ and λ0 or ν and ν0 belongto different spins.2. Derivation of a GW self-energy for excitonsWe start from the definition of the exciton Schwinger-Keldysh Green function:Kð1; 10Þ ¼ 1iℏ⟪hk1ðt1Þep1ðt1Þe0†p01ðt01Þh0†k01ðt01Þ⟫T; ðA5Þwhich describes the creation of an electron-hole pair at timet01 and annihilation at time t1. The notation ⟪ · ⟫T is shortfor the statistical average⟪ÂðtÞ⟫T ¼ hTC½SCÂðtÞ�ihSCiðA6Þwith the time evolution operatorSC ¼ TC exp�−iℏZCdt0Hextðt0Þ�ðA7Þand TC being the time-ordering operator on the Keldyshcontour, which is depicted in Fig. 4. The Keldysh timecoordinate t includes the branch index n ¼ � of the con-tour as well as the physical time coordinate t. The super-index 1 stands for electron momentum p1, hole momentumk1, electron band e, hole band h, and Keldysh time t1.The band index also contains the spin quantum number. InEq. (A5), the external HamiltonianHextðtÞ ¼eAXk;k0;a;a0Va0;aext ðk;k0; tÞρa0;aðk0;−k; tÞ ðA8Þcouples the time-dependent external potential Vext to theparticle-density operatorρa0;aðk0;k; tÞ ¼ saa0†k0þkðtÞak0 ðtÞ: ðA9ÞHere, we introduce the convention to assign a positive signsa ¼ þ1 to electrons (a ¼ e) and a negative sign sa ¼ −1to holes (a ¼ h). Whenever densitylike operator pairs ofthe type in Eq. (A9) appear, it is understood that both bandindices describe the same carrier species.In the following, we derive an equation of motion (EOM)for the exciton Green function (A5) by means ofHeisenberg’s equation for operators in the interactionpicture with respect to HextðtÞ:iℏ∂∂tÂðtÞ ¼ ½ÂðtÞ; H�; ðA10Þwhere the Hamiltonian describes the interacting gas ofelectrons and holes:H ¼Xp;eεepe†pep þXk;hεhkh†khk þ 12Xk;k0;qXe;e0;a;a0sesaVe;a;a0;e0k;k0;k0þq;k−qe†ka†k0ak0þqe0k−qþ 12Xk;k0;qXh;h0;a;a0shsaVh;a;a0;h0k;k0;k0þq;k−qh†ka†k0ak0þqh0k−q: ðA11ÞHole energies and wave functions are obtained from the valence-band quantities as εhk ¼ −εvk and chα;k ¼ ðcvα;kÞ�. The EOMfor creation and annihilation operators can be derived using½Â B̂; Ĉ� ¼ ½B̂; Ĉ�þ − ½Â; Ĉ�þB̂; ðA12Þwith the anticommutator ½·; ·�þ. To proceed, we apply the time-ordering operator in Eq. (A5):iℏ∂∂t1Kð1; 10Þ ¼ ∂∂t1nθCðt1; t01Þ⟪hk1ðt1Þep1ðt1Þe0†p01ðt01Þh0†k01ðt01Þ⟫þ θCðt01; t1Þ⟪h0†k01ðt01Þe0†p01ðt01Þep1ðt1Þhk1ðt1Þ⟫o; ðA13Þwhere the step function on the Keldysh contour is introduced as θCðt1; t2Þ ¼ δn1;− þ n1δn1;n2θðt1 − t2Þ to take all fourcombinations of branch indices n1 and n2 into account correctly. We use brackets ⟪ · ⟫ to denote the statistical average (A6)FIG. 4. The Keldysh contour C starts at initial time t0, extendsto ∞ on the upper branch (þ), and returns to t0 on the lowerbranch (−).EXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-9without time ordering. The time derivative of expectation values is evaluated by splitting the exponential in SC according tothe Keldysh time ordering. Then, we obtainiℏ∂∂t1Kð1; 10Þ ¼ n1δn1;n01δðt1 − t01ÞFð1; 10Þ þ ðεep þ εhkÞKð1; 10Þ þ eAXk�Xe00Ve;e00ext ðp1 þ k;k; t1ÞKð1; 10Þ����p1þk;e00−Xh00Vh;h00ext ðk1 þ k;k; t1ÞKð1; 10Þ����k1þk;h00�−Xq;h00;e00Ve;h;h00;e00p1;k1;k1þq;p1−qKð1; 10Þ����k1þq;h00p1−q;e00þXk;q;a;a0;h00shsaVh;a;a0;h00k1;k;kþq;k1−qrð1; 10;k;q; a; a0; tþ1 Þ����k1−q;h00þXk;q;a;a0;e00sesaVe;a;a0;e00p1;k;kþq;p1−qrð1; 10;k;q; a; a0; tþ1 Þ����p1−q;e00ðA14Þwith the phase-space filling factorFð1; 10Þ ¼ ⟪F̂ð1; 10Þ⟫ ¼ δk1;k01δp1;p01δh;h0δe;e0 − δp1;p01δe;e0⟪h0†k01ðt1Þhk1ðt1Þ⟫ − δk1;k01δh;h0⟪e0†p01ðt1Þep1ðt1Þ⟫ ðA15Þand the three-particle Green functionrð1; 10;k;q; a; a0; tÞ ¼ 1iℏ⟪a†kðtþÞa0kþqðtÞhk1ðt1Þep1ðt1Þe0†p01ðt01Þh0†k01ðt01Þ⟫T: ðA16ÞHere, tþ denotes a time that is infinitesimally later on the Keldysh contour than t. We also introduce the notationKð1; 10Þjk1þq;h00ðp1þq;e00Þ for a Green function where the hole state (electron state) of argument 1 is changed. By introducingthe free inverse exciton Green functionK−10 ð1; 10Þ ¼ n1δn1;n01δðt1 − t01Þ�iℏ∂∂t1δk1;k01δp1;p01δh;h0δe;e0 −H0ð1; 10Þ − Vextð1; 10; t1Þ�ðA17Þwith the effective exciton HamiltonianH0ð1; 10Þ ¼ ðεep þ εhkÞδk1;k01δp1;p01δh;h0δe;e0 −Xq;h00;e00Ve;h;h00;e00p1;k1;k1þq;p1−qδk1þq;k01δp1−q;p01δh00;h0δe00;e0 ðA18Þand the external potentialVextð1; 10; t1Þ ¼eAXk�Xe00Ve;e00ext ðp1 þ k;k; t1Þδk1;k01δp1þk;p01δh;h0δe00;e0 −Xh00Vh;h00ext ðk1 þ k;k; t1Þδk1þk;k01δp1;p01δh00;h0δe;e0�;ðA19Þwe bring the EOM (A14) to its integral form:ZCdt2Xk2;p2;h2;e2K−10 ð1; 2ÞKð2; 10Þ ¼ n1δn1;n01δðt1 − t01ÞFð1; 10Þ þXk;q;a;a0;h00shsaVh;a;a0;h00k1;k;kþq;k1−qrð1; 10;k;q; a; a0; tþ1 Þ����k1−q;h00þXk;q;a;a0;e00sesaVe;a;a0;e00p1;k;kþq;p1−qrð1; 10;k;q; a; a0; tþ1 Þ����p1−q;e00: ðA20ÞThe exciton Green function is coupled via Coulomb interaction to three-particle Green functions, which, in turn, arecoupled to four-particle Green functions and so forth. To truncate this hierarchy of EOM, we eliminate the three-particleGreen function in Eq. (A20) by using the functional derivative with respect to the external potential, which is contained inthe time evolution operator:ALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-10δKð1; 10ÞδVb0;bext ðq0;q; tÞ ¼eA�sb1iℏrð1; 10;q0 − q;q; b0; b; tÞ − Kð1; 10Þdb0;bðq0;−q; tÞ�ðA21Þwith db0;bðq0;−q; tÞ ¼ ð1=iℏÞ⟪ρb0;bðq0;−q; tÞ⟫. The integral EOM (A20) is then written asZCdt2Xk2;p2;h2;e2K−1effð1;2ÞKð2;10Þ¼n1δn1;n01δðt1− t01ÞFð1;10Þþ iℏAeXk;q;a;a0(Xe00Ve;a;a0;e00p1;k;kþq;p1−qδKð1;10Þjp1−q;e00δVa;a0ext ðkþq;q;tþ1 Þ−Xh00Vh;a;a0;h00k1;k;kþq;k1−qδKð1;10Þjk1−q;h00δVa;a0ext ðkþq;q;tþ1 Þ):ðA22ÞThe operator K−1effð1; 2Þ follows from K−10 ð1; 2Þ if the external potential Vext is replaced by the effective potential Veff thatalso contains the electrostatic potential of the mean charge density:Va;a0eff ðkþ q;q; t1Þ ¼ Va;a0ext ðkþ q;q; t1Þ þAeXk0;b;b0Va;b;b0;a0k;k0;k0−q;kþqρb;b0 ðk0 − q;q; t1Þ: ðA23ÞWe exploit the chain rule for functional derivatives to reformulate the interaction term on the rhs of Eq. (A22):δKð1; 10ÞδVa;a0ext ðkþ q;q; tþ1 Þ¼ZCdt2Xk2;q2;a2;a02δKð1; 10ÞδVa2;a02eff ðk2 þ q2;q2; t2ÞδVa2;a02eff ðk2 þ q2;q2; t2ÞδVa;a0ext ðkþ q;q; tþ1 Þ: ðA24ÞThe second factor is identified as the inverse dielectric functionε−1ðk2 þ q2;q2; a2; a02; t2;kþ q;q; a; a0; t1Þ ¼δVa2;a02eff ðk2 þ q2;q2; t2ÞδVa;a0ext ðkþ q;q; t1Þ¼ n1δn1;n2δðt1 − t2Þδk;k2δq;q2δa;a2δa0;a02þXk0;b;b0Va2;b;b0;a02k2;k0;k0−q2;k2þq2nLðk0 − q2;q2; b; b0; t2;kþ q;−q; a; a0; t1Þ− iℏdb;b0 ðk0 − q2;q2; t2ÞÞda;a0 ðkþ q;−q; t1ÞoðA25Þwith the density-density expectation valueLðk0 − q2;q2; b; b0; t2;kþ q;−q; a; a0; t1Þ ¼1iℏ⟪ρa;a0 ðkþ q;−q; t1Þρb;b0 ðk0 − q2;q2; t2Þ⟫T: ðA26ÞAlso, by defining the inverse exciton Green function viaZCdt2Xk2;p2;h2;e2K−1ð1; 2ÞKð2; 10Þ ¼ n1δn1;n01δðt1 − t01Þδk1;k01δp1;p01δh;h0δe;e0 ; ðA27Þwe find for the polarization functionΠð3; 10;k4 þ q4;q4; a4; a04; t4Þ ¼δKð3; 10ÞδVa4;a04eff ðk4 þ q4;q4; t4Þ¼ eAZCdt1ZCdt2Xk1;p1;h1;e1Xk2;p2;h2;e2Kð3; 1ÞΓð1; 2;k4 þ q4;q4; a4; a04; t4ÞKð2; 10Þ ðA28Þwith the vertex functionΓð1; 2;k4 þ q4;q4; a4; a04; t4Þ ¼ −AeδK−1ð1; 2ÞδVa4;a04eff ðk4 þ q4;q4; t4Þ: ðA29ÞEXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-11Combining Eqs. (A22), (A24), (A25), and (A28), we arrive at the Dyson equation for the exciton Green function:ZCdt2Xk2;p2;h2;e2K−1effð1; 2ÞKð2; 10Þ¼ n1δn1;n01δðt1 − t01ÞFð1; 10Þ þ iℏXk;q;a;a0ZCdt2Xk2;q2;a2;a02ε−1ðk2 þ q2;q2; a2; a02; t2;kþ q;q; a; a0; t1Þ×ZCdt3Xk3;p3;h3;e3�Xe00Ve;a;a0;e00p1;k;kþq;p1−qKð1; 3Þjp1−q;e00 −Xh00Vh;a;a0;h00k1;k;kþq;k1−qKð1; 3Þ����k1−q;h00�×ZCdt4Xk4;p4;h4;e4Γð3; 4;k2 þ q2;q2; a2; a02; t2ÞKð4; 10Þ¼ n1δn1;n01δðt1 − t01ÞFð1; 10Þ þZCdt4Xk4;p4;h4;e4Σð1; 4ÞKð4; 10Þ ðA30Þwith the self-energy Σð1; 10Þ. Equation (A30) is an exact EOM for the exciton Green function, where the many-bodyinteraction effects are contained in the self-energy or, more specifically, in the inverse dielectric function and the vertexfunction. To solve the many-body problem for the dense exciton gas, approximations have to be applied to the self-energy.We derive the lowest approximation to the vertex function (A29) by neglecting the density terms in the phase-spacefilling factor Fð1; 10Þ and inverting Eq. (A30) to obtain an EOM for the inverse exciton Green function:K−1ð1; 10Þ ¼ K−1effð1; 10Þ − Σð1; 10Þ: ðA31ÞThe derivative of K−1effð1; 10Þ with respect to Veff then yields a deltalike vertex function:Γð3; 4;k2 þ q2;q2; a2; a02; t2Þ ≈ n3δn3;n4δðt3 − t4Þn3δn3;n2δðt3 − t2Þ×nδk3;k4δp3þq2;p4δk2;p3δh3;h4δa2;e3δa02;e4 − δp3;p4δk3þq2;k4δk2;k3δe3;e4δa2;h3δa02;h4o: ðA32ÞInserting this into Eq. (A30), we obtain the self-energy in random-phase approximation (RPA):Σð1; 10Þ ¼ iℏXk;q;q2;a;a0�Xēε−1ðp01;q2; ē; e0; t01;kþ q;q; a; a0; tþ1 Þ×�Xe00Ve;a;a0;e00p1;k;kþq;p1−qKð1; 10Þ����p1−q;e00p01−q2;ē−Xh00Vh;a;a0;h00k1;k;kþq;k1−qKð1; 10Þ����k1−q;h00p01−q2;ē þXh̄ε−1ðk01;q2; h̄; h0; t01;kþ q;q; a; a0; tþ1 Þ×�Xh00Vh;a;a0;h00k1;k;kþq;k1−qKð1; 10Þ����k1−q;h00k01−q2;h̄−Xe00Ve;a;a0;e00p1;k;kþq;p1−qKð1; 10Þ����p1−q;e00k01−q2;h̄ �: ðA33ÞThe matrix elements V contain background screening fromthe filled valence bands as well as from the dielectricenvironment. On the other hand, the inverse dielectricfunction ε−1 describes screening from photoexcited elec-trons and holes. Therefore, the (matrix) product ε−1V canbe identified as the fully screened Coulomb interaction W,which gives the self-energy (A33) the well-known GWform. As main difference to the standard GW self-energy,its excitonic version contains four interaction terms be-tween the electrons and holes constituting two interactingexcitons.Note thatCoulomb interaction between the electronand hole within the same exciton is already included in thefree inverse exciton Green function (A17). However, asdiscussed in Ref. [34], the RPA self-energy does not describeeffects due to the exchange of fermions between two excitons.In the following subsection, we again follow Ref. [34] toderive such fermionic corrections to the GW self-energy.3. Exchange corrections to the GW self-energyInstead of eliminating the three-particle Green function rby means of functional derivatives from the EOM of theexciton Green function (A14), we derive an EOM for r itselfusing an alternative decoupling mechanism. To this end, wemake use of the free inverse exciton Green function (A17):ALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-12ZCdt2Xk2;p2;h2;e2rð1;2;k;q; a;a0; tþ1 ÞK−10 ð2;10Þ¼ZCdt2Xk2;p2;h2;e2rð1;2;k;q; a;a0; tþ1 Þn2δn2;n01δðt2 − t01Þ�iℏ∂∂t2δk2;k01δp2;p01δh2;h0δe2;e0 −H0ð2;10Þ−Vextð2;10; t2Þ�: ðA34ÞSince K−10 acts to the left, we have to apply the adjoint operators in each term. The time derivative−iℏ∂∂t01rð1; 10;k;q; a; a0; tþ1 Þ ¼ −∂∂t01nθCðt1; t01Þ⟪a†kðt1Þa0kþqðt1Þhk1ðt1Þep1ðt1Þe0†p01ðt01Þh0†k01ðt01Þ⟫þ θCðt01; t1Þ⟪e0†p01ðt01Þh0†k01ðt01Þa†kðt1Þa0kþqðt1Þhk1ðt1Þep1ðt1Þ⟫oðA35Þis evaluated along the same lines as for the two-particle Green function (A13). Note that, in Eq. (A35), the time derivative isacting on creation instead of annihilation operators. Multiplying Eq. (A34) with K0ð10; 3Þ from the right, integrating over thesuperindex 10, and using the relationZCdt2Xk2;p2;h2;e2K−10 ð1; 2ÞK0ð2; 10Þ ¼ n1δn1;n01δðt1 − t01Þδk1;k01δp1;p01δh;h0δe;e0 ; ðA36Þwe arrive at the following explicit expression for r:rð1; 3;k;q; a; a0; tþ1 Þ ¼ZCdt01Xk01;p01;h0;e0r0ð1; 10;k;q; a; a0; tþ1 ÞK0ð10; 3Þ þZCdt01Xk01;p01;h0;e0Xk0;q0;a00;a000×�Xh00sh0sa00Vh00;a00;a000;h0k01þq0;k0;k0þq0;k01sð1; 10;k;q; a; a0; t1;k0;q0; a00; a000; t01Þ����k01þq0;h00×Xe00se0sa00Ve00;a00;a000;e0p01þq0;k0;k0þq0;p01sð1; 10;k;q; a; a0; t1;k0;q0; a00; a000; t01Þ����p01þq0;e00�K0ð10; 3Þ: ðA37ÞHere, r0 results from the time derivative of the theta functions in Eq. (A35), for which we findr0ð1; 10;k;q; a; a0; tþ1 Þ ¼ n1δn1;n01δðt1 − t01Þ⟪½a†kðt1Þa0kþqðt1Þhk1ðt1Þep1ðt1Þ; e0†p01ðt1Þh0†k01ðt1Þ�⟫¼ n1δn1;n01δðt1 − t01Þn⟪a†kðt1Þa0kþqðt1ÞF̂ð1; 10Þ⟫þ δp01;kþqδe0;a0⟪a†kðt1Þh0†k01ðt1Þhk1ðt1Þep1ðt1Þ⟫þ δk01;kþqδh0;a0⟪e0†p01ðt1Þa†kðt1Þhk1ðt1Þep1ðt1Þ⟫o:ðA38ÞThe three-particle Green function couples to four-particle Green functions of the typesð1; 10;k;q; a; a0; t;k0;q0; a00; a000; t0Þ ¼ 1iℏ⟪a†kðtþÞa0kþqðtÞhk1ðt1Þep1ðt1Þe0†p01ðt01Þh0†k01ðt01Þa00†k0 ðt0þÞa000k0þq0 ðt0Þ⟫T: ðA39ÞThe expression (A37) is inserted in the EOM for the two-particle Green function (A20), which yieldsZCdt2Xk2;p2;h2;e2K−10 ð1; 2ÞKð2; 10Þ ¼ n1δn1;n01δðt1 − t01ÞFð1; 10Þ þZCdt2Xk2;p2;h2;e2Mð1; 2ÞK0ð2; 10Þ ðA40ÞwithEXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-13Mð1; 10Þ ¼Xk;q;a;a0�Xh̄sashVh;a;a0;h̄k1;k;kþq;k1−q�r0ð1; 10;k;q; a; a0; tþ1 Þ����k1−q;h̄þXk0;q0;a00;a000�Xh00sh0sa00Vh00;a00;a000;h0k01þq0;k0;k0þq0;k01sð1; 10;k;q; a; a0; t1;k0;q0; a00; a000; t01Þ����k1−q;h̄k01þq0;h00þXe00se0sa00Ve00;a00;a000;e0p01þq0;k0;k0þq0;p01sð1; 10;k;q; a; a0; t1;k0;q0; a00; a000; t01Þ����k1−q;h̄p01þq0;e00� þXēsaseVe;a;a0;ēp1;k;kþq;p1−q�r0ð1; 10;k;q; a; a0; tþ1 Þ����p1−q;ēþXk0;q0;a00;a000�Xh00sh0sa00Vh00;a00;a000;h0k01þq0;k0;k0þq0;k01sð1; 10;k;q; a; a0; t1;k0;q0; a00; a000; t01Þ����p1−q;ēk01þq0;h00þXe00se0sa00Ve00;a00;a000;e0p01þq0;k0;k0þq0;p01sð1; 10;k;q; a; a0; t1;k0;q0; a00; a000; t01Þ����p1−q;ēp01þq0;e00� �: ðA41ÞEquation (A40) is not a full Dyson equation, since only the free exciton Green functionK0ð2; 10Þ appears on the rhs. It ratherrepresents the first skeleton diagram of an expansion of the full exciton Green function with respect to M. As discussed inRef. [34], it is meaningful to interpretM as a self-energy and upgrade K0ð2; 10Þ to Kð2; 10Þ. Indeed, it turns out that the four-particle Green functions can be factorized into two-particle Green functions such that the GW self-energy obtained fromfunctional derivatives is reproduced. For example, the first four-particle term in Eq. (A41) can be rewritten asXk;q;a;a0Xk0;q0;a00;a000Xh̄Xh00sasa00Vh;a;a0;h̄k1;k;kþq;k1−qVh00;a00;a000;h0k01þq0;k0;k0þq0;k01sð1; 10;k;q; a; a0; t1;k0;q0; a00; a000; t01Þ����k1−q;h̄k01þq0;h00¼ 1iℏXk;q;a;a0Xk0;q0;a00;a000Xh̄Xh00sasa00Vh;a;a0;h̄k1;k;kþq;k1−qVh00;a00;a000;h0k01þq0;k0;k0þq0;k01× ⟪a†kðtþ1 Þa0kþqðt1Þh̄k1−qðt1Þep1ðt1Þe0†p01ðt01Þh00†k01þq0 ðt01Þa00†k0 ðt0þ1 Þa000k0þq0 ðt01Þ⟫T≈1iℏXk;q;a;a0Xk0;q0;a00;a000Xh̄Xh00sasa00Vh;a;a0;h̄k1;k;kþq;k1−qVh00;a00;a000;h0k01þq0;k0;k0þq0;k01× ⟪a†kðtþ1 Þa0kþqðt1Þa00†k0 ðt0þ1 Þa000k0þq0 ðt01ÞiiT⟪h̄k1−qðt1Þep1ðt1Þe0†p01ðt01Þh00†k01þq0 ðt01Þ⟫T¼ 1iℏXk;q;a;a0Xk0;q0;a00;a000Xh̄Xh00Vh;a;a0;h̄k1;k;kþq;k1−qVh00;a00;a000;h0k01þq0;k0;k0þq0;k01× iℏLðk0 þ q0;−q0; a00; a000; t01;kþ q;−q; a; a0; t1ÞiℏKð1; 10Þ����k1−q;h̄k01þq0;h00≈ iℏXk;q;q0;a;a0Xh̄;h00Vh;a;a0;h̄k1;k;kþq;k1−qðε−1ðk01;−q0; h00; h0; t0þ1 ;kþ q;q; a; a0; t1Þ− n1δn1;n01δðt1 − t01Þδk01;kþqδq;−q0δh00;aδh0;a0 ÞKð1; 10Þ����k1−q;h̄k01þq0;h00: ðA42ÞIn the fourth line, we use the definition (A26) of thedensity-density expectation value. In the fifth line, wereplaced the latter by the inverse dielectric functionaccording to Eq. (A25), neglecting the term quadratic inthe carrier density. The first term in the fifth line reproducesthe hole-hole interaction term in theGW self-energy (A33),while the second cancels the corresponding instantaneous(Fock-like) contribution. We, therefore, note that the four-particle terms in Mð1; 10Þ contain the retarded part of theexcitonicGW self-energy, while any instantaneous Hartree-and Fock-type interaction are described by the r0 terms.Unlike in Ref. [34], we keep the retarded part of the GWself-energy as it is and focus on instantaneous correctionterms contained in r0.ALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-14We evaluate the r0 terms in Eq. (A41) with Eq. (A38) by replacing the general band summation over a and a0 by electronand hole bands and bringing all terms to normal order. Neglecting all expectation values with four electron or hole operators,respectively, the resulting instantaneous exciton self-energy is given byMδð1; 10Þ ¼ n1δn1;n01δðt1 − t01Þ�Mδ;ð1Þð1; 10Þ þMδ;ð2Þð1; 10Þ þMδ;ð3Þð1; 10Þ þMδ;ð4Þð1; 10Þ�ðA43ÞwithMδ;ð1Þð1; 10Þ ¼Xk�Xh00;h000Vh;h00;h000;h0k1;k;kþk1−k01;k01δp1;p01δe;e0⟪h00†kðt1Þh000kþk1−k01ðt1Þ⟫−Xe00;e000Vh;e00;e000;h0k1;k;kþk1−k01;k01δp1;p01δe;e0⟪e00†kðt1Þe000kþk1−k01ðt1Þ⟫−Xh00;h000Ve;h00;h000;e0p1;k;kþp1−p01;p01δk1;k01δh;h0⟪h00†kðt1Þh000kþp1−p01ðt1Þ⟫þXe00;e000Ve;e00;e000;e0p1;k;kþp1−p01;p01δk1;k01δh;h0⟪e00†kðt1Þe000kþp1−p01ðt1Þ⟫�; ðA44ÞMδ;ð2Þð1; 10Þ ¼ −Xk�Xh00;h000Vh;h00;h0;h000k1;k;k01;kþk1−k01δp1;p01δe;e0⟪h00†kðt1Þh000kþk1−k01ðt1Þ⟫ −Xe00Vh;e00;e0;h0k1;k;p01;k01δk;p01þk01−k1⟪e00†kðt1Þep1ðt1Þ⟫−Xh00Ve;h00;h0;e0p1;k;k01;p01δk;k01þp01−p1⟪h00†kðt1Þhk1ðt1Þ⟫þXe00;e000Ve;e00;e0;e000p1;k;p01;kþp1−p01δk1;k01δh;h0⟪e00†kðt1Þe000kþp1−p01ðt1Þ⟫�;ðA45ÞMδ;ð3Þð1; 10Þ ¼Xq�Xh00;h000Vh;h00;h0;h000k1;k01−q;k01;k1−q⟪h00†k01−qðt1Þe0†p01ðt1Þep1ðt1Þh000k1−qðt1Þ⟫−Xe00;h000Vh;e00;e0;h000k1;p01−q;p01;k1−q⟪h0†k01ðt1Þe00†p01−qðt1Þep1ðt1Þh000k1−qðt1Þ⟫−Xh00;e000Ve;h00;h0;e000p1;k01−q;k01;p1−q⟪h00†k01−qðt1Þe0†p01ðt1Þe000p1−qðt1Þhk1ðt1Þ⟫þXe00;e000Ve;e00;e0;e000p1;p01−q;p01;p1−q⟪h0†k01ðt1Þe00†p01−qðt1Þe000p1−qðt1Þhk1ðt1Þ⟫�; ðA46ÞandMδ;ð4Þð1; 10Þ ¼ −Xq�Xh00;h000Vh;h00;h000;h00k1;k01−q;k1−q;k01⟪h00†k01−qðt1Þe0†p01ðt1Þep1ðt1Þh000k1−qðt1Þ⟫−Xk;e00;e000;h000Vh;e00;e000;h000k1;k;kþq;k1−qδp1;p01δe;e0⟪h0†k01ðt1Þe00†kðt1Þe000kþqðt1Þh000k1−qðt1Þ⟫−Xk;h00;h000;e000Ve;h00;h000;e000p1;k;kþq;p1−qδk1;k01δh;h0⟪h00†kðt1Þe0†p01ðt1Þe000p1−qðt1Þh000kþqðt1Þ⟫þXe00;e000Ve;e00;e000;e0p1;p01−q;p1−q;p01⟪h0†k01ðt1Þe00†p01−qðt1Þe000p1−qðt1Þhk1ðt1Þ⟫�: ðA47ÞThe terms Mδ;ð1Þ and Mδ;ð2Þ, being proportional to electron and hole occupancies, correspond to the single-particle Hartreeand Fock self-energies, respectively. As discussed in detail in Ref. [34], the single-particle occupancies can beapproximately projected to the phase space of excitons using the relationsEXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-15h†kh0k0 ≈Xp;ee†ph†kh0k0ep; e†pe0p0 ≈Xk;he†ph†khke0p0 : ðA48ÞA similar projection technique has also been used inRefs. [18,72]. The resulting four-particle expectation val-ues, along with the corresponding terms in Mδ;ð3Þ andMδ;ð4Þ, can be interpreted as exciton Green functionsevaluated on the time diagonal according toiℏKð1; 10Þ���t01¼tþ1¼ ⟪h†k01ðtþ1 Þe0†p01ðtþ1 Þep1ðt1Þhk1ðt1Þ⟫:ðA49ÞWe note that the above projection in phase space assumesthat any photoexcited electron or hole is part of a boundexciton, while no higher-order complexes such as biexci-tons contribute beyond a pure factorization in pairs ofexcitons. This is an assumption that becomes more accu-rate for lower densities. It is consistent with a treatmentof renormalization effects by means of a Dyson-typeequation for the exciton Green function with a self-energyformulated solely in terms of excitons [upgraded Eq. (A40)].The Dyson-type equation establishes a self-consistencywithin the exciton gas but does not introduce any higher-order complexes.4. Real-time self-energy in the exciton representationWe now switch to an exciton representation of theinstantaneous self-energy by introducingGðν1;Q1; t1; ν01;Q01; t01Þ ¼Xk1;p1;k01;p01;e;h;e0;h0Φ�ν1;Q1ðe;p1; h;k1ÞKð1; 10ÞΦν01;Q01ðe0;p01; h0;k01Þ ðA50Þwith exciton wave functions as solutions to the Bethe-Salpeter equation that corresponds to the effective excitonHamiltonian (A18):ðεep þ εhkÞΦν;Qðe;p; h;kÞ −Xq;h0;e0Ve;h;h0;e0p;k;kþq;p−qΦν;Qðe0;p − q; h0;kþ qÞ ¼ Eν;QΦν;Qðe;p; h;kÞ: ðA51ÞThe quantum number of the relative electron-hole motion is denoted by ν, while the total exciton momentum isQ ¼ kþ p.Exciton energies are denoted by Eν;Q. Since the wave functions form a basis of the electron-hole-pair Hilbert space, theexciton representation can be inverted:Kð1; 10Þ ¼Xν1;ν01;Q1;Q01Φν1;Q1ðe;p1; h;k1ÞGðν1;Q1; t1; ν01;Q01; t01ÞΦ�ν01;Q01ðe0;p01; h0;k01Þ: ðA52ÞWe obtain the exciton Hartree self-energy from Eq. (A44) and the first and fourth terms in Eq. (A45) asΣHðν1;Q1; t1; ν01;Q01; t01Þ ¼ iℏn1δn1;n01δðt1 − t01ÞXν2;ν02;Q2;Q02Ṽν1;ν02;ν2;ν01Q1;Q02;Q2;Q01Gðν2;Q2; t1; ν02;Q02; t1Þ ðA53Þwith the effective exciton-exciton interaction matrix elementsṼν1;ν2;ν3;ν4Q1;Q2;Q3;Q4¼ Ṽν1;ν2;ν3;ν4Q1;Q2;Q¼Q3−Q2¼Xk1;k2;h1;h2;e1;e2Φ�ν1;Q1ðe1;h1;k1ÞΦν3;Q2þQðe2;h2;k2Þ�XēΦ�ν2;Q2ðē;h2;k2Þ×�Xe0Ve1;ē;e2;e0Q1−k1;Q2−k2;Q2−k2þQ;Q1−k1−QΦν4;Q1−Qðe0;h1;k1Þ−Xe0Ve1;ē;e0;e2Q1−k1;Q2−k2;Q1−k1−Q;Q2−k2þQΦν4;Q1−Qðe0;h1;k1Þ−Xh0Vh1;ē;e2;h0k1;Q2−k2;Q2−k2þQ;k1−QΦν4;Q1−Qðe1;h0;k1−QÞ þXh̄Φ�ν2;Q2ðe2; h̄;k2−QÞ�Xh0Vh1;h̄;h2;h0k1;k2−Q;k2;k1−QΦν4;Q1−Qðe1;h0;k1−QÞ−Xh0Vh1;h̄;h0;h2k1;k2−Q;k1−Q;k2Φν4;Q1−Qðe1;h0;k1−QÞ−Xe0Ve1;h̄;h2;e0Q1−k1;k2−Q;k2;Q1−k1−QΦν4;Q1−Qðe0;h1;k1Þ �¼ ṼðDÞ;ν1;ν2;ν3;ν4Q1;Q2;Q¼Q3−Q2− ṼðXÞ;ν1;ν2;ν3;ν4Q1;Q2;Q¼Q3−Q2: ðA54ÞALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-16Note that we make use of the momentum conservationrelation p ¼ Q − k implied in the exciton wave functionsto eliminate electron momenta. Terms number 1, 3, 4, and 6stem from the single-particle Hartree term and correspondto the semiclassical dipolar Coulomb interaction betweenexcitons [27], which is described by ṼðDÞ. We find tworepulsive and two attractive contributions caused by theopposite charge of electrons and holes. Terms number 2and 5, however, are nonclassical terms due to the fermionicnature of the exciton constituents: These fermioniccorrection terms describe the exchange of either an electronor a hole between two excitons leading to an attractiveinteraction proportional to ṼðXÞ. Notably, besides the signthe only difference between “direct” and “exchange” termsinside the effective matrix element is the swap of two statesin the elementary Coulomb matrix element.In a similar way, an exciton Fock self-energy is obtainedfrom Eq. (A46) as well as the first and fourth termsin Eq. (A47):ΣFðν1;Q1; t1; ν01;Q01; t01Þ ¼ iℏn1δn1;n01δðt1 − t01ÞXν2;ν02;Q2;Q02Ṽν1;ν02;ν01;ν2Q1;Q02;Q01;Q2Gðν2;Q2; t1; ν02;Q02; t1Þ: ðA55ÞThe exciton Fock self-energy can be interpreted as resulting from the bosonic exchange of two excitons, with only twostates in the effective matrix element swapped with respect to the exciton Hartree term. As the latter, the Fock term containsfermionic exchange terms.The remaining four terms of the instantaneous self-energy, which are the second and third terms in Eqs. (A45) and (A47),respectively, can be grouped together similar to the exciton Hartree and Fock self-energies:ΣPBðν1;Q1; t1; ν01;Q01; t01Þ ¼ iℏn1δn1;n01δðt1 − t01ÞXν2;ν02;Q2;Q02�ṼPB;ν1;ν02;ν2;ν01Q1;Q02;Q2;Q01þ ṼPB;ν1;ν02;ν01;ν2Q1;Q02;Q01;Q2�Gðν2;Q2; t1; ν02;Q02; t1Þ: ðA56ÞWe introduce new effective matrix elementsṼPB;ν1;ν2;ν3;ν4Q1;Q2;Q3;Q4¼ ṼPB;ν1;ν2;ν3;ν4Q1;Q2;Q̄¼Q4−Q2¼Xk1;k2;h1;h2;e1;e2Φ�ν1;Q1ðe1; h1;k1Þ×�Φν4;Q2þQ̄ðe2; h2;k2ÞXh̄;e00Vh1;e00;e2;h2k1;Q2þQ̄−k1;Q2þQ̄−k2;k2Φν3;Q1−Q̄ðe1; h̄;k1 − Q̄ÞΦ�ν2;Q2ðe00; h̄;k1 − Q̄ÞþΦν3;Q1−Q̄ðe2; h2;k2 − Q̄ÞXē;h00Vh00;e1;e2;h2k1−Q̄;Q1−k1;Q1−k2;k2−Q̄Φν4;Q2þQ̄ðē; h1;k1ÞΦ�ν2;Q2ðē; h00;k1 − Q̄Þ�: ðA57ÞThe self-energy contribution (A56) appears on the samelevel as the Hartree-Fock energy renormalizations, is propor-tional to electron-hole interaction matrix elements, and hasthe positive sign of a repulsive interaction. It is, therefore,plausible to attribute the self-energy to Pauli blocking causedby phase-space filling from the fermionic constituents. In aBethe-Salpeter or Wannier equation for the exciton, phase-space filling appears as a prefactor similar to Eq. (A15) of theelectron-hole term, thereby reducing the exciton bindingenergy [42]. Expanding the densitylike terms in the phase-space filling factor according to Eq. (A48) and adding thecorresponding bosonic exchange term, ΣPB can be derived.We note that, unlike the fermionic exchange corrections inΣH and ΣF, ΣPB cannot simply be obtained by swapping twostates in a Coulomb matrix element. Both the fermioniccorrection terms and the Pauli-blocking terms are beyond theRPA self-energy derived in the previous subsection viafunctional derivative technique.We now bring the retarded part of the exciton GW self-energy to a form similar to the above exciton Hartree andFock self-energies. Starting from Eq. (A42), we derive anexpression that contains the fully screened Coulombinteraction W. Based on the arrangement of operators inthe definition (A25), the elements of the dielectric matrixcan be identified asε−1ðk2 þ q2;q2; a2; a02; t2;kþ q;q; a; a0; t1Þ ¼ hψa2k2jhψa0kþqjε−1ðr; r0; t2; t1Þjψakijψa02k2þq2i≈ ε−1;a2;a0;a;a02k2;kþq;k;k2þq2ðt2; t1Þδq;q2δk;k2; ðA58ÞEXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-17in analogy to Coulomb matrix elements. The first Kronecker delta takes momentum conservation into account, while thesecond is an additional approximation. Expanding the dielectric and Coulomb matrix elements in terms of Wannierfunctions and making use of the completeness of Bloch states, we findiℏXk;q;q0;a;a0Xh̄;h00Vh;a;a0;h̄k1;k;kþq;k1−qε−1ðk01;−q0; h00; h0; t0þ1 ;kþ q;q; a; a0; t1ÞKð1; 10Þ����k1−q;h̄k01þq0;h00¼ iℏXq;a;a0Xh̄;h00Vh;a;a0;h̄k1;k01−q;k01;k1−qε−1;h00;a0;a;h0k01−q;k01;k01−q;k01ðt0þ1 ; t1ÞKð1; 10Þ����k1−q;h̄k01−q;h00¼ iℏXqXh̄;h00Xα;α0ðch00α0;k01−qÞ�ðchα;k1Þ�Wα0α−qðt0þ1 ; t1Þch̄α;k1−qch0α0;k01Kð1; 10Þ����k1−q;h̄k01−q;h00¼ iℏXqXh̄;h00Wh00;h;h̄;h0k01−q;k1;k1−q;k01ðt0þ1 ; t1ÞKð1; 10Þ����k1−q;h̄k01−q;h00ðA59Þwith the fully screened Coulomb matrixWαβq ðt; t0Þ ¼Xγε−1;αγq ðt; t0ÞVγβq : ðA60ÞAnalogous calculations are carried out for the three remaining four-particle terms in Eq. (A41). Then, we arrive at theexciton GW self-energy:ΣGWðν1;Q1; t1; ν01;Q01; t01Þ ¼ iℏXν2;ν02;Q2;Q02W̃ν1;ν02;ν01;ν2Q1;Q02;Q01;Q2ðt01; t1ÞGðν2;Q2; t1; ν02;Q02; t01Þ ðA61Þwith the effective matrix elementW̃ν1;ν2;ν3;ν4Q1;Q2;Q¼Q3−Q2ðt01; t1Þ ¼Xk1;k2;h1;h2;e1;e2Φ�ν1;Q1ðe1; h1;k1ÞΦν3;Q2þQðe2; h2;k2Þ×�XēΦ�ν2;Q2ðē; h2;k2Þ�Xe0Wē;e1;e0;e2Q2−k2;Q1−k1;Q1−k1−Q;Q2−k2þQðt01; t1ÞΦν4;Q1−Qðe0; h1;k1Þ−Xh0Wē;h1;h0;e2Q2−k2;k1;k1−Q;Q2−k2þQðt01; t1ÞΦν4;Q1−Qðe1; h0;k1 −QÞ þXh̄Φ�ν2;Q2ðe2; h̄;k2 −QÞ�Xh0Wh̄;h1;h0;h2k2−Q;k1;k1−Q;k2ðt01; t1ÞΦν4;Q1−Qðe1; h0;k1 −QÞ−Xe0Wh̄;e1;e0;h2k2−Q2;Q1−k1;Q1−k1−Q;k2ðt01; t1ÞΦν4;Q1−Qðe0; h1;k1Þ �: ðA62ÞIn summary, we obtain a Dyson-type equation for the exciton Green function by upgrading Eq. (A40):ZCdt2Xk2;p2;h2;e2K−10 ð1; 2ÞKð2; 10Þ ¼ n1δn1;n01δðt1 − t01ÞFð1; 10Þ þZCdt2Xk2;p2;h2;e2Mð1; 2ÞKð2; 10Þ: ðA63ÞThe exciton self-energyMð1; 2Þ contains Hartree- and Fock-type contributions including fermionic corrections and a Pauli-blocking term as well as an RPA-type correlation term stemming from the factorization of four-particle Green functions.ALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-18Consistent with the derivation of the RPAvertex function in Eqs. (A31) and (A32), we drop the densitylike terms in Fð1; 10Þ.Then, using the Bethe-Salpeter equation (A51), Dyson’s equation can be expressed in the exciton representation:ZCdt2Xν2;Q2��iℏ∂∂t1− Eν1;Q1 δν1;ν2δQ1;Q2n1δn1;n2δðt1 − t2Þ− ΣHðν1;Q1; t1; ν2;Q2; t2Þ − ΣFðν1;Q1; t1; ν2;Q2; t2Þ − ΣPBðν1;Q1; t1; ν2;Q2; t2Þ�Gðν2;Q2; t2; ν01;Q01; t01Þ¼ n1δn1;n01δðt1 − t01Þδν1;ν01δQ1;Q01þZCdt2Xν2;Q2�ΣGWðν1;Q1; t1; ν2;Q2; t2Þ − ΣGW;δðν1;Q1; t1; ν2;Q2; t2Þ�Gðν2;Q2; t2; ν01;Q01; t01Þ: ðA64ÞHere, ΣGW;δ denotes the instantaneous part of the exciton-GW self-energy emerging from the factorization of four-particleGreen functions according to Eq. (A42). As described in Refs. [34,40,73], we transform Dyson’s equation to physical timesby unfolding the Keldysh contour ðRC dt2 ¼ Pn2R∞−∞ dt2Þ, dropping the external potential and using the Langreth-Wilkinstheorems. In close analogy to single-particle Green functions, we introduce greater and lesser Green functionsg>ðν1;Q1; t1; ν01;Q01; t01Þ ¼ Gn1¼−;n01¼þðν1;Q1; t1; ν01;Q01; t01Þ;g<ðν1;Q1; t1; ν01;Q01; t01Þ ¼ Gn1¼þ;n01¼−ðν1;Q1; t1; ν01;Q01; t01Þ ðA65Þand retarded Green functionsgretðν1;Q1; t1; ν01;Q01; t01Þ ¼ θðt1 − t01Þ�g>ðν1;Q1; t1; ν01;Q01; t01Þ − g<ðν1;Q1; t1; ν01;Q01; t01Þ�: ðA66ÞRetarded self-energies additionally contain an instantaneous contribution:Σretðν1;Q1; t1; ν01;Q01; t01Þ ¼ Σδðν1;Q1; t1; ν01;Q01; t01Þ þ θðt1 − t01Þ�Σ>ðν1;Q1; t1; ν01;Q01; t01Þ − Σ<ðν1;Q1; t1; ν01;Q01; t01Þ�:ðA67ÞIn the following, we assume that Green functions and self-energies are diagonal in the exciton representation. Equal-timeGreen functions, which appear in the instantaneous self-energy, correspond to g< functions, consistent with the replacement(A49). We further assume that the exciton gas is in a quasiequilibrium state, where g< is a time-independent distributionfunction. We again discard Pauli-blocking nonlinearities by describing excitons as ideal bosons, which impliesiℏg<ðν1;Q1Þ ¼ NXν1;Q1¼exp�Eν1;Q1− μXkBT − 1�−1ðA68Þwith the temperature T and the exciton chemical potential μX. The exciton density is given by nX ¼ ð1=AÞPν;Q NXν;Q.According to the bosonic commutator rules, it isiℏg>ðν1;Q1Þ ¼ 1þ iℏg<ðν1;Q1Þ: ðA69ÞIn the quasiequilibrium case, evaluating iℏð∂=∂t1Þgretðν1;Q1; t1Þ with Eqs. (A66) and (A64), an algebraic equation for theretarded exciton Green function can be derived in the frequency domain:�ℏω − Eν1;Q1− ΣHðν1;Q1Þ − ΣFðν1;Q1Þ − ΣPBðν1;Q1Þ − ΣMW;retðν1;Q1;ωÞ�gretðν1;Q1;ωÞ ¼ 1: ðA70ÞThe retarded Green function describes the spectral properties of the dense exciton gas. Many-body interaction effects aretaken into account via the self-energy ΣðωÞ ¼ ΣH þ ΣF þ ΣPB þ ΣMW;retðωÞ, which acts as a frequency-dependent operator.We introduce a Montroll-Ward self-energy ΣMW;retðωÞ for excitons, which contains only the retarded part of the GW self-energy (A61). Assuming that excitons can be described as quasiparticles with a renormalized energy Ẽν1;Q1and aEXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-19broadening Γν1;Q1, the exciton Montroll-Ward self-energy can be derived along the same lines as in the single-particlecase [8,73]:ΣMW;retðν1;Q1;ωÞ ¼ iℏXν01;Q01Z∞−∞dω02πð1þ NXν01;Q01ÞW̃>;ν1;ν01;ν1;ν01Q1;Q01;Q1;Q01ðω0Þ − NXν01;Q01W̃<;ν1;ν01;ν1;ν01Q1;Q01;Q1;Q01ðω0Þℏω − Ẽν01;Q01þ iΓν01;Q01− ℏω0 : ðA71ÞThe plasmon propagators W̃>=<ðωÞ are given by theeffective fully screened matrix elements (A62), whichare linear combinations of Coulomb matrix elements Win the Bloch basis. Using the unitary transformation toWannier orbitals as shown in Eq. (A59), the excitonMontroll-Ward self-energy can be expressed in terms ofplasmon propagators in the Wannier basis. The latter fulfillthe Kubo-Martin-Schwinger relation [73]W>αβ;qðωÞ ¼ eℏω=kBTW<αβ;qðωÞ; ðA72Þwhich in combination with 2iImWretαβ;qðωÞ ¼ W>αβ;qðωÞ −W<αβ;qðωÞ allows one to directly relate the propagators to theretarded Coulomb interaction:W>αβ;qðωÞ ¼ ½1þ nBðωÞ�2iImWretαβ;qðωÞ;W<αβ;qðωÞ ¼ nBðωÞ2iImWretαβ;qðωÞ ðA73Þwith the Bose distribution function nBðωÞ. The retardedfully screened Coulomb matrix is obtained using theinverse dielectric matrix for photoexcited carriers accordingto Eq. (A60):Wretαβ;qðωÞ ¼Xγε−1;ret;αγq ðωÞVγβq : ðA74ÞThe dielectric matrix itself is given byεret;αβq ðωÞ ¼ δα;β −XγVαγq Πret;γβq ðωÞ ðA75Þand is discussed in the following subsection.The instantaneous part of the self-energy is derived fromEqs. (A53), (A55), and (A56) using the approximation forthe equal-time (lesser) Green functions (A68) and thesplitting of effective matrix elements (A54) into directand exchange parts:ΣHðν1;Q1Þ ¼Xν01;Q01Ṽν1;ν01;ν01;ν1Q1;Q01;Q01;Q1NXν01;Q01¼ ΣH;ðDÞðν1;Q1Þ þ ΣH;ðXÞðν1;Q1Þ; ðA76aÞΣH;ðDÞðν1;Q1Þ ¼Xν01;Q01ṼðDÞ;ν1;ν01;ν01;ν1Q1;Q01;Q01;Q1NXν01;Q01; ðA76bÞΣH;ðXÞðν1;Q1Þ ¼ −Xν01;Q01ṼðXÞ;ν1;ν01;ν01;ν1Q1;Q01;Q01;Q1NXν01;Q01; ðA76cÞΣFðν1;Q1Þ ¼Xν01;Q01Ṽν1;ν01;ν1;ν01Q1;Q01;Q1;Q01NXν01;Q01¼ ΣF;ðDÞðν1;Q1Þ þ ΣF;ðXÞðν1;Q1Þ; ðA77aÞΣF;ðDÞðν1;Q1Þ ¼Xν01;Q01ṼðDÞ;ν1;ν01;ν1;ν01Q1;Q01;Q1;Q01NXν01;Q01; ðA77bÞΣF;ðXÞðν1;Q1Þ ¼ −Xν01;Q01ṼðXÞ;ν1;ν01;ν1;ν01Q1;Q01;Q1;Q01NXν01;Q01; ðA77cÞΣPBðν1;Q1Þ ¼Xν01;Q01�ṼPB;ν1;ν01;ν01;ν1Q1;Q01;Q01;Q1þ ṼPB;ν1;ν01;ν1;ν01Q1;Q01;Q1;Q01�NXν01;Q01: ðA78ÞBecause of the Coulomb singularity at long wavelength [Q ¼ 0 in Eq. (A54)], the Hartree interaction requires a separatetreatment in the Wannier representation. For a charge-neutral system, the macroscopic (leading) term drops out, and onlymicroscopic contributions to Hartree interaction remain. Following the procedure in Ref. [74], we calculate Hartree-typematrix elements by setting the macroscopic eigenvalue of the Coulomb matrix to zero before transforming the matrixALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-20element to the Wannier representation. This means that thedipole-dipole interaction is not screened by the dielectricenvironment but only by the polarizability of the bilayeritself. This observation is consistent with the modeldeveloped in Ref. [18].The quasiparticle approximation for exciton Green func-tions implies that the exciton self-energy Σðν1;Q1;ωÞ¼ΣHðν1;Q1ÞþΣFðν1;Q1ÞþΣPBðν1;Q1ÞþΣMW;retðν1;Q1;ωÞhas to be evaluated self-consistently:Ẽν1;Q1¼ Eν1;Q1þ ReΣðν1;Q1; Ẽν1;Q1=ℏÞ;Γν1;Q1¼ −ImΣðν1;Q1; Ẽν1;Q1=ℏÞ þ Γ0: ðA79ÞSince the instantaneous part of the self-energy is real valued,quasiparticle broadening induced by exciton-exciton inter-action stems only from the Montroll-Ward self-energy. Weadditionally introduce a phenomenological broadening Γ0that takes into account scattering of excitons with phononsand defects. We choose Γ0 ¼ 5 meV independent of thetemperature.Finally, we note that a static limit can be applied to theself-energy similar to the single-particle case [75], whichresults in a screened-exchange-Coulomb-hole (SXCH)self-energy for excitons:ΣF;ðDÞðν1;Q1Þ þ ΣMW;retðν1;Q1; Ẽν1;Q1=ℏÞ ≈ ΣSXCHðν1;Q1Þ¼Xν01;Q01W̃ν1;ν01;ν1;ν01Q1;Q01;Q1;Q01ðω ¼ 0ÞNXν01;Q01þ 12Xν01;Q01�W̃ν1;ν01;ν1;ν01Q1;Q01;Q1;Q01ðω ¼ 0Þ − ṼðDÞ;ν1;ν01;ν1;ν01Q1;Q01;Q1;Q01�: ðA80ÞFrom the structure of the exciton Montroll-Ward self-energy, it follows that (static) screening due to photo-excited carriers is applied to the bosonic exchangeinteraction described by the exciton Fock self-energy.However, the fermionic correction terms to the Fock self-energy are not screened. Also, the exciton Hartree andPauli-blocking contributions remain unscreened. Besidesscreening of excitonic exchange, a Coulomb hole self-energy arises and leads to a redshift of energies.Throughout this work, we take into account the fullfrequency dependence of screening, which is discussed indetail in the following.5. Excitonic screening in RPAConsistent with the inverse dielectric function (A25), thedielectric function itself is given byεðkþ q;q; a; a0; t1;k2 þ q2;q2; a2; a02; t2Þ¼ δVa;a0ext ðkþ q;q; t1ÞδVa2;a02eff ðk2 þ q2;q2; t2ÞðA81Þwith the relation between Vext and Veff given in Eq. (A23).Using the definition of the particle-density operator (A9),we obtainεðkþ q;q; a; a0; t1;k2 þ q2;q2; a2; a02; t2Þ¼ n1δn1;n2δðt1 − t2Þδk;k2δq;q2δa;a2δa0;a02−AeXk0δδVa2;a02eff ðk2 þ q2;q2; t2Þ�Xe;e0Va;e;e0;a0k;k0;k0−q;kþq⟪e†k0 ðt1Þe0k0−qðt1Þ⟫ −Xh;h0Va;h;h0;a0k;k0;k0−q;kþq⟪h†k0 ðt1Þh0k0−qðt1Þ⟫�:ðA82ÞConsistent with the derivation of Hartree-type self-energy terms, we proceed by expanding the single-particle densities interms of exciton densities according to Eq. (A48) and interpreting the resulting four-particle expectation values as excitonGreen functions. This means that Eq. (A82) contains functional derivatives of the exciton Green function with respect to theeffective potential, which we have identified before as the polarization function; see Eq. (A28). Thus, the dielectric functionhas the well-known form ε ¼ 1 − VΠ with an excitonic polarization function Π. The polarization function is calculated inRPA by inserting the deltalike vertex function (A32):Πð3; 4;k2 þ q2;q2; a2; a02; t2Þ ≈eA�Xk1;h1Kð3; a2;k2; h1;k1; t2ÞKða02;k2 þ q2; h1;k1; t2; 4Þ−Xp1;e1Kð3; e1;p1; a2;k2; t2ÞKðe1;p1; a02;k2 þ q2; 4Þ�: ðA83ÞEXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-21Since we consider only two-particle Green functions K that describe electron-hole pairs, the band index a2 has to be anelectron index in the first line and a hole index in the second line in Eq. (A83), respectively. We now evaluate Eq. (A82)using the RPA polarization function, introducing the exciton representation (A52) and assuming that Green functions arediagonal in the exciton basis, which results inεðkþ q;q; a; a0; t1;k2 þ q2;q2; a2; a02; t2Þ¼ n1δn1;n2δðt1 − t2Þδk;k2δq;q2δa;a2δa0;a02 − iℏXν1;ν2;Q1;Q2Gðν1;Q1; t1; ν1;Q1; t2ÞGðν2;Q2; t2; ν2;Q2; tþ1 Þ×Xk0� Xe;e0;h00;k00Va;e;e0;a0k;k0;k0−q;kþq×�Xh1;k1Φν1;Q1ðe0;k0 − q; h00;k00ÞΦ�ν1;Q1ðe000;k2; h1;k1ÞΦν2;Q2ðe0000;k2 þ q2; h1;k1ÞΦ�ν2;Q2ðe;k0; h00;k00Þδa2;e000δa02;e0000−Xe1;p1Φν1;Q1ðe0;k0 − q; h00;k00ÞΦ�ν1;Q1ðe1;p1; h000;k2ÞΦν2;Q2ðe1;p1; h0000;k2 þ q2ÞΦ�ν2;Q2ðe;k0; h00;k00Þδa2;h000δa02;h0000 −Xh;h0;e00;p00Va;h;h0;a0k;k0;k0−q;kþq×�Xh1;k1Φν1;Q1ðe00;p00; h0;k0 − qÞΦ�ν1;Q1ðe000;k2; h1;k1ÞΦν2;Q2ðe0000;k2 þ q2; h1;k1ÞΦ�ν2;Q2ðe00;p00; h;k0Þδa2;e000δa02;e0000−Xe1;p1Φν1;Q1ðe00;p00; h0;k0 − qÞΦ�ν1;Q1ðe1;p1; h000;k2ÞΦν2;Q2ðe1;p1; h0000;k2 þ q2ÞΦ�ν2;Q2ðe00;p00; h;k0Þδa2;h000δa02;h0000 �:ðA84ÞSimilar to the inverse dielectric matrix (A58), the dielectric matrix elements are identified asεðkþ q;q; a; a0; t1;k2 þ q2;q2; a2; a02; t2Þ ≈ εa;a02;a2;a0k;k2þq2;k2;kþqðt1; t2Þδq;q2δk;k2¼Xα;βðcaα;kÞ�ðca02β;k2þq2Þ�εαβ−qðt1; t2Þca2β;k2ca0α;kþqδq;q2δk;k2; ðA85Þassuming again diagonality in the k index. We can now expand the Coulomb matrix elements in terms of Wannierfunctions and make use of the completeness of Bloch states to derive the dielectric matrix in the Wannierrepresentation:εαβ−qðt1; t2Þ ¼Xa;a0;a2;a02caα;kca02β;kþqεa;a02;a2;a0k;kþq;k;kþqðt1; t2Þðca2β;kÞ�ðca0α;kþqÞ�¼ n1δn1;n2δðt1 − t2Þδα;β − iℏXν1;ν2;Q1;Q2Fðν1;Q1; t1; ν2;Q2; t2ÞXδVαδ−qΠ̃δβ−qðν1; ν2;Q1;Q2Þ ðA86ÞwithFðν1;Q1; t1; ν2;Q2; t2Þ ¼ Gðν1;Q1; t1; ν1;Q1; t2ÞGðν2;Q2; t2; ν2;Q2; tþ1 Þ: ðA87ÞUsing the momentum conservation implied by the exciton wave functions, the matrix elements of Π̃ can be explicitlycalculated asΠ̃δβ−qðν1; ν2;Q1;Q2Þ ¼ δQ2;Q1þqMδν1;ν2;Q1ð−qÞðMβν1;ν2;Q1ð−qÞÞ� ðA88ÞALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-22with the wave function overlapMδν1;ν2;Q1ðqÞ ¼Xk;e;hΦν1;Q1ðe; h;kÞ�Xe0ce0δ;Q1−k−qðceδ;Q1−kÞ�Φν2;Q1−qðe0; h;kÞ −Xh0ch0δ;k−qðchδ;kÞ�Φν2;Q1−qðe; h0;k − qÞ��:ðA89ÞFinally, we calculate the retarded dielectric matrix by transforming the Keldysh dielectric matrix (A86) to physicaltimes as discussed in the previous chapter. Assuming stationary exciton distribution functions, one can show thatiℏFretðν1;Q1; ν2;Q2;ωÞ ¼NXν2;Q2− NXν1;Q1Eν2;Q2− Eν1;Q1þ ℏωþ iγ: ðA90ÞWith this, we arrive at the final result [see Eq. (A75)]:εret;αβq ðωÞ ¼ δα;β −XδVαδq Πret;δβq ðωÞ ðA91Þwith the polarization matrixΠret;δβq ðωÞ ¼Xν1;ν2;Q1NXν2;Q1−q − NXν1;Q1Eν2;Q1−q − Eν1;Q1þ ℏωþ iγMδν1;ν2;Q1ðqÞðMβν1;ν2;Q1ðqÞÞ�: ðA92ÞWe, thus, derive a microscopic dielectric function thatdescribes excitonic screening in RPA corresponding to thebubble-type polaripzation shown in Fig. 5. Our result is ageneralization of the bound-state dielectric functionderived by Röpke and Der [76], which has also beenused in Ref. [8]. In fact, the Röpke-Der dielectric functioncan be considered a macroscopic limit of our result. Weemphasize that a microscopic treatment of screeningencoded in the matrix form of εret;αβq ðωÞ is essential tocapture local-field effects arising due to the layeredstructure of the crystal unit cell. The interplay of interlayerand intralayer interactions in the TMD bilayer is, thereby,naturally taken into account. For numerical calculations,we use a phenomenological damping γ ¼ min(10 meV,ℏω) in Eq. (A92) to ensure the correct analytic behavior inthe static limit ω → 0.6. Numerical details for the exciton Dyson equationThe exciton energy renormalizations (2) are evaluatedself-consistently with the self-energy (3) based on excitoneigenstates from the Bethe-Salpeter equation (BSE) (1). Weuse a Brillouin zone sampling with 48 × 48 × 1 k points,limiting the Brillouin zone to the regions with radius 2.3around the K and K0 points. The two highest valence andtwo lowest conduction bands are considered. For every totalexciton momentum Q, 48 exciton eigenstates from the BSEare taken into account. Thus, for the given sampling of theBrillouin zone, an energy range between the lowest boundstate at about −150 meV and about 100 meV above thecontinuum edge is covered. This is particularly important toconverge the excitonic polarization function that determinesthe dynamical excitonic screening. Even at elevated densitiesup to 1012 cm−2, only a small fraction of these states isoccupied. The larger fraction of unoccupied states is thenneeded to properly describe the excitonic polarization. Thefrequency integration in the Montroll-Ward self-energy(A71) extends from −500 to 500 meV using a samplingwith 80ω points. We have checked that exciton energyrenormalizations are converged to within 1 meV.FIG. 5. Feynman diagram for the excitonic RPA polarizationΠδβq ðωÞ inserted between two Coulomb interaction lines inWannier representation. The interaction vertices depicted in graycorrespond to the wave function overlaps defined in Eq. (A89).EXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-237. Numerical results for spin-triplet excitonsTo complement the numerical results for spin-singletexcitons shown in Fig. 2, we provide data for spin-tripletexcitons in Fig. 6. The larger population of triplet statesleads to an increase of fermionic and bosonic exchangeeffects among exciton triplets. As a result, the interactionbetween exciton triplets is slightly less repulsive thanbetween singlets.APPENDIX B: EXPERIMENTAL DETAILS1. Estimation of electron-hole pair densityThe injected electron-hole pair density is estimated usingthe following equation:neh ¼Pαfrepπr2EphðB1Þwith P the laser power (ranging from 250 nW to 25 μW),α ¼ 0.11 absorption, frep ¼ 80 MHz laser repetition rate,Eph ¼ 1.63 eV photon energy (resonant with the MoSe2A∶1s state), and r ¼ ð0.7� 0.1Þ μm radius of the focusedlaser spot (according to our previous work [20]). Theabsorption value is estimated from transfer matrix analysisof reflectance measurements on the heterobilayer; seeFig. 7. For estimating the peak density, the radius is chosenas r ¼ FWHM=ð2 ffiffiffiffiffiffiffiln 2p Þ so that the effective density isequal to the maximum density of the spot center,neh × πr2 ¼ N0. It is assumed that the majority of thecreated excitons form interlayer excitons due to the ultrafastcharge transfer.2. Time-dependent energy shifts at higher temperatureA streak camera image acquired at 70 K is presented inFig. 8(a). It is taken in the spectral region of neutralinterlayer triplets, demonstrating the time-dependentchange of the emission peak. As the time after the(b) 5K100K200K300KX energy shift (meV)−202e-h-pair density (1011 cm-2)0 5 10 15e-h-pair density (1011 cm-2)dip-dip+ ferm X + PSF+  bos X+ screen(a)X energy shift (meV)01020300 5 10 15FIG. 6. Exciton energy renormalizations induced by exciton-exciton interaction. (a) Cumulative density-dependent renormalization ofthe zero-momentum “gray” 1s-exciton (spin-triplet) energy at a temperature T ¼ 100 K, subsequently adding dipole-dipole interaction(dip-dip), fermionic exchange interaction and phase-space filling (þ ferm Xþ PSF), bosonic exchange interaction (þ bos X), andscreened bosonic exchange (þ screen). The latter represents the result of the full calculation. (b) Calculated temperature and densitydependence of energy renormalization for the zero-momentum 1s-exciton triplet. The result for T ¼ 5 K is obtained from extrapolatingthe high-temperature data.MeasFitFIG. 7. (Upper) Reflectance contrast derivative of MoSe2=WSe2 heterobilayer and respective fit from transfer matrix analysis.The extracted parametrized dielectric function is used to obtain theeffective absorption spectrum (lower). Spectrum of the pump laser(λ ¼ 761 nm) in red, overlapping with the absorption peak of theMoSe2 A∶1s state. Adapted from Ref. [20].ALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-24excitation increases, the exciton density and the strength ofthe PL signal decreases, as illustrated by the transient. Forthe analysis, we estimate the initial exciton density from thepump fluence and absorption and set the changes of thisdensity proportional to the decrease of the PL signal.Spectra taken at different times are shown in Fig. 8(b)including Gaussian fits to extract peak energies. The datademonstrates the time- and thus density-dependent shift ofthe emission to lower energies with a weak reversalobservable at later times toward higher energies.[1] H. Fang, C. Battaglia, C. Carraro, S. Nemsak, B. Ozdol, J. S.Kang, H. a. Bechtel, S. B. Desai, F. Kronast, A. a. Unal, G.Conti, C. Conlon, G. K. Palsson, M. C. Martin, A. M.Minor, C. S. Fadley, E. Yablonovitch, R. Maboudian, andA. Javey, Strong interlayer coupling in van der Waalsheterostructures built from single-layer chalcogenides,Proc. Natl. Acad. Sci. U.S.A. 111, 6198 (2014).[2] P. Rivera, J. R. Schaibley, A. M. Jones, J. S. Ross, S. Wu, G.Aivazian, P. Klement, K. Seyler, G. Clark, N. J. Ghimire, J.Yan, D. G. Mandrus, W. Yao, and X. Xu, Observation oflong-lived interlayer excitons in monolayer MoSe2-WSe2heterostructures, Nat. Commun. 6, 6242 (2015).[3] P. Rivera, H. Yu, K. L. Seyler, N. P. Wilson, W. Yao, and X.Xu, Interlayer valley excitons in heterobilayers of transitionmetal dichalcogenides, Nat. Nanotechnol. 13, 1004 (2018).[4] P. Merkl, F. Mooshammer, P. Steinleitner, A. Girnghuber,K.-Q. Lin, P. Nagler, J. Holler, C. Schüller, J. M. Lupton, T.Korn, S. Ovesen, S. Brem, E. Malic, and R. Huber,Ultrafasttransition between exciton phases in van der Waals hetero-structures, Nat. Mater. 18, 691 (2019).[5] S. Ovesen, S. Brem, C. Linderälv, M. Kuisma, T. Korn, P.Erhart,M. Selig, and E.Malic, Interlayer exciton dynamics invan derWaals heterostructures, Commun. Phys. 2, 1 (2019).[6] N. Peimyoo, T. Deilmann, F. Withers, J. Escolar, D. Nutting,T. Taniguchi, K. Watanabe, A. Taghizadeh, M. F. Craciun,K. S. Thygesen, and S. Russo, Electrical tuning of opticallyactive interlayer excitons in bilayer MoS2, Nat. Nano-technol. 16, 888 (2021).[7] E. Barré, O. Karni, E. Liu, A. L. O’Beirne, X. Chen,H. B. Ribeiro, L. Yu, B. Kim, K. Watanabe, T. Taniguchi,K. Barmak, C. H. Lui, S. Refaely-Abramson, F. H. daJornada, and T. F. Heinz, Optical absorption of interlayerexcitons in transition-metal dichalcogenide heterostruc-tures, Science 376, 406 (2022).[8] A. Steinhoff, M. Florian, M. Rösner, G. Schönhoff, T. O.Wehling, and F. Jahnke, Exciton fission in monolayertransition metal dichalcogenide semiconductors, Nat.Commun. 8, 1166 (2017).[9] J. Wang, J. Ardelean, Y. Bai, A. Steinhoff, M. Florian, F.Jahnke, X. Xu, M. Kira, J. Hone, and X.-Y. Zhu, Opticalgeneration of high carrier densities in 2D semiconductorheterobilayers, Sci. Adv. 5, eaax0145 (2019).[10] D. J. Morrow and X. Ma, Trapping interlayer excitons invan der Waals heterostructures by potential arrays, Phys.Rev. B 104, 195302 (2021).[11] W. Li, X. Lu, S. Dubey, L. Devenica, and A. Srivastava,Dipolar interactions between localized interlayer excitons invan der Waals heterostructures, Nat. Mater. 19, 624 (2020).[12] E. C. Regan, D. Wang, E. Y. Paik, Y. Zeng, L. Zhang, J. Zhu,A. H. MacDonald, H. Deng, and F. Wang, Emerging excitonphysics in transition metal dichalcogenide heterobilayers,Nat. Rev. Mater. 7, 778 (2022).[13] H. Guo, X. Zhang, and G. Lu, Tuning moiré excitons inJanus heterobilayers for high-temperature Bose-Einsteincondensation, Sci. Adv. 8, abp9757 (2022).(a) (b)FIG. 8. Time-dependent energy shift of neutral interlayer triplets at 70 K. (a) Streak camera image of the spectrally and time-resolved PLat a peak density of 1012 cm−2. The spectra are normalized at each time step, and the color scale is logarithmic. Corresponding PL transientis presented in the right. (b) Extracted PL spectra at different times after the excitation following the decrease of the PL signal proportionalto the exciton density. Included are Gaussian fits, and the dashed line indicates the peak energy immediately after the excitation.EXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-25https://doi.org/10.1073/pnas.1405435111https://doi.org/10.1038/ncomms7242https://doi.org/10.1038/s41565-018-0193-0https://doi.org/10.1038/s41563-019-0337-0https://doi.org/10.1038/s42005-018-0096-2https://doi.org/10.1038/s41565-021-00916-1https://doi.org/10.1038/s41565-021-00916-1https://doi.org/10.1126/science.abm8511https://doi.org/10.1038/s41467-017-01298-6https://doi.org/10.1038/s41467-017-01298-6https://doi.org/10.1126/sciadv.aax0145https://doi.org/10.1103/PhysRevB.104.195302https://doi.org/10.1103/PhysRevB.104.195302https://doi.org/10.1038/s41563-020-0661-4https://doi.org/10.1038/s41578-022-00440-1https://doi.org/10.1126/sciadv.abp9757[14] X. Sun, Y. Zhu, H. Qin, B. Liu, Y. Tang, T. Lü, S. Rahman,T. Yildirim, and Y. Lu, Enhanced interactions of interlayerexcitons in free-standing heterobilayers, Nature (London)610, 478 (2022).[15] J. Wang, Q. Shi, E. M. Shih, L. Zhou, W. Wu, Y. Bai, D.Rhodes, K. Barmak, J. Hone, C. R. Dean, and X. Y. Zhu,Diffusivity reveals three distinct phases of interlayer ex-citons inMoSe2=WSe2 heterobilayers, Phys. Rev. Lett. 126,106804 (2021).[16] L. Yuan, B. Zheng, J. Kunstmann, T. Brumme, A. B.Kuc, C. Ma, S. Deng, D. Blach, A. Pan, and L. Huang,Twist-angle-dependent interlayer exciton diffusion inWS2-WSe2 heterobilayers, Nat. Mater. 19, 617 (2020).[17] Z. Sun, A. Ciarrocchi, F. Tagarelli, J. F. Gonzalez Marin, K.Watanabe, T. Taniguchi, and A. Kis, Excitonic transportdriven by repulsive dipolar interaction in a van der Waalsheterostructure, Nat. Photonics 16, 79 (2022).[18] D. Erkensten, S. Brem, R. Perea-Causín, and E. Malic,Microscopic origin of anomalous interlayer exciton trans-port in van der Waals heterostructures, Phys. Rev. Mater. 6,094006 (2022).[19] J. Choi, J. Embley, D. D. Blach, R. Perea-Causín, D.Erkensten, D. S. Kim, L. Yuan, W. Y. Yoon, T. Taniguchi,K. Watanabe, K. Ueno, E. Tutuc, S. Brem, E. Malic, X. Li,and L. Huang, Fermi pressure and Coulomb repulsiondriven rapid hot plasma expansion in a van der Waalsheterostructure, Nano Lett. 23, 4399 (2023).[20] E. Wietek, M. Florian, J. Göser, T. Taniguchi, K. Watanabe,A. Högele, M. M. Glazov, A. Steinhoff, and A. Chernikov,Nonlinear and negative effective diffusivity of interlayerexcitons in moiré-free heterobilayers, Phys. Rev. Lett. 132,016202 (2024).[21] B. Miller, A. Steinhoff, B. Pano, J. Klein, F. Jahnke, A.Holleitner, and U. Wurstbauer, Long-lived direct andindirect interlayer excitons in van der Waals heterostruc-tures, Nano Lett. 17, 5229 (2017).[22] P. Nagler, G. Plechinger, M. V. Ballottin, A. Mitioglu, S.Meier, N. Paradiso, C. Strunk, A. Chernikov, P. C. M.Christianen, C. Schüller, and T. Korn, Interlayer excitondynamics in a dichalcogenide monolayer heterostructure,2D Mater. 4, 025112 (2017).[23] W. Li, X. Lu, J. Wu, and A. Srivastava, Optical control ofthe valley Zeeman effect through many-exciton interactions,Nat. Nanotechnol. 16, 148 (2021).[24] L. V. Butov, A. Zrenner, G. Abstreiter, G. Böhm, and G.Weimann, Condensation of indirect excitons in coupledAlAs=GaAs quantum wells, Phys. Rev. Lett. 73, 304(1994).[25] L. V. Butov, C. W. Lai, A. L. Ivanov, A. C. Gossard, andD. S. Chemla, Towards Bose-Einstein condensation ofexcitons in potential traps, Nature (London) 417, 47 (2002).[26] A. L. Ivanov, Quantum diffusion of dipole-oriented indirectexcitons in coupled quantum wells, Europhys. Lett. 59, 586(2002).[27] C. Schindler and R. Zimmermann, Analysis of the exciton-exciton interaction in semiconductor quantum wells, Phys.Rev. B 78, 045313 (2008).[28] B. Laikhtman and R. Rapaport, Exciton correlations incoupled quantum wells and their luminescence blue shift,Phys. Rev. B 80, 195313 (2009).[29] S. V. Lobanov, N. A. Gippius, and L. V. Butov, Theory ofcondensation of indirect excitons in a trap, Phys. Rev. B 94,245401 (2016).[30] L. A. Jauregui, A. Y. Joe, K. Pistunova, D. S. Wild, A. A.High, Y. Zhou, G. Scuri, K. De Greve, A. Sushko, C.-H. Yu,T. Taniguchi, K. Watanabe, D. J. Needleman, M. D. Lukin,H. Park, and P. Kim, Electrical control of interlayer excitondynamics in atomically thin heterostructures, Science 366,870 (2019).[31] D. Unuchek, A. Ciarrocchi, A. Avsar, Z. Sun, K. Watanabe,T. Taniguchi, and A. Kis, Valley-polarized exciton currentsin a van der Waals heterostructure, Nat. Nanotechnol. 14,1104 (2019).[32] Y. E. Lozovik and O. L. Berman, Phase transitions in asystem of two coupled quantum wells, JETP Lett. 64, 573(1996).[33] L. V. Butov, A. A. Shashkin, V. T. Dolgopolov, K. L.Campman, and A. C. Gossard, Magneto-optics of thespatially separated electron and hole layers inGaAs=AlGaAs coupled quantum wells, Phys. Rev. B 60,8753 (1999).[34] V. May, F. Boldt, and K. Henneberger, Many-body theoryfor the dense exciton gas of direct semiconductors I.General considerations, Phys. Status Solidi (b) 129, 717(1985).[35] M. Combescot and O. Betbeder-Matibet, The effectivebosonic Hamiltonian for excitons reconsidered, Europhys.Lett. 58, 87 (2002).[36] X. Sun, Y. Zhu, H. Qin, B. Liu, Y. Tang, T. Lü, S. Rahman,T. Yildirim, and Y. Lu, Enhanced interactions of interlayerexcitons in free-standing heterobilayers, Nature (London)610, 478 (2022).[37] N. Götting, F. Lohof, and C. Gies, Moiré-Bose-Hubbardmodel for interlayer excitons in twisted transition metaldichalcogenide heterostructures, Phys. Rev. B 105, 165419(2022).[38] F. Lohof, J. Michl, A. Steinhoff, B. Han, M. v. Helversen,S. Tongay, K. Watanabe, T. Taniguchi, S. Höfling, S.Reitzenstein, C. Anton-Solanas, C. Gies, and C.Schneider, Confined-state physics and signs of fermioniza-tion of moiré excitons in WSe2=MoSe2 heterobilayers, 2DMater. 10, 034001 (2023).[39] M. Katzer, M. Selig, L. Sigl, M. Troue, J. Figueiredo, J.Kiemle, F. Sigger, U. Wurstbauer, A. W. Holleitner, and A.Knorr, Exciton-phonon scattering: Competition between thebosonic and fermionic nature of bound electron-hole pairs,Phys. Rev. B 108, L121102 (2023).[40] F. Boldt, K. Henneberger, and V. May, Many-body theoryfor the dense exciton gas of direct semiconductors II.Calculation of exciton level shift and damping in depend-ence on exciton density, Phys. Status Solidi (b) 130, 675(1985).[41] H. Haug and S. W. Koch, Quantum Theory of the Opticaland Electronic Properties of Semiconductors, 5th ed.(World Scientific, Singapore, 2009).[42] A. Steinhoff, M. Rösner, F. Jahnke, T. O. Wehling, and C.Gies, Influence of excited carriers on the optical andelectronic properties of MoS2, Nano Lett. 14, 3743 (2014).[43] A. Chernikov, C. Ruppert, H. M. Hill, A. F. Rigosi, and T. F.Heinz, Population inversion and giant bandgap renormal-ALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-26https://doi.org/10.1038/s41586-022-05193-zhttps://doi.org/10.1038/s41586-022-05193-zhttps://doi.org/10.1103/PhysRevLett.126.106804https://doi.org/10.1103/PhysRevLett.126.106804https://doi.org/10.1038/s41563-020-0670-3https://doi.org/10.1038/s41566-021-00908-6https://doi.org/10.1103/PhysRevMaterials.6.094006https://doi.org/10.1103/PhysRevMaterials.6.094006https://doi.org/10.1021/acs.nanolett.3c00678https://doi.org/10.1103/PhysRevLett.132.016202https://doi.org/10.1103/PhysRevLett.132.016202https://doi.org/10.1021/acs.nanolett.7b01304https://doi.org/10.1088/2053-1583/aa7352https://doi.org/10.1038/s41565-020-00804-0https://doi.org/10.1103/PhysRevLett.73.304https://doi.org/10.1103/PhysRevLett.73.304https://doi.org/10.1038/417047ahttps://doi.org/10.1209/epl/i2002-00144-3https://doi.org/10.1209/epl/i2002-00144-3https://doi.org/10.1103/PhysRevB.78.045313https://doi.org/10.1103/PhysRevB.78.045313https://doi.org/10.1103/PhysRevB.80.195313https://doi.org/10.1103/PhysRevB.94.245401https://doi.org/10.1103/PhysRevB.94.245401https://doi.org/10.1126/science.aaw4194https://doi.org/10.1126/science.aaw4194https://doi.org/10.1038/s41565-019-0559-yhttps://doi.org/10.1038/s41565-019-0559-yhttps://doi.org/10.1134/1.567264https://doi.org/10.1134/1.567264https://doi.org/10.1103/PhysRevB.60.8753https://doi.org/10.1103/PhysRevB.60.8753https://doi.org/10.1002/pssb.2221290232https://doi.org/10.1002/pssb.2221290232https://doi.org/10.1209/epl/i2002-00609-3https://doi.org/10.1209/epl/i2002-00609-3https://doi.org/10.1038/s41586-022-05193-zhttps://doi.org/10.1038/s41586-022-05193-zhttps://doi.org/10.1103/PhysRevB.105.165419https://doi.org/10.1103/PhysRevB.105.165419https://doi.org/10.1088/2053-1583/acd265https://doi.org/10.1088/2053-1583/acd265https://doi.org/10.1103/PhysRevB.108.L121102https://doi.org/10.1002/pssb.2221300231https://doi.org/10.1002/pssb.2221300231https://doi.org/10.1021/nl500595uization in atomically thinWS2 layers, Nat. Photonics 9, 466(2015).[44] S. Ulstrup, A. G. Čabo, J. A. Miwa, J. M. Riley, S. S.Grønborg, J. C. Johannsen, C. Cacho, O. Alexander, R. T.Chapman, E. Springate, M. Bianchi, M. Dendzik, J. V.Lauritsen, P. D. C. King, and P. Hofmann, Ultrafast bandstructure control of a two-dimensional heterostructure,ACS Nano 10, 6315 (2016).[45] S. Zhao, Z. Li, X. Huang, A. Rupp, J. Göser, I. A. Vovk,S. Y. Kruchinin, K. Watanabe, T. Taniguchi, I. Bilgin, A. S.Baimuratov, and A. Högele, Excitons in mesoscopicallyreconstructed moiré heterostructures, Nat. Nanotechnol.18, 572 (2023).[46] K. Wagner, J. Zipfel, R. Rosati, E. Wietek, J. D. Ziegler, S.Brem, R. Perea-Causín, T. Taniguchi, K. Watanabe, M. M.Glazov, E. Malic, and A. Chernikov, Nonclassical excitondiffusion in monolayer WSe2, Phys. Rev. Lett. 127, 076801(2021).[47] A. Segura, L. Artús, R. Cuscó, T. Taniguchi, G. Cassabois,and B. Gil, Natural optical anisotropy of h-BN: Highestgiant birefringence in a bulk crystal through the mid-infrared to ultraviolet range, Phys. Rev. Mater. 2,024001 (2018).[48] A. Arora, M. Koperski, K. Nogajewski, J. Marcus, C.Faugeras, and M. Potemski, Excitonic resonances in thinfilms ofWSe2: From monolayer to bulk material, Nanoscale7, 10421 (2015).[49] M. Selig, G. Berghäuser, A. Raja, P. Nagler, C. Schüller,T. F. Heinz, T. Korn, A. Chernikov, E. Malic, and A. Knorr,Excitonic linewidth and coherence lifetime in monolayertransition metal dichalcogenides, Nat. Commun. 7, 13279(2016).[50] F. Cadiz, E. Courtade, C. Robert, G. Wang, Y. Shen, H. Cai,T. Taniguchi, K. Watanabe, H. Carrere, D. Lagarde, M.Manca, T. Amand, P. Renucci, S. Tongay, X. Marie, and B.Urbaszek, Excitonic linewidth approaching the homo-geneous limit in MoS2-based van der Waals heterostruc-tures, Phys. Rev. X 7, 021026 (2017).[51] G. Gupta and K. Majumdar, Fundamental exciton linewidthbroadening in monolayer transition metal dichalcogenides,Phys. Rev. B 99, 085412 (2019).[52] Y.-h. Chan, J. B. Haber, M. H. Naik, J. B. Neaton, D. Y. Qiu,F. H. da Jornada, and S. G. Louie, Exciton lifetime andoptical line width profile via exciton-phonon interactions:Theory and first-principles calculations for monolayerMoS2, Nano Lett. 23, 3971 (2023).[53] N. Kumar, Q. Cui, F. Ceballos, D. He, Y. Wang, and H.Zhao, Exciton-exciton annihilation in MoSe2 monolayers,Phys. Rev. B 89, 125427 (2014).[54] S. Mouri, Y. Miyauchi, M. Toh, W. Zhao, G. Eda, andK. Matsuda, Nonlinear photoluminescence in atomicallythin layered WSe2 arising from diffusion-assistedexciton-exciton annihilation, Phys. Rev. B 90, 155449(2014).[55] L. Yuan, T. Wang, T. Zhu, M. Zhou, and L. Huang, Excitondynamics, transport, and annihilation in atomically thintwo-dimensional semiconductors, J. Phys. Chem. Lett. 8,3371 (2017).[56] Y. Hoshi, T. Kuroda, M. Okada, R. Moriya, S. Masubuchi,K. Watanabe, T. Taniguchi, R. Kitaura, and T. Machida,Suppression of exciton-exciton annihilation in tungstendisulfide monolayers encapsulated by hexagonal boronnitrides, Phys. Rev. B 95, 241403(R) (2017).[57] A. Raja, L. Waldecker, J. Zipfel, Y. Cho, S. Brem, J. D.Ziegler, M. Kulig, T. Taniguchi, K. Watanabe, E. Malic, T. F.Heinz, T. C. Berkelbach, and A. Chernikov, Dielectricdisorder in two-dimensional materials, Nat. Nanotechnol.14, 832 (2019).[58] L. Sigl, F. Sigger, F. Kronowetter, J. Kiemle, J. Klein, K.Watanabe, T. Taniguchi, J. J. Finley, U. Wurstbauer, andA.W. Holleitner, Signatures of a degenerate many-bodystate of interlayer excitons in a van der Waals heterostack,Phys. Rev. Res. 2, 042044(R) (2020).[59] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevX.14.031025 for additionaldata on observations demonstrating the excitation densitydependent energy shift of interlayer excitons.[60] P. Giannozzi et al., QUANTUM ESPRESSO: A modular andopen-source software project for quantum simulations ofmaterials, J. Phys. Condens. Matter 21, 395502 (2009).[61] P. Giannozzi et al., Advanced capabilities for materialsmodelling with QUANTUM ESPRESSO, J. Phys. Condens.Matter 29, 465901 (2017).[62] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalizedgradient approximation made simple, Phys. Rev. Lett. 78,1396(E) (1997).[63] M. J. van Setten, M. Giantomassi, E. Bousquet, M. J.Verstraete, D. R. Hamann, X. Gonze, and G. M.Rignanese, The PseudoDojo: Training and grading a 85element optimized norm-conserving pseudopotential table,Comput. Phys. Commun. 226, 39 (2018).[64] R. Gillen and J. Maultzsch, Interlayer excitons inMoSe2=WSe2 heterostructures from first principles, Phys.Rev. B 97, 165306 (2018).[65] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistentand accurate ab initio parametrization of density functionaldispersion correction (DFT-D) for the 94 elements H-Pu,J. Chem. Phys. 132, 154104 (2010).[66] K. Nakamura, Y. Yoshimoto, Y. Nomura, T. Tadano, M.Kawamura, T. Kosugi, K. Yoshimi, T. Misawa, and Y.Motoyama, RESPACK: An ab initio tool for derivation ofeffective low-energy model of material, Comput. Phys.Commun. 261, 107781 (2021).[67] G. Schönhoff, M. Rösner, R. E. Groenewald, S. Haas, andT. O. Wehling, Interplay of screening and superconductivityin low-dimensional materials, Phys. Rev. B 94, 134504(2016).[68] C. Steinke, T. O. Wehling, and M. Rösner, Coulomb-engineered heterojunctions and dynamical screening intransition metal dichalcogenide monolayers, Phys. Rev.B 102, 115111 (2020).[69] M. Rösner, E. Şaşioğlu, C. Friedrich, S. Blügel, and T. O.Wehling, Wannier function approach to realistic Coulombinteractions in layered materials and heterostructures,Phys. Rev. B 92, 085102 (2015).[70] M. Florian, M. Hartmann, A. Steinhoff, J. Klein, A. W.Holleitner, J. J. Finley, T. O. Wehling, M. Kaniber, and C.Gies, The dielectric impact of layer distances on exciton andtrion binding energies in van der Waals heterostructures,Nano Lett. 18, 2725 (2018).EXCITON-EXCITON INTERACTIONS IN VAN DER WAALS … PHYS. REV. X 14, 031025 (2024)031025-27https://doi.org/10.1038/nphoton.2015.104https://doi.org/10.1038/nphoton.2015.104https://doi.org/10.1021/acsnano.6b02622https://doi.org/10.1038/s41565-023-01356-9https://doi.org/10.1038/s41565-023-01356-9https://doi.org/10.1103/PhysRevLett.127.076801https://doi.org/10.1103/PhysRevLett.127.076801https://doi.org/10.1103/PhysRevMaterials.2.024001https://doi.org/10.1103/PhysRevMaterials.2.024001https://doi.org/10.1039/C5NR01536Ghttps://doi.org/10.1039/C5NR01536Ghttps://doi.org/10.1038/ncomms13279https://doi.org/10.1038/ncomms13279https://doi.org/10.1103/PhysRevX.7.021026https://doi.org/10.1103/PhysRevB.99.085412https://doi.org/10.1021/acs.nanolett.3c00732https://doi.org/10.1103/PhysRevB.89.125427https://doi.org/10.1103/PhysRevB.90.155449https://doi.org/10.1103/PhysRevB.90.155449https://doi.org/10.1021/acs.jpclett.7b00885https://doi.org/10.1021/acs.jpclett.7b00885https://doi.org/10.1103/PhysRevB.95.241403https://doi.org/10.1038/s41565-019-0520-0https://doi.org/10.1038/s41565-019-0520-0https://doi.org/10.1103/PhysRevResearch.2.042044http://link.aps.org/supplemental/10.1103/PhysRevX.14.031025http://link.aps.org/supplemental/10.1103/PhysRevX.14.031025http://link.aps.org/supplemental/10.1103/PhysRevX.14.031025http://link.aps.org/supplemental/10.1103/PhysRevX.14.031025http://link.aps.org/supplemental/10.1103/PhysRevX.14.031025http://link.aps.org/supplemental/10.1103/PhysRevX.14.031025http://link.aps.org/supplemental/10.1103/PhysRevX.14.031025https://doi.org/10.1088/0953-8984/21/39/395502https://doi.org/10.1088/1361-648X/aa8f79https://doi.org/10.1088/1361-648X/aa8f79https://doi.org/10.1103/PhysRevLett.78.1396https://doi.org/10.1103/PhysRevLett.78.1396https://doi.org/10.1016/j.cpc.2018.01.012https://doi.org/10.1103/PhysRevB.97.165306https://doi.org/10.1103/PhysRevB.97.165306https://doi.org/10.1063/1.3382344https://doi.org/10.1016/j.cpc.2020.107781https://doi.org/10.1016/j.cpc.2020.107781https://doi.org/10.1103/PhysRevB.94.134504https://doi.org/10.1103/PhysRevB.94.134504https://doi.org/10.1103/PhysRevB.102.115111https://doi.org/10.1103/PhysRevB.102.115111https://doi.org/10.1103/PhysRevB.92.085102https://doi.org/10.1021/acs.nanolett.8b00840[71] Z. Gong, G.-B. Liu, H. Yu, D. Xiao, X. Cui, X. Xu, and W.Yao, Magnetoelectric effects and valley-controlled spinquantum gates in transition metal dichalcogenide bilayers,Nat. Commun. 4, 2053 (2013).[72] F. Katsch, M. Selig, A. Carmele, and A. Knorr, Theory ofexciton-exciton interactions in monolayer transition metaldichalcogenides, Phys. Status Solidi (b) 255, 1800185 (2018).[73] D. Kremp, M. Schlanges, and W.-D. Kraeft, QuantumStatistics of Nonideal Plasmas (Springer, Berlin, 2005).[74] A. Schobert, J. Berges, E. G. C. P. van Loon, M. A. Sentef,S. Brener, M. Rossi, and T. O. Wehling, Ab initioelectron-lattice downfolding: Potential energy landscapes,anharmonicity, and molecular dynamics in charge densitywave materials, SciPost Phys. 16, 046 (2024).[75] D. Erben, A. Steinhoff, C. Gies, G. Schönhoff, T. O.Wehling,and F. Jahnke, Excitation-induced transition to indirect bandgaps in atomically thin transition-metal dichalcogenidesemiconductors, Phys. Rev. B 98, 035434 (2018).[76] G. Röpke and R. Der, The influence of two-particle states(excitons) on the dielectric function of the electron-holeplasma, Phys. Status Solidi (b) 92, 501 (1979).ALEXANDER STEINHOFF et al. PHYS. REV. X 14, 031025 (2024)031025-28https://doi.org/10.1038/ncomms3053https://doi.org/10.1002/pssb.201800185https://doi.org/10.21468/SciPostPhys.16.2.046https://doi.org/10.1103/PhysRevB.98.035434https://doi.org/10.1002/pssb.2220920220