k.p theory for two-dimensional transition metal dichalcogenide semiconductors

We present $\mathbf{k}\cdotp\mathbf{p}$ Hamiltonians parametrised by {\it ab initio} density functional theory calculations to describe the dispersion of the valence and conduction bands at their extrema (the $K$, $Q$, $\Gamma$, and $M$ points of the hexagonal Brillouin zone) in atomic crystals of semiconducting monolayer transition metal dichalcogenides. We review the parametrisation of the essential parts of the $\mathbf{k}\cdotp\mathbf{p}$ Hamiltonians for MoS$_2$, MoSe$_2$, WS$_2$, and WSe$_2$, including the spin-splitting and spin-polarisation of the bands, and we discuss the vibrational properties of these materials. We then use $\mathbf{k}\cdotp\mathbf{p}$ theory to analyse optical transitions in two-dimensional transition metal dichalcogenides over a broad spectral range that covers the Van Hove singularities in the band structure (the $M$ points). We also discuss the visualisation of scanning tunnelling microscopy maps.

Like graphene, the group-IVB monolayer TMDCs MX 2 (where M=Mo or W and X=S or Se) considered in this work have hexagonal lattice structures, and the extrema (valleys) in the dispersion relations of both the valence and conduction bands (VB and CB) can be found at the K and −K points of the hexagonal Brillouin zone (BZ). Unlike graphene, however, these 2D crystals do not have inversion symmetry. The minimalistic approach to the theoretical modelling of monolayer TMDCs is therefore based on mimicking them as graphene with a staggered sublattice potential that breaks inversion symmetry [10,25,26]. This approach captures some optical and transport effects related to the valley degree of freedom of the electrons [26,27,28,29,30,31]. The staggered graphene analogue [26] has also been generalised to the tight-binding (TB) description of TMDCs [26,49,53,54,55,56,57,39], but this approach suffers from the large number of atomic orbitals that have to be included on each site and the need for beyond-nearest-neighbour hopping to account for the variation of the weight of individual atomic orbitals in the band wave functions across the BZ revealed by detailed density functional theory (DFT) modelling (see, e.g., Figure 3). As a result, the accumulation of experimental data [32,33,34,35] and the drive towards the implementation of monolayer TMDCs in practical devices call for more detailed theoretical descriptions of these materials and more compact models of their electronic properties.
In this Review, we describe two complementary theoretical approaches that have recently been used to achieve a detailed description of the electronic properties of these materials. One consists of ab initio DFT modelling of the band structure, which has the potential to be accurate. DFT can be combined with transport codes [17,36,37,38,39,40,41,42] or used to calculate optical spectra [43,44,45,46], but ab initio calculations are prohibitively expensive for many practical problems focused on modelling devices and studies of quantum dots [47]. Moreover, magnetic-field effects [47,48,49,50,51,52] and certain questions regarding neutral and charged excitons [155] cannot easily be addressed by DFT-based techniques. The second approach uses the k· p methodology [59,60,156], which exploits the symmetries of the system. This approach provides an accurate characterisation of the dispersion of the valence and conduction bands in the vicinity of the K and −K points of the BZ in terms of a relatively small number of parameters [61]. Magnetic-field and spin-orbit coupling effects can also be taken into account in a straightforward way [47]. In contrast to DFT modelling, this method is only valid in the vicinity of certain high-symmetry k-space points; however, for those intervals, it enables one to quantify all the essential features of the electronic properties. One can also relate a k· p Hamiltonian to a particular TB model [26,49,54], although it is not necessary to set up a TB model in order to derive a k· p Hamiltonian. Here we present phenomenological k· p Hamiltonians derived for all extrema of the bands (at the K, Q, Γ, and M points of the BZ) using the symmetry properties of TMDC atomic crystals, with specific material parameters obtained by fitting them to the DFT band structures of MoS 2 , MoSe 2 , WS 2 and WSe 2 . The DFT calculations discussed in this Review were performed using the vasp [62] and fleur [63] codes. The robustness of our results is well illustrated by the close agreement between the results obtained from these two different first-principles codes and through comparison to all available experimental results.
This Review is organised as follows. Section 2 is devoted to the crystalline lattice parameters and vibrational properties of TMDCs. Sections 3 and 4 discuss spinsplitting due to spin-orbit coupling (SOC) and band width [relevant for angle-resolved photoemission spectroscopy (ARPES) studies of TMDCs]. Sections 5, 6, 7, and 8 describe the structure and parametrisation of k· p Hamiltonians for K, Q, Γ, and M points of the BZ, respectively. Finally, we draw our conclusions in Section 9.

Lattice parameters, band-structure calculations and vibrational properties
The crystal structure of each MX 2 monolayer considered in this work consists of three atomic layers, X-M-X. Within each layer the M or X atoms form a 2D hexagonal lattice: see Figure 1. The M atoms in the middle plane are surrounded by three nearestneighbour X atoms in both the bottom and the top layer so that the crystal has D 3h symmetry. The crystal structure is characterised by the in-plane lattice constant a 0 and the distance d X−X between the two chalcogen planes. It has already been noted [65] that certain details of the band structure obtained from DFT calculations depend rather sensitively on a 0 and d X−X . Indeed, we have also found that agreement with the available experimental results regarding, e.g., the effective mass m vb Γ at the Γ point of the BZ or the energy difference E KΓ between the top of the VB at the K and Γ points can only be achieved if the values of a 0 and d X−X fall in a rather narrow range. As a first step, we have used two approaches to calculate the basic lattice parameters a 0 and d X−X . The first approach used vasp [62]. The vasp geometries were calculated using the Heyd-Scuseria-Ernzerhof 2006 (HSE06) exact-exchange density functional [71]. The plane-wave cutoff energy was set to 600 eV and the BZ was sampled by a 12 × 12 × 1 Monkhorst-Pack grid. The vertical separation between the layers was set to 20Å to make the interaction between the repeated images of the layer in the three-dimensional cell negligible. Optimisation was carried out until atomic forces fell below 0.005 eV/Å. The second approach used the full-potential linearised augmented plane-wave (FLAPW) method as implemented in the fleur code [63]. The FLAPW method is an all-electron method within DFT. The fleur code allows 2D systems to be studied without constructing slabs in three-dimensionally periodic cells and the resulting electronic spectra are free of plane-wave continua. All our fleur calculations were carried out with a cut-off k max of 10.6 eV −1 for the plane-wave basis set and 144 k points corresponding to a 12 × 12 × 1 Monkhorst-Pack grid in the irreducible wedge of the BZ. Muffin-tin radii of 1, 1.21, 1.27, and 1.27Åwere used for S, Se, Mo, and W, respectively. We note that considering local orbitals for Mo (s, p), Se (s, p, d), and W (s, p, f ) to improve the linearised augmented plane-wave basis proved to be crucial for a correct description of the excited states. We used the Perdew-Burke-Ernzerhof (PBE) generalised gradient approximation [72] to the exchange-correlation potential. The structures were relaxed (with the effects of SOC included) until the forces were less than 0.0005 eV/Å.
The calculated values of a 0 and d S−S for monolayer TMDCs are shown in Table  1 and compared to measured values for the corresponding bulk materials. The lattice parameters obtained from the first of the DFT approaches described above are shown in the rows labelled by "(HSE)", the ones from the second approach are in the rows labelled by "(PBE)". "(Exp)" indicates experimental results found in the literature. Although Table 1. Lattice vector a 0 and chalcogen-chalcogen distance d X−X as obtained from DFT calculations. Experimental values for the corresponding bulk material are shown in rows labelled by "Exp".  [73]. b [76]. c [74]. d [77]. e [78]. f [79], measurement at 293 K. g [75]. h [80].
there is some scatter in the experimental data, Table 1 suggests that using the HSE06 functional to relax the monolayer crystal structure leads to a good agreement with the room-temperature empirical bulk a 0 values. On the other hand, the PBE functional seems to slightly overestimates a 0 . However, the situation is less clear in the case of d X−X . We note that both the HSE06 and the PBE results are in good agreement with Reference [81]. Recent experiments show that the energy of the photoluminescence peak is quite sensitive to the temperature [3,82,83], which can be understood in terms of the dependence of the band structure on a 0 and d X−X . Indeed, a recent computational study [84] was able to qualitatively reproduce the redshift of the photoluminescence peak of MoS 2 as a function of temperature by assuming a thermal expansion of the lattice. The good agreement between the calculated lattice parameters and the corresponding experimental ones suggests that, interestingly, the predictions based on our DFT results are expected to be most accurate at room temperature (except for the band gap, which is known to be underestimated by DFT). To our knowledge systematic measurements of the temperature-dependence of the lattice parameters of bulk MX 2 have not been performed, except for MoS 2 [79].
As in the case of the lattice parameters, we have used both the vasp and the fleur codes to calculate the band structures of monolayer TMDCs. For the vasp calculations we used the HSE lattice parameters as input. The band structures were calculated in the local density approximation (LDA). SOC was taken into account in the non-collinear magnetic structure approach with the symmetry turned off. The charge density was obtained self-consistently using a 12 × 12 × 1 k-point grid and a 600 eV cutoff energy. The results obtained by this method are shown in rows denoted by "(HSE,LDA)" in Tables 2-9 below. For the fleur calculations the charge densities obtained from the geometry relaxation calculations (see Section 2) were used for further calculation of the band structure and spin expectation values. SOC in fleur is included within the second variational method for the valence electrons, whereas the core electrons are treated fully relativistically. These results are in rows denoted by "(PBE,PBE)" in Tables 2-9 below. One possibility, which we did not explore, is to use the HSE lattice parameters and the HSE06 functional for band-structure calculations, as in Reference [81]. We note that the results of Reference [81] seem to indicate that the HSE06 functional gives larger VB spin-splittings than found experimentally.
In addition to the band structure of the TMDCs, which is our main focus in this work, electron-phonon coupling is also essential in order to understand transport [36,37,38] and relaxation [122] processes. For completeness, we give a brief review of the vibrational characteristics of monolayer TMDCs. Ab initio lattice-dynamics calculations indicate that single layers of the TMDCs MoS 2 , MoSe 2 , WS 2 , and WSe 2 are dynamically stable [123,124,125], in agreement with experiments.
A comprehensive group-theory analysis of the different polytypes and stacking arrangements of few-layer TMDCs is presented in Reference [126]. The symmetry of few-layer structures determines which phonon modes are Raman-active, and therefore provides an important means of characterising samples. As mentioned earlier, monolayer MX 2 has D 3h point-group symmetry (see Table 10 for the character table and irreducible  representations). The six zone-centre optical phonon modes may be classified according to the irreducible representations under which their eigenvectors transform: in the twofold-degenerate E modes the metal atom remains stationary while the chalcogen atoms vibrate in opposite in-plane directions; in the twofold-degenerate E modes the chalcogen atoms vibrate together in-plane in the opposite direction to the metal atom; in the non-degenerate A 1 mode the metal atom remains stationary while the chalcogen atoms vibrate in opposite out-of-plane directions; finally, in the non-degenerate A 2 mode the chalcogen atoms vibrate together out-of-plane in the opposite direction to the metal atom. Of these vibrations, all but the A 2 mode are Raman-active. Only the E and A 2 modes are infrared-active.
DFT-LDA and DFT-PBE results for the phonon frequencies are summarised in Table 1 of Reference [127]. There is a reasonable degree of agreement between the LDA and PBE results, suggesting that the DFT phonon frequencies are accurate. Subsequent theoretical studies [123,124,125] have reproduced the results of Reference [127] for the monolayer. In experimental studies of thin films of WS 2 , WSe 2 , and MoS 2 it is found that modes that were Raman inactive in the bulk become active in thin films and that there are small shifts in the phonon frequencies on going from the bulk to a thin film [128,129,130]. Where comparison is possible, the experimental Raman frequencies of thin films are in agreement with the DFT results.

