Optical excitations in 2D semiconductors

Similar to the previously developed Hartree–Fock theory [18, 19], the Kohn–Sham scheme of DFT is an effective one-particle description in the mean field of all other particles. Hohenberg and Kohn have proven [20] that the ground-state energy and the wave function are related to the density $\rho(\vec{r}) = \sum_i^\text{occ} |\varphi_i(\vec{r})|^2$ by a unique functional. By varying them, the ground state density $\rho_0(\vec{r})$ can be found for which the functional derivative becomes zero. Eventually, the effective one-particle equation for the ith particle (wave function $\varphi_{i}(\vec{r})$ and energy ε_i ) is

$$\begin{equation} \left(\frac{\hat{\vec{p}}^{\,2}}{2 m} + \hat{V}\,^{\text{eff}}[\rho(\vec{r})]\right) \varphi_{i}(\vec{r}) = \epsilon_{i}\:\varphi_{i}(\vec{r}), \end{equation} \tag{ 1 }$$

which is the so-called Kohn–Sham equation [21]. The effective potential mediating the interaction to all other electrons is given by the Coulomb interaction to all electrons, the interaction with the external potential (of the nuclei), and the functional derivative of the exchange-correlation (xc) potential

$$\begin{equation} \hat{V}\,^{\text{eff}}[\rho(\vec{r})] = \frac{e^2}{4 \pi \varepsilon_0} \int d^3r^{^{\prime}} \frac{\rho(\vec{r}\,^{^{\prime}})}{|\vec{r}-\vec{r}\,^{^{\prime}}|} + {V}_{\text{ext}}(\vec{r}) + \frac{\delta E_{\text{xc}}[\rho(\vec{r})]}{\delta \rho(\vec{r})}. \end{equation} \tag{ 2 }$$

The solution has to be found self-consistently, since the input density is composed of the resulting wave functions. The quality of the formally exact DFT depends on the quality of the functional. For solid state systems, the local density approximation (LDA) [22, 23] or the generalized gradient approximation (GGA) [24, 25] are typically chosen since they systematically exploit the dependencies of $\rho(\vec{r})$ and $\nabla \rho(\vec{r})$. While ε_i are strictly the Lagrange multipliers, in practice the interpretation as single-particle excitation energies and band structures is often a reasonable first approximation. The latter can be improved (somewhat empirically) by mixing the non-local exchange potential with the local GGA potential leading to the so-called hybrid functionals [26]. For HSE06 [27], for example, we have observed an improvement of about 20% for the band gaps of 2D materials [3]. Unfortunately, HSE06 still underestimates the band gaps by more than 20%, and further systematic improvements of the functional are challenging. Another drawback is the lack of weak but long-range van der Waals-type interactions. Methods to correct for these interactions are summarized, e.g. in [28].

In general, DFT has been very successful in predicting structural and thermodynamic properties of materials while its description of excited states has not been satisfactory, and can most be taken as a starting-point for more advanced treatments such as MBPT. The calculations of optical properties from standard DFT, however, do not lead to reasonable predictions. To achieve this, the time-dependent generalization (TD-DFT) can be used [29–31], in which the adiabatic LDA (ALDA) is routinely applied for the exchange-correlation potential. Unfortunately, several shortcomings remain, especially for extended systems [16, 32]. These can be overcome using many-body perturbation techniques.

In many-body perturbation theory (MBPT), the discussion starts from the many-body ground state $| N \rangle$ [33]. Using the shorthand notation for position, spin, and time $(\vec{x}t) = (\vec{r},\sigma,t) = 1$ and atomic units (i.e. $\hbar = 1$, $\hbar^2/2 m = 1$, and $e^2 = 2$ leads to lengths in a_B and energies in Ry), the causal one-particle Green function is given by

$$\begin{equation} G(1,2) = -i \langle N | \hat{T}(\hat{\psi}(1)\hat{\psi}^\dagger(2)) | N \rangle. \end{equation} \tag{ 3 }$$

It describes the creation of an electron at 1, its propagation and annihilation at 2, or vice versa. $\hat{T}$ is the time-order operator which adds a negative sign when $t_1 \gt t_2$ and $\hat{\psi}^{(\dagger)}(1)$ are field operators for the annihilation (creation) of an electron that can be expressed as $\hat{\psi}^\dagger(\vec{x},t) = e^{i \hat{H} t} \sum_i \hat{a}^\dagger_i \varphi_i^*(\vec{x}) e^{-i \hat{H} t}$. $\hat{a}^\dagger_i$ denotes the creation of an electron in the ith state and ϕ_i the corresponding wave function.

The equation of motion for the Green function

$$\begin{align} \left(i\frac{\partial}{\partial t} + \nabla^2 - V_\text{ext}(\vec{r}) - V_\text{Coul}(\vec{r})\right) G(\vec{x}t,\vec{x}^{^{\prime}}t^{^{\prime}}) +i\int d\vec{x}^{^{\prime\prime}} \frac{e^2}{|\vec{r}-\vec{r}^{^{\prime\prime}}|} G_2(\vec{x}^{^{\prime\prime}}t,\vec{x}^{^{\prime\prime}}t,\vec{x}t,\vec{x}^{^{\prime}}t^{^{\prime}}) = \delta(\vec{x},\vec{x}^{^{\prime}}) \delta(t,t^{^{\prime}}), \end{align} \tag{ 4 }$$

where V_Coul is the Hartree potential, links the one-particle Green function to the two-particle Green function G₂, and we eventually find an infinite hierarchy of equations. The hierarchy can be broken by introducing the self-energy $\Sigma$

$$\begin{align} \left(i\frac{\partial}{\partial t} + \nabla^2 - V_\text{ext}(\vec{r}) - V_\text{Coul}(\vec{r})\right) G(\vec{x}t,\vec{x}^{^{\prime}}t^{^{\prime}}) -\int d\vec{x}^{^{\prime\prime}}dt^{^{\prime\prime}} \Sigma(\vec{x}t,\vec{x}^{^{\prime\prime}}t^{^{\prime\prime}}) G(\vec{x}^{^{\prime\prime}}t^{^{\prime\prime}},\vec{x}^{^{\prime}}t^{^{\prime}}) = \delta(\vec{x},\vec{x}^{^{\prime}}) \delta(t,t^{^{\prime}}) \,. \end{align} \tag{ 5 }$$

In principle, the self-energy can be obtained self-consistently by five equations, the so-called 'Hedin's equations' [33]:

$$\begin{equation} G(12) = G^0(12) + \int d(34) G^0(13) \Sigma(13) G(42) \end{equation} \tag{ 6 }$$

$$\begin{equation} \Sigma(12) = i \int d(34) G(14) W(1^+3) \Gamma(423) \end{equation} \tag{ 7 }$$

$$\begin{equation} W(12) = v(12) + \int d(34) v(42)P(34)W(13) \end{equation} \tag{ 8 }$$

$$\begin{equation} P(12) = -i \int d(34) G(23) G(42) \Gamma(341) \end{equation} \tag{ 9 }$$

$$\begin{equation} \Gamma(123) = \delta(12)\delta(13) + \int d(4567) \frac{\delta \Sigma(12)}{\delta G(45)} G(46) G(75) \Gamma(673). \end{equation} \tag{ 10 }$$

Equation (6) is the Dyson equation that links the Green function to the non-interacting G⁰ and the self-energy Σ. The latter is determined (equation (7)) by the screened Coulomb interaction W and the vertex function Γ ($1^+$ means the evaluation at $t_1+0^+$). W can be calculated from the bare Coulomb interaction $v(12) = \frac{\delta(t_1-t_2)}{|\vec{r}_1-\vec{r}_2|}$ and an additional term including the polarizability P. P in turn depends on G and Γ. The most complicated form is obtained for the vertex function (10), which also includes a functional derivative of the self-energy Σ.

Although Hedin's equations provide a formally exact framework for obtaining Σ, practical calculations require approximations. Often, the GW approximation is applied. This approximation starts from $\Sigma = 0$ in equations (6) and (10) and the self-energy simplifies to

$$\begin{equation} \Sigma(12) = i G(12) W(1^+2) \,. \end{equation} \tag{ 11 }$$

W can be calculated via the dielectric function in the random phase approximation $W(12) =$ $\int d(3) v(13) \varepsilon^{-1}(32)$, with $\varepsilon(12) = \delta(12) - \int d(3) v(13)P(32)$. For the frequency integration (after a Fourier transformation from the time domain), a simple plasmon-pole model [34, 35] is used. Finally, the wave function $\psi_{m\vec{k}}^\text{QP}(\vec{x})$ and energy $E_{m\vec{k}}^\text{QP}$ of the quasi-particle (m includes its spin) can be calculated by

$$\begin{equation} \left( -\nabla^2 + V_\text{ext}(\vec{r}) + V_\text{Coul}(\vec{r}) \right) \psi_{m\vec{k}}^\text{QP}(\vec{x}) + \int d\vec{x}^{^{\prime}} \Sigma(\vec{x},\vec{x^{^{\prime}}},E_{m\vec{k}}^\text{QP}) \psi_{m\vec{k}}^\text{QP}(\vec{x}^{^{\prime}}) = E_{m\vec{k}}^\text{QP} \psi_{m\vec{k}}^\text{QP}(\vec{x}) \,. \end{equation} \tag{ 12 }$$

Since Σ depends on $\psi_{m\vec{k}}^\text{QP}(\vec{x})$ and $E_{m\vec{k}}^\text{QP}$, solving equation (12) would require another self-consistency. However, due to the simplified vertex in the GW approximation self-consistency may actually worsen the results [36, 37], so we stick to updating the wave functions and energies based on DFT. The energy dependence of the self-energy is approximated as

$$\begin{equation} \Sigma_{mm^{^{\prime}}}(E) = \Sigma_{mm^{^{\prime}}}(\bar{E}) + \Sigma^{^{\prime}}_{mm^{^{\prime}}}(\bar{E}) \cdot (E - \bar{E}) \end{equation} \tag{ 13 }$$

