Phonon-assisted Photoluminescence from Dark Excitons in Monolayers of Transition Metal Dichalcogenides

TL;DR

This paper develops a microscopic theory to explain phonon-assisted photoluminescence from dark excitons in transition metal dichalcogenides, clarifying emission features below the bright exciton line in tungsten and molybdenum-based materials.

Contribution

It introduces a comprehensive microscopic model for phonon-assisted dark exciton recombination, explaining experimental PL spectra in TMD monolayers.

Findings

WSe₂ shows visible low-temperature PL from dark exciton recombination.

Different configurations of bright and dark states explain spectral differences.

The theory aligns well with experimental observations.

Abstract

The photoluminescence (PL) spectrum of transition metal dichalcogenides (TMDs) shows a multitude of emission peaks below the bright exciton line and not all of them have been explained yet. Here, we study the emission traces of phonon-assisted recombinations of momentum-dark excitons. To this end, we develop a microscopic theory describing simultaneous exciton, phonon and photon interaction and including consistent many-particle dephasing. We explain the drastically different PL below the bright exciton in tungsten- and molybdenum-based materials as result of different configurations of bright and dark states. In good agreement with experiments, we show that WSe$_2$ exhibits clearly visible low-temperature PL signals stemming from the phonon-assisted recombination of momentum-dark excitons.

Equations

$$H=\sum$$

$$\frac{d}{dt}{n}_{\mathbf{k}}$$

$$i\hslash \frac{d}{dt}{\mathcal{S}}_{\mathbf{k}}^{\mu }$$

$$i\hslash \frac{d}{dt}{\mathcal{U}}_{\mathbf{kq}}^{\nu ,\pm }$$

$$i\hslash \frac{d}{dt}{\mathcal{C}}_{\mathbf{kq}}^{\nu \mu ,\pm }$$

$$\displaystyle I_{\sigma}(\omega)=\dfrac{2}{\hbar}\sum_{\mu}\dfrac{|M^{\mu}_{\sigma}|^{2}}{(E^{\mu}_{\mathbf{0}}-\omega)^{2}+(\gamma^{\mu}_{\sigma\mathbf{0}}+\Gamma^{\mu}_{\mathbf{0}})^{2}}\bigg{(}\gamma^{\mu}_{\sigma\mathbf{0}}N^{\mu}_{\mathbf{0}}+\sum_{\nu\mathbf{q},\alpha\pm}|D^{\mu\nu}_{\alpha\mathbf{q}}|^{2}N^{\nu}_{\mathbf{q}}\eta^{\pm}_{\alpha\mathbf{q}}\dfrac{\Gamma^{\nu}_{\mathbf{q}}}{(E^{\nu}_{\mathbf{q}}\pm\Omega^{\alpha}_{\mathbf{q}}-\omega)^{2}+(\Gamma^{\nu}_{\mathbf{q}})^{2}}\bigg{)}.$$

$$I_{\sigma }(\omega ){\text{arrowvert}}_{\omega \approx E_0}=\frac{2}{\hslash }\sum _{\mu }\frac{|M_{\sigma }^{\mu }{|}^2({\gamma }_{\sigma 0}^{\mu }+{\Gamma }_0^{\mu })}{(E_0^{\mu }-\omega {)}^2+({\gamma }_{\sigma 0}^{\mu }+{\Gamma }_0^{\mu }{)}^2}{N}_0^{\mu }.$$

$$H$$

$$\frac{d}{dt}{n}_{\mathbf{k}}$$

$$i\hslash \frac{d}{dt}{\mathcal{S}}_{\mathbf{k}}$$

$$i\hslash \frac{d}{dt}{\mathcal{U}}_{\mathbf{kq}}^{\zeta }$$

$$i\hslash \frac{d}{dt}{\mathcal{C}}_{\mathbf{Qq}}^{\zeta }$$

$$\frac{d}{dt}{N}_{\mathbf{Q}}$$

$$⟨c_{\mathbf{k}}^{†}{b}_{\mathbf{q}}^{†}{b}_{{\mathbf{q}}^{\prime }}{a}_{\mathbf{q}-{\mathbf{q}}^{\prime }}⟩\approx ⟨c_{\mathbf{k}}^{†}{a}_{\mathbf{q}-{\mathbf{q}}^{\prime }}⟩⟨b_{\mathbf{q}}^{†}{b}_{{\mathbf{q}}^{\prime }}⟩+⟨c_{\mathbf{k}}^{†}{b}_{\mathbf{q}}^{†}{b}_{{\mathbf{q}}^{\prime }}{a}_{\mathbf{q}-{\mathbf{q}}^{\prime }}{⟩}^{\text{corr.}}$$

$$i\hslash \frac{d}{dt}{X}_{\mathbf{kq}{\mathbf{q}}^{\prime }}^{\zeta {\zeta }^{\prime }}=(E_{\mathbf{q}+{\mathbf{q}}^{\prime }}-\zeta {\Omega }_{\mathbf{q}}-{\zeta }^{\prime }{\Omega }_{{\mathbf{q}}^{\prime }}-{\omega }_{\mathbf{k}})X_{\mathbf{kq}{\mathbf{q}}^{\prime }}^{\zeta {\zeta }^{\prime }}+D_{{\mathbf{q}}^{\prime }}^{\ast }{\eta }_{{\mathbf{q}}^{\prime }}^{{\zeta }^{\prime }}{\mathcal{U}}_{\mathbf{kq}}^{\zeta }+\dot{X}{|}_{\text{source}}+\dot{X}{|}_{\text{corr.}}.$$