Band-edge energy differences and spin-splittings
Detailed discussion of the conduction and valence band dispersions in the vicinity of the k-space points of interest (K, Q, Γ, and M ) will be given in Sections 5, 6, 7, and 8. In this section we briefly introduce the various band-splittings and band-edge energy differences that we use to characterise the band structure. An overview of the band structure obtained from DFT calculations is shown in Figure 2. The direct band gap E bg of monolayer TMDCs can be found at the K and −K points of the BZ. Due to the lack of inversion symmetry, all bands are split by the intrinsic SOC except at the time-reversal invariant points M and Γ. We denote by 2∆ vb and 2∆ cb the spin-splitting of the VB and CB, respectively. There are another six minima in the CB that might be important, e.g., for transport or relaxation processes in certain compounds. We denote these points by Q i , i = 1 . . . 6. They can be found roughly half way between the K (−K) and the Γ points. The spin-splitting of the CB at Q i given by 2∆ Q . The importance of the Q i points depends, amongst other things, on the energy difference between the bottom of the CB at the K and Q i points. This energy difference is denoted by E KQ . Looking at the VB now, the energy difference between the top of the VB at K and Γ is denoted by E KΓ . Finally, since it is directly available in recent ARPES measurements [85,86,87], we also record the width of the VB, which we define as the energy difference between the maximum of the VB at K and the minimum that can be found on the Γ-K line.
Certain properties of TMDCs are easier to understand if one considers which atomic orbitals contribute to a given band at a given k-space point. For example, as pointed out in References [26] and [88], the different atomic orbital composition can explain the difference in the spin-splitting magnitude of the CB and VB at the K point. The contribution of individual atomic orbitals to a given band is shown in Figure 3 for the d orbitals of the metal atoms and the p orbitals of the chalcogens (the weights of other atomic orbitals are much smaller). Comparing Figures 3(a) and (b) we find that in general more than one type of atomic orbital contributes to both the CB and the  VB and the weight of the atomic orbitals changes throughout the BZ. Setting up a consistent tight-binding model for TMDCs is therefore more difficult than is the case for, e.g., graphene.

Valence band width D vb
An observable that can be directly compared to experimental ARPES measurements [85,86,87] is the width of the VB D vb . In order to be able to compare the experimental and theoretical results, we define D vb to be the difference between the top of the VB at the K point and the minimum, which lies between the Γ and K points: see Figure  2. (Note that the absolute minimum of the VB is not at this k-space point. However, Reference [86] shows the dispersion only between Γ and K; therefore we use the definition of D vb given above.) Comparison between the calculated and experimental values is given in Table 2.  [85], exfoliated samples on a SiO substrate. b [87], samples grown by chemical vapour deposition on a highly oriented pyrolytic graphite substrate. c [86], samples grown by molecular beam epitaxy on bilayer graphene on top of SiC (0001).
In the case of MoS 2 the experimental values [85,87] suggest that the VB is narrower than the calculated one by ≈ 10%, whereas for MoSe 2 [86] the opposite is true. Since the orbital composition of the VB away from the K point is not purely of Mo d orbital type but also p orbitals of X atoms are admixed (see Figure 3), D vb can be sensitive to interactions with substrates, which are not considered in our calculations and which might explain some of the differences with respect to measurements. 5. Effective model at the K and −K points

K and −K points
The physics around the K and −K points has attracted the most attention both experimentally and theoretically so far. This is mainly due to the exciting optical properties of these materials at the direct band gap, which can be found at the K and −K points. Moreover, it turns out that the effect of SOC is strong at this BZ point, leading to spin-split and spin-polarized bands. Since the K and −K points are connected by time-reversal symmetry, the polarization of the bands has to be opposite at K and −K, i.e., the spin and the valley degrees of freedom are coupled [26]. We start our discussion in Section 5.2 with a basic characterization of the band structure in terms of effective masses and spin-splittings. Then, in Section 5.3, a detailed k · p theory is presented which captures the salient features of the DFT band structure and allows us to understand the results of recent experiments [32,33,34,35].

