Excitons versus electron-hole plasma in monolayer transition metal dichalcogenide semiconductors

TL;DR

This paper investigates the conditions under which excitons or electron-hole plasma dominate in monolayer transition metal dichalcogenides, revealing a tunable Mott transition influenced by dielectric environment and doping.

Contribution

It introduces a comprehensive many-body theoretical framework to predict and control the exciton-to-plasma transition in 2D semiconductors, including the Mott transition point.

Findings

Exciton-to-plasma ratio can be tuned by dielectric screening and doping.

A Mott transition occurs at specific excitation densities.

Predicted observable effects in photoemission and tunneling microscopy.

Abstract

When electron-hole pairs are excited in a semiconductor, it is a priori not clear if they form a fermionic plasma of unbound particles or a bosonic exciton gas. Usually, the exciton phase is associated with low temperatures. In atomically thin transition metal dichalcogenide semiconductors, excitons are particularly important even at room temperature due to strong Coulomb interaction and a large exciton density of states. Using state-of-the-art many-body theory including dynamical screening, we show that the exciton-to-plasma ratio can be efficiently tuned by dielectric substrate screening as well as charge carrier doping. Moreover, we predict a Mott transition from the exciton-dominated regime to a fully ionized electron-hole plasma at excitation densities between $3\times10^{12}$ cm$^{-2}$ and $1\times10^{13}$ cm$^{-2}$ depending on temperature, carrier doping and dielectric environment. We propose the observation of these effects by studying excitonic satellites in photoemission spectroscopy and scanning tunneling microscopy.

Tables

Table 1. Table 1: Parameters describing the Coulomb interaction and the static screening in the semiconducting TMDCs M ​ X 2 𝑀 subscript 𝑋 2 MX_{2} .

	MoS2	MoSe2	WS2	WSe2
lattice constant $a$(Å)	3.180	3.320	3.191	3.325
bare Interaction $U$
$\gamma$ (Å)	1.932	2.232	2.130	2.297
$\delta$ (Å2)	0.395	-0.356	0.720	0.174
A (Å2)	8.758	9.546	8.818	9.574
$U_2$ (eV)	0.810	0.837	0.712	0.715
$U_3$ (eV)	0.367	0.376	0.354	0.360
screening $\epsilon$
$a$ (1/Å2)	2.383	2.856	3.947	2.430
$b$	17.836	11.635	29.931	20.764
$c$ (Å)	5.107	1.979	5.440	5.761
$h$ (Å)	2.740	4.298	1.578	2.489
$e$	5.739	6.303	4.497	5.305
${\epsilon }_2$	3.077	3.148	2.979	3.028
${\epsilon }_3$	2.509	2.510	2.494	2.481

Equations

$$n_a({\mu }_a,T)=\frac{i\hslash }{\mathcal{A}}{\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }\sum _{\mathbf{k}\sigma }{f}^a(\omega )A_{\mathbf{k}\sigma }^a(\omega ).$$

$$G_{\mathbf{k}\sigma }^{\text{ret},a}(\omega )=\frac{1}{\hslash \omega -{\epsilon }_{\mathbf{k}\sigma }^a-{\Sigma }_{\mathbf{k}\sigma }^{\text{ret},a}(\omega )}.$$

$$\begin{array}{r}{n}_a=n_{\text{free}}^a+n_X,\end{array}$$

$${\alpha }_a=\frac{n_{\text{free}}^a}{n_a}$$

$${\mu }_{X,\nu }={\mu }_e+{\mu }_h-E_{\nu },$$

$$\frac{n_X}{n_e{n}_h}\propto e^{\beta I_{\nu }^{\text{eff}}}.$$

$$\begin{array}{r}S=N{k}_B\left[\text{ln}\left(\frac{V}{N}c(T)\right)+\frac{5}{2}\right],\end{array}$$

$$A_{\mathbf{k}\sigma }^a(\omega )=2i\text{Im}\frac{1}{\hslash \omega -{\epsilon }_{\mathbf{k}\sigma }^a-{\Sigma }_{\mathbf{k}\sigma }^{\text{ret},a}(\omega )}.$$

$$\begin{array}{rl}& n_a({\mu }_a,T)\\ & =\frac{1}{\mathcal{A}}\sum _{\mathbf{k}\sigma }{f}^a(E_{\mathbf{k}\sigma }^a)\\ & -\frac{1}{\mathcal{A}}\sum _{\mathbf{k}\sigma }{\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }\frac{2}{\hslash }\text{Im}{\Sigma }_{\mathbf{k}\sigma }^{\text{ret},a}(\omega )\left[f^a(E_{\mathbf{k}\sigma }^a)-f^a(\omega )\right]\\ & \frac{d}{d\omega }\frac{\mathcal{P}}{\omega -E_{\mathbf{k}\sigma }^a/\hslash }\\ & =n_a^{\text{QP}}+n_a^{\text{corr}},\end{array}$$

$$\begin{array}{rl}{A}_{\mathbf{k}\sigma }^a(\omega )=& -2\pi i\delta (\hslash \omega -E_{\mathbf{k}\sigma }^a)(1-Z_{\mathbf{k}\sigma }^a)\\ & -2\pi i{\Gamma }_{\mathbf{k}\sigma }^a(\omega ),\end{array}$$

$$\begin{array}{rl}& \text{Im}{\Sigma }_{\mathbf{k}\sigma }^{\text{ret},a}(\omega )\\ & =\frac{1}{\mathcal{A}}\sum _{\mathbf{k}b}{V}_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{k}{\mathbf{k}}^{\prime }}^{ab}\text{Im}{\epsilon }_{\mathbf{k}-{\mathbf{k}}^{\prime }}^{-1,\text{ret}}(\omega -E_{{\mathbf{k}}^{\prime }\sigma }^b/\hslash )\\ & \left[f^b(E_{{\mathbf{k}}^{\prime }\sigma }^b)-n^B(E_{{\mathbf{k}}^{\prime }\sigma }^b-\hslash \omega )\right]\\ & +\frac{1}{\mathcal{A}}\sum _{\mathbf{k}b{\sigma }^{\prime }}\text{Im}{T}_{\mathbf{k}{\mathbf{k}}^{\prime }\sigma {\sigma }^{\prime }}^{\prime \prime ,\text{ret},ab}(\omega +E_{{\mathbf{k}}^{\prime }{\sigma }^{\prime }}^b/\hslash )\\ & \left[f^b(E_{{\mathbf{k}}^{\prime }{\sigma }^{\prime }}^b)+n_{ab}^B(\hslash \omega +E_{{\mathbf{k}}^{\prime }{\sigma }^{\prime }}^b)\right].\end{array}$$

$$\begin{array}{rl}& V_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{k}{\mathbf{k}}^{\prime }}^{\text{S},<,ab}(\omega )=n^B(\omega )V_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{k}{\mathbf{k}}^{\prime }}^{ab}2i\text{Im}{\epsilon }_{\mathbf{k}-{\mathbf{k}}^{\prime }}^{-1,\text{ret}}(\omega ),\\ & V_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{k}{\mathbf{k}}^{\prime }}^{\text{S},>,ab}(\omega )=(1+n^B(\omega ))V_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{k}{\mathbf{k}}^{\prime }}^{ab}2i\text{Im}{\epsilon }_{\mathbf{k}-{\mathbf{k}}^{\prime }}^{-1,\text{ret}}(\omega ).\end{array}$$

$$T^{\text{ret},ab}(\omega )=V^{\text{S},\text{ret},ab}+i\hslash V^{\text{S},\text{ret},ab}{\mathcal{G}}^{\text{ret},ab}(\omega )T^{\text{ret},ab}(\omega ),$$

