Dipolar Coupling at Interfaces of Ultrathin Semiconductors, Semimetals, Plasmonic Nanoparticles, and Molecules

Recent progress in growth techniques has enabled the fabrication of stacks of transition metal dichalcogenide monolayers combined with different nanostructures ranging from other 2D layers over dye molecules to even plasmonic nanoparticles. Such structures promise to combine the optoelectric properties of the constituents allowing to design structures with desired properties. For all of these examples, a detailed knowledge of the coupling among different constituents is crucial. In this article, a unified description is presented based on Maxwell Bloch equations to describe dipolar interactions among different types of heterostructures. Exemplary, Förster‐type energy transfer from dye molecules to MoS2 monolayers, strong coupling at MoSe2–metal nanoparticle interfaces, Meitner–Auger‐like interlayer coupling in WSe2–graphene stacks, and relaxation processes of hot interlayer excitons in MoSe2–WSe2 heterobilayers are discussed.


Introduction
Monolayers of transition metal dichalcogenides (TMDCs) exhibit an extremely strong light-matter interaction [1] with absorption rates up to 10% in the visible spectrum. [2,3]This fact is even more remarkable considering their 2D nature with less than 1 nm thickness.6][7][8][9][10][11][12][13] Moreover, the atomical thickness provides a strong sensitivity to the environment: atomically thin materials can be manipulated by the choice of the embedding material and its geometry [14] as well as by defects [15] or functionalization. [16]When coated in leaky cavities, 2D TMDCs exhibit an interesting interplay of phonon and photon dissipation, leading to anomalous dispersion and negative mass effects. [17]Due to the variety of feasible material combinations, functionalization proves to be a powerful tool for tailoring electronic and optical properties on the nanoscale.The scope of this work is to present a theoretical framework to analyze hybrids of a TMDC layer parallel to the xy-plane functionalized with other nanomaterials as displayed in Figure 1.We focus on interactions that are mediated by the electric near-field, starting from the setup of the general light-matter Hamiltonian valid for all discussed hybrids in Section 2. The corresponding matrix elements are specified for molecules, metals, graphene, and TMDCs in Section 2.1-2.3,where we also develop the Heisenberg equations of motion for excitations in these TMDCs as the overall substrate.Subsequently, we apply the formalism on hybrids of a TMDC monolayer functionalized with molecules, metal nanoparticles (MNP) or graphene, and a TMDC heterobilayer.In this work, we focus on coupling effects that are mediated by the mutual electric fields between the constituents: [29][30][31][32] 2) MNPs provide an impressive amplification of the electric near-field, [33,34] Section 4. The interaction of TMDC excitons with the MNP plasmons features exciton localization, Section 4.1, and allows to enter the strong coupling regime, [35][36][37][38][39][40][41][42][43][44] Section 4.2.