$$X_{\mathbf{kq}{\mathbf{q}}^{\prime }}^{\zeta {\zeta }^{\prime }}\approx -\frac{D_{{\mathbf{q}}^{\prime }}^{\ast }{\eta }_{{\mathbf{q}}^{\prime }}^{{\zeta }^{\prime }}{\mathcal{U}}_{\mathbf{kq}}^{\zeta }}{E_{\mathbf{q}+{\mathbf{q}}^{\prime }}-\zeta {\Omega }_{\mathbf{q}}-{\zeta }^{\prime }{\Omega }_{{\mathbf{q}}^{\prime }}-{\omega }_{\mathbf{k}}-i{\delta }^{+}}\approx -\frac{D_{{\mathbf{q}}^{\prime }}^{\ast }{\eta }_{{\mathbf{q}}^{\prime }}^{{\zeta }^{\prime }}{\mathcal{U}}_{\mathbf{kq}}^{\zeta }}{E_{\mathbf{q}+{\mathbf{q}}^{\prime }}-E_{\mathbf{q}}-{\zeta }^{\prime }{\Omega }_{{\mathbf{q}}^{\prime }}-i{\delta }^{+}},$$

$$E_{\mathbf{Q}}⟶{\tilde{E}}_{\mathbf{Q}}=E_{\mathbf{Q}}+{\Sigma }_{\mathbf{Q}}^{\text{ph}};{\Sigma }_{\mathbf{Q}}^{\text{ph}}=\sum _{\zeta \mathbf{q}}\frac{|D_{\mathbf{q}}{|}^2{\eta }_{\mathbf{q}}^{\zeta }}{{\tilde{E}}_{\mathbf{Q}}-{\tilde{E}}_{\mathbf{Q}+\mathbf{q}}+\zeta {\Omega }_{\mathbf{q}}},$$

$$E_0⟶E_0-i{\gamma }_0;{\gamma }_0=\frac{1}{i}\sum _{\mathbf{k}}\frac{|M_{\mathbf{k}}{|}^2}{E_0-{\omega }_{\mathbf{k}}-i\kappa }\approx \pi \sum _{\mathbf{k}}|M_{\mathbf{k}}{|}^2\delta (E_0-{\omega }_{\mathbf{k}}),$$

$$G_{\mathbf{Q}}(\omega )=(E_{\mathbf{Q}}-\omega {)}^{-1};{\mathcal{G}}_{\mathbf{Q}}(\omega )=({\tilde{E}}_{\mathbf{Q}}-\omega {)}^{-1},$$

$${\mathcal{S}}_{\mathbf{k}}$$

$${\Sigma }^{\text{ph}}(\omega )=\sum _{\zeta \mathbf{q}}{\mathcal{G}}_{\mathbf{q}}(\omega +\zeta {\Omega }_{\mathbf{q}})|D_{\mathbf{q}}{|}^2{\eta }_{\mathbf{q}}^{\zeta }\text{and}\frac{G_0(\omega +i{\gamma }_0)}{1+{\Sigma }^{\text{ph}}(\omega )G_0(\omega +i{\gamma }_0)}\approx {\mathcal{G}}_0(\omega ),$$

$$\displaystyle\mathcal{G}_{\mathbf{q}}(\omega+\zeta\Omega_{\mathbf{q}})\mathcal{G}_{\mathbf{q}}(\tilde{E}^{\ast}_{0}+\zeta\Omega_{\mathbf{q}})=-\mathcal{G}^{\ast}_{\mathbf{0}}(\omega)\bigg{(}\mathcal{G}_{\mathbf{q}}(\omega+\zeta\Omega_{\mathbf{q}})-\mathcal{G}_{\mathbf{q}}(\tilde{E}^{\ast}_{0}+\zeta\Omega_{\mathbf{q}})\bigg{)}$$

$$\frac{d}{dt}{n}_{\mathbf{k}}=\frac{2}{\hslash }|M_{\mathbf{k}}{|}^2$$

$$\frac{d}{dt}{N}_0{\text{arrowvert}}_{\text{phon.}}$$

$${\Sigma }_{\mathbf{Q}}^{\text{ph}}\approx -i{\Gamma }_{\mathbf{Q}}=-i\pi \sum _{\mathbf{q},\zeta }|D_{\mathbf{q}}{|}^2{\eta }_{\mathbf{q}}^{\zeta }\delta (E_{\mathbf{Q}}-E_{\mathbf{Q}+\mathbf{q}}+\zeta {\Omega }_{\mathbf{q}}),$$

$$\displaystyle\dfrac{d}{dt}n_{\mathbf{k}}=\dfrac{2}{\hbar}\dfrac{|M_{\mathbf{k}}|^{2}}{(E_{\mathbf{0}}-\omega_{\mathbf{k}})^{2}+(\gamma_{\mathbf{0}}+\Gamma_{\mathbf{0}})^{2}}\bigg{(}\gamma_{\mathbf{0}}N_{\mathbf{0}}+\sum_{\mathbf{q}\zeta}|D_{\mathbf{q}}|^{2}N_{\mathbf{q}}\eta^{\zeta}_{\mathbf{q}}\dfrac{\Gamma_{\mathbf{q}}}{(E_{\mathbf{q}}+\zeta\Omega_{\mathbf{q}}-\omega_{\mathbf{k}})^{2}+(\Gamma_{\mathbf{q}})^{2}}\bigg{)}$$

$$D_{\lambda \mathbf{q}}^{\nu \mu }=\sum _{\mathbf{k},\alpha =e,h}{g}_{\lambda \mathbf{q}}^{\alpha \mathbf{k}}{\Phi }_{\nu }^{\ast }(\mathbf{k}){\Phi }_{\mu }(\mathbf{k}+{\mathbf{q}}_{\alpha }).$$

$${\gamma }_{\sigma 0}^{\nu }=\pi \sum _{\mathbf{k}}|M_{\sigma \mathbf{k}}^{\nu }{|}^2\delta (E_0^{\nu }-{\omega }_{\mathbf{k}})=\frac{\hslash e_0^2}{2m_0^2{ϵ}_{0}n{c}_0}|{\mathbf{m}}^{vc}\cdot {\mathbf{e}}_{\sigma }{|}^2\frac{|{\tilde{\Phi }}_{\nu }(r=0){|}^2}{E_0^{\nu }},$$