$$\begin{array}{rl}& 0=(E_{\mathbf{k}\sigma }^a+E_{\mathbf{k}+\mathbf{Q}{\sigma }^{\prime }}^b)G_{2,\mathbf{k},\mathbf{k}+\mathbf{Q}\sigma {\sigma }^{\prime }}^{\text{ret},ab}\\ -& (1-f_{\mathbf{k}\sigma }^a-f_{\mathbf{k}+\mathbf{Q}{\sigma }^{\prime }}^b)\frac{1}{\mathcal{A}}\sum _{{\mathbf{k}}^{\prime }}{V}_{\mathbf{k}+\mathbf{Q},{\mathbf{k}}^{\prime },\mathbf{k},{\mathbf{k}}^{\prime }+\mathbf{Q}}^{\text{S},ab}{G}_{2,{\mathbf{k}}^{\prime },{\mathbf{k}}^{\prime }+\mathbf{Q}\sigma {\sigma }^{\prime }}^{\text{ret},ab}.\end{array}$$

$$\begin{split}&\psi^{\sigma\sigma^{\prime}}_{\nu\mathbf{Q}}(\mathbf{k})=\big{<}\mathbf{k}\mathbf{k}^{\prime}\sigma\sigma^{\prime}ab\big{|}\nu\sigma\sigma^{\prime}\mathbf{Q}\big{>}\delta_{\mathbf{k}^{\prime},\mathbf{k}+\mathbf{Q}}\,,\end{split}$$

$$\begin{array}{rl}& T^{\text{ret},ab}(\omega )=\\ & \sum _{\nu \sigma {\sigma }^{\prime }\mathbf{Q}}{N}_{ab}^{-1}(E_{\text{kin}}-\hslash \omega )\frac{|\nu \sigma {\sigma }^{\prime }\mathbf{Q}⟩⟨\nu \widetilde{\sigma {\sigma }^{\prime }}\mathbf{Q}|}{\hslash \omega -E_{\nu \mathbf{Q}}^{\sigma {\sigma }^{\prime }}}(E_{\text{kin}}-E_{\nu \mathbf{Q}}^{\sigma {\sigma }^{\prime }})\end{array}$$

$$\begin{array}{rl}{n}_a({\mu }_a,T)& =\frac{1}{\mathcal{A}}\sum _{\mathbf{k}\sigma }{f}^a(E_{\mathbf{k}\sigma }^a)\\ & +\frac{1}{{\mathcal{A}}^2}\sum _{\mathbf{k}{\mathbf{k}}^{\prime }b{\sigma }^{\prime }}{\int }_{-\infty }^{E_{\text{Gap}}}\frac{d\omega }{\pi }{n}_{ab}^B(\omega )\\ & \text{Im}{T}_{\mathbf{k}{\mathbf{k}}^{\prime }\sigma {\sigma }^{\prime }}^{\text{ret},ab}(\omega )\frac{d}{d\omega }i\hslash {\mathcal{G}}_{\mathbf{k}{\mathbf{k}}^{\prime }\sigma {\sigma }^{\prime }}^{\text{ret},ab}(\omega )\\ & +n_{\text{scatt}}.\end{array}$$

$$\begin{split}&\frac{1}{\mathcal{A}^{2}}\sum_{\mathbf{k}\mathbf{k}^{\prime}b\sigma^{\prime}}\textrm{Im}\,T^{\textrm{ret},ab}_{\mathbf{k}\mathbf{k}^{\prime}\sigma\sigma^{\prime}}(\omega)\frac{d}{d\omega}i\hbar\mathcal{G}^{\textrm{ret},ab}_{\mathbf{k}\mathbf{k}^{\prime}\sigma\sigma^{\prime}}(\omega)\Big{|}_{\textrm{bound}}\\ &=\pi\frac{1}{\mathcal{A}}\sum_{\nu\sigma^{\prime}\mathbf{Q}}\delta(\hbar\omega-E^{\sigma\sigma^{\prime}}_{\nu\mathbf{Q}})\,.\end{split}$$

$$\begin{array}{rl}{n}_a({\mu }_a,T)& =\frac{1}{\mathcal{A}}\sum _{\mathbf{k}\sigma }{f}^a(E_{\mathbf{k}\sigma }^a)\\ & +\frac{1}{\mathcal{A}}\sum _{b\ne a}\sum _{\sigma {\sigma }^{\prime }}\sum _{\nu \mathbf{Q}}{n}_{ab}^B(E_{\nu \mathbf{Q}}^{\sigma {\sigma }^{\prime }})\\ & =n_{\text{free}}^{\text{GW},a}+n_X.\end{array}$$

$$\begin{array}{rl}{E}_{\mathbf{k}\sigma }^a& ={\epsilon }_{\mathbf{k}\sigma }^a+\text{Re}{\Sigma }_{\mathbf{k}\sigma }^{\text{GW},\text{ret},a}(E_{\mathbf{k}\sigma }^a)+\text{Re}{\Sigma }_{\mathbf{k}\sigma }^{\text{T},\text{ret},a}(E_{\mathbf{k}\sigma }^a).\end{array}$$

$$\begin{array}{rl}& {\Gamma }_{\mathbf{k}\sigma }^a(\omega )=\frac{1}{\mathcal{A}}\sum _{b\ne a}\sum _{\nu {\sigma }^{\prime }\mathbf{Q}}\sum _{{\mathbf{k}}^{\prime }}\\ & (|{\psi }_{\nu \mathbf{Q}}^{\sigma {\sigma }^{\prime }}(\mathbf{k}){|}^2{\delta }_{a,h}{\delta }_{{\mathbf{k}}^{\prime },\mathbf{k}+\mathbf{Q}}+|{\psi }_{\nu \mathbf{Q}}^{\sigma {\sigma }^{\prime }}(\mathbf{k}-\mathbf{Q}){|}^2{\delta }_{a,e}{\delta }_{{\mathbf{k}}^{\prime },\mathbf{k}-\mathbf{Q}})\\ & \delta (\hslash \omega +E_{{\mathbf{k}}^{\prime }{\sigma }^{\prime }}^b-E_{\nu \mathbf{Q}}^{\sigma {\sigma }^{\prime }})[n_{ab}^B(E_{\nu \mathbf{Q}}^{\sigma {\sigma }^{\prime }})+f^b(E_{{\mathbf{k}}^{\prime }{\sigma }^{\prime }}^b)].\end{array}$$

$$\begin{split}\Sigma_{\mathbf{k}\sigma}^{\textrm{F},a}=&i\hbar\frac{1}{\mathcal{A}}\sum_{\mathbf{k}^{\prime}}V^{aa}_{\mathbf{k}\mathbf{k}^{\prime}\mathbf{k}\mathbf{k}^{\prime}}G^{<,a}_{\mathbf{k}^{\prime}\sigma}\\ =-&i\hbar\frac{1}{\mathcal{A}}\sum_{\mathbf{k}^{\prime}}V^{aa}_{\mathbf{k}\mathbf{k}^{\prime}\mathbf{k}\mathbf{k}^{\prime}}\int\frac{d\omega}{2\pi}f^{a}(\omega)A_{\mathbf{k}^{\prime}\sigma}^{a}(\omega)\\ =-&\frac{1}{\mathcal{A}}\sum_{\mathbf{k}^{\prime}}V^{aa}_{\mathbf{k}\mathbf{k}^{\prime}\mathbf{k}\mathbf{k}^{\prime}}\Big{(}f^{a}(E^{a}_{\mathbf{k}^{\prime}\sigma})+f^{a,\textrm{bound}}_{\mathbf{k}^{\prime}\sigma}\Big{)}\,~{}.\end{split}$$