Basic characterization and material parameters
The aim of this section is twofold. First, we want to point out that there is a difference between the MoX 2 and WX 2 materials regarding the sign of the SOC constant in the CB (for a microscopic explanation see References [47], [54], and [56]). This difference is important for the interpretation of experiments in which properties of A and B excitons [27,30] are compared (for introduction to exciton physics see e.g., [58]). Second, we report effective masses and spin-splittings extracted from our DFT calculations and compare them to experimental results, where available; see Tables 3 and 4.
One of the phenomena that first sparked strong interest in monolayer TMDCs was the pronounced effect of SOC on the VB around the K and −K points. SOC leads to the spin-splitting and spin-polarization of the VB and the energy scale associated with SOC is several hundreds of meVs: see Table 4. SOC in the VB was first studied using DFT calculations [69,88,89,90], but it can be readily understood using, e.g, a tight-binding model and first-order perturbation theory [26,56,57]. An experimental signature of the spin-splitting of the VB is the energy difference of the A and B excitons [27,30].
SOC also affects the CB. This was initially neglected, mainly because in MoS 2 , which is the most widely studied of the TMDCs, it is indeed a small effect and it was assumed that the situation would be similar in other monolayer TMDCs. In general the magnitude of the spin-splitting of the CB is only 7-10% of that of the VB, with the exception of MoS 2 , where it is only ≈ 2%: see Table 3. However, in absolute terms it is an energy scale that can be important at low temperatures and in ballistic samples. Note that the SOC in the CB at the K point is a more subtle effect than in the VB. In the simplest theoretical approximation, which assumes that it is sufficient to consider only the d z 2 atomic orbitals of the metal atoms, the SOC vanishes. DFT calculations, on the other hand, indicate that there is a finite spin-splitting in the CB at the K point [61,69,88,89,90].
As it turns out, the SOC in the CB can be understood in terms of a competition between two contributions [47,54,56,61,91]: i) a first-order contribution from the chalcogen atoms, which have a small, but finite weight [54,55] and ii) a second-order contribution due to the coupling to other bands [47,54,56,61], where the d xz and d yz atomic orbitals have large weights; see Figure 3. Due to this competition the spinpolarisation of the spin-split CBs is different in MoX 2 and WX 2 . Our latest results were obtained using the fleur code, which allows the explicit calculation of the spin expectation value s z in a given band. We find that the spin-split CB with s z > 0 ( s z < 0) is higher (lower) in energy in MoX 2 , while the opposite is true for WX 2 : see Figures 4(a) and 4(b), in which the CBs of MoSe 2 and WSe 2 are shown, respectively. As a consequence, the band with the lighter effective mass is lower in energy for MoX 2 leading to band crossing of the two spin-split bands in the vicinity of the K and −K points [54,47,56]. By contrast, as shown in Figures 4(c) and 4(d), in the VB the sign of s z is the same for both MoX 2 and WX 2 . These differences notwithstanding, there  is a spin-valley coupling in the CB similar to the VB. In Figure 4 we also introduce the notation K vb ) for the higher-in-energy (lower-in-energy) spin-split VB, and similarly for the CB. As a consequence of the spin polarisation of the bands in optical experiments the lowest-energy spin-allowed transition is K cb for WX 2 . Very recently, the first spin-resolved ARPES measurement on bulk WSe2 has appeared [99] and seems to indicate an out-of-plane spin polarisation of the spin-split VB around K and K points. Assuming that the measurements predominantly probe the top layer [99], i.e., effectively a monolayer sample, they are in agreement with the DFT calculations presented here.
The dispersion around the K and −K points is not simply parabolic [61], which has to be borne in mind when fitting the band structure to obtain the effective masses and other band parameters. This can already be appreciated in Figures 2(b) and (c), where a trigonal warping (TW) of the dispersion around the K and −K points can clearly be seen. The TW is more pronounced in the VB than in the CB. In the simplest approximation this can be taken into account by a cubic term in the dispersion. Therefore the dispersion of each spin-split band in the VB and the CB can be described by where the wavevector q = (q x , q y ) is measured from the K point, ϕ q = arctan(q y /q x ), m eff is the effective mass of the given band, and C 3w is a parameter describing the TW. The derivation of E K (q) based on a multi-band k · p model is presented in Section 5.3 and Appendix A. We note that a similar model was recently used in Reference [96].
The values of the m eff and C 3w that we have extracted from our DFT calculations for each band and material are given in Tables 3 and 4. We note that several works have already presented tables of, e.g., effective masses [41,43,64,65,66,67,68] for different monolayer TMDCs. However, the effects of SOC have often been neglected leading to, e.g., the conclusion that the effective masses of the spin-split VBs are the same. Recent experimental evidence shows that this is not the case [86]. Moreover, due to the presence of the TW, some care has to be taken when defining the effective mass and, especially, when choosing the fitting range that is used to obtain it from a DFT band structure. All our DFT band-structure calculations were performed along the Γ-K-M line in the BZ. We first fitted m eff , i.e., we set C 3w = 0 in Equation (1). The fitting range corresponded to 5% of the Γ-K distance. The dispersion over such range was considered to be isotropic and the difference in the effective masses along K-Γ and K-M was neglected. Therefore the effective masses shown in Tables 3 and  4 characterise, strictly speaking, a rather narrow vicinity of the band edge. The nonparabolicity of the band structure and the trigonal distortion of the constant energy contours, described by the second term in Equation (1), was taken into account in a second step, whereby Equation (1) was fitted over a wider range (typically ≈ 10% of the Γ-K distance), but m eff , obtained in the previous step, was kept fixed. This twostep fitting was needed to obtain coherent parameter sets between the simple approach outlined here and a more accurate model presented in Section 5.3. Further details of the fitting procedure are discussed in Appendix B. Looking at Tables 3 and 4 one can see that the effective masses and spin-splittings obtained from the two different DFT calculations are in almost perfect agreement, while there are some differences in the extracted values of C 3w .
Considering first the CB, the extracted band parameters and SOC splittings 2∆ cb for different monolayer TMDCs are shown in Table 3. In addition we show the charge density n cb at which the upper spin-split CB K (1) cb starts to be populated. This charge density is calculated using the effective mass of the K (2) cb band given in Table  3 and assuming a simple parabolic dispersion (i.e., neglecting C 3w ), which is a good approximation in the CB. Note that typical charge densities achieved by gating in Table 3. Band dispersion parameters and spin-splittings at the K and −K points in the CB from DFT calculations. m cb ) band, and similarly for C . m e is the free electron mass. n cb is the electron density above which the upper spin-split CB starts to fill. MoS 2 are reported to be ∼ 4 · 10 12 cm −2 -3.6 · 10 13 cm −2 [92] a few times 10 12 cm −2 for monolayer samples [93] and few-layer samples [94]), and up to 10 14 cm −2 in few-layer WS 2 using ionic liquid gating [95]. To our knowledge there are no direct measurements of ∆ cb or m cb yet. Turning now to the VB, the band parameters and SOC splitting 2∆ vb obtained from our DFT calculations are shown in Table 4. In the case of MoSe 2 , very recent highresolution ARPES measurements [86] allow for a direct comparison with the calculations, because the difference between the effective masses of K (1) vb and K (2) vb could be directly observed. We show two theoretical values for the effective masses in the VB of MoSe 2 . The first one is obtained using the fitting procedure described above, i.e., by averaging the values along the K-Γ and K-M directions. The second value, shown in parenthesis, is obtained by following the fitting procedure that was used for the experimental data [97]. This latter procedure involves fitting only along the K-Γ direction, and a fitting range of ≈ 13% of the K-Γ distance. One can see that the theoretical and experimental effective masses that were obtained using the same fitting range are in good agreement. Moreover, the calculated value of 2∆ vb also corresponds rather well to the measured one. MoS 2 is the only other monolayer TMDC for which ARPES measurements are Table 4. Effective masses and spin-splittings at the K point in the VB from DFT calculations. m vb ) band, and similarly for C . m e is the free electron mass. The values in brackets were obtained using a slightly different fitting range, as explained in the text. Experimental values are shown in rows denoted by "Exp". available to extract the effective mass. However, the ARPES data of Reference [98] do not resolve K (1) vb and K (2) vb separately; therefore the reported effective mass is the average of m (1) vb and m (2) vb . Taking into account the experimental uncertainty, our results are in reasonable agreement with the measurements of Reference [98]. The available data for MoS 2 and MoSe 2 suggest that DFT can capture the VB effective masses quite well even without GW corrections, such as those found in Reference [45]. For monolayer WS 2 and WSe 2 , to our knowledge, no ARPES data are yet available.
In optical experiments the difference of the A and B exciton energies are usually identified with 2∆ vb providing the results shown in Table 4. We note that there are two assumptions behind the identification of the A and B exciton energy difference with 2∆ vb : i) that the spin-splitting in the CB is negligible and ii) that the binding energies of the A and B excitons are the same. Regarding i), one can see in Table 3 that ∆ cb is small, but finite, and for quantitative comparisons between theory and experiment it should not be neglected. As for ii), we note that the binding energy of the A and B excitons depends on their reduced mass, which, according to Table 4, should be different for the different exciton species. With these caveats the agreement between the calculations and the experiments is qualitatively good, especially for MoS 2 and MoSe 2 .
Comparing the DFT-calculated effective masses in Tables 3 and 4 for the VBs and CBs that have the same spin-polarisation, one can observe that there is no electronhole symmetry in the band structure. The first experimental evidence to support this observation, coming from magnetoluminescence experiments, has appeared very recently [33,34,35]. Regarding the experimental relevance of TW, it has been argued [32] that it leads to measurable effects in the polarisation of electroluminescence in p-n junctions. We note that due to the heavier effective mass in the VB and the larger values of C 3w , the TW is more pronounced in the VB than in the CB. In the latter a simple parabolic approximation is often adequate.
We finish Section 5.2 with a brief discussion of the quasiparticle band gap E bg,K , which we define as the difference between the maximum of the K (1) vb and K (2) cb bands at the K and −K points. DFT calculations for monolayer TMDCs underestimate the band gap (see Table 5) and its evaluation requires the use of GW methodology [69,43,65,44,45]. Experimental evidence that supports the conclusions of the GW calculations is now also emerging. Apart from its fundamental importance, the main reason for discussing E bg,K and showing our DFT results is that E bg,K enters into the fitting procedure that we use to obtain the parameters of the k · p Hamiltonian that describes the dispersion in the vicinity of the band edge. The details of the k · p model and the fitting procedure are given in Section 5.3 and Appendix B. As one can see, our DFT calculations significantly underestimate the experimental quasiparticle band gaps. We also note that in heavily doped samples, which were used in the ARPES measurements [98,86], the observed band gap is reduced with respect to results obtained by other methods [70,80,24,104], hinting at the crucial importance of screening in monolayer TMDCs.