and $Z_{m\vec{k}} = 1/\lbrace 1-\Sigma^{^{\prime}}_{mm}(\vec{k})\rbrace$ denotes the renormalization constant including the derivative of the self-energy with respect to the energy. $Z_{m\vec{k}}$ describes the spectral weight of the corresponding quasi-particle and its value is typically around 0.7 [38]. Transformed to $\vec{k}$ space, we can finally calculate

$$\begin{equation} \tilde{H}_{nn^{^{\prime}}}(\vec{k},E) = E^\text{DFT}_{n\vec{k}}\delta_{nn^{^{\prime}}} + (\Sigma(E)-V_\text{xc})_{nn^{^{\prime}}}(\vec{k}) \end{equation} \tag{ 14 }$$

and diagonalize $\tilde{H}(\vec{k},E^\text{QP}_{m\vec{k}}) D_{m\vec{k}} = E^\text{QP}_{m\vec{k}} D_{m\vec{k}}$. The solution is found when $E^\text{QP}_{m\vec{k}}$ agree on the left and right sides of this equation and lets us identify the corresponding wave function $\Psi_{m{\vec{k}}}^\text{QP} = \sum_n D^n_{m{\vec{k}}} \Psi_{n{\vec{k}}}^\text{DFT}$. While this treatment is important in many heterogeneous systems [39, 40], in other cases non-diagonal elements are less crucial and the solution simplifies to $E_{m\vec{k}}^\text{QP} = E_{m\vec{k}}^\text{DFT} + \langle \psi_{m\vec{k}}^\text{DFT} | \Sigma(E_{m\vec{k}}^\text{QP}) - V_\text{xc} | \psi_{m\vec{k}}^\text{DFT} \rangle$ with $\psi_{m\vec{k}}^\text{QP} = \psi_{m\vec{k}}^\text{DFT}$. In this case the linearized form is given as $E_{m\vec{k}}^\text{QP} = E_{m\vec{k}}^\text{DFT} + Z_{m\vec{k}} \langle \psi_{m\vec{k}}^\text{DFT} | \Sigma(E_{m\vec{k}}^\text{DFT}) - V_\text{xc} | \psi_{m\vec{k}}^\text{DFT} \rangle$. We note that the vertex corrections in general have relatively small effect on the QP band gap while they tend to up-shift the absolute band energies (relative to the average potential) [41, 42].

Equation (14) illustrates that only the difference of the self-energy and the exchange-correlation potential is important. Based on the observation that a fictitious metallic screening reproduces the LDA $V^\text{LDA}_\text{xc} \approx i G W_\text{metal}$ for many systems, an alternative approximation is possible [43, 44]

$$\begin{equation} \hat{H}^\text{QP} \approx \hat{H}^\text{LDA} + \overbrace{iG\underbrace{(W - W_\text{metal})}_{dW}}^{\Delta \Sigma}. \end{equation} \tag{ 15 }$$

This $GdW(\text{LDA})$ approximation allows to focus on the changed screening with respect to the LDA. Such changes are usually more than an order of magnitude smaller than the original quantities. This allows us to avoid the time-consuming RPA calculation of the dielectric functions and apply RPA-parametrised model functions [45]. In particular, we have found good agreement for 2D materials [44, 46], and we refer the reader to these works for further details.

Until now, we have focused on single excitations that can be described by the Green function $G\equiv G_1$. To understand the optical properties of materials, a detailed understanding of the correlated electron–hole excitations is necessary [16, 47], i.e. G₂ and its equation of motion, the Bethe–Salpeter equation (BSE) introduced by Strinati [48, 49]. Unlike the one-particle Green function, which mediates between the N and $(N\pm 1)$ particle system, G₂ contains electron–hole pairs (we neglect possible hole-hole or electron-electron pairs). The ground state $|N,0\rangle$ can therefore be excited by the creation and annihilation operators of G₂ to a state S, denoted $|N,S\rangle$, without changing the number of particles. The two-particle Green function is defined analogously to equation (3) as

$$\begin{equation} G_2(12;1^{^{\prime}}2^{^{\prime}}) = - \langle N | \hat{T}(\hat{\psi}(1)\hat{\psi}(2)\hat{\psi}^\dagger(2^{^{\prime}})\hat{\psi}^\dagger(1^{^{\prime}})) | N \rangle. \end{equation} \tag{ 16 }$$

For excitons, we can assume that an electron–hole pair is created at times t₁ and $t_1^{^{\prime}} = t_1+0^+$ and annihilated at t₂ and $t_2^{^{\prime}} = t_2+0^+$. We note that G₂ and the corresponding correlation function L are four-point functions

$$\begin{equation} L(12;1^{^{\prime}}2^{^{\prime}}) = -G_2(12;1^{^{\prime}}2^{^{\prime}}) + G_1(1,1^{^{\prime}}) G_1(2,2^{^{\prime}}) . \end{equation} \tag{ 17 }$$

Its equation of motion, the BSE is given by

$$\begin{equation} L(12;1^{^{\prime}}2^{^{\prime}}) = L_0(12;1^{^{\prime}}2^{^{\prime}}) + \int d(3456) L_0(14;1^{^{\prime}}3) K(35;46) L(62;52^{^{\prime}}), \end{equation} \tag{ 18 }$$

where $L_0(12;1^{^{\prime}}2^{^{\prime}}) = G_1(1,2^{^{\prime}}) G_1(2,1^{^{\prime}})$ describes the non-interacting case and $K(35;46)$ is the interaction kernel (electron–hole interaction). After the transformation to the energy domain and to two effective particles, the BSE can be rewritten as a generalized eigenvalue problem

$$\begin{align} (E_c-E_v)A_{vc}^S + \sum_{v^{^{\prime}}c^{^{\prime}}} K_{vc,v^{^{\prime}}c^{^{\prime}}}^{AA}(\Omega_S) A_{v^{^{\prime}}c^{^{\prime}}}^S + \sum_{v^{^{\prime}}c^{^{\prime}}} K_{vc,v^{^{\prime}}c^{^{\prime}}}^{AB}(\Omega_S) B_{v^{^{\prime}}c^{^{\prime}}}^S = \Omega_S A_{vc}^S \nonumber\\\sum_{v^{^{\prime}}c^{^{\prime}}} K_{vc,v^{^{\prime}}c^{^{\prime}}}^{BA}(\Omega_S) A_{v^{^{\prime}}c^{^{\prime}}}^S + (E_c-E_v) B_{vc}^S + \sum_{v^{^{\prime}}c^{^{\prime}}} K_{vc,v^{^{\prime}}c^{^{\prime}}}^{BB}(\Omega_S) B_{v^{^{\prime}}c^{^{\prime}}}^S = -\Omega_S B_{vc}^S . \end{align} \tag{ 19 }$$

$(E_c-E_v)$ are the energy differences of valence and conduction bands and K are the corresponding matrix elements of the interaction kernel

$$\begin{equation} K_{vc,v^{^{\prime}}c^{^{\prime}}}^{AA}(\Omega_S) = i \int d(3456) \psi_v(\vec{x}_4) \psi_c^*(\vec{x}_3) K(35,46;\Omega_S) \psi_{v^{^{\prime}}}^*(\vec{x}_5) \psi_{c^{^{\prime}}}(\vec{x}_6) \end{equation} \tag{ 20 }$$

(for K^AB , K^BA and K^BB the structure is similar, but the arguments are interchanged). Employing the GW approximation, the exchange part is the bare Coulomb interaction while the direct part is screened by the dielectric function $K(34;56) = K^x(35;46) + K^d(35;46) = -i\delta(3,4)\delta(5^-,6)v(3,6) + i\delta(3,6)\delta(4,5)W(3^+,4)$. By solving equation (19) we can find the exciton energies $\Omega_S$ and the amplitudes A^S and B^S of the excitations and de-excitations $\chi_S(\vec{x},\vec{x}^{^{\prime}}) = \sum_v^\text{occ} \sum_c^\text{emp} (A^S_{vc} \psi_c(\vec{x}) \psi^*_v(\vec{x}^{^{\prime}}) +B^S_{vc} \psi_v(\vec{x}) \psi^*_c(\vec{x}^{^{\prime}}))$.

Many studies have found that neglecting K^AB (i.e. setting $K^{AB} = 0 = K^{BA}$) leads to reasonable results for extended semiconductors; this procedure is known as the Tamm–Dancoff approximation [50]. Following this route leads to the more compact form

$$\begin{equation} (E_c-E_v)A^S_{vc} + \sum_{v^{^{\prime}}c^{^{\prime}}}K^{AA}_{vc,v^{^{\prime}}c^{^{\prime}}}(\Omega_S)A^S_{v^{^{\prime}}c^{^{\prime}}} = \Omega_S A^S_{vc} , \end{equation} \tag{ 21 }$$

where the dimension of the eigenproblem can still become large, given by the product of the hole and electron bands and the number of $\vec{k}$ points considered. Since spin–orbit coupling plays an important role in many 2D materials, the Hamiltonian does not have a block structure and does not decomposes into singlet and triplet states [16, 47].

The optical properties of a material are measured via the macroscopic dielectric function [16], which depends on the energy ω and the direction of the limit $\vec{q}\to 0$

$$\begin{equation} \varepsilon_M(\omega,\measuredangle\vec{q}) = \lim_{\vec{q}\to 0} \frac{1}{\varepsilon^{-1}_{\vec{G} = 0,\vec{G}^{^{\prime}} = 0}(\vec{q},\omega)}. \end{equation} \tag{ 22 }$$

For ω = 0 the dielectric constant ε₀ is obtained. Because the absorption spectrum is given by the imaginary part ($\varepsilon_M = \varepsilon_1 + i\varepsilon_2$), we will focus on $\varepsilon_2(\omega)$ in this work. To evaluate the coupling to the light field, we calculate the imaginary part

$$\begin{equation} \varepsilon_2(\omega) = \frac{16\pi e^2}{\omega^2} \sum_S \left|\frac{\vec{A}}{|\vec{A}|} \langle 0 | \vec{v} | S \rangle\right|^2 \delta(\omega-\Omega_S) \end{equation} \tag{ 23 }$$