$$\begin{split}E^{a}_{\mathbf{k}\sigma}&\approx\varepsilon_{\mathbf{k}\sigma}^{a}+\textrm{Re}\,\Sigma_{\mathbf{k}\sigma}^{\textrm{GW},\textrm{ret},a}(E^{a}_{\mathbf{k}\sigma})\Big{|}_{\textrm{free}}\\ &=\varepsilon_{\mathbf{k}\sigma}^{a}+\Sigma_{\mathbf{k}\sigma}^{\textrm{F},a}\Big{|}_{\textrm{free}}+i\hbar\int_{-\infty}^{\infty}\frac{d\omega^{\prime}}{2\pi}\\ \sum_{\mathbf{q}}&\frac{(1-f^{a}(E^{a}_{\mathbf{q},\sigma}))V^{\textrm{S},>,aa}_{\mathbf{k}\mathbf{q}\mathbf{k}\mathbf{q}}(\omega^{\prime})+f^{a}(E^{a}_{\mathbf{q},\sigma})V^{\textrm{S},<,aa}_{\mathbf{k}\mathbf{q}\mathbf{k}\mathbf{q}}(\omega^{\prime})}{E^{a}_{\mathbf{k}\sigma}-E^{a}_{\mathbf{q},\sigma}-\hbar\omega^{\prime}}\,.\end{split}$$

$$\begin{array}{rl}& {\epsilon }_{\mathbf{q}}^{\text{ret}}(\omega )={\epsilon }_{\mathbf{q}}^{\text{ret},\text{RPA}}(\omega )\\ & -V_{\mathbf{q}}\frac{1}{\mathcal{A}}\sum _{\nu {\nu }^{\prime }\sigma {\sigma }^{\prime }\mathbf{Q}}{\left|M_{\nu {\nu }^{\prime }\mathbf{Q}}^{\sigma {\sigma }^{\prime }}(\mathbf{q})\right|}^2\frac{n_{ab}^B(E_{\nu \mathbf{Q}}^{\sigma {\sigma }^{\prime }})-n_{ab}^B(E_{{\nu }^{\prime }\mathbf{Q}-\mathbf{q}}^{\sigma {\sigma }^{\prime }})}{E_{\nu \mathbf{Q}}^{\sigma {\sigma }^{\prime }}-E_{{\nu }^{\prime }\mathbf{Q}-\mathbf{q}}^{\sigma {\sigma }^{\prime }}-\hslash \omega -i\gamma }\end{array}$$

$$\begin{split}&M^{\sigma\sigma^{\prime}}_{\nu\nu^{\prime}\mathbf{Q}}(\mathbf{q})=\\ &\frac{1}{\mathcal{A}}\sum_{\mathbf{p}}\psi^{\sigma\sigma^{\prime}}_{\nu\mathbf{Q}}(\mathbf{p})\Big{[}(\psi^{\sigma\sigma^{\prime}}_{\nu^{\prime}\mathbf{Q}-\mathbf{q}}(\mathbf{p}+\mathbf{q}))^{*}-(\psi^{\sigma\sigma^{\prime}}_{\nu^{\prime}\mathbf{Q}-\mathbf{q}}(\mathbf{p}))^{*}\Big{]}\end{split}$$

$$\begin{array}{rl}{\hat{V}}_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{k}{\mathbf{k}}^{\prime }}^{\text{S},ab}(\omega )& =V_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{k}{\mathbf{k}}^{\prime }}^{\text{S},>,ab}(\omega )-V_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{k}{\mathbf{k}}^{\prime }}^{\text{S},<,ab}(\omega )\\ & =V_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{k}{\mathbf{k}}^{\prime }}^{ab}2i\text{Im}{\epsilon }_{\mathbf{k}-{\mathbf{k}}^{\prime }}^{-1,\text{ret}}(\omega ).\end{array}$$

$$\displaystyle\mathbf{U}_{\mathrm{diag}}(q)=\left(\begin{array}[]{ccc}U_{1}(q)&0&0\\ 0&U_{2}&0\\ 0&0&U_{3}\end{array}\right),$$

$$U_i=⟨e_i⟩\mathbf{U}{e}_i$$

$$e_1=\frac{1}{\sqrt{3}}\left(\begin{array}{c}1\\ 1\\ 1\end{array}\right),e_2=\frac{1}{\sqrt{6}}\left(\begin{array}{c}+2\\ -1\\ -1\end{array}\right),e_3=\frac{1}{\sqrt{2}}\left(\begin{array}{c}0\\ +1\\ -1\end{array}\right).$$

$$U_1(q)=\frac{3e^2}{2{\epsilon }_{0}A}\frac{1}{q(1+\gamma q+\delta q^2)},$$

Full text

Excitons play a prominent role in the optical properties of atomically thin transition metal dichalcogenide (TMDC) semiconductors due to electron-hole Coulomb interaction being strongly enhanced by carrier confinement and reduced dielectric screening. This suggests an interpretation of experimental results as well as theoretical prediction in terms of excitons rather than unbound electrons and holes. Zhang et al. (2014); Berghäuser and Malic (2014); Sun et al. (2014); Schmidt et al. (2016); Selig et al. (2016) On the other hand, it is well-known that at a certain excitation density of electron-hole pairs the Mott transition is observed. Shah et al. (1977); Manzke et al. (2012); Chernikov et al. (2015) Here a phase where excitons and unbound carriers can coexist evolves into a fully ionized electron-hole plasma, which shows the significance of plasma effects to the physics of excited semiconductors.

Since excitons are more or less neutral compound particles, many-particle renormalization and screening effects in an exciton gas are very different from those in a plasma of unbound electrons and holes, which we refer to in the following as quasi-free carriers. For this reason it is highly desirable to quantify the relative importance of excitonic and plasma effects over a wide range of electron-hole excitation densities and for various material and externally controllable parameters. It has already been suggested to tune exciton binding energies by electrical doping Chernikov et al. (2015) and some effort has been devoted to study the influence of dielectric screening on excitons in TMDC semiconductors. Lin et al. (2014); Ugeda et al. (2014); Latini et al. (2015); Trolle et al. (2017)

In the past, a very powerful scheme has been developed to theoretically describe the so-called ionization equilibrium between the formation of bound particles and their dissociation into quasi-free particles Kremp et al. (1984a, b); Zimmermann and Stolz (1985); Kremp et al. (1993, 2005); Semkat et al. (2009), with applications to atomic plasmas and highly excited semiconductors. In general, a dominance of free-carrier plasma over excitons is expected at very low excitation density due to entropy ionization of excitons Mock et al. (1978) and above the so-called Mott density due to pressure ionization of excitons. The fraction of excitons reaches a maximum at elevated densities still below the Mott transition and can assume values of practically $100\%$ in TMDCs, though being sensitive to various externally controllable parameters. The scheme relies on the assumption of a quasi-equilibrium between plasma and excitons being established before electron-hole recombination sets in. Ultrafast equilibration is facilitated by efficient carrier-carrier Steinhoff et al. (2016) and carrier-phonon interaction Selig et al. (2016) as well as exciton formation Ceballos et al. (2016); Steinleitner et al. (2017) after optical excitation, see Ref. 23 for a review.

Experimental verification of the ionization equilibrium has been achieved in GaAs quantum wells using THz spectroscopy to probe transitions between 1s- and 2p-exciton states. Kaindl et al. (2003); Koch et al. (2006) A similar technique in the mid-infrared range has been applied recently to monolayer WSe2. Steinleitner et al. (2017) Alternatively, the fractions of excitons and plasma can be determined from their contributions to photoluminescence (PL) spectra Chatterjee (2004) in combination with additional PL simulations. As we suggest below, other ways to quantify the degree of exciton formation are angular-resolved photoemission spectroscopy (ARPES) and scanning tunneling spectroscopy (STS).

In this paper, we combine calculations of material-realistic band structures and Coulomb matrix elements of the monolayer TMDC materials MX2 (M$=$W,Mo and X$=$S,Se) with the state-of-the-art theory of ionization equilibrium, which we briefly introduce in the following section before discussing the results in detail. The theory is based on GW- and T-matrix self-energies describing the excited carriers and the effect of frequency-dependent screening. We also include excitonic screening that we find to be relevant and that is usually not taken into account. On these grounds, we study the influence of experimental and device-relevant parameters like dielectric screening, temperature and carrier doping on the ionization equilibrium and the Mott density.