k · p Hamiltonian
We now present a low-energy effective k · p Hamiltonian that describes the coupled dynamics of the VB and CB. Part of the theory was previously discussed in References [61] and [47]; in the present work we both overview and extend these earlier results. To obtain a model that captures the most important features of the dispersion of the VB and CB one can start from a seven-band model, which was introduced in Reference [47]; a brief recap is given in Appendix A. An effective low-energy Hamiltonian can be derived by systematically eliminating all degrees of freedom other than the ones corresponding to the VB and CB using Löwdin partitioning [110]. We keep terms up to third order in the off-diagonal coupling elements of the original seven-band model and use the spinful basis {|Ψ vb , s , |Ψ cb , s }, where |Ψ vb (|Ψ cb ) are spinless Bloch wave functions in the VB (CB) and |Ψ b , s = |Ψ b ⊗ |s , with b = {cb, vb} and s = {↑, ↓} denoting the band spin degree of freedom, respectively. One finds that the low-energy effective Hamiltonian is the sum of the following terms: i) The free-electron term H 0 = 2 q 2 2me (1 2 ⊗s z ), where 1 2 is a unit matrix in the electronhole space, s z is a spin Pauli matrix, and m e is the free electron mass. Here and in Equations (4b)-(4f) the wavevector q = (q x , q y ) is measured from the K or −K points. We note that H 0 is usually neglected in the GaAs literature on account of the light effective mass in this material, but here we want to keep it.
ii) The SOC Hamiltonian H τ,s so , which contains the diagonal and q-independent contributions of the SOC. It reads i.e., it is diagonal in spin space and is proportional to the Pauli matrix s z (for further details see Appendix A). H τ,s so describes the spin-splittings of the CB and the VB, which are due to the absence of inversion symmetry in monolayer TMDCs. Since H τ,s so is diagonal, one can also write it in terms of the eigenvalues s = ±1 of s z ; we will use the two notations interchangeably. Moreover, the index τ = 1 (τ = −1) denotes the valley K (−K). Where it is more convenient, we will also use the matrix τ z which acts in the valley space. In the VB, the parameter ∆ vb that describes the strength of the SOC can always be taken to be positive. As explained in Section 5.2, the situation is more complicated in the CB [54,47,56], because DFT calculations show that, in the case of MoX 2 , the spin-split bands cross close to the K and −K points, while there is no such band crossing for WX 2 . This can be understood in terms of ∆ cb having opposite signs in MoX 2 and WX 2 .
iii) Finally, the k · p Hamiltonian H τ,s k·p in Equation (2) is given by where Here q ± is defined as q ± =q x ± iτ q y , ϕ q = arctan(q y /q x ), ε vb and ε cb are band-edge energies, γ τ,s , α τ,s , β τ,s , κ τ,s , η τ,s , and ω s are material parameters discussed below. H τ,s k·p is a generalisation of the results given in Reference [61] for the case in which the material parameters depend on the SOC.
In general all off-diagonal material parameters appearing in H τ,s k·p are complex numbers such that for τ = −1 (−K valley) they are the complex conjugate of the τ = 1 (K valley) values. Concrete values of the material parameters for each MX 2 material can be obtained by, e.g., fitting a DFT band structure, see Table 6. Note however, that the fitting procedure (see Appendix B) yields real numbers for each parameter. We now briefly discuss each of the terms [Equations (4b)-(4f)].
i) Terms up to linear in order inq + andq − can be found in Equation (4b). H τ,s D is basically the massive Dirac fermion model introduced in Reference [26]. It describes an isotropic dispersion around the band edge and it does not break the electron-hole symmetry. The value of γ s,τ also depends on the SOC, but the Löwdin-partitioning calculations suggest that this dependence should be weak. This is indeed what we have found from fits to the DFT band structure. Therefore in the following we suppress both the spin index s and, since γ is taken to be a real number, the valley index τ .
iii) Off-diagonal terms quadratic inq + andq − are given in Equation (4d). H τ,s 3w , in combination with H τ,s D , gives a contribution to the TW of the energy contours that can be observed in Figures 2(b) and (c). [For further details see Equation (B.1) in Appendix B]. The TW is expected to play an important role in the explanation of recent electroluminescence experiments [32]. Moreover, it was observed in ARPES measurements [87]. iv) Off-diagonal terms cubic inq + andq − appear in Equation (4e). H τ,s cub,1 is important for obtaining a good fit to the DFT band structure away from the K point in a two-band model that describes the coupled dynamics of the VB and CB. In combination with H τ,s D they contribute to the TW of the bands; see Equation (B.1) in Appendix B.
v) Diagonal terms cubic inq + andq − . In some cases it is more convenient to work with a model that gives the dispersions of the VB and CB separately. Cubic terms in q are needed to capture the non-parabolicity of the bands, and such a model is given by Equation (1). It can easily be obtained by applying Löwdin partitioning to Equation (2) and eliminating either the electron or the hole degrees of freedom. In this case, for consistency, the term H τ,s cub,2 in Equation (4f) also has to be taken into account.
The eigenstates and eigenvalues of the k · p Hamiltonian (2) can also be used as a starting point for analytical calculation of the Berry curvature [108]. The Berry curvature is relevant for the quantum transport characteristics of TMDCs, such as the valley Hall effect [26] and weak localisation [157], while a related quantity, the spin Berry curvature [109], gives rise to a finite spin Hall conductivity for moderate hole doping.
Finally, we show the k · p parameters obtained from fitting of the DFT band structure (see Table 6) using the model that explicitly contains the coupling between the VB and the CB. In this case the diagonal cubic term [Equation (4f)] is not important for obtaining a good fit to the band structure and therefore the ω s parameter is not shown. Close to the band edge the k · p parameters given in Table 6 reproduce the effective masses shown in Tables 3 and 4. The details of the fitting procedure are given in Appendix B. Since the effective masses and C 3w parameters extracted from the (HSE,LDA) and (PBE,PBE) approaches are rather similar, we only show results that are based on (HSE,LDA) DFT band-structure calculations. Due to the SOC all parameters, with the exception of γ, are different for different spin indices s. Since the Hamiltonian of Equation (4a) is diagonal in the spin space, i.e., it describes the coupled dynamics of the VB and CB having the same spin, it is convenient to introduce the notation s =↑ (s =↓) for s = 1 (s = −1). Regarding the correspondence between the notation used in Section 5.2 and here, note that the order of the bands with ↑ and ↓ polarisation in the CB is different for MoX 2 and WX 2 . Therefore in the VB the upper index (1) ( (2)) is equivalent to ↓ (↑), but in the CB the relation depends on which material is considered. We note that the parameter γ can, in principle, also be obtained directly as a momentum matrix element between the Kohn-Sham wave functions of the VB and CB. For these calculations we used the castep code [111], where the necessary plane-wave coefficients of the wave functions at the band edges are readily accessible. These values are denoted by |γ KS | in Table 6. On the one hand, the good agreement between |γ| and |γ KS | indicates the consistency of our fitting procedure. This is not trivial, because the fitting involves a non-linear function of the k · p parameters. On the other hand, one has to bear in mind that |γ| is obtained such that it would give the best fit to the DFT band structure over a certain range in the BZ. Therefore it may differ from the value of |γ KS | that is calculated at a single point of the BZ. The valley index τ is suppressed in Table 6 because, as mentioned above, from the fitting procedure we obtain real numbers for the off-diagonal terms.
As explained in Appendix B, our fitting procedure involves the quasiparticle band gap E bg,K . For this reason two sets of k · p parameters are reported in Table 6: one in which we used E bg,K values obtained from our DFT calculations and one in which we used E bg,K values found in GW calculations; see Table 5. In the latter case we make the assumption that the bands above the Fermi energy are rigidly shifted upwards in energy such that the effective masses and the TW in the VB and CB remain the same. We believe that this is a reasonable assumption because the available experimental evidence (see Tables 4 and 9) suggests that, at least in the VB, the effective masses are captured quite well by the DFT calculations. In addition to the K and −K points, there are six other minima in the CB which may be important for, e.g., relaxation processes. We denote the BZ points where these minima are located by Q i , (i = 1 . . . 6); they are also known as Λ points [see Figures 2 and 5(b)]. We note that phonon scattering between the K and −K points and Q i points is symmetry-allowed [122] and that, depending on the energy difference E KQ (see Figure 2), the electron mobility may be significantly affected by these scattering processes [36,37,40]. However, as we will show, understanding the SOC at the Q i points is also important when considering the possible scattering processes, a fact which seems to have been overlooked in some recent publications. We start in Section 6.2 with a basic characterisation of the band structure in terms of the effective masses and point out an important effect of SOC on the spin polarisation of the bands. A detailed k · p theory is given in Section 6.3.

