Microscopic Coulomb interaction in transition-metal dichalcogenides

For a strictly 2D semiconductor, the basic system Hamiltonian is given by

$$\begin{equation*}H={H}_{0}+{H}_{\mathrm{I}}+{H}_{\mathrm{C}}\end{equation*}$$

where

$$\begin{equation*}{H}_{0}=\sum _{\alpha {\mathbf{k}}_{{\Vert}}}{{\epsilon}}_{\alpha {\mathbf{k}}_{{\Vert}}}{c}_{\alpha {\mathbf{k}}_{{\Vert}}}^{{\dagger}}{c}_{\alpha {\mathbf{k}}_{{\Vert}}}\end{equation*}$$

describes the single-particle part,

$$\begin{equation*}{H}_{\mathrm{I}}=\frac{e}{{m}_{0}c}\sum _{\alpha {\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}\mathbf{A}\cdot {\mathbf{p}}_{\alpha {\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}{c}_{\alpha {\mathbf{k}}_{{\Vert}}}^{{\dagger}}{c}_{{\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}},\end{equation*}$$

contains the light–matter interaction, and

$$\begin{equation*}{H}_{\mathrm{C}}=\frac{1}{2}\sum _{\begin{subarray}{c}\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }\\ {\mathbf{k}}_{{\Vert}}{{\mathbf{k}}_{{\Vert}}}^{\prime }{\mathbf{q}}_{{\Vert}}\ne 0\end{subarray}}{V}_{{\mathbf{q}}_{{\Vert}};{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}}}^{\alpha \beta {\beta }^{\prime }{\alpha }^{\prime }}{c}_{\alpha {\mathbf{k}}_{{\Vert}}-{\mathbf{q}}_{{\Vert}}}^{{\dagger}}{c}_{\beta {{\mathbf{k}}_{{\Vert}}}^{\prime }+{\mathbf{q}}_{{\Vert}}}^{{\dagger}}{c}_{{\beta }^{\prime }{{\mathbf{k}}_{{\Vert}}}^{\prime }}{c}_{{\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}.\end{equation*}$$

represents the Coulomb interaction, respectively. Here, ℏ k_∥ is the in-plane crystal momentum, α combines spin and band index, ${{\epsilon}}_{\alpha {\mathbf{k}}_{{\Vert}}}$ contains the effective single-particle band dispersion, and ${\mathbf{p}}_{\alpha {\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}$ and ${V}_{{\mathbf{q}}_{{\Vert}};{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}}}^{\alpha \beta {\beta }^{\prime }{\alpha }^{\prime }}={V}_{{\mathbf{q}}_{{\Vert}}}\langle \alpha {\mathbf{k}}_{{\Vert}}-{\mathbf{q}}_{{\Vert}}\vert {\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}\rangle \langle \beta {{\mathbf{k}}_{{\Vert}}}^{\prime }+{\mathbf{q}}_{{\Vert}}\vert {\beta }^{\prime }{{\mathbf{k}}_{{\Vert}}}^{\prime }\rangle$ denote the momentum and Coulomb matrix elements, respectively. Furthermore, ${V}_{{\mathbf{q}}_{{\Vert}}}$ is the Fourier transform of the 2D Coulomb potential. In the presence of a dielectric environment, the environmental screening of the Coulomb potential can be included into the definition of the 'bare' Coulomb potential.

The optical response is then given by $\mathbf{j}=-\frac{e}{{m}_{0}}{\sum }_{\alpha {\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}{\mathbf{p}}_{\alpha {\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}\langle {c}_{\alpha {\mathbf{k}}_{{\Vert}}}^{{\dagger}}{c}_{{\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}\rangle$ and can be computed from the Heisenberg equations of motion for the transition amplitudes ${P}_{{\mathbf{k}}_{{\Vert}}}^{\alpha {\alpha }^{\prime }}=\langle {c}_{\alpha {\mathbf{k}}_{{\Vert}}}^{{\dagger}}{c}_{{\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}\rangle$:

$$\begin{align}\hfill \mathrm{i}\hslash \enspace \frac{\mathrm{d}}{\mathrm{d}t}\enspace {P}_{{\mathbf{k}}_{{\Vert}}}^{\alpha {\alpha }^{\prime }}& =\left({{\epsilon}}_{{\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}-{{\epsilon}}_{\alpha {\mathbf{k}}_{{\Vert}}}-\mathrm{i}\gamma \right){P}_{{\mathbf{k}}_{{\Vert}}}^{\alpha {\alpha }^{\prime }}+{\tilde {Q}}_{{\mathbf{k}}_{{\Vert}}}^{\alpha {\alpha }^{\prime }}-{\left({\tilde {Q}}_{{\mathbf{k}}_{{\Vert}}}^{{\alpha }^{\prime }\alpha }\right)}^{{\ast}}\hfill \\ \hfill & \quad +{\left.\mathrm{i}\hslash \enspace \frac{\mathrm{d}}{\mathrm{d}t}\enspace {P}_{{\mathbf{k}}_{{\Vert}}}^{\alpha {\alpha }^{\prime }}\right\vert }_{\text{corr}}.\hfill \end{align} \tag{ 1 }$$

Here,

$$\begin{equation}{\tilde {Q}}_{{\mathbf{k}}_{{\Vert}}}^{\alpha {\alpha }^{\prime }}=\sum _{\beta }{P}_{{\mathbf{k}}_{{\Vert}}}^{\alpha \beta }{\tilde {{\Omega}}}_{{\mathbf{k}}_{{\Vert}}}^{\beta {\alpha }^{\prime }}\end{equation} \tag{ 2 }$$

$$\begin{align}\hfill {\tilde {{\Omega}}}_{{\mathbf{k}}_{{\Vert}}}^{\beta {\alpha }^{\prime }}& =\frac{e}{{m}_{0}c}\mathbf{A}\cdot {\mathbf{p}}_{{\alpha }^{\prime }\beta {\mathbf{k}}_{{\Vert}}}-\sum _{\begin{subarray}{c}\gamma {\gamma }^{\prime }\\ {{\mathbf{k}}_{{\Vert}}}^{\prime }\ne {\mathbf{k}}_{{\Vert}}\end{subarray}}{V}_{{\mathbf{k}}_{{\Vert}}-{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}};{{\mathbf{k}}_{{\Vert}}}^{\prime }}^{{\alpha }^{\prime }\beta \gamma {\gamma }^{\prime }}{P}_{{{\mathbf{k}}_{{\Vert}}}^{\prime }}^{\gamma {\gamma }^{\prime }}\hfill \end{align} \tag{ 3 }$$

contains the sources and the renormalizations of the single particle energies and internal field. For the interband polarization ${P}_{{\mathbf{k}}_{{\Vert}}}^{\nu c}$, we explicitely have

$$\begin{align*}\hfill {\tilde {Q}}_{{\mathbf{k}}_{{\Vert}}}^{\nu c}& -{\left({\tilde {Q}}_{{\mathbf{k}}_{{\Vert}}}^{c\nu }\right)}^{{\ast}}=\left({\tilde {{\Omega}}}_{{\mathbf{k}}_{{\Vert}}}^{cc}-{\tilde {{\Omega}}}_{{\mathbf{k}}_{{\Vert}}}^{\nu \nu }\right){P}_{{\mathbf{k}}_{{\Vert}}}^{\nu c}+\left({f}_{{\mathbf{k}}_{{\Vert}}}^{\nu }-{f}_{{\mathbf{k}}_{{\Vert}}}^{c}\right){\tilde {{\Omega}}}_{{\mathbf{k}}_{{\Vert}}}^{\nu c},\hfill \\ \hfill {\tilde {{\Omega}}}_{{\mathbf{k}}_{{\Vert}}}^{cc}& =\frac{e}{{m}_{0}c}{\mathbf{p}}_{cc{\mathbf{k}}_{{\Vert}}}\cdot \mathbf{A}-\sum _{{{\mathbf{k}}_{{\Vert}}}^{\prime }\ne {\mathbf{k}}_{{\Vert}}}\hfill \\ \hfill & \quad {\times}\left({V}_{{\mathbf{k}}_{{\Vert}}-{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}};{{\mathbf{k}}_{{\Vert}}}^{\prime }}^{cccc}{f}_{{{\mathbf{k}}_{{\Vert}}}^{\prime }}^{c}+{V}_{{\mathbf{k}}_{{\Vert}}-{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}};{{\mathbf{k}}_{{\Vert}}}^{\prime }}^{cc\nu \nu }{f}_{{{\mathbf{k}}_{{\Vert}}}^{\prime }}^{\nu }\right)+\mathrm{N}\mathrm{R},\hfill \\ \hfill {\tilde {{\Omega}}}_{{\mathbf{k}}_{{\Vert}}}^{\nu c}& =\frac{e}{{m}_{0}c}{\mathbf{p}}_{\nu c{\mathbf{k}}_{{\Vert}}}\cdot \mathbf{A}\hfill \\ \hfill & \quad -\sum _{{{\mathbf{k}}_{{\Vert}}}^{\prime }\ne {\mathbf{k}}_{{\Vert}}}{V}_{{\mathbf{k}}_{{\Vert}}-{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}};{{\mathbf{k}}_{{\Vert}}}^{\prime }}^{c\nu \nu c}{P}_{{{\mathbf{k}}_{{\Vert}}}^{\prime }}^{\nu c}+\mathrm{N}\mathrm{R},\hfill \end{align*}$$

and a similar expression for ${\tilde {{\Omega}}}_{{\mathbf{k}}_{{\Vert}}}^{\nu \nu }$. Here, NR refers to nonresonant Auger and pair-creation contributions. As can be recognized, Coulomb renormalizations of the single particleenergies are contained in ${\tilde {{\Omega}}}_{{\mathbf{k}}_{{\Vert}}}^{cc}$ and ${\tilde {{\Omega}}}_{{\mathbf{k}}_{{\Vert}}}^{\nu \nu }$ whereas the attractive electron–hole interaction responsible for the formation of bound excitons is contained in ${\tilde {{\Omega}}}_{{\mathbf{k}}_{{\Vert}}}^{\nu c}$. Correlation effects beyond the mean field approximation are contained in

$$\begin{equation}{\left.\mathrm{i}\hslash \enspace \frac{\mathrm{d}}{\mathrm{d}t}\enspace {P}_{{\mathbf{k}}_{{\Vert}}}^{\alpha {\alpha }^{\prime }}\right\vert }_{\text{corr}}.\end{equation} \tag{ 4 }$$

It can be shown that the dominant correlation effect is the replacement of the 'bare' Coulomb potential in equation (3) by its dynamically screened counterpart, yielding the screened Hartree–Fock approximation [6], whereas if correlation effects are neglected, the renormalizations of the single particle energies correspond to the Hartree–Fock approximation.

If the system is excited with an optical pulse of central frequency ω_L , the induced optical current is dominated by transition amplitudes ${P}_{{\mathbf{k}}_{{\Vert}}}^{\alpha {\alpha }^{\prime }}$ corresponding to dipole allowed transitions with transition energies ${{\epsilon}}_{{\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}-{{\epsilon}}_{\alpha {\mathbf{k}}_{{\Vert}}}\approx \hslash {\omega }_{L}$ that are resonant or nearly resonant with the optical frequency. If the relevant bands are sufficiently isolated, one can restrict the microscopic analysis to the resonant transition amplitudes and occupation numbers ${f}_{{\mathbf{k}}_{{\Vert}}}^{\alpha }={P}_{{\mathbf{k}}_{{\Vert}}}^{\alpha \alpha }$ of those bands and the only coupling to remote bands is via screening of the Coulomb potential. Treating screening within the RPA approximation, we can make the separation

$$\begin{equation}{{\Pi}}_{{\mathbf{q}}_{{\Vert}}}\left(\omega \right)={{\Pi}}_{{\mathbf{q}}_{{\Vert}}}^{\text{GS}}\left(\omega \right)+{\Delta}{{\Pi}}_{{\mathbf{q}}_{{\Vert}}}\left(\omega \right)\end{equation} \tag{ 5 }$$

$$\begin{equation}{W}_{{\mathbf{q}}_{{\Vert}}}\left(\omega \right)=\frac{{V}_{{\mathbf{q}}_{{\Vert}}}}{1-{V}_{{\mathbf{q}}_{{\Vert}}}{{\Pi}}_{{\mathbf{q}}_{{\Vert}}}\left(\omega \right)}=\frac{{W}_{{\mathbf{q}}_{{\Vert}}}^{\text{GS}}}{1-{W}_{{\mathbf{q}}_{{\Vert}}}^{\text{GS}}{\Delta}{{\Pi}}_{{\mathbf{q}}_{{\Vert}}}\left(\omega \right)}.\end{equation} \tag{ 6 }$$

Here, ${W}_{{\mathbf{q}}_{{\Vert}}}^{\text{GS}}$ includes all ground state and environmental screening contributions,

$$\begin{align*}\hfill {\Delta}{{\Pi}}_{{\mathbf{q}}_{{\Vert}}}\left(\omega \right)& =\sum _{\alpha {\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}\frac{{\Delta}{f}_{{\mathbf{k}}_{{\Vert}}-{\mathbf{q}}_{{\Vert}}}^{\alpha }-{\Delta}{f}_{{\mathbf{k}}_{{\Vert}}}^{{\alpha }^{\prime }}}{\hslash \omega +\mathrm{i}{\gamma }_{\mathrm{T}}+{{\epsilon}}_{\alpha {\mathbf{k}}_{{\Vert}}-{\mathbf{q}}_{{\Vert}}}-{{\epsilon}}_{{\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}}\hfill \\ \hfill & \quad {\times}{\left\vert \langle \alpha {\mathbf{k}}_{{\Vert}}-{\mathbf{q}}_{{\Vert}}\vert {\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}\rangle \right\vert }^{2},\hfill \end{align*}$$

contains the optically induced deviations from the ground state polarization function only, and γ_T is the triplet dephasing rate, respectively. Hence, if the single particle dispersions, optical dipole and Coulomb matrix elements as well as the ground state screening of the Coulomb potential are known from ab initio calculations, the SBE provide a very efficient scheme to compute the optical response since they only need to be solved for a few bands and that part of k_∥ space in which populations and/or polarizations are optically induced.

Specifically, for TMDCs it is often sufficient to include the two spin–split valence and conduction bands near the K-points where the single particle dispersion displays a direct gap. Here, the single particle dispersion can be approximated by the relativistic dispersion

$$\begin{equation}{{\epsilon}}_{i\mathbf{K}+\mathbf{k}}^{c/v}={E}_{\mathrm{F},i}{\pm}\frac{1}{2}\sqrt{{{\Delta}}_{i}^{2}+{\left(2\hslash {v}_{\mathrm{F},i}k\right)}^{2}},\end{equation} \tag{ 7 }$$

that results from the widely used massive Dirac–Fermion (MDF) model Hamiltonian^?. Here, i = sτ combines the spin and valley index, Δ_i , v_F,i and E_F,i are the spin and valley dependent gap, Fermi velocity and Fermi level, respectively. The occurring parameters can be adjusted to reproduce the DFT bandstructure near the K-points. The overlap matrix elements resulting from the MDF Hamiltonian are

$$\begin{align*}\hfill \langle c\mathbf{K}+{\mathbf{k}}_{{\Vert}}\vert c\mathbf{K}+{{\mathbf{k}}_{{\Vert}}}^{\prime }\rangle & =\langle \nu \mathbf{K}+{{\mathbf{k}}_{{\Vert}}}^{\prime }\vert \nu \mathbf{K}+{\mathbf{k}}_{{\Vert}}\rangle \hfill \\ \hfill & ={u}_{{\mathbf{k}}_{{\Vert}}}{u}_{{{\mathbf{k}}_{{\Vert}}}^{\prime }}+{v}_{{\mathbf{k}}_{{\Vert}}}{v}_{{{\mathbf{k}}_{{\Vert}}}^{\prime }}\enspace {\mathrm{e}}^{-\mathrm{i}{\theta }_{{\mathbf{k}}_{{\Vert}}-{{\mathbf{k}}_{{\Vert}}}^{\prime }}}\hfill \\ \hfill \langle c\mathbf{K}+{\mathbf{k}}_{{\Vert}}\vert \nu \mathbf{K}+{{\mathbf{k}}_{{\Vert}}}^{\prime }\rangle & =\left({u}_{{\mathbf{k}}_{{\Vert}}}{v}_{{{\mathbf{k}}_{{\Vert}}}^{\prime }}\enspace {\mathrm{e}}^{-\mathrm{i}{\theta }_{{{\mathbf{k}}_{{\Vert}}}^{\prime }}}-{v}_{{\mathbf{k}}_{{\Vert}}}{u}_{{{\mathbf{k}}_{{\Vert}}}^{\prime }}\enspace {\mathrm{e}}^{-\mathrm{i}{\theta }_{{\mathbf{k}}_{{\Vert}}}}\right)\hfill \end{align*}$$

with ${u}_{{\mathbf{k}}_{{\Vert}}}^{2}=1-{v}_{{\mathbf{k}}_{{\Vert}}}^{2}=\frac{\sqrt{{{\Delta}}_{i}^{2}+{\left(2\hslash {v}_{\mathrm{F},i}k\right)}^{2}}+{{\Delta}}_{i}}{2\sqrt{{{\Delta}}_{i}^{2}+{\left(2\hslash {v}_{\mathrm{F},i}k\right)}^{2}}}$.

Although the material excitations in quasi-2D materials are confined to a region |z| ≲ d/2 ⩽ L/2, the Coulomb potential is long ranged and the interaction among particles confined in the layer is sensitive to the dielectric environment. Furthermore, for the evaluation of the opto-electronic system response, one needs the screened Coulomb potential that not only contains screening contributions from the dielectric environment but also all ground state contributions from the material itself.

Within DFT, the macroscopic (3D) dielectric constant is obtained as

$$\begin{equation*}{{\epsilon}}_{\mathrm{M}}^{-1}\left(\mathbf{q}\right)=\frac{{\langle {\tilde {V}}_{\mathbf{q}}\left(\mathbf{r}\right)\rangle }_{{\Omega}}}{{V}_{0}}={{\epsilon}}_{00}^{-1}\left(\mathbf{q}\right),\end{equation*}$$

where _q (r) is a lattice periodic function and ⟨...⟩_Ω denotes a spatial average over the elementary unit cell. If applied to an artificial 3D crystal consisting of parallel monolayers separated by large vacuum regions, the averaging over a region much larger than the extension of the electron density yields a dielectric function that decreases with the size of the artificial unit cell and hence, cannot be used to construct the screened 2D Coulomb potential [7]. To avoid this complication, we developed in a previous publication [6] an electrostatic model that is based on the bulk dielectric functions which can be directly computed from 3D DFT calculations, thus avoiding any difficulties related to the artificial insertion of large vacuum layers.

Whereas we refer to reference [6] for a detailed derivation of our electrostatic model, we briefly summarize the most important steps with respect to the screened interaction for further use in this paper (figure 1). As visualized in figure 1, we considered TMDCs embedded in different dielectric surroundings, such as in vacuum, on a substrate and embedded in two different dielectric media.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To account for the environmental and ground state screening, we consider a slab geometry with

$$\begin{equation*}{{\epsilon}}_{{\Vert}}\left(z\right)=\begin{cases}{{\epsilon}}_{{\Vert}}^{\mathrm{T}} \hfill & z{< }-L/2,\hfill \\ {{\epsilon}}_{{\Vert}} \hfill & \vert z\vert {< }L/2,\hfill \\ {{\epsilon}}_{{\Vert}}^{\mathrm{B}} \hfill & L/2{< }z\hfill \end{cases}\quad {{\epsilon}}_{\perp }\left(z\right)=\begin{cases}{{\epsilon}}_{\perp }^{\mathrm{T}} \hfill & z{< }-L/2,\hfill \\ {{\epsilon}}_{\perp } \hfill & \vert z\vert {< }L/2,\hfill \\ {{\epsilon}}_{\perp }^{\mathrm{B}} \hfill & L/2{< }z.\hfill \end{cases}\end{equation*}$$

and determine the interaction potential as solution of Poisson's equation. The resulting interaction potential within the slab is given by

$$\begin{align}\hfill {\phi }_{{\mathbf{q}}_{{\Vert}}}\left(z,{z}^{\prime }\right)& =\frac{2\pi }{\kappa q}\left({\mathrm{e}}^{-\tilde {q}\vert z-{z}^{\prime }\vert }\right.\hfill \\ \hfill & \quad +{c}^{+-}\enspace {\mathrm{e}}^{-\tilde {q}\left(L+z+{z}^{\prime }\right)}+{c}^{-+}\enspace {\mathrm{e}}^{-\tilde {q}\left(L-z-{z}^{\prime }\right)}\hfill \\ \hfill & \quad +\left.{c}^{--}\enspace {\mathrm{e}}^{-\tilde {q}\left(2L-z+{z}^{\prime }\right)}+{c}^{--}\enspace {\mathrm{e}}^{-\tilde {q}\left(2L+z-{z}^{\prime }\right)}\right)\hfill \\ \hfill & =\frac{2\pi }{\kappa q}\enspace {\mathrm{e}}^{-\tilde {q}\vert z-{z}^{\prime }\vert }+{\Delta}{\phi }_{{\mathbf{q}}_{{\Vert}}}\left(z,{z}^{\prime }\right)\hfill \end{align} \tag{ 8 }$$

with $\tilde {q}=\sqrt{\frac{{{\epsilon}}_{{\Vert}}}{{{\epsilon}}_{\perp }}}q$, $\kappa =\sqrt{{{\epsilon}}_{{\Vert}}{{\epsilon}}_{\perp }}$ and similar for κ_T/B and

$$\begin{equation*}{c}^{\eta {\eta }^{\prime }}=\frac{\left(\kappa +\eta {\kappa }_{\mathrm{T}}\right)\left(\kappa +{\eta }^{\prime }{\kappa }_{\mathrm{B}}\right)}{\left(\kappa +{\kappa }_{\mathrm{T}}\right)\left(\kappa +{\kappa }_{\mathrm{B}}\right)-\left(\kappa -{\kappa }_{\mathrm{T}}\right)\left(\kappa -{\kappa }_{\mathrm{B}}\right){\mathrm{e}}^{-2{\tilde {q}}_{{\Vert}}L}},\end{equation*}$$

from which the 2D Coulomb potential is obtained as ${\phi }_{{\mathbf{q}}_{{\Vert}}}^{2\text{D}}={\phi }_{{\mathbf{q}}_{{\Vert}}}\left(0,0\right)$. In the last line of equation (8), the first term describes the direct interaction of two electrons located at z and z' whereas the second term describes the interaction with image charges. From equation (8), the 'bare' Coulomb interaction ${V}_{{\mathbf{q}}_{{\Vert}}}\left(z,{z}^{\prime }\right)$, containing environmental screening only, is obtained by setting _∥ = _⊥ = 1 within the slab, whereas the true vacuum Coulomb interaction ${V}_{{\mathbf{q}}_{{\Vert}}}^{\text{vac}}\left(z,{z}^{\prime }\right)$ is obtained by additionally setting κ_T = κ_B = 1. Hence, equation (8) suggests that the Coulomb potential without intrinsic screening can be expressed as ${V}_{{\mathbf{q}}_{{\Vert}}}\left(z,{z}^{\prime }\right)={V}_{{\mathbf{q}}_{{\Vert}}}^{\text{vac}}\left(z,{z}^{\prime }\right)+{\Delta}{V}_{{\mathbf{q}}_{{\Vert}}}\left(z,{z}^{\prime }\right)$, where the last term describes the interaction via image charges. Similarly, one obtains the screened Coulomb interaction in bulk, ${W}_{{\mathbf{q}}_{{\Vert}}}^{\text{bulk}}\left(z,{z}^{\prime }\right)$, by inserting the bulk parameters for the in- and out-of-plane dielectric constants ${{\epsilon}}_{{\Vert}}^{\mathrm{B}}$ and ${{\epsilon}}_{\perp }^{\mathrm{B}}$ both within the slab and for the top and bottom environment.

Whereas both ${V}_{{\mathbf{q}}_{{\Vert}}}^{\text{vac}}\left(z,{z}^{\prime }\right)$ and ${W}_{{\mathbf{q}}_{{\Vert}}}^{\text{bulk}}\left(z,{z}^{\prime }\right)$ are independent of the slab thickness L, the interaction with image charges in the presence of environmental screening introduces a dependence on the slab thickness L both for the 'bare' and screened monolayer interaction potential. To model the monolayer potential, we shall assume a slab thickness L = D where D = c/2 is the natural layer-to-layer distance in the naturally stacked bulk crystal. In the long wavelength limitq_∥ D → 0 one obtains ${W}_{{\mathbf{q}}_{{\Vert}}}^{\text{GS,}2\text{D}}=\frac{2\pi }{q\left(\left({\kappa }_{\mathrm{T}}+{\kappa }_{\mathrm{B}}\right)/2+{r}_{0}q\right)}$, corresponding to the widely used Rytova–Keldysh potential with a screening length ${r}_{0}=D\left(2{\kappa }^{2}-{\kappa }_{\mathrm{T}}^{2}-{\kappa }_{\mathrm{B}}^{2}\right)/4{{\epsilon}}_{\perp }$ that depends on the dielectric contrast between the TMDC material and dielectric environment. In particular, if the dielectric environment is the bulk TMDC crystal itself, the screening length vanishes. For the short wavelengths q_∥ L/2 → ∞, the screened potential approaches $2\pi /\kappa q\enspace {\mathrm{e}}^{-\tilde {q}\vert z-{z}^{\prime }\vert }$, independently of the dielectric environment. Comparison with the bulk Coulomb potential allows us to identify ${{\epsilon}}_{{\Vert}}={{\epsilon}}_{{\Vert}}^{\mathrm{B}}$ and ${{\epsilon}}_{\perp }={{\epsilon}}_{\perp }^{\mathrm{B}}$ for the ground state screening even in the case of a monolayer. Hence, we can write for the Coulomb potential including environmental and ground state screening

$$\begin{equation}{W}_{{\mathbf{q}}_{{\Vert}}}^{\text{GS}}\left(z,{z}^{\prime }\right)={W}_{{\mathbf{q}}_{{\Vert}}}^{\text{bulk}}\left(z,{z}^{\prime }\right)+{\Delta}{W}_{{\mathbf{q}}_{{\Vert}}}\left(z,{z}^{\prime }\right),\end{equation} \tag{ 9 }$$

where ${\Delta}{W}_{{\mathbf{q}}_{{\Vert}}}\left(z,{z}^{\prime }\right)$ describes the screened interaction via image charges. Equation (9) should be compared with the division W = W^ML + ΔW that has been used in the GΔW approach where ΔW contains the changes in the screened potential as compared to the suspended monolayer [7].

Even monolayer TMDCs have an intrinsic thickness since they consist of three atomic layers and their atomic orbitals have a finite extension perpendicular to the layer plane. Consequently, the Coulomb interaction in these materials differs both from the exact 2D and the three-dimensional cases. Taking the in-plane periodicity and the finite out-of-plane extension into account, a 2D ansatz gives the Bloch states,

$$\begin{equation}{\phi }_{\alpha {\mathbf{k}}_{{\Vert}}}\left(\mathbf{r}\right)=\frac{1}{\sqrt{A}}\enspace {\mathrm{e}}^{\mathrm{i}{\mathbf{k}}_{{\Vert}}\cdot \mathbf{r}}{u}_{\alpha {\mathbf{k}}_{{\Vert}}}\left(\mathbf{r}\right)\end{equation} \tag{ 10 }$$

with strictly 2D crystal momenta but 3D spatial coordinates r. Hence, the Coulomb Hamiltonian contains the quasi-2D matrix elements

$$\begin{align*}\hfill {V}_{{\mathbf{q}}_{{\Vert}};{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}}}^{\alpha \beta {\beta }^{\prime }{\alpha }^{\prime }}& ={\int }_{ec}{\mathrm{d}}^{3}r{\int }_{ec}{\mathrm{d}}^{3}{r}^{\prime }\enspace {u}_{\alpha {\mathbf{k}}_{{\Vert}}-{\mathbf{q}}_{{\Vert}}}^{{\ast}}\left(\mathbf{r}\right){u}_{\beta {{\mathbf{k}}_{{\Vert}}}^{\prime }+{\mathbf{q}}_{{\Vert}}}^{{\ast}}\left({\mathbf{r}}^{\prime }\right)\hfill \\ \hfill & \quad {\times}{V}_{{\mathbf{q}}_{{\Vert}}}\left(z,{z}^{\prime }\right){u}_{{\beta }^{\prime }{{\mathbf{k}}_{{\Vert}}}^{\prime }}\left({\mathbf{r}}^{\prime }\right){u}_{{\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}\left(\mathbf{r}\right),\hfill \end{align*}$$

and a similar expression holds for the matrix elements of the screened Coulomb potential.

According to equation (8), the Coulomb interaction differs from the exact two-dimensional potential only in the exponential terms. Hence, we can write for the screened Coulomb potential

$$\begin{equation}{W}_{{\mathbf{q}}_{{\Vert}};{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}}}^{\alpha \beta {\beta }^{\prime }{\alpha }^{\prime }}={W}_{{\mathbf{q}}_{{\Vert}}}^{\mathrm{2}\mathrm{D}}{F}_{{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}}}^{\alpha \beta {\beta }^{\prime }{\alpha }^{\prime }}\left({\mathbf{q}}_{{\Vert}}\right)\end{equation} \tag{ 11 }$$

$$\begin{equation}\enspace \quad -{\Delta}{W}_{{\mathbf{q}}_{{\Vert}}}^{\mathrm{2}\mathrm{D}}{\Delta}{F}_{{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}}}^{\alpha \beta {\beta }^{\prime }{\alpha }^{\prime }}\left({\mathbf{q}}_{{\Vert}}\right).\end{equation} \tag{ 12 }$$

Here,

$$\begin{align}\hfill {F}_{{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}}}^{\alpha \beta {\beta }^{\prime }{\alpha }^{\prime }}\left({\mathbf{q}}_{{\Vert}}\right)& ={\int }_{ec}{\mathrm{d}}^{3}r{\int }_{ec}{\mathrm{d}}^{3}{r}^{\prime }\enspace {u}_{\alpha {\mathbf{k}}_{{\Vert}}-{\mathbf{q}}_{{\Vert}}}^{{\ast}}\left(\mathbf{r}\right){u}_{\beta {{\mathbf{k}}_{{\Vert}}}^{\prime }+{\mathbf{q}}_{{\Vert}}}^{{\ast}}\left({\mathbf{r}}^{\prime }\right)\hfill \\ \hfill & \quad {\times}{\mathrm{e}}^{-{\tilde {q}}_{{\Vert}}\vert z-{z}^{\prime }\vert }{u}_{{\beta }^{\prime }{{\mathbf{k}}_{{\Vert}}}^{\prime }}\left({\mathbf{r}}^{\prime }\right){u}_{{\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}\left(\mathbf{r}\right)\hfill \end{align} \tag{ 13 }$$

is a form factor and

$$\begin{equation}{\Delta}{F}_{{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}}}^{\alpha \beta {\beta }^{\prime }{\alpha }^{\prime }}\left({\mathbf{q}}_{{\Vert}}\right)={F}_{{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}}}^{\alpha \beta {\beta }^{\prime }{\alpha }^{\prime }}\left({\mathbf{q}}_{{\Vert}}\right)-{f}_{{\mathbf{k}}_{{\Vert}}}^{\alpha {\alpha }^{\prime }}\left({\mathbf{q}}_{{\Vert}}\right){f}_{{{\mathbf{k}}_{{\Vert}}}^{\prime }}^{\beta {\beta }^{\prime }}\left({\mathbf{q}}_{{\Vert}}\right)\end{equation} \tag{ 14 }$$