I Spectral Functions and Exciton Satellites.

To examine the equilibrium properties of excited carriers in TMDCs, we follow the approach developed in Refs. 13; 14; 15; 16 and reviewed in Ref. 18. We use the quantum-statistical expression for the carrier density $n_{a}$ of the species $a$, which can be electrons or holes, as a function of temperature $T$ and chemical potential $\mu_{a}$ as a starting point:

[TABLE]

$f^{a}(\omega)$ denotes the Fermi distribution function depending on $\mu_{a}$ and $T$, $\mathcal{A}$ is the crystal area and $A_{\mathbf{k}\sigma}^{a}(\omega)=2i\textrm{Im}G^{\textrm{ret},a}_{\mathbf{k}\sigma}(\omega)$ is the spectral function of the single-particle state $\ket{\mathbf{k}\sigma a}$ related to the retarded single-particle Green’s function

[TABLE]

The self-energy $\Sigma_{\mathbf{k}\sigma}^{\textrm{ret},a}(\omega)$ accounts for many-particle effects giving rise to renormalizations of the single-particle band structure $\varepsilon_{\mathbf{k}\sigma}^{a}$ as well as contributions of bound states. For a given self-energy, the inversion of Eq. (1) yields the chemical potential $\mu_{a}(n_{a},T)$ for each species and therefore any thermodynamic property of the system in the grand canonical description.

As we describe in detail in the Methods section, by using a T-matrix self-energy in screened ladder approximation and assuming small quasi-particle damping, we obtain a spectral function $A^{a}(\omega)$ in the so-called extended quasi-particle approximation. It exhibits poles for quasi-free and bound carriers as shown in Fig. 1. From the multiple valleys in the single-particle band structure of electrons and holes (Fig. 1(a)), a rich spectrum of bound states emerges (Fig. 1(b)), that contains a variety of dark excitons with large total momentum $Q$ besides the bright $K$-valley excitons commonly referred to as A and B. The dark excitons, though playing a minor role in optical experiments, are essential to the description of the ionization equilibrium. Various bound states are reflected in the low-energy satellites of the single-particle spectral function. This shows that excitonic contributions are expected to be observed in experiments that are sensitive to these spectral properties. In ARPES Damascelli et al. (2003), momentum-resolved images of the electron spectral function comparable to Fig. 1(c) are obtained, weighted with Fermi distribution functions that are defined by the chemical potential $\mu_{e}$ and temperature $T$. For a fixed quasi-momentum state as shown in Fig. 1(e), this is typically referred to as energy distribution curve. On the other hand, STS Fischer et al. (2007) probes the local density of states (LDOS) and thus momentum-averaged spectral functions of electrons and holes, which are displayed in Figs. 1(c) and (d). We therefore propose to use these experimental techniques to spectrally distinguish between excitons and quasi-free carriers. This opens the possibility to quantify the degree of exciton ionization and verify the results we present in this paper.

II Degree of Ionization and Mott Transition.

According to Eq. (1), the spectral function can be used to separate the total electron and hole density ($a=e,h$),

[TABLE]

into contributions from quasi-free carriers and from carriers bound as excitons, where the excitons are approximately treated as bosons. Hence the properties of the excited semiconductor at a given temperature and excitation density are defined by the density of electrons $n_{e}$, the density of holes $n_{h}$ and the density of excitons $n_{X}$. A certain degree of ionization of the excited carriers

[TABLE]

will be established that is determined by the ionization equilibrium between electrons, holes and excitons. While for optical excitation, equal densities of electrons and holes are generated, we distinguish here between electron and hole ionization to also include the effect of carrier doping where electron and hole densities are different.

Using single-particle band structures and bound-state spectra, which are determined as discussed in the Methods section, we solve Eq. (3) numerically to obtain the degree of ionization $\alpha_{a}$ in various TMDC materials under different experimental conditions. The results are collected in Fig. 2 and exhibit the behaviour of ionization degree as a function of excitation density. There are different regimes of ionization to be observed. At high excitation densities between $3\times 10^{12}$ cm*-2* and $1\times 10^{13}$ cm*-2*, depending on experimental parameters, efficient screening and many-particle renormalizations lead to a full ionization of excited carriers, which is known as Mott effect. At lower densities around $n_{a}=1\times 10^{12}$ cm*-2*, excitons dominate the physical properties of TMDCs for the parameters studied here due to the large exciton binding energies and density of states that is mostly given by dark excitons. Bright excitons with very small momenta that are optically active make up only a tiny fraction of the total exciton density. As Fig. 2(b) shows, an efficient tuning knob for the degree of ionization is the dielectric screening by the environment, which can change over a wide range depending on the experimental situation or device realization in which the TMDC monolayer is the active material. The reason is the strong impact of dielectric screening on the exciton binding energies. Typical examples for substrates are Borofloat ($\varepsilon=2$), SiO2 ($\varepsilon=3.9$) and sapphire ($\varepsilon=10$). The dielectric constant of the environment on top of the monolayer is often given by the vacuum value. On the other hand, in devices the TMDC monolayer is fully encapsulated by dielectric material. As an example we consider a full dielectric environment with $\varepsilon=10$, which might be either sapphire or additional layers of TMDC material in a vertical heterostructure whose main influence on the excitons is the dielectric screening. Andersen et al. (2015) We find that the minimal degree of ionization can be tuned from below $0.1\%$ ($99.9\%$ excitons) at weak dielectric screening to about $30\%$ at strong screening, while the Mott density is lowered at the same time by roughly a factor of $3$. The second important parameter that is relevant to applications of TMDC monolayers is the doping with additional carriers which might be either intrinsic or induced by external electric fields in a capacitor structure. Here the fractions of ionized electrons and holes, $\alpha_{e}$ and $\alpha_{h}$, are discussed separately as the densities of the species are not equal anymore. We consider hole doping of WS2, but similar results are expected in case of electron doping. Looking at Figs. 2 (e) and (f), we find that even at weak doping the minority carriers are practically all bound as excitons below the Mott transition. On the other hand, at higher doping levels an increasing fraction of majority carriers exists as quasi-free plasma due to missing partners for exciton formation. As a function of minority-carrier density the Mott transition is lowered by roughly the density of doped excess carriers. From this we conclude that at doping levels above $10^{13}$ cm*-2* neither dark nor bright excitons will exist in any case. This is supported by the experimental estimate for doping-induced ionization at several $10^{13}$ cm*-2*. Chernikov et al. (2015) Another crucial parameter is the temperature, see Fig. 2(c), which can vary in experiments or devices due to heating of the active material under strong optical or electrical pumping. This effect has been taken into account to explain the observed exciton-to-plasma ratio in monolayer WSe2 in Ref. 22. At room temperature and even at elevated temperatures up to $700$ K excitons clearly dominate below the Mott transition. At the same time, the Mott density slightly increases with temperature due to weaker renormalizations of the quasi-particle gap. It turns out that strain is no efficient tuning knob as both bright and dark excitons contribute to the ionization equilibrium, although bright excitons are preferred in moderately tensile-strained TMDCs, see the discussion in the Methods section. A comparison of different TMDC materials shows that excitons are slightly more important in molybdenum- than in tungsten-based TMDCs due to the larger binding energies, which leads to larger Mott densities.