as the sum of all possible excitonic states whose weights are determined by the normalized polarization of the light $\vec{A}/|\vec{A}|$ (we assume a linear polarization) and the velocity matrix elements. These elements can be calculated as $\langle 0 | \vec{v} | S \rangle = \sum_v^\text{occ} \sum_c^\text{emp} A^S_{vc} \langle v | \vec{v} | c \rangle$, therefore both required values (the excitation energies $\Omega_S$ and the electron–hole amplitudes $A^S_{vc}$) are known after diagonalizing the BSE Hamiltonian.

After the theoretical derivations in the last sections, the actual calculations have to be performed numerically on a computer.

The first approximation that is typically applied is the introduction of pseudopotentials (PS) [51–54] and corresponding wave functions. The basic idea is to add the inner electrons to the potential of the nuclei (because they are not relevant for chemical bonds etc) and treat only the valence electrons explicitly. By this procedure, equation (1) can be rewritten as

$$\begin{equation} \left(\frac{\hat{\vec{p}}^{\,2}}{2 m} + \hat{V}^{\text{eff,PS}}[\rho(\vec{r})]\right) \varphi^\text{PS}_{i}(\vec{r}) = \epsilon^\text{PS}_{i}\:\varphi^\text{PS}_{i}(\vec{r}) . \end{equation} \tag{ 24 }$$

Moreover, this allows the spin–orbit coupling [55] to be taken into account without having to solve the Dirac equation for the entire system.

In addition, the wave functions have to be treated numerically. For most calculations we use atom-centered Gaussian orbitals. The advantage of Gaussian orbitals is their short-range behavior in real and reciprocal space and the analytical solution of their integrals. Each function consists of several orbitals with decay constants γ (denoted by the index α including their s-, p-, d-character, etc) at $\vec{\tau}_\nu$

$$\begin{align} \varphi_{n,\vec{k}}(\vec{r})& = \sum_{\alpha, \nu} c_{n,\vec{k}}^{\alpha,\nu} \psi_{n,\vec{k}}^{\alpha,\nu}(\vec{r})\quad\text{with}\quad \nonumber\\ \psi_{n,\vec{k}}^{\alpha,\nu}(\vec{r})& = \frac{1}{N_0}\sum_{\vec{R}_j} \mathrm{e}^{i\vec{k}(\vec{r}_j + \vec{\tau}_\nu)} \phi^{\alpha,\nu}(\vec{r} - \vec{R}_j - \vec{\tau}_\nu). \end{align} \tag{ 25 }$$

For details and an efficient numerical implementation, we refer to [56, 57]. The coefficients $c_{n,\vec{k}}^{\alpha,\nu}$ are obtained by the solution of the eigenvalue problem

$$\begin{equation} \sum_{\alpha^{^{\prime}}, \nu^{^{\prime}}} \left\lbrace H_{\alpha,\nu}^{\alpha^{^{\prime}},\nu^{^{\prime}}}(\vec{k}) - E_{n}(\vec{k}) S_{\alpha,\nu}^{\alpha^{^{\prime}},\nu^{^{\prime}}}(\vec{k}) \right\rbrace c_{n}^{\alpha^{^{\prime}},\nu^{^{\prime}}}(\vec{k}) = 0 , \end{equation} \tag{ 26 }$$

involving the overlap matrix S in addition to the Hamilton matrix H.

The two-point quantities in equation (11) have to be calculated with an auxiliary basis β, e.g. $W_{\beta,\beta^{^{\prime}}}(\vec{q},\omega)$. The required three-center integrals

$$\begin{equation} M^{mn}_{\beta}(\vec{k},\vec{q}) = \int d\vec{r} \: \varphi^*_{m\vec{k}}(\vec{r}) \chi^*_{\beta}(\vec{q},\vec{r}) \varphi_{n,\vec{k}+\vec{q}}(\vec{r}) \end{equation} \tag{ 27 }$$

use the Gaussian orbitals $\varphi_{n,\vec{k}}$ from DFT and again Gaussian orbitals or plane waves for $\chi_{\beta}(\vec{q},\vec{r})$ [58]. In particular, the GdW allows for a small number of auxiliary basis functions (plane waves) and enables many of the calculations discussed in the following. For example, periodic supercell calculations up to about 100 atoms are possible. Another critical aspect is the convergence of the supercell size when employing the GW/GdW approximation in a 2D material (in a 3D periodic code). Due to the long-range character, an inversely proportional extrapolation to infinite cell sizes is required; for details and the discussion of convergence see [44].

In addition, we use GPAW [59, 60] for several calculations. This code is based on the Projector Augmented Wave (PAW) method [61] which is as a generalized pseudopotential method. Already at the DFT level, we employ a plane wave basis. For further details we refer to [3].

The theoretical approach described in section 2.2 has been developed for neutral semiconductors. In the case of doping, this is no longer strictly correct. The electron and the corresponding hole created by the excitation with light can bind to a third particle (an additional electron or hole in n- and p-doped systems, respectively).

The treatment of three-particle states requires the description of three independent quasi-particles. The description as an exciton with a 'small' perturbation is often insufficient because both electrons or holes are indistinguishable. In principle, this would lead to the three-particle Green function G₃ and the solution of its equation of motion. To the best of our knowledge, the complete formal derivation has not been presented in the literature and seems to be out of reach. Here we will approximate the three-particle Hamiltonian similarly to the BSE. This approach has been adopted by Tempelaar and Berkelbach [62], Torche and Bester [63], and Zhumagulov et al [64].

In second quantization the many-body Hamiltonian is given as

$$\begin{equation} \hat{H} = \sum_{ij} h_{ij} \hat{a}^\dagger_{i} \hat{a}_{j} + \frac{1}{2} \sum_{ijmn} V_{ij,mn}~ \hat{a}^\dagger_{i} \hat{a}^\dagger_{j} \hat{a}_{n} \hat{a}_{m} \, , \end{equation} \tag{ 28 }$$

where the matrix elements of the single-particle terms (kinetic energy and external potential) are h_ij and the matrix elements of the Coulomb interaction are $V_{ij,mn}$. In a periodic system each index has to be understood including the momentum $i = (i,{\vec{k}}_i)$. Within Hartree–Fock theory, the single-particle states $\phi_{n}(x)$ are the solutions of the Hartree–Fock equations, and the matrix elements can be calculated by

$$\begin{equation} V_{ij,mn} = \int \phi_{i}^\ast(x) \phi_{j}^\ast(y) V(x,y) \phi_{m}(x) \phi_{n}(y) dx \, dy , \end{equation} \tag{ 29 }$$

where $x = (\vec{r},\sigma)$ and the bare Coulomb interaction is $V(x,x^{^{\prime}}) = e^2/|\vec{r}-\vec{r}^{^{\prime}}|$. For an electron–hole pair configuration $|vc\rangle : = \hat{a}^\dagger_{c} \hat{a}_{v} |0\rangle$ the matrix element is

$$\begin{equation} \begin{aligned} \langle vc | \hat{H} | v^{^{\prime}}c^{^{\prime}} \rangle =(E_0^{(HF)}+\epsilon_c^{(HF)}-\epsilon_v^{(HF)})\delta_{cc^{\prime}}\delta_{vv^{\prime}} - (V_{v^{^{\prime}}c,vc^{^{\prime}}} - V_{v^{^{\prime}}c,c^{^{\prime}}v}) , \end{aligned} \end{equation} \tag{ 30 }$$

with $\epsilon^{(\text{HF})}_n$ being the band-structure energy and $E_0^{(\text{HF})}$ the ground-state energy. Note that equation (30) is not diagonal and the diagonalization of the Hamiltonian

$$\begin{equation} \sum_{v^{^{\prime}}c^{^{\prime}}} \langle vc | \hat{H} | v^{^{\prime}}c^{^{\prime}} \rangle A^{S}_{v^{^{\prime}}c^{^{\prime}}} = \Omega^{S} A^{S}_{vc} \end{equation} \tag{ 31 }$$

leads to the excitation energies $\Omega^{S}$ and the excited states, which are linear combinations of the free electron–hole pair configurations $|S\rangle = \sum_{vc} A^{S}_{vc} | vc \rangle$. In the last section we found a very similar form using many-body perturbation theory, i.e. after the Tamm–Dancoff approximation in equation (21). However, two important changes can be concluded: First, the energy levels from Hartree–Fock are replaced by those from GW and second, the direct interaction is screened by the dielectric function of the system $W(\vec{r},\vec{r}^{^{\prime}}) = \int \epsilon^{-1}(\vec{r},\vec{r}^{^{\prime\prime}}) V(\vec{r}^{^{\prime\prime}},\vec{r}^{^{\prime}}) d^3r^{^{\prime\prime}}$.

For the three-particle states, we take a similar route. For the configurations of two electrons and one hole $| vc_1c_2 \rangle : = \hat{a}^\dagger_{c_2} \hat{a}^\dagger_{c_1} \hat{a}_{v} |0\rangle$ ($c_1 \lt c_2$), we can derive the matrix elements within Hartree–Fock [65, 66]

$$\begin{align} \langle vc_1c_2 | \hat{H} | v^{^{\prime}}c^{^{\prime}}_1c^{^{\prime}}_2 \rangle& = \left(E_0^{(\text{HF})} + \epsilon^{(\text{HF})}_{c_1} + \epsilon^{(\text{HF})}_{c_2} - \epsilon^{(\text{HF})}_{v}\right) \delta_{c_1c_1^{^{\prime}}} \delta_{c_2c_2^{^{\prime}}} \delta_{vv^{^{\prime}}} - (V_{v^{^{\prime}}c_1,vc^{^{\prime}}_1} - V_{v^{^{\prime}}c_1,c^{^{\prime}}_1v}) \delta_{c_2,c_2^{^{\prime}}}\nonumber\\ & \quad - (V_{v^{^{\prime}}c_2,vc^{^{\prime}}_2} - V_{v^{^{\prime}}c_2,c^{^{\prime}}_2v}) \delta_{c_1,c_1^{^{\prime}}} + (V_{v^{^{\prime}}c_1,vc_2^{^{\prime}}} - V_{v^{^{\prime}}c_1,c_2^{^{\prime}}v}) \delta_{c_2,c_1^{^{\prime}}} + (V_{v^{^{\prime}}c_2,vc_1^{^{\prime}}} - V_{v^{^{\prime}}c_2,c_1^{^{\prime}}v}) \delta_{c_1,c_2^{^{\prime}}}\nonumber\\ & \quad+ (V_{c_1c_2,c^{^{\prime}}_1c^{^{\prime}}_2} - V_{c_1c_2,c^{^{\prime}}_2c^{^{\prime}}_1}) \delta_{v,v^{^{\prime}}} . \end{align} \tag{ 32 }$$