$$\begin{equation}{f}_{{\mathbf{k}}_{{\Vert}}}^{\alpha {\alpha }^{\prime }}\left({\mathbf{q}}_{{\Vert}}\right)={\int }_{ec}{\mathrm{d}}^{3}r\enspace {u}_{\alpha {\mathbf{k}}_{{\Vert}}-{\mathbf{q}}_{{\Vert}}}^{{\ast}}\left(\mathbf{r}\right){\mathrm{e}}^{-{\tilde {q}}_{{\Vert}}z}{u}_{{\alpha }^{\prime }{\mathbf{k}}_{{\Vert}}}\left(\mathbf{r}\right)\end{equation} \tag{ 15 }$$

is its correlated part, respectively. In the 2D limit, the exponential term in the form factor approaches unity and the 2D Coulomb matrix element is recovered due to the wavefunctions orthonormality. Furthermore, ${F}_{{{\mathbf{k}}_{{\Vert}}}^{\prime };{\mathbf{k}}_{{\Vert}}}^{\alpha \beta {\beta }^{\prime }{\alpha }^{\prime }}\left({\mathbf{q}}_{{\Vert}}=0\right)={\delta }_{\alpha {\alpha }^{\prime }}{\delta }_{\beta {\beta }^{\prime }}={f}_{{\mathbf{k}}_{{\Vert}}}^{\alpha {\alpha }^{\prime }}\left({\mathbf{q}}_{{\Vert}}=0\right){f}_{{{\mathbf{k}}_{{\Vert}}}^{\prime }}^{\beta {\beta }^{\prime }}\left({\mathbf{q}}_{{\Vert}}=0\right)$, while for large scattering vectors q_∥ → ∞ the form factor approaches 0. Since the potential of the image charges is particularly important in the long wavelength limit where the correlated part of the form factor is negligible, the Coulomb matrix elements of the quasi-2D TMDC structures can be expressed in a good approximation by the exact 2D term modified by the form factor.

