Calculating excitons, plasmons, and quasiparticles in 2D materials and van der Waals heterostructures

Atomically thin two-dimensional (2D) materials host a rich set of electronic states that differ substantially from those of their bulk counterparts due to quantum confinement and enhanced many-body effects. This Topical Review focuses on the theory and computation of excitons, plasmons and quasiparticle band structures in 2D materials and their heterostructures. The general theory is illustrated by applications to various types of 2D materials including transition metal dichalcogenides, graphene, phosphorene, and hexagonal boron nitride. The weak and highly non-local dielectric function of atomically thin crystals is shown to be responsible for many of the unique properties exhibited by the 2D materials such as the formation of strongly bound, non-Hydrogenic excitons, large band gap renormalization effects, and the different signatures of excitons and plasmons in electron energy loss spectroscopy (EELS). Among other topics covered are spin–orbit coupling, trions, interlayer excitons, exciton dissociation, and environmental screening. Technical issues associated with the application of the many-body GW method and the Bethe–Salpeter equation (BSE) to 2D materials are also discussed. A combined quantum/classical method is introduced and used throughout to account for dielectric screening and self-energy effects from substrates and van der Waals heterostructures including the difficult case of non-matching lattices.

2 K S Thygesen sensitive to screening from the environment, e.g. substrates [29,30]. This state of affairs makes it essential to go beyond the single-particle approximation and explicitly account for many-body effects when modelling the electronic and optical properties of 2D materials.
Although the 2D materials are interesting in their own right, an even larger potential lies in the possibility of reassembling different 2D crystals into new layered compounds [31]. Such designer materials have been coined van der Waals heterostructures (vdWHs) with reference to the weak van der Waals forces holding the 2D crystal planes together. The concept provides an ideal platform for controlling the unique 2D electronic states with atomic scale precision, and offers an enormous amount of possibilities to broaden the versatility of 2D materials and to achieve superior and unusual material properties. In contrast to conventional heterostructure fabrication methods, which involve complex and expensive ultrahigh-vacuum techniques to epitaxially grow the single-crystalline layers, vdWHs can be stacked in ambient conditions with no requirements of lattice matching. Furthermore, the weak interlayer bonding in vdWHs leads to clean and atomically well defined interfaces, thus reducing detrimental interface scattering effects [32,33].
The vast size of the vdWH compound space can be illustrated by considering a heterostructure containing 20 atomic layers (which is similar in size to the most complex structures made to date [34]). If each layer is picked from a pool of ten different 2D mat erials this yields 10 20 different stacking combinations. In reality, the space of presently available 2D mat erials comprises several hundreds of compounds [35] and the possibilities for combining them are essentially endless. To harness the full potential of vdWHs it is essential to develop efficient and reliable computational methods that can scan the vdWH compound space and guide experimentalists towards the most interesting materials before undertaking the time consuming task of synthesising specific heterostructures.
The computational modelling of vdWHs, or substrate supported 2D materials, is complicated by the incommensurable nature of the interfaces. Solid state computational codes apply periodic boundary conditions which imply the use of large super cells to avoid straining the different 2D layers of the heterostructure. However, such an approach quickly becomes prohibitive for anything but the simplest few-layer structures. Consequently, reliable modelling of realistic, incommensurable heterostructures requires novel approaches that combine the quantum description of the individual layers with a more coarse grained description of the effect of interlayer interactions. Recently, a few such methods have appeared addressing different aspects of the interactions including the effect of interlayer hybridisation [36] and di electric screening [37].
This review focuses on the theory and computational methods used to describe elementary electronic excita-tions in 2D materials and heterostructures. To limit the scope, only pristine crystals are considered, although defects are known to play an important role for the chemical and physical properties of 2D materials [38][39][40]. The review is organised into four main sections devoted to the topics: dielectric screening (section 2), quasiparticle band structures (section 3), excitons (section 4), and plasmons (section 5). Each section is structured such that the first part focusses on isolated 2D materials while the last part extends the discussion to van der Waals heterostructures and effects of substrate interactions with an emphasis on the role of the di electric screening.