In addition to the independent motion of each of the three particles, the two body interactions between each pair of particles are included. In analogy to the BSE, we replace V → W in the direct term of the electron–hole interaction and for the interaction of both electrons (last terms). To ensure invariance when changing c₁ and c₂, both terms are screened. Finally, we find the matrix elements

$$\begin{align} \langle vc_1c_2 | \hat{H}^\text{trion} | v^{^{\prime}}c^{^{\prime}}_1c^{^{\prime}}_2 \rangle& = (\epsilon_{c_1} + \epsilon_{c_2} - \epsilon_{v}) \delta_{c_1c_1^{^{\prime}}} \delta_{c_2c_2^{^{\prime}}} \delta_{vv^{^{\prime}}} - (W_{v^{^{\prime}}c_1,vc^{^{\prime}}_1} - V_{v^{^{\prime}}c_1,c^{^{\prime}}_1v}) \delta_{c_2,c_2^{^{\prime}}}\nonumber\\ & \quad - (W_{v^{^{\prime}}c_2,vc^{^{\prime}}_2} - V_{v^{^{\prime}}c_2,c^{^{\prime}}_2v}) \delta_{c_1,c_1^{^{\prime}}} + (W_{v^{^{\prime}}c_1,vc_2^{^{\prime}}} - V_{v^{^{\prime}}c_1,c_2^{^{\prime}}v}) \delta_{c_2,c_1^{^{\prime}}} + (W_{v^{^{\prime}}c_2,vc_1^{^{\prime}}} - V_{v^{^{\prime}}c_2,c_1^{^{\prime}}v}) \delta_{c_1,c_2^{^{\prime}}}\nonumber\\ & \quad + (W_{c_1c_2,c^{^{\prime}}_1c^{^{\prime}}_2} - W_{c_1c_2,c^{^{\prime}}_2c^{^{\prime}}_1}) \delta_{v,v^{^{\prime}}} . \end{align} \tag{ 33 }$$

and by diagonalizing the corresponding Hamiltonian

$$\begin{equation} \sum_{v^{^{\prime}}c_1^{^{\prime}}c_2^{^{\prime}}} \langle vc_1c_2 | \hat{H}^\text{trion} | v^{^{\prime}}c^{^{\prime}}_1c^{^{\prime}}_2 \rangle A^{(T,\vec{K})}_{v^{^{\prime}}c^{^{\prime}}_1c_2^{^{\prime}}} = E^{(T,\vec{K})} A^{(T,\vec{K})}_{vc_1c_2} \end{equation} \tag{ 34 }$$

we can calculate the energy of the trions $E^{(T,\vec{K})}$ and their wave functions

$$\begin{equation} \Phi^{(T,\vec{K})}(x_h,x_1,x_2) = \frac{1}{\sqrt{2}} \sum_{vc_1c_2} A^{(T,\vec{K})}_{vc_1c_2} \phi_{v}^\ast(x_h) \left\lbrace \phi_{c_1}(x_1) \phi_{c_2}(x_2) - \phi_{c_2}(x_1) \phi_{c_1}(x_2) \right\rbrace \end{equation} \tag{ 35 }$$

for trions T with momentum $\vec{K}$. Because of momentum conservation, the wave numbers of $vc_1c_2$ have to satisfy $\vec{k}_v = \vec{k}_{c_1} + \vec{k}_{c_2} - \vec{K}$. Note that the trion in an optical spectrum is not found at the trion energy $E^{(T,\vec{K})}$, but at the transition energy

$$\begin{equation} \Omega( |T,\vec{K}\rangle\leftrightarrow |c\vec{K}\rangle ) = E^{(T,\vec{K})} - \epsilon_{c\vec{K}} \end{equation} \tag{ 36 }$$

as an electron remains in the conduction band $|c\vec{K}\rangle$ (or is already present before absorption). Note that the $\vec{k}$ points in equation (33) are responsible for both, the convergence as well as the applied doping. As discussed in previous works, our results should be understood in the limit of vanishing doping [67].

The general treatment of electro-magnetic fields in a periodic calculation is challenging. Corresponding to the electric and magnetic fields $\vec{E}$ and $\vec{B}$, the (vector) potentials ϕ and $\vec{A}$ can be determined by $\vec{E} = - \nabla \varphi - \frac{\partial \vec{A}}{\partial t}$ and $\vec{B} = \nabla \times \vec{A}$. These potentials ϕ and $\vec{A}$ enter the Hamiltonian, and if they are not only varying vertical to the 2D material, the lattice periodicity is destroyed. To overcome this in magnetic fields, Kohn [68] and Roth [69] have reformulated the theory around the 1960s, as did Chang and others later using Berry phase techniques [70–74]. Although several effects of a magnetic field on optical absorption have been measured and theoretically modeled over the years, a predictive theoretical method was lacking. In the following, we show a first step in this direction [75].

Building on the ideas of Kohn and Roth, we solve the Hamiltonian without a field $\hat{H}_0^\text{eff} | \Psi_{n\vec{k}}^0 \rangle = E_{n\vec{k}} |\Psi_{n\vec{k}}^0 \rangle$ (see section 2.1). In a next step, we include the magnetic field (along the out-of-plane z-axis) as a perturbation

$$\begin{align} \hat{H}^\text{eff}& = \hat{H}_0^\text{eff} + \frac{e}{2m_e}(\hat{L}_z + g_e \hat{S}_z) B_z + \frac{e^2}{2m_e}(x^2+y^2)B_z^2\nonumber\\& = : \hat{H}_0^\text{eff} - \hat{m}_z B_z + \mathcal{O}(B_z^2) . \end{align} \tag{ 37 }$$

Here $\hat{L}$ and $\hat{S}$ are the angular momentum and spin operators, respectively, where we assume that $\hat{L}$ can be evaluated at the position of our local basis functions. For small fields, we only need to consider the linear term in B (we omit z for brevity), and the magnetic moment is given by an orbital and a spin contribution

$$\begin{align} \hat{m}& = : \hat{m}^\text{orb} + \hat{m}^\text{spin} \nonumber\\ m_{n\vec{k}}& = \langle \Psi_{n\vec{k}}^0 | \hat{m}^\text{orb} | \Psi_{n\vec{k}}^0 \rangle + \langle \Psi_{n\vec{k}}^0 | \hat{m}^\text{spin} | \Psi_{n\vec{k}}^0 \rangle . \end{align} \tag{ 38 }$$

The destroyed lattice periodicity is reflected in the spatial dependency of the angular momentum operator $\hat{L}_z = x\hat{p}_y - y\hat{p}_x$ which prohibits a direct evaluation. Therefore, we follow the derivation of Chang and Niu [70] for the magnetic moment of a wave packet

$$\begin{equation} m_{n\vec{k}}^\text{orb} = \mu_\text{B} \text{Im} \left\langle \frac{\partial u_{n\vec{k}}}{\partial k_x} \middle| \hat{H}^\text{eff}(\vec{k}) - E_{n\vec{k}} \middle| \frac{\partial u_{n\vec{k}}}{\partial k_y} \right\rangle \,, \end{equation} \tag{ 39 }$$

where $u_{n\vec{k}}$ is the lattice periodic part of the wave function and its derivative is evaluated as finite difference. Since the spin part of equation (38) can be calculated straight forward, this allows the study of the changes of the electronic states for small fields.

To investigate the changes of the optical properties, we have to take into account the wave function of the excitons. This can be evaluated from the BSE (19). Consequently, the shift m^S of each exciton S can be evaluated by

$$\begin{align} m^{S}& = \sum_{vc\vec{k}} |A^S_{vc\vec{k}}|^2 (m_{c\vec{k}} - m_{v\vec{k}}) \quad \text{and} \nonumber\\ g^{S}& = -2 \sum_{vc\vec{k}} |A^S_{vc\vec{k}}|^2 (m_{c\vec{k}} - m_{v\vec{k}}) / \mu_\text{B} . \end{align} \tag{ 40 }$$

The g factor of the exciton S is often related to the difference of corresponding experimental measurements with $\sigma^{+/-}$ circularly polarized light $g \mu_B B : = \Omega_{\sigma^+} - \Omega_{\sigma^-}$ [76]. Because the effects are reversed for different helicities, the factor −2 in equation (40) is obtained. If excitonic effects were neglected, the g factor of the transition from ($v\vec{k}$) to ($c\vec{k}$) could be approximated by

$$\begin{equation} g^{(v\vec{k}\to c\vec{k})}_\text{band} = -2 (m_{c\vec{k}} - m_{v\vec{k}}) / \mu_\text{B}. \end{equation} \tag{ 41 }$$

The theory presented in this section focused on the case of weak magnetic fields, and will be applied in section 8. For stronger magnetic fields, the diamagnetic term and the formation of Landau levels, i.e. effects beyond linear response, become important.

Electronically the TMDCs are semiconductors with indirect band gap in the bulk or for several layers, while they have a direct band gap in the monolayers. Figure 2 shows the band structures of MoS₂ and WSe₂ along a path through the high-symmetry points Γ, M, and K. The dashed lines show the results of DFT-LDA, that clearly underestimate the band gaps compared to our MBPT calculations within the GdW approximation. The latter reveal direct band gaps of 2.89 and 2.44 eV at the K point. Both values are in reasonable agreement to experimental measurements of 2.72 and 2.56 eV [79]. The large difference to the DFT band gaps of 1.78 and 1.65 eV (quasi-particle correction) is the result of the reduced screening in two-dimensional materials (in bulk, the direct gap in MoS₂ is 2.19 eV and the quasi-particle correction is less than half the size [80, 81]). The smaller the screening, the larger the inaccuracy of the DFT-LDA which is based on the homogeneous electron gas. In figure 2, the red and blue colors denote the spin polarization $\langle S_z\rangle$ and exhibit a clear out-of-plane polarization (i.e. orthogonal to the layers). The highest valence bands (VBs) and the lowest conduction bands (CBs) are polarized by more the 90% (often nearly 100%) at the band extrema. In particular, the VBs are distinctly split, which can be traced back to the strong spin–orbit coupling. Due to the heavier elements in WSe₂, the VBs are split by 570 meV, while the value for MoS₂ is only 180 meV. The CBs show an interesting characteristic close to the K point. They approach each other, and in the case of Mo-based TMDCs, they even change order. A corresponding schematic sketch is shown in the right panel of figure 2.

Zoom In Zoom Out Reset image size

The low-energy part of the optical spectrum is mostly governed by the transitions from the highest VBs to the lowest CBs. Due to the Coulomb interaction between the excited electron and the generated hole, peaks with strong optical weight can be found below the corresponding band gaps in figure 3. The lowest peaks are called A excitons and have, with respect to the band gaps, exciton binding energy of 0.76 eV and 0.68 eV, respectively. For bright excitons the spin of electron and hole have to be the same, i.e. in the absorption spectrum only transitions between bands with the same color in figure 2 can occur. Therefore, the next bright exciton B stems from the bands with different spin. The splitting between A and B follows the large spin-orbit splitting of the VBs and is thus stronger in WSe₂ compared to MoS₂. Furthermore, peaks with lower optical weight corresponding to excited states (A$_\text{2s}$, B$_\text{2s}$) are observed. Similar to possible excitations in a hydrogen atom, excitons can also form a Rydberg-like series [83, 84].

Zoom In Zoom Out Reset image size

A detailed understanding of the excitations is possible by analyzing their wave functions. In figure 4(a), the A exciton in MoS₂ is shown in real and reciprocal space. When the hole is fixed at the central Mo atom, the probability to find an electron is significant in a circular region of about 20 Å in diameter. The Fourier transformation relates this to the reciprocal contributions which are significant around the ±K points. For WSe₂ the comparison between dark and bright excitons is shown in figure 4(b). While their spatial extent is similar (but not identical, not shown here), the corresponding conduction bands have different spin character, and the selection rules forbid that spin flip processes are bright. Due to the energetic position of the bands in WX₂, the appearance of D as the lowest state is not surprising. Interestingly, also in MoS₂ the transition to the second-lowest conduction band leads to a dark exciton D below A. We find that the reason for this is a combination of three tiny differences: (1) The effective mass of the bands is different, (2) the exchange interaction is absent for D, and (3) the slightly different exciton wave functions also lead to small differences of the direct interaction [82]. Experiments [85] have proven that dark excitons are the ground state for MoS₂, WS₂, and WSe₂. For MoSe₂ the optically active and dark states are close in energy. Further details of the calculated energetic alignment can be found in table 1 of [82].

Zoom In Zoom Out Reset image size

After the optically bright exciton and the spin-forbidden dark states have been introduced, we now discuss the characteristic of excitons whose momentum is much larger than that of a photon. Consequently, all transitions are dark in the absence of additional momentum transfer by phonons or the like.

In figure 5, we examine possible excitations for different momenta. We restrict ourselves to Q pointing in the y direction in figure 5(a). These momenta can be used in the general BSE (21), in which all conduction band properties and wave functions have to be evaluated at $\vec{k}+\vec{Q}$. For further details we refer to our work [86]. The resulting spectra for 16 different Q are shown in figure 5(b). While the lowest one has already been discussed in figure 3, we find that the intensities of the peaks for larger Q decrease rapidly before increasing again at about half of the probed window. Figure 5(c) is more useful to follow the energies of the excitons independently of their oscillator strength. Here the color denotes the absolute square of the dipole operator assuming complete momentum transfer Q (due to phonons or other processes)

$$\begin{equation} \langle 0 | e^{-i \vec{Q} \vec{r}} \hat{\vec{p}} | S,\vec{Q} \rangle = \sum_{vc\vec{k}} A^{(S,\vec{Q})}_{vc\vec{k}} \sum_\sigma \int \phi_{v,\vec{k}}^*(\vec{r}, \sigma) e^{-i \vec{Q} \vec{r}} \hat{\vec{p}} \phi_{c,\vec{k}+\vec{Q}}(\vec{r}, \sigma) d\vec{r}. \end{equation} \tag{ 42 }$$

As discussed above, we find a dark exciton below A (at Q = 0). At slightly higher energy, this characteristic is repeated for the B exciton. As Q increases, we find increasing exciton energies [87], while for even larger Q there are two more local minima at Q ≈ 0.65 Å⁻¹ and Q ≈ 1.3 Å⁻¹. These correspond to transitions between Λ or the high-symmetry points (K$\to\Lambda$ and K'→K), respectively. The red color (or the peaks in figure 5(b), respectively) can be interpreted physically as follows: If the system can provide the momentum Q during absorption, a peak in the optical spectrum with the calculated dipole strength may be observed. We note that subsequent works [88, 89] have found similar physical trends. In experiment electron energy-loss spectroscopy (EELS) or neutron scattering could be applied. Measurement for small momentum transfer have been reported by Hong et al [90].

Zoom In Zoom Out Reset image size

In experiments 2D materials are either deposited on a substrate or encapsulated e.g. in hBN. Even if the 2D compounds are chemically inert, the substrate can change several properties significantly. In the previous section, it was argued that the large exciton binding energies in freestanding monolayers are a consequence of quantum confinement and weak intrinsic screening. The weak intrinsic screening implies that there is plenty of room for enhancing the screening, and thereby influence the excitations, via the dielectric environment. Here we demonstrate this effect by means of GdW calculations [44, 91] for supported monolayer TMDCs.

Figure 6 compares the vacuum results discussed above with MoS₂ placed on the semiconductor SiO₂ and the metallic Au(111) surface. The enhanced screening significantly reduces the band gap E_g from 2.89 to 2.25 eV. At the same time, the exciton binding energies E_b shrink, leading to a much smaller decrease (red-shift) of the exciton energies $\Omega^\text{A}$ of only 50 meV compared to the changes of the band gap. The physical origin of the reduction is the image charge effect [92–94]. When an electron is removed from the system (valence band), the remaining hole gets screened. The screening reduces the energy cost of removing the electron and thus the ionisation potential decreases, i.e. the valence band moves up. Similarly, the electron affinity decreases due to the screening, i.e. the conduction band moves down and consequently the gap reduces. However, the polarization similarly reduces the electron–hole binding, resulting in a small red-shift (sketched in figure 7(a)). Mathematically, this can be formulated by the change of the screened interaction W, that enters in both $\Sigma = iGW$ and the electron–hole attraction. Thus the large changes of the electronic properties are largely cancelled for the optical properties. We note that a similar result is found when studying MoS2-2H on MoS2-1T (metallic phase) [95]. The exciton binding energy is 0.10 eV compared to the freestanding monolayer of 0.50 eV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A more detailed analysis is shown in figures 7(b) and (c). In addition to the shift of the A peak, we find that B shifts by similar values. The shift of the 2s state is much stronger, which is due to the larger extent in real space (i.e. the screening acts over a larger area) and the much stronger shift of the band gap (i.e. 2s should be found at a similar fraction between 1s and band gap). A further red-shift can be observed when the 2D material is moved closer to the substrate compared to the 3 Å used previously. The closer the substrate, the stronger the dielectric modification and thus the lower the exciton energy.

In summary, the dielectric screening from the substrate greatly reduces the QP gap of the 2D semiconductor while the lowest lying excitons are only slightly red-shifted.

In figure 8(a) the optical absorption of a neutral MoS₂ monolayer is compared with that of a negatively doped system. The black curve is repeated from the previous section and shows the A, B, and corresponding 2s excitons. For the trions, we find similar resonances at these positions, which are slightly broadened (unlike the bound states, convergence is much more challenging for these unbound trions and is reflected in their broadening [78]). In addition, clear peaks are found at slightly lower energies A$^-_\text{1,2,3}$ and B⁻. The red-shift results from the additional binding energy of the third particle. This can be understood similarly to an H⁻ ion, in which the second electron can bind because it polarizes the system. Note, however, that the three particles here are quasi particles, i.e. they are dressed by their polarisation cloud governed by the dielectric screening in the undoped system. Due to the additional repulsion of the third particle, trions are more extended in real space (see discussion below). In our numerical results, both the convergence and the doping concentration depend on the $\vec{k}$ mesh used. Therefore, a converged result should be understood in the limit of vanishing doping, and the intensities between excitons and trions cannot be compared directly.

Zoom In Zoom Out Reset image size

A similar picture arises for the monolayer of WSe₂. In addition to a resonant feature close to the A exciton, we find three different optically bright trions. We also observe two dark trions below the energy of D. To disentangle their properties, the contributions of holes and electrons are shown in figure 8(c). Similar to figure 4 for the neutral excitons, we find that the two dark states have electrons and holes in bands at K with different spins. The splitting of D$_1^-$ and D$_2^-$ is related to the electron at −K stemming from the lowest or second lowest conduction band. Three different states are possible for the bright trions (three further degenerate states exists with momentum −K instead of K). A$_1^-$ and A$_i^-$ are transitions of an electron from the VB into CB+1 (similar to the A exciton) with an additional electron in CB at −K. For A$_2^-$ the probability of the additional electron in CB+1 is highest at −K. The same characteristic is present for MoS₂, but with a smaller splitting of the states. The origin of these small splittings is related to the details of the Coulomb interaction between the corresponding quasi particles. In particular, the exchange interaction changes if both electrons are in bands with the same or different spin character. A detailed comparison of our results can be found in [82]. Meanwhile, several experimental works have reported the presence of two or more bright trions and a dark trion [85, 102].

We note that similar to the neutral exciton, the trion is also influenced by a substrate. While we find a binding energy of 43 meV in vacuum, this is reduced to 23 meV on a Au(111) surface.

Signatures of trions have been first reported in the 1980s in bulk Ge [103], Si [104], and CuCl [105]. Due to their small trion binding energies of about 1 meV, the study of ground state trions has already been challenging. Since that time, most studies have focused on the investigation of the energetically lowest red-shifted peak, even though the binding energies are much larger, e.g. 130 meV for carbon nanotubes [106]. Therefore, trions have been expected only as the lowest energy state.

On the other hand, figure 8 already suggest that trions may also be possible for energetically higher excitons. Due to the larger spin-orbit splitting (i.e. the B exciton is split off by more than 400 meV), we focus on the study of the WS₂ monolayer. Figure 9 compares experimental and theoretical findings. The characteristics for the lowest states are very similar to WSe₂ discussed above. Below the band gap, further distinct peaks are observed in the experiment: two small peaks in front of the dominant peak with about 24% absorption at about 2 eV and two further peaks a few tens of meV below the step-like feature of the band gap. Our calculations show a similar structure as discussed in figure 8. We find three bright trions T$_\text{A}^{1s}$ below the exciton X$_\text{A}^{1s}$ and similarly T$_\text{A}^{2s}$ below X$_\text{A}^{2s}$. Note that the same characteristic is also present for positive trions (dashed red line, only one trion peak), which we will not further discuss. To verify the character of the excitations in our calculation, we investigate their wave functions in figure 9(c). The top view shows a circle for the 1s states, while a nodal line indicates the 2s character for the excitons in real space. To compare exciton and trion, we show a slice in the right panel. The trion is slightly broader due to the additional repulsion, but the character is similar for both. Therefore, T$_\text{A}^{2s}$ can be identified as excited state trion. Because of additional relaxation channels, the lifetime/width is larger and the observation of two (or three) distinct levels is not possible. In addition to this clear 2s peak, a tiny step also appears at the exciton peak (figure 9(a)). We attribute this to a continuum of scattering states (bound exciton plus one extra electron), i.e. a 'trion continuum', and can conclude that the entire optical spectrum is modified by doping.

Zoom In Zoom Out Reset image size

In addition to our study on WS₂, Wagner et al have recently been able to measure positive and negative 2s trions with similar characteristics in WSe₂ [107].

In figure 10 we show the absorption spectra of a WSe₂ monolayer when applying uniaxial strain. Panel (a) shows the experimental results measured from top to bottom and reveals that the strain can be applied reversible. While this is expected theoretically, it is not a priori clear that no bonds are broken and the phase remains the same. In addition, the key issue is the transfer of the strain from the substrate to the 2D material. The inert 2D compound (which typically only binds via weak van der Waals interactions) has to be connected sufficiently strong to the substrate and, at the same time, the 2D material has to stay intact.

Zoom In Zoom Out Reset image size

After the successful strain transfer is proven, we can track the resonances and their strain-dependent shifts (gauge factors). For the ground-state A exciton and its excited state A$_\text{2s}~-54$ and −55 meV/% are observed, respectively, which is similar to the B exciton (−50 meV/%). The broad C and D features show shifts of +17 and −22 meV/%, respectively. Our calculations in figure 10(b) show similar trends, leading to −43 meV/%, −44 meV/%, −36 meV/%, +10 meV/%, and −19 meV/%, respectively. The physical explanation for the almost identical values of the 1s and 2s states of the A exciton is that they stem from the same bands. And since the spin-orbit split bands have a similar (but not exactly identical) orbital character, the B exciton also shows a similar gauge factor as well. The main contributions from the C and D resonances stem from bands with different character. Especially for the C exciton we find contributions not only at K but also around $\Gamma$. After averaging all contributions between 2.4 and 2.6 eV by their optical weight, we find that the gauge factor has changed its sign. Note that this mixture is also the reason why the excitation energy for C and D is not perfectly linear in our calculations.

Due to the different extent of the wave function of the 1s and 2s excitons in real space, their contributions in reciprocal space are also different, i.e. more localized at ±K. The almost identical gauge factors indicate that the shifts of the excitons can be explained by shifted bands, while the electron–hole interaction remains the same. To verify this, we employ WS₂ including excitons with different momenta. In figure 11 the dashed lines indicate the band gaps at K or between K and $\Lambda$. For a more detailed discussion of the band structure, we refer to section 3.2. Furthermore, the excitons with momentum Q ≈ 0 (optical bright), $Q\approx 1/3$, and $Q = 2/3$ are shown. In particular, for the transition to $\Lambda$ ($Q\approx 1/3$) a nearly perfect agreement of the changes of electronic and optical gap is found. For direct transitions at K, this trend is similar, but the energy difference is slightly smaller for the excitons. Thus, the exciton binding energy changes by about 25 meV/%.

Zoom In Zoom Out Reset image size

In both examples discussed above, we find that more than 90% of the strain-dependent effects can already be understood by the shifts on the single particle level. If the contributions to the relevant excitons are known, the investigation of the quasi-particle band structure is mostly sufficient. Even the use of MBPT theory is not mandatory, since the qualitative effects are well described within DFT (not shown here).

In the previous section the optical properties of neutral excitons in TMDCs have indicated that strain-dependent effects can predominantly be described by the changes of the band structure. Here this conclusion is challenged with the anisotropic 2D material black phosphorus (bP) and, moreover, is also checked for charged trions. In its mono- or few-layer form bP has been discussed frequently since 2014, e.g. its anisotropy or its layer dependence [111–113]. We have investigated the optical properties of excitons and trions in [46]. The energetically low-lying excitations can be found at $\Gamma$ below the band gap of 2.09 eV. The lowest bright exciton has an energy of 1.54 eV and thus a strong exciton binding energy of 0.55 eV. With a trion binding energy of about 80 meV, the (positive) ground state trion resides at 1.46 eV. For further details, we refer the reader to [46]. Figure 12 shows the shift of the GW band gap compared to the exciton and trion energies. For biaxial strain from −2% to 3% our calculations show a nearly linear trend with a shift of about 130 meV/% also for this anisotropic material. For excitons and trions the gradient is minimally smaller. These differences are shown in figure 12(b) as exciton and trion binding energies. The binding energy shifts by about 20 meV/% for excitons, confirming that the main contribution (again nearly 90%) stems from the shifts of the bands. The same holds for trions in which the absolute deviation is similar. However, the relative deviation is consequently larger, and a complete evaluation including many-body effects is desirable.

Zoom In Zoom Out Reset image size

In summary, we find that strain shifts optical excitations very selectively. If the exciton composition is known, the changes of the corresponding band structure on the DFT level are often a reasonable approximation. For better accuracy (e.g. required for trions), the full BSE/trion calculation is required.

O₂ has a triplet ground state. The $\pi^*$ orbitals are both occupied once ($^3\Sigma_\text{g}$). Thus, the system is spin-polarized. Even though it is quite reactive (and thus not welcome in many experiments), it represents one of the simplest cases to study the physical properties of a hybrid system of a 2D material and an adsorbate. Figure 13(a) shows the HOMO and LUMO (doubly degenerated $\pi$ and $\pi^*$) employing different methods. While LDA drastically underestimates the gap with 1.96 eV, GW and GdW (11.80 and 11.62 eV, respectively) show good agreement compared to the experimental ionization potential and electron affinity. The reason for this large correction is very low screening by the molecule itself. On the other hand, this can also lead to large changes when the adsorbate is brought closer to a substrate. In figure 13(b) the results are shown in GdW, in which we model the surface by its dielectric response. Similar to the discussion in section 3.3, we find that the band gap is inversely proportional to the distance due to the image charge effect.

Zoom In Zoom Out Reset image size

For a reasonable description of the heterosystem of O₂ and MoS₂, the van der Waals interaction is crucial, and we optimize the structure using PBE+D3 [25, 59, 60, 117] in the $2\times 2$ cell. The final structure is shown in figure 13(c) with a Mo-O distance of 4.56 Å and $d_0 = 2.98$ Å. The binding energy is 94 meV and quite unspecific with respect to the exact position, as can be seen in the right panel of (d). The electronic properties of the hybrid system are shown in figure 14. The spin-polarized character of O₂ is transferred to the complete system. In panel (a) the band structure is colored in red, overlaid by the one of the pure MoS₂ in dashed black. Compared to the previous results the bands are folded back because of the $2\times 2$ cell. While the O₂ HOMO in DFT-LDA is about 1 eV below the valence band maximum, the LUMO is a mid-gap state. Applying many-body perturbation theory drastically changes this picture. While the electronic gaps in both subsystems change, this is much stronger for O₂ and the HOMO is 4 eV below the MoS₂ VBM, while the LUMO is about 1 eV above the MoS₂ CBM. Even if the interaction is weak, it is essential to allow changes of the quasi-particle wave functions compared to DFT (equation (14)). Furthermore, the band structure shows some more subtle changes. On the one hand, the spin-polarized ground state of O₂ is also reflected in MoS₂. The comparison of −K and +K shows small changes also in MoS₂-dominated bands depending on the parallel or antiparallel alignment of the spin to O₂ (up to 2 meV). On the other hand, the screening of O₂ slightly reduces the gap of the monolayer by about 20 meV. Even though both effects are small, they strongly depend on the adsorbate density and the exact distance to the surface (for more details see [39]).

Zoom In Zoom Out Reset image size

The corresponding optical properties are shown in figure 14(c). For clearer presentation we use $d = d_0-0.8$ Å to enlarge the interaction. In contrast to the spectrum of the free monolayer (figure 3), which is overall very similar, O₂:MoS₂ shows slightly split peaks. They originate from the different coupling to the adsorbate. The contributions to A₁ and A₂ stem from the ±K points (see panel (d)). Therefore, the different electronic coupling is directly reflected in the optical properties. The strength with respect to the distance for O₂ is shown in the inset. Although the strength is small for this test system, larger molecules with stronger magnetic coupling could enhance the discussed effects.

The last section discussed adsorbates on MoS₂, which is a common situation in experiment. Defects are another usual observation. Here we investigate the sulfur vacancy and the possibility that adsorbates can attach to it. Many further defects have been observed [118, 119] and can potentially allow the binding of several adsorbates. Figure 15 compares (a) pure $3\times 3$ MoS₂ with the properties of (b) an empty sulfur vacancy, (c) NH, or (d) trans-aminoazobenzene (AAB) attached. The band structures show that the vacancy gives rise to an unoccupied mid-gap state and a shallow level that is partially occupied. Note that the employed cell is not large enough to decouple the vacancies completely and the corresponding bands are left with some dispersion. The corresponding optical properties show several peaks from about 1.0 eV to more than 3.0 eV. When larger cells are used, these peaks shift upward by about 0.6 eV (not shown). Compared to pure MoS₂ with distinct A and B excitons, the optical response is washed out. If either NH or trans-AAB binds to the vacancy, the defect levels are strongly shifted (i.e. they change their character due to the chemical bond) and no longer appear close to the band edges.

Zoom In Zoom Out Reset image size

Consequently, the low-energy optical properties of MoS₂ can be nearly restored by the present unpolarized molecules. On the other hand, spin-polarization or further selected properties might be introduced by choosing other adsorbates.

In section 3, we have discussed the properties of a monolayer MoS₂ in detail. To investigate the changes in stacked materials, figure 16 shows the distance-dependent changes of the optical properties. The excitation peaks discussed for $d = \infty$ (monolayer) shift to lower energies if the distance of the layers is reduced. Furthermore, an additional peak becomes visible at $d_0+0.6$ Å. In panel (b), we follow the positions and also include (dark) charge transfer states that are absent in the monolayer and shift strongly with the distance. When the B and charge transfer (CT) excitons are close in energy, two bright peaks can be observed. To gain further insight, we investigate the exciton wave function in figures 16(c)–(f) where the hole is fixed in layer 1. While the electron in the CT state is completely located on layer 2, for the intralayer exciton the electron is also at layer 1. Both the second and third excitons have a mixed character. A detailed analysis of the contributions shows that the transitions making up the A exciton are between the VB and CB at the K-point (similar to the monolayer). On the other hand, both mixed interlayer excitons have contributions from VB−1 to CB+1 and from the VB of the other layer to the same CB+1. Note that this mixing is only possible because the spin of both bands is the same.

Zoom In Zoom Out Reset image size

To describe this mixing between the intralayer (B) and CT (A) excitons, we obtain a $4\times 4$ matrix:

$$\begin{equation} H = \left(\begin{array}{cccc} E_{11} & 0 & 0 & \vartheta_h(d) \\ 0 & E_{22} & \vartheta_h(d) & 0 \\ 0 & \vartheta_h(d) & E_{12}(d) & 0 \\ \vartheta_h(d) & 0 & 0 & E_{21}(d) \end{array}\right) \begin{array}{cc}\small \begin{array}{c} (\leftrightarrow \Psi_{11}) \\ (\leftrightarrow \Psi_{22}) \end{array} \\ \begin{array}{c}\small (\leftrightarrow \Psi_{12}) \\ \small(\leftrightarrow \Psi_{21}) \end{array} \end{array}~. \end{equation} \tag{ 43 }$$

On the diagonal, the exciton energies $E_{11} = E_{22} = E^\text{B}$ (intralayer B excitons) and $E_{12} = E_{21} = E^\text{CT} + \frac{C}{d}$ (CT states) enter, where d is the distance and C a constant. These states only couple due to the non-vanishing interaction of the valence bands, denoted as tunneling ϑ_h . The corresponding four parameters are extracted from our calculations and the results of the diagonalization of H are given in figure 16(h). Each solution is two-fold degenerate and nicely reproduces the data points from our BSE calculations. At about 2 Å, both curves get energetically closer and an avoided crossing is visible. Due to the latter the intralayer and CT states mix and both peaks are optically active and share the large optical amplitude from the original B transition. It should be noted that similar models have been previously used to describe bilayer TMDCs [121, 122] as well as few-layer systems [123, 124].

While the nature of the excitations has been disentangled by the first principles calculations, distance-dependent measurements to verify these predictions are complicated. For instance, pressure would change d to closer distances while also modifying the in-plane properties at the same time. An alternative way to distinguish between intra- and interlayer contributions is to apply a perpendicular electric field. While the physics in one layer remains almost unchanged, the different potential between the layers modifies interlayer properties linearly with the field strength (Stark effect). Figure 17 shows the absorption spectra for different E_z up to 0.05 eV Å⁻¹. At zero field all peaks are two-fold degenerate due to symmetry. Either intralayer excitons have the same energy in layer 1 and 2, or interlayer states are from layer 1 to 2 and 2 to 1, respectively. For non-zero field strength the A exciton remains almost unchanged. On the other hand, each mixed interlayer exciton splits in two peaks which shift in different directions. The physical reason for these shifts can be understood from the corresponding bands. Both VB and CB of layer 2 are rigidly shifted. Hence, the degeneracy of the transitions $\text{VB}1\to \text{CB}2$, $\text{VB}2\to \text{CB}1$ is broken, and they change by the potential difference $\pm \Delta V$, respectively. The inset in figure 17 shows the difference between both peaks. A shift of about 46 meV for 0.05 eV Å⁻¹ is found for the mixed interlayer peaks.

Zoom In Zoom Out Reset image size

In a joint work with experimental colleagues from Exeter [125], we have been able to demonstrate the expected characteristic field dependence. A schematic sketch of the device is shown in figure 18(a). The bilayer MoS₂ is encapsulated by hBN, which is placed between two graphene electrodes. The latter are electrically connected by gold. Panel (b) shows the actual device where each layer is marked. The active part of the device has a large area of about $15\times 25\,\mu$m [125]. The results of the field-dependent reflectance contrast measurements are summarized in figures 18(c) and (d). For all E_z a clear peak at the same energy is found slightly above 1.9 eV which is the A exciton. Likewise, a broad B-like excitation is observed. In between an X-shaped peak structure appears. From the zoom in the difference of the gradients of about 0.4 nm can be read, which is the effective interlayer distance perpendicular to the electric field/layers. In theory, we can estimate a similar value. Considering that 35% (65%) are localized on the same (other) layer and the layer distance is 0.6 nm, we find $0.65\times 0.6\,\text{nm} \approx 0.4\,\text{nm}$.

Zoom In Zoom Out Reset image size

For a better comparison with our theoretical results, we model the dielectric environment of the device in figure 19. Compared to the previous results in figure 17, we find two main differences: i) The shift of the first peak (IL$_{1/2}$) is stronger compared to the second one (B$_{1/2}$). The reason for the different mixture is the energetic alignment of the initial CT and intralayer B exciton. Due to the additional screening the CT state is lower in energy and the fraction on inter- and intralayer character is no longer symmetric. ii) Furthermore, we find the 2s peak below B$_{1/2}$. As discussed in section 3.3, both the band gap and the exciton binding energy are strongly modified. In contrast to the 1s peak, which is only slightly red-shifted, the distinct energy shift of A$_\text{2s}$ is another manifestation of the previously discussed effects. Figure 19(b) compares the calculations with the second derivative of the reflectance contrast. Due to the unknown exact device parameters, we have rigidly shifted the energies and scaled the electric fields. Overall a good agreement is found for A and A$_\text{2s}$, as well as the IL$_{1/2}$ and B$_{1/2}$. In experiment, the B peak is relatively broad and prohibits an exact comparison.

Zoom In Zoom Out Reset image size

While it may not be surprising that two subsystems have charge transfer states, the appearance in bulk materials has been unexpected. With the same model for the hybridisation as discussed in equation (43), we can expect that mixed IL states will form in multilayer and bulk if intralayer and CT states are close in energy. In particular, the model Hamiltonian (43) is also valid for the bulk (2H phase) with two layers in the unit cell. In several works [81, 126–128] we have investigated the characteristics in bulk MoS₂, MoSe₂, MoTe₂, WS₂, and WSe₂. While for Mo-based TMDCs interlayer states have been predicted and observed experimentally, the energy alignment (too close to A, small oscillator strength) in W-based materials hinders their detection.

In most homobilayers the energetically lowest state is an intralayer state. The reason for this is that the Coulomb interaction is distance dependent and thus smaller between two layers compared to the layer itself. In the last section we discussed the change of the initial band alignment by an electric field. If two different layers are combined, the bands are typically not aligned due to the different band gaps and work functions. In such a heterobilayer the exciton ground state may be an interlayer transition with a particularly long lifetime.

Prototypical examples are again TMDCs. If we stack two of these layers, the system is no longer commensurate due to the different lattice constants, or can even have rotated lattices. Although these structural differences can lead to many interesting phenomena, we will focus on the properties of commensurable systems which are computationally manageable with our methods. In figure 20, we examine the MoS₂/WS₂ bilayer. The lattice constant of WS₂ is $a_\text{lat}^{\text{WS}_2} = 3.155$ Å [130] which is only 0.2% smaller compared to MoS$_2~a_\text{lat}^{\text{MoS}_2} = 3.160$ Å [131]. Here we employ a = 3.160 Å and neglect the small effects due to strain (see section 5). MoS₂/WS₂ forms a type II bands alignment, i.e. the VBM stems from WS₂, while the CBM is located in MoS₂. Note that the strongly hybridized bands around $\Gamma$ and halfway between K and Γ may seem to prohibit such a definition. However, for optics we are mainly interested in physics close to K, where the attribution to individual layers is still possible. Thus, the optical properties of the heterobilayer share several similarities with the monolayers. The first strong absorption peaks can be found as A- and B-like excitons. Figure 20(d) shows that both peaks have predominantly intralayer character, i.e. the probability to find the electron and hole in the respective layer is largest. In addition, a tiny peak can be observed about 200 meV below A. From its exciton wave function with hole and electron on different layers, we can clearly identify this as an interlayer state. Although it might be challenging to measure such small oscillator strength in absorption, clear evidence has been found in photoluminescence due to the high occupation [132, 133]. We find an oscillator strength of about 1% compared to the A$_{\text{WS}_2}$ in our calculations. The intensity strongly depends on the distance between the layers in figure 20(e). In a Moiré structure this may suppress the oscillator strength in regions where the atoms are not aligned [134].