$$|{\mathbf{m}}^{vc}\cdot {\mathbf{e}}_{\sigma }{|}^2=\frac{1}{2}[\frac{a_0{m}_{0}t}{\hslash }(1+\sigma ){]}^2.$$

$${ϵ}_{scr}(q)={\kappa }_1\tanh (\frac{1}{2}[{\alpha }_{1}dq-\ln (\frac{{\kappa }_1-{\kappa }_2}{{\kappa }_1+{\kappa }_2})]),$$

$$({\epsilon }_{c,\mathbf{k}+\alpha \mathbf{Q}}-{\epsilon }_{v,\mathbf{k}-\beta \mathbf{Q}}){\Phi }_{\nu \mathbf{Q}}(\mathbf{k})-\sum _{\mathbf{q}}{W}_{\mathbf{q}}{\Phi }_{\nu \mathbf{Q}}(\mathbf{k}+\mathbf{q})=E_{\nu \mathbf{Q}}{\Phi }_{\nu \mathbf{Q}}(\mathbf{k}).$$

$$\frac{{\hslash }^2{k}^2}{2m_{\text{r}}}{\Psi }_{\nu }(\mathbf{k})-\sum _{\mathbf{q}}{W}_{\mathbf{q}}{\Psi }_{\nu }(\mathbf{k}+\mathbf{q})=E_{\nu }^{\text{bind}}{\Psi }_{\nu }(\mathbf{k}).$$

Full text

The cryogenic photoluminescence (PL) spectrum of two dimensional semiconductors provides a powerful tool to study intriguing quantum phenomena invisible at room temperature or in linear optical experiments. The large variety of low temperature emission features at energies below the bright exciton resonance indicates the existence of bound exciton configurations, such as trions, biexcitons Mak et al. (2013); Wang et al. (2014); Courtade et al. (2017); Ye et al. (2018a); Barbone et al. (2018) and trapped excitons Tonndorf et al. (2015); Zhang et al. (2017), but could also result from the indirect recombination of intrinsically dark exciton states Lindlau et al. (2017); Koperski et al. (2017); Mueller and Malic (2018); Malic et al. (2018); Selig et al. (2016); Feierabend et al. (2017). In particular, the radiative decay of a momentum indirect electron-hole pair requires the simultaneous interaction with a phonon to fulfill the momentum conservation and is therefore very inefficient compared to the direct recombination of an exciton with zero center-of-mass momentum. However, in indirect semiconductors, where the momentum indirect exciton is located below the bright state, the low temperature emission can exhibit strong phonon-assisted signals due to a large population of dark states, as depicted in Fig. 1. Recently, the PL spectrum of hexagonal boron nitride has been shown to be dominated by phonon assisted processes Cassabois et al. (2016) resulting from an indirect bandgap. Several theoretical and experimental studies have demonstrated that in tungsten-based monolayer materials intervalley excitons are located below the optically bright exciton Zhang et al. (2015); Echeverry et al. (2016); Selig et al. (2016); Berghäuser et al. (2017) and recent experimental PL studies on hBN-encapsulated tungsten diselenide have revealed a multitude of low-temperature emission peaks whose microscopic origin still needs to be clarified Courtade et al. (2017); Lindlau et al. (2017); Ye et al. (2018a); Barbone et al. (2018).

While the radiative decay of bright excitons has been extensively studied Singh et al. (2015); Palummo et al. (2015); Steinhoff et al. (2015); Pöllmann et al. (2015); Robert et al. (2016); Wang et al. (2016), the theoretical description of simultaneous exciton, phonon and photon interactions remains challenging because of the non-markovian nature of these processes. In previous studies, specific cases, such as optical phonon replicas of bright excitons, have been theoretically investigated within the polaron-frame Feldtmann et al. (2009, 2010) and non-markovian treatments of the full density matrix Chernikov et al. (2012); Shree et al. (2018). However, there is so far no general theoretical framework for the phonon-assisted exciton recombination including both optical and acoustic phonons as well as intra- and intervalley recombination channels on a microscopic footing. Moreover, the underlying dephasing rates shaping the phonon-assisted emission signals have not been derived in literature yet. Here, we present a fully microscopic model of the phonon-assisted PL from dark excitons. Based on the fundamental equations-of-motion of the many-particle density matrix we find an analytical formula which allows to calculate temperature-dependent photoluminescence spectra for thermalized exciton gases. We apply the derived formula to calculate the luminescence spectrum of hBN-encapsulated TMD monolayers. For tungsten diselenide (WSe2) we are able to explain so far observed but unidentified low temperature emission features between $50$ and $80\,\mathrm{meV}$ below the bright exciton Courtade et al. (2017); Lindlau et al. (2017); Ye et al. (2018b), by assigning them to the phonon assisted recombination of momentum-dark excitons. In contrast, for molybdenum diselenide (MoSe2) we find no additional peaks and an opposite asymmetric broadening of the bright exciton resonance compared to WSe2, which can be explained by the absence of lower lying dark excitons. Overall, our work provides new insights into phonon-assisted exciton luminescence and can be applied to determine emission spectra of arbitrary semiconducting materials.

Theoretical approach: The PL signal can be derived from the many-particle density matrix of an interacting system of electrons, phonons and photons. Throughout this work we assume a low excitation density, where the fermionic substructure of excitons can be neglected, and excitons are well described as a gas of non-interacting bosons. In this regime the many-particle physics of an undoped monolayer can be described by the excitonic Hamiltonian Toyozawa (1958); Ivanov and Haug (1993); Chernikov et al. (2012); Katsch et al. (2018), which in rotating frame reads

[TABLE]