Basic characterisation and material parameters
Let us consider the Q 1 minimum, which can be found along the Γ-K direction [see Figure  2(c)]. We choose k x to be parallel to the Γ-K direction, while k y is perpendicular to it. Neglecting SOC for a moment, our DFT calculations show that, close to the Q 1 point, the energy contours are to a good approximation ellipses whose axes are parallel to k x and k y [see Figure 5 is quadratic with different effective masses m 0 Q,x and m 0 Q,y along k x and k y : where the wavenumbers q x and q y are measured from the energy minimum of the dispersion (see Section 6.3 for details). As one can see in Table 7, the ratio of the effective masses is m 0 Q,y /m 0 Q,x ≈ 2 for MoX 2 and m 0 Q,y /m 0 Q,x ≈ 1.3-1.8 for WX 2 . SOC has two major effects [see Figure 5 Similarly to Section 5.2, we introduce the notation Q i ) for the higher-inenergy (lower-in-energy) spin-split CB at the Q i point [see Figure 5(b)]. The basic characterisation of these bands therefore requires two effective masses for each of the two spin-split bands and the spin-splitting energy 2∆ Q . In addition, it is also important to know the energy difference E KQ between the band extrema of the Q (2) i and K (2) cb bands. These parameters, obtained by fitting to our DFT band structures, are shown in Table  7. The fitting range we used was ≈ ±7.5% of the Γ-K distance around the Q point in the k x direction and roughly half of that in the k y direction. Looking at Table 7 one can see that the effective masses obtained from the two DFT calculations are again in good agreement, while small differences can be seen in the results for 2∆ Q . However, there are noticeable differences in the energy separation E KQ between the bottom of the CB at the Q and the K points, which we ascribe to the different lattice constants used in Table 7. Material parameters at the Q point. n Q is the electron density above which the carriers start to populate the Q valleys.   [98]. Note, however, that the ARPES measurements in Reference [98] were performed on potassium-intercalated samples, and the effects of the intercalation on the band structure of TMDCs have not yet been studied in detail. We also note that E KQ , in contrast to the band gap E bg , appears to be less sensitive to GW corrections [114] if the latter calculations are well converged, and therefore we expect that our DFT results are at least qualitatively correct.
We have also calculated the carrier density n Q at which the Fermi energy, measured from the bottom of the K-point valley in the CB, reaches the bottom of the Q-point valley; see Table 7. In our calculation we assumed a simple parabolic dispersion for the CB in the vicinity of K, where the effective masses of K (1) cb and K (2) cb are given in Table 3. Our results suggest that for MoX 2 it will be difficult to achieve the doping levels needed to populate the Q (1) i valleys, but for WX 2 the required doping levels are more realistic. As noted in Reference [57], the valley-spin coupling is present not only in the K and −K valleys, but also in the Q i valleys, and this may have experimental consequences. The calculated spin polarisation of the CB between the K and the Q 1 point is shown in Figure 6 for MoSe 2 and WSe 2 . One finds that despite the band crossing(s) between the K and Q 1 points, for MoX 2 the spin-polarisation of the Q 1 band is the same as the spin polarisation of the K (2) cb band [see Figure 6(a)]. For WX 2 , however, due to the band crossings, the spin-polarisation of K (2) cb is opposite to the spin polarisation of Q cb to Q 1 , Q 3 , and Q 5 is allowed, while scattering to Q

, and Q
(1) 6 is, strictly speaking, also allowed but should be suppressed with respect the former processes due to the relatively large spin-splitting 2∆ Q .

k · p Hamiltonian
Due to the low symmetry of the Q i points in the BZ and because there are many nearby bands in energy, there is a large number of band-overlap parameters that would need to be taken into account in a detailed multi-band k · p model. Therefore it is more difficult to develop such a theory and it would offer less insight. Nevertheless, a low-energy effective k · p Hamiltonian can be derived with the help of the theory of invariants (for a recent discussion see, e.g., References [115] and [116]). The pertinent symmetry group is C 1h ; for convenience, its character table is shown in Table 8 [59].
As in the case of the K and −K valleys, the Q i -point minima are pairwise connected by time-reversal symmetry and to describe this one can introduce the matrix τ z , whose eigenvalues, τ = ±1 label individual members of the pairs of valleys. As an example, let us consider the Q 1 (τ = 1) and Q 4 (τ = −1) minima, which can be found along the Γ-K and Γ-(−K) directions, respectively [see Figure 2(c)]. This direction is parallel to the k x component of k. Including up to second-order terms in k and taking SOC into account the effective Hamiltonian reads where k x and k y are measured from the Γ point, 1/m τ,s (Q,x,y,xy) are effective masses, s = ±1 are the eigenvalues of the spin Pauli matrix s z . Furthermore, E Q is the bandedge energy if SOC is neglected, ∆ Q is the spin-splitting at Q, and a 1,2 and b 1,2 are material parameters to be discussed later. Since we are going to develop a theory in which the dispersion is parabolic, in contrast to Section 5, we will not keep track of the free-electron contribution explicitly.
To simplify the discussion, let us first neglect the spin degree of freedom. Then ∆ Q = a 1,2 = 0 and m τ,s (Q,x,y,xy) = m 0 (Q,x,y,xy) . Since close to Q 1 (Q 4 ) the energy contours are, to a good approximation, ellipses whose axes are in the k x and k y directions [see Figure 5(a)], one finds that 1/m 0 (Q,xy) ≈ 0 and introducing the wavenumbers q x and q y measured from k = (τ k Q , 0), i.e., from the Q 1 (Q 4 ) point, one can write whereẼ KQ measures the energy difference with respect to the K point in the absence of SOC. The effective masses m 0 Q,x and m 0 Q,y are in general different. Taking SOC into account, H τ,s Q [Equation (6)] can be re-written in the following form: where E KQ is defined in Figure 2. One can see that SOC has the following effects: i) it splits the bands and opens a gap ∆ Q between the spin-up and spin-down bands; ii) it shifts the minima of the spin-split bands off from the k = (τ k Q , 0) point; iii) it makes the effective masses of the spin-polarised bands different, so they are given by 1/m τ,s (Q,x,y) = 1/m 0 (Q,x,y) − τ s/δm (Q,x,y) .  Figure 5(b) taking WSe 2 as an example, where these effects are most clearly seen. The material parameters m τ,s Q,x,y , q Q,x , q Q,y , ∆ Q , and E KQ can be obtained from, e.g., DFT calculations; see Section 6.2. We find that q Q,y is zero within the precision of our calculations and q Q,x is always very small.