Zoom In Zoom Out Reset image size

In addition to the neutral excitons, trions could also be present in the case of a heterobilayer. Figures 21(a) and (b) show a zoom in at low energies for positively and negatively charged trions. The lowest peaks can be found at a trion binding energy of $E_\text{b}^{\text{T}^+_\text{IL}} = 18$ meV and $E_\text{b}^{\text{T}^-_\text{IL,1}} = 28$ meV, respectively. Interestingly, we find that the nature of this quantum mechanical state can be counterintuitive. Both particles with the same charge can be found on the same layer, i.e. both holes on WS₂ or both electrons on MoS₂, respectively. The trion with one hole on both layers $\dot{\text{T}}_\text{IL,(1/2)}$ can be found close to the interlayer exciton and appears as its scattering state (which is not fully converged). This assignment is further demonstrated from our ab initio results in figures 21(c) and (d). Although not shown in detail here, the negative T$^-_\text{IL}$ further splits in three states (like for the intralayer trions), while T$^+_\text{IL}$ remains a single peak. Note that the signature of a trion can generally appear several times in the optical spectra (figure 21(d)). This is a result of detection at $\Omega_{T\to \nu} = E^T - \epsilon_\nu$, where $\nu$ can be different final states. For example, $\dot{\text{T}}_\text{IL,(1/2)}$ is responsible for the peak discussed above, and another one is present slightly above the MoS₂ exciton (not shown).

Zoom In Zoom Out Reset image size

Calman et al have experimentally measured interlayer trions in MoSe₂/WSe₂ with a binding energy of 28 meV [135] and thus have validated our prediction.

As for conventional 3D materials, the number of hypothetical 2D materials is essentially unlimited. Thus, it is important to select reasonable candidates and start with a good guess of their lattice and atomic structure. The most important classification scheme is the symmetry. Figure 24(a) shows the five different 2D Bravais lattices. Together with the basis, this defines different materials. Some examples are shown in figure 24(b), e.g. graphene [7], MoS₂ [10], WTe₂ [154], CrI₃ [155], TiS₃ [156], FeSe [157], or FeOCl. For our screening study, we use already known 2D materials as prototypes and substitute their atoms by chemically similar ones. For example, we replace carbon with silicon, germanium, tin, etc in graphene, or molybdenum with other transition metals and sulfur with other chalcogenides in MoS₂. Due to this structure generation strategy, most of the candidates are quite simple, i.e. they are binary or ternary compounds, and the number of atoms in the unit cell is limited to a small number.

Zoom In Zoom Out Reset image size

With these structural guesses, we begin our structure workflow, which is shown in figure 25. After an initial relaxation which preserves the symmetry, we check that the system is still one two-dimensional layer and that it is not already included in the C2DB. If it is, we continue with the classification including its magnetic state. To check the stability, three different criteria are evaluated: 1) The formation energy needs to be less than zero which means that the total energy of the atoms reduces (relative to the standard states of the constituent elements) due to the interaction with other atoms. If this is not case, a small perturbation or even a temperature larger than 0 K would be sufficient to degenerate the material. We also compare the total energy to those of other phases, e.g. for graphene we compare to diamond etc. If the energy of the current 2D configuration is much higher, it is most probably unstable. 2) We evaluate the stiffness tensor and its eigenvalues. If any of them is negative, this is a clear sign that the structural relaxation has not found a minimum on the energy surface. 3) The dynamic matrix is set up at the high-symmetry points of the $2\times 2$ supercell, and we calculate its eigenvalues. If an eigenvalue is below zero, this indicates that the system would further relax in a larger cell and is therefore unstable in the current configuration. If the material is unstable, we store this result in the database and do not consider the material further. If, on the other hand, the material is stable, we continue to investigate its properties in more detail. Materials of particular interest for opto-electronics are transition metal dichalcogenides for which we find about 70 candidates with different compositions that are stable (see MoS₂ and WTe₂ in figure 24).

Zoom In Zoom Out Reset image size

Once the structure has been determined, we are able to evaluate a large set of properties automatically by our ab initio property workflow, which includes Bader charges, Born charges, out-of-plane dipoles, work functions, exchange couplings, magnetic anisotropies, deformation potentials, plasma frequencies, piezoelectric tensors, spontaneous polarizations, band structures (within different approximations) and density of states, effective masses, Fermi surfaces, topological invariants, infrared polarizabilities, optical polarizabilities, optical absorbances, Raman spectra, and second harmonic generation. For a more detailed discussion of all these properties, we refer to [3, 150] in the appendix. Here we briefly discuss the statistics of band gaps and band positions. Figure 26 shows the resulting band gaps for more than 250 semiconducting materials for which we used different approximations: PBE [25] is a widely used generalised gradient approximation, HSE [27] replaces 25% of the local exchange by the non-local exchange potential for the short-range part of the Coulomb interaction, and GLLBSC [158] is an approximation to the optimised effective exchange potential that includes an estimate of the derivative discontinuity for the band gap. The band gaps obtained with these approximations are plotted in blue, orange, and green with respect to the $G_0W_0$ results. The left panel shows the band gaps and indicates that PBE and HSE clearly underestimate the $G_0W_0$ values on average by about 45% and 22% in 2D materials. GLLBSC provides a much better estimate for the band gap, but the position with respect to vacuum (right panel) and distances to other bands (not shown here) are problematic. On the other hand, the predicted band positions in PBE and HSE clearly follow the trends of the $G_0W_0$ results. While we perform $G_00W_0$ and BSE calculations for semiconductors with less than six atoms per unit cell, all other properties are calculated using the PBE functional.

Zoom In Zoom Out Reset image size

While DFT has become a standard method and is relatively easy to use, MBPT is much more involved. This is due to several parameters which have to be converged simultaneously, such as the additional basis representing the two-point functions, or the $\vec{q}$ integral (over the first Brillouin zone) and the energy dependence of the self-energy. Moreover, MBPT in the GW approximation scales an order of magnitude slower, i.e. $\mathcal{O}(N^4)$ where N is the number of electrons (compared to $\mathcal{O}(N^3)$ in DFT), and the prefactor is much larger. Depending on the exact implementation, this only allows for small systems of a few dozen of atoms or less. Therefore, high-throughput calculations are particularly challenging [38].

As a first step towards fully automated GW calculations, we have analysed a set of 370 GW calculations for 2D materials to establish the validity and limitations of various commonly used approximations and extrapolation schemes [38]. As an example, we briefly mention one aspect pertaining to the renormalization factor Z. The left panel of figure 27 shows the GW band structure of monolayer MoS₂. While most points are reasonable, the third conduction band at K shows a drastic and unphysical jump. This can be traced back to the moderately converged derivative of the self-energy (equation (13)), which enters by the quasi-particle weights $Z = \lbrace 1-\frac{\partial \Sigma}{\partial\omega}\rbrace^{-1}$. This behavior is observed for other materials as well. In general such cases can be identified by the condition $Z \not\in [0.5, 1.0]$, which indicates either a numerical/convergence issue or a break down of the quasi-particle approximation. We refer to such states as quasi-particle inconsistent states. While for a single material more rigorous convergence can be helpful, this is impractical in a database. Another possibility is to set Z to the average value 0.75. In the right panel of figure 27 the resulting changes of the quasi-particle shift are shown. While the distribution is slightly broader, the absence of outliers reduces the mean absolute error to 0.09 eV. An even better result can be obtained if Z = 0.75 is used only when $Z \not\in [0.5, 1.0]$ (MAE becomes 0.06 eV, not shown here).

Zoom In Zoom Out Reset image size

A novel class of 2D materials that we have investigated are the so-called Janus monolayers. These are materials, for instance TMDCs, in which the atoms at the top and bottom of the monolayer are not identical such that the out-of-plane mirror symmetry is broken. Lu et al [159] have fabricated and measured MoSSe in 2017. Its structure (see figure 28(a)) is similar to that of the TMDCs discussed previously. The lattice constant of 3.241 Å is in between those of MoS₂ and MoSe₂ (3.16 and 3.30 Å). Due to the different work functions above and below the layer, the potential differs by 0.75 eV, so the material has a build-in dipole.

Zoom In Zoom Out Reset image size

The band structure in figure 28(b) shows the typical direct gap at the K point. The gap of 2.64 eV is (again) in between those of 2.89 eV for MoS₂ and 2.57 eV for MoSe₂. This is also true for the spin-orbit splittings of the VBs and the CBs (210 meV and 30 meV), compared to MoS₂ (180 meV and 15 meV) and MoSe₂ (230 meV and 42 meV). Similarly, the lowest CB has opposite spin polarization compared to the highest VB. But unlike MoS₂ and MoSe₂, we do not find that the CBs cross. Instead, the character changes due to an avoided crossing.

The optical absorption shows a similar trend (figure 28(c)). Two clear excitons A and B originate from direct transition between the spin-orbit-split bands. In addition, excited states (2s, etc) and a broad peak structure C (originating from transitions in different parts in the Brillouin zone) reside below the band gap. The first bright exciton can be found at 1.88 eV, and the first dark exciton is about 25 meV lower in energy. Compared to previous results (25 meV in MoS₂) [82], we expect that the lowest excitation is dark in MoSSe similar to MoS₂. MoSSe is only one possible candidate. Employing tellurium or changing the transition metal (to tungsten or similar) suggests many further interesting Janus materials. E.g. we found semiconducting materials with direct and indirect band gaps ranging from 0.7 eV for CrSTe to 3.0 eV for BiClS. These have finite out-of-plane dipoles with a difference of the vacuum level up to 1.5 eV as well as large Rashba-splittings up to 75 meV.