Here, ${\mathbf{k}_{\parallel}}$ denotes the component of the vector $\mathbf{k}$ parallel to the monolayer plane and the summation is performed over all appearing quantum numbers. We use annihilation (creation) operators $a^{(\dagger)}_{\mu{\mathbf{Q}_{\parallel}}}$, $b^{(\dagger)}_{\alpha{\mathbf{q}_{\parallel}}}$ and $c^{(\dagger)}_{\sigma\mathbf{k}}$ for excitons in the state $\mu$, phonons in the mode $\alpha$, and photons with the polarization $\sigma$, respectively, while the vectorial quantum numbers denote the momentum of the particle. Moreover, the corresponding dispersions for the three particle species are given by $E^{\mu}_{{\mathbf{Q}_{\parallel}}}$, $\Omega^{\alpha}_{{\mathbf{q}_{\parallel}}}$ and $\omega^{\sigma}_{\mathbf{k}}$. Here, the excitonic bandstructure $E^{\mu}_{{\mathbf{Q}_{\parallel}}}$ is decomposed into valleys around energetic extrema, where excitons can be described as free particles with effective valley masses. The exciton index $\mu$ acts as a compound index containing main, angular, spin and valley quantum numbers. The second line of Eq. 8 describes the conversion of excitons to photons and vice versa under conservation of the in-plane momentum, whose probability is determined by the exciton-photon matrix element $M^{\mu}_{\sigma\mathbf{k}}$. Finally, excitons can scatter from the state $(\mu,{\mathbf{Q}_{\parallel}})$ to $(\nu,{\mathbf{Q}_{\parallel}}+{\mathbf{q}_{\parallel}})$ by emitting or absorbing a phonon, guided by the exciton-phonon matrix element $D^{\nu\mu}_{\alpha{\mathbf{q}_{\parallel}}}$ Selig et al. (2016, 2018); Brem et al. (2018). For convenience, we will drop the index $\parallel$ in the following keeping in mind that exciton and phonon momenta are two dimensional. Moreover, we suppress the photon and phonon mode index for a better readability.

In this framework, the PL can be determined from the temporal change of the photon density $n_{\mathbf{k}}=\langle c^{\dagger}_{\mathbf{k}}c^{\phantom{\dagger}}_{\mathbf{k}}\rangle$ Kira et al. (1999). Applying the Heisenberg equation of motion we find coupled differential equations for the photon density, the polarization $\mathcal{S}^{\mu}_{\mathbf{k}}=\langle c^{\dagger}_{\mathbf{k}}a^{\phantom{\dagger}}_{\mu{\mathbf{k}_{\parallel}}}\rangle$, the phonon-assisted polarization $\mathcal{U}^{\mu,\pm}_{\mathbf{kq}}=\langle c^{\dagger}_{\mathbf{k}}b^{(\dagger)}_{\mp\mathbf{q}}a^{\phantom{\dagger}}_{\mu{\mathbf{k}_{\parallel}}-\mathbf{q}}\rangle$ and the exciton-phonon correlation $\mathcal{C}^{\nu\mu,\pm}_{\mathbf{kq}}=\langle a^{\dagger}_{\nu{\mathbf{k}_{\parallel}}}a^{\phantom{\dagger}}_{\mu{\mathbf{k}_{\parallel}}-\mathbf{q}}b^{(\dagger)}_{\mp\mathbf{q}}\rangle$:

[TABLE]

where $\eta^{\pm}_{\mathbf{q}}=1/2\mp 1/2+\langle b^{\dagger}_{\mathbf{q}}b^{\phantom{\dagger}}_{\mathbf{q}}\rangle$ denotes the relevant phonon occupation factor for absorption/emission, $N^{\mu}_{\mathbf{Q}}=\langle a^{\dagger}_{\mu\mathbf{Q}}a^{\phantom{\dagger}}_{\mu\mathbf{Q}}\rangle$ represents the exciton occupation and $Q^{\nu\mu,\pm}_{\mathbf{kq}}=\eta^{\mp}_{\mathbf{q}}N^{\mu}_{{\mathbf{k}_{\parallel}}-\mathbf{q}}-\eta^{\pm}_{\mathbf{q}}N^{\nu}_{{\mathbf{k}_{\parallel}}}$ is the source of the corresponding exciton-phonon correlation.

To obtain the above set of equations we have factorized appearing many-particle expectation values according to the cluster expansion scheme Kira and Koch (2006). Here, we neglect coherent quantities and non-linear density dependencies, assuming a fast decoherence after the initial laser excitation and low exciton densities. Moreover, we disregard contributions connected to multi-phonon processes, cf. supplementary. Important higher order correlations have been abbreviated with $(\cdot)\rvert^{\text{phon}}_{\text{deph.}}$. In the supplementary material we show that these correlations give rise to a center-of-mass dependent phonon dephasing, which in the Born-Markov approximation reads $\Gamma^{\mu}_{\mathbf{Q}}=\pi\sum_{\nu\mathbf{q},\pm}|D^{\mu\nu}_{\mathbf{q}}|^{2}\eta^{\pm}_{\mathbf{q}}\delta(E^{\mu}_{\mathbf{Q}}\pm\Omega_{\mathbf{q}}-E^{\nu}_{\mathbf{Q+q}}).$ Moreover, additional correlations $(\cdot)\rvert^{\text{rad}}_{\text{deph.}}$ induced by the exciton-photon interaction yield a radiative dephasing for exciton states within the light cone $\gamma^{\mu}_{\mathbf{{\mathbf{k}_{\parallel}}}}=\pi\sum_{\mathbf{k}^{\prime}}|M^{\mu}_{\mathbf{k}^{\prime}}|^{2}\delta_{{\mathbf{k^{\prime}}_{\parallel}}{\mathbf{k}_{\parallel}}}\delta(E^{\mu}_{\mathbf{k}_{\parallel}}-\omega_{\mathbf{k}^{\prime}}).$

After the initial thermalization of hot excitons, the occupation numbers $\eta$ and $N$ change slowly compared to the dephasing times and the derived set of equations can be solved in the adiabatic limit (cf. supplementary material). It is important to note that a straightforward numerical evaluation of the obtained equations can yield unphysical negative PL signals, when dephasing and density scattering contributions are treated in different approximations. It is crucial to prevent double counting of the appearing processes as dephasing and as phonon-assisted decay. A detailed discussion is provided in the supplementary material. For a thermalized exciton distribution, we find the following analytic expression for the $\sigma$-polarized photon flux emitted in perpendicular direction with respect to the monolayer:

[TABLE]

The first term is analogue to the markovian result for the direct recombination, with the small difference, that the phonon dephasing here only appears in the denominator. However, for energies $\omega\approx E_{0}$ the second term containing the scattering contributions, can be rewritten as the phonon-guided in-scattering to the bright state, which in thermal equilibrium is equal to the out-scattering $2\Gamma^{\mu}_{\mathbf{0}}N^{\mu}_{\mathbf{0}}$. Therefore, in the resonant case, Eq. 26 can be well approximated with the Elliot formula Koch et al. (2006)

[TABLE]

The second part in Eq. 26 containing phonon-assisted decay channels agrees with luminescence formulas derived in previous studies Feldtmann et al. (2009, 2010) for optical phonon replicas in the polaron picture. However, it is here generalized to arbitrary phonon modes, intervalley scattering and a consistent description of the phonon-induced dephasing of direct and indirect transitions. Note, that the phonon-assisted recombination corresponds to additional luminescence decay channels involving exciton densities of momentum dark states. In contrast, in absorption experiments, phonon-sidebands arise from a frequency dependent dephasing and energy renormalization, which asymmetrically shape the optical response function of the bright exciton Christiansen et al. (2017). The above presented procedure can in principle also be applied to obtain transition rates for other crossed interactions, such as defect-scattering assisted tunnelling.

Phonon-assisted photoluminescence in TMDs: Now, we apply Eq. 26 to study the luminescence of TMD monolayers encapsulated with hBN. To this end, we use ab initio parameters from literature for the electronic bandstructure, phonon dispersion, dielectric constants and electron-phonon coupling Kormányos et al. (2015); Jin et al. (2014); Laturia et al. (2018). The excitonic properties of the system are derived in effective mass approximation by solving the Wannier equation Koch et al. (2006). Details about the used parameters, the applied screening model and the form of the excitonic matrix elements are given in the supplementary material.

Figure 2 shows calculated PL spectra for hBN-encapsulated monolayer WSe2 at four different temperatures. In addition to the full spectra calculated from Eq. 26 (color shaded), we also present the spectra in the case of $\Gamma^{\nu}_{\mathbf{q}}\longrightarrow 0$, corresponding to the Fermi’s golden rule (solid line). The sharp steps in the solid line resulting from the strict energy conservation in the Fermi’s golden rule conveniently illustrate the multitude of exciton valleys and phonon modes contributing to the overall PL signal. Here, each step in the spectrum corresponds to a transition from the energetic minimum of a certain valley while absorbing/emitting a phonon.

At 300K and 150K the PL spectrum is dominated by the bright exciton resonance, but already here phonon-assisted indirect recombinations have a significant impact resulting in an asymmetric resonance. At energies larger than $E_{0}$ the signal is mainly shaped by recombinations of small momentum excitons in the K-K valley, which scatter to virtual states in the light cone via low energy acoustic phonons. For energies smaller then $E_{0}$, much larger broadening is observed. This can be ascribed to K-K excitons emitting optical phonons and - more importantly - to intervalley transitions from the K-$\Lambda$ and the K-K’ state (hole located at the K and electron at the $\Lambda$ and the K’ valley, respectively). For the considered dielectric environment our Wannier model predicts that the latter are located about $34\,\mathrm{meV}$ and $46\,\mathrm{meV}$ below the bright exciton, respectively.

At low temperatures, the optical response shifts towards multiple indirect PL peaks below the bright exciton stemming from the phonon-assisted recombination of the energetically lowest K-K’ exciton. Although scattering to virtual states becomes more improbable with increasing detuning from the bright state, the large population ratio between dark and bright states gives rise to a strong indirect signal at low temperatures (Fig. 1). The predicted low temperature PL signals at about $-60\,\mathrm{meV}$ (K-K’ - acoustic phonon assisted) and $-75\,\mathrm{meV}$ (K-K’ - optical phonon assisted) correspond well to experimentally observed PL peaks below the trion resonance. In agreement with our theory, these peaks are not visible in reflection/absorption spectra Christiansen et al. (2017). Although clearly visible in PL spectra, so far these peaks were to a large extent ignored in literature due to their unclear origin. Our work reveals that these peaks stem from indirect phonon-assisted transitions from lower lying momentum-dark excitons in tungsten-based TMDs. We have also performed calculations for WS2 (not shown) obtaining similar low-temperature PL features.

In Fig. 3 we directly compare an experimentally measured PL spectrum with our simulation. The blue shaded curve shows a spectrum measured at T=15 K for hBN-encapsulated WSe2 at charge neutrality. Details about the experiment can be found in Ref. Courtade et al. (2017). The solid and dashed lines show calculated spectra for 15K phonon temperature and three different effective exciton temperatures. For a better comparisson the calculated spectra where convoluted with a 1 meV broad Gaussian to mimic weak disorder. Note that the calculated spectra in Fig. 2 correspond to stationary PL signals of exciton distributions which are in thermal equilibrium with the lattice. In Fig. 3 we introduce different phonon and exciton temperatures to take into account that in experiments performed at low temperatures the optically injected hot excitons decay before they have reached thermal equilibrium with the lattice. Figure 3 illustrates that the theoretically predicted luminescence features from dark excitons agree well with the experimentally observed resonances between 45 and 80meV below the bright exciton denoted with L, P and Q. The two pronounced peaks around 40 and 30 meV below the bright state (denoted with D and T) potentially stem from spin-forbidden states and residual trions Ye et al. (2018a), which are not included in our model. We find the best agreement between experiment and theory for an effective exciton temperature of about 80K, indicating that the luminescence stems from an exciton distribution that is much hoter than the lattice. In the supplementary material we have summarized several independent PL measurements on hBN-encapsulated WSe2 Courtade et al. (2017); Ye et al. (2018a); Barbone et al. (2018) in which the above discussed phonon-assisted peaks have been observed. It is important to note, that the peak denoted as P in Fig. 3, which we attribute to the acoustic phonon assisted decay of K-K’ excitons, in experiments performed at 4 K shows a double peak structure. In the supplemetary, we have included a discussion about this splitting.