In the following, we will assume that the single-particle band dispersion and wave functions can reliably be computed from DFT and use the DFT parameters as input for the SBE/DBE evaluations. From a fundamental point, this is in fact a nontrivial assumption. As is well known, DFT is based on the assumption that the two-particle interaction in any given system is a universal functional of the electron density while all system specific contributions to the many-body Hamiltonian are contained in the external potential. Using the dielectric model presented in section 2.2, it becomes clear that the substrate not only changes the screened Coulomb potential but also modifies the 'vacuum' potential that now contains interactions with image charges. Thus, in the presence of external screening, all particles (electrons and ions) interact via a modified Coulomb potential and hence, the electron–electron interaction can no longer be considered as universal.

However, analytical estimates based on the expression for the 'bare' Coulomb potential presented in section 2.2 show that the influence of the dielectric screening on the DFT single particle energies should be small for particles confined to a region |z|≲ d/2 < D/2. Indeed, rigorous treatments of monolayer TMDCs embedded in different dielectric environments have shown that the single-particle bandstructure in the proximity of the direct band gap at the K point remains unchanged for different environments [10]. Therefore, the single-particle bandstructure near the K-point can be obtained from an artificial 3D supercell calculation even in the presence of external screening.

The finite extension of the electron density in the direction perpendicular to the TMDC plane is illustrated in figure 2(a) for the example of MoS₂. Here, the normalized density distribution for both valence and conduction band at the K-point is plotted for different cuts through the unit cell.