Heterostructure Hamiltonian
The general heterostructure Hamiltonian considered throughout this work reads with annihilation (creation) operators i ð †Þ for electrons with compound quantum number i, which will be specified in the next subsections to distinguish the different applications (i-iv).The first term in Equation (1) accounts for the dispersion of the electrons of all constituents of the hybrid.The second term incorporates monopole contributions of electron-electron interactions V ijj 0 i 0 ¼ Z d 3 rΨ Ã i ðrÞΨ Ã j ðr 0 ÞVðr, r 0 ÞΨ j 0 ðr 0 ÞΨ i 0 ðrÞ (2)   with Coulomb potential Vðr, r 0 Þ and electron wave functions Ψ i ðrÞ.The third term accounts for the electric field-matter coupling in dipole approximation with the Rabi energy including the charge q of the carriers and the electric field Eðr, tÞ.
In course of evaluating the Rabi energy for the different materials, we approximate the electric field to be constant within one unit cell.This allows to separate electronic-and electric-field coordinates and to identify electronic intra-and inter-band transitions.We specify the matrix elements V ijj 0 i 0 and Ω ii 0 , for the molecule, metal, and semiconductor/semimetal electrons in Section 2.1-2.3.The framework allows to derive diverse energy-transfer mechanisms between the different material systems on the same footing by considering the respective contributions in the Hamiltonian when deriving Heisenberg's equation of motion as illustrated in Sections 3-6.For each hybrid, the electric field is determined by Maxwell's equations, which can be accessed in the Green's function formalism [16,62] with the 2D polarization densities P l Q k for layer l and the incident external electric field E 0 Q k as a source.To achieve this algebraic form for the electric field, we assumed the individual layers to act as infinitesimally thin. [63]The Green's dyadic G Q k ðz, z l , ωÞ transfers the polarization densities located at z l to the observers' position z.It depends on the scalar Green's function The Green's functions for the chosen geometries (additional dielectric interfaces in z-direction) are given in the appendix.Each considered heterostructure, see Figure 1, consists of two layers ðl ¼ 0, 1Þ at z ¼ z 0 , z 1 , respectively.To solve the electric near-field, that is responsible for the coupling between the layers, we apply the quasi-static approximation, setting ω !0 in the Green's dyadic, Equation (5).This corresponds to neglecting propagation effects and is valid as long as the layer distance Δz ¼ jz 0 À z 1 j ( λ is small compared to the wavelength of the exciting light field E 0 Q k . [64]In contrast, we can simplify Q k !0 to derive the far-field of the heterostructure for perpendicular plane-wave excitation.

Molecules
In the case of molecules, the compound quantum number i in Equation (1) accounts for the molecular orbitals λ ¼ H, L (highest occupied molecular orbital (HOMO), lowest unoccupied molecular orbital (LUMO)).The single-particle energies are already treated as effective excitonic energies.The corresponding Hamiltonian reads The first term represents the energy ε λ of the molecular orbitals, [65] described as fermions, a sufficient approach in linear optics. [66]The molecular transition energy reads For molecules, the electron-electron interactions in the general light-matter Hamiltonian, Equation (1), are already included in the effective two-level energies ε λ .The second term accounts for the light-matter coupling where we have assumed a point-like molecule, located at the center of charge r ¼ r 1 , exhibiting a dipole moment d λλ .For the plot, we assume exemplary a dichloromethane (DCM) molecule.The strength of the molecular dipole moment, d λλ ¼ d HL is taken from ab initio calculations. [67]

Metals
Metals have a partially occupied conduction band that allows for intraband as well as interband transitions.Depending on the spectral range of interest, the optical response in metal structures is either dominated by intraband transitions, known as plasmons, or interband transitions. [68]Accordingly, the compound quantum number i in the Hamiltonian in Equation (1) contains the band index λ, constituting for valence ðvÞ and conduction ðcÞ band, and the electronic momentum k.The Hamiltonian is given by The band structure is accounted for by the first term, which includes the electronic dispersion ε λk , renormalized due to electron-electron interactions.The second term covers interband transitions between different bands λ 6 ¼ λ, with the transition probability determined by the optical dipole element d λλk .The last term describes intraband transitions where momentum Q is transferred from the electric field to an electron within one band.The probability of this process is given by the Fourier transform of the electric field E Q and the electron charge q ¼ Àe normalized by the volume V.Only the conduction band (λ ¼ c) needs to be considered for intraband transitions as the valence band is fully occupied.
The permittivity of gold, ϵ Au , can be derived using the Heisenberg equation of motion framework.The intraband contribution results in the Drude susceptibility for a quasi-free electron plasma in the gold conduction band, [69,70] while the interband transitions can, in its most simple form, be approximated by Lorentz-shaped contributions around the respective transition energies.In gold, interband transitions become relevant above 2.4 eV [68] and while the correct allocation of the visible interband transitions in the gold band structure is still controversial, there are several models that fit the experimental data in refs.75][76] In this paper, we use the permittivity provided by ref. [74].
The first two terms resemble a standard Drude model that incorporates dielectric screening to account for strongly polarized d-orbitals.Damping rates Γ, Γ a are added phenomenologically to fit experimental curves measured from bulk material.All parameters are given in the appendix.

Plasmonic Crystals
The optical response of MNPs is given in Mie theory similar to the electric field of a point dipole. [77,78]In the case of a 2D plasmonic crystal (PC) consisting of periodically arranged MNPs, the particle and lattice geometry is as important for the optical properties as the material.The macroscopic polarization density is given by [79,80] with an effective PC polarizability α Ã ðωÞ, taking MNP interactions and Umklapp processes into account. [79]For noninteracting MNPs, Equation ( 9) is also applicable with MNP polarizability αðωÞ in Mie-Gans theory [78,81] that depends on the surrounding and MNP permittivity ϵ Au ðωÞ, Equation (8). [74]We provide the MNP polarizability in the appendix.

2D Materials
In the case of 2D semiconductors and semimetals, for example, TMDCs and graphene, the compound quantum numbers of the electrons constitute the layer l, the band λ ¼ c, v (conduction band, valence band) and the momentum ξ þ k where ξ denotes the considered high symmetry points, in our case K=K 0 with the 2D wave number k therein, that is, i ¼ ðl, λ, ξ, kÞ.The semiconductor/semimetal Hamiltonian takes the form points; their electronic dispersion is approximated parabolic.For simplicity, the TMDC is illustrated for one spin direction.
The first term accounts for the energy of electrons with the dispersion ε lξ λk .The second term describes monopole, Coulomb electron-electron interactions. [82]As we are only interested in energy-transfer mechanisms, Equation (10) disregards electronic tunnel processes between the constituents.The third term accounts for light-matter interaction.It is treated in thin film [63] and low wave number approximation, such that the momentum dependence of the dipole moment d lξ λλ drops. [83]urther, E Q k ðz l , tÞ ¼ ∫ d 2 re ÀiQ k ⋅r k Eðr k , z l , tÞ accounts for the Fourier component of the electric field inside the semiconductor/semimetal layer z ¼ z l with respect to the in-plane component.Equation ( 10) is applicable for mono-and multilayer.We concentrate on the description of monolayer (MoSe 2 , MoS 2 ) and heterobilayer (WSe 2 /graphene, MoSe 2 /WSe 2 ), where the layer number l also specifies the layer material.

Graphene
In the case of graphene, we concentrate on the reciprocal space near the K points where the valence and conduction bands touch and form the Dirac cones.Since graphene does not exhibit a bandgap, excitonic effects are typically of minor importance in the optical range and we use the free particle band structure.As displayed in Figure 2, the graphene dispersion is linear close to the K point: ε k ¼ AEℏv F jkj for the conduction (þ) or valence ðÀÞ bands, respectively. [51]