Finally, we compare our results obtained for WSe2 (representing tungsten-based 2D materials) with molybdenum-based TMDs, cf. Fig. 4. Here, we show a continuous temperature study of the PL spectrum for (a) WSe2 and (b) MoSe2. In tungsten-based TMDs the lower lying K-$\Lambda$ and K-K’ valley provide a continuous density of (momentum-dark) states below the bright exciton, yielding an asymmetric broadening of the bright exciton peak towards lower energies. In contrast, for MoSe2 the bright state constitutes the global minimum of the exciton dispersion calculated in our Wannier model. Therefore there are no phonon-assisted PL peaks below the bright state, and the main peak is asymmetrically broadened towards higher energies, cf. Fig. 4(b). Here the broadening of the low energy side of the main line is predominantly given by the energy uncertainty of the bright state giving rise to a Lorentzian shape. In contrast, the high energy side is additionally shaped by the decay of small momentum excitons in the K-K valley assisted by long range acoustic phonons, which agrees with the findings in Ref. Shree et al., 2018.

In summary, we have presented an analytical expression for the phonon-assisted exciton photoluminescence from momentum-dark excitons which allows to model the emission spectrum of an arbitrary semiconducting material. When applying our model to the luminescence of WSe2, we find in good agreement with experiment clear PL signals stemming from dark intervalley excitons. This explains the origin of so far observed but unidentified PL peaks. In contrast, MoSe2 does not exhibit indirect PL peaks and shows an opposite asymmetric broadening which is consistent with the absence of lower lying dark exciton states. Our work will trigger future experimental studies on dark exciton PL and in particular time-resolved PL measurements addressing exciton thermalization.

Appendix A Equations of Motion

In the low excitation regime the fermionic substructure of excitons as well as the exciton-exciton interaction can be neglected, so that the properties of the TMD monolayer are described by the following excitonic hamiltonian in rotating frame:

[TABLE]

For definitions see main text. For convinience we will in the following supress the photon, phonon and exciton mode index. The derivation shown below can be performed analogously for multiple exciton states and photon/phonon modes, yielding and additive solution. To further reduce the amount of indices in the equations we consider the case of perpendicular emission to the monolayer plane, removing the in-plane momentum transfer due to interaction with photons, whose consideration is however straight forward. Applying the Heisenberg equation of motion $i\hbar\partial_{t}\langle\mathcal{\cdot}\rangle=\langle[\cdot,H]\rangle$ to calculate the evolution of the photon density $n_{\mathbf{k}}=\langle c^{\dagger}_{\mathbf{k}}c^{\phantom{\dagger}}_{\mathbf{k}}\rangle$ we find a coupled system of equations involving the polarization $\mathcal{S}_{\mathbf{k}}=\langle c^{\dagger}_{\mathbf{k}}a^{\phantom{\dagger}}_{0}\rangle$, the phonon-assisted polarization $\mathcal{U}^{\pm}_{\mathbf{kq}}=\langle c^{\dagger}_{\mathbf{k}}b^{(\dagger)}_{\mp\mathbf{q}}a^{\phantom{\dagger}}_{-\mathbf{q}}\rangle$ (with phonon creation for the $+$ case), the exciton-phonon correlation $\mathcal{C}^{\pm}_{\mathbf{Qq}}=\langle a^{\dagger}_{\mathbf{Q}}a^{\phantom{\dagger}}_{\mathbf{Q}-\mathbf{q}}b^{(\dagger)}_{\mp\mathbf{q}}\rangle$ and the exciton density $N_{\mathbf{Q}}=\langle a^{\dagger}_{\mathbf{Q}}a^{\phantom{\dagger}}_{\mathbf{Q}}\rangle$:

[TABLE]

Here $\eta^{\pm}_{\mathbf{q}}=1/2\mp 1/2+\langle b^{\dagger}_{\mathbf{q}}b^{\phantom{\dagger}}_{\mathbf{q}}\rangle$ denotes the relevant phonon occupation factor for absorbtion/emission, which are assumed to be time independent (bath approximation). Moreover, $Q^{\pm}_{\mathbf{Q}\mathbf{q}}=\eta^{\mp}_{\mathbf{q}}N_{\mathbf{Q}-\mathbf{q}}-\eta^{\pm}_{\mathbf{q}}N_{\mathbf{Q}}$ is the source of the corresponding exciton-phonon correlation and $n_{\mathbf{kk^{\prime}}}=\langle c^{\dagger}_{\mathbf{k}}c^{\phantom{\dagger}}_{\mathbf{k^{\prime}}}\rangle$ stands for the photonic density matrix.

To obtain the above set of equations we have factorized appearing many-particle expectation values acording to the cluster expansion scheme [cite], viz. for example

[TABLE]