Γ point
Next we consider the band structure at the Γ point. There are three main motivations to include the Γ point in our work: i) there is a local maximum in the VB at the Γ point, which could be observed in recent ARPES measurement [85,86] and therefore it is of interest to compare the experimental and calculated effective masses; ii) it was theoretically predicted [45] and recently verified experimentally [70,106] that in addition to excitons at the K point, in MoS 2 excitons can also be formed in the vicinity of the Γ point. It was argued that these "C-excitons" are qualitatively different from the ones at the K point because they arise from an effectively one-dimensional energy minimum [70] in the "optical band structure" (for the exact definition see below); and iii) finally, understanding of the VB and CB behaviour at the Γ point is important in the interpretation of scanning tunnelling microscopy (STM) experiments [80,119].

Basic characterisation and material parameters
We start the discussion with the VB maximum at the Γ point (VBMG). The spinsplitting is at Γ point is zero, and in the VB it remains negligible over a considerable region of k space [see Figure 2(a)]. Moreover, to a good approximation the dispersion around the VBMG is isotropic (see Figure 2(b) and Section 7.3) and parabolic. Therefore it can be described by which is characterised by a single effective mass m vb Γ and the energy E KG , which is the energy difference between the maximum of the K (1) vb band at the K point and the VBMG. The values of m vb Γ obtained from fitting the results of our DFT calculations are given in Table 9 and experimental results, where available, are also shown. The effective masses m vb Γ were obtain by fitting the band structure along the Γ-K direction in a range of ≈ 21% of the Γ-K distance. The calculated m vb Γ values are in reasonable agreement with the available experimental results. There are, however, noticeable differences between the theoretical and experimental E KG values, which we attribute to substrate effects. Note that the weight of the chalcogen p orbitals is substantial at the Γ point (see Figure  3 or Reference [55]) so that one can expect a stronger interaction between the substrate and the electronic states. Interestingly, the m vb Γ parameter does not seem to be affected as strongly as E KG by the substrate. Comparing Tables 4 and 9, one would expect the VBMG to be the most important for the transport properties of MoS 2 because it is probably quite close in energy to the maximum of the K vb band (which we denote by VBMK1; we also introduce the label VBMK2 for the maximum of the K (2) vb band), thus facilitating scattering processes in the VB that do not require spin-flips. Furthermore, according to our DFT calculations VBMG lies in between VBMK1 and VBMK2 for MoS 2 and WS 2 , while it lies below VBMK2 for MoSe 2 and WSe 2 .
Very recently, there has been considerable interest in optical transitions near the Γ point for MoS 2 [45,70,117]. The exciton species associated with this transition is named "C-exciton" and its existence is explained by the presence of a minimum in the optical band structure around the Γ point. Here, following References [70] and [118], the optical band structure is defined as the difference between the dispersions of the CB and VB.
To obtain an insight into possible optical transitions at higher energies than the ones at the fundamental band gap, we plot the optical band structure for monolayer TMDCs in Figures 7 and 8 over the whole BZ. For simplicity, SOC is neglected in these calculation. A clear "gear-shaped" minimum [70] is noted both for MoS 2 and WS 2 around the Γ point [Figures 7(a) and 8 (a)]. For a more quantitative understanding of the interband transitions one also needs to consider the optical density of states (optical DOS), which is defined as the density of states of the optical band structure (the terminology "joint density of states" is also used; see, e.g., [118]). A peak in the optical DOS corresponding to the minimum in the optical band structure is present at  Figure 8(c)]. However, other peak(s) and a wide shoulder extending into higher energies can also be seen in the optical DOS. We attribute these features to saddle points in the optical band structure, which can be observed, e.g., along the Γ-K line for WS 2 and to the saddle points at the M point in the optical band structure of all four MX 2 materials. These observations motivate us to have a closer look at the band structure at the M point as well, which is presented in Section 8. Finally, we emphasise that for a quantitative understanding of the optical band structure and the interband optical transitions the effects of SOC are also important. In general, they lead to spin-splitting of the bands (e.g., along the Γ-K line), or splitting of the Van Hove singularity (see Section 8.1) and the energies of these splittings may be comparable to or larger than the linewidth of the optical transitions.

k · p Hamiltonian
As in previous sections, we use group theory to obtain effective k·p Hamiltonians for the VB and CB. Similarly to the K point, it is possible to set up a multi-band k · p model. We have found, however, that the number of necessary bands, even if one neglects SOC, is quite large and, as will be shown later, terms up to fourth order in k need to be taken into account in order to capture the features of the band structure related to the C-exciton terms. Therefore we present here only a simplified discussion of the problem; a more complete theory is left for a future work. As we will show, important insight can be gained from the spinless case, i.e., in the discussion that follows we will neglect SOC.
The bands of interest are the VB, the (doubly degenerate) CB [shown by green lines in Figure 9(a)] and the first (doubly degenerate) band above the CB, which we denote by CB+1 [black lines in Figure 9(a)]. We will rely on group-theoretical arguments, which are very convenient at the Γ point, where, as mentioned above, several atomic orbitals contribute with significant weight to each band. The pertinent symmetry group is D 3h and the character table is shown in Table 10.
Symmetry analysis of the contributing atomic orbitals implies (see, e.g., Table IV in Reference [61]) that the VB at the Γ point belongs to the A 1 irreducible representation of D 3h . As already given in Equation (9), up to second order in q, the dispersion is parabolic and isotropic [see Figure 2 Table 9. Along the Γ-K direction, the spin-splitting of the VB is small up to wavevectors corresponding to about half of the Γ-Q distance. This is due to the fact that in the vicinity of Γ the d z 2 atomic orbitals of the metal and the p z atomic orbitals of the chalcogen atoms contribute with large weight to the VB (see Figure 3 and [54,55]). Along Γ-M all bands remain spin-degenerate due to symmetry; see Section 8.1. The spin-splitting of the VB is therefore suppressed around the Γ point.
Turning now to the CB [shown by green lines in Figure 9(a)], at the Γ point it is doubly degenerate and antisymmetric with respect to the horizontal mirror plane σ h of the crystal lattice. In group theoretical terms, it corresponds to the E irreducible representation (irrep) of D 3h . Since the VB is symmetric with respect to σ h , one can show using group-theoretical arguments that the optical matrix element between the VB, which has A 1 symmetry, and the CB, which has E symmetry, is zero at the Γ point.
However, as shown in Figure 9, due to band crossings one of the degenerate CB+1 bands becomes the CB at some distance from Γ. The doubly degenerate CB+1 band belongs to the 2D E irreducible representation of D 3h . This irrep is symmetric with respect to σ h , and optical transitions between bands of A 1 and E symmetries are allowed. Therefore as a starting point for studying the optical transitions in the vicinity of the Γ point one has to describe the CB+1 bands. Up to second-order terms in the wavevector k, the effective Hamiltonian describing the E bands in the vicinity of Γ reads: where 1 2 is a 2 × 2 unit matrix, and φ k is the argument of k x + ik y (here the wavevector components k x and k y are measured from Γ). We also keep explicit the free-electron describes the coupling of the CB+1 bands to other remote bands with the same E symmetry, while H (2) wr captures the coupling of the CB+1 bands to other remote bands with A 1 symmetry. In contrast to the VB, one can see that H (2) wr leads to a hexagonal distortion of the energy contours of the CB+1 bands already in second order of k. Looking at Equation (10c) one can also note that, e.g., along the Γ-K line the off-diagonal and one of the diagonal terms become zero. Therefore Equation (10c) alone would suggest that one of the E bands is dispersionless. Since the dispersion of the higher-in-energy E band is indeed very flat along Γ-K [see Figure 9(a)], we expect that H (2) d largely cancels H 0 . SOC, as illustrated in Figure 9(b), has two main effects: i) At the Γ point it leads to a splitting of the otherwise degenerate states. Therefore, instead of four-fold degeneracies, which would follow from taking into account the spin but not the SOC there are only two-fold degeneracies.
[For the E bands in MoS 2 the splitting is too small to be seen on the scale of Figure 9(b)].
ii) Close to the Γ point the band crossings between the E and E bands are turned into avoided crossings.
One can observe, however, that beyond these avoided crossings the dispersion of the spin-split CB follows that of the spinless CB quite closely. This is remarkable for the following reason: it has been argued [70] that the existence of the C-exciton is related to a minimum in the optical band structure. Looking along the Γ-K or Γ-M lines, the minimum in the optical band structure can be found for k values where the spinful CB closely follows the lower-in-energy E band. We expect therefore that, theoretically, the starting point for describing the C-exciton physics in an effective-mass approximation would be to extend the model shown in Equations (9) and (10a) by terms that contain higher powers of k, especially for the E bands, where these corrections become important closer to the Γ point than is the case in the VB. Neglecting the coupling between the two E bands and considering only the one lower-in-energy band that becomes the CB, the terms up to fourth order in k that need to be added to the dispersion are (11) where the constants C (3) and C (4) 1,2,3 can be obtained from fitting the band structure. Looking at Figures 7 and 8, this approach appears to be most useful for MoS 2 and WS 2 , where a clear minimum in the optical band structure in the vicinity of Γ can be seen. However, the exact location of the minimum in the optical band structure may also depend on the SOC, which was not taken into account in Figures 7 and 8. This introduces additional complexity into the problem because, as noted in Section 6, the results at the Q point are somewhat sensitive to the computational method that is used. A detailed discussion of the optical band structure is therefore left for a future work.