Zoom In Zoom Out Reset image size

For MoS₂, the wave functions are linear combinations of the transition metal's d-orbitals that are strongly localized in between the chalcogenite layers. While the conduction-band wave function is dominated by the d_z -orbital, the wave function of the valence band is a linear combination of the in-plane d_xy - and ${d}_{{x}^{2}{y}^{2}}$-orbitals [19]. This is different in tungsten based materials where the transition metal's s-orbital and the chalcogen atom's p_z -orbital both have a finite contribution to the conduction or valence band respectively. These features result in a relatively higher density distribution in the tungsten layer for the conduction band whereas for the valence band the density distribution extends farther in the direction towards the chalcogen atoms, respectively.

For the optical response, the properties of the form factor are of interest especially close to the direct band. Therefore, we compute the form factor for scattering vectors $\vert {\mathbf{k}}_{{\Vert}}-{\mathbf{k}}_{{\Vert}}^{\prime }\vert$ with k_∥ fixed at the K-point and ${\mathbf{k}}_{{\Vert}}^{\prime }$ modified along different paths through the first Brillouin zone. The band structure is similar for all paths close to the K-point but changes significantly farther away, particularly for the conduction band where in the K + 60° direction a secondary minimum can be seen that forms the indirect band gap in the reference bulk structure.