Which means that we neglect coherent quantities, e.g. $\langle a\rangle$, assuming a fast decoherence after the initial laser excitation, and discard pure photon-phonon correlations. Moreover, we have introduced the notation $\langle\cdot\rangle^{\text{corr.}}$, which in the upper case quantifies the correlation between phonon density and polarization. These higher order correlations are abbreviated in Eq. 11 and 12 via $X^{\zeta\zeta^{\prime}}_{\mathbf{kqq^{\prime}}}=\langle c^{\dagger}_{\mathbf{k}}b^{(\dagger)}_{-\zeta\mathbf{q}}b^{(\dagger)}_{-\zeta^{\prime}\mathbf{q}^{\prime}}a^{\phantom{\dagger}}_{-\mathbf{q}-\mathbf{q^{\prime}}}\rangle^{\text{corr.}}$ and $Y^{\zeta\zeta^{\prime}}_{\mathbf{Qqq^{\prime}}}=\langle b^{(\dagger)}_{-\zeta\mathbf{q}}b^{(\dagger)}_{-\zeta^{\prime}\mathbf{q}^{\prime}}a^{\dagger}_{\mathbf{Q}}a^{\phantom{\dagger}}_{\mathbf{Q}-\mathbf{q}-\mathbf{q^{\prime}}}\rangle^{\text{corr.}}$, while the case $\zeta=+1$ again corresponds to the phonon creation. Similar as the coupling of $\mathcal{S}$ to $\mathcal{U}$ gives rise to a phonon induced energy renormalization as well as additional phonon assisted transitions, the coupling of $\mathcal{U}$ to $X$ introduces two-phonon-assisted transitions as well as an energy renormalization for $\mathcal{U}$. The microscopic treatment of higher order phonon-emission sidebands goes beyond the scope of this work, however, we will include the energy renormalizations induced by the 2-phonon-correlations.

Appendix B Many-Particle Dephasing

To obtain the many-particle dephasing induced by electron-phonon interaction, we additionally consider the equations of motion of $X$ and $Y$. In the following we only show the derivation for $X$, since $Y$ can be treated in complete analogy. Applying the same approximations as above we find:

[TABLE]

Additional source terms $\dot{X}\rvert_{\text{source}}$ following from the exciton-photon interaction give rise to two-phonon-assisted photoemission processes, which will not be considered here. Moreover, higher order exciton-phonon correlations $\dot{X}\rvert_{\text{corr.}}$ give rise to equivalent energy renormalizations as $X$ induces for $\mathcal{U}$, denoted with $\delta^{+}$ which will be taken into accoutn below via a self-consistent renormalization. In the static limit we obtain,

[TABLE]

where in the second step we assume that $\mathcal{U}^{\zeta}_{\mathbf{kq}}$ only gives a non-zero response for photon energies close to its resonance, viz. $\omega_{\mathbf{k}}\approx E_{\mathbf{q}}-\zeta\Omega_{\mathbf{q}}$ which is analogous to the treatment of scattering contributions in second order Born-Markov-approximation [cite]. Plugging Eq. 16 into 11 , we find that the higher order exciton-phonon correlations give rise to the polaronic energy renormalization

[TABLE]

Where $\delta^{+}$ is assumed to yield the same energy renormalization for $X$ as for all other orders, giving rise to the self-consistent polaron self-energy $\Sigma^{\text{ph}}$. The same treatment of $Y$ yields an analogous result in Eq. 12. Moreover, the coupling of $\mathcal{S}$ to $n_{\mathbf{kk^{\prime}}}$ and the backcoupling of $\mathcal{C}$ to $\mathcal{U}$ can be similarly shown to yield a radiative dephasing for exciton states in the lightcone,

[TABLE]

where we assume short photon lifetimes $1/\kappa$ and a weakly varying coupling element $M_{\mathbf{k}}$ in vicinity of $E_{0}$.

Appendix C Photoluminescence in Static Limit

Using the derived renormalizations, $E_{\mathbf{Q}}\longrightarrow\tilde{E}_{\mathbf{Q}}=E_{\mathbf{Q}}+\Sigma^{\text{ph}}_{\mathbf{Q}}-i\gamma_{\mathbf{0}}\delta_{\mathbf{Q,0}}$, we now solve Eq. 9 - 12 by assuming slowly varying densities allowing to find adiabatic solutions from the static limit of our equations of motion. Introducing the Greensfunctions

[TABLE]

we find:

[TABLE]

Now we can identify,

[TABLE]

where we consistent with Eq. 16 evaluate the selfenergy at resonance. Next we recast

[TABLE]

Hence, pluging Eq.20 in 9 yields,

[TABLE]

When considering the equation of motion for the exciton density in the lightcone, the phonon-scattering contribution in the static limit can be writen as

[TABLE]

resembling the semi-classical Boltzman scattering-equations. This temporal change of the exciton occupation exactly corresponds to the second contribution in the second line of Eq. 23. Since we have assumed slowly varying occupations we have to neglect this contribution to stay consisent with our adiabatic solution.

Finally, the contribution proportional to $N_{0}$ of the remaining term in the second line of Eq. 23 can again acording to Eq. 21 be rewritten as $-|\mathcal{G}_{0}(\omega_{\mathbf{k}})|^{2}\Im\text{m}\{\Sigma^{\text{ph}}(\omega_{\mathbf{k}})\}N_{0}$, which cancels the imaginary part of the polaron self-energy in the numerator of the first line in Eq. 23. Throughout this work we neglect polaron shifts and only consider the imaginary part of the self-energy (dephasing),

[TABLE]

Hence, we obtain

[TABLE]

The generalization to several exciton,phonon and photon modes is straight forward and can be found in the main text.

Appendix D Excitonic Matrix Elements and Wavefunctions

To evaluate the above derived PL formula we make use of the Mott-Wannier model of excitons. In this framework the exciton-phonon matrix element is determined by the electron-phonon couplings for electrons and holes $g^{e/h}$ and the exciton wavefunction in momentum space $\Phi(\mathbf{k})$ via Selig et al. (2018); Brem et al. (2018)

[TABLE]

Here the momentum transferred to the electron (hole) is denoted by $q_{e(h)}$, when the exciton gains a center-of-mass momentum $q=q_{e}+q_{h}$. The exciton index $\nu$ here acts as a compound index containing principal quantum number, angular momentum, electron/hole valley and spin configuration. For the carrier-phonon coupling $g$ we use the deformation potential approximations for acoustic and optical phonons deduced from density functional perturbation theory (DFPT) in Ref. Jin et al. (2014). Furthermore, the radiative dephasing is determined via

[TABLE]

where the interband momentum matrix element $\mathbf{m}^{vc}$ is derived from a two band $k\cdot p$ Hamiltonian, which in vicinity of the K point yields Xiao et al. (2012)

[TABLE]