When approaching the Mott density from the low-density side, many-particle renormalizations become increasingly important, cf. Eq. (23) in the Methods section. Exchange interaction and efficient screening due to free carriers as well as excitons reduce the quasi-particle band gap and the exciton binding energies. More and more excitons are ionized, which leads to an increase of efficient free-carrier screening and thereby to a self-amplification of the ionization effect until all excitons are dissociated into an electron-hole plasma and the degree of ionization becomes $\alpha_{a}=1$. Note that $\alpha_{a}$ includes not only bright but also dark excitons with large total momentum for example between $K$ and $\Sigma$ valleys. They may have larger binding energies and be slightly more stable against ionization than bright excitons visible in an optical experiment. Fig. 3 shows an illustration of the Mott effect in terms of the spectral functions in extended quasi-particle approximation, which contain both exciton and quasi-free-particle signatures. At low excitation densities, the only spectral contribution stems from quasi-free carriers at the band edge. With increasing density, the quasi-particle peak is shifted to lower energies due to many-particle renormalizations. At the same time, spectral weight is transferred from the quasi-particle to the bound-state peaks as exciton populations increase, see the explicit expression of the spectral function in Eq. (21). The appearence of several exciton satellites in the hole spectral function is due to different bound states involving electrons either in the $K$- or $\Sigma$-valleys, see Fig. 1 (b). The spectral position of a bound-state peak in the spectral function of carrier $a$ is given by the difference of the corresponding exciton energy $E^{ab}$ and the energy of the second carrier $b$ involved in the bound state. The bound resonance might therefore be interpreted as an effective ionization energy of the actual carrier $a$ with respect to its energy in the quasi-particle band structure. This underlines the fact that excitonic signatures are observable in experiments that are sensitive to the single-particle spectral function, such as ARPES and STS. Note that while the amplitudes of bound-state resonances in the spectral functions are relatively small, observables like the carrier density (1) and the photoemission intensity involve weighting with a Fermi function that strongly favors the low-energy resonances over the quasi-free contribution. With increasing excitation density, quasi-particle and excitonic resonances approach until above the Mott density all excitons are ionized and only a quasi-particle peak remains. Fig. 4 shows the reduction of the quasi-particle gap until the Mott transition appears around $n_{a}=8\times 10^{12}$ cm*-2*.

An alternative picture of the interacting electrons and holes, that is consistent with the extended quasi-particle approximation, is the so-called chemical picture in which excitons are considered as a new particle species besides electrons and holes. Kremp et al. (2005); Semkat et al. (2009) They are characterized by a chemical potential

[TABLE]

with bound-state energies $E_{\nu}$ that are given by the relative motion of electron and hole, and an ideal Bose distribution function. In the chemical picture, solving Eq. (3) corresponds to an adaption of the chemical potentials of the different particle species, namely electrons, holes and excitons, as in a chemical reaction. These considerations are consistent with the theory based on spectral functions, that we use to obtain all numerical results presented in this paper. Only for the purpose of illustration, we simplify the theory considering the nondegenerate case ($f^{a}(E^{a}_{\mathbf{k}\sigma})\ll 1$) and a single band-structure valley for electrons and holes each. Then a Saha equation can be formulated that determines the degree of ionization:

[TABLE]

In analogy to the usual mass action law, $I_{\nu}^{\textrm{eff}}=\Delta E_{\textrm{Gap}}-E_{\nu}$ can be interpreted as an effective ionization potential of excitons that corresponds to the exciton binding energy, see also the inset in Fig. 4. It is obvious from Saha’s equation, that a large exciton binding energy and low temperature favor the formation of excitons versus the dissociation into an unbound electron-hole plasma. The ionization potential depends on excitation density as a consequence of the excitation-induced lowering of the band continuum edge given by $\Delta E_{\textrm{Gap}}$ and the shift of the bound-state energy $E_{\nu}$. The bound-state shift on the other hand is a net result of band-gap shrinkage, screening of exciton binding energy and Pauli blocking Steinhoff et al. (2014) and is much weaker than the band-gap shift due to compensation effects. In the end, the ionization potential is lowered with increasing excitation density until at $I_{\nu}^{\textrm{eff}}=0$ the bound state vanishes and merges with the continuum edge, which is the Mott effect.

A striking observation in Fig. 2 is the degree of ionization approaching unity at low excitation densities, which is somewhat counter-intuitive but can be understood from a thermodynamical point of view. The potential that is minimized by the many-particle system is the free energy $F=U-TS$. At low densities and fixed temperature, the entropy $S$ gained by a dissociation of an exciton into two separate particles overcompensates the reduction of internal energy $U$ by the exciton binding energy $E_{B}$. Hence the so-called entropy ionization already discussed by Mock et al. Mock et al. (1978) is connected to the huge phase space available for quasi-free carriers in the low-density limit. We may clarify this using the entropy of an ideal gas with $N$ particles in a volume $V$ as given by the Sackur-Tetrode equation:

[TABLE]

where $c(T)$ is a temperature-dependent parameter. Obviously, the dissociation of an exciton gas ($N$ particles) into a free electron-hole plasma ($2N$ particles) yields the entropy $\Delta S=Nk_{B}\textrm{ln}\left(n^{-1}\right)$ with $n=N/V$ up to some additive constant. It follows that the critical density $n_{\textrm{crit}}$ below which the free energy is dominated by entropy essentially scales as $\textrm{exp}\left(-E_{B}/k_{B}T\right)$ with temperature.

Although the extended quasi-particle approximation and the chemical picture are very descriptive, we have to be aware of their limitations. The approach relies on the assumption of a quasi-equilibrium of both types of carriers, which tends to overestimate the fraction of bound carriers as exciton formation takes a certain time after excitation of electron-hole pairs. Nevertheless, we expect the approach to give a reasonable quantitative description of the ionization equilibrium, particularly of the trends that can be expected under variation of external experimental parameters. A rather fundamental discussion is concerned with the Mott transition as a first-order phase transition between an exciton gas and a fully ionized electron-hole plasma. Kremp et al. (2005); Semkat et al. (2009) The phase transition would be connected to an instability of thermodynamic functions that manifests itself in an ambiguity of $\alpha_{a}$ in a certain region below the Mott density. Due to excitation-induced broadening of the two-particle states, which is assumed small in our approach, and the shrinkage of the ionization potential towards the Mott transition, quasi-free and bound carriers cannot really be separated in this density regime. We avoid this regime as a more sophisticated theory including full spectral functions and exciton-exciton interaction would be required. Also, screening in a correlated many-particle system near the Mott transition is an intricate problem Kraeft et al. (1986). However, we observe that taking into account excitonic screening is necessary to obtain meaningful results around the minima of $\alpha_{a}$ in the exciton-dominated regime. Otherwise, coming from the low-density side of the ionization curves in Fig. 2, there would be no mechanism to efficiently break up the excitonic binding and trigger the transition to an ionized plasma. We discuss the contribution of excitons to screening in the following section.

Another prominent feature of TMDC semiconductors is the formation of trions, which could in principle be included as additional particle species in the spirit of the chemical picture. Kremp et al. (1993) In practice, obtaining bound-trion spectra on the same footing as excitons is a very challenging task of its own that is beyond the scope of this paper. Due to the relatively small trion binding energies, we assume that their influence on our results are negligible.

III Excitonic Screening.

In the spirit of the extended quasi-particle approximation to the spectral function, there are two types of contributions to excited-carrier screening of the Coulomb interaction, the metal-like free-carrier screening and dipolar screening due to bound excitons.

The screening can be characterized by the plasmon spectral function, see Eq. (26), that contains excitations in the interacting electron-hole plasma as poles in the $q$-$\omega$-plane, see Fig. 5. In the exciton-dominated regime shown in Fig. 5(a), besides the usual 2-d free-carrier contribution at small energies and small momenta a broad resonance above $150$ meV appears. It stems from transitions between 1s- and 2s-like exciton states, see Fig. 1(b), and also from comparable transitions between exciton states with large momenta. There are contributions at smaller energies as well that can not be as easily distinguished from free-carrier screening. At large densities beyond the Mott transition the plasmon spectral function shows a pronounced peak structure with a square-root-like behaviour at small momenta which has been discussed for TMDCs in Ref. 31 and which is typical for a two-dimensional electron gas. Haug and Koch (1993) Although excitons are expected to be much less polarizable than a free electron-hole plasma and hence contribute less to screening, at elevated excitation densities with a large fraction of carriers bound as excitons, their contribution can be essential. From the many-particle renormalization caused by free-carrier and excitonic screening, we deduce that in monolayer TMDC semiconductors excitonic screening is less efficient by two to three orders of magnitude for comparable excitation densities. Nevertheless, in a regime where more than $99\%$ of carriers are bound as excitons, excitonic screening still yields a significant contribution. As the plasmon spectral function is directly observable by electron energy loss spectroscopy (EELS) Kogar et al. (2016), we suggest to use this technique to explore exciton signatures in the dielectric function eperimentally.