In Figure 3, we show a density plot of the form factor for inter- and intraband interaction in the region around the K-points for the example of MoS₂. Despite the anisotropy of the band structure, the MoS₂ form factor of the interband interactions F_{c v} decreases isotropically around the K-point up to a range of $\vert {\mathbf{k}}_{{\Vert}}-{\mathbf{k}}_{{\Vert}}^{\prime }\vert \approx$ 0.45 Å⁻¹. The same is true for tungsten based materials, but only up to $\vert {\mathbf{k}}_{{\Vert}}-{\mathbf{k}}_{{\Vert}}^{\prime }\vert \approx$ 0.30 Å⁻¹. Thereafter, the gradient in K-direction is smaller than in Γ-direction, whereby the absolute differences are less than 0.06. Path dependent differences are larger with regard to the intraband interactions. Here, the gradient is smaller in the direction towards the Γ-point for the valence and towards the M-points for the conduction band, respectively. However, but even here, the form factor is isotropic in a range of $\vert {\mathbf{k}}_{{\Vert}}-{\mathbf{k}}_{{\Vert}}^{\prime }\vert \approx$ 0.15 Å⁻¹ and absolute differences are always less than 0.12. Thus, the differences occurring for the intraband form factor depending on different orbital compositions of the states seem to balance out in the integration over both valence and conduction band states in the interband form factor. In the following, averages over calculations along different paths are shown and used in subsequent calculations.

In order to be useful, the introduced form factor should be insensitive to the periodic boundary conditions applied in the DFT code and it should be applicable not only to model the Coulomb interaction in monolayers but also in bulk structures. To check these properties, we calculated the form factor for different extensions of the vacuum region included in the supercell.

Our results show that as long as we properly relax the structure before calculating the form factor we obtain similar results independent of the used super cell size. The same is true even if we consider the transition to a quasi-bulk structure differing from the common crystalline structure of MoS₂ in stacking order (usually 2H-/3R-phases) and in the interlayer distance by about 0.6 Å. On this basis, we conclude that we can use the determined form factors to model a wide range of mono- and multi-layer systems in different dielectric environments.

To illustrate the effect of the form factor in real space, we compare in figure 4 the exact 2D and the quasi-2D Coulomb potentials for MoS₂. As we can see, the form factor effectively removes the Coulomb singularity at the origin but leaves the potential unchanged for larger distances, i.e., for r_∥ ⪆ 7 Å. the slope of the quasi-2D potential approaches that of the original 2D potential.

Zoom In Zoom Out Reset image size

In figure 5, we show examples of the momentum dependent form factors for inter- and intraband transitions that were numerically evaluated using the full wavefunctions. To simplify detailed calculations on the SBE/DBE level, it is useful to develop analytical approximations for the form factor. For this purpose, we replace the full wavefunctions by their respective values at the K-point and incorporate the orbital dependent overlap matrix elements that occur within the MDF model.

Zoom In Zoom Out Reset image size

As shown by the green curve in figure 5, this procedure captures the main features of the full form factors for small scattering vectors. However, for scattering vectors larger than 0.1 Å⁻¹ or 0.2 Å⁻¹ clear deviations are seen such that additional corrections are needed.

In some of our previous works (see e.g. references [6, 20]), we adopted a semi-empirical Ohno potential to account for the finite thickness of the TMDC layers. This approximation has also been used to describe the Coulomb potential of molecules and nanotubes (see references [21, 22]) and was successfully applied in calculations of optical properties of graphene and monolayer TMDCs.

By using the Ohno potential, one regularizes the singularity of the Coulomb potential by introducing an effective thickness parameter d, which in reciprocal space occurs in an additional prefactor ${\text{e}}^{-{\mathbf{q}}_{{\Vert}}d}$ of the Coulomb potential. Even though this ansatz reduces the Coulomb scattering for larger scattering vectors, the detailed analysis shows that it underestimates the form factor for small scattering vectors. Therefore, we additionally introduce an orbital independent exponential function of the form $F\left({\mathbf{q}}_{{\Vert}}\right)={\text{e}}^{-\left({\sigma }^{2}{{q}_{{\Vert}}}^{2}/2+{q}_{{\Vert}}d\right)}$, while the orbital dependency is given by the MDF overlap matrix elements. With help of this function, the form factor can be described properly and physically correctly as it approximates 0 for large scattering vectors. Since the orbital dependence is contained in the MDF overlap matrix elements, the orbital independent part can be combined with the strictly 2D Coulomb potential to define the quasi-2D Coulomb potential ${V}_{{\mathbf{q}}_{{\Vert}}}=F\left({\mathbf{q}}_{{\Vert}}\right){V}_{{\mathbf{q}}_{{\Vert}}}^{\mathrm{2}\mathrm{D}}$.

A physical interpretation of the quadratic contribution can be obtained by considering the real space representation of the quasi-2D Coulomb potential:

$$\begin{align}\hfill V\left({\mathbf{r}}_{{\Vert}}\right)& =\sum _{{\mathbf{q}}_{{\Vert}}}{\mathrm{e}}^{\mathrm{i}{\mathbf{q}}_{{\Vert}}\cdot {\mathbf{r}}_{{\Vert}}}F\left({\mathbf{q}}_{{\Vert}}\right){V}_{{\mathbf{q}}_{{\Vert}}}^{\mathrm{2}\mathrm{D}}\hfill \\ \hfill & =\frac{1}{2\pi {\sigma }^{2}}\int {\mathrm{d}}^{2}{r}_{{\Vert}}^{\prime }\enspace {\mathrm{e}}^{-\vert {\mathbf{r}}_{{\Vert}}-{\mathbf{r}}_{{\Vert}}^{\prime }{\vert }^{2}/2{\sigma }^{2}}\frac{1}{\sqrt{{r}_{{\Vert}}^{\prime 2}+{d}^{2}}}.\hfill \end{align} \tag{ 16 }$$

Hence, the linear term in the exponent of the form factor changes the pure 2D Coulomb potential into the Ohno potential with thickness d, whereas the quadratic term lead to an additional convolution with a Gaussian of with σ.

In figure 6, we show a comparison of the interband form factor for five commonly discussed examples of semiconducting TMDCs. We notice that for tungsten based materials the form factor decreases slightly faster than for molybdenum based materials. Even more significant modifications occur for different chalcogen atoms. The comparison in figure 6 reveals a more rapid form factor decay when changing from S to Se to Te demonstrating that the gradient becomes steeper with increasing atom size. Since the microscopic distance between the chalcogenide sheets increases with the size of the involved atoms, the systems become slightly more 3D which is reflected in the form factor being a measure for the influence of the finite thickness.

Zoom In Zoom Out Reset image size

As a first application, we calculate the renormalization of the direct K point band gap for different dielectric environments. For the unexcited monolayer, the renormalized band gap is given by

$$\begin{equation}{{\Delta}}_{\mathbf{K}}={{\Delta}}^{\text{DFT}}+\sum _{{\mathbf{q}}_{{\Vert}}}\left({W}_{{\mathbf{q}}_{{\Vert}};\mathbf{K};\mathbf{K}+{\mathbf{q}}_{{\Vert}}}^{\nu \nu \nu \nu }-{W}_{{\mathbf{q}}_{{\Vert}};\mathbf{K};\mathbf{K}+{\mathbf{q}}_{{\Vert}}}^{cc\nu \nu }\right)\end{equation} \tag{ 17 }$$

In figure 7, we show the dependency of the gap on the environmental screening for five widely studied semiconducting monolayers (upper part). Since we plot the gap against the inverse value of the substrate dielectric constant, the figure covers the whole range of {1, ∞} and {−∞, −10}. For all materials investigated, we find a similar gap reduction for increased environmental screening. In the lower part of figure 7, we compare the gap shift of the different materials relative to the respective suspended monolayer with and without the form factor included. We note that the shift is almost independent of the form factor. To understand this behavior, we write the equation for the renormalized gap as

$$\begin{equation}{E}_{\mathrm{G}}={{\Delta}}^{\text{DFT}}+\left({W}_{\mathbf{K};\mathbf{K}}^{\nu \nu \nu \nu }\left({\mathbf{r}}_{{\Vert}}=0\right)-{W}_{\mathbf{K};\mathbf{K}}^{cc\nu \nu }\left({\mathbf{r}}_{{\Vert}}=0\right)\right)\end{equation} \tag{ 18 }$$

and make use of equation (16).

Zoom In Zoom Out Reset image size

Thus, for a strictly 2D system, the gap renormalization is determined by the screened Coulomb matrix elements at r_∥ = 0 and z = z' = 0. However, the form factor replaces z = z' = d in the value of the screened Coulomb potential and additionally averages the Coulomb potential over a region around the origin using a Gaussian weight function. Both effects have a significant impact on the direct interaction contained in W^{2D bulk}(r_∥) = 1/κr_∥ that varies strongly in the region around the origin. In contrast, changes contained in ${\Delta}W\propto 1/\sqrt{{{r}_{{\Vert}}}^{2}+{L}^{2}}$ vary only weakly within the region around the origin and hence are less affected by the form factor.

Furthermore, a comparison of the relative changes for the different materials in the regime of positive dielectric constants shows that these are slightly smaller for Mo rather than for W based materials with maximum differences in the range of few meV. Again, the influence of the incorporated chalcogenide atoms is more significant. Here, the maximum relative differences increase from Te to Se to S based compounds by approximately 80 meV and 45 meV, respectively.

In table 1, we give an overview of the computed band gaps for a variety of TMDCs and different dielectric environments. Comparison with available experimental data shows excellent agreement for all three Mo-based material systems. E.g. for MoS₂ we find E_G = 2.025 eV for the direct gap at the K-points in bulk, whereas we find E_G = 2.473 eV for the freely suspended, E_G = 2.279 eV on top of a fused silicon substrate, and E_G = 2.160 eV for an hBN encapsulated monolayer. The gap on a metal substrate (e.g. gold) can be estimated from the limit 1/_S → 0⁻, yielding E_G = 1.944 eV. These values are in excellent agreement with available experimental data, listed also in table 1. In contrast, the gaps for the W-based materials appear to be 130–150 meV below reported experimental values. Comparison with the MoX₂ systems shows that the predicted gaps for the corresponding WX₂ and MoX₂ are very similar, whereas the experimentally reported WX₂ gaps are always above those of the corresponding MoX₂ systems under comparable conditions. Since the environmentally induced band gap changes are captured very well by our analysis, we speculate that the origin for the systematic underestimation of the gap in the W-based systems is most likely at the DFT level.

For the suspended monolayers, we can also compare our results with other ab initio methods. Reported quasi-particle gaps based on the GW/GdW approach vary over a wide range of up to 400 meV and we find our results at the lower end of this spectrum. Regarding the relative differences in the gaps of the W- and Mo-based monolayers, the LDA and PBE based GW results show the same trend as our results, namely almost equal gaps for the sulfides and the selenides, thus supporting the assumption that the systematic discrepancy for the WX₂ originates at the DFT level. Regarding the shifts of the band gap with increased dielectric screening, our results agree well with the available GW-based predictions.

As a second application, we calculate the band gap renormalization resulting from finite carrier densities in different TMDCs. In the presence of excited carriers, the renormalized band gap

$$\begin{align}\hfill {{\Delta}}_{\mathbf{K}}& ={{\Delta}}^{\text{DFT}}+\sum _{{\mathbf{q}}_{{\Vert}}}\left({W}_{{\mathbf{q}}_{{\Vert}};\mathbf{K};\mathbf{K}+{\mathbf{q}}_{{\Vert}}}^{\nu \nu \nu \nu }-{W}_{{\mathbf{q}}_{{\Vert}};\mathbf{K};\mathbf{K}+{\mathbf{q}}_{{\Vert}}}^{cc\nu \nu }\right)\hfill \\ \hfill & \quad {\times}\left({f}_{{\mathbf{k}}_{{\Vert}}-{\mathbf{q}}_{{\Vert}}}^{\nu }-{f}_{{\mathbf{k}}_{{\Vert}}-{\mathbf{q}}_{{\Vert}}}^{c}\right).\hfill \end{align} \tag{ 19 }$$

contains the dynamically screened Coulomb matrix elements ${W}_{{\mathbf{q}}_{{\Vert}}}\left(\hslash \omega -{\tilde {{\epsilon}}}_{c\mathbf{K}+{\mathbf{k}}_{{\Vert}}-{\mathbf{q}}_{{\Vert}}}+{\tilde {{\epsilon}}}_{\nu \mathbf{K}+{\mathbf{k}}_{{\Vert}}}\right)$, that, in turn, contain the renormalized single particle dispersions ${\tilde {{\epsilon}}}_{c/\nu {\mathbf{k}}_{{\Vert}}}$. These equations are solved self consistently together with equation (6), assuming thermal carrier distributions at 300 K within the renormalized bands. Here, we distinguish between optically excited carrier densities (equal numbers of electrons and holes) and carrier densities due to doping, where we restrict the analysis to electron doping only. In the presence of excited carriers, both, screening of the Coulomb potential by the excited carriers and phase space filling contribute to the conduction and valence band renormalization. In contrast, for the case of electron doping, the valence band renormalization is solely due to screening effects corresponding to the Coulomb hole, whereas the phase space filling effects additionally contribute to the conduction band renormalization.

The interplay between phase space filling and screening contributions sensitively depends on the employed screening model. In particular, the approximation of static screening leads to an overestimation of screening effects. To show the importance of dynamical screening, we compare the band gap of a SiO₂ supported MoS₂ monolayer in figure 8 for different carrier excitation conditions. Using the static approximation for the screened Coulomb interaction, the gap decreases very rapidly with increasing carrier density. Here, the decrease depends only little on the fact whether the carriers are generated by optical excitation or doping, indicating that the gap renormalization is completely dominated by the Coulomb hole. In contrast, for the case of dynamic screening, the gap decreases more slowly with increasing carrier density and is larger for the symmetrical, optically induced electron–hole densities than for electron densities alone.

Zoom In Zoom Out Reset image size

Using the dynamical screening calculations, we compare in figure 9 the density dependent gap for the five investigated monolayer materials assuming optically induced electron–hole densities (solid lines) and electron densities only (dashed lines). Whereas the overall behavior is similar in all materials, the effects are somewhat more pronounced in Mo than in based systems, also increasing from S to Te. Furthermore, the renormalization depends on the distribution of carrier densities, thus the band gap changes are weaker due to doping densities than due to optically excited carriers. This effect is more pronounced in than in Mo based materials.

Zoom In Zoom Out Reset image size

Finally, we compare the influence of excited carriers on the band gap of MoS₂ for different dielectric environments in figure 10. As we can see, the environmentally induced band gap offset vanishes rapidly with increasing carrier density, showing that density dependent screening effects dominate the band gap over dielectric screening effects for moderate to large carrier densities.

Zoom In Zoom Out Reset image size

As third and final application, we calculate the resonance energies of the A exciton series for different TMDCs in various dielectric environments. Similar to the band gap renormalization, we find that the exciton binding energy decreases due to enhanced screening of the Coulomb interaction in media with increased dielectric constant. For the 1s exciton, the decrease in the band gap renormalization and binding energy nearly cancel. As a net result, the 1s exciton shows only a small red shift with increased dielectric screening. In contrast, for the higher exciton states, the band gap renormalization dominates over the decreased binding, yielding a pronounced red shift of the excited states with enhanced dielectric screening.

As representative examples, we show in figure 11 the excitonic resonances generated by solving the Dirac–Wannier equation for MoS₂ and WS₂ on a quartz substrate and for an hBN encapsulated configuration. Here, the four energetically lowest resonances are shown together with the quasiparticle band gap. As mentioned earlier, the gap of the W-based systems is systematically underestimated such that we have to shift the DFT gap by 149 meV for WS₂ while keeping the dipole matrix elements constant BEFORE we compute the renormalizations and exciton binding energies. While the band gap decreases from the SiO₂ to the hBN encapsulated sample by about 110 meV, the 1s resonance energy decreases only by less than a fourth of this amount. Thus, the binding energy is reduced by about 90–100 meV for WS₂ and MoS₂ respectively. In contrast, the binding energies of the 2s, 3s, 4s excitons are reduced by 40 − 20 meV only.

Zoom In Zoom Out Reset image size