TMDCs as Substrates
In the monolayer limit, TMDCs facilitate an optically accessible direct bandgap at the K=K 0 points.In their vicinity, we approximate the dispersion ε ξ λk to be parabolic, parameterized from DFT calculations, [84] see Figure 2. Furthermore, TMDCs feature a spin-splitting of the electronic valence and conduction bands.Their dispersion ε ξ λk depends on the valley but also on the spin index s which we absorb in a compound quantum number ξ for valley and spin.The optical dipole moment d ξ cv of the TMDC inter-band transition is taken from ab initio electronic structure calculations. [85]For the hybrid materials considered in this manuscript, see Figure 1, a TMDC monolayer serves as the substrate.Therefore, we derive already at this point the Heisenberg equation of motion for the excitonic transition [86] which all hybrids have in common given by the expectation value of the exciton operator denotes the excitonic wave function in real and momentum space, respectively.α ξ and β ξ are the ratios of the effective masses for the respective excitonic configuration ξ as defined in ref. [86].
The relative motion q k , arising from the electron-electron Coulomb interactions in Equation ( 1), is responsible for the formation of bound exciton states with index ν.The remaining degree of freedom is the excitonic center-of-mass motion Q k .Using the center-of-mass momentum Q k instead of the full set of electron/hole momenta, see the k-axis in Figure 2, corresponds to the exciton picture.As we perform the calculations in the exciton picture, we also give schematic dispersions in further sections in the Q k -space.For the derivation of the Heisenberg equation, we restrict ourselves to low exciton densities, that is, to a regime where Pauli blocking and exciton-exciton scattering is negligible and excitons can be treated as bosonic particles.Screening would occur on the same order as the neglected fermionic contributions for the exciton operators pξν Q k , [86,87] which modify the linear equations of motion.[90] The main effect of screening would be a reduction of the interaction strength between and in the layers.We arrive at the excitonic Bloch equation [82,86,91] ðℏω À E ξν A phenomenological damping constant γ ξν mainly stemming from electron-phonon scattering is introduced as a dephasing. [92]The excitonic dispersion reads in parabolic approximation with exciton transition energy E ξν and M ξ as the exciton mass at valley ξ.The spin splitting of the electronic bands results in two types of excitons with different transition energies corresponding to their spins, namely the A and the B exciton.The macroscopic polarization within the TMDC layer l in the excitonic picture reads The macroscopic polarization density connects the exciton dynamics to the generated electric field E Q k via Equation (4), that mediates the interaction between the constituents of the hybrids and that has to be solved self-consistently as illustrated in the following sections.
The excitonic transition p ξν Q k is the crucial microscopic quantity to describe linear optical experiments as in Sections 3 and 4, for example, direct transmission, reflection, and absorption.However, if it comes to photoluminescence or nonlinear optical spectroscopy, one further needs to consider the excitonic occupation N ξν Q k , that is defined as the correlated (c) expectation value of an electron-hole pair [93,94] N ξν From the Heisenberg equation of motion, it is possible to derive the Boltzmann equation for the excitonic occupations. [28]ere, we provide a general form of the Boltzmann equation for low exciton densities N ξν Q k ( 1, allowing to neglect contributions that are nonlinear in N ξν Equation ( 17) describes scattering that can, for example, emerge due electromagnetic coupling, see Section 5, or exciton-phonon coupling, see Section 6.The decay of the excitonic occupation with center-of-mass momentum Q k depends on the total decay rate Γ ξν,out Q k , which is given by a summation over all possible target states, in example, center-of-mass momenta Q k þ K k , of the microscopically calculated out-scattering rates The buildup of the excitonic occupation in Equation ( 17) is caused by the total in-scattering rate Γ ξν,in Qk .This rate typically depends on the respective scattering process and on other occupations, such as occupations in adjacent materials, see Section 5, or excitonic occupations with different momenta N ξν Section 6, as well as excitonic transitions It is furthermore common to theoretically divide N ξν Q k into an incoherent and a coherent part: which allows to couple it to Equation ( 13).However, in most cases, it is sufficient to take into account p ξν Q k ¼0 (lightcone).This allows to accurately describe excitation entering the equation by penetration with light, both resonantly [7,28] and off-resonant. [12]

Dipolar Coupling in TMDC Molecule Heterostructures
We assume one molecular layer l ¼ 1 at position z ¼ z 1 on top of a monolayer TMDC ðl ¼ 0Þ, see Figure 3a.We omit the respective layer index l in the following.We assume a dielectric environment as given in Figure 3a.For the TMDC, placed on a hBN substrate, we assume a dielectric constant due to the off-resonant transitions [95] and we consider the molecules to be in vacuum.All parameters can be found in the appendix.

Bloch Equations
In frequency domain ω, from Equations ( 4), (6), and (13), we can derive a set of Bloch equations for the interaction between molecule and TMDC mediated by the electric field, for the HOMO-LUMO transition σ HL ¼ H † L h iin the molecule and the excitonic transition p ξν Q k .We compute the TMDC Bloch equation as given in Equation ( 13), also analogously for the molecule, and selfconsistently solve the set of equations with the electric field, Equation ( 4), for details see ref. [16].We find a coupled set of equations Figure 3. a) Sketch of a TMDC (MoS 2 ) monolayer at z 0 with thickness Δ TMDC % 0.6nm, which is screened by its off-resonant transitions ϵ T ¼ 13.36, [95] on a hBN substrate with permittivity ϵ 1 ¼ 4.5 [95] covered by a molecular layer at z 1 in vacuum environment with ϵ 2 ¼ 1. b) Illustration of a Förster-type energy transfer between TMDC excitons (illustration in the exciton picture) and molecular excitation.c) Förster rate for a DCM molecule on a monolayer MoS 2 , at a distance of Δz ¼ 1nm, with an hBN substrate underneath, at low temperatures, that is, γ ξν % 0. The absorption of the pristine TMDC is given in light gray for comparison.Figure adapted from ref. [16].
The coupling between Equations ( 19) and ( 20) occurs in the near-field via the stationary version (ω !0) of the Green´s dyadic in Equation ( 5), that yields, (and V M!T Q k ξν analogously). [16]2.Förster Rates From Equations ( 19) and ( 20), a rate for the Förster process can be derived occurring as a decay channel for excitations σ HL of the molecular HOMO-LUMO transition.In frequency space, we find The rate is given by the imaginary part of the self-energy γ F ¼ ImðΣ F ðωÞÞ, and can be computed analytically in the limit γ ξν !0, which nicely approximates more exact results with γ ξν 6 ¼ 0. [92] In this limit, the rate depends on the transition energy E HL of the molecule with respect to the exciton energy E ξν , and is only nonzero when the energy of the molecular transition is larger than the bright TMDC 1s resonance (energy conservation) [16] In Equation (22), the in-plane momentum resembles the momentum for which the exciton E ξν Q k is in resonance with the molecular transition energy E HL , see Figure 3b.Obviously, at this point, the rate strongly depends on the molecular transition energy inducing a spectral dependence in the rate.The dielectric environment is taken into account via a Rytova-Keldysh-type approach, [96,97] and is expressed here with the abbreviations δ ϵ T þϵ 2 , respectively.d k stands for the absolute value of the in-plane part of the optical dipole moment of the TMDC. [85]t can be seen from Equation ( 22) that the Förster transfer depends on the squares of the optical dipole matrix elements d k φ ξν r k ¼0 , d HL of both constituents, as it is typical for this energy transfer. [50,98,99]Screening due to the dielectric environment occurs due to the corrections δ 1 , δ 2 .This includes the case of homogeneous dielectric environment ϵ 1 ¼ ϵ 2 ¼ ϵ T ≡ ϵ, where δ 1 ¼ 0 ¼ δ 2 results in the usual 1 ϵ -dependence of the Förster rate.For small distances, the rate shows a combination of exponential and power law decay as a function of the distance Δz between molecule and TMDC sheet.100] Figure 3c shows the Förster rate for negligible γ ξν % 0, which is a good approximation for low temperatures below the phonon threshold [92] for DCM molecules evaporated over a MoS 2 monolayer placed on a hBN substrate.The rate is plotted over the molecular transition energy E HL .It is strongest directly above the energy A and B resonances in the TMDC.As seen from Equation ( 22), finite-exciton momenta (constituting dark excitons) are necessary to allow for the energy-transfer process (although this is smeared out slightly for nonzero temperatures).Due to the ability of this near-field effect occurring from the momentum distribution of the molecular dipole which activates excitonic states with nonzero momentum, it becomes possible to resolve these momentum-dark states of the TMDC in the farfield.The signatures of dark excitons are directly experimentally accessible via the photoluminescence of the functionalizing molecules, see ref. [16].