Γ point wave functions and STM measurements
The shape and extent of the VB and CB wave functions at the Γ point can also play an important role in the interpretation of STM measurements. Since there is a growing experimental interest [80,119,120,121] in STM studies of monolayer TMDCs, we give a brief account of calculations that can be used to interpret STM measurements.
We first focus on the STM maps that one can obtain using a tip with a curvature radius larger than atomic distances at scanning distances comparable to or larger than the lattice constant. In this case the current is dominated by electrons tunnelling from the metal with the largest k z momentum component at the energy given by the scanning voltage. Therefore the in-plane momentum components of the tunnelling electron can be neglected: k x,y → 0. On one hand, this implies that the real space 2D maps of the tunnelling current should reflect the vertical extent of the Γ-point wave functions in the CB and VB. On the other hand, electron tunnelling into the band edges, which are at the K and −K points, can take place as a two-step process: first, the electron tunnels into a virtual state close to the Γ point in the corresponding band, then it emits a BZ-corner phonon to scatter into the final state near the band edge. The expected I-V characteristic therefore should have a tunnelling gap in the current of magnitude of the phonon energy, counted from the Fermi level in a doped 2D semiconductor or the band edge in an undoped one.
STM images of bulk MX 2 can be simulated from first-principles. Since in STM measurements one detects the tail of either the VB or CB wave function, depending on whether electrons or holes are injected into the material, one has to determine the decay rate of the wave function of the relevant states in different parts of the unit cell. This can be achieved by numerical differentiation of the logarithm of the square modulus of the band-decomposed wave function along the z direction (i.e., perpendicular to the sheet). For this purpose we use a trilayer geometry to model the surface of the bulk material, since we do not expect the inter-layer interaction to affect the tail of the wave functions severely. In these calculations we used the optB88 van der Waals density functional [131]. This functional should provide a significantly better description of the interlayer interaction than the LDA or PBE functionals. Figure 10 illustrates the decay rate of the VB of MoS 2 . Three points are highlighted: the position of the metal atom (M), the chalcogen atom (X), and the centre of the hexagon formed by three M and three X atoms on the surface (H). Large tunnelling currents occur when the decay rate is low. For example, in the case of the VB, the tunnelling current is dominated by the contribution of the sublattice where the metal atom is located in the VB. Note that the centre of the hexagon is also quite bright; this is due to constructive interference between the p x and p y orbitals of the chalcogen atoms. Table 11 summarises the decay rates at the three notable positions in the unit cell for the four MX 2 materials studied in this work. It can be used to explain which sublattice is expected to dominate the tunnelling current in a particular MX 2 material.
One can also envisage an alternative STM arrangement where the tunnelling current  is determined by coupling with a single atomic site at the end of the tip brought to atomic/subatomic distances from the 2D material. In this case, momentum transfer and momentum conservation are not problems for the tunnelling electron; hence, the tunnelling spectrum may reflect the structure of the electronic wave function at the band edges in the BZ corners. However, in this case, the actual current maps would be affected by the form of the atomic orbital of the last atom in the tip and analysis of such details lies outside the scope of this review.

The M point: spin-orbit splitting of the Van Hove singularity
In Section 7 we have already shown that optical transitions in monolayer TMDCs are expected to occur not only at the K and −K points, but at other points in the BZ as well. Indeed, a strong light-matter interaction was observed in Reference [20] and attributed to Van Hove singularities in the electronic density of states. Moreover, strong absorption beyond the energy range of visible light has been found in MoS 2 [106] and high-energy optical transitions (in the range of 1.5-9 eV) have been studied using ellipsometry in Reference [107]. Motivated by these observations and by the fact that, according to Figures 7 and 8, the optical DOS is finite at energies that correspond to transitions at the M point of the BZ, we briefly discuss the dispersion of the VB and CB at the M point. Optical transitions associated with a saddle point in the dispersion have also been studied recently in monolayer graphene [132].  Figure 11. Band structure of a monolayer TMDC at the M point obtained from DFT calculations a) without taking SOC into account; and b) with SOC. In a) the labels above the bands denote the pertinent irreps of the group C 2v . The actual calculations were performed for WS 2 using the (PBE,PBE) approach.

Basic characterisation and k · p Hamiltonian
monolayer TMDCs considered here the energy difference between the CB and the CB+1 (VB and VB−1) bands is rather small compared to the band gap; the difference between the CB and CB+1 is around 0.5 eV and the difference between the VB and VB−1 is 0.15-0.3 eV. Therefore, regarding optical transitions, the situation at the M point is different from the K point, where the CB and the VB are well separated in energy from all other bands. It is also different from the situation encountered at the Γ point, where the CB was antisymmetric, while the VB and the CB+1 bands were symmetric with respect to the horizontal mirror plane. Here all four bands are symmetric and in-plane polarised electromagnetic radiation can, in principle, induce transitions between them. The discussion of the band dispersion at the M point is simplified if one introduces the local coordinate system shown in Figure 12. Here both q x and q y are measured from the M point, the former being parallel to the K-M direction, the latter to the Γ-M direction. Similarly to the Q point, we content ourselves with the construction of a k · p Hamiltonian based on the theory of invariants. The six M i points in the BZ are pairwise connected by time-reversal symmetry. To describe this one can introduce the matrix τ z , whose eigenvalues, τ = ±1 label individual members of the pairs of M points. In the simplest approximation, the Hamiltonian of all four bands of interest is One can see that the dispersion is parabolic and characterised by different effective masses m τ,s M,x and m M,y along the M -K and M -Γ directions, respectively. To understand the implications of Equation (12), let us first neglect the SOC, i.e., we set ∆ M = 0 and m τ,s M,x = m M,x . Looking at Figure 11(a) one can notice that for the CB (denoted by B 2 ) the effective masses m M,x and m M,y have the same sign. For the VB (denoted by A 1 ), however, their sign is different. Similar conclusions hold for the CB+1 and the VB−1 bands. Therefore in the optical band structure one has a saddle point in the dispersion and consequently a Van Hove singularity.
SOC, as shown in Figure 11(b), has two main effects: i) It leads to a linear-in-q x splitting of the bands along M -K, while the bands remain spin-degenerate along M -Γ. This means that the saddle point in the optical DOS of any two bands will also be split along M -K. In addition, the effective mass m τ,s M,x becomes spin dependent.
ii) It turns band crossings into avoided crossings.
As shown in Figure 11(b), the linear-in-q x splitting of the bands is a rather good approximation close to M . However, the situation is complicated by the fact that the CB+1 (VB−1) band is quite close in energy to the CB (VB). SOC couples the CB and CB+1 (VB and VB−1) bands and leads to avoided crossings between them. Therefore a more complete description would require a model similar to the K point, where the coupling of nearby bands is explicitly taken into account. We leave the construction of such a model to a future work. The Hamiltonian of Equation (12) can be constructed using the theory of invariants. The symmetry group at the M point (and along Γ-M ) is C 2v , which includes the following symmetry operations: a twofold rotation C 2 around the Γ-M direction, a reflection σ v with respect to the q y -q z plane, and the reflection σ h with respect to the q x -q y plane. The character table of C 2v is shown in Table 12, where the relevant basis functions, in the chosen coordinate system, are also given. The Hamiltonian, which Table 12. Character table for the group C 2v . Basis functions for a given irrep are also shown. R x,y,z denotes the angular momentum components.
can be constructed with the help of Table 12 and which is at most second order in the wavenumbers q x and q y , is given in Equation (12). The symmetries of the individual bands at the M point, which can be deduced by, e.g., considering which atomic orbitals contribute to a certain band, are indicated in Figure 11(a).
It is important to note that there is another possible optical transition, which may have a similar energy to the one between the VB and the CB at the M point. This transition can occur between the upper spin-split VB and the lower spin-split CB+2 (see Figure 13). Our DFT results shown in Table 13 suggest that for MoX 2 the transition at the M point has lower energy, while for WX 2 they are nearly degenerate. This prediction does not take into account excitonic effects, which are also expected to be important and may determine which transition actually has the lower energy, because the exciton binding energies at K and M may be different. We expect that, e.g., the polarisation of the photoluminescence can give important information about these transitions. If the incident light is circularly polarised, then the photoluminescence related to the transition VB→CB+2, which takes place at the K and −K points, should also be circularly polarised, as is the case for the well known VB→CB+2 transition. Since the local symmetry at the M point is different, we do not expect that the photoluminescence due to the VB→CB transition at the M points is circularly polarised.