The next neighbor hopping integral $t=\hbar/a_{0}\sqrt{E_{g}/(m_{e}+m_{h})}$ is determined by the effective masses $m_{e/h}$ of electrons and holes and the single particle bandgap $E_{g}$ at the K-point, while $\sigma=\pm 1$ for left-(right-)handed circularly polarized light. To calculate the exciton binding energies in a monolayer we use an approach analog to the Keldish screening for charges in a thin film of thickness d surrounded by a dielectric environment. We explicitly take into account anisotropic dielectric tensors. Solving the Poisson equation for the boundary conditions of an encapsulated monolayer yields $W_{q}=V_{q}/\epsilon_{scr}(q)$, with the bare 2D-Fourier transformed Coulomb potental $V_{q}$ and the non-local screening,

[TABLE]

where $\kappa_{i}=\sqrt{\epsilon^{\parallel}_{i}\epsilon^{\bot}_{i}}$ and $\alpha_{i}=\sqrt{\epsilon^{\parallel}_{i}/\epsilon^{\bot}_{i}}$ account for the parallel and perpendicular component of the dielectric tensor $\epsilon$ of the monolayer ($i=1$) and the environment ($i=2$), which can be found in refs. Laturia et al. (2018); Geick et al. (1966).

Finally, to calculate the excitonic wavefunctions we numerically solve the the Wannier equation,

[TABLE]

Here we use relative (k) and center-of mass coordinates (Q) with $\alpha(\beta)=m_{c(v)}/(m_{c}+m_{v})$ assuming effective electron ($m_{c}$) and hole ($m_{v}$) masses. Within the vicinity of minima and maxima of valence and conduction band, we approximate the dispersions quadratically, which allows us to separate relative and center of mass motion. When $\mathbf{K}_{c}$ denotes the conduction band valley and $\mathbf{K}_{v}$ the valence band valley, we find $\Phi_{\nu\mathbf{Q}}(\mathbf{k})=\Phi_{\nu}(\mathbf{k})=\Psi_{\nu}(\mathbf{k}-\alpha\mathbf{K}_{v}-\beta\mathbf{K}_{c})$, with $\Psi$ obeying the effective electron-hole Schroedinger equation,

[TABLE]

where $m_{\text{r}}=(m_{c}m_{v})/(m_{c}+m_{v})$ is the reduced exciton mass for the corresponding valley. Furthermore, the parabolic approximation yields the center-of-mass dispersion $E_{\nu\mathbf{Q}}=E^{\text{bind}}_{\nu}+\hbar^{2}(\mathbf{Q}-[\mathbf{K_{c}}-\mathbf{K_{v}}])^{2}/(2[m_{c}+m_{v}])+\varepsilon_{c\mathbf{K_{c}}}-\varepsilon_{v\mathbf{K_{v}}}$. Note, that exciton wavefunctions with different valley configurations are centered at different momenta. All necessary electronic band parameters can be found in Kormányos et al. (2015).

Appendix E Experimental Observations

Recent experimental PL studies on hBN-encapsulated tungsten diselenide have revealed a multitude of low-temperature emission peaks whose microscopic origin has not been fully clarified yet Courtade et al. (2017); Lindlau et al. (2017); Ye et al. (2018a); Barbone et al. (2018). Most studies have so far focused on the impact of bound exciton configurations, such as trions, biexcitons and trapped excitons, while the potential influence of indirect phonon-assisted recombination of intrinsically dark exciton states has been ignored to a large extend. In Fig. 5 we have summarized the experimental observations of three independent measurements of the cryogenic luminescence from hBN-encapsulated WSe2.

Fig. 5a) shows PL spectra of different samples at charge neutrality. All three measured samples denoted with S1, S2 and H7 (same sample as used in main text) consistently exhibit several peaks below the bright exciton. In particular, the theoretically predicted signal denoted with P in the main text, which we attribute to the acoustic phonon assisted recombination of K-K’ excitons, is systematically observed in form of two peaks denoted with P1 and P2 in Fig. 5a)-c). The exact shape and intensity of the peaks P1 and P2 depends on excitation conditions (laser energy and power density) and temperature and we observe slight variations even for the same sample during different experiments. The same qualitative behaviour for the peaks P1 and P2 is observed in the measurements on very similar samples shown in Fig. 5b) and d) which have been performed by other groups namely Ye et.al. Ye et al. (2018a) and Barbone et.al. Barbone et al. (2018), respectively. Note, that Fig. 5b) presents a continuous study of the impact of the applied gate voltage, whereas d) is measured without gate voltage. Although the shape and relative intensity of the two peaks P1 and P2 seems to vary on different samples, they are clearly visible in all three measurements at a consistent energetic position, which agrees with the theoretically predicted peak stemming from acoustic phonon-assisted recombination of K-K’ excitons, cf. main text.

However, the signal in all three measurements appears as two separate peaks with a splitting of about 5 meV. Furthermore, it is striking that the energetic distance between P1 and P2 to the peak denoted as “T” in Fig. 5a) corresponds approximately to the energy of KTO and KLO phonon, respectively. This suggests that these peaks could be phonon side bands of the “T” peak whose origin is still under debate. However, Fig. 5c) shows a polarization-resolved PL spectrum at the charge neutrality point measured by Ye et.al. (cf. supplementary of ref.Ye et al., 2018a). Here a clear polarization dependence of those two peaks is visble, i.e. in cross-polarized case, only a single peak is observed. This polarization dependence indicates that the two peaks do not stem from the same exciton state, but might rather be a result of e.g. Coulomb exchange-mediated splitting of the two spin configurations of momentum-dark excitons. In particular, a mixing of momentum-dark states and direct K-K excitons could give rise to a Coulomb-exchange induced splitting of the K-K’ and K’-K exciton, similar as in the case of spin-forbidden dark states Robert et al. (2017). In summary, further research is needed to understand the origin of the “T” peak and the possible impact of Coulomb exchange coupling on the splitting and hybridization of bright and dark exciton states.