Dipolar Coupling in TMDC MNP Heterostructures
[39][40][41][42][43][44] These observations motivate the development of a microscopic theory and its means to explore exciton localization near the nanoparticles.A schematic of a TMDC-MNP and a TMDC-PC heterostructure is given in Figures 4a and 5a, respectively.All parameters can be found in the appendix.As we are interested in the basics of resonant 1s exciton-plasmon interaction, we neglect higher excitonic states in the excitonic Bloch Equation ( 13), setting ν ¼ 1s.This is valid since the excitonic states are energetically separated [6] further than the plasmon broadening.
This way, we obtain a direct as well as a plasmon-mediated inter-valley coupling between K and K 0 -excitons in the excitonic Bloch equation.In the eigenbasis of the in-plane contribution of the Green's dyadic, Equation ( 5), with momentum-dependent eigenvectors e U Q k , e V Q k , e z , the inter-valley coupling, induced by the vector-mixing self-consistent electric field, is decoupled. [101]oreover, the electric near-field is only driving the V-component of the excitonic transition p The first term in the parentheses describes the exciton dispersion and includes direct exciton-exciton (exchange) interactions, X Q k .The second term incorporates effective, plasmon-induced exciton-exciton interactions, V Q k Q 0 k .The external field E 0 Q k at the TMDC as well as the MNP positions ðz 0 =z 1 Þ is included in the source term S Q k on the right-hand side of Equation ( 23), providing a direct and a plasmon-mediated excitation for TMDC excitons.

Localization of Excitons
Considering the coupled Equation ( 23), we find that a set of eigenfunctions ψ R,μ Q k for the center-of-mass motion can be defined via to extract the plexcitonic eigenvalues E μ .It describes the excitonic motion under the influence of an exchange coupling X Q k [82] in the plasmon-induced potential V Q k Q 0 k .For further investigation, the full excitonic Bloch Equation ( 23) will be mapped on these eigenstates using the projection where the wave functions are left-and right-handed solutions of the non-Hermitian eigenvalue problem in Equation (24).These eigenfunctions can be used to derive the polarization induced by the plasmon in real space that we display in Figure 4b.

