Optical excitations in 2D semiconductors

Two-dimensional (2D) materials have revealed many fascinating physical and chemical properties. Due to the quantum confinement and enhanced many-body effects especially the optical properties are altered compared to their bulk counterparts. The optics of 2D materials can easily be modified by various means, e.g. the substrate, doping, strain, stacking, electric or magnetic fields. In this review we focus on the theoretical description of the excited states and optical properties of 2D semiconductors paying particular attention to the current challenges and future opportunities. While the presented methodology is completely general and applicable to any 2D material, we discuss results for the transition metal dichalcogenides, their heterostructures, and some novel materials from the computational 2D materials database.


Introduction
Two-dimensional (2D) materials are crystalline, covalently bonded, and chemically inert sheets of atoms that measure just a few angstroms in thickness, but can have lateral extensions in the micro-or millimeter range. Such materials can be exfoliated from layered bulk materials or synthesised bottom-up using e.g. chemical vapour deposition or pulsed laser deposition. In 2D materials, the confinement of electrons to the atomically thin layers result in highly anisotropic properties. For example, in-plane bonds are typically strong and covalent or ionic, while weak van der Waals bonds act in the out-of-plane direction. This anisotropy is also reflected in the electronic and optical properties and eventually leads to many interesting phenomena such as valley polarization, charged excitations, quantum confinement, hyperbolic response, etc [1]. There are many two-dimensional materials; more than 50 have been synthesised or exfoliated in monolayer form while many hundreds are predicted to be exfoliable from already known layered bulk compounds [2,3]. Several materials have been used for a long time in their bulk forms in which hundreds of layers are stacked on top of each other. For instance, minerals like graphite or MoS 2 have been known since thousands of years [4] and their atomic structures have been clarified several decades ago [5,6]. Although techniques for cleaving and thinning are well established, it took until 2004 for graphite to be first observed in its monolayer form graphene [7][8][9]. In the following years, many further materials such as the monolayers of several transition metal dichalcogenides (e.g. MoS 2 [10]) or black phosphorus [11] have been observed experimentally.
In particular, conventional semiconductors like silicon or gallium arsenide have been studied intensively long before two-dimensional materials, e.g. due to their relevance for modern computer technology, integrated electronics, or solar cells. Silicon has been used in most information technology devices in daily life for more than 50 years. Because of its indirect band gap, efficient coupling of electronic and optical properties is challenging and silicon has a vanishing gap in its two-dimensional phase [3,12]. However, this long-standing attempt promises easier interfacing with today's fiber optic communications and faster and more parallel information processing. Due to the reduced dimensionality, light-matter interactions in many two-dimensional materials are stronger than in their three-dimensional counterparts [1,13]. This also makes them promising candidates for renewable energy and sustainability applications (e.g. solar cells or light-emitting diodes).
A fundamental understanding of light-matter interactions at the atomic length scale is mandatory for targeted development of novel and improved opto-electronic devices. In this work, we will theoretically investigate the optical properties of various two-dimensional materials and different ways to modify them. All investigations will be performed using ab initio techniques which allow for a parameter-free prediction of the corresponding physical properties and even the prediction of novel hypothetical materials. These highly accurate and system-independent methods are accompanied by a high numerical cost, i.e. highly efficient implementations in computer programs for supercomputers are required. They develop their full potential gaining new physical insights in combination with experimental observations. The detailed theoretical analysis allows clarification of the physical mechanisms involved, and at the same time the experiments verify the theoretical predictions.
We start this work with a brief overview of the state-of-the-art ab initio techniques in section 2. The starting point is the well established density functional theory (DFT), explicitly corrected by van der Waals interactions, which are not well described by the standard functionals. The variational principle can be used to determine the lowest total energy and the corresponding density for the ground state of a system. The minimal energy with respect to the atomic positions determines the optimal structure. The electronic and optical properties require a more accurate description including excited states and are today mostly based on many-body perturbation theory. Single particle processes such as band structures can be approximated very reliably by single-particle Green functions and their equation of motion in the GW approximation. Optical properties are described by the transition between two different electronic states preserving the number of particles. Consequently, their description requires the two-particle Green function and its equation of motion (the Bethe-Salpeter equation). Due to the electron-hole interaction, bound excitons can form and significantly change the spectrum. Moreover, three-particle excitations (trions) are important for many naturally doped 2D materials and lead to additional red-shifted peaks. Their theoretical description requires a rigorous quantum mechanical treatment, but is often missed in simulations. We conclude the section with novel theoretical approaches to evaluate the influence of magnetic fields on electronic and optical properties.
In the following sections we investigate the optical properties and different ways to modify them in detail. Figure 1 summarizes several external stimuli which can be used in two ways: On the one hand, these modifications are suitable to design the optical properties, and on the other hand, they provide insight into the physical nature of the underlying excitations. In section 3, we begin the discussion of monolayers of transition metal dichalcogenides and examine their characteristic optical spectra. Due to the electronic band structure, the low-energy excitations mostly originate from transitions close to the ±K valleys. In addition to the large exciton binding energies, which include a Rydberg series of excitons, the strong spin orbit coupling leads to two excitons with similar amplitudes. We take a closer look at possible excitations which do not couple directly to optical measurements. Such excitations may be dark due to the character of the wave function, e.g. if the spin of the electron and hole are reversed. Especially in transition metal dichalcogenides with their large spin orbit coupling of the valence bands, but energetically very close conduction bands, there is a close competition which state is the energetically lowest. If the difference of the momenta of both particles, electron and hole, is larger than the momentum of the photon, no optical transition is possible in the first place. However, the momentum can be transferred to a phonon and an assisted optical transition is possible, for example in photoluminescence experiments. In this work we focus on the calculation of the momentum-dependent exciton dispersions. We also investigate the modification due to a possible substrate. While such a substrate is artificial in our calculations, it is always present in experiments, and thus modifying the substrate is the most fundamental way to change the optical response. Section 4 introduces another common perturbation in 2D materials, the doping. In many experiments, samples are unintentionally doped with additional electrons. On the other hand, the doping can also be tuned intentionally if a voltage is applied. In our theoretical calculations we assume small doping concentrations, which allows us to neglect the influence on the band structure and limit ourselves to three-particle correlations. A bound state (trion) can appear below the corresponding exciton. This is not only true for the lowest energy transition, but we have also been able to discover trions for the second Rydberg state.
Due to mismatch of the lattice constants or imperfections in the sample preparations 2D materials may naturally be strained locally or on a larger scale. Strain can also be selectively tuned by bending the substrate. In section 5 we investigate the influence of the electronic and optical properties from first principles by modifying the structural properties. Biaxial strain can be simulated by changing the lattice constant, while uniaxial strain also distorts the structure and breaks the hexagonal symmetry. The resulting structural changes shift the bands according to the composition and symmetry of the corresponding wave functions. At small strain values an almost linear response can be observed, which is characteristic for each band. As a consequence, the excitation energies of the excitons exhibit characteristic shifts that can be used to identify them. On the other hand, a sufficiently large strain may even be utilized to change the energetic order of the excitons.
Two-dimensional materials are almost only surface. Like any other material, they have a certain concentration of defects. Additional atoms and molecules can adsorb both at the defect-free areas and at defects sites. Section 6 examines the properties of some exemplary molecules. On the one hand, we have chosen an oxygen molecule which has a triplet ground state, making it one of the simplest examples of a magnetic molecule. Even if the material is inert, a weak van der Waals bond can form, and a certain overlap is observed. This changes the electronic and optical properties of the monolayer and depends on the spin character of the molecule. On the other hand, vacancies and the adsorption of different molecules therein are studied. Several molecules are able to recover the defect-free properties of the 2D material very well. Section 7 introduces the stacking of individual layers to form homo-or heterostructures. Several 2D materials exist in such naturally layered three-dimensional structures, e.g. graphene, hexagonal boron nitride or MoS 2 . In addition, it is possible to build artificial structures like twisted homobilayers or heterostructures. Even if the layers are connected only by van der Waals interaction, several interesting phenomena can be observed. For example, a monolayer of transition metal dichalcogenides shows a direct band gap, while multiple layers have indirect band gaps. Besides the resulting changes of the optical spectrum, spatially indirect states are present on different layers which have no counterpart in a monolayer. Such interlayer excitations consist of a hole and an electron on two different layers and are the analog to charge-transfer states in molecules. In particular, in systems where interlayer excitons are the ground state, their long lifetime might be promising. The level alignment of the different layers can also be utilized to distinguish interlayer states from other excitations by employing a perpendicular electric field. A further interesting question concerns the changes of the excitations due to additional doping, i.e. interlayer trions.
In section 8 we calculate the response of the electronic and optical properties due to an external magnetic field. While the use of electric fields orthogonal to the layer is conceptual simple in the supercell method, magnetic fields break the translational symmetry in the layer and seem to prohibit the applicability of our periodic calculations. However, by using Berry phase techniques the problem can be reformulated and the band-and ⃗ k-dependent shifts can be evaluated. By employing the Bethe-Salpeter equation, we can connect such shifts to the shifted exciton energies and eventually to the changed optical spectrum. Even if magnetic fields probably have the smallest absolute impact, they can still be used to split excitons with different character, e.g. the ±K valleys in transition metal dichalcogenides.
In section 9 we consider the central object, the two-dimensional material itself, from a broader perspective. After presenting a high-throughput screening approach for the discovery of novel hypothetical materials, we address the statistics of some material properties like band gaps or exciton binding energies. The large set of systems with very diverse physical and chemical properties is also an ideal test set to compare different methods, such as those discussed in the first section.