Dielectric screening
The macroscopic dielectric constant of a conventional bulk material is defined as the limiting value of ( ) ω ε q q, as → q q 0 0. For a 2D material this definition cannot be straightforwardly adopted since ( ) . In fact, the dielectric function of a 2D material is strongly q q-dependent and a more elaborate treatment of screening is required [41]. In the simplest approximation, the 2D material can be modelled as an infinitely thin isotropic slab. This approximation leads to a dielectric function of the form where α is the 2D polarizability of the layer [42]. Using that the 2D Fourier transform of the 1/r Coulomb interaction equals / π q 2 , the screened interaction can be written 2D 2D (2) This simple form of the screened interaction has the merit of leading to an analytical expression in real space and has been used by several authors to study excitons and charged impurities in various 2D semiconductors [23,[42][43][44]. We note that the form / π q 2 for the Coulomb interaction and equation (1) for the dielectric function are consistent approximations which both become exact in the limit of vanishing thickness of the material, i.e. the strict 2D limit. In this section, the strict 2D screening model is first compared to more realistic models taking the finite thickness of the material into account. Next, we consider the effects of support or van der Waals heterostructuring on the dielectric properties.

Quasi-2D Coulomb interaction
To understand the nature of dielectric screening in atomically thin materials, we begin by considering the form of the bare Coulomb interaction. It is instructive to compare the case of a strict 2D material to the physically relevant case of a 2D material with a finite thickness. To this end we consider the interaction energy between two normalized line charges, ρ 1 and ρ 2 , located at the in-plane positions ∥ r r ,1 and ∥ r r ,2 , respectively, and extending over the 2D material in the out-of-plane direction (see figure 1(b)), r r r r r r r r r r r r r r r r , dd .
The (2D) Fourier transform of the quasi-2D interaction becomes Physically, V q q Q2D ( ) ∥ equals the Coulomb energy of a quasi-2D charge density wave with wave vector ∥ q q and thickness d. We note that for in-plane wavelengths much larger than the thickness of the material ( ) equals the strict 2D potential while for smaller wave lengths it equals the 3D potential, i.e.

Dielectric function: 2D versus 3D
For a material of any dimensionality, the (inverse) microscopic dielectric function gives the total potential to first order in the applied external potential, r r r r r r r r d , , tot 1 ext (6) where we have neglected the frequency dependence for ease of notation. For periodic systems it is often more convenient to use a plane wave representation. Taking q q to represent a wave vector in the first Brillouin zone (BZ) and G G a reciprocal lattice vector, the dielectric matrix of a crystal evaluated within the random phase approximation (RPA), takes the form where ( ) + V q q G G is the 3D Fourier transform of the Coulomb potential and χ 0 the non-interacting response function, r r e , r q q q q tot i (9) where ˜( ) V r r q q is a lattice periodic function. Since usually we are interested in macroscopic fields, we define the 3D macroscopic dielectric function as where ⟨ ⟩ Ω ... denotes a spatial average over a unit cell.
due to local field effects [46].
First-principles codes developed for bulk crystals use super cells with a finite vacuum region of thickness L to separate the periodically repeated 2D layers. When such codes are used direct application of equation (10) leads to ( ) . This is a consequence of an averaging region much larger than the effective extension of the electron density around the material. The standard definition in equation (10) becomes meaningless in this case, which is the reason why relatively different values for ε D 3 have been reported for monolayer materials in the recent literature [47][48][49].
From the first equality in equation (10), it is natural to substitute an average over the entire unit cell in the out of plane direction with an average over a confined region corresponding to the actual extension of the electronic density. In practice, we average the in-plane coordinates ( ) ∥ r r over the unit cell area, A, and the outof-plane coordinate from z 0 − d/2 to z 0 + d/2, where z 0 denotes the center of the material and d its width. The macroscopic dielectric function then becomes: according to the RPA expression in equation (7).
We stress that it is essential to use a truncated Coulomb potential in equation (7) in order to decouple the layers in neighboring supercells [41]. Note that we used the label Q2D in equation (11) since this definition of the macroscopic dielectric function is consistent with the quasi-2D picture. As a rule of thumb we choose d to be the distance between the layers in the layered bulk, but the results are not very sensitive to this choice [41].  The calculated ∥ q q -dependent (static) dielectric functions of monolayer hBN and MoS 2 are illustrated in figure 2. The dielectric function approaches unity for → ∥ q q 0 0 meaning that long range screening is completely absent. Mathematically, this follows from equation (7) and the fact that in the small q q limit, ( ) χ ω ∝ q q q, 0 00 0 0 2 while the truncated Coulomb interaction, / ∝ V q 1 trunc . In addition to the Q2D dielectric function, figure 2 shows the linear expression in equation (1) as well as the di electric function of the 3D layered crystals. As was the case for the bare quasi-2D Coulomb interaction, see equation (5), the quasi-2D dielectric function follows the strict 2D result for ∥ q d 1 and the bulk result for ∥ q d 1.

Dielectric screening in van der Waals heterostructures
The incommensurate nature of vdWHs presents a great challenge for first-principles calculations because it is generally not possible to represent the heterostructure in a computational cell that is small enough to allow the calculation to be performed without straining one or more of the layers and thereby alter its electronic properties. The problem is particularly severe for manybody calculations for which the computational cost grows rapidly with system size. For example, the heterobilayer MoS 2 /WSe 2 , which is studied in greater detail in section 4.6, requires a unit cell of at least 100 atoms to reduce the strain to 1% rendering standard plane wave based many-body calculations highly challenging even for this simple system. To overcome the problem, the quantum electrostatic heterostructure (QEH) model computes the dielectric function of a general vdWHs from the ab initio response functions of the constituent 2D layers coupled only via their electrostatic interactions, see figure 3. A minimal representation of the response functions of the 2D layers (the dielectric building blocks), consisting of the ( ) ∥ ω q q , -resolved monopole and dipole components (out of plane), is employed without compromising the accuracy. As shown in later sections, the QEH model allows us to evaluate the screened interaction entering the theories of excitons, plasmons, and quasiparticle band structures, for general incommensurate vdWHs. A database containing the pre-calculated dielectric building blocks of a collection of 2D materials is available on-line (https://cmr. fysik.dtu.dk). From here the data files can be downloaded together with a Python module for constructing the dielectric building block of any 2D mat erial and solving the coupled electrostatic equations for the heterostructure dielectric function. QEH model calculations for vdWHs containing a few hundred layers can be performed on a standard PC. Details on the QEH model can be found in [37].
As a simple illustration of the QEH model, figure 4 shows the calculated static dielectric function of multilayer MoS 2 slabs of varying thickness as a function of in-plane wave vector. The dielectric function measures the ratio between an applied potential of the form ⋅ V r r q q exp i 0 ( ) ∥ and the resulting total potential averaged over the slab, i.e. the macroscopic dielectric function of the multilayer as defined in equation (11). For large ∥ q the dielectric functions show similar behavior. However, whereas ( ) = ε 0 14 for the bulk, the dielectric functions of the slabs decrease sharply to 1 for small ∥ q . This demonstrates that the dielectric properties of a vdWH of thickness d are 2D-like for / ∥ q d 1 and 3D-like for / ∥ q d 1 . The QEH model describes the change in the dielectric function from monolayer to bilayer very accurately in spite of the well known differences between their electronic band structures [47]. This shows that hybridisation driven band structure effects, i.e. quantum confinement, have rather weak influence on the dielectric properties of a vdWH and is the main reason for the success of the QEH model.

Quasiparticle band structures
The electronic band structure of a 2D material can be routinely calculated using density functional theory (DFT) adopting one of the many existing approximations for the exchange-correlation (xc) functional. It is well known that semi-local xc-functionals underestimate band gaps of semiconductors (including 2D semi-conductors [50]) while hybrid functionals, such as the HSE [51], perform better albeit at higher computational cost. It is perhaps less well known that the single-particle DFT band structure, even when calculated with the exact xc-functional, is not supposed to be exact. The exact band gap can in principle be obtained by adding the socalled derivative discontinuity, ∆ xc , to the exact DFT single-particle gap, however, most functionals do not evaluate ∆ xc or incorrectly predict it to be zero. One exception is the GLLBSC functional, which has been shown to predict fairly accurate band gaps for a range of materials [52].

GW approximation
The true band energies of a solid can be obtained as the poles of the single-particle Green function. For weakly correlated materials, these poles are well defined, i.e. the so-called renormalization factor Z is close to unity, and can be obtained by solving the quasiparticle (QP) eigenvalue equation where the xc-potential of the DFT Hamiltonian is replaced by the non-local and energydependent self-energy operator, ( ) ω Σ ′ r r r r , , . Physically, the QP energies describe the energy cost of adding electrons/holes to the neutral ground state of the material. They are thus defined in terms of the manybody total energies of the neutral and charged systems, ,for empty states The QP energies can be measured by photo-emission or transport spectroscopy which probe the transition energies between states with N and ± N 1 electrons. The QP gap is defined as the difference between the highest occupied state (ionisation potential) and lowest unoccupied state (electron affinity). It differs from the optical gap by the exciton binding energy which can be substantial for atomically thin materials, see section 4.
A highly successful approach to predict QP energies is the so-called GW method, which approximates Σ to first order in the screened interaction, W. The GW approximation was introduced by Hedin [53] in 1965 and first applied to real solids in an ab initio framework by Hybertsen and Louie [54] and Godby et al [55]. Compared to mean field methods, like the DFT Kohn-Sham scheme, the GW method is significantly more involved both in terms of numerical implementation and the required computation time. In real space the GW self-energy takes the form where G is the single-particle Green function. Within the simplest G 0 W 0 approach, the Green function and  [45].
To obtain the QP energies one treats Σ − V xc using firstorder perturbation theory which leads to the expression where Z nk k is the renormalisation factor. For more details and background theory on the GW method we refer to the extensive literature on the topic.

GW calculations of freestanding 2D materials
The performance of the GW approximation has been thoroughly established for bulk materials [56][57][58] and more recently also for molecules [59,60]. In comparison, systematic studies assessing the accuracy and numerical convergence of GW calculations for 2D materials are still few. Nevertheless, it has become clear that (i) it is essential to use a truncated Coulomb interaction when evaluating W in order to avoid long range screening between periodically repeated layers which otherwise leads to an unphysical reduction of the QP band gap, and (ii) when a truncated Coulomb interaction is used, the convergence of the GW calculation with respect to the number of in-plane k-points becomes much slower than is the case for similar bulk systems [41,50,61,62].
The slow k k-point convergence of the GW band structure is a direct consequence of the strong wave vector-dependence of the dielectric function of a 2D semiconductor, see figure 2, which implies that the number of q q-points required to obtain a proper sampling of the screened interaction ( ) W q q over the BZ is much higher for the 2D material than what would be anticipated from the 3D case. For example, the band gap of bulk 2H-MoS 2 is converged to within 0.1 eV with an in-plane k k-point grid of × 12 12 while the same accuracy for monolayer MoS 2 requires a grid of × 40 40 when standard BZ sampling schemes are applied.
The k k-point convergence can be greatly improved by using an analytical expression for the small q q-limit of the 2D response function when performing the BZ integral over the critical region around = q q 0 0 [63]. In brief, an expression for the dielectric function similar to equation (1) is used to obtain a closed expression for ( ) ′ W q q G GG G which can be analytically averaged around = q q 0 0. As seen in figure 6(a) the analytical expression for ( ) W q q 0 00 0 is valid for surprisingly large in-plane q q-vectors. The use of the analytical model drastically reduces the requirements on the k k-point mesh. Figure 6(b) shows that the band gap of MoS 2 converges much faster than the traditional sampling scheme. The latter underestimates the screening around = q q 0 0 because ( ) = = ε q q 0 0 1, and consequently no screening is picked up at this sampling point which results in a large overestimation of the gap. Table 1 shows the calculated QP band gaps of mono layer hBN, MoS 2 , and phosphorene. These calculations were performed using a truncated Coulomb interaction, extrapolated as /N 1 G G to the infinite plane wave basis limit, and converged with respect to k k-point grid. Spin-orbit interactions were not included in the reported values. Inclusion of spin-orbit interactions has no effect for hBN and phosphorene, but split the valence band of MoS 2 at the K point by 0.15 eV thereby lowering the QP gap by around 0.07 eV [50,64] (see  table 2 for the spin orbit splittings in the TMDs). For MoS 2 the converged G 0 W 0 @PBE band gap of 2.54 eV agrees well with our previously reported value of 2.48 eV (with spin-orbit coupling) obtained using a Wigner-Seitz truncated Coulomb interaction and × 36 36 k k-points [50]. Other reported G 0 W 0 gaps range from 2.40-2.82 eV [47,48,49,65,66,67]. However, these calcul ations were performed (i) without the use of a truncated Coulomb interaction (ii) using relatively small k k-point grids of × 6 6 to × 18 18, and (iii) using different in-plane lattice constants vary-ing between 3.15 and 3.19 Å. These different settings can affect the band gap by as much as 0.5 eV [41], and therefore we refrain from providing a detailed comparison of our result to these earlier calculations. An overview of previous GW results for MoS 2 can be found in [41,62].
For many applications, not only the distance between the occupied and unoccupied bands, i.e. the band gap, but also the absolute position of the band edges relative to vacuum are of interest. For a 2D mat erial the abso- lute band energies are obtained by referring them to the asymptotic value of the Hartree potential in the vacuum region. For bulk materials, the problem is much more delicate as it requires the use of thick slabs to represent both the bulk interior and its surface, and because the Hartree potential depends on the surface dipoles (on both sides of the slab), which makes the band alignment dependent on the surface termination. These problems are obviously not present for the 2D materials making them ideal as benchmark systems for the band alignment problem. Figure 5 shows the positions of the valence band maximum (VBM) and conduction band minimum (CBM) relative to the vacuum level for 51 MX 2 monolayer metal-dichalcogenides and -oxides evaluated at both the LDA and G 0 W 0 levels. For all materials, the effect of the G 0 W 0 correction is to shift the conduction band up and the valence band down with respect to the LDA values. There are clear correlations between the band energies and chemical composition of the materials with the oxides showing the largest gaps and deepest lying band edges followed by the sulphides, selenides, and tellurides in that order. These trends can be explained in terms of the relative size and electronegativity of the O, S, Se, and Te atoms [50]. In section 4.6, the band alignment at the MoS 2 /WSe 2 heterobilayer is analysed in greater detailed.

Spin-orbit interactions
Several of the interesting 2D materials, in particular the TMDs and the novel class of 2D metal-organic hybrid perovskites [20,21], contain heavy elements for which spin-orbit (SO) interactions can have significant effects on the band structure. The SO coupling enters the Hamiltonian via the term [68] ( ) where c is the speed of light, m the electron mass, V the electron potential, and S S and L L are the electron spinand angular momentum operators, respectively. The above form assumes a spherical potential. Even in a crystal lattice this is usually an excellent approximation because the SO coupling is completely dominated by the spherically symmetric regions close to the nuclei.
To compute the effect of SO coupling on the band structure, the single-particle Hamiltonian (including the H SO term) must be diagonalised within the space of spin-up and spin-down eigenstates, { } ψ σ nk k, . Table 2 shows the spin-orbit induced splittings at the valence band maximum and conduction band minimum for some monolayer TMDs and TMOs. The splitting has been calculated by diagonalizing the LDA or the G 0 W 0 @ LDA Hamiltonian, respectively. Clearly, the SO splitting is rather insensitive to the underlying band structures (LDA or GW) with large large differences occurring only in the cases where LDA and GW predict different locations of the band extrema within the BZ. Note that only monolayer TMDs in the 2H phase appear in table 2. This is because the combination of inversion symmetry and time reversal symmetry in the 1T phase inhibits spin splitting of the bands. The absence of band splitting does, however, not imply that SO coupling can be neglected for the 1T structures.
For the MoX 2 and WX 2 monolayers in the 2H-phase, the SO-induced splitting of the valence band at the K and ′ K points of the BZ, leads to a splitting of the lowest exciton into two excitons with opposite spin (the A and B excitons in figure 9). Since time reversal symmetry maps ⟩ σ | K , to σ |− ′ K , 〉, the energetic ordering of the spin states is reversed in the K and ′ K valleys. This makes it possible to selectively excite excitons in one of the two valleys using circularly polarised light tuned to the Table 1. Calculated band gaps (in eV) from PBE, G 0 W 0 @PBE and partially self-consistent GW 0 @PBE for three monolayers in PBE-relaxed structures. The GW calculations were performed using a truncated Coulomb interaction and analytic integration of ) ( W q q around = q q 0 0. Spin-orbit interactions were not included. The k k-point grids were × 18 18 (h-BN), × 18 18 (MoS 2 ) and × 10 14 (phosphorene).  Table 2. Spin-orbit induced splittings at the valence band maximum and conduction band minimum for some monolayer transition metal dichalcogenides and -oxides. The band splittings have been calculated by full diagonalization of the LDA and G 0 W 0 @ LDA Hamiltonians, respectively [73]. The location of the band extrema in the BZ is indicated in parenthesis. Note that this can be different in LDA and GW. Table reproduced from [50].
frequency of the A or B exciton [69][70][71]. These discoveries have recently fueled the idea of using monolayer TMDs for valleytronics [72], i.e. technologies that exploit control over the valley degree of freedom.

Environmental screening
As described in the previous section, the QP energies depend on the screened interaction, = − ε W V 1 , via the electron self-energy. Owing to the long range nature of the Coulomb interaction and the weak intrinsic screening in 2D materials, the screened interaction, and consequently the QP energies, can be strongly affected by the dielectric environment outside the 2D material. For example, it was shown using non-linear optical spectroscopy and many-body calculations that the band gap of monolayer WS 2 can vary by 0.3 eV depending on the choice of substrate [29]. This effect can be utilised to control the band structure of 2D materials via the nonlocal screening by the substrate/embedding as demonstrated by GW calculations for both lateral [74] and vertical [75] heterostruture designs. Figure 7 shows the band gap of an hBN layer as a function of distance, L, to a single graphene sheet. The G 0 W 0 gap follows a 1/L dependence and at the equilibrium distance around L = 3.5 Å, the band gap is 1 eV lower than that of an isolated hBN layer despite the fact there is essentially no overlap between the wave functions of the two systems at the considered separations. This phenomenon is known as the image charge effect. It has been studied extensively for molecules on metal surfaces, where it can shift molecular energy levels by several electron volts [76,77]. For metal-moleculemetal junctions a proper incorporation of the image charge effect is essential for correct prediction of energy level alignment and transport properties [78,79].
In the case of a molecule on a metal surface, the origin of the image charge effect is easy to understand. Indeed, according to equations (12) and (13), the QP energies represent electron addition/removal energies.
When an electron is added or removed, the molecule is left in a ±1 charge state. As a result of the electrostatic interaction between the charged molecule and its image charge in the metal surface, the QP energies will be lowered relative to those of the gas-phase molecule. Due to the different sign of occupied and unoccupied states, see equations (12) and (13), the image charge effect shifts occupied (empty) states up (down) in energy thereby closing the gap. For a 2D material, one might expect the image charge effect to vanish because the electrostatic interaction between a completely delocalised charge distribution of a Bloch state and its equally delocalised image charge, scales as 1/A where A is the area of the 2D material. However, this picture is in fact incorrect because the electron is not delocalised -only its probability density. Thus the image charge effect is a pure correlation effect.
Returning to figure 7, it can be seen that the LDA, in addition to underestimating the hBN band gap by 3 eV, does not capture the image charge effect. This is a common failure of all semi-local and hybrid xc-functionals [77] which are unable to detect changes in the singleparticle energies when the (occupied) wave functions of the two subsystems do not overlap. This shows that quantitatively accurate modelling of band structures of supported 2D materials or vdWHs must employ manybody methods that are able to capture the non-local screening effects.
It is interesting to note that the effect of environmental screening on the QP band gap is highly mat erial dependent. For example, GW calculations predict a reduction of the QP gap of a benzene molecule, an hBN monolayer, and an MoS 2 monolayer when placed on a graphite surface to be 3 eV, 1 eV, and 0.2 eV, respectively [30,76,80]. These differences cannot be explained by a simple image charge model that only accounts for the interaction between the bare electron/hole in the 2D material and its image charge in the substrate. Indeed, such a model would yield the same reduction of the QP gap for all three systems. The reason is that the change in QP band gap is determined by the relative change in screened interaction induced by the environment. Since the dielectric function of the 2D material has the form πα = + q q q 1 2 ( ) ε one expects that materials with larger polarizability (α) are less affected by the environ mental screening. Indeed, since the internal screening increases from benzene, to hBN to MoS 2 , this argument provides an explanation for the observed differences in QP gap reduction.
The dependence of the band gap reduction on the intrinsic polarisability of the 2D material is clearly illustrated by figure 8 which shows the change in band gap (∆E gap ) for six 2D semiconductors (SC) in four different heterostructure configurations: on top of a single graphene sheet, sandwiched between two graphene sheets, and the same configurations with graphene replaced by the metallic 1T phase of MoS 2 . The change in band gap due to the environmental screening has been calculated using the G∆W method. In this method, the screened potential in the 2D material (layer i of a heterostructure (HS)) is split into two parts, i i i ,HS (18) where W i is the screened potential of the freestanding monolayer and ∆W i is the change due to the surrounding layers. The latter term can be evaluated efficiently using the QEH method described in section 2.3. By combining equation (18) with equations (14) and (16) and assuming the Green function of the embedded layer to be the same as that of the freestanding monolayer (which is a good approximation due to the weak interlayer hybridisation), the change in QP energy can be readily computed. For details on the G∆W method see [75]. Returning to figure 8, it is seen that the band gap correction is larger for the 1T-MoS 2 substrate than for graphene reflecting the increased screening from a metallic compared to a semi-metallic substrate. Moreover, as anticipated, the magnitude of the band gap correction is seen to scale with the inverse (static) polarisability of the 2D semiconductor.
The results of this section show that dielectric engineering can be used to manipulate the band gap of a 2D material and ultimately the vertical band structure profile of a vdWH. As another interesting application, it has been proposed that lateral control over the band edges within a homogeneous 2D material may be achieved by using a laterally structured dielectric environment [74].

Excitons
One of the hallmarks of the 2D semiconductors and insulators is the existence of strongly bound excitons with binding energies reaching up to 30% of the band gap. These excitons couple strongly to light and lead to a substantial modification of the optical spectrum both below and above the QP band gap. The large binding energies of excitons in 2D materials was initially predicted theoretically for monolayer hBN [81], graphane [23] and TMD monolayers [41,82,83], and has subsequently been confirmed experimentally [29,[84][85][86]. The reason for the large exciton binding energies in 2D materials is a combination of geometric confinement [87] and reduced dielectric screening [23,41,45]. The latter implies that exciton binding energies in 2D materials can be controlled via the dielectric environment, e.g. by encapsulation or deposition on substrates. Furthermore, by stacking 2D semiconductors with staggered band alignment, interlayer excitons with the electron and hole residing on different layers, can be formed. Such heterostructures provide an ideal platform for studying charge separation-and recombination processes in atomically well defined crystal environments.

Excitons in freestanding 2D materials
The most accurate method for calculating interband excitations in semiconductors and insulators is the Bethe-Salpeter equation (BSE) based on many-body perturbation theory [46,88]. The excitation energies of crystal momentum q q are found by solving an eigenvalue problem of the form where ( ) H q q is the BSE two-particle Hamiltonian eva luated in a basis of e-h states ψ = r r r r , The kernel contains two terms, namely the electronhole exchange interaction (V ) and the direct screened electron-hole interaction (W ), The factor 2 accounts for spin. Neglecting the direct interaction (W) results in the random phase approximation (RPA). For a 2D material it can be shown that the exchange term becomes negligible in the small q q limit (see section 5.6). Thus the e-h exchange only reduces the binding energy of the lowest zero momentum exciton of monolayer hBN and MoS 2 by 0.08 eV and 0.02 eV, respectively, which amounts to less than 5% of the total binding energy.
To illustrate the important role of the direct e-h interaction for the optical properties of 2D materials, figure 9 shows the absorbance of a single layer of MoS 2 calculated from the BSE with and without inclusion of W. Since the e-h exchange term is negligible in the optical limit ( = q q 0 0), the RPA spectrum is essentially equal to the non-interacting result consisting of the bare e-h transitions. Inclusion of the direct e-h term changes the spectrum completely, introducing distinct peaks below and just above the QP band gaps due to absorption by excitons. The A and B excitons (and their Rydberg states A' and B') are composed of transitions from the spin-orbit split valence band to the conduction band at the K-point of the BZ. The C exciton involves transitions around the Γ-point [83]. The spectrum has been numerically broadened to simulate the effect of coupling to phonons at finite temperatures.
In the following we specialise to = q q 0 0 and thus neglect ′ V SS . Furthermore, we specialise to a single valence and conduction band, i.e.
( / ) = S v c k k , . This two-band approximation is typically excellent for the lowest bound excitons as long as the valence and conduction bands are separated from other bands at the relevant points of the BZ. Assuming that the wave functions are of the form ( ) , the direct screened interaction becomes, In analogy to the 3D case it can be shown for a 2D semiconductor that the so-called exciton envelope function, ( ) ( ) F k k are the BSE expansion coefficients in equation (19), satisfies the 2D Schrödinger equation Here E B is the exciton binding energy given as the difference between the QP band gap and the lowest optical transition, and µ ex is the exciton effective mass, calculated from the hole and electron masses according to  probability amplitude for the distance between the electron and hole. The screened interaction can be evaluated at different levels of approximations. First, the strict 2D approximation with the screened potential in equation (2) admits a closed analytical expression in real space where H 0 (x) and N 0 (x) are the Struve and Neumann functions, respectively. Next, taking the finite width of the electron and hole charge densities into account by approximating the out-of-plane wave functions by step functions yields the quasi-2D screened potential [45] ( ) ( ) ( ) where ε Q2D was defined in equation (11). Finally, the out-of-plane variation of the electron and hole wave functions can be calculated from the exact wave functions at the relevant point of the BZ, . However, as shown below, the exciton binding energies obtained from the effective Hamiltonian equation (23) are rather insensitive to the precise form chosen for the out-of-plane electron and hole densities. Table 3 shows the binding energy of the lowest exciton in monolayer MoS 2 and hBN obtained from equation (23) with W evaluated in three different ways, namely (i) from equation (22) with / φ v c given by the exact wave functions at the direct band gap (K point for MoS 2 and Γ point for hBN) (ii) from equation (25) corre sponding to / φ v c being step functions, and (iii) equation (24) corre sponding to a strict 2D Coulomb potential. It can be seen that all the results obtained from the effective Hamiltonian agree within a few percent. The reason for this is that the Fourier transformed exciton wave function, ( ) ∥ F k k , is localised in the small ∥ k k regime where the linear approximation, ( ) , is valid. For comparison, the results obtained from the BSE are listed in the first column. It should be stressed that BSE calcul ations for 2D semiconductors, performed with a truncated Coulomb interaction to avoid screening between periodically repeated layers, converge very slowly with respect to the number of k k-points, see figure 10. The reason is again the strong variation of the screening around = q q 0 0 as discussed in detail in section 3 in the context of GW calculations [41,62].

Screened 2D hydrogenic model
Apart from the dimensionality, the main difference between the hydrogen model of 3D excitons and the effective Hamiltonian equation (23) is that the screened interaction entering the latter is not simply /εr 1 . This implies that a closed expression for the exciton energies and radii cannot be directly obtained in the 2D case. However, as shown below it is possible to define an effective 2D dielectric constant and thereby recover a 2D Hydrogenic model which admits analytical solutions. The effective dielectric constant is obtained by averaging ( ) ε q q over the q q-vectors up to /a 1 where a is the radius of the exciton in real space, q q q q , , d The rationale behind this definition is that the q q-vectors coming into play when evaluating the average interaction energy ( ) ( ) , are mainly those fulfilling | | < a q q 1/ . For the 2D Hydrogen atom the Bohr radius is given by /( ) µ = ε a 2 ex and equation (26) has to be solved self-consistently for ε given an expression for ) ( ε q q . In a strictly 2D system the screening is linear in q q and equation (26) can be solved to yield [90] Using that the Hydrogenic binding energy in 2D is a factor of four larger than in 3D [87], we obtain This expression constitutes a long-sought-for 2D analog of the 3D Hydrogenic exciton model. A remarkable property of the expression (28) is that it becomes independent of the effective mass if the polarizability is sufficiently large. More precisely It may come as a surprise that the binding energy becomes independent of mass, since a large mass gives rise to a localized exciton and the binding energy typically increases with localization. However, in 2D, short range interactions are screened more effectively than long range interactions. Thus, there are two opposing effects of the exciton mass and for large polarizabilities, the binding energy becomes independent of mass. In order to assert the applicability of the expressions (28) and (29), we have calculated the effective masses and static polarizabilities (in the random phase approximation) of the 51 semiconducting TMDs shown in figure 5. In figure 11 we compare the model binding energies with the full solution of equation (23). Using the expression (28), the agreement is seen to be on the order of 10%. With the  (29), we also obtain excellent agreement for binding energies up to ∼0.5 eV, whereas the binding energies are somewhat underestimated for the more strongly bound excitons. Recently, first principles calculations have indicated that exciton binding energies in 2D materials scale linearly with the band gaps [91]. In the present model, this behavior comes out naturally since the in-plane component of the polarizability, α, is roughly inversely proportional to the band gap as shown in [90]. Combining this with equation (29) thus gap . However, in the present model the scaling originates solely from the screening and not the effective mass as previously proposed [91]. In fact, there is no direct correlation between the exciton binding energies and effective masses for this class of materials. The screened Hydrogen model is not restricted to the lowest excitonic state, but can be used to account for the entire exciton series in 2D materials. In [92], the exciton series of graphene derivatives was predicted to deviate from the 2D Rydberg series and in [93] this was experimentally demonstrated for a monolayer of WS 2 on SiO 2 . The reason is that the effective screening depends on the n quantum number due to the increasing spatial extent of higher lying excitons. In [93] the authors defined an n-dependent effective dielectric constant, ε n , which was determined by fitting each term in the Rydberg series to a 2D Hydrogen model. The Rydberg series is thus written We have previously shown that the Rydberg series of monolayer WS 2 can be accurately reproduced by the effective Hamiltonian equation (23) [37]. Thus the full solution of equation (23) is used as a reference for the Hydrogenic model in the following. We calculate the n-dependent effective screening by replacing a in equation (26) by an n-dependent characteristic extension of the state. To this end we note that for l = 0, where ˆˆ= + r x y 2 2 . With these more general definitions, the a defined previously is given by , the effective dielectric constant for state n then becomes Using a 2D polarizability of α = 5.25 Å and µ = m 0.19 ex 0, both obtained from first principles, we obtain the Rydberg series of a freestanding WS 2 monolayer shown in figure 12. Clearly the analytical results from the screened Hydrogenic model agree very well with the full solution of equation (23). In contrast, the pure 2D Hydrogen model with a fixed state-independent dielectric constant significantly underestimates the binding energies of higher lying states, since the decreased screening of extended states is not taken into account. We also note that the model binding energies of the n = 1 state agree very well with a full solution of the Bethe-Salpeter equation which yields an exciton binding energy of 0.54 eV [66]. The effective dielectric constants from equation (32) are shown in the inset and are in excellent agreement with the values fitted to experiments [93] (not shown).
It is straightforward to generalize the above expressions to ≠ l 0 [87], which results in a larger value of the effective radius, a nl , and thus < > ′ ′ l l for n l n l , , ε ε . The energy is still given by equation (30) and for a given n, the higher angular momentum excitons will therefore have a larger binding energy. This effect was indeed observed experimentally for WS 2 monolayers [84].
Finally, it should be noted that, in addition to the nonlocal screening captured by the effective Hamiltonian equation (23) and the screened Hydrogen model, there are a number of other effects that can lead to deviations from the Hydrogenic exciton series. These include the break down of the effective mass approx imation, the two-band assumption as well as non-zero Berry curvature of the Bloch bands [94,95]. All these effects are, however, related to the properties of the band structure rather than the dimensionality of the material and thus are not genuine 2D effects.

Finite momentum excitons
The above discussion has focussed on excitons with zero center of mass momentum, i.e. = q q 0 0. However, excitons with finite momentum, although optically inaccessible, are essential to describe exciton dynamics and the response of the electron system to momentum carrying probes such as x-rays or fast charged particles [96,97]. In fact, the exciton q q-dispersion can be readily accessed by inelastic x-ray scattering or momentum resolved electron energy loss spectroscopy as discussed further in section 5.5.
In 2D materials, the electron-hole (e-h) exchange interaction entering the Bethe-Salpeter equation becomes proportional to q q ∥ for small momentum transfer, see discussion around equations (43)- (45). This has been shown to result in an anomalous linear (light-like) dispersion for excitons in some of the 2D semiconductors including MoS 2 , phosphorene and hBN [98][99][100]. Only the singlet exciton exhibits linear dispersion because the exchange term vanishes for the triplet exciton where the electron and hole have opposite spin. Interestingly, the standard parabolic exciton dispersion that is inherited directly from the single-particle band energies, was found for the lowest singlet exciton of graphane. This indicates a suppression of the e-h exchange arising from a small spatial overlap of the electron and hole wave functions in this particular material [43]. It has been proposed that the difference in exciton dispersion, i.e. whether linear or parabolic, could be used to distinguish between Wannier-and Frenkel-like excitons in 2D semiconductors where the usual distinction based on the binding Figure 12. Rydberg series of a monolayer of WS 2 calculated with the generalized hydrogen model with linear screening (equations (30) and (32)) and from the solution of the effective 2D Schrödinger equation (23). The results are compared with the 2D Hydrogen model with a fixed state-independent dielectric constant. The inset shows the effective n-dependent dielectric constants obtained from equation (32). Figure reproduced from [90]. energy strength becomes difficult to maintain due to the universally large binding energies of 2D excitons [99].

Trions and biexcitons
Beyond excitons, the rich excitation spectrum of 2D materials features electron-hole complexes involving more than two particles. In doped semiconductors, free electrons (in n-doped regions) or holes (in pdoped regions) can interact with the neutral excitons forming three-particle complexes known as trions. Experimental binding energies of trions (relative to the free particle and the exciton) in monolayer TMDs are in the range 20-30 meV [24,44,101] while binding energies as high as 150 meV have been reported for phosphorene [102]. The size of the trion binding energy is thus larger than or comparable to the thermal energy at room temperature suggesting that trions should play active roles in the photophysics of monolayer semiconductors. For comparison, trion binding energies in conventional semiconductor quantum wells are in the range 1-5 meV [103].
The exciton model Hamiltonian equation (23) can be straightforwardly generalised to an N-particle 2D excitonic system where ∥ ( ) r r i and * m i are the position and effective mass of particle i. By adopting relative coordinates the problem reduces to an (N − 1)-particle problem. For trions (N = 3), this type of Hamiltonian, with the analytical 2D potential of equation (24) and parameters (effective masses and screening lengths) obtained from DFT, has been solved for the MoX 2 and WX 2 monolayer TMDs by variational optimisation of simple trial wave functions [44], path-integral Monte Carlo [104], and explicitly correlated Gaussians [105]. In all cases, reported binding energies for the lowest trion lie within a factor of two of the experimental values. Using similar methods, biexciton (N = 4) binding energies of around 20 meV were predicted, in reasonable agreement with experiments (50-70 meV) [106,107].
Recently, a first-principles approach to trion excitations based on an extension of the BSE to three-particle states was developed and applied to carbon nanotubes [108]. The significantly larger size of the configurational space of trions as compared to excitons prohibits standard diagonalization of the three-body Hamiltonian. Fortunately, all Coulomb interaction terms are of two-particle nature and the resulting Hamiltonian is very sparse, allowing for iterative solutions of the corresponding eigenvalue problem [108].

Excitons in van der Waals heterostructures
It was shown in the previous section that the simple 2D model for dielectric screening, ( ) πα = + ε q q q 1 2 2D , gives an accurate description of excitons in isolated atomically thin materials. The reason is that typical exciton radii in monolayer semiconductors and insulators are larger than the thickness of the material such that the screening is mainly governed by small q q-vectors where ( ) ε q q 2D is accurate, see figure 2. This condition is bound to break down for thicker structures or supported/encapsulated 2D materials. In such cases the full non-linear q q-dependence of the dielectric function must be taken into account. Figure 13 shows the static dielectric function of monolayer MoS 2 deposited onto (upper panel) and embedded inside (lower panel) an N-layer stack of hBN. The dielectric function of the MoS 2 monolayer has been calculated using the QEH model, see [45] for details. It is clear that the region of ∥ q q -vectors where the strict 2D screening model (dashed lines) is valid shrinks as the number of hBN layers increases. The exciton binding energies and radii obtained from the effective Hamiltonian equation (23) with a screened interaction evaluated using (i) the exact ∥ q q -dependent di electric function and (ii) the linear 2D screening model, are shown in figure 14. The 2D model underestimates the binding energy and overestimates the radius significantly for more than just a few hBN layers and the results diverge in the limit of infinite hBN support. The reason is that the size of the exciton calculated using method (i), and indicated by the vertical grey lines, does not change much with the number of hBN layers and consequently the range of ∥ q q -vectors that are important for the screening also remains fairly constant.
As another example, we study the 2D to 3D transition of the exciton in MoS 2 . In layered bulk materials, the exciton equation can be written as follow: In layered materials like MoS 2 , the exciton mass in the out of plane direction is typically much higher than in the in-plane directions ( ∥ µ µ ⊥ ). Consequently, we can neglect the out of plane component of the kinetic energy and be left with the 2D exciton model. Moreover, the in-plane effective mass does not differ considerably between monolayer and bulk MoS 2 [109]. Therefore, one would expect the main difference between the physics of excitons in monolayer and layered bulk to be governed by the screened potential rather than the geometric confinement. To test this hypothesis we consider a multilayer MoS 2 structure and calculate the binding energy of an exciton localized in the central MoS 2 layer using the effective Hamiltonian equation (23) with the screened potential calculated from the QEH model. The results for the binding energy as a function of the number of MoS 2 layers are plotted in figure 15. As expected, the reduction of the exciton binding energy is larger when the monolayer is embedded in MoS 2 compared to hBN (see figure 14 panel(b)). Remarkably, the binding energy converges towards a value of 0.16 eV only 0.03 eV higher than previously reported first principles BSE results for bulk MoS 2 [82]. This shows that the large difference in exciton binding energies of 2D and 3D (layered) bulk semiconductors comes mainly from the difference in screening while quantum confinement plays a minor role.

Interlayer excitons
In addition to the intralayer excitons discussed in previous sections, vdWHs can host more complex types of excitons where the electron and hole reside in different layers, so-called (spatially) indirect excitons or interlayer excitons. Because of the spatial separation of the electron and hole, interlayer excitons posses longer lifetimes [112,113] than intralayer excitons, which make them ideal candidates for realization of exotic many-particle states such as room temperature Bose-Einstein condensates [114]. Indirect excitons also play a central role in the important charge separation process of excitonic solar cells [115,116]. In this case, the binding energy of the indirect exciton determines the minimum band offset required to separate the electron and hole and thus the achievable output voltage.
A realistic model of interlayer excitons in vdWHs must (i) start from an accurate description of the band alignment at the heterointerface and (ii) account for the non-local screening of the electron-hole interaction by the 2D host materials. Below the importance of these effects are illustrated using the heterostructure MoS 2 / WSe 2 as an example. As demonstrated experimentally in [110], by intercalating hBN spacer layers between the optically active TMD layers it is possible to vary the electron-hole separation with atomic precision making this system ideally suited for exploring the properties of interlayer excitons.
It was demonstrated in section 3, that QP band structures of isolated monolayers can be accurately calculated by the GW method. Moreover, it was shown how the effect of interlayer screening on the band structure of vdWHs can be included within a GW framework by utilizing the QEH model to evaluate the change in screened interaction (∆W ) due to the surrounding layers. Figure 16(a) shows the G 0 W 0 @LDA calculated band alignment for the MoS 2 /WSe 2 bilayer. The coloured bars represent the vacuum level aligned monolayer band structures (same as figure 6) while the black lines show the result when interlayer screening is included via the G∆W method. The effect is rather small due to the modest screening provided by a single neighbouring layer and is comparable to the band gap reduction of ∼0.2 eV found for MoS 2 on graphene in figure 8. It should be noted that DFT calculations performed for realistic MoS 2 /WSe 2 heterostructure configurations employing large supercells to account for the  lattice mismatch, show that the hybridisation between the two layers is negligible, in particular around the relevant band extrema at the K-point [111].
An interlayer exciton can be modelled as a bound electron-hole pair with the particle motion confined to two normally displaced 2D planes. This implies that the binding energy of the interlayer exciton can be computed from the usual 2D exciton Hamiltonian equation (23) provided that the exciton mass is calculated from the effective electron and hole masses of the relevant layers and that the screened interaction takes the (constant) out-of-plane separation of the electron and hole into account. Figure 16(a) shows the binding energies of the intra-and interlayer excitons of the MoS 2 /WSe 2 bilayer calculated in this way. The effective masses were obtained from the monolayer band structures and the screened interaction was calculated using the QEH model, see [111]. The binding energy of the interlayer exciton (0.28 eV) is only around 0.1 eV smaller than the binding energy of the intralayer excitons (0.39-0.42 eV). Upon insertion of hBN spacer layers, the binding energy of the interlayer exciton decreases as ∼1/d as the distance between the electron and hole (d) increases, while the intralayer excitons are only weakly affected [111].
In [110] the photoluminescence (PL) due to exciton recombination in MoS 2 /WSe 2 -based heterostructures was measured. The measured PL spectra are reproduced in figure 16(b) for the individual MoS 2 and WSe 2 monolayers (left panel) and the hBN intercalated heterostructures (right panel). The position of the lowest optical transition calculated by subtracting the exciton binding energy from the QP band gap calculated using different methods is shown by the symbols. When the G 0 W 0 band gap is used, the calculated PL peaks are in excellent agreement with the experiments for both the intralayer and interlayer excitons. The small underestimation (∼0.1 eV) of the interlayer exciton peak is also found for the WSe 2 intralayer exciton indicating that the origin is a too high positioning of the WSe 2 valence band by the G 0 W 0 method.

Field-induced exciton dissociation
As illustrated by the previous sections, excitons in atomically thin 2D semiconductors can be extremely strongly bound. While this has some exciting fundamental perspectives [114], it is problematic for applications such as photodetectors and solar cells which rely on efficient conversion of photons into electrical currents. Indeed, photocurrent measurements on suspended MoS 2 show that large voltages are required to dissociate the excitons into free carriers [85]. On the other hand, the previous section showed that van der Waals heterostructuring could be employed as a strategy to weaken the electron-hole interaction and thereby lower the exciton binding energy without altering the band structure of the 2D semiconductor (apart from a reduction of the band gap due to image charge screening as described in section 3.4). As shown in the following, encapsulation of a single MoS 2 layer by just two layers of hBN is sufficient to increase exciton dissociation rates by an order of magnitude [117].
When a constant electric field is applied to an exciton, it will eventually decay into a free electron and hole [118,119]. The application of a constant electric field changes the exciton from a bound state to a resonance with a finite spectral line width, , where τ is the resonance life time or equivalently the inverse exciton dissociation rate, see figure 17. There are generally two approaches used to compute resonances. The so-called indirect methods identify resonances as the poles of the scattering amplitude analytically extended into the complex energy plane [120], while the direct methods obtain the resonance states directly as the eigenstates of a complex scaled non-hermitian Hamiltonian [121]. Figure 18 shows the dissociation rates of the lowest exciton in MoS 2 as a function of in-plane field strength with and without hBN encapsulation. The rates have been obtained by complex scaling of the exciton Hamiltonian equation (23) using the QEH model to include the screening from hBN. As expected, larger fields lead to shorter lifetimes and the dependence of the rate on the field strength is approximately Figure 17. Illustration of a 2D exciton in the absence (left) and presence (right) of a constant in-plane electric field. The exciton potential is shown in blue, the exciton wave function is sketched in green, and the energy is shown in red. When an electric field is applied, the energy of the exciton shifts down and the sharp energy peak is broadened due to the coupling to the continuum of states. Figure reproduced from [117]. considered field strengths. It can also be seen that the dissociation rate can be tuned to a high degree by encapsulating the MoS 2 layer in hBN. This is as expected, since the dielectric screening from the hBN weakens the e-h interaction which makes the exciton dissociate more readily.
In a real device, the field-induced dissociation of excitons competes with other decay mechanisms such as direct radiative recombination [122], defect-assisted non-radiative recombination [66], and exciton-exciton annihilation [123]. The relative importance of these effects is highly dependent on temperature, defect concentration, and exciton density.
At very low temperatures, the direct radiative decay of zero momentum excitons dominates, with a characteristic lifetime of ∼200 fs. [122,124] At room temperature, most of the excitons have non-vanishing momenta, and the radiative recombination lifetime is on average ∼1 ns [66,122]. Under such conditions, defect-assisted recombination becomes the dominant mechanism reducing exciton lifetimes to 2-5 ps [66,125]. Exciton-exciton annihilations become important only when the density of excitons in a sample is large. At a density of × 1 10 12 cm −2 , the effective lifetime due to annihilation is around 10 ps [123].
The calculations show that for field strengths larger than 100 V μm −1 , the dissociation lifetime is shorter than one picosecond for all the systems considered. This is shorter than the smallest characteristic lifetimes of the alternative decay channels at room temperature (indicated by the grey shaded region in figure 18) and thus field-induced dissociation should dominate. We note that a potential gradient of 100 V μm −1 is not unlikely to exist in the metal-MoS 2 contact region where charge transfer and interface dipole formation driven by Fermi level mismatch can lead to significant variations in the potential and band energies even in the absence of an applied bias voltage. Recently, a chemical treatment has successfully been used to eliminate the influence of defects on exciton decay, leaving radiative decay as the main decay channel and leading to exciton lifetimes of 10 ns [40]. Under such conditions, field dissociation should become the dominant decay mechanism for field strengths as small as 10 V μm −1 . Recently, the complex scaling method was also used to study field-induced exciton dissociation in layered bulk TMDs [126] and electroabsorption in monolayer TMDs [127].

Plasmons
Plasmons in metals are quantised versions of the collective oscillatory modes of the electron gas against the background of positively charged ions [128]. By their collective nature, plasmons couple strongly to electromagnetic fields. This makes it possible to utilize metallic structures to confine, enhance, or guide light fields and thereby enable a range of applications from single-molecule spectroscopy [129,130] to energy conversion [131,132] and imaging [133,134]. In this context, 2D materials present a number of unique properties which could lead to improved performances or new functionalities of plasmonic devices. For example, graphene plasmons are long lived, can be easily tuned via electrostatic or chemical doping, and can be used to confine light to extremely small volumes [135,136]. In addition, van der Waals stacking of different 2D metals and dielectrics opens up the possibility of controlling plasmons on the atomic length scale and for creating deep-subwavelength metamaterials with unique optical properties that cannot be found in nature or achieved with conventional heterostructures.

The Drude model and its limitations
The semi-classical Drude model, which describes the electron system as a gas of classical particles subject to a constant empirical friction force, is the workhorse for describing electron dynamics in metal structures. Within the Drude model, the dielectric function of a metal takes the form Figure 18. Calculated field-induced dissociation rate of the lowest exciton in a single MoS 2 layer in isolation, deposited on a single hBN layer or sandwiched between two hBN layers. The range of decay rates due to other recombination mechanisms at room temperature are indicated by the grey shaded area. They span from non-radiative defect assisted recombination (upper limit) to the direct radiative recombination (lower limit). Figure reproduced from [117].
where / ω π = * n m 4 p is the bulk plasma frequency (n the electron density and m * the effective carrier mass) and τ is the empirical relaxation time. For simple bulk metals and surfaces, equation (35) provides an excellent description of the optical properties for frequencies below the onset of interband transitions where v F is the Fermi velocity. The latter condition expresses that the wave length of the driving field is either longer than the distance travelled by an electron during one period of the field or the distance travelled by an electron between two collisions, i.e. the mean free path, τ v F . This condition defines the regime where the q q-dependence of the dielectric function, also known as spatial dispersion or non-local response, can be neglected. For finite nano structures, the optical response naturally involves large q q vectors and non-local effects become important [137][138][139]. Furthermore, surface scattering and spill-out effects must be taken into account in this regime [140]. For noble-and transition metals, the semi-classical picture must be modified to account for the less mobile d-electrons as well as the additional damping due to s − d interband transitions [141]. For semiconductors, a full quantum theory is required to describe interband plasmons which do not have a direct classical analogue.
Quantum mechanically, the plasmon is a collective excited state which can be represented as a coherent superposition of electron-hole transitions (intraband or interband or a mix thereof). As such the plasmon state, in particular for interband plasmons, is very similar to the exciton state. However, while the energy of the exciton is lower than the energy of its constituent electron-hole pairs, the plasmon energy is higher. Thus one way to differentiate between the two types of excitations is to say that the collective nature of the states, i.e. the coherent rather than incoherent superposition of electron-hole pairs, leads to an increased (decreased) Coulomb energy for plasmons (excitons). As will be shown in section 5.6 the response of the two types of collective excitations to external probes is very different in 2D materials.

Plasmons in 2D metals
For an isotropic metallic film of charge density n and thickness d, solution of Maxwell's equations in the non-retarded regime (speed of light → ∞ c ) yields two plasmons with energy dispersion [142] ( ) ω ω = ± − e q q 2 1 p qd (36) The two solutions correspond to tangential and normal oscillations of the electron gas. For qd 1, the energy of both plasmon modes converge towards the surface plasmon energy of a semi-infinite metal, / ω 2 p . This happens because the penetration depth of the field associated with the plasmon becomes smaller for larger q q, and consequently the density oscillations on the two surface become decoupled in this regime. For an atomically thin conductor, we can assume qd 1, and the tangential solution becomes (37) In contrast to 3D metals where the plasmon energy remains constant for → q q 0 0, the 2D plasmon vanishes in this limit.
To understand the different behaviour of plasmons in 2D and 3D, it is instructive to consider the Dyson equation that links the long range ( = G G 0 0) components of the irreducible and reducible response functions via the Coulomb interaction, The irreducible response function includes local field effects and e-h exchange-correlation effects. It satisfies the Dyson equation where matrix multiplication in the G G and ′ G G indices is implied. χ 0 is the non-interacting response function from equation (8), is the short range part of the Coulomb interaction, and ( ) ω f q q, xc is the xc-kernel. For metals, the simplest adiabatic local density approximation (ALDA) to f xc is often sufficient for describing the plasmonic properties and generally improves the RPA results.
Physically, the difference between χ 0 00 0 and χ 0 00 0 irr is that the former includes long-range screening, i.e. the effect of potentials arising from variations in the induced density on length scales larger than the unit cell. In 3D metals, these are the potentials leading to plasmon formation, but in 2D, they are too weak to sustain the plasmon in the small q q limit. In general, i.e. independent of system dimensionality, ( ) χ ω ∝ q q q, 0 00 0 irr 2 in the small q limit (for metals this holds for ω v q F ). For a 2D material, since ∝ V q q q 1 2D ( ) / for → q 0, it follows from equation (38) that This means that also for metals, for → q q 0 and finite ω. The fact that there is no difference between the total and external macroscopic fields means that there is no distinction between the loss spectrum and optical absorption spectrum of a 2D material.
It should be stressed that, in contrast to the metallic intraband plasmons, the formation of interband plasmons can be driven by short range Coulomb interactions and thus may be present, although weak, in the → q q 0 0 limit. An example illustrating the different behavior of intraband and interband plasmons in 2D is provided by doped graphene. The doping introduces free carriers in the π or π * bands which leads to the formation of a metallic intraband plasmon with the q-dispersion characteristic of a 2D metal in the long wavelength limit [26,143]. In contrast, the π π − * interband plasmon approaches a finite energy for → q q 0 0 [144]. Interband plasmons are discussed further in section 5.6.

Plasmons in metallic transition metal dichalcogenides
A highly attractive feature of the metallic plasmon in doped graphene is that its frequency can be tuned by controlling the density of free carriers e.g. via electrostatic gating. On the other hand, the energy of the metallic plasmon is restricted by the achievable carrier concentration and this limits applications of graphene plasmonics to the terahertz and infrared regimes. Metallic 2D materials have much higher charge carrier densities leading to plasmon energies up to 1 eV or above and could potentially extend 2D plasmonics to the near infrared or even optical regimes.
In [27] the (intraband) plasmons of six different metallic monolayer TMDs were explored by first-principles RPA calculations. As an example, figure 19 shows the plasmon dispersion of TaSe 2 . The grey shaded areas indicate the dissipative regions of Landau damping formed by inter-and intraband transitions, respectively. The important role of the interband transitions for the plasmon energy and strength (indicated by the dot size) is evident from the large difference between the results of a full calculation based on the true response function which includes all possible transitions (blue dots), and a restricted calculation based on a response function which includes only intraband transitions (red dots). Clearly, the pure intraband plasmon is well described by the 2D free electron model equation (37) while the true plasmon is significantly redshifted by the coupling to interband transitions.
In addition to the extended sheet plasmons, certain 2D materials have been predicted to host one-dimensional (1D) plasmons. For example, first-principles calcul ations revealed the existence of such 1D plasmons along the metallic edges of MoS 2 nanoribbons [145]. Furthermore, as the metallic edge states have been proposed to be a general feature of polar 2D materials like the TMDs [146], 1D plasmons may be found at edges or domain boundaries of many other 2D materials.

Plasmons in graphene/hBN heterostructures
The prospects of graphene plasmonics were long obstructed by the strong damping of the its plasmons. Recently, it was demonstrated that by encapsulating graphene by two hBN films the losses can be significantly reduced, presumably down to the intrinsic limit set by Landau damping, electron-phonon and electronelectron interactions within the graphene itself [136]. Below we explore some properties of plasmons in graphene/hBN/graphene heterostructures with a focus on the possibility of tuning the plasmon energy via plasmon-plasmon hybridisation and dielectric screening.
The close to perfect lattice match between graphene and hBN allows full first-principles calculations to be performed for the thinnest heterostructures, see figure 20. Here we use doped graphene that has a finite density of states at the Fermi level, giving rise to metallic plasmons with energies in the range 0-2 eV depending in the in-plane momentum transfer. Note that such high plasmon energies require extremely high doping levels that are probably not achievable in practice. The plasmon energies go to zero in the optical limit as is characteristic for plasmons in 2D metals, see equation (36). We calculate the effect of hBN on the plasmons using the QEH model for up to 20 layers of hBN and compare to full first-principles calculations for 1-3 layers of hBN.
To identify the plasmons of the heterostructure we follow [147]. In brief, we compute the eigenvalues, ( ) ω ε n , of the heterostructure dielectric function for each frequency point and identify a plasmon energy, ħω p , from the condition The corre sponding eigenvector, ( ) φ ω n p , represents the potential associated with the plasmon oscillation (not shown). This analysis identifies two plasmons corre sponding to the symmetric (++) and antisymmetric (+−) combinations of the graphene plasmons as previously found for two freestanding graphene sheets [148]. More details on these calculations can be found in [37].
The middle panel of figure 20 shows the calculated dispersion of the symmetric (red) and antisymmetric (blue) plasmons for 1 and 3 layers of hBN, respectively. For ∥ q d 1, the charge density oscillations in the two graphene sheets do not interact and both plasmon modes converge to that of a single freestanding graphene sheet. The right panel shows the evolution of the plasmon energies at a fixed in-plane momentum transfer as function of the distance between the two graphene sheets. The full (dashed) lines represent the case of vacuum (hBN) filling the gap. As expected, the hBN dielectric redshifts the plasmon energies by screening the interaction between the two charge densities. For 1-3 hBN layers, the QEH model perfectly reproduces the first-principles results for the plasmon energy.

Electron energy loss spectroscopy
Low loss electron energy loss spectroscopy (LL-EELS) is a commonly used technique to probe plasmons and interband transitions. In angular resolved LL-EELS the momentum loss is simultaneously recorded which allows the dispersion of the plasmon to be determined. As a complementary tool, spatially resolved LL-EELS can provide information about the energy loss at a specific point in space. Today, state of the art scanning transmission electron microscopes (STEM) provide loss spectra with simultaneous spatial and energy resolutions below 1 Å and 10 meV, respectively [149].
From Fermi's Golden rule and the spectral representation of the reducible response function χ, the power absorbed by the material when exposed to an external potential of the form ( ) where a mixed ( ) ∥ z q q , -representation of χ has been used. Discussions of the interpretation and calculation of electron energy loss spectroscopy for 2D materials can be found in [150,151].

Interband plasmons versus excitons
The similar structure of interband plasmons and excitons has previously been alluded to. Below the difference between these two types of collective excitations is explored in more detail. It turns out that 2D materials present some novelty in this regard compared to the well   (21) it follows that the collective nature of the excitations, i.e. the mixing of the e-h pairs, is governed by the xc-kernel consisting of the e-h exchange (V) and direct screened interaction (W). These terms respectively increase and decrease the energy of the collective excitation relative to that of the constituent e-h pairs. Roughly speaking, this means that the e-h exchange interaction is responsible for plasmon formation while the direct e-h interaction is responsible for exciton formation. Consider first the relevant e-h exchange matrix element in the BSE Hamiltonian (where k k and q q are vectors in the 2D Brillouin zone): This equation shows that the long range ( = G G 0 0) component of each of the wave functions ψ S and ψ ′ S are proportional to q in the long wave length limit. In the same limit the (quasi) 2D Coulomb interaction takes the form 1/q. Therefore the long range part of the Coulomb interaction does not contribute to ( ) ′ V q q SS for small q q. On the other hand, short range Coulomb interactions generally have little effect on the excitations due to the rapid decay of the Coulomb interaction with wave vector, see equation (5). We thus conclude that for a 2D material the e-h exchange, which drives plasmon formation, is very weak in the small q q limit. Note that this is very different from the 3D case where the long range part of the Coulomb interaction dominates the e-h exchange for small q q.
In contrast to the e-h exchange operator, the direct e-h interaction is relatively independent of q q. In particular, the contribution to ( ) ′ W q q SS from the long range component of the Coulomb interaction does not vanish in the → q q 0 0 limit but in fact dominates the matrix element. Thus we must expect that the collective excitations of a 2D material are of excitonic nature for small q q while plasmons become increasingly important for larger q q. The example given below shows that this is indeed the case. Figure 21(a) shows the loss spectrum of 1-3 layer MoS 2 for four different values of the in-plane momentum transfer q q. All spectra are normalized by the number of layers. The response function has been calculated at the RPA level using DFT-PBE wave functions for the multilayer films (full lines) or using the QEH model to include interlayer interactions (dashed lines). The excellent agreement between the full calculations and the QEH model demonstrates that interlayer hybridisation has negligible effect on the loss spectra.
A consistent description of plasmons and excitons should ideally be based on the BSE rather than the RPA. Indeed, the RPA neglects the direct e-h interaction (W) and thus cannot account for excitons. However, the computational cost of the BSE makes it impossible to describe the excitation spectrum accurately up to energies relevant for the plasmons (around 5-20 eV for MoS 2 ). On the other hand we argue that the low energy part of the RPA spectrum behaves qualitatively similar to the BSE spectrum with respect to the q q-dependence.
Returning to figure 21(a) it can be seen that for = q q 0 0 the loss spectrum is independent of film thickness. This is a result of the absence of long range screening which leads to an effective decoupling of the layers. In fact, the RPA spectrum consists essentially of uncoupled single-particle transitions because the e-h exchange coupling vanishes (almost) for = q q 0 0. The peak at ∼ 3 eV originates from transitions at the Γ-point of the BZ. These are the same transitions that form the C exciton giving rise to the peak at 2.5 eV in the BSE spectrum (see figure 9) [83]. Thus the 3 eV peak in the RPA spectrum represents an uncorrelated version of the C exciton. For finite q q, plasmon peaks around 8 eV and 14 eV emerge and completely dominate the loss spectrum. The plasmon peaks are sensitive to both q q and film thickness which again underlines the strongly q q-dependent screening in 2D (an effect also referred to as spatial dispersion or non-local response). We conclude that excitonic states dominate the loss spectra for small q q while plasmonic states dominate for larger q q. Figure 21(b) shows the local loss spectrum of mono layer MoS 2 (blue) obtained from equation (42). For comparison the = q q 0 0 spectrum is shown. The local spectrum involves a summation over all q q-vectors and consequently exhibits peaks due to both excitons and plasmons. Physically, the localised electron beam probes the collective excitations at all momenta; it thus represents a sensitive but not very selective probe. Figure 21(c) shows experimentally measured local loss spectra obtained at different positions of an MoS 2 sample of varying thickness. The spectra increase in intensity as the beam is moved from 50 nm outside the sample over regions with increasing number of MoS 2 layers. Figure 21(d) shows the calculated local, i.e. q q-summed, loss spectra for 1-10 layer MoS 2 calculated using the QEH model. The agreement with experiments is reasonable although the RPA calcul ations and experiments show the opposite trend regarding the relative intensity of the C-exciton and the 8 eV plasmon peaks. Interestingly, the experimental spectra recorded up to 50 nm away from the MoS 2 sample are dominated by the 3 eV excitonic peak rather than the plasmons. A possible explanation is that when the beam is far away from the sample only the long wave length (small q q) components of the probing beam potential reaches the MoS 2 . Therefore, by varying the position of the electron beam relative to the sample it is possible, to some extent, to selectively excite collective excitations with different q q.

Conclusions
This topical review has discussed some of the concepts, theories, and computational methods used to describe elementary electronic excitations in atomically thin 2D materials and layered van der Waals heterostructures. It was shown that the weak and highly non-local nature of the dielectric screening in atomically thin crystals is responsible for a number of unique properties exhibited by the 2D materials. The non-conventional physics of excitons in 2D materials was illustrated and it was shown how the non-local screening can be incorporated into a 2D Hydrogenic model providing analytical expressions for the exciton properties with quantitative accuracy and simple explanations for the non-Hydrogenic Rydberg series, non-degeneracy of excitonic states with different l quantum numbers, and the (almost) linear dependence of the exciton binding energy with band gap. The problem of exciton dissociation in 2D materials by an in-plane electric field was discussed in the context of the complex scaling method. It was shown that the dissociation rate can be greatly enhanced by encapsulating the 2D material in a dielectric which screens the e-h interaction. For plasmons, the absence of long range screening in 2D leads to acoustic dispersion relations for the metallic intraband plasmons and suppresses the formation of plasmonic states in general for small momentum transfers, q q. Consequently, it was shown that the loss spectrum of a 2D material is dominated by excitonic states for small q q and plasmonic states for larger q q. The weak intrinsic screening in the 2D materials make them highly sensitive to the environment such as substrate or encapsulation. It was demonstrated that dielectric screening from the environment can shift the collective and single-particle excitation energies by up to 1 eV. The fact that standard density functional theory (DFT) methods are unable to account for these screening effects highlights the need to go beyond the independent particle approximation for band structure calculations. The application of the many-body GW method to 2D materials was briefly reviewed, and a method to overcome the slow k k-point convergence of such calculations was described. The quantum electrostatic heterostructure (QEH) model was presented as a method for computing the dielectric properties of general, incommensurate van der Waals heterostructures with first-principles accuracy and minimal computational cost. The ability of the QEH model to provide an accurate and seamless connection between the screening properties of 2D and layered 3D materials was used to illustrate the transition of the exciton in MoS 2 from the monolayer to the bulk. These calculations demonstrated that the large exciton binding energy in 2D materials is a result of weak dielectric screening rather than quantum confinement.