Strong Coupling
To connect the exciton localization to related optical experiments that observe strong exciton-plasmon coupling, [41,[102][103][104][105] transmission T, reflection R, and true absorption A of TMDC-PC hybrids can be evaluated for a 2D PC, see Figure 5a, specified  for the case of perpendicular incident plane-wave excitation In the course of this, an extension of the quasi-static approximation is incorporated. [106]The absorption of a MoSe 2 -PC hybrid at room temperature is displayed in Figure 5b for different TMDC-PC distances Δz.Since the distance changes the interaction strength between excitons and plasmons, Figure 5b displays the transition between weak (red) and strong (blue) coupling.
For a close stack Δz ¼ r z þ 1nm, Figure 5b shows a typical strong coupling spectrum with a Rabi-splitting Ω that is several tens of meV.The splitting, which is a measure of the coupling strength, shrinks significantly with increasing distance between TMDC and PC.The trend is elaborated in the inset of Figure 5b, showing a Ω ∝ ðΔzÞ À2 dependency.This behavior differs from expectations for two spatially localized coupled dipoles, but is characteristic for a single dipole interacting with the polarization of a 2D plane.However, when maximizing the coupling strength in an experimental setup, for interlayer gaps ðΔz À r z Þ significantly smaller than 1 nm, [42] one has to additionally consider tunneling processes respectively charge transfer between the constituents which would suppress the strong exciton-plasmon interaction.
In summary, the developed framework allows to estimate the Rabi-splitting occurring for strong exciton-plasmon coupling.Since the equations are based on a microscopic description of the exciton dynamics in TMDCs, we can connect the strong coupling with the appearance of exciton localization.Thus, the theory complements previous models based on quasi-normal modes, [35,107,108] that focus more on a detailed description for the plasmonic nanoparticles photonics modes.

Dipole-Monopole Coupling of TMDC Graphene Heterostructures
So far, we investigated the dipolar coupling of atomically thin semiconductors with molecules and MNPs.In this section, we focus on the coupling between TMDC and graphene, illustrated in Figure 6a.The dielectric environment of the TMDC-graphene stack is assumed as displayed in Figure 6a.All parameters can be found in the appendix.We note that in addition to energy transfer between the sheets due to Coulomb coupling, also charge transfer via electron tunneling between the layers due to finite-electronic-wave function overlaps occurs.This electron tunneling is assisted by phonon scattering to match energy and momentum conservation during the tunneling event.The time scale for the tunneling transfer is on the order of about 100 fs. [51]As a new effect beyond the electron-tunneling process, the combination of a gapped semiconductor and a semimetal leads to the intriguing interfacial Meitner-Auger energy transfer corresponding to a dipole-monopole coupling, as observed in ref. [51].The process we address is sketched in Figure 6b.Here, the excitation energy of an optically excited exciton in the TMDC sheet is transferred to the graphene layer inducing an intraband transition, which creates hot-hole distributions.This kind of coupling is only possible for either transiently excited graphene leading to vacancies slightly below the Fermi level or p-doped graphene.We consider our description as a first step to illustrate the basic transfer mechanism.In principle, for large hole occupation difference, screening must be taken into account.While we discuss here the intraband excitation in the valence band, the mirrored process is possible also in the conduction band in case of n-doped graphene.
To study the Meitner-Auger energy transfer, Figure 6b, we derive the equation of motion for the incoherent (Q k 6 ¼ 0) exciton occupation in the K-valley N Q k ≡ N K1s Q k , see Equation ( 16) [51] ∂ The graphene valence band occupation is denoted by f k with momentum k, the linear graphene dispersion by ε k and the parabolic exciton dispersion by E Q k as introduced in Section 2.3.Equation ( 26) consists of an out-scattering from the TMDC, an in-scattering of graphene and the matrix element The occurrence of a single-dipole moment and the Coulomb potential V Q k , see Equation (10), shows that a single-dipole-monopole transition occurs.For a predominantly TMDC excitation as the initial occupation, comparing Equation (26) with the Boltzmann Equation ( 17) allows to identify the decay rate of the excitonic occupation as where the additional factor of 4 accounts for the valley and spin degree of freedom in graphene.Evaluating the energy conserving delta function, we find that k accounts for electrons close to the Dirac point and k À Q k for electrons deep in the valence band.To obtain an analytical expression, we assume that electrons close to the Dirac point have much smaller momenta than electrons deep in the valence band, that is, k ( Q k and jk À Q k j % Q k .A subsequent analytical evaluation of the sum yields with the density of states in graphene We can see that the rate depends on the matrix element W Q k of the Meitner-Auger transfer, on the density of states in graphene (via the k-sum), and on the occupation difference of the occupied states f Q k in graphene, which accounts for Pauli blocking.The Heaviside function accounts for the fact that during the interlayer transfer a minimal center-of-mass momentum is required to fulfill the energy and momentum conservation between the strongly different band curvature of a TMDC and graphene.
Figure 7a shows the calculated Meitner-Auger rate for a WSe 2 /graphene interface as function of center-of-mass momentum Q k , Equation ( 27), for different chemical potentials of the graphene occupation.As can be observed from Figure 7a due to the matrix element W Q k , the interfacial coupling vanishes for zero-excitonic center-of-mass momentum.Therefore, only incoherent excitons N Q k 6 ¼0 contribute to this type of coupling mechanism in contrast to coherent excitons p Q k ¼0 , Equation (13) [28,92,109,110] optically generated at Q k ¼ 0. For momenta close to 1.5 nm À1 , we observe a substantial transition rate of 1 meV, which increases with increasing chemical potential.A decreasing chemical potential leads to larger hole occupations close to the Dirac point, favoring the occupation difference in Equation ( 28), and therefore an increase of the transition rate.As discussed at the beginning, the chemical potential can either be set initially by doping or induced transiently due to the pulsed excitation.
Finally, we want to have a closer look at the minimal center-ofmass momentum necessary to obtain a finite-transfer rate.The amount of minimal required center-of-mass momentum can be controlled via the steepness of the graphene dispersion, that is, the Fermi velocity, and the TMDC optical bandgap, that is, the exciton resonance.For the chosen parameters (see appendix), Figure 7a shows a minimum required center-of-mass momentum of about 1.5 nm À1 .This corresponds to a kinetic energy of 100 meV, which far exceeds the mean kinetic energy at room temperature.However, the required momentum can also be provided by a hot TMDC exciton occupation N Q k 6 ¼0 induced for nonresonant optical pumping of the TMDC exciton, see Equation ( 26). [12]With increasing detunings above the exciton energy, the amount of injected excitons decreases, but they occupy larger momenta due an efficient scattering with acoustic and optical phonons. [12]Also for detuning below the optical bandgap, a substantial amount hot excitons with large centerof-mass momentum is formed, which can contribute to the interfacial Meitner-Auger energy transfer. [12,51]dditionally, we investigate the influence of the interlayer distance to the Meitner-Auger transfer and compare with the Förster coupling.To evaluate the z-dependence of the Meitner-Auger rate, we perform a thermal average [50] Γ the out-scattering rate.To evaluate the scattering rate, Equation ( 28), we set the out-scattering and in-scattering occupations to one and zero, respectively.Although, assuming thermalized electron occupations in graphene for a transient coupling mechanism is a strong approximation, this approach enables analytic access to the interlayer distance dependence of the Meitner-Auger coupling mechanism.Figure 7b shows the Δzdependence of Meitner-Auger transfer and for comparison also the result for the Förster transfer.First of all, we see that both coupling mechanisms decrease exponentially in the near-field.However, in the far-field, the Meitner-Auger transfer continues to decrease exponentially, while the Förster transfer changes to a Δz À4 dependence.Interestingly, the different far-field behavior can be traced back to the linear band structure of graphene, reflected by the Heaviside function in the scattering rate, and not the difference in the type of dipole-dipole/dipole-monopole interaction.
In this section, we presented a, to our knowledge, so far unstudied excitation transfer mechanism (Meitner-Auger process) uniquely occurring at semiconductor-semimetal interfaces. [51]In contrast to well-known energy-transfer mechanisms such as Förster or Dexter transfer or carrier transfer such as tunneling, the interfacial Meitner-Auger energy transfer enables an excitation transfer from the semiconductor to the semimetal inducing hot-carrier distributions in the semimetal.This unique coupling could advance conceptual devices as hot-carrier-based photocells. [111] Interlayer Exciton Relaxation in a MoSe 2 -WSe 2 Heterostructure In this section, we investigate the relaxation process of hot (Q k 6 ¼ 0) interlayer excitons in a TMDC heterostructure, which consists of a MoSe 2 layer and a WSe 2 layer vertically stacked on top of a SiO 2 substrate, see Figure 8a.All parameters can be found in the appendix.Hot interlayer excitons typically emerge by either pumping the MoSe 2 intralayer excitonic transition (favoring a subsequent hole transfer) [60] or by pumping the WSe 2 excitonic transition (electron transfer) [61] due to the type-II band alignment, see Figure 8b.They also form a microscopic dipolar field due to the interlayer electron-hole separation.
We employ the semiconductor Hamiltonian from Equation (10) starting from an initially injected hot exciton occupation and add the phonon Hamiltonian [112] adjusted to interlayer excitons Here, the first term denotes the free phonon Hamiltonian with the phonon dispersion ℏΩ l K k α with momentum K k , mode α, layer l, and the phonon annihilation (creation) operators b ð †Þ .The second term denotes the exciton-phonon interaction with the excitonic transition operators given in Equation (11)  and the exciton-phonon matrix element G l K k α , given by Here, g cðvÞ,MðWÞ K k α denotes the interaction matrix elements for electron/hole scattering (c=v).The compound index IL denotes As depicted in Figure 8b, we consider interlayer excitons with holes located at the K valley in the WSe 2 layer and electrons located at the K valley in the MoSe 2 layer.The phonon-assisted scattering process of electrons and holes as constituents of interlayer excitons is schematically displayed in Figure 8b.The excitonic wave functions φ IL are obtained by solving the Wannier equation for a static Coulomb potential V q k in Equation ( 10) in a five-layer geometry, see Figure 8a, as in ref. [113].With Equations ( 10) and (30), we are able to obtain the Boltzmann equation for the excitonic interlayer occupations [93] see Equation (17).The in-scattering rates are given by where n l K k α denotes the phonon occupation.The signs reflect phonon emission (þ) and absorption (À) processes.The outscattering rates can be obtained via the symmetry relation  [114] on top of a SiO 2 substrate (ϵ 1 ¼ 3.9).The gap and the region above the heterostructure are assumed as vacuum (ϵ 2 ¼ ϵ 3 ¼ 1).Δz denotes the distance between both layers.b) Sketch of the phonon-assisted relaxation processes of interlayer occupations (electron in MoSe 2 , hole in WSe 2 ), described by Equation (30).The parabolas denote the electron (À) and hole (þ) dispersion at the K valley in the respective layer over their momenta k.Hot excitons, that are, Coulomb-bound electron-hole pairs at high momenta, scatter down via optical and acoustic phonons to cold excitons, that are, Coulomb-bound electron-hole pairs at low momenta.c) Similar process as in (a), illustrated in the exciton picture for interlayer excitons.
We solve Equation ( 32) by assuming a Boltzmann distribution at the initial condition , where N IL 0 is the total exciton density within the linear regime, where phonon-assisted scattering dominates.Since we want to model hot interlayer excitonic occupations, we set the initial exciton temperature to T 0 ¼ 300K and assume a lattice temperature of T ¼ 77K.This temperature mismatch leads to a relaxation process mediated by optical and acoustic phonons until the temperature of the interlayer occupations equals the lattice temperature in equilibrium.In Figure 9a, we show the thermalization process of initially injected hot excitonic occupations.Within about two picoseconds (bright lines in Figure 9a), the excitonic occupations at higher momenta are depleted down to momenta or kinetic energies of about 1 nm À1 or 40 meV, respectively.Then, the scattering process significantly slows down until a Boltzmann distribution, whose temperature matches the lattice temperature of 77 K, is formed (dark line in Figure 9a).The entire thermalization process takes about 20-40 ps.The peculiar characteristics of the thermalization process will be examined in the following.
In Figure 9b, we display the decay rates of the thermalization process of the excitonic interlayer occupations N IL Q k in Figure 9a.We observe two features.
First, the decay rates are strongly non-monotonic and can be divided into two parts depending on the center-of-mass momentum Q k or the kinetic energy 2M IL : at high center-of-mass momenta Q k , the decay rate exhibits a moderate value, which increases with decreasing momentum until about 1 nm À1 or 40 meV, where a sharp edge at a high value is formed, after that the decay rate exhibits a substantial drop.Then, the rate increases again up to a moderate value at zero momentum.This edge constitutes an optical phonon bottleneck, that is, states of momenta smaller than the edge location that cannot be reached by optical phonons due to their large energy transfer.Thus, the scattering processes on the high-momentum side of the edge are dominated by optical phonons (fast downscattering in Figure 9a), whereas the scattering in the low-momentum range of the edge is solely caused by acoustic phonons, which are slow and exhibit a very small energy transfer (slow thermalization in Figure 9a).
The layer spacing Δz, see Figure 8a, is an important feature also for TMDC heterostructures.Since the exciton-phonon matrix elements in Equation ( 31) depend on the excitonic interlayer wave functions φ IL , the Δz dependence enters through the localization of the wave functions in q k -space: a decreasing layer spacing Δz weakens the localization of the interlayer wave function φ IL q k and vice versa, since an increasing layer spacing leads to an increasing charge separation and thus to a weakening of the Coulomb bonding of the interlayer exciton.In Figure 9b, we observe that the overall magnitude of the decay rate decreases monotonous by increasing layer spacing Δz.The cause is the localization of the excitonic interlayer wave function in momentum space: an increasing localization, caused by a growing layer spacing, sharpens the exciton-phonon matrix element in momentum space, see Equation (31), and therefore reduces the range of scattering partners.Thus, the larger the layer spacing, the longer the duration of the thermalization process displayed in Figure 9a becomes.
In this section, we discussed the layer spacing-dependent cooling and thermalization time of hot interlayer excitons, where the phonon-assisted decay time of hot interlayer excitons increases with the layer spacing.This can be important, if the TMDC heterostructure sample exhibits additional layer spacer.In addition, the relevant exciton-phonon coupling determines the relaxation of interlayer excitons into their ground state after their formation due to electron/hole transfer from intralayer excitons.

Conclusion
In conclusion, we have presented a general Maxwell Bloch equation framework to describe dipolar interactions at interfaces of Q k ranging from 0 fs to 40 ps, displaying the thermalization process of hot excitons, see Equation (32), for different interlayer spacings.The initial interlayer occupation is assumed as a hot Boltzmann distribution with a temperature of T 0 ¼ 300K, whereas the lattice temperature is assumed as T ¼ 77K.b) Decay rates Γ out Q k ¼ 33) and ( 34), of the thermalization process of the interlayer occupation in a MoSe 2 -WSe 2 heterostructure.
TMDC monolayers with a functionalization, for example, dye molecules, plasmonic nanoparticles, arrays of plasmonic nanoparticles, graphene, and other TMDC monolayers on the same footing.To focus on the electromagnetic near-field coupling, electron/hole-tunneling processes are not taken into account.
We have demonstrated that Förster interaction manifests an important relaxation pathway of optical dye-molecule excitation: tightly bound excitons in the TMDC monolayer modulate the dependence on the molecular transition energy.In particular, momentum-dark excitons induce sharp spectral features in the Förster rate, and therefore the functionalization with molecules makes these states observable in the far-field.
We further have revealed the appearance of bound excitonplasmon states in stacks of TMDC monolayers and plasmonic nanoparticles which arise due to dipole-dipole coupling of excitons and plasmons, so-called plexcitons.The framework, based on microscopic electron dynamics, allows us to describe the exciton localization near the plasmonic nanoparticles.At the same time, strong coupling in optical spectra is a fingerprint for exciton localization.Moreover, in stacks of TMDC monolayers with 2D PCs, we find strong coupling leading to a hybridization of exciton and plasmon dispersion with mode splittings on the order of some tens of meV.Due to the periodicity of the structure, this splitting can be observed in far-field detection.
Next, we have studied the energy transfer from WSe 2 excitons to graphene plasmon hole excitations, where we have proposed a new type of interlayer transfer, which is based on the Meitner-Auger interlayer coupling: here, an exciton in the WSe 2 layer decays via excitation of an intraband hole plasmon in graphene.Following our microscopic evaluation, this transfer process is fast compared to conventional energy-transfer mechanisms such as Dexter and Förster coupling in the considered structure and allows an excitation transfer from the TMDC inducing hotcarrier distributions in the graphene layer.
Last, we study the phonon-assisted relaxation dynamics of hot interlayer excitons in a TMDC heterobilayer with type-II band alignment: we find that the duration of the thermalization process of hot interlayer excitons grows by an increasing layer spacing.This behavior can be traced back to the layer spacingdependent Coulomb bonding of interlayer excitons and, thus, can be of crucial importance in heterostructure samples with additional layer spacer.
In summary, we have presented a general theoretical framework for the investigation of electromagnetic interactions in hybrid structures and we provided examples to demonstrate its capabilities.While the structures presented here are by no means exhaustive, they illustrate the potential of the framework to explore hybrid systems from different perspectives.Future work will involve more detailed investigations of the structures presented here, as well as the application of the framework to different hybrid structures or excitons at elevated densities.

Figure 2 .
Figure2.Approximated valence (blue) and conduction (red) bands ε λ¼v=c of molecules, graphene, and a TMDC.The molecular transition energy E HL is given by the energy of the HOMO and the LUMO.In graphene, valence and conduction bands touch at the K points and are approximately linear in their vicinity.TMDC monolayers feature a direct bandgap at the K=K ' points; their electronic dispersion is approximated parabolic.For simplicity, the TMDC is illustrated for one spin direction.

Figure 4 .
Figure 4. a) Sketch of a single MNP with semiaxis ðr x , r y , r z Þ at z 1 in a surrounding permittivity ϵ 2 on top of a TMDC monolayer at z 0 placed on a substrate with ϵ 1 .A spacer dielectric with permittivity ϵ 1 is inserted between the TMDC and the MNP.b) Plasmon-and light-induced polarization within TMDC layer in real space.

Figure 6 .
Figure 6.a) Sketch of a TMDC-graphene heterostructure with interlayer distance Δz and environment dielectric constants ϵ 1 , ϵ 2 , ϵ 3 .b) Dispersion of the TMDC in the exciton picture depending on the excitonic center-ofmass momentum Q k (left) and electronic graphene dispersion over the electron momentum k (right).The decay of an exciton in the TMDC induces an intraband transition in the graphene valence band.

Figure 7 .
Figure 7. a) Meitner-Auger transfer rate for different transient chemical potentials at room temperature.b) Distance dependence of Meitner-Auger transfer and Förster.We use v F ¼ 1.8 nm=f s and E 1s ¼ 1.65 eV and find a minimum required center-of-mass momentum of around 1.5 nm À1 with an effective exciton mass of M ¼ 0.65m e .

Figure 8 .
Figure 8. a) Sketch of a TMDC heterobilayer consisting of one MoSe 2 layer (ϵ M ¼ 17.4) and one WSe 2 layer (ϵ W ¼ 15.6)[114] on top of a SiO 2 substrate (ϵ 1 ¼ 3.9).The gap and the region above the heterostructure are assumed as vacuum (ϵ 2 ¼ ϵ 3 ¼ 1).Δz denotes the distance between both layers.b) Sketch of the phonon-assisted relaxation processes of interlayer occupations (electron in MoSe 2 , hole in WSe 2 ), described by Equation(30).The parabolas denote the electron (À) and hole (þ) dispersion at the K valley in the respective layer over their momenta k.Hot excitons, that are, Coulomb-bound electron-hole pairs at high momenta, scatter down via optical and acoustic phonons to cold excitons, that are, Coulomb-bound electron-hole pairs at low momenta.c) Similar process as in (a), illustrated in the exciton picture for interlayer excitons.

Figure 9 .
Figure 9. a) Snapshots of the interlayer occupations N ILQ k ranging from 0 fs to 40 ps, displaying the thermalization process of hot excitons, see Equation(32), for different interlayer spacings.The initial interlayer occupation is assumed as a hot Boltzmann distribution with a temperature of T 0 ¼ 300K, whereas the lattice temperature is assumed as T ¼ 77K.b) Decay rates Γ out Q k ¼