Theoretical methods
A quantum mechanical system can be described by the Schrödinger equation or, more generally, by the Dirac equation. In the time-independent case this is given byĤΨ = EΨ, where E are the energies and Ψ are the wave functions. While the HamiltonianĤ is known exactly even for systems with macroscopic extent, the direct solution for O (10 23 ) particles is out of scope (even after the ionic motion is decoupled via the Born-Oppenheimer approximation). Different concepts have been developed to rewrite the problem, followed by different approximations, of which the two most commonly used are presented here: The DFT, or more precisely the Kohn-Sham scheme of DFT, models the interacting electron system as a non-interacting system moving in an effective potential, and is well suited for calculating the ground state and related properties. Many-body perturbation theory (MBPT) provides a systematic way of accounting for the two-body Coulomb interaction beyond the mean-field level, by summing certain subsets in the perturbation series to infinite order. In this way, it is possible to obtain an accurate description of the optical excitations. Thus, the combination of both techniques is a suitable set of tools to describe and predict the physical properties in this work.
After a brief summary of the two techniques in sections 2.1 and 2.2 and the practical aspects in section 2.3 (for a more detailed discussion, we refer the reader e.g. to [14][15][16][17]), the new developments for doping (trions) are discussed in section 2.4 and for magnetic fields (Zeeman shifts) in section 2.5.