IV Conclusion.

The exciton ionization equilibrium in monolayer TMDC semiconductors has been studied for various material as well as experimentally and device-relevant external parameters on the basis of an ab initio description of the electronic band structure and Coulomb interaction. We observe entropy ionization of excitons at low excitation densities and a Mott transition to a fully ionized plasma at high densities between $3\times 10^{12}$ cm*-2* and $1\times 10^{13}$ cm*-2* depending on experimental parameters. Below the Mott transition, excitons become dominant in all cases with maximal fractions of excitons between $70\%$ and more than $99.9\%$. The most efficient tuning knobs are dielectric screening of the Coulomb interaction via the choice of dielectric environment and carrier doping that can induce complete ionization above a level of $10^{13}$ cm*-2*. We suggest that fingerprints of excitonic contributions can be observed in ARPES and STS experiments, which are sensitive to the single-particle spectral functions. Moreover, we find that excitonic screening, although two to three orders of magnitude less efficient than free-carrier screening at comparable excitation densities, plays an important role in the description of ionization equilibrium. Exciton signatures in the dielectric function suggest EELS as another way to study the ionization equilibrium in excited semiconductors.

V Methods.

V.1 Theory of Ionization Equilibrium.

We start from the general expression for the carrier density (1) and the spectral function

[TABLE]

In the limit of small quasi-particle damping $\textrm{Im}\,\Sigma^{\textrm{ret},a}\ll\textrm{Re}\,\Sigma^{\textrm{ret},a}$, the spectral function can be expanded in linear order of $\textrm{Im}\,\Sigma^{\textrm{ret},a}$ yielding the carrier density in so-called extended quasi-particle approximation

[TABLE]

where the quasi-particle energy $E^{a}_{\mathbf{k}\sigma}$ is given by $E^{a}_{\mathbf{k}\sigma}=\varepsilon_{\mathbf{k}\sigma}^{a}+\textrm{Re}\,\Sigma_{\mathbf{k}\sigma}^{\textrm{ret},a}(E^{a}_{\mathbf{k}\sigma})$ and $\mathcal{P}$ denotes the Cauchy principal value. Kremp et al. (1984a); Semkat et al. (2009) The total density is divided into contributions from quasi-free particles and correlated particles, the latter being either in bound or scattering many-particle states.

The spectral function in extended quasi-particle approximation corresponding to this separation into free and correlated carriers is given by

[TABLE]

with $\Gamma_{\mathbf{k}\sigma}^{a}(\omega)=\textrm{Im}\Sigma_{\mathbf{k}\sigma}^{\textrm{ret},a}(\omega)\frac{1}{\pi}\frac{d}{d\hbar\omega}\frac{\mathcal{P}}{\hbar\omega-E^{a}_{\mathbf{k}\sigma}}$ and the renormalization factor $Z_{\mathbf{k}\sigma}^{a}=\int d\hbar\omega\Gamma_{\mathbf{k}\sigma}^{a}(\omega)$. The first term describes quasi-free particles at renormalized energies. Their spectral weight is reduced according to the renormalization factor to account for correlated carriers, which are spectrally described by the second term.

To evaluate the expressions (9) and (10), we have to choose an approximation for the self-energy $\Sigma^{\textrm{ret},a}(\omega)$. The real and imaginary parts of $\Sigma$ determine the quasi-particle energies and the correlated part of the carrier density, respectively. An appropriate choice is the screened ladder approximation Stolz and Zimmermann (1979); Kremp et al. (1984a, 1993) $\Sigma(\omega)=\Sigma^{\textrm{H}}+\Sigma^{\textrm{GW}}(\omega)+\Sigma^{\textrm{T}}(\omega)$ that takes into account screening of Coulomb interaction due to excited carriers as well as the formation of bound two-particle states and consists of Hartree, GW and T-matrix contributions. We assume that renormalizations due to the Hartree self-energy are small compared to exchange and correlation effects. In the T-matrix contribution, we neglect exchange terms and assume static screening so that the T-matrix depends only on one instead of three frequency arguments. Thus we obtain for the imaginary part of the self-energy using the generalized Kadanoff-Baym ansatz Kremp et al. (2005):

[TABLE]

Here we applied thermal equilibrium relations for the screened Coulomb interaction Kremp et al. (2005):

[TABLE]

$\varepsilon^{-1,\textrm{ret}}_{\mathbf{q}}(\omega)$ is the longitudinal dielectric function describing screening due to excited carriers and $n^{B}(\omega)$ is the Bose distribution function of the elementary plasma excitations called plasmons. $V^{ab}$ denotes Coulomb matrix elements between species $a$ and $b$ which contain dielectric screening due to carriers in the ground state and due to the environment but no screening due to excited carriers. $T^{\prime\prime}$ denotes the T-matrix with the two lowest-order terms subtracted from the ladder expansion and is discussed in the following subsection.

V.2 T-matrix and Bound Carriers.

The T-matrix in statically screened ladder approximation describing bound and scattering two-particle states between carrier species $a$ and $b$ obeys a Lippmann-Schwinger equation (LSE)

[TABLE]

where $\mathcal{G}^{\textrm{ret},ab}(\omega)$ is the free two-particle Green’s function in the particle-particle channel. The corresponding interacting two-particle Green’s function $G^{\textrm{ret},ab}_{2}(\omega)$ fulfills a Bethe-Salpeter equation, that has been discussed in detail in Kraeft et al. (1986); Bornath et al. (1999) and is equivalent to the LSE. We will exploit this fact later when solving the LSE and evaluating the T-matrix self-energy. In its homogeneous form, the BSE in static ladder approximation is given by

[TABLE]

Diagonalization yields bound states $\ket{\nu\sigma\sigma^{\prime}\mathbf{Q}}$ and eigenenergies $E^{\sigma\sigma^{\prime}}_{\nu\mathbf{Q}}$. We drop the indices $a$ and $b$ here, assuming that only two-particle states between different carrier species are involved. Due to the translational invariance of the crystal, the bound states can be classified by the total exciton momentum $\mathbf{Q}$ as discussed in Qiu et al. (2015). Here we neglect the effect of electron-hole exchange interaction that leads to a fine-structure splitting of excitons and trions Qiu et al. (2015); Jones et al. (2016); Plechinger et al. (2016) in the meV range, which is small compared to the exciton binding energies of several hundred meV. As a consequence, electron and hole spins, which are already good quantum numbers in monolayer TMDC materials due to crystal symmetry, also classify the bound states. For each total momentum and spin combination a series of excitons exists, which is labeled here by $\nu$, analogue to the angular momentum states of Hydrogen-like Wannier excitons. Due to the two-dimensional nature of monolayer TMDCs and the related strong momentum dependence of dielectric screening, nontrivial exciton series deviating from a Hydrogen-like spectrum are observed. Qiu et al. (2013); Chernikov et al. (2014); Berghäuser and Malic (2014) The eigenenergies decompose into a part from the relative motion of electron and hole and a kinetic part depending on the total momentum: $E^{\sigma\sigma^{\prime}}_{\nu\mathbf{Q}}=E^{\sigma\sigma^{\prime}}_{\textrm{rel},\nu}+E^{\sigma\sigma^{\prime}}_{\textrm{kin},\mathbf{Q}}$. We can use Bloch basis functions to find a representation of the bound states corresponding to exciton wave functions