Conclusions
In this short review we have focused on the band structure of monolayer TMDCs. Our aim has been to discuss all the details of the band structure that we believe are relevant for transport and relaxation processes and optical transitions. The two main tools that we used were the (local) symmetries of the BZ (an essential ingredient of the k · p expansion) and DFT calculations. The first of these tools allowed us to capture general features of the band structure. Material parameters, such as effective masses, spin-splittings, and band edge energy differences depend on the chemical composition of particular TMDCs and are important for quantitative predictions. For this reason we also performed extensive DFT calculations which can, in many cases, predict material parameters accurately. From a theoretical point of view, an important aspect of the approach used in this work is that it leads to explicit k · p Hamiltonians that can be used to address a variety of problems. In particular, they are expected to be accurate when external perturbations vary on length scales much larger than the interatomic distances. Therefore we believe our results will help to develop further (semi)analytical approaches to study, e.g., exciton physics [83,133,134,135,136,137,138], plasmons [139,140], diffusive transport [141], spin [142,145], noise [150], topological properties [149,143], valley-currents [96,146,151], proximity effect [147,148], electron-electron interaction [152], and quantum dots [47,153]. On the other hand, TB-based methods are probably more appropriate for studying the effects of point or line defects [39,154] on transport.
The picture that emerges from this study is that monolayer TMDCs in the ballistic limit should display a remarkable variety of electronic properties, many of which are yet to be verified experimentally. In this section we give a brief account of the seven-band k· p model, which lies behind the effective Hamiltonian of Equation (2). A somewhat longer version of the following discussion has already been published in Reference [47], but we think that it is helpful to present a shorter version here to make this work self-contained.
Our seven-band model (without spin) contains every band between the third band below the VB (which we denote by VB−3) and the second band above the CB (denoted by CB+2), i.e., we take the basis vb, cb, cb + 1, cb + 2} denotes the band, the lower index µ indicates the pertinent irreducible representation of the point group C 3h , which gives the symmetry of the bands at the K point of the BZ. The spinful symmetry basis functions are introduced by |Ψ b µ , s = |Ψ b µ ⊗ |s , where s = {↑, ↓} denotes the spin degree of freedom. An important symmetry of the system is that is has a horizontal mirror plane. As a consequence, the basis states can be grouped into two groups: the first group contains states whose orbital part is symmetric with respect to the mirror operation wherep ± =p x ± ip y are momentum operators. Non-zero matrix elements only exist between states |Ψ b µ , s and |Ψ b µ , s that are either both symmetric or both antisymmetric with respect to the mirror operation σ h . Since Equation (A.1) does not contain spinoperators, the matrix elements are diagonal in the spin-space. The matrix elements H K k·p calculated at the K point of the BZ are shown in Table A1, where the diagonal elements are the band-edge energies. The matrix elements at the −K point can be obtained with the substitutions γ i → γ * i and q ± → −q ∓ . Concrete values for the parameters γ i can be obtained for each material by, e.g., directly evaluating the matrix elements Ψ b µ |p ± |Ψ b µ using Kohn-Sham orbitals. We used this method to calculate the values denoted by γ KS in Table 6.

Appendix A.2. Spin-orbit coupling
In the atomic approximation the SOC is given by the Hamiltonian where V (r) is the spherically symmetric atomic potential,L is the angular-momentum operator, andŜ = (s x , s y , s z ) is a vector of spin Pauli matrices s x , s y , and s z (with eigenvalues ±1). Note thatL·Ŝ =L z s z +L + s − +L − s + , whereL ± =L x ± iL y and s ± = 1 2 (s x ± is y ). The task is then to calculate the matrix elements of Equation (A.2) in the basis introduced earlier in this section.
The non-zero matrix elements H so can be obtained by considering the transformation properties of the basis states and angular-momentum operators with respect to the mirror operation σ h and the rotation by 2π/3 around an axis perpendicular to the plane of the monolayer TMDC. Note, for example, in contrast to the Hamiltonian in Table (A1), the SOC Hamiltonian shown in Table (A2) has non-zero matrix elements between symmetric and antisymmetric basis states as well This is due to the fact that theL ± operators are themselves antisymmetric with respect to σ h . The full SOC Hamiltonian at K is shown in Table ( The change of the wavefunction symmetry notation follows from the assumption that orbital wave functions at K and −K are connected by time-reversal symmetry, i.e., |Ψ b µ (K) =K 0 |Ψ b µ (−K) , whereK 0 denotes complex conjugation.

Low-energy effective Hamiltonian
The low-energy effective Hamiltonian of Equation (4a) can be obtained from H 0 +H k·p + H so by means of Löwdin partitioning (see, e.g., Reference [110]) by considering terms up to third order in various off-diagonal couplings.
First, one can determine |γ|, α s , and β s in the following way. For small enough q, the largest energy scale under the square root in Equation (B.1) is the band gap (for a given spin s) E can be obtained by fitting the CB and VB in the vicinity of the K point with a parabola, Equations (B.4)-(B.5) constitute four equations for five unknown parameters |γ|, α s , and β s . As explained in Section 5.2, the fitting around K was done in a range that corresponds to 5% of the Γ-K distance. The dispersion over this range can be considered to be isotropic and the difference in the effective masses along K-Γ and K-M can be neglected. Over the same range in q, one can also fit the functionε vb +ε cb 2 + c i.e., we have obtained four more equations for |γ|, α s , and β s . Using Equations (B.4)-(B.7) one finds eight equations for the five unknown parameters |γ|, α s , and β s , which can be solved as a linear least-squares problem. The solution, however, depends on the value of the quasiparticle band gap E (s) bg /2 used in the least-squares problem. As shown in Table 5 this is significantly underestimated in DFT calculations. Therefore we have performed the fitting using both the DFT band gap and the GW gap. Note that in order to find E (s) bg one has to add to the E bg values shown in Table 5 the relevant spin-splitting energies, which can be found in Tables 3 and 4. The two approach lead to two sets of parameters. In both cases the same effective masses m (s) cb and m (s) vb , obtained from our DFT calculations, were used. Since the available experimental results suggest that DFT can capture the effective masses quite well, at least in the VB (see Tables 4  and 9), we think that this is a reasonable approach to take into account the results of GW calculations.
Finally the four remaining parameters κ s and η s were determined in following way. Similarly to the previous step, a function of the form ε vb +ε cb 2 + c   4 were required to give the best fit simultaneously to both the CB and the VB. The fitting was performed along the Γ-K-M directions around K and the fitting range corresponded to ≈ 16% of the Γ-K distance. Note, that cos 3φ q = −1 (cos 3φ q = 1) along Γ-K (K-M ). Since |γ|, α s , and β s are already know by this step, κ s and η s is calculated as κ s = c