DFT
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 ρ(⃗ r) = occ i |φ i (⃗ r)| 2 by a unique functional. By varying them, the ground state density ρ 0 (⃗ r) can be found for which the functional derivative becomes zero. Eventually, the effective one-particle equation for the ith particle (wave function φ i (⃗ r) and energy ϵ i ) is 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 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 ρ(⃗ r) and ∇ρ(⃗ 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][30][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.

Many-body perturbation theory
In many-body perturbation theory (MBPT), the discussion starts from the many-body ground state |N⟩ [33]. Using the shorthand notation for position, spin, and time (⃗ xt) = (⃗ r, σ, t) = 1 and atomic units (i.e. h = 1,h 2 /2m = 1, and e 2 = 2 leads to lengths in a B and energies in Ry), the causal one-particle Green function is given by It describes the creation of an electron at 1, its propagation and annihilation at 2, or vice versa.T is the time-order operator which adds a negative sign when t 1 > t 2 andψ ( †) (1) are field operators for the annihilation (creation) of an electron that can be expressed asψ â † 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 where V Coul is the Hartree potential, links the one-particle Green function to the two-particle Green function G 2 , and we eventually find an infinite hierarchy of equations. The hierarchy can be broken by introducing the self-energy Σ In principle, the self-energy can be obtained self-consistently by five equations, the so-called 'Hedin's equations' [33]: Equation (6) is the Dyson equation that links the Green function to the non-interacting G 0 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) = δ(t1−t2) |⃗ r1−⃗ r2| 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 Σ = 0 in equations (6) and (10) and the self-energy simplifies to W can be calculated via the dielectric function in the random phase approximation (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 ψ QP m ⃗ k (⃗ x) and energy E QP m ⃗ k of the quasi-particle (m includes its spin) can be calculated by Since Σ depends on ψ 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 and } denotes the renormalization constant including the derivative of the self-energy with respect to the energy. Z m ⃗ k describes the spectral weight of the corresponding quasi-particle and its value is typically around 0.7 [38]. Transformed to ⃗ k space, we can finally calculatẽ and 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 QP In this case the linearized form is given as E QP 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 LDA xc ≈ i GW metal for many systems, an alternative approximation is possible [43,44] This GdW(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 ≡ 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 2 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 ± 1) particle system, G 2 contains electron-hole pairs (we neglect possible hole-hole or electron-electron pairs). The ground state |N, 0⟩ can therefore be excited by the creation and annihilation operators of G 2 to a state S, denoted |N, S⟩, without changing the number of particles. The two-particle Green function is defined analogously to equation (3) as For excitons, we can assume that an electron-hole pair is created at times t 1 and t ′ 1 = t 1 + 0 + and annihilated at t 2 and t ′ 2 = t 2 + 0 + . We note that G 2 and the corresponding correlation function L are four-point functions Its equation of motion, the BSE is given by where L 0 (12; 1 ′ 2 ′ ) = G 1 (1, 2 ′ )G 1 (2, 1 ′ ) 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 (E c − E v ) are the energy differences of valence and conduction bands and K are the corresponding matrix elements of the interaction kernel (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; (3,6)δ(4, 5)W(3 + , 4). By solving equation (19) we can find the exciton energies Ω S and the amplitudes A S and B S of the excitations and de-excitations ). 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 where the dimension of the eigenproblem can still become large, given by the product of the hole and electron bands and the number of ⃗ 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 ⃗ q → 0 For ω = 0 the dielectric constant ε 0 is obtained. Because the absorption spectrum is given by the imaginary part (ε M = ε 1 + iε 2 ), we will focus on ε 2 (ω) in this work. To evaluate the coupling to the light field, we calculate the imaginary part as the sum of all possible excitonic states whose weights are determined by the normalized polarization of the light ⃗ A/| ⃗ A| (we assume a linear polarization) and the velocity matrix elements. These elements can be calculated as ⟨0|⃗ v|S⟩ = occ v emp c A S vc ⟨v|⃗ v|c⟩, therefore both required values (the excitation energies Ω S and the electron-hole amplitudes A S vc ) are known after diagonalizing the BSE Hamiltonian.

Some practical aspects
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][52][53][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 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 ⃗ τ ν For details and an efficient numerical implementation, we refer to [56,57]. The coefficients c α,ν n, ⃗ k are obtained by the solution of the eigenvalue problem 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.
use the Gaussian orbitals φ n, ⃗ k from DFT and again Gaussian orbitals or plane waves for χ β (⃗ q,⃗ 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].

Three-particle states
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 3 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 aŝ 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, ⃗ k i ). Within Hartree-Fock theory, the single-particle states ϕ n (x) are the solutions of the Hartree-Fock equations, and the matrix elements can be calculated by where x = (⃗ r, σ) and the bare Coulomb interaction is V(x, x ′ ) = e 2 /|⃗ r −⃗ r ′ |. For an electron-hole pair configuration |vc⟩ :=â † câv |0⟩ the matrix element is with ϵ (HF) n being the band-structure energy and E (HF) 0 the ground-state energy. Note that equation (30) is not diagonal and the diagonalization of the Hamiltonian (31) leads to the excitation energies Ω S and the excited states, which are linear combinations of the free electron-hole pair configurations |S⟩ = vc A S vc |vc⟩. 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 For the three-particle states, we take a similar route. For the configurations of two electrons and one hole , we can derive the matrix elements within Hartree-Fock [65,66] 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 1 and c 2 , both terms are screened. Finally, we find the matrix elements and by diagonalizing the corresponding Hamiltonian we can calculate the energy of the trions E (T, ⃗ K) and their wave functions for trions T with momentum ⃗ K. Because of momentum conservation, the wave numbers of vc 1 c 2 have to Note that the trion in an optical spectrum is not found at the trion energy E (T, ⃗ K) , but at the transition energy as an electron remains in the conduction band |c ⃗ K⟩ (or is already present before absorption). Note that the ⃗ 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].

Magnetic fields
The general treatment of electro-magnetic fields in a periodic calculation is challenging. Corresponding to the electric and magnetic fields ⃗ E and ⃗ B, the (vector) potentials φ and ⃗ A can be determined by These potentials φ and ⃗ 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][71][72][73][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Ĥ eff 0 |Ψ 0 In a next step, we include the magnetic field (along the out-of-plane z-axis) as a perturbationĤ HereL andŜ are the angular momentum and spin operators, respectively, where we assume thatL 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 The destroyed lattice periodicity is reflected in the spatial dependency of the angular momentum operator L z = xp y − yp x which prohibits a direct evaluation. Therefore, we follow the derivation of Chang and Niu [70] for the magnetic moment of a wave packet where u n ⃗ 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 The g factor of the exciton S is often related to the difference of corresponding experimental measurements with σ +/− circularly polarized light gµ B B := Ω σ + − Ω σ − [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 ⃗ k) to (c ⃗ k) could be approximated by 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.

Excitations of free-standing 2D materials and substrate screening
A particularly well-suited class for the investigations of opto-electronics properties are the group-6 transition metal dichalcogenides (TMDCs) MX 2 (M = Mo/W, X = S/Se/Te) for which many experimental results are available [77]. For this reason we will mainly focus on the TMDCs in the following sections. Some electronic and optical properties of other 2D materials will be studied in section 9.
TMDCs consist of three atoms per unit cell: one transition metal in the central plane and two chalcogenides above and below. A side view of the hexagonal structure is shown in figure 1 for MoS 2 . The following discussions are based on calculations employing the experimental lattice constants, which are close to the theoretically optimized ones [78]. After we have discussed the ground-and excited-state properties of free-standing monolayers in vacuum, the influence of the substrate is analysed and is shown to provide a simple, yet efficient means to tune the properties of the 2D material.

Free-standing monolayers of transition metal dichalcogenides
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 2 and WSe 2 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 2 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 ⟨S z ⟩ 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 2 , the VBs are split by 570 meV, while the value for MoS 2 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.
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 2 compared to MoS 2 . Furthermore, peaks with lower optical weight corresponding to excited states (A 2s , B 2s ) are observed. Similar to possible excitations in a hydrogen atom, excitons can also form a Rydberg-like series [83,84].
A detailed understanding of the excitations is possible by analyzing their wave functions. In figure 4(a), the A exciton in MoS 2 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 2 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 2 , the appearance of D as the lowest state is not surprising. Interestingly, also in MoS 2 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 2 , WS 2 , and WSe 2 . For MoSe 2 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].

Finite-momentum exciton landscape
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 ⃗ k + ⃗ 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) 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 Å −1 and Q ≈ 1.3 Å −1 . These correspond to transitions between Λ or the high-symmetry points (K→ Λ 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   (111) substrate. The substrate distance is fixed to 3 Å. A sketch of the corresponding systems, the band gap Eg, the excitation energy Ω A , and the resulting binding energy E b is given below. Reproduced from [78]. CC BY 4.0.
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].

Substrate screening
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 2 placed on the semiconductor SiO 2 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 Ω 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][93][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 Σ = 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.
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.

Three-particle excitations in doped materials
In the last section we have discussed several properties of neutral 2D materials. However, it has been observed that an additional peak (smaller in intensity and slightly red-shifted to the exciton) is typically present in grown and exfoliated samples [96,97]. This peak has been attributed to the natural doping in the system and the formation of excitations consisting of three quasi-particles, i.e. negatively charged trions consisting of two electrons and one hole. In recent years, it has become possible to build gated devices (see, among others, [98,99] and many more) in which the doping concentration can be tuned.
Theoretically, trions have often been modeled using parameter-based approaches [100,101]. Having introduced an ab-initio method in carbon nanotubes [65], we use an extended version that includes spin-orbit interactions (equation (33)) to evaluate trionic properties in 2D materials.

Trions from first principles
In figure 8(a) the optical absorption of a neutral MoS 2 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 − 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 ⃗ 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.
A similar picture arises for the monolayer of WSe 2 . 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 2 , 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.

Excited-state trions
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 2 monolayer. Figure 9 compares experimental and theoretical findings. The characteristics for the lowest states are very similar to WSe 2 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 1s A below the exciton X 1s A and similarly T 2s A below X 2s A . 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 2s A 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.
In addition to our study on WS 2 , Wagner et al have recently been able to measure positive and negative 2s trions with similar characteristics in WSe 2 [107].

Strain-dependent effects
One of the conceptually simplest approaches to probe and to tune the properties of a material is to change its structure by applying strain [108,109]. In section 3, we have studied the band structure of TMDCs and found that it can be easily modified, e.g. by the dielectric environment. Therefore, a similarly strong dependence on the internal structure can be expected in 2D materials in which it is easy to apply strain.

Strain-dependent excitons in transition metal dichalcogenides
In figure 10 we show the absorption spectra of a WSe 2 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.
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 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 Γ. After averaging all contributions between 2.4 and 2.6 eV by their optical weight, we find that  [110]. © IOP Publishing Ltd. All rights reserved. Figure 11. Strain-dependent changes of the quasi-particle gaps at K (red dashed line) and between K-Λ (blue dashed line) in monolayer WS2. The solid lines show the excitons at Q = 0, 1/3, and 2/3 reciprocal units (red, blue, and black dots, respectively). The minimal direct gap and the A exciton are set to zero for Q = 0. Note that the curves partially overlap (the lines are linear fits). Reproduced from [86]. © IOP Publishing Ltd. All rights reserved. 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 2 including excitons with different momenta. In figure 11 the dashed lines indicate the band gaps at K or between K and Λ. 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 ≈ 1/3, and Q = 2/3 are shown. In particular, for the transition to Λ (Q ≈ 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/%.
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 The strain is applied biaxially from −2% to +3%. Here we employ the GW approximation to evaluate the quasi-particle gap and electron-hole matrix elements. (b) Energy shift of the exciton and trion binding energy on an enlarged scale with respect to (a). Reproduced from [46]. © IOP Publishing Ltd. All rights reserved. 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).

Strain-dependent excitons and trions in black phosphorus
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][112][113]. We have investigated the optical properties of excitons and trions in [46]. The energetically low-lying excitations can be found at Γ 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.
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.

Adsorbates on 2D materials
Due to their large surface to volume ratio, 2D materials are very sensitive for adsorbates and allow for several ways of tuning. If the adsorbates are not radicals, they can either attach to the pure materials by van der Waals interactions, or, if present, form a chemical bond to defects. Their interaction will influence the adsorbed molecule and the 2D material electronically and optically. Following our previous works [39,114], we will discuss both aspects for an oxygen molecule adsorbed on MoS 2 and a sulfur vacancy in MoS 2 with NH or azobenzene attached.

Oxygen on MoS 2
O 2 has a triplet ground state. The π * orbitals are both occupied once ( 3 Σ 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 π and π * ) 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 The binding energy of O2 at various positions above MoS2 is shown on the right. We employ PBE+D3 to capture the essential van der Waals interactions, which leads to a maximal binding energy above the Mo atoms. Reprinted figure with permission from [39], Copyright (2020) by the American Physical Society.
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.
For a reasonable description of the heterosystem of O 2 and MoS 2 , the van der Waals interaction is crucial, and we optimize the structure using PBE+D3 [25,59,60,117] in the 2 × 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 2 is transferred to the complete system. In panel (a) the band structure is colored in red, overlaid by the one of the pure MoS 2 in dashed black. Compared to the previous results the bands are folded back because of the 2 × 2 cell. While the O 2 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 2 and the HOMO is 4 eV below the MoS 2 VBM, while the LUMO is about 1 eV above the MoS 2 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 2 is also reflected in MoS 2 . The comparison of −K and +K shows small changes also in MoS 2 -dominated bands depending on the parallel or antiparallel alignment of the spin to O 2 (up to 2 meV). On the other hand, the screening of O 2 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]).
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 2 :MoS 2 shows slightly split peaks. They originate from the different coupling to the adsorbate. The contributions to A 1 and A 2 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 2 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.

Electronic repair of MoS 2
The last section discussed adsorbates on MoS 2 , 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 × 3 MoS 2 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 2 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.
Consequently, the low-energy optical properties of MoS 2 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.

Interlayer excitations in stacked materials
After we have discussed several ways to modify the monolayers themselves, we investigate what happens if two or more layers are joined into homo-or heterostructures. In the simplest case two identical monolayers are aligned and stacked in their natural ordering. For instance, MoS 2 bulk is found in the 2H phase, meaning that the unit cell consists of two layers, one being rotated by 180 • . In experiment several of such stacked samples arise each time a bulk crystal is cleaved. But there are many further phases possible, e.g. the 3R phase, in which the unit cell consists of three layers that are not rotated but translated relative to each other. Eventually, the layers can also be rotated with respect to each other. A further possibility is the stacking of two (or more) different materials. The resulting heterostructures may be very versatile, complex, and form Moiré patterns.

Bilayer MoS 2 in electric fields
In section 3, we have discussed the properties of a monolayer MoS 2 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 = ∞ (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.
To describe this mixing between the intralayer (B) and CT (A) excitons, we obtain a 4 × 4 matrix: On the diagonal, the exciton energies E 11 = E 22 = E B (intralayer B excitons) and E 12 = E 21 = E CT + 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 Å −1 . 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 VB1 → CB2, VB2 → CB1 is broken, and they change by the potential difference ±∆V, respectively. The inset in figure 17 shows the difference between both peaks. A shift of about 46 meV for 0.05 eV Å −1 is found for the mixed interlayer peaks.
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 2 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 × 25 µ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 × 0.6 nm ≈ 0.4 nm.
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 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 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.
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][127][128] we have investigated the characteristics in bulk MoS 2 , MoSe 2 , MoTe 2 , WS 2 , and WSe 2 . 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.

Heterobilayers
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 2 /WS 2 bilayer. The lattice constant of WS 2 is a WS2 lat = 3.155 Å [130] which is only 0.2% smaller compared to MoS 2 a MoS2 lat = 3.160 Å [131]. Here we employ a = 3.160 Å and neglect the small effects due to strain (see section 5). MoS 2 /WS 2 forms a type II bands alignment, i.e. the VBM stems from WS 2 , while the CBM is located in MoS 2 . Note that the strongly hybridized bands around Γ 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 WS2 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].
In addition to the neutral excitons, trions could also be present in the case of a heterobilayer. = 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 2 or both electrons on MoS 2 , respectively. The trion with one hole on both layersṪ 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 − IL further splits in three states (like for the intralayer trions), while T + 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 Ω T→ν = E T − ϵ ν , where ν can be different final states. For example,Ṫ IL,(1/2) is responsible for the peak discussed above, and another one is present slightly above the MoS 2 exciton (not shown).
Calman et al have experimentally measured interlayer trions in MoSe 2 /WSe 2 with a binding energy of 28 meV [135] and thus have validated our prediction.

Probing excitations in magnetic fields
In addition to the electric fields discussed in the last section, magnetic fields are frequently employed to probe optical properties [136]. In hydrogen atoms the response to magnetic fields is described by the normal and anomalous Zeeman effect for small fields and the diamagnetic effect for larger fields. Since excitons are conceptually similar to hydrogen, one would also expect their response to magnetic fields to be similar. While effects of a magnetic field on the optical absorption in semiconductors have been measured experimentally and modeled theoretically for various systems in previous decades, a predictive theoretical technique has been missing so far. In [75] we have worked out a first step towards the systematic determination of exciton g factors. Section 2.5 recapitulates the required theoretical framework, which we apply to 2D materials in this section.
We start our investigations for one of the well-known transition metal dichalcogenides namely MoSe 2 . The calculation of the magnetic moment m orb n ⃗ k in a periodic system is given by equation (39). Before we evaluate this, we take a look at the local limit Here χ ⃗ k α ′ µ ′ (⃗ r) are the lattice periodic basis functions which lead to a double sum over the lattice vectors and an integral over the unit cell. Because x and y appear in addition to the two localized functions, the common transformation to the relative distance (used for other matrix elements) is no longer valid. To check the importance of the double sum, we calculate ⟨L z ⟩ loc, ⃗ k αµ,α ′ µ ′ and thus m orb,loc n ⃗ k , which is plotted in figure 22. In panel (a), all values m orb,loc n ⃗ k are between ±2 µ B which is consistent with l ⩾ |m loc | for l ⩽ 2 (s, p, and d electrons). In the topmost valence bands and lowest conduction band close to K point, we find a value of m orb,loc VB,K = 1.93 µ B for both spin-orbit split bands and m orb,loc CB,K = −0.09 µ B for the conduction bands. These values can be understood by the character of the wave functions and the spherical harmonics. The VBM has contributions of 80% of Mo d orbitals with Y 2,±2 character (d x 2 −y 2 and d xy orbitals) and 18% of p x/y orbitals of Se. Close to the CBM, the contribution of Mo Y 2,0 (d z 2 orbital) is about 55% and the magnetic moment is much closer to zero.  band which resemble the main contribution of the bright A and B transitions (e.g., for MoS2 g A band = 2(mCB,K − mVB,K) and g B band = 2(mCB+1,K − m VB−1,K )), and the resulting g factors from equation (40). For comparison, several experimental measurements of the g factor of the A and B exciton are listed. Note that we expect an error of about ±0.05 µB for our magnetic moments m and about ±0.1 for the g factors.  [137], −3.0 [138], −3.8 [138], −4.0 [76], −4.2 [139], −4.6 [ [141,142], −4.1 [143], −4.2 [139], −4.3 [138], −4.4 [ [139], −3.94 [76], −4.0 [138], −4.25 [145], −4.35 [ [148], −3.7 [142], −3.8 [139], −4.37 [149] Going beyond the local approximation [68,69] The final step to calculate the optical response is to use the weights of the corresponding excitations in equation (40). To this end, we solve the Bethe-Salpeter equation to find the exciton weights on a mesh on ⃗ k points, similar to what we have previously discussed for MoS 2 in figure 4. Figure 23 shows the differences of the magnetic moments of the CB and the VB, and CB+1, respectively. Panels (a) and (b) are MoSe 2 , (c) and (d) show WS 2 . Due to the different ordering of the spin character in both materials, the values close to ±K are ±2 µ B for the CB in MoSe 2 but for the CB+1 in WS 2 . The contributions of the dark (spin-forbidden) state are about two time larger. The much smoother character for WS 2 stems from the missing band crossings in contrast to Mo-based materials.
The exact values are listed in table 1. Besides the magnetic moments at K and their differences g band (equation (41)), we list the g factors of the A and B excitons evaluated with equation (40). We find that g is more than 30% smaller compared to g band due to the spatial structure of the excitons. For instance, for MoSe 2 the interband transition exactly at K changes by the magnetic moment of −2.3 µ B (reversed at −K), and thus g A band = 2 · (−2.3) = −4.6. Away from ±K the absolute value of the magnetic moment decreases and thus g A is only −3.2. Overall, the theoretical and experimental g factors agree well (table 1), even if an exact comparison is difficult due to the varying experimental values. The g factor of the dark exciton is reduced from g D band = 2 · (−4.6) = −9.2 to g D = −7.4 respecting the excitonic nature. In addition to momentum direct excitons discussed previously, we can also use equation (40) to evaluate indirect excitons [86]. For example, we find g = 12.0 (g band = 15.6) for the lowest energy transitions for K→K' and g = −6.0 (g band = −8.6) for K→ Λ.

High-throughput exploration of hypothetical two-dimensional materials
Two-dimensional materials have emerged as a new class of materials in the last decade. In the last sections we have mostly discussed the properties of transition metal dichalcogenides and black phosphorus. However, since the observation of graphene [7] more than fifty monolayer compounds have been synthesised or exfoliated [3]. This asks for a more systematic study of stable 2D materials and their properties. In this section we will discuss a set of more than 4000 materials which we have studied and stored in the 'Computational 2D Materials Database (C2DB)' . To be able to compute such a large set of materials, a well-defined and automated workflow is essential. This high-throughput approach is implemented in python scripts employing the ASE and GPAW code [59,60].
In this large set of materials, we find a great variability of properties published in [3,38,150,151]. For instance, several TMDC show an anisotropy like ReSe 2 [152], several materials have a magnetic ground state like CrI 3 (which we have analysed in detail in the heterobilayer with WSe 2 [40]), or they belong to the class of Janus materials, i.e. compounds with different top and bottom layers, that have a built-in dipole moment. Further high-throughput studies suggest many more new materials [153].
After a short discussion of the workflow with special attention to GW calculations, we investigate the Janus material MoSSe in more detail. Many further results are accessible through the website https://cmrdb. fysik.dtu.dk/c2db.

Structure workflow
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 2 [10], WTe 2 [154], CrI 3 [155], TiS 3 [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 2 . 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.
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 × 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 2 and WTe 2 in figure 24).

Properties and statistics
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  . The structure workflow: starting from an initial guess, the relaxed structure is analysed and the stability is checked (see main text). Depending on the result, the candidate is stored in the database with or without further properties. Reproduced from [150]. © IOP Publishing Ltd. All rights reserved. Figure 26. Band gaps and gap centers (average of VBM and CBM) of more than 250 materials calculated with PBE, HSE, and GLLBSC (blue, yellow, and green dots). Reproduced from [3]. © IOP Publishing Ltd. All rights reserved. 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 0 W 0 results. The left panel shows the band gaps and indicates that PBE and HSE clearly underestimate the G 0 W 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 0 W 0 results. While we perform G 0 0W 0 and BSE calculations for semiconductors with less than six atoms per unit cell, all other properties are calculated using the PBE functional.

Automated GW calculations
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 ⃗ 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. O(N 4 ) where N is the number of electrons (compared to 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 2 . 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 = {1 − ∂Σ ∂ω } −1 . This behavior is observed for other materials as well. In general such cases can be identified by the condition Z ̸ ∈ [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 ̸ ∈ [0.5, 1.0] (MAE becomes 0.06 eV, not shown here).

Janus MoSSe
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 2 and MoSe 2 (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.
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 2 and 2.57 eV for MoSe 2 . This is also true for the spin-orbit splittings of the VBs and the CBs (210 meV and 30 meV), compared to MoS 2 (180 meV and 15 meV) and Band structure in the LDA (grey dashed) and the GdW(LDA) approximation (black). The band gap is 2.64 eV, and the VBM are set to zero. Note that the spin-orbit-split CBs close to K do not cross. (c) Optical absorption spectrum of MoSSe. The main peaks are labeled. The energetically lowest bright peak A is at 1.88 eV, and the first (dark) excitation is 25 meV below. Reprinted with permission from [151]. Copyright (2019) American Chemical Society. MoSe 2 (230 meV and 42 meV). Similarly, the lowest CB has opposite spin polarization compared to the highest VB. But unlike MoS 2 and MoSe 2 , 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 2 ) [82], we expect that the lowest excitation is dark in MoSSe similar to MoS 2 . 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.

Conclusions and outlook
In this review we summarize several investigations of the electronic and optical properties of 2D semiconductors which are obtained by employing and extending ab initio techniques. We examine several 2D materials and ways to design the optical excitation in them. In addition to exploring the excitations in the intrinsic, freestanding monolayer, we find potentially strong modification by the substrate on which the material is placed, by free carrier doping which leads to three-particle states, or by structural deformations when applying strain. In addition, defects and adsorbates can strongly alter the excitation spectrum, as well as the stacking of different 2D materials, which allows for excitations between the sub-systems, most notably interlayer excitons and trions. Furthermore, electric and magnetic fields can be applied to study and tune the response of the optical properties. The main computational methods applied are the DFT, the many-body perturbation theory in the GW and GdW approximation, and the Bethe-Salpeter equation which are extended by complementary developments, e.g. to capture the three-particle effects or the magnetic response of excitons.
The emphasis of the results presented is on the comparison to experiments. Due to the large number of measurements for the transition metal dichalcogenides, we have often focused one of the five materials MoS 2 , MoSe 2 , MoTe 2 , WS 2 , or WSe 2 in our calculations. For several effects like substrates, momentum dependent excitons, or the magnetic response our results provide new fundamental insight. For other topics such as excited state trions, interlayer excitons, or interlayer trions, our predictions have been verified by corresponding experiments. Many of the results obtained for the TMDCs are of general nature, i.e. related to the 2D nature rather the specifics of materials, and are thus expected to carry over to other 2D semiconductors.
The described methods in this review and their combination can also provide helpful ways for experimental studies of further physical effects. For example, the usage of electric fields is predestined to distinguish intra-and interlayer excitations. Doping or magnetic fields can provide further insights into the character of the involved quantum mechanical states.
In the last section, we have presented a few results from high-throughput studies of thousands of 2D materials. Besides the opto-electronic properties, a couple of systems have interesting topological [160] or magnetic properties [40]. For instance several materials show out-of-plane magnetization (like CrI 3 ), while other are polarized in-plane (and couple antiferromagnetically to the neighboring layers, like CrSBr [161]). In particular, the possibility to stack these layers promises novel physical phenomena and applications.
From a methodological point of view, ab-initio techniques are well-established for the investigation of opto-electronic properties. But there are still many open questions that need to be addressed. For instance, the introduced methodology to calculate trions can (in principle) be straightforwardly extended to evaluate four-or five-particle complexes, etc which have been reported in the literature [162,163]. To gain deeper insights into the lifetimes and dynamic processes of excited states, the coupling of electronic and vibronic degrees of freedom and the consideration of the temperature is essential [164][165][166]. Their careful treatment is a difficult task which is still beyond state-of-the-art.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.