[TABLE]

where $\mathbf{k}$ conventionally denotes the hole momentum, while the electron momentum is fixed via the total momentum.

An explicit expression for the T-matrix can be obtained by writing the LSE (13) in the basis of two-particle eigenstates $\ket{\nu\sigma\sigma^{\prime}\mathbf{Q}}$ as shown in detail in Ref. 17. Since the BSE represents a generalized eigenvalue problem, the eigenstates form a biorthogonal basis. The procedure yields a spectral representation of the T-matrix in operator form that is referred to as bilinear expansion:

[TABLE]

with the Pauli blocking factor $N_{ab}=1-f^{a}-f^{b}$, the operator of kinetic energy $E_{\textrm{kin}}$ and the eigenstate of the adjoint BSE $\bra{\nu\widetilde{\sigma\sigma^{\prime}}\mathbf{Q}}$. The bilinear expansion is used in the following to evaluate the imaginary part of the self-energy (11) and thereby the contribution of correlated carriers.

V.3 Separation of Bound and Quasi-Free Carriers.

Inserting Eq. (11) into Eq. (9) and noting that neither the GW self-energy nor the two lowest T-matrix terms contribute to the carrier density Stolz and Zimmermann (1979), we obtain Kremp et al. (1993); Semkat et al. (2009)

[TABLE]

n^{B}_{ab}(\omega)=\big{[}\textrm{exp}(\beta(\hbar\omega-\mu_{a}-\mu_{b}))-1\big{]}^{-1} is the Bose distribution function depending on the chemical potentials of both carrier species. Eq. (17) contains contributions of both bound two-particle states (below the single-particle gap $E_{\textrm{Gap}}$) and scattering two-particle states, the latter being explicitely given in Refs. 14; 15; 16. The renormalization factor $Z_{\mathbf{k}\sigma}^{a}$ of the quasi-particle resonance in the spectral function (10) enters the contribution of correlated carriers as Pauli-blocking factor and as correction to the two-particle scattering spectrum. To simplify the following discussion, we neglect the contribution $n_{\textrm{scatt}}$ of scattering states beyond the quasi-free carriers and consider only the bound-state contribution given by the real-frequency poles of the T-matrix Kremp et al. (1984b):

[TABLE]

Using Eq. (18), we arrive at the final expression for the carrier density:

[TABLE]

The total carrier density separates into contributions from quasi-free carriers and from carriers bound as excitons according to the two poles in the spectral function $A^{a}(\omega)$. For a specific material, the ionization equilibrium has to be computed numerically. The electron and hole chemical potentials are determined by adapting the Fermi functions $f^{a}(E^{a}_{\mathbf{k}\sigma})$ of electrons and holes to a given density of quasi-free carriers at the quasi-particle energies $E^{a}_{\mathbf{k}\sigma}$. As the chemical potentials also enter the bound-carrier density via the Bose function $n^{B}_{ab}$, Eq. (19) represents an implicit equation for the fraction of quasi-free carriers $\alpha_{a}=n_{\textrm{free}}^{a}/n_{a}$, that has to be solved self-consistently with the quasi-particle energies in GW approximation, see Eq. (23), and the bound-state energies $E^{\sigma\sigma^{\prime}}_{\nu\mathbf{Q}}$. To simplify the procedure, we exploit the fact that shifts of excitonic resonances are naturally much smaller than band-gap shifts, which is due to compensation effects between gap shrinkage and binding-energy reduction. Bornath et al. (1999); Steinhoff et al. (2014) Hence we assume that the exciton spectrum depends only weakly on the excitation density so that we can limit ourselves to the BSE (14) in the limit of zero excitation density.

Consistent with the imaginary part of the self-energy (11), the quasi-particle energies $E^{a}_{\mathbf{k}\sigma}$ contain GW- and T-matrix contributions:

[TABLE]

The GW self-energy is separated into the Fock term and the so-called Montroll-Ward term containing all contributions beyond bare exchange interaction. The T-matrix contribution is explicitely given in Ref. 34 and leads to a blue shift of single-particle energies that is in the nondegenerate case ($f^{a}(E^{a}_{\mathbf{k}\sigma})\ll 1$) caused by the bound-carrier population. At the same time, the Fock self-energy contains exchange interaction with both quasi-free and bound carriers via the extended spectral functions that leads to a lowering of single-particle energies. This can be seen by using the T-matrix self-energy in Eq. (11) to obtain an excitonic contribution to the spectral function (10) given by

[TABLE]

It yields a sharp resonance for each bound state weighted by its Bose population function and the exciton wave functions at the corresponding position in k-space. Note that the spectral positions of the resonances are not given by the bound-state energies $E^{\sigma\sigma^{\prime}}_{\nu\mathbf{Q}}$, which are two-particle quantities, but by an effective binding energy of the carrier in state $\ket{\mathbf{k}\sigma a}$, as $\Gamma$ represents a single-particle spectral function. The Fock self-energy Kremp et al. (2005) can then be expressed in terms of the spectral function using the Kubo-Martin-Schwinger relation for the propagators $G^{<}(\omega)$ in thermal equilibrium:

[TABLE]

The first contribution to the Fock self-energy scales, besides the Coulomb matrix elements, with the free-carrier density while the second contribution scales with the density of bound carriers. It turns out that similar to exchange interaction with free carriers, bound-carrier exchange leads to k-dependent renormalizations according to the exciton wave functions and populations that are contained in the population factor $f^{a,\textrm{bound}}_{\mathbf{k}^{\prime}\sigma}$. As a conclusion, the real part of the self-energy (20) contains quasi-particle renormalizations due to exciton populations via the T-matrix in two different places that act in opposite directions. We assume that these renormalizations cancel to a large degree and focus on the free-carrier contributions in accordance with Refs. 16; 18. Then we obtain for the quasi-particle energies:

[TABLE]

In a similar manner as for the Fock self-energy, the full spectral functions could be used to evaluate the Montroll-Ward self-energy in Eq. (23). Due to the spectral structure of the self-energy, however, renormalizations of the single-particle band structure caused by bound carriers involve a denominator of the order of the exciton binding energy, which is very off-resonant. Therefore the Montroll-Ward term is evaluated using spectral functions for quasi-free carriers.

V.4 Screening due to Excited Carriers.

In the spirit of the extended quasi-particle approximation, screening of the Coulomb interaction due to both free carriers and bound excitons is taken into account. The free-carrier screening is treated in RPA with a Lindhard dielectric function Steinhoff et al. (2014), while the excitonic polarizibilities are calculated as described in Röpke and Der (1979); Röpke et al. (1981):

[TABLE]

with matrix elements

[TABLE]

and exciton wave functions $\psi^{\sigma\sigma^{\prime}}_{\nu\mathbf{Q}}(\mathbf{p})$ as defined above. The momentum and frequency dependence of screening is characterized by the plasmon spectral function

[TABLE]

V.5 Ab initio Single- and Two-Particle Properties.

For a material-realistic description of the strong Coulomb interaction in atomically thin systems, we use band structures $\varepsilon^{a}$ and Coulomb matrix elements $V^{ab}$ from material-realistic G0W0 calculations for freestanding slabs of various monolayer TMDCs as a basis, see the Supporting Information and Ref. 29. The matrix elements already take into account dielectric screening due to charge carriers in the ground state of the system. Furthermore, we account for additional screening provided by a possible dielectric environment, like a substrate, as described in the Supporting Information. The main effect of the additional screening is a reduction of exciton binding energies, which are key quantities for the exciton-plasma balance. The corresponding reduction of quasi-particle band gaps on the other hand plays no role here, as all bound-state energies and excited-carrier chemical potentials are measured relative to the gap. By convention, we choose the gap between the spin-up conduction and valence bands at $K$ as reference which belong to the A-exciton resonance.

The calculation of ionization equilibrium requires knowledge of two-particle energies, which are obtained by solving the Bethe-Salpeter equation (14). Fig. 1 (b) shows the full exciton spectrum for monolayer WS2 on a SiO2 substrate, for which we assume a dielectric constant $\varepsilon=3.9$. Besides the $K$-excitons at $Q=0$ the spectra exhibit a rich structure of various inter-valley excitons with large total momenta. Even for direct-band-gap materials, inter-valley excitons between $K$- and $\Sigma$-valleys at $Q\approx 6$ nm*-1* may be lower than intra-valley $K$-excitons due to a larger exciton binding energy. As can be seen in Fig. 6 (a) the binding energies of the lowest $Q=0$-excitons (“1s”) are comparable for all considered materials, the main difference being the splitting between excitons involving different conduction bands, which is larger for tungsten-based materials. For MoS2 bound states belonging to the B exciton are visible at $-175$ meV corresponding to the valence-band splitting of about $140$ meV. Higher exciton states (“2s”) appear as well. Fig. 6 (b) shows that dielectric screening from the environment has a strong impact on the exciton spectrum as it can reduce exciton binding energies significantly Latini et al. (2015). This has a major influence on the ionization equilibrium, as we discuss in the main text. Another tunable parameter in experiments is the application of strain to the TMDC monolayer, which leads to changes in the band structure Conley et al. (2013). This is reflected by MoS2 changing from an indirect to a direct semiconductor in the exciton picture under tensile strain, see Fig. 6 (c).

V.6 Numerical Details.

We calculate the ionization equilibrium from the fraction of quasi-free carriers as root of the implicit equation (19). The two highest valence and two lowest conduction bands are considered to cover all excitons that are relevant in a quasi-equilibrium situation. Likewise, we limit the Brillouin zone to circles with radius $2.7$ nm*-1* around the $K$,$K^{\prime}$,$\Sigma$,$\Sigma^{\prime}$ and $\Gamma$ points using a Monkhorst-Pack mesh with $30$ mesh points along $\Gamma$-$M$, which yields reasonable convergence of all results. The frequency integrals involved in the Montroll-Ward self-energy (23) are extended from $-600$ to $600$ meV exploiting the relation $V^{\textrm{S},<,ab}_{\mathbf{k}\mathbf{k}^{\prime}\mathbf{k}\mathbf{k}^{\prime}}(-\omega)=V^{\textrm{S},>,ab}_{\mathbf{k}\mathbf{k}^{\prime}\mathbf{k}\mathbf{k}^{\prime}}(\omega)$. For simplicity, we use a dielectric function (24) which is isotropic in momentum by evaluating its dependence on $|\mathbf{q}|$ along the contour $\Gamma$-$K$ and using Coulomb matrix elements $V_{|\mathbf{q}|}$ that are averaged over Wannier orbitals. Both the Lindhard and the excitonic dielectric function (24) are evaluated using ground-state energies and extrapolated to the limit of vanishing phenomenological quasi-particle broadening $\gamma$. For a sum over bound states, only those states that are below the local (in $\mathbf{k}$-space) renormalized band structure measured at the maximum of the corresponding exciton wave functions are taken into account.

Acknowledgement We acknowledge financial support from the Deutsche Forschungsgemeinschaft (JA 14-1 and RTG 2247 ”Quantum Mechanical Materials Modelling”) and the European Graphene Flagship as well as resources for computational time at the HLRN (Hannover/Berlin). We thank Michael Lorke, Christopher Gies and Dirk Semkat for fruitful discussions.

Supporting Information available, containing a description of the ab initio procedures to obtain band structures and screened Coulomb matrix elements.

VI Supporting Information.

Ab-initio based parametrization of the TMDCs

To get an analytic description of the Coulomb interaction in semiconducting TMDC monolayers, we apply the approach as presented in Steinhoff et al. (2014); Liu et al. (2013); Schönhoff et al. (2016) for MoS2 also to MoSe2, WS2, and WSe2. We start with ab initio calculations for density-density like bare $U_{\alpha\beta}(q)$ and screened $V_{\alpha\beta}(q)$ Coulomb interaction matrix elements in the Wannier basis (with $\alpha,\beta\in[d_{z^{2}},d_{xy},d_{x^{2}-y^{2}}]$ ) for the freestanding, undoped TMDC slabs using the FLEUR and SPEX codes FLE (2014); Friedrich et al. (2009, 2010) on discrete $18\times 18\times 1$ q grids. To interpolate the resulting $3\times 3$ matrices analytically we make use of the (sorted) eigenbasis of the bare Coulomb interaction $\mathbf{U}$ by diagonalizing it

[TABLE]

where the diagonal matrix elements are given by

[TABLE]

using the eigenvectors of $\mathbf{U}$ in their long-wavelength limits

[TABLE]

While the leading eigenvalue $U_{1}(q)$ is a function of $q$ the other two eigenvalues can be readily approximated as constants (see Fig. S7). For the analytic description of the former we use

[TABLE]

where $A=\frac{\sqrt{3}}{2}a^{2}$ is the area of the hexagonal unit cell, $a$ is the lattice constant, and $\varepsilon_{0}$ is the vacuum permittivity.

The matrix elements of the screened interaction $\mathbf{V}(q)$ in the eigenbasis of the bare interaction $\mathbf{U}(q)$ are then obtained via

[TABLE]

where $\varepsilon_{i}(q)$ accounts for both the material-specific internal polarizability and the screening by the environment. Its diagonal representation is given by

[TABLE]

Once again, the leading eigenvalue $\varepsilon_{1}(q)$ is a function of $q$ while the other elements are described sufficiently well as constants. Here, the former can be expressed by

[TABLE]

which describes the macroscopic dielectric function of a two-dimensional semiconductor. The parameters $\beta_{i}$ include the screening effects of substrates (see Ref. Rösner et al. (2015)) via

[TABLE]

where the dielectric constants of the substrate ($i=1$) and the superstrate ($i=2$) are introduced. In order to describe the original ab initio data as close as possible we fit $\varepsilon_{\infty}(q)$ using

[TABLE]

and $\varepsilon_{\mathrm{sub},1}=\varepsilon_{\mathrm{sub},2}=1$ as the ab intio calculations were performed for freestanding layers.

As soon as all fitting parameters are obtained (see Tab. 1) the screening of a dielectric environment can be included by choosing $\varepsilon_{\mathrm{sub},1}$ or $\varepsilon_{\mathrm{sub},2}$ correspondingly. Thus we have a closed analytic description of the screened Coulomb interaction $\mathbf{V}(q)$ in the eigenbasis of the bare interaction $\mathbf{U}(q)$ at arbitrary momenta $q$ in the first Brillouin zone. In order to transform it to the original orbital basis we make use of the eigensystem given in Eq. (S32). In fact, this model is appropriate for every two-dimensional semiconductor.

Besides the analytical description of the screened Coulomb matrix elements we make use of a tight-binding Hamiltonian to describe the electronic band structure (as obtained from $G_{0}W_{0}$ calculations) of the TMDC slab. To this end, we utilize the same Wannier basis as described before (see Ref. Groenewald et al. (2016)) and derive a minimal three-band model describing the highest valence and two lowest conduction bands using the Wannier90 package Mostofi et al. (2014). Thereby we solely disentangle our target bands from the rest without performing a maximal localization in order to preserve the original transition metal d-orbital characters. The latter is crucial for the subsequent addition of first and second order Rashba spin-orbit coupling following Ref. Liu et al. (2013), which takes into account the large spin-orbit splitting in the conduction- and the valence-band K valleys.

Throughout the whole model-building process we assume that the substrate mainly affects the internal Coulomb interaction and neglect its influence on the band structure. We have used relaxed lattice constants as given in Tab. 1. The values of the bare and screened Coulomb interaction were extrapolated from vacuum heights of $16\,$Å to $32\,$Å.