Theory of optically induced F\"orster Coupling in van der Waals coupled Heterostuctures

TL;DR

This paper models optically induced F"orster coupling in van der Waals heterostructures, predicting rapid energy transfer rates that align with experimental exciton lifetimes, advancing understanding of exciton dynamics in 2D materials.

Contribution

It provides a theoretical framework for F"orster coupling in heterostructures, predicting dephasing and energy transfer rates consistent with experimental data.

Findings

Dephasing rates and energy transfer times are in the picosecond range.

Predicted transition rate of thermalized excitons is about 4 ps$^{-1}$ at room temperature.

Theoretical results agree with recent experimental measurements.

Abstract

We investigate the impact of optically induced F\"orster coupling in van der Waals heterostructures consisting of graphene and a monolayer transition metal dichalcogenide (TMD). In particular, we predict the corresponding dephasing rates and a fast energy transfer between the TMD layer and graphene being in the picosecond range. Exemplary we find a transition rate of thermalized excitons of about 4 ps$^{-1}$ in a MoSe$_2$-graphene stack at room temperature. This timescale is in good agreement with the recently measured exciton lifetime in this heterostructure.

Tables

Table 1. Table 1: Parameters used in the computation. ∗ determined numerically by using the method given in Berghäuser and Malic ( 2014 ); Selig  et al. ( 2016 )

Param.		Param.		Ref.
$\hslash$	0.658 eV fs	$d_G$	0.25 e nm	Malic et al. (2011)
$e$	1 e	$v_F$	1 nm fs-1	Castro Neto et al. (2009)
${ϵ}_0$	5.5$\cdot$10-2 e2eV-1nm-1	$d_{W{S}_2}$	0.4 e nm	Selig et al. (2016)
$k_B$	8.6$\cdot$10-5 eV K-1	$\left|{\phi }_{W{S}_2}\left(𝐫=0\right)\right|$	0.49 nm-1	∗
${ϵ}_{Si{O}_2}$	3.9	$E_{W{S}_2}^{1s}$	2.0 eV	Christiansen et al. (2017)
		$M_{W{S}_2}$	3.5 eVfs2nm-2	Kormanyos et al. (2015)

Table 2. Table 2: Parameters for MoSe 2 as TMD material used in the computation. ∗ determined numerically by using the method given in Berghäuser and Malic ( 2014 ); Selig  et al. ( 2016 )

Param.		Ref.
$d_{MoS{e}_2}$	0.25 e nm	Selig et al. (2016)
$\left|{\phi }_{MoS{e}_2}\left(𝐫=0\right)\right|$	0.64 nm-1	∗
$E_{MoS{e}_2}^{1s}$	1.6 eV	Christiansen et al. (2017)
$M_{MoS{e}_2}$	6.1 eVfs2nm-2	Kormanyos et al. (2015)

Equations

$$H=\sum _{\mathbf{k},{\mathbf{k}}^{\prime },\lambda ,{\lambda }^{\prime }}\sum _{\mathbf{q},{\mathbf{q}}^{\prime },\nu ,{\nu }^{\prime }}{V}_{\mathbf{kq}{\mathbf{k}}^{\prime }{\mathbf{q}}^{\prime }}^{\lambda \nu {\lambda }^{\prime }{\nu }^{\prime }}{a}_{\mathbf{k}}^{†\lambda }{a}_{\mathbf{q}}^{†\nu }{a}_{{\mathbf{q}}^{\prime }}^{{\nu }^{\prime }}{a}_{{\mathbf{k}}^{\prime }}^{{\lambda }^{\prime }},$$

$$V_{\mathbf{kq}{\mathbf{k}}^{\prime }{\mathbf{q}}^{\prime }}^{\lambda \nu {\lambda }^{\prime }{\nu }^{\prime }}=⟨{\Psi }_{\mathbf{k}}^{\ast \lambda }(\mathbf{r}){\Psi }_{\mathbf{q}}^{\ast \nu }({\mathbf{r}}^{\prime })V(\mathbf{r}-{\mathbf{r}}^{\prime }){\Psi }_{{\mathbf{q}}^{\prime }}^{{\nu }^{\prime }}({\mathbf{r}}^{\prime }){\Psi }_{{\mathbf{k}}^{\prime }}^{{\lambda }^{\prime }}(\mathbf{r})⟩$$

$$V_{\mathbf{kq}{\mathbf{k}}^{\prime }{\mathbf{q}}^{\prime }}^{\lambda \nu {\lambda }^{\prime }{\nu }^{\prime }}=\frac{{\Omega }^2}{4\pi {ϵ}_0ϵA^2}\sum _{{\mathbf{R}}_{\mathbf{n}},{\mathbf{R}}_{\mathbf{n}}^{\prime }}{e}^{-i{\mathbf{R}}_{\mathbf{n}}(\mathbf{k}-{\mathbf{k}}^{\prime })}{e}^{-i{\mathbf{R}}_{\mathbf{n}}^{\prime }(\mathbf{q}-{\mathbf{q}}^{\prime })}\times$$

$$\times \left(\frac{{\mathbf{d}}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{\lambda {\lambda }^{\prime }}\cdot {\mathbf{d}}_{\mathbf{q}{\mathbf{q}}^{\prime }}^{\nu {\nu }^{\prime }}}{|{\mathbf{R}}_{\mathbf{n}}-{{\mathbf{R}}_{\mathbf{n}}}^{\prime }{|}^3}-3\frac{{\mathbf{d}}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{\lambda {\lambda }^{\prime }}\cdot ({\mathbf{R}}_{\mathbf{n}}-{\mathbf{R}}_{\mathbf{n}}^{\prime }){\mathbf{d}}_{\mathbf{q}{\mathbf{q}}^{\prime }}^{\nu {\nu }^{\prime }}\cdot ({\mathbf{R}}_{\mathbf{n}}-{\mathbf{R}}_{\mathbf{n}}^{\prime })}{|{\mathbf{R}}_{\mathbf{n}}-{\mathbf{R}}_{\mathbf{n}}^{\prime }{|}^5}\right).$$

$$V_{\mathbf{kq}{\mathbf{k}}^{\prime }{\mathbf{q}}^{\prime }}^{\lambda \nu {\lambda }^{\prime }{\nu }^{\prime }}=\frac{\Omega }{4\pi {ϵ}_0ϵA}\sum _{\mathbf{s},\mathbf{Q}}{e}^{-i\mathbf{Q}\cdot (\mathbf{s}+\mathbf{z})}{\delta }_{\mathbf{Q},\mathbf{k}-{\mathbf{k}}^{\prime }}{\delta }_{\mathbf{Q},\mathbf{q}-{\mathbf{q}}^{\prime }}\times$$

$$\displaystyle\times\bigg{(}\frac{\mathbf{d}^{\lambda\lambda^{\prime}}_{\mathbf{k\,k-Q}}\cdot\mathbf{d}^{\nu\nu^{\prime}}_{\mathbf{q\,q+Q}}}{|\mathbf{s}+\mathbf{z}|^{3}}-3\frac{\mathbf{d}^{\lambda\lambda^{\prime}}_{\mathbf{k\,k-Q}}\cdot(\mathbf{s}+\mathbf{z})\mathbf{d}^{\nu\nu^{\prime}}_{\mathbf{q\,q+Q}}\cdot(\mathbf{s}+\mathbf{z})}{|\mathbf{s}+\mathbf{z}|^{5}}\bigg{)}.$$

$$H_F=\sum _{\mathbf{k},\mathbf{q},\mathbf{Q}}\frac{1}{4\pi {ϵ}_0ϵA}{d}_{\mathbf{k}+\frac{1}{2}\mathbf{Q}k-\frac{1}{2}\mathbf{Q}}^{cv}{d}_{\mathbf{q}-\frac{1}{2}\mathbf{Q}q+\frac{1}{2}\mathbf{Q}}^{vc}a(\mathbf{Q},\mathbf{z})\times$$

$$\times a_{\mathbf{k}+\frac{1}{2}\mathbf{Q}}^{†c}{a}_{\mathbf{q}-\frac{1}{2}\mathbf{Q}}^{†v}{a}_{\mathbf{q}+\frac{1}{2}\mathbf{Q}}^c{a}_{\mathbf{k}-\frac{1}{2}\mathbf{Q}}^v+h.c.$$

$$H_F=$$

$$\times {\phi }_{\mathbf{k}}^{\ast \mu }{\phi }_{\mathbf{q}}^{\kappa }{P}_{\mathbf{Q}}^{†\mu }{R}_{\mathbf{Q}}^{\kappa }+h.c..$$

$$H_F=\sum _{\mu ,\mathbf{q},\mathbf{Q}}V(\mathbf{Q},\mathbf{z},\mu )P_{\mathbf{Q}}^{†\mu }{R}_{\mathbf{Q}}^{\mathbf{q}}+h.c.$$

$$i\hslash {\partial }_t{P}_{\mathbf{Q}}^{\nu }$$

$$i\hslash {\partial }_t{R}_{\mathbf{Q}}^{\mathbf{q}}$$

$${\partial }_t{P}_{\mathbf{Q}}^{\nu }=-\frac{i}{\hslash }(E_{\mathbf{Q}}^{\nu }+{\Delta }_{\mathbf{Q}}^{\nu }-i{\Gamma }_{\mathbf{Q}}^{\nu })P_{\mathbf{Q}}^{\nu }.$$

$${\Gamma }_{\mathbf{Q}}^{\nu }=|V(\mathbf{Q},\mathbf{z},\nu ){|}^2\pi \sum _{\mathbf{q}}\delta (E^{\mathbf{q}}-E_{\mathbf{Q}}^{\nu })$$

$$=\frac{d_T^2{d}_G^2|{\phi }^{\nu }(\mathbf{r}=0){|}^2{E}_{\mathbf{Q}}^{\nu }}{64{\pi }^2{ϵ}_0^2{ϵ}^2{\hslash }^2{v}_F^2}|a(\mathbf{Q},\mathbf{z}){|}^2,$$

$${\Delta }_{\mathbf{Q}}^{\nu }$$

$$=-\frac{d_T^2{d}_G^2|{\phi }^{\nu }(\mathbf{r}=0){|}^2{q}_{max}}{32{\pi }^3{ϵ}_0^2{ϵ}^2\hslash v_F}\times$$

$$\times \left(1+\frac{E_{\mathbf{Q}}^{\nu }}{2\hslash v_F{q}_{max}}ln|1-\frac{2\hslash v_F{q}_{max}}{E_{\mathbf{Q}}^{\nu }}|\right)|a(\mathbf{Q},\mathbf{z}){|}^2$$

$$\frac{1}{x^k}=\frac{1}{\Gamma (k)}{\int }_0^{\infty }dt{t}^{k-1}{e}^{-tx},$$

$$a(\mathbf{Q},z)=2\pi \frac{{\mathbf{e}}_1\cdot \mathbf{Q}{\mathbf{e}}_2\cdot \mathbf{Q}}{Q}{e}^{-Qz}$$

$$a(Q,z)=\pi Q{e}^{-Qz}{\mathbf{e}}_1\cdot {\mathbf{e}}_2$$

$$|a(Q,z){|}^2={\pi }^2{Q}^2{e}^{-2Qz}.$$

$$\displaystyle|a(T,z)|^{2}=\frac{\pi^{2}}{2}\frac{1}{\lambda^{6}}\bigg{(}2\lambda(T)^{2}(\lambda(T)^{2}+z^{2})$$

$$\displaystyle-e^{z^{2}/\lambda(T)^{2}}\lambda\sqrt{\pi}z(3\lambda(T)^{2}+2z^{2})erfc(z/\lambda(T))\bigg{)}$$

Full text

Since the discovery of graphene in 2004, a new research field on atomically thin quasi two dimensional (2D) materials has been established. In particular, monolayers of transition metal dichalcogenides (TMDs) became one of the most investigated materials beyond graphene. They exhibit technologically promising characteristics, such as a direct band gap, tightly bound electron-hole pairs (excitons)Mak et al. (2010); Ramasubramaniam (2012); Chernikov et al. (2014), strong light-matter interactionLi et al. (2014); Moody et al. (2015), and valley-selective circular dichroism Cao et al. (2012); Schmidt et al. (2016). Recent progress in growth techniques has also enabled the production of van der Waals heterostructures by vertically stacking atomically thin monolayers Shi et al. (2012); Geim and Grigorieva (2013). Consequently, heterostructures consisting of graphene and a monolayer TMD have been investigated theoretically Trushin (2018) and experimentally Georgiou et al. (2012); Britnell et al. (2013); He et al. (2014); Lin et al. (2014); Zhang et al. (2014); Hill et al. (2017); Froehlicher et al. (2018). Performing differential reflection measurements, J. He and co-workers found the relaxation of optically injected carriers from tungsten disulfide to graphene to occur within 1 psHe et al. (2014). In recent photoluminescence measurements G. Froelicher and co-workers found found a shortened exciton lifetime of about 1 ps in a coupled MoSe2-graphene structure at room temperature Froehlicher et al. (2018). In another study, Hill and co-workers found an additional broadening of 5 meV and a redshift of 23 meV of the excitonic resonance in a linear optics experiment due to the coupling between monolayer WS2 and grapheneHill et al. (2017). The experimental data have not been complemented by microscopic theory yet: Possible coupling mechanisms constitute of electronic tunneling, Dexter- and Förster energy transfer Richter et al. (2006); Rozbicki and Machnikowski (2008); Specht et al. (2015); Ovesen et al. (2018). Since both materials exhibit strong optical dipole moments, Förster coupling is expected to have a significant impact on the excitation transfer in these heterostructures, if the barrier potential is large enough to suppress electronic overlap Specht et al. (2015), cf. Fig. 1.

Förster-induced relaxation dynamics were investigated recently in heterostructures of different dimensions, namely between quantum dotsRichter et al. (2006); Rozbicki and Machnikowski (2008), between graphene and attached moleculesMalic et al. (2014) and between two quantum wellsBatsch et al. (1993). Here, we present a microscopic approach based on the Heisenberg equation of motion formalism allowing us to derive an analytic description of the Förster coupling in a van der Waals heterostructure consisting of graphene and a TMD monolayer. We predict the impact of the Förster coupling on the optical response of the heterostructure including the linewidth and spectral shift of excitonic resonances. In principle, also charge transfer (tunneling and Dexter processes) between both monolayers could affect the spectral width and the spectral position of the exciton line in the TMD monolayer. However, charge transfer does strongly dependent on the overlap of the wavefunctions of electrons in the TMD and graphene. For electronically decoupled monolayers, the Förster coupling can be expected to dominate tunnel and Dexter coupling Specht et al. (2015). For closely stacked heterostuctures however further investigation are required to evaluate the actual strength of the Dexter- and tunnel- coupling. This requires first principle methods and is addressed in future work. On the other hand, pure Förster coupling provides at least a limiting case that can be discussed in the analysis of experiments without or involving electronic overlap.

I Theoretical Model

We start with defining the many-particle Hamilton operator describing the Coulomb interaction between an exemplary TMD layer and graphene

[TABLE]

with electron annihilation (creation) operators $a^{(\dagger)\lambda}_{\mathbf{k}}$ acting on states with the band index $\lambda=c,v$ and the momentum $\mathbf{k}$. We denote electrons in the TMD monolayer by $(\mathbf{k}^{(^{\prime})},\lambda^{(^{\prime})})$ and in graphene by $(\mathbf{q}^{(^{\prime})},\nu^{(^{\prime})})$. The appearing Coulomb matrix element inducing the Förster transfer reads

[TABLE]

with the single particle electronic wavefunctions $\Psi_{\mathbf{k}}^{\lambda}$. We consider both materials to be aligned in the x-y-plane. Further, $V(\mathbf{r}-\mathbf{r^{\prime}})=\frac{e^{2}}{4\pi\epsilon_{0}\epsilon}\frac{1}{|\mathbf{r}-\mathbf{r^{\prime}}|}$ denotes the three dimensional Coulomb potential, with $\epsilon_{0}$ the vacuum permittivity and $\epsilon$ the mean dielectric constant of the surrounding material. This assumption is valid, since it differs only weakly from the exact static Coulomb potential for a dielectric surrounding typically found in experiments Hill et al. (2017), i.e. a TMD monlayer on a quartz substrate. This was checked carefully by explicitly evaluating the Possion equation Ovesen et al. (2018). Next, we decompose the space coordinates $\mathbf{r}=\mathbf{R_{n}}+\mathbf{r_{n}}$ in both constituents in $\mathbf{R_{n}}$ pointing to the center of the n-th unit cell (uc) and $\mathbf{r_{n}}$ describing a point within the unit cell. Thus, the integral transforms to $\int_{\mathbb{R}^{3}}d^{3}r=\sum_{\mathbf{R_{n}}}\int_{\text{uc}}d^{3}\mathbf{r_{n}}$. Here, the sum over the unit cells $\mathbf{R_{n}}$ is restricted to the monolayers, whereas the $\mathbf{r_{n}}$ integrations are performed in the 3 dimensional unit cells, cf. figure 1 (c). Now, we perform a Taylor expansion for the intra cell coordinate $\mathbf{r_{n}}$ and assume $|\mathbf{r_{n}-r_{n}^{\prime}}|\ll|\mathbf{R_{n}-R_{n}^{\prime}}|$ which can only be fulfilled, if the wavefunctions of the electrons in graphene and the TMD layer do not overlap. Here, $min|\mathbf{R_{n}}-\mathbf{R_{n}^{\prime}}|=z$ denotes the distance between both constituents. The monopole-monopole contribution of the Taylor expansion can be shown to contribute to the diagonal part of the Hamiltonian, providing an energy renormalization with respect to the uncoupled TMD-graphene heterostructure Richter et al. (2006). This energy renormalization vanishes in the linear optics limit and is therefore neglected in the following. The monopole-dipole can be neglected within the rotating-wave approximation of excitation with optical frequencies Specht et al. (2015) and therefore the dominating part is the dipole-dipole contribution of the Coulomb potential $\frac{1}{|\mathbf{R_{n}}+\mathbf{r_{n}}-\mathbf{R_{n}^{\prime}}-\mathbf{r_{n}^{\prime}}|}=\frac{\mathbf{r_{n}}\cdot\mathbf{r_{n}^{\prime}}}{|\mathbf{R_{n}}-\mathbf{R_{n}^{\prime}}|^{3}}-3\frac{\mathbf{r_{n}}\cdot(\mathbf{R_{n}}-\mathbf{R_{n}^{\prime}})\mathbf{r_{n}^{\prime}}\cdot(\mathbf{R_{n}}-\mathbf{R_{n}^{\prime}})}{|\mathbf{R_{n}}-\mathbf{R_{n}^{\prime}}|^{5}}.$

Assuming electronic Bloch functions in both structures, we obtain for the coupling element

[TABLE]

$\mathbf{d}^{\lambda\lambda^{\prime}}_{\mathbf{kk^{\prime}}}=\frac{e}{\Omega}\int_{\text{uc}}d^{3}r\,u^{*\lambda}_{\mathbf{k}}(\mathbf{r})\,\mathbf{r}\,u^{\lambda^{\prime}}_{\mathbf{k^{\prime}}}(\mathbf{r})$ denotes the dipole moment with the lattice-periodic Bloch factor $u^{\lambda}_{\mathbf{k}}(\mathbf{r})$ and the area of the unit cell $\Omega$. The Bloch factors further contain envelope functions ensuring the confinement in $z$-direction. Since they are normalized with respect to the $z$-integration, they drop in the computation.

To further evaluate the sums over $\mathbf{R_{n}}^{(^{\prime})}$, we introduce center-of-mass coordinates $\mathbf{S}=\frac{1}{2}(\mathbf{R_{n}}+\mathbf{R_{n}^{\prime}}+\mathbf{z})$ and $\mathbf{s}=\mathbf{R_{n}}-\mathbf{R_{n}^{\prime}}-\mathbf{z}$, with $\mathbf{z}$ being the vector pointing from the TMD layer to the graphene layer. Then, we obtain the matrix element:

[TABLE]

Without loss of generality, we write the dipole moments as $\mathbf{d}^{\lambda\lambda^{\prime}}_{\mathbf{k,k-Q}}=d^{\lambda\lambda^{\prime}}_{\mathbf{k,k-Q}}\mathbf{e}^{\lambda\lambda^{\prime}}_{1}$ and $\mathbf{d}^{\nu\nu^{\prime}}_{\mathbf{q,q+Q}}=d^{\nu\nu^{\prime}}_{\mathbf{q,q+Q}}\mathbf{e}^{\nu\nu^{\prime}}_{2}$ with $\mathbf{e}_{1/2}$ denoting the directions of the dipole moments in the TMD layer and in graphene, respectively. Before we provide the final Hamiltonian for the Förster coupling, we evaluate the sums over the electronic bands. We find, that contributions with $\lambda=\lambda^{\prime}$ and $\nu=\nu^{\prime}$ do not induce an interlayer energy transfer but an energy renormalization Specht et al. (2015). This renormalization vanishes in the limit of linear optics and therefore is dropped from the Hamiltonian. Additionally we perform a rotating frame approximation to remove non-energy conserving terms for excitations with optical frequencies from the Hamiltonian. Thus we obtain the Hamiltonian

[TABLE]

with $a(\mathbf{Q},\mathbf{z})=\int d^{2}s\frac{e^{-i\mathbf{Q}\cdot(\mathbf{s}+\mathbf{z})}}{|\mathbf{s}+\mathbf{z}|^{5}}(|\mathbf{s}+\mathbf{z}|^{2}\mathbf{e}^{cv}_{1}\cdot\mathbf{e}^{vc}_{2}-3\mathbf{e}^{cv}_{1}\cdot(\mathbf{s}+\mathbf{z})\mathbf{e}^{vc}_{2}\cdot(\mathbf{s}+\mathbf{z})).$

The next step is to introduce electron-hole pair polarizations $P^{\dagger cv}_{\mathbf{k}+\frac{1}{2}\mathbf{Q},\mathbf{k}-\frac{1}{2}\mathbf{Q}}=a^{\dagger c}_{\mathbf{k}+\frac{1}{2}\mathbf{Q}}a^{v}_{\mathbf{k}-\frac{1}{2}\mathbf{Q}}$ in the TMD and $R^{\dagger cv}_{\mathbf{q}+\frac{1}{2}\mathbf{Q},\mathbf{q}-\frac{1}{2}\mathbf{Q}}=a^{\dagger c}_{\mathbf{q}+\frac{1}{2}\mathbf{Q}}a^{v}_{\mathbf{q}-\frac{1}{2}\mathbf{Q}}$ in graphene to further simplify the Hamilton operator. Performing a transformation to center-of-mass coordinates and projecting the relative coordinate to exciton wavefunctions $\varphi^{\kappa}_{\mathbf{q}}$ in graphene and $\varphi^{\mu}_{\mathbf{k}}$ in the TMD layer with quantum numbers $\kappa$ and $\mu$, we obtain for the Hamiltonian

[TABLE]

Now, we assume that the dipole moments do not depend on the momentum at the band minimum ($d^{cv}_{\mathbf{k}+\frac{1}{2}\mathbf{Q}\,k-\frac{1}{2}\mathbf{Q}}\equiv d_{T}$ and $d^{vc}_{\mathbf{q}-\frac{1}{2}\mathbf{Q}\,q+\frac{1}{2}\mathbf{Q}}\equiv d_{G}$) and exploit the fact that at optical frequencies there are no bound excitons in graphene, thus we approximate $\varphi^{\kappa}_{\mathbf{q}}=\delta^{\kappa}_{\mathbf{q}}$ in the free particle limit. In contrast, for TMDs with strongly bound excitons, we use the relation $\sum_{\mathbf{k}}\varphi^{\mu}_{\mathbf{k}}=\sqrt{A}\varphi^{\mu}(\mathbf{r}=0)$ (where $\varphi^{\mu}(\mathbf{r})$ is the TMD exciton wavefunction in realspace). The TMD wavefunction is obtained as eigenvector of the Wannier equation Haug and Koch ; Berghäuser and Malic (2014), where we treated the appearing Coulomb potential within the Keldysh approach to account for the finite width of the TMD monolayerKeldysh (1978); Berkelbach et al. (2013). We obtain the Förster coupling Hamilton operator

[TABLE]

with the center of mass dependent Förster coupling element $V(\mathbf{Q},\mathbf{z},\mu)=\frac{a(\mathbf{Q},\mathbf{z})d_{T}d_{G}\varphi^{\mu}(\mathbf{r}=0)}{4\pi\epsilon_{0}\epsilon\sqrt{A}}$. The Förster Hamiltonian, eq. 7, can be interpreted as the annihilation of a bound exciton in the TMD monolayer and the creation of a free electron hole pair in the graphene monolayer, cf. figure 1 (b).

Now, we have all ingredients to define the Bloch equations for the exciton polarizations $P^{\nu}_{\mathbf{Q}},R^{\mathbf{q}}_{\mathbf{Q}}$ by exploiting the Heisenberg equation of motion

[TABLE]

The first term in both equations accounts for the free energy of the polarizations $E^{\nu}_{\mathbf{Q}}=E^{\nu}+\frac{\hbar^{2}\mathbf{Q}^{2}}{2M}$ with $E^{\nu}$ being the spectral position of the TMD exciton and $M=m_{e}+m_{h}$ as the total mass of the exciton. In graphene, we assume the energy to be independent of the center-of-mass momentum $E_{\mathbf{Q}}^{\mathbf{q}}=\hbar v_{F}(|\mathbf{q}+\frac{1}{2}\mathbf{Q}|+|\mathbf{q}-\frac{1}{2}\mathbf{Q}|)\approx 2\hbar v_{F}|\mathbf{q}|$ with $v_{F}$ as the Fermi velocity, since typically it holds $|\mathbf{q}|\gg|\mathbf{Q}|$ and the dispersion of graphene is linear in the investigated $\mathbf{q}$-region for coherent or even thermalized TMD excitons. The remaining terms in eqs. (8) and (9) result from the Förster coupling between both monolayers. The coupled eqs. (8) and (9) describe new quasiparticles from electron-hole pairs in graphene and excitons in TMDs. The most direct approach in the weak coupling limit is to perform a Born-Markov approximationKuhn and Rossi (1992); Schilp et al. (1994); Scully and Zubairy (1997) for the polarization in graphene $R^{\mathbf{q}}$. This approximation is valid as long as the Förster induced interlayer coupling is small compared to the spectral bandwidth provided by the broad electronic energy distribution in graphene. The Born Markov approximation carried out is formally equivalent to a Wigner-Weisskopf approximation where the graphene continuum acts as a bathScully and Zubairy (1997). This way, we have access to the dephasing rate and the energy renormalization for excitons in the TMD due to Förster coupling with the graphene layer:

[TABLE]

For the broadening we obtain

[TABLE]

while the energy renormalization reads

[TABLE]

Here, we neglect the Förster interaction mediated influence of different TMD excitons on each other. For an optical excitation, we have $E^{\nu}_{\mathbf{Q}}\approx\hbar\omega_{opt}$, where $\omega_{opt}$ is the incident laser frequency. Note that we added a factor of $2$ to take account to the spin degree of freedom in graphene. Note that the momentum sum appearing in the expression for the lineshift, equation 13, in general diverges. This problem is also known from the computation of the Lambshift, where the atomic states couple to the mode continuum of the radiational field which results in a divergent self energy contributionScully and Zubairy (1997). Similar, in the current case, the energy renormalization does depend on the choice of a physically motivated momentum cutoff $q_{max}$. Therefore we do not compute the energy renormalization explicitly but will discuss it qualitatively in the following. We find that the dephasing rate, eq. (12), and the energy renormalization, eq. (14), show the same momentum and layer separation dependence, since both are proportional to $a(\mathbf{Q},\mathbf{z})|^{2}$.

The appearing function $|a(\mathbf{Q},\mathbf{z})|^{2}$ in eqs. 12 and 14 has to be evaluated. $a(\mathbf{Q},\mathbf{z})$ can be integrated by Schwinger parameterizing the denominator of the integrand Prausa (2017)

[TABLE]

with $x,k\in\mathbb{R}$ and $\Gamma(k)$ denoting the Gamma function. Inserting this expression in the integral $a(\mathbf{Q},z)$ the $d^{2}s$ integration can be performed straight forward. Also the $dt$ integration over hermite gaussian polynomials

[TABLE]

Our result coincides nicely with the result given in reference Batsch et al., 1993. Note that we obtain a different expression compared to reference Tomita et al., 1996. Since we are interested in the general behavior of the Förster rate, the dependence on the angle of the center of mass momentum is averaged out, which yields

[TABLE]

The last step is to square this expression and sum it over $K$ and $K^{\prime}$ point in graphene. The latter is implicitly included in the summation over the momentum $\mathbf{q}$ in equation 13 and 11. Note that we have to sum over two orthogonal dipole moments in graphene Stroucken et al. (2011), and therefore the angular dependence in the final result drops. We obtain the final expression

[TABLE]

Using eqs. 12 and 14, the corresponding dephasing rates and energy renormalizations can be evaluated quantitatively.

II Results

We exploit the derived equations to calculate Förster induced dephasing rate for the case of a van der Waals heterostructure consisting of graphene and monolayer tungsten disulfide (WS2) as an exemplary TMD on a quartz substrate ($\epsilon=$3.9). In particular, we exploit Eqs. (12) with parameters given in Table 1 and compute the dephasing rate in the WS2 monolayer in the presence of the graphene layer.

In Figure 2 (a), we show a contour plot of the Förster induced dephasing rate as a function of the center-of-mass momentum $\mathbf{Q}$ of the TMD exciton and the distance between the TMD and the graphene layer. We find dephasing rates ranging from more than 5 meV at $\mathbf{Q}=$2 nm*-1* and $z=$0.5 nm (relevant for the generation of incoherent excitons for photoluminescence Thränhardt et al. (2000)), to 0 meV for $\mathbf{Q}=0$ (relevant for coherent optical absorption) and $z=$0.5 nm (corresponding to the situation where the WS2 layer and the graphene layer are closely stackedHill et al. (2017)). As expected, the dephasing rate decreases as a function of the interlayer spacing.

To investigate the dependence on the center-of-mass momentum in more detail, we show cuts at different interlayer separations of 0.5 nm, 1.0 nm and 2 nm in Fig. 2 (b). We find for all interlayer separations an initial increase of the Förster rate followed by an exponential decay, which can be directly extracted from eq. (16). Excitons, visible in coherent optical experiments (without any scattering or thermalization), exhibit vanishing center-of-mass momenta $\mathbf{Q}\approx 0$, due to momentum conservation Kira and Koch (2006); Haug and Koch ; Berghäuser and Malic (2014). Our computations predict no impact of the Förster coupling on the linewidth and resonance position of coherent excitons at $\mathbf{Q}=0$. In reference Hill et al., 2017, H. Hill and co-workers found an additional broadening of the WS22 resonance in a coherent reflectance measurement of about 5 meV compared to the monolayer. Since the Förster dephasing rate vanishes for $\mathbf{Q}=0$ it can be ruled out as a direct source of the peak broadening which requires ongoing investigations.

In our computations, we find a strong impact of the dielectric constant of the enviroment due to the screening of the Coulomb potential and exciton wavefunction. The latter effect turns out to be less significant. To illustrate this, we discuss briefly two limits. Exemplary, for a free standing heterostructure ($\epsilon=$1) with an interlayer separation of 0.5 nm we find a Förster induced dephasing of about 10 meV at $Q=$2 nm*-1*. In contrast, for a heterostructure encapsulated by hexagonal boron nitride ($\epsilon=$7) we find a Förster rate of approximately 1.5 meV at $Q=$2 nm*-1*.

Next we investigate the dependence of the Förster induced dephasing rate on the interlayer separation. In Fig. 2 (c) we show the dephasing rate as a function of the interlayer distance for three different center-of-mass momenta $Q$. First, we find for all $Q$ a decreasing behavior. For elevated center of mass momenta, the Förster induced dephasing rate exhibits an exponential decay which is depicted in figure 2 (c).

To get a first impression of the Förster transition rate of a thermalized exciton distribution in the WS2 monolayer into the graphene layer, we apply a thermal average. In the low excitation limit the equation for the incoherent exciton density $N_{\mathbf{Q}}=\delta\langle P^{\dagger}_{\mathbf{Q}}P_{\mathbf{Q}}\rangle$ Selig et al. (2018) reads : $\partial_{t}N_{\mathbf{Q}}=-2\Gamma_{\mathbf{Q}}N_{\mathbf{Q}}$. Here, supported by the fast relaxation within the graphene layer Winzer et al. (2010), the occupation in graphene was assumed to vanish. Under the assumption of an initially thermalized exciton Boltzmann distribution $N_{\mathbf{Q}}=Ne^{-\frac{E_{\mathbf{Q}}}{k_{B}T}}$ in the TMD, the temperature dependent transition rate can be obtained by thermally averaging the momentum dependent transition rate $\Gamma_{T}=2\langle\Gamma_{\mathbf{Q}}^{1s}\rangle_{T}$. This assumption is justified by the fast exciton-phonon scattering rates which mediate the thermalization, being in the order of some 10 fsSelig et al. (2018). The corresponding thermal average of the function $|a(T,z)|^{2}=\langle|a(\mathbf{Q},z)|^{2}\rangle_{T}$ reads

[TABLE]

with the error function $\rm{erfc}(z/\lambda(T))$ and the thermal wavelength $\lambda(T)=\frac{\hbar}{\sqrt{2Mk_{B}T}}$. Interestingly, in the second term a high numerical accuracy is required in the evaluation since the function $e^{a^{2}}$ is fast increasing and the function $erfc(a)$ is fast decreasing as a function of $a$ Tomita et al. (1996).

Figure 3 shows the transition rate as a function of temperature for three selected interlayer separations $z$. We find for all $z$ an increase of the rate with increasing temperature. Exemplary, we predict a rate of 5 ps*-1* for a closely stacked heterostructure ($z=$0.5 nm) at room temperature. A similar time scale was reported for the transition of carriers from WS2 to graphene for a closely stacked heterostructure from non-linear differential reflectance contrast measurements (pump-probe)He et al. (2014). The incoherent exciton occupation is visible in coherent reflectance contrast measurements as a bleaching of the excitonic polarization $i\hbar\partial_{t}P_{\mathbf{0}}=EP_{\mathbf{0}}+(1-2\sum_{\mathbf{Q}}\Xi_{\mathbf{Q}}N_{\mathbf{Q}})\mathbf{d}\cdot\mathbf{E}_{T}$, with the excitonic form factor $\Xi_{\mathbf{Q}}=\sum_{\mathbf{q}}|\varphi_{\mathbf{q}}|^{2}|\varphi_{\mathbf{q+\frac{1}{2}Q}}|^{2}$ during and after the thermalization Hawton and Nelson (1998); Katsch et al. (2018). Our computation reveals a value of 4 ps*-1* for a heterostructure consisting of monolayer MoSe2 and graphene (obtained with parameters from table 2). This value also coincides qualitativley with the recently reported exciton lifetime of about 1 ps in a MoSe2-graphene stackFroehlicher et al. (2018).

In figure 3 (a), the increasing Förster induced transition rate as a function of temperature can be explained as follows: At low temperatures, the exciton thermalizes in narrow Boltzmann distributions at very small momenta with a width of $kT$. For this distribution we find a vanishing Förster rates, because the Förster coupling as a function of the center of mass momentum $\mathbf{Q}$ vanishes for $\mathbf{Q}=0$. Increasing the temperature leads to a broadening of the exciton distribution in momentum space. The leads to the occupation of exciton states with larger Förster rates, cf. Fig. 2 (a) and (b). This results in a increase of the Förster transition rate at elevated temperatures. However, since even at room temperature most of the excitons are located at low momenta, the transition rate is still increasing at elevated temperatures.

Finally, in Fig. 3 (b), we show the transition rate as a function of the interlayer separation for three different temperatures. For small distances we find a $e^{-z}$ law for all investigated temperatures, which is consistent with the observations for the dephasing rate. At larger interlayer spacings we find a $z^{-4}$ law, which is consistent with the limits of the function in equation 19. The observed $z$ dependence is consistent with the observations for Förster rates in porpyrin-functionalized graphene. Malic et al. (2014) The developed theoretical approach can be applied to other van der Waals heterostructures including multilayers of the same TMD or stacks of different TMD monolayersOvesen et al. (2018).

III Conclusion

In conclusion, we presented a simple analytic model describing the Förster mechanism in van der Waals heterostructures consisting of graphene and a monolayer TMD. We predict Förster rates leading to a fast energy transfer between the TMD and graphene on a picosecond timescale. Exemplary the Förster induced transition rate of thermalized WS2 excitons is about 5 ps*-1* at room temperature. This was found to nicely coincide with recent pump probe He et al. (2014) and photoluminescence Froehlicher et al. (2018) measurements. So far, we have not included recently investigated dark exciton states with momenta far beyond the light coneQiu et al. (2015); Wu et al. (2015); Selig et al. (2016, 2018), which will be addressed in future work.

Acknowledgement

In particular we thank the referee of the manuscript for a fruitful discussion. The authors gratefully acknowledge stimulating discussions with Gunnar Berghäuser (Chalmers, Göteborg), Manuel Kraft and Florian Katsch (TU Berlin). This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), Projektnummer 182087777–SFB 951 (projects A8 and B12, N.K. and A.K.) and Projektnummer 43659573–SFB 787 (project B1, M.S.). This project has also received funding from the European Unions Horizon 2020 research and innovation programme under Grant Agreements No. 696656 (Graphene Flagship, E.M.) and No. 734690 (SONAR, A.K.). Further we acknowledge funding from the Swedish Research Council and Stiftelsen Olle Engkvist (E.M.). This work was partially supported by the National Research Foundation of Korea (NRF) Grants funded by the Korean Government (MSIP) (2018R1A2B6001449).