Vibrational and optical properties of MoS: from monolayer to bulk
Abstract
Molybdenum disulfide, MoS, has recently gained considerable attention as a layered material where neighboring layers are only weakly interacting and can easily slide against each other. Therefore, mechanical exfoliation allows the fabrication of single and multilayers and opens the possibility to generate atomically thin crystals with outstanding properties. In contrast to graphene, it has an optical gap of 1.9 eV. This makes it a prominent candidate for transistor and optoelectronic applications. Singlelayer MoS exhibits remarkably different physical properties compared to bulk MoS due to the absence of interlayer hybridization. For instance, while the band gap of bulk and multilayer MoS is indirect, it becomes direct with decreasing number of layers. In this review, we analyze from a theoretical point of view the electronic, optical, and vibrational properties of singlelayer, fewlayer and bulk MoS. In particular, we focus on the effects of spinâorbit interaction, number of layers, and applied tensile strain on the vibrational and optical properties. We examine the results obtained by different methodologies, mainly ab initio approaches. We also discuss which approximations are suitable for MoS and layered materials. The effect of external strain on the band gap of singlelayer MoS and the crossover from indirect to direct band gap is investigated. We analyze the excitonic effects on the absorption spectra. The main features, such as the double peak at the absorption threshold and the highenergy exciton are presented. Furthermore, we report on the phonon dispersion relations of singlelayer, fewlayer and bulk MoS. Based on the latter, we explain the behavior of the Ramanactive and modes as a function of the number of layers. Finally, we compare theoretical and experimental results of Raman, photoluminescence, and opticalabsorption spectroscopy.
1 Introduction
For many layered materials, it has been established that the fewlayer or monolayer phases have distinct properties with respect to their bulk counterparts. Often these properties are even changing between the mono, bi and, trilayer phases. Within the layers, the atoms are held together by strong covalent bonds while the interlayer bonds are rather weak and mostly due to van der Waals interaction. As a consequence, the layers can easily be separated by mechanical exfoliation and single, quasi twodimensional (2D), and fewlayer systems of various materials can easily be produced.Novoselov2005 () Some examples are graphene, hexagonal boron nitride (BN), semiconducting transition metal dichalcogenides (M = Mo, W, Ta, and X = S, Se, Te), Wilson1969 () the superconducting metal , or the elemental 2D systems silicene, germanene, and phosphorene Jiang2015 ().
Many of these materials have potential for novel technological functionalities. Graphene is the most prominent singlelayer material Novoselov2004 (). It does not only have outstanding physical properties such as high conductivity, flexibility, and hardness Geim2009 (), but it is also a benchmark for fundamental physics. E.g., it displays an anomalous halfinteger Quantum Hall effect due to the quasirelativistic behavior (linear crossing in the bandstructure) of the electronsKatsnelson2006 (); Katsnelson2012 (). The fascinating properties of graphene have paved the way for intense investigations of alternative layered materials.Novoselov2005 ()
Electronics and optical applications often require materials with a sizeable band gap. For instance, the channel material in fieldeffect transistors must have a sufficient band gap to achieve high on/off ratios Lembke2015 (). In this respect, the semiconducting transition metal dichalcogenides (TMDs) can complement or substitute the zeroband gap material grapheneRadisavljevicB.2011 (). Singlelayer is an appealing alternative for optoelectronic applications with an optical gap of 1.81.9 eV, high quantum efficiencyMak2010 (); Splendiani2010 (), an acceptable value for the electron mobilityLembke2012 (); Baugher2013 (), and a low power of dissipationRadisavljevic2011 (). It has potential application in nanoscale transitors RadisavljevicB.2011 (); Zhang2012 (); Lembke2012 (); Lembke2015 (), photodetectors LopezSanchez2013 (); Zhang2015 (), and photovoltaics applications Fontana2013 (); Furchi2014 (); Pospischil2015 (). Other TMDs such as singlelayer also exhibit high photoluminescence yield Gutierrez2013 ().
In this stimulating scenario, TMDs are being intensively investigated. Fabrication techniques such as the mechanical exfoliation Li2012c (); CastellanosGomez2012b () and the liquid exfoliation Coleman2011 () produce single and multilayer crystals with high crystalline quality at low cost. This has increased notably the amount of research groups working in both fundamental and applied aspects of TMDs. Concerning the electrical and optical properties of singlelayer, multilayer and bulk , extensive experimental investigations have been carried out within the last few years. The most important techniques are photoluminescence, optical absorption, and electroluminescence spectroscopy Mak2010 (); Splendiani2010 (); Korn2011 (); Sundaram2013 (). It is widely accepted that singlelayer has a direct band gap that transforms into an indirect gap with increasing number of layers. Similarly, bandgap engineering is possible by applying strain. The application of strain drives a directtoindirect band gap transition in singlelayer Scalise2012 (); Scalise2014 (); He2013 (); Conley2013 (); Dong2014 (); Guzman2014 (). Moreover, suitable hydrostatic pressure reduces the band gap of single and multilayer resulting in a phase transition from semiconductor to metal Nayak2014 (); Nayak2015 (). The group symmetry and the spinorbit interaction in also raise interesting properties. The control of the valley polarization of the photogenerated electronhole pairs paves the way for using in applications related to nextgeneration spin and valleytronics Xiao2012 (); Mak2012 (); Zeng2012 (); Cao2012 (); Kumar2014 (). Further studies dealing with charged exciton complexes (trions)Mak2013 (); Plechinger2015 () or with second harmonic generation have also been published Kumar2013 (); Malard2013 ().
Many challenges remain to be solved in the field of TMDs. The problem of obtaining high hole mobility in singlelayer hinders the realization of pn diodes. A proposed solution is using a monolayer diode, in which the pn junction is created electrostatically by means of two independent gate voltages Pospischil2015 (); Baugher2015 (); Ross2015 (). Another active research field is the design of Van der Waals heterostructures. Assembling atomically thin layers of distinct 2D materials allows to enrich the physical properties geim2013 (). Techniques like chemical vapor deposition and wet chemical approaches are triggering the fabrication of heterostructures Huang2014 (); Gong2014 (). For example, flexible photovoltaic devices of TMDs/graphene layers exhibit quantum efficiency above 30% Britnell2013 (). A photovoltaic effect has also been achieved using a pn heterojunction Furchi2014 (). The different stacking configurations and the band alignments are important aspects in bilayer heterostructureshe:prb:89 (); debbichi:prb:89 (); liu:natcommun:5 ().
Another important activity in the field of TMDs is the characterization of the vibrational properties of . Earlier studies of bulk using Raman and infrared spectroscopy Verble1970 (); Wieting1971 () as well as Neutron scatteringWakabayashi1975 () and electronenergyloss spectroscopyBertrand1991 () had already well characterized the phonons at and the phonon dispersion. In the recent years, a large number of Raman studies on mono and fewlayer systems has emerged Lee2010 (); Windom2011 (); Li2012 (); Plechinger2012 (); Boukhicha2013 (); Luo2013 (); Zhao2013 (); Zhang2013 (); ZhangX2015 (). The Raman frequencies are correlated with the number of layers which allows their unequivocal identification. The trend of the Raman modes (inplane mode) and (outofplane mode) with the number of layers has been intensively discussed, both theoretically Ataca2011 (); MolinaSanchez2011 (); Luo2013 (); Terrones2014 () and experimentally Lee2010 (); Rice2013 (); Terrones2014 (). The mode follows a predictable behavior. Its frequency grows with increasing number of layers, due to the interlayer interaction. The mode shows the opposite trend, i. e., decreasing in frequency for an increasing number of layers.
The experimental findings are accompanied by a vast theoretical literature. The characteristic stacking of ultrathin layers of adds new challenges to the theoretical approaches. For example, the layer thickness influences the dielectric constant which becomes strongly anisotropic. This enhances the Coulomb interaction between carriers. The calculation of the excitations has to include these dimensional effect for applying accurately the GW method and the BetheSalpeter equation. For instance, one consequence is an exciton binding energy in some layered materials of hundreds of meV, much larger than in bulk semiconductors. Also, new models have been developed to explain the interplay of the spinorbit interaction and the crystal symmetry (which is layer dependent), and its consequences, like the valleyHall effect. Moreover, the interlayer interaction has demanded the improvement of the modelling of the van der Waals interaction in extended systems. The precision of the Raman spectroscopy has allowed to evaluate the accuracy of ab initio methods for calculating phonons, and it proves how useful is the simulation of the vibrational properties for understanding the interlayer interaction and the chemical bonding. Therefore, the research on layered materials has contributed to the appeareance of new methods and to the reformulation of existing ones. In this review, we give an overview of the challenges in the modelling of the spectroscopic properties of and the solutions proposed. The discussion of the literature results is complemented by additional calculations. We believe the topics discussed here will be also useful in the modelling and understanding of other twodimensional materials.
2 Structural properties
Bulk molybdenum disulfide () belongs to the class of transition metal dichalcogenides (TMDs) that crystallize in the characteristic 2H polytype. The corresponding Bravais lattice is hexagonal and the space group of the crystal is ( nonsymmorphic group). The unit cell is characterized by the lattice parameters (inplane lattice constant) and (outofplane lattice constant). The basis vectors are
(1) 
The unit cell contains 6 atoms, two Mo atoms are located at the Wyckoff sites and four S atoms at the Wyckoff sites. With the internal parameter , the positions, expressed in fractional coordinates, are (1/3, 2/3,1/4) for the Mo atoms and (2/3, 1/3,) as well as (2/3, 1/3,1/2) for the S atoms.
The singlelayer contains one Mo and two S atoms. In this case, the inversion symmetry is broken and the space group (more precisely, layer group) is ( symmorphic group). The doublelayer is constructed by adding another SMoS layer, having now the layer group ( symmorphic group). Consequently, an odd number of layers has the same symmetry as the singlelayer (absence of inversion symmetry), whereas an even number has the symmetry of a doublelayer (with inversion symmetry).
(Å)  (Å)  (GPa)  

Ref. dickinson:jacs:45 ()  3.160  12.294  0.621  3.890  
Ref. schoenfeld:acb:39 ()  3.161  12.295  0.627(5)  3.890  
Ref. alhilli:jcg:15 ()  3.140  12.327  3.926  53.41.0  
Ref. petkov:prb:65 ()  3.168(1)  12.322(1)  0.625  3.890 
The crystal structure of can be specified as a stacking of quasitwodimensional (2D) SMoS layers along the direction. Within each layer, Mo atoms are surrounded by 6 S atoms in a trigonal prismatic geometry as illustrated in Fig. 1. The bonding type is predominantly covalent within the atomically thin SMoS layers, whereas the layers themselves are weakly bound by Van der Waals (VdW) forces in the crystal. The inherent weakness of the interlayer interactions can result in different stacking sequences and therefore in different polytypisms as shown in Ref. he:submitted ().
Defining the optimized geometry is the first step for any calculation of the phonon spectra and/or the band structure. Most of the previous investigations used densityfunctional theory (DFT) on the level of the localdensity approximation (LDA) or the generalizedgradient approximation (GGA) Ataca2011 (). We want to emphasize that in DFT the accuracy of the calculated quantities is determined by the treatment of the exchange correlation (XC) energy given by the XC functional. However, the standard local (LDA) and semilocal (GGA) XC functionals do not account for the longrange van der Waals interactions, which are responsible for the stable stacking of the layers and thus particularly relevant in twodimensional materials. Nevertheless, the wellknown LDA overestimates the (weak) covalent part of the interlayer bonding and compensates thus the missing vdW forces yielding a bound ground state for most layered materials. This explains the success of LDA in obtaining the geometry of many layered materials such as graphite Wirtz2004 (), boron nitride Kern1999 (); Serrano2007 () or graphene on different substrates Allard2010 (); Fromm2013 (); Endlich2013 (). The good performance of LDA in layered materials (although fortuitous) has made this approximation widely used in the calculations of structural properties.
In the present work, the calculations were partly performed with the Vienna ab initio simulation package (VASP)kresse:cms:6 (); kresse:prb:54 () utilizing the projector augmented plane wave (PAW) method bloechl:prb:50 (); kresse:prb:59 () to describe the corevalence interaction. PAW potentials with nonlocal projectors for the molybdenum (Mo) 4, 4, 4, 5 as well as sulfur (S) 3, and 3 valence states were generated to minimize errors arising from the frozen core approximation. The valence electrons were treated by a scalarrelativistic Hamiltonian and spin orbit coupling (SOC) was selfconsistently included in all VASP calculations as described elsewhere kim:prb:80 (). VASP uses DFT with a variety of XC functionals ranging from LDA to different types of GGAs, to hybrid functionals, and VdW density functionals. Furthermore, VASP has an implementation of many body perturbation theory such as the approximation ranging from the singleshot shishkin:prb:74 () to a selfconsistent (sc) approximationshishkin:prb:75 (); shishkin:prl:99 (). Concerning standard DFT results presented in this work, the XC energy was treated within the LDALDA:prl:45 () and the GGA. For the latter, the parametrization of Perdew, Burke, and Ernzerhof (PBE), in particular the PBEsol functional PBEsol:prl:100 () was used.
In order to improve the theoretical lattice parameters calculated within DFTLDA/GGAAtaca2011 (), we have also studied the structural properties including VdW interactions starting from the experimentally observed structural parameters summarized in Tab. LABEL:tab:lattice. For this purpose we used the optB86bVdW functional, recently implemented in VASP klimes:prb:83 (). The optB86b functional is a nonlocal correlation functional that approximately accounts for dispersion interactions. It is based on the VdWDF proposed by Dion et al.dion:prl:92 (), but employs an accurate exchange functional particularly optimized for the correlation part. klimes:prb:83 ()
The structural optimization of hexagonal bulk requires a 2D energy minimization, since the ground state energy depends on two degrees of freedom, i. e., the volume and the ratio. The experimentally observed structural parameters summarized in Tab. LABEL:tab:lattice have been used as starting point for the calculation of the electronic properties of bulk as well as single layer (1L) within DFT and methods that go beyond as described in detail below. This minimization was performed manually using LDA and PBEsol as well as the optB86bVdW functional, for comparison by varying the unit cell volume of bulk within of the experimentally observed equilibrium volume (). For each chosen volume the ratio was first optimized by fitting the energy dependence on to the Murnaghan equation of state (EOS). murnaghan:pnas:30 () The final DFToptimized unit cell volume was obtained by subsequently fitting to the Murnaghan EOS. Note that in each single optimization step the atomic positions were also relaxed by minimizing the total forces on the atoms until they were converged to 0.05 eV/Å.
In order to avoid effects from the changes in size of the basis set due to changes in the unit cell volume , the kinetic energy cutoff has been increased to 350 eV. Convergence with respect to sampling within the Brillouin zone was reached with centered meshes in case of optB86bVdW and with centered meshes for LDA and PBEsol. The manually performed structural optimization was cross checked with VASP calculations employing minimization algorithms parallel for the atomic positions and the ratio for selected volumes in the range of and one subsequent Murnaghan EOS fit. From these calculations the Bulk modulus is obtained from
(2) 
In Table 2 the results of the structure optimization corresponding to different functionals are summarized.
(Å)  (Å)  (GPa)  (Å)  VdW gap (Å)  

LDA  3.120  12.09  0.1214  3.87  4043  6.039  2.933 
PBEsol  3.138  12.60  0.1264  4.01  1821  6.305  3.188 
optB86bVdW  3.164  12.40  0.1236  3.92  3940  6.203  3.068 
When using LDA or optB86bVdW functionals, the theoretical values of for bulk agree well with the experimental value given in Table LABEL:tab:lattice. Concerning lattice parameters, we observe a small underestimation of the inplane parameter , both in LDA (1.3 %) and PBEsol (0.7 %). Contrary, the parameter is underestimated by 1.6% and overestimated by 2.4% in LDA and PBEsol, respectively. Compared to LDA, the PBEsol functional yields the larger deviation of the resulting ratio. The improvement after including the VdW interactions is substantial, with an error of only 0.09 % (see Fig. 2 for a comparison with the experimental results).
In conclusion, van der Waals functionals give the most accurate results for lattice parameters and the bulk modulus. LDA tends to underestimate the interlayer distance and the parameter, but in average gives acceptable results and it can be trusted in the prediction of structural properties.
Based on the ground state structures summarized in the bulk charge density was calculated on a centered 12123 mesh by converging the total energy to 0.1 meV using a kinetic cutoff energy of 350 eV and a Gaussian smearing with a smearing width of 50 meV. Tests with 18183 point grids have shown that the electronic band gaps are converged within 20 meV compared to the results obtained with the 12123 grid.
The singlelayer structure has been constructed from the optimized bulk structure (Tab. 2) by selecting only the bottom SMoS layer and adding 20 Å vacuum along direction. The atomic positions in the slab geometry have again been relaxed (force convergence criterion of 0.05 eV/Å) before calculating the band structure on centered 12121 point grids. Convergence tests of the eigenvalues as a function of the vacuum space between repeated layers has been performed up to an accuracy of 15 meV, establishing a convergence distance of 20 Å.
As we will demonstrate later, the main features of the band structure of critically depend on the lattice optimization and even small differences can induce significant changes. In the case of the phonon band structure, deviations are reflected in a rigid shift of the phonon frequencies but the main trends are less affected than in the case of the electronic band structure.
3 Lattice dynamics of MoS
The knowledge of the phonon dispersion in a material is indispensable for the understanding of a large number of macroscopic properties such as the heat capacity, thermal conductivity, (phononlimited) electric conductivity, etc. Vibrational spectroscopy (Raman spectroscopy and Infrared absorption spectroscopy)Wieting1971 () give access to the phonons at the Brillouin zone center ( point). Inelastic neutron scattering citeWakabayashi1975 allows to measure (almost) the full phonon dispersion. Precise semiempirical modeling of the phonon dispersion and abinitio calculations in comparison to experimental data are a challenge by itself. However, precise modeling is also required because details in the vibrational spectra may also carry some information about the number of layers and the underlying substrate. For graphene, this has been widely explored: the socalled 2D line in the spectra splits into subpeaks when going from the single to the multilayer caseFerrari2006 (); Graf2007 (). Last but not least, the 2Dline also changes position as a function of the underlying substrateBerciaud2009 (); Forster2013 (); Starodub2011 (); Endlich2013 (). All these features are related to the doubleresonant natureThomsen2000 () of Raman scattering in graphene and on the dependence of the highest optical mode on the screening. For MoS (and related semiconducting transitionmetal dichalcogenides), the layerdependence of the vibrational spectra is less spectacular than for graphene. Nevertheless a clear trend in the lower frequency interlayer shear and breathing modes with increasing layer number can be observedPlechinger2012 (); Zhang2013 (); Boukhicha2013 () (similar to graphene Tan2012 ()). Also in the highfrequency optical modes at , a clear trend from single to multilayer MoS has been observedLee2010 () and reproduced in other dichalcogenidesBerkdemir2013 ().
Interestingly, the behaviour of the phonon frequencies as a function of the number of layers does not always follow the intuitive trend. For instance, the frequency of the Raman active inplane phonon decreases with the increment of the number of layers. This is contrary to the expectation that the weak interlayer forces should increase the restoring forces and consequently result in an increase of the frequency. The Raman active does follow the expected behaviour, increasing the frequency with the number of layers. Several attempts have been done to explain this trend Lee2010 (); MolinaSanchez2011 (); Luo2013 (). This will be critically reviewed in this section.
General Theory
Before discussing the phonons of bulk and fewlayer MoS, we briefly review the abinitio calculation of phonons. (A complete discussion of the theory of phonons can be found, e.g., in Ref. Bruesch1982 ().) Starting from the equilibrium geometry of the system (see details in Section 2), one obtains the phonon frequencies from the solution of the secular equation
(3) 
where is the phonon wavevector, and and are the atomic masses of atoms and . The dynamical matrix is defined as
(4) 
where denotes the displacement of atom in direction . The second derivative of the energy in Eq. 4 corresponds to the change of the force acting on atom in direction with respect to a displacement of atom in direction Bruesch1982 (). The elements of the dynamical matrix at a given wavevector can be obtained from an abinitio total energy calculation with displaced atoms in a correspondingly chosen supercell (that needs to be commensurate with ). Another approach consists in the use of density functional perturbation theory (DFPT)Baroni2001 (); Gonze1997 (), where atomic displacements are taken as a perturbation potential and the resulting changes in electron density and energy are calculated selfconsistently through a system of KohnSham like equations. Within this approach the phonon frequency can be obtained for any while using the primitive unitcell. Another way to obtain the dynamical matrix is to use empirical interatomic potentials or force constants. A decent fit of the MoS phonon dispersion from StillingerWeber potentials has recently been suggested by Jiang et al.Jiang2013 (). Such a fit allows to study much larger systems such as nanoribbons, nanotubes or heterostructures. The advantage of DFPT is, however, the higher accuracy and the automatic inclusion of mid and longrange interactions. For later reference, we note that realspace force constants can be obtained from the dynamical matrix by discrete Fourier Transform:
(5) 
where is the number of unitcells (related to the density of the point sampling). The physical meaning of the is the force acting on atom in direction as atom in a unit cell at is displaced in direction and all other atoms in the crystal are kept constant.
It is important to note that for polar systems, the phonons in the longwavelength limit can couple to macroscopic electric fields. In mathematical terms, this means that the dynamical matrix in the limit can be written as the sum of the dynamical matrix at zero external field and a “nonanalytic” part that takes into account the coupling to the electric field and depends on the direction in which the limit is taken:
(6) 
The nonanalytic part contains the effect of the longrange Coulomb forces and is responsible for the splitting of some of the longitudinal optical (LO) and transverse optical (TO) modes. Baroni2001 (); Gonze1997 (): Its general form is as follows:
(7) 
where is the volume of the unit cell, stands for the Born effective charge tensor of atom and is the dielectric tensor. Since the dielectric tensor is fairly large in bulk MoS (, the effect of LOTO splitting is visible, but not very pronounced ( 2.6 cm).
We will discuss in the following first the phonon dispersions of bulk and singlelayer MoS. Afterwards, we will discuss in detail the symmetry of bulk and singlelayer phonons at . We will single our the Raman and infrared (IR) active modes and discuss how their frequencies evolve with the number of layers.
Phonon dispersion
The phonon dispersions of singlelayer and bulk MoS are shown in Figure 3. Overall, the singlelayer and bulk phonon dispersions have a remarkable resemblance. In the bulk, all singlelayer modes are split into two branches but since the interlayer interaction is weak, the splitting is very low (similar to the situation in graphite and graphene.Wirtz2004 () The only notable exception from this is the splitting of the acoustic modes of bulk MoS around .
We have also depicted the experimental data obtained with neutron inelastic scattering spectroscopy for bulk Wakabayashi1975 () as well as the result of IR absorption and Raman scattering at . The overall agreement between theory and experiment is rather good, even for the interlayer modes. This confirms our expectation that the LDA describes reasonably well the interlayer interaction (even though not describing the proper physics of the interlayer forces, as discussed in Section 2).
The bulk phonon dispersion has three acoustic modes. Those that vibrate inplane (longitudinal acoustic, LA, and transverse acoustic, TA) have a linear dispersion and higher energy than the outofplane acoustic (ZA) mode. The latter displays a dependence analogously to that of the ZA mode in graphene (which is a consequence of the pointgroup symmetry Saito1998 ()). The low frequency optical modes are found at 35.2 and 57.7 cm and correspond to rigidlayer shear mode and layerbreathing mode (LBM) respectively (see left panels of Fig. 5) in analogy with the low frequency optical modes in graphiteWirtz2004 ()). It is worth to mention the absence of degeneracies at the high symmetry points and and the two crossings of the LA and TA branches just before and after the point. The high frequency optical modes are separated from the low frequency modes by a gap of 49 cm.
The singlelayer phonon dispersion is very similar to the bulk one. The number of phonon branches is reduced to nine. At low frequencies, the shear mode and layerbreathing mode are absent. At higher energies, very little difference between bulk and singlelayer dispersion can be seen. This is due to the fact that the interlayer interaction is very weak. The subtle splitting and frequency shifts of zonercenter modes in gerade and ungerade modes (as going from single layer to the bulk) will be discussed below.
The densities of states (DOS) of singlelayer and bulk are represented in the right panels of Fig. 3. The differences between singlelayer and bulk DOS are minimal, except a little shoulder around 60 cm in the bulk DOS due to the low frequency optical modes. In both cases the highest peaks are located close to the Raman active modes and and they are due to the flatness of the bands around .
The density of states can be partially measured in 2nd and higherorder Raman spectra. We have represented in Fig. 4 the Raman spectrum of MoS bulk of Ref. Windom2011 (), obtained by exciting the sample with a laser frequency of 632 nm. Similar spectra can be found in older studies Chen1974 (); Stacy1985 () and other recent studies Livneh2010 (); Livneh2014 (). First order Raman peaks are due to the excitation of zeromomentum phonons. We can identify in the spectrum the modes and . Moreover, working in conditions of resonant Raman scattering, secondorder Raman modes can be obtained. These modes come from the addition or subtraction of modes with opposite momentum, , together with the resonance of an intermediate excited electronic stateBerkdemir2013 (). The result is a rich combination of Raman modes, as shown in Fig. 4. The identification can be done with the help of the density of states. We have calculated the 1phonon and the 2phonon density of states for MoS bulk, represented by the solid blue line and the dashed green line, respectively.
A careful examination of the 2phonon density of states tells us which phonons are participating in the Raman spectra. Thus, the longitudinal acoustic branch at M, denoted as , couples to the modes , , and , resulting in overtones in the Raman spectrum. Other combinations include or , always with momentum . The secondorder Raman modes are much more restrictive than firstorder modes. The concurrence of phonon modes and electronic levels is needed. Such alignment depends strongly on the electronic structure. Consequently, the secondorder Raman spectrum has revealed useful to establish the fingerprints of singlelayer systems with respect to the bulk, as discussed in Ref. Berkdemir2013 ().
Symmetry analysis of phonon modes
We have drawn in Figs. 5 and 6 the atomic displacements (eigenvectors) of optical phonon modes of bulk and singlelayer MoS at . Group theoretical analysisVerble1970 (); Wieting1971 (); MolinaSanchez2011 (); Luo2013 (); RibeiroSoares2014 () yields for the 15 optical modes of bulk MoS (D symmetry) the following decomposition in irreducible representations: . The , and modes are Raman active and the and modes are infrared active. For the 15 optical modes of doublelayer MoS (D symmetry), the decomposition is: where the gerade modes are Raman active and the ungerade ones are IR active. For the 6 optical modes of the singlelayer, one obtains the following irreducible representations: . The and modes are Raman active, the mode is both IR and Raman active.
The attribution of the different symmetries to the phonon modes in Figs. 5 and 6 can be understood quite intuitively (taking into account that the drawings are simplified 2D versions of the real 3D modes). All modes are doubly degenerate and correspond thus to inplane vibrations of Mo and/or S atoms because a rotation by any angle yields another phononmode with the same frequency. The nondegenerate and modes must therefore correspond to perpendicular movement of the atoms. For bulk MoS, (space group , point group ), there is an inversion center halfway between two S atoms of neighboring layers. We can thus distinguish between gerade (g) and ungerade (u) modes which are symmetric/antisymmetric with respect to inversion. The gerade modes are those where atoms in neighboring layers move with a phase shift of while the ungerade modes correspond to inphase movement. All phonon modes of bulk MoS thus come in pairs of gerade and ungerade modes which are close in frequency. This can be clearly seen for the modes in panels 1 to 4 of Fig. 5. Furthermore, the “shear mode” at 35.2 cm is the gerade counterpart of the inplane acoustic mode (not shown) and the “layerbreathing mode” (LBM) at 55.7 cm is the gerade counterpart of the outofplane acoustic mode. In almost all cases, the gerade mode is higher in frequency than the ungerade mode. This is because the weak (VanderWaals like) bond between S atoms of neighboring sheets is elongated and squeezed in the gerade mode (thus gives rise to an additional restoring force) but kept constant in the ungerade mode. The notable exception is the case of the modes in panel 2 where the ungerade mode is higher in energy. We will come back to this important case in the next subsection. One can easily see that only ungerade modes can be IR active: for a mode to be IR active, a net dipole must be formed through the displacement of positive charges in one direction and negative charges in the opposite direction. However, in gerade modes, the dipoles formed on one layer are canceled out by the oppositely oriented dipoles on the neighboring layer. Since in systems with inversion symmetry, a phonon mode cannot be both IR and Raman active, only the gerade modes can be Raman active in bulk MoS.
The distinction between and modes is made by rotating the crystal by around the principal rotation axis. This rotation is a nonsymmorphic symmetry, i. e., it has to be accompanied by a translation normal to the layerplane in order to map the crystal into itself. In our reduced 2D representation of the vibrational modes this corresponds to a translation of the 3 atoms of the upper layer onto the 3 atoms of the lower layer. If the arrows change direction, the mode is , otherwise . Finally, for the singly degenerate modes, the subscript () stand for modes that are symmetric (antisymmetric) with respect to rotation around a C axis crossing an Mo atom perpendicularly to the 2D plane of projection. For the doubly degenerate modes, it is the other way around.
For even Nlayers of MoS, the spacegroup symmetry is and the assignment of the phononmode symmetries has to be done according to the pointgroup symmetry. Since inversion symmetry is present, the mode assignment is very similar to the one of bulk MoS. For the doubly degenerate modes (see Fig. 5), the subscripts 1 and 2 are dropped. All modes are IR active and all modes are Raman active. Out of the perpendicularly polarized modes, the inactive mode turns into an IR active mode, the inactive modes turns into a Raman active modes. Notably, the layered breathing mode (LBM) is, in principle, Raman active. Indeed, for double and 4layer MoS, this mode has been detected in Raman measurements, albeit with small amplitudeZhang2013 ().
For the single layer and for addnumbered multilayers, the space group is and the corresponding pointsymmetry group is D. Since inversion symmetry is absent in this group, there is no distinction between gerade and ungerade modes. Instead, modes that are symmetric under (reflection at the xyplane) are labeled with a prime and antisymmetric modes with a double prime (Fig. 6).
The experimental and theoretical frequencies of all phonon modes of singlelayer and bulk MoS at are summarized in Table 3. For the IR active modes of bulk MoS, we give both the values for longitudinaloptical (LO) and transverseoptical (TO) modes. The LOTO splitting is calculated from the nonanalytic part of the dynamical matrix (Eq. 7) which only affects the IR active modes.
single layer ( symmetry)  bulk ( symmetry)  
Pol.  Atoms  Irrep.  Char.  Freq. (cm)  Irrep.  Char.  Freq. (cm)  
Calc.MolinaSanchez2011 ()  Exp.Lee2010 ()  Calc.MolinaSanchez2011 ()  Exp.Wieting1971 ()  
Mo+S  a  0.0  a  0.0  
R  35.2  33  
Mo+S  a  0.0  a  0.0  
i  55.7  
S  R  289.2  i  287.1  
R  288.7  287  
Mo+S  R+IR()  391.7  384.3  R  387.8  383  
IR()  TO: 388.3  384  
LO: 391.2  387  
S  R  410.3  403.1  i  407.8  
R  412.0  409  
Mo+S  IR()  476.0  IR()  TO: 469.4  470  
LO: 472.2  472  
i  473.2 
Anomalous Davydov splitting
As mentioned above, in bulk MoS, all modes appear as pairs of gerade and ungerade modes (Fig. 5 and Table 3). This is because the unit cell of bulk MoS contains 6 atoms while the singlelayer unit cell only contains 3. The frequency splitting of gerade and ungerade modes in layered materials and molecular crystals is known as “Davydov splitting” or “factor group splitting” Davydov1969 (); Dawson1975 (). The model of Davydov has been used in particular to explain the splitting of modes that are IR and R active for the singlelayer into a Raman active mode and an IR active mode of the bulk (mode No. 2 in Figs. 5 and 6). It was recognized early Wieting1971 (); Kuroda1979 () that neither the weak VanderWaals coupling between neighboring layers nor a simple model of dipolar couplings matches the experimental observation that for some layered materials and, in particular, MoS. A model involving quadrupole interaction was proposed by Ghosh et al.Ghosh1976 (); Ghosh1983 () but could not be underpinned by numerical calculations.
The explanation of the “normal” Davydov splitting in vanderWaals bonded layered materials is straightforward. As can be seen in Fig. 5, the weak interlayer bonding can be viewed as an additional (weak) spring constant acting between sulfur atoms from neighboring layers (red dashed lines). For the ungerade modes, the Satoms are moving in phase and the additional spring thus is not “used”. However, for the gerade modes, where the Satoms are vibrating with a phase shift of , the additional spring is elongated and compressed and thus yields an additional restoring force. This leads, in general, to an upshift of the frequencies of the gerade modes. Since the interaction is weak, the frequency shift is small ( 5cm). Furthermore, the effect is more pronounced for the perpendicularly polarized modes than for the inplane modes (Fig. 5 and Table 3). The only exception from the “normal” Davydov splitting is thus the mode No. 2. One might argue that this case is exceptional, because the LOshift of the mode makes its frequency higher then the one of the mode. However, even without the LOshift, experimentsWieting1971 () and calculationsMolinaSanchez2011 () agree that .
Since our abinitio calculations reproduce the anomalous Davydov splitting for mode No. 2, we can use the interatomic force constants, derived from the calculations, in order to find the origin of this seemingly anomalous behaviour. The situation is depicted in Fig. 7 (a) and (b). We analyze in detail the force constants between S atoms of neighboring layers (blue springs) and also the force constants between S atoms on one layer with Mo atoms on the neighboring layer (red springs). In the mode it is the sum of all the S–S spring constants that leads to an additional restoring force and thus to an upshift (with respect to the same mode in the fictitious isolated layer). However, for the mode, an additional restoring force arises as well, this time due to the Mo–S interactions. As it turns out, their effect is stronger than the ones of the S–S springs. This follows from the numerical values of the horizontal components of the S–S and Mo–S force constants. We present the values in Fig. 7 (c). For a given S atom, we calculate the sum of the (horizontal) force constants over all nearest, nextnearest, …, S and Mo atoms of the adjacent layer. Negative sign implies restoring force (the S atoms is pushed back to the left when displaced to the right). In Fig. 7 (c) one can see that the interaction of S with the three nextnearest S atoms of the adjacent layer is stronger than the interaction of S with the closest Mo atom. However, the cumulative effect of the S–Mo interactions is larger than the one of the S–S interaction. This explains why for this mode pair the sign of the Davydov splitting is negative. The dominance of the interlayer S–Mo interaction over the interlayer S–S interaction was already invoked in the forceconstant model of Luo et al.Luo2013 () In that model, all the interaction was ascribed to the closest atom pairs of adjacent layers. This renders the semiempirical model simple and quantitatively successful. In semiempirical calculations using the code GULPGale1997 (); Gale2003 (), we have verified that it is possible to inverse the frequencies of the and modes by changing the relative strength of the crosslayer S–S and S–Mo interaction (see Fig. 7 (d)). However, the abinitio calculations demonstrate that the physical reality is more complex. The 2ndnearest neighbor interaction between sulfur atoms across layers is even repulsive. Thus, the correct balance between S–S and S–Mo interaction is only found by summing over all interactions. It turns out that the force constants decrease rather quickly with increasing distance and the Coulomb contribution (from the effective charges) is rather small.
Note that the situation is different for mode No. 1 in which the Mo atoms do not move and the interlayer S–Mo interaction thus plays a negligible role for the Davydov splitting. In this case, the Raman active mode is slightly higher in frequency than the inactive mode, following the intuitive expectation and yielding the “normal” sign of the Davydov splitting.
We will see in the following that also the dependence of the frequencies of the on the number of layers follows an unexpected trend which can be used for the determination of the number of layers via Raman spectroscopy.
Dependence of Raman active modes on number of layers
Since the beginning of the research on flakes, the Raman modes have been used to identify the number of layersLee2010 (); Plechinger2012 (); Zhao2013 (); Zhang2013 (). The correspondence between frequency and number of layers has been done by comparing with other techniques such as atomic force microscopy or optical contrast. The phonon modes used to characterize the thickness are usually the highfrequency Raman modes and (see Fig. 5) and the breathing and shearing modes at lowfrequency Zhang2013 (). We will summarize in the following the results and analyze the physics of the frequency trends.
Highfrequency modes
In the single layer, the high frequency modes and collapse into the mode . (From Fig. 6 it is evident that with increasing interlayer distance, the modes and acquire the same frequency.) Interestingly, as measured in Ref. Lee2010 () and indicated in Table 3 (see also Figs. 5 and 6), the bulk mode is lower in frequency than the singlelayer mode. This contradicts the expectation that the additional interlayer interaction should increase the frequency from singlelayer to bulk. Even the bulk mode (which is higher in frequency than the mode due to the anomalous Davydov splitting) is slightly lower than the single layer mode. The same behaviour (that the bulk modes are lower in frequency than the singlelayer mode) can be observed for the inplane mode No. 1 ( and in bulk versus in the single layer) and for the outofplane mode No. 4 ( and in bulk versus in the single layer). Only the outofplane mode (No. 3) follows the expected trend that the interlayer interaction increases the frequency with respect to the singlelayer mode .
Figure 8 shows the frequency of the and modes as a function of layer number. We compare LDA calculations (circles comes from Ref. MolinaSanchez2011 () and triangles from Ref.Luo2013 ()) with the experimental data of Refs. Lee2010 () and Luo2013 (). The outofplane mode increases in frequency with increasing number of layers. The explanation lies in the interlayer interaction, caused by weak vanderWaals bonding of the sulphur atoms of neighbouring layers. The stacking of layers can thus be seen as the addition of a spring between the sulfur atoms of neighboring layers. Within the picture of connected harmonic oscillators, this results in an increase of the frequency. The LDA calculations reproduces well this trend, as can be seen in Fig. 8. The small disagreement between theory and experiment can be attributed to the inaccuracy of the interlayer interaction given by LDA.
The inplane mode displays the opposite trend, decreasing in frequency by about 4 cm from singlelayer to bulk (Fig. 8, lower panel). This is  at first sight  unexpected, because the additional “spring” between the sulfur atoms should lead to an increased restoring force and thus to a frequency increase as in the case of the mode. Several attempts have been made in the past to explain this anomalous behaviour, ascribing it to longrange Coulomb interactionsLee2010 (). In our previous previous publicationMolinaSanchez2011 (), we have investigated how the dielectric screening in the bulk environment reduces the longrange (Coulomb) part of the force constants. However, the longrange part plays only a minor role. We have verified this by performing an abinitio phonon calculation of the mode of singlelayer MoS sandwiched between graphenelayers. If the distance between the sulfur atoms and the graphene layer is higher than 6Å, there is no “chemical” interaction between the different layers and the graphene just enhances the dielectric screening of the MoS layer. Since the remains unaffected, we conclude that the longrange Coulomb effect can be discarded as a possible effect for the anomalous frequency trend.
The solution to the problem has been given by Luo et al.Luo2013 () and is related to a weakening of the nearest neighbor Mo–S forceconstant in the bulk environment. To be precise, one has to compare the Real space force constants for the force and displacement parallel to the layer. (See Eq. (5) for the definition of the force constant). This is the dominant force constant determining the frequency of the mode as becomes immediately clear from the displacement pattern in Fig. 6. Force constants between more distant atoms play a very small role. There are two reasons why is larger for the singlelayer than for bulk. One reason is that in the singlelayer, the Mo–S bond is slightly shortened. However, even without this change in bondlength, the Mo–S force constant is slightly larger in the singlelayer than in bulk. This can be obtained in an abinitio calculation of the singlelayer using the interatomic distances from bulk. The results of our calculations are shown in Table 4 and are very similar to the values of Fig. 3 in Ref. Luo2013 (). The frequency of the mode is proportional to the forceconstant: . Even though the differences seem tiny, they explain the calculated and observed frequency differences between singlelayer and bulk in table 3. The finding that the Mo–S force constant is smaller in bulk than in singlelayer, even at equal bondlength, is related to a (slight) redistribution of the charge density when a neighboring layer is present.
bulk  SL  

bulk geom.  opt. geom.  
0.1102  0.1119  0.1127 
The fact that the and the modes move in opposite directions with increasing number of layers, makes the distance between the two corresponding peaks in the Raman spectra a reliable measure for the layer numberLee2010 (); Luo2013 (). But this is not the only way to detect the layer number in Raman spectroscopy. The low frequency Raman active modes display an even stronger dependence as explained in the following.
The shearing mode (C), denoted in bulk as , is the rigidlayer displacement inplane. This mode is Raman active in bulk, as indicated in Table 3. The layerbreathing mode (LBM) corresponds to vertical rigidlayer vibrations, in the case of bulk, where is has symmetry, it is a silent mode. However, in the bilayer case it has symmetry and is (weakly) visible. Several groups have investigated the low frequency behaviour of fewlayer MoS Plechinger2012 (); Zhang2013 (); Zhao2013 (). The frequency trends as a function of the layer number can be explained via a simple analytical model that was first developed to explain the corresponding modes in fewlayer grapheneTan2012 (); Michel2012 (). In this model, layers with a mass per unit area are connected via harmonic springs. One distinguishes between force constants (per unit area) for displacement perpendicular and parallel to the layer, respectively. Mathematically, the model is equivalent to a linear chain of atoms. Assuming a time dependence Assuming a time dependence of for all the atoms, Newton’s equation of motion lead to the secular equation
(8) 
where for the layerbreathing modes and for the shear modes. The frequency of the th phonon mode is
(9) 
For one obtains the acoustic (translational) mode, yields the lowest nonvanishing frequency of the mode with a nodeless envelope function, and yields the highest frequency mode where neighboring layers are vibrating with a phase shift of (almost) .
Fig. 9 shows the Raman spectra published in Ref. Zhang2013 (). The number of layers ranges from 1 to in the case of odd number of layers (ONL) and from 2 to in the case of even number of layers (ENL). Evidently, the singlelayer Raman spectra has no signature of lowfrequency modes. The peaks at cm are due to the Brillouin scattering of the LA mode from the silicon substrate. One can see that some peaks stiffen (dashed lines) for increasing thickness and others are softened (dotted lines). Fig. 9(c) shows the shear (C) and breathing (LBM) mode as a function of the number of layers. We can see also in Fig. 9(d) that the full width at half maximum (FWHM) behaves in a different way for the C or the LBM modes. In the case of the LBM mode (blue dots in in Fig. 9 (a,b)), the branch with the largest weight is the one with . According to Eq. 9, the frequency of this branch as a function of layer number goes like
(10) 
where . This is the functional form of the blue diamonds in Fig. 9 (c). For , the LBM increases in frequency according to
(11) 
and approaches, for , the value of the bulk mode at 55.7 cm. However, since the bulk mode is not Raman active, the intensity of this mode quickly decreases with increasing and the mode is already almost invisible for . For intermediate values of , side branches of the LBM appear and are clearly visible in the Raman spectra of Fig. 9.
The same analysis can be done for the shear (C) mode. In this case, it is the branch that has the dominant weight and converges towards the bulk shear mode with increasing :
(12) 
This is the functional form of of the red squares in Fig. 9 (c). The frequency ratio of bulk and double layer shear modes is which is verified by the experimental dataPlechinger2012 (); Zhang2013 (). Some sidebranches for are visible in the spectra as well, however with lower intensity than the branch.
Due to the strong layer dependence of the frequencies, the shear and compression mode are a very sensitive tool for the determination of layerthickness by Raman spectroscopy Plechinger2012 (); Zhang2013 (). The monoatomic chain model is able to explain the main physics of these modes. Small deviations from the analytic result have been observed Plechinger2012 () and might be due to anharmonic effectsBoukhicha2013 (). The disadvantage of the use of these modes lies in the detection of frequencies close to those of the Brillouin scattering.
4 Electronic structure
This section is devoted to the main properties of the electronic structure of . Based on first principles calculations, the characteristics of the band structures of singlelayer, doublelayer and bulk are discussed. In particular, we analyze the electronic band gaps and show ab initiohow the band structures depend on the unit cell parameters and the structural optimization. We explain the reasons for the differences in the results obtained through the different computational approaches.
Historically, TMDs are used in the field of tribology as lubricants. The attention given to TMDs decades ago has lead to several theoretical studies of the band structure of in singlelayer and bulk formsMattheiss1973 (); Mattheiss1973a (). These studies were complemented more recently with angleresolved photoelectron spectroscopy measurements for bulk accompanied by ab initio calculations Boker2001 ().
The current interest in RadisavljevicB.2011 (); Peelaers2012 (), the availability of highquality singlelayer flakes Coleman2011 (), and the improvement of experimental results have prompted new theoretical studies in the past 5 years. Regarding the electronic structure, the most efficient approach with respect to computational cost and accuracy is the use of DFTLDA/GGA. Due to the potential application of singlelayer in transistors, most calculations are focused on the band gap. By using LDA, singlelayer MoS is determined as a direct semiconductor with a band gap of 1.78 eV at the point of the Brillouin zone Lebegue2009 (). The transition from a direct to an indirect gap semiconductor with increasing number of layers is also confirmed Scalise2012 (); Scalise2014 (). The extensive use of standard DFT in MoS (and the comparison with more advanced methods) has demonstrated that this approach is adequate for investigating trends. However, an inherent problem of DFTLDA/GGA is the underestimation of the band gap which is related to the lack of the derivative discontinuity in (semi)local exchangecorrelation potentialsJones1989 ().
In a first attempt to improve the electronic band gaps and band dispersions at low computational cost, the modified BeckeJohnson (MBJ) potential becke:jcp:124 () combined with LDA correlation was applied. The MBJLDA approach was developed by F. Tran and P. Blaha tran:prl:102 () and implemented in VASP Kim2010 (). The MBJ potential is a local approximation to an atomic exact exchange potential plus a screening term (screening parameter ) and is given by
(13) 
In this expression, denotes the electron density, the kinetic energy density, and the BeckeRoussel (BR) potential becke:pra:39 (). In the parameterfree MBJLDA calculation employed in this study, the parameter is chosen to depend linearly on the square root of the average of over the unit cell volume and is selfconsistently determined.
Alternatively, the screened hybrid functional HSE heyd:jcp:118 (); paier:jcp:122 (); krukau:jcp:125 () and the improved HSEsol functional schimka:jcp:134 () were tested. The success of HSE in predicting band gaps with accuracy comparable to that of schemes based on many body perturbation theory ( methods) but significantly reduced computational cost is multiply demonstrated in the work of G. Scuseria (see the review in Ref. janesko:pccp:11 ()) as well as in independent studies including a variety of materials and properties paier:jcp:124 (); paier:jcp:125:erratum (); hummer:prb:75 (); paier:prb:79 (); kim:prb:80 (). In general, hybrid functionals are constructed by using the DFT correlation energy and adding an exchange energy that consists of 25% HartreeFock (HF) exchange and 75% DFT exchange. Furthermore, in the concept of the screened HSE functional heyd:jcp:118 () the expensive integrals of the slowly decaying longranged part of the HF exchange are avoided by further separating the into a short (SR) and longranged (LR) term, where the latter is replaced by its DFT counterpart. An additional parameter defines the rangeseparation note:HSE03 (). The HSEsol functional is analogous to the HSE functional, but is based on the PBEsol functional for all DFT parts according to
(14)  
Concerning the electronic properties, the highest level of accuracy has been achieved by calculations. In this work the band structures were studied using the singleshot () and the selfconsistent (sc) approximation. In both approaches, the dynamically (frequency dependent) screened Coulomb interaction and the self energy were calculated using the DFTLDA wave functions. In the sc case, both, the quasiparticle (QP) energies (one electron energies) and the one electron orbitals (wave functions) are updated in and . Note that the attractive electronhole interactions were not accounted for by vertex corrections in . Thus the calculated band gaps are expected to be slightly overestimated compared to experiment shishkin:prl:99 (). In the particular case of , the method has been used in many flavours, yielding a significant spread of values for the band gap, to be discussed below.
4.1 Characterization of the band structure of singlelayer and bulk MoS
First of all, we discuss the main features of the electronic structure of . Figure 10 shows the band structure for singlelayer (1L), doublelayer (2L) and bulk MoS. The vdWDF optimized lattice constant (Section 2) is used in all calculations. The suitability of using either the experimental or the theoretical lattice parameter in the calculations is still controversial. For the calculations presented in this review, we have chosen the latter, which guarantees a strainfree structure. Thereby, we will further be able to investigate strain effects on the electronic structure.
Figure 10 shows the relevant features of the bulk band structure:

two distinct valence band edges (VBEs) located at and ,

three conduction band extrema (CBEs) at (half way between and ), , and (half way between and ),

the valence band maximum (VBM) located at and the conduction band minimum (CBM) at ,

a fundamental electronic band gap of indirect nature that is defined by the energy difference ,

the splitting of the VBM at into states and due to interlayer interaction,

twofold spin degeneracy of all states due to inversion symmetry, and

nearly parabolic band dispersions at , , and .
On the level of accuracy the CBM in bulk ( point) is 0.4 eV lower in energy than the CBE at and the VBM at is favored by 0.5 eV energy difference with respect to the band edge at . In bulk , the valence band splitting at amounts to roughly 240 meV.
In contrast to bulk and multilayer , the main attribute of the singlelayer band structure is the direct fundamental band gap located at since both, the VBM and the CBM, are at . This direct band gap at has been clearly demonstrated by several preceding works on singlelayer Komsa2012 (); Ramasubramaniam2012 (); Yun2012 (); Cheiwchanchamnangij2012 (); Shi2013 (); MolinaSanchez2013 () and explains the significant increase of photoluminescence yield in singlelayer Mak2010 (). The other important feature in the bandstructure of 1L is the splitting of the valence band maximum which amounts to 150 meV. Since interlayer interactions are absent, this splitting must have a origin different from the splitting in bulk. Indeed, due to missing inversion symmetry, the spindegeneracy of the bands is lifted, resulting in a particularly large spinorbit splitting at zhu:prb:84 (); Cheiwchanchamnangij2012 (). This splitting explains the doublet structure of the strong peak in the absorption spectrum of 1L Mak2010 ().
The CBM in 1L is also located at but the splitting due to SOC is one order of magnitude weaker than the splitting of the VBM. Its absolute value is strongly affected by the exchange correlation functional used in the calculations Kuc2015 () as will be discussed later. Both, the valence and conduction bands exhibit nearly parabolic dispersion at this point, which explains the small effective charge carrier masses and indicates promising conductivity. However, compared to bulk the valence band at is considerably flattened. This flattening results in a much higher effective hole mass in 1L, which was also observed in AngleResolved Photoemission Spectroscopy (ARPES) experiments Jin2013 (). A second local conduction band minimum close in energy is observed at . The relative energy positions of the states and and the location of the VBM (either or ) determine, whether the material is a direct or an indirect semiconductor. We observed that in 1L the energy difference is very sensitive to (i) the structural optimization, (ii) the applied inplane strain, and (iii) the accuracy (around 0.05 eV). We discuss these issues in more detail later in this section.
We now describe the changes stemming from quasiparticle corrections in the bandstructures of bulk, single, and doublelayer . The most notable change is the sizable increase of the band gap on the level of the method. Also the valence band width increases slightly. Note that for the singlelayer, the VBM at is shifted downwards as compared to the VBM at . The conduction band is even more profoundly affected. The upshift of the CBM at is larger than the one of the secondary CBM . In the singlelayer, this results in the lowering of the energy relative to and thus in the reduction of the energy difference by 60% compared to LDA turning the material almost indirect on the level. In bulk the CBM is lower in energy than (due to interlayer interaction) already on the DFT level. On the level, this trend is amplified even more. In both, doublelayer and bulk , the corrections
Besides that, one has also to consider that the correction to the band structures slightly depends on the number of layers in multilayer systems. We find that the band gap correction at decreases in doublelayer and bulk with respect to the case of singlelayer . The larger number of layers is associated with an increase of the dielectric screening, which results in a smaller correction Wirtz2006 (); Komsa2012 ().
In order to better understand, the origin of the differences between single, double, and bulk band structures, we have summarized the orbital composition of the highest valence and lowest conduction bands at points of interest in the Brillouin zone in Table 5. In singlelayer , the S and Mo orbitals dominate the composition of the wave functions, with a minor contribution of the orbitals. The conduction band edge at is mainly a Mo (86 %) and the remaining part shared between the S and Mo orbital. The valence states at and are predominantly Mo (80 %) and 20 % of S without any contribution from orbitals. The wave function at the local minimum at has a more complex composition, typical for points of low symmetry, as summarized in Table 5. These findings qualitatively reproduce previous DFTPBE results (e.g., Fig. 4 in Ref. Scalise2014 (), Fig. 5 in Ref. Guzman2014 ()) as well as the tightbinding (TB) model of Liu and coworkers Liu2013 (); Liu2015 (). The latter, however, suggest also a significant contribution of the Mo states to the valence band edge at which we contribute to a deficiency of the TB model using only three Mo bands. The composition of the and states in bulk is very similar to the singlelayer values. The bulk and states however, are now predominantly composed of Mo and the valence band states at change the weight of the orbital of sulphur atoms. The latter increase is related to the bonding between sulphur atoms of different layers, which produces the interlayer coupling Cappelluti2013 ().
Singlelayer  

Atom  Sulphur  Molybdenum  
Orbital  
  0.54        0.46        
    0.23  0.02      0.75      
  0.09    0.05      0.86      
  0.20              0.80  
  0.20              0.80  
0.03  0.22  0.06    0.54    0.12    0.01  
Bulk  
Atom  Sulphur  Molybdenum  
Orbital  
  0.53        0.47        
0.07    0.30  0.03      0.60      
  0.09  0.05        0.86      
  0.21      0.79        
  0.18      0.82        
0.07  0.18  0.06    0.52    0.14    0.03 
A correct description of the electronic properties requires (i) the inclusion of Molybdenum semicores states ( and orbitals) in the basis set, (ii) a plane wave cutoff of 350 eV, (iii) at least a () centered mesh for bulk (1L) , and (iv) the explicit inclusion of the spinorbit interaction Sangalli2012 (). We interpolate the band structure to a finer grid using the WANNIER90 code mostofi:cpc:178 () and the VASP2WANNIER90 interface vaspwannier (). With respect to calculations, it is important to mention that (i) solely including valence electrons leads to an erroneous wavevector dependence of the correction MolinaSanchez2013 (), (ii) the convergence with respect to virtual states when calculating is particularly slow for 1LQiu2013 (), and (iii) the default value for the number of quasiparticle energies that are calculated and updated in the sc VASP calculation must be substantially increased. While taking into account 500 virtual states is sufficient to converge the quasiparticle gaps of bulk within 20 meV, more than 1000 bands are required for 1L gaps to be stable within 40 meV. The logarithmic scaling of the direct and indirect gap in 1L with respect to the number of bands included in the calculation is shown in Fig. 11.
Concerning the number of quasiparticle energies that have been updated in the sc calculations (NBANDSGW parameter in VASP), it is emphasized that more than 200 are required to converge the quasiparticle gaps. In particular the conduction band extremum at point strongly depends on this parameter.
4.2 Dependence on the crystal structure
The analysis of the preceding paragraphs underlines the importance of accurately calculating the energy difference between the conduction band minima at and . This is a challenge for the different theoretical approaches mentioned before, because these quantities also sensitively depend on the details of the crystal structure. In order to discuss this, we focus on singlelayer and bulk but the conclusions can be extended to multilayer .
One source of controversy between the results of several calculations performed for 1L Komsa2012 (); Ramasubramaniam2012 (); Cheiwchanchamnangij2012 (); Shi2013 () could be the underlying crystal structure. In particular, the lattice constant, interlayer distance (relevant for multilayer and bulk ), and atomic positions defining the interatomic distance (MoS bond length and SMoS bond angle) may significantly affect the energy gaps and band dispersion. The dependence of the band structure on the details of the crystal structure has not been addressed so far and will be elucidated in the following.
Most calculations reported so far, have used the experimental room temperature lattice constant of bulk , dickinson:jacs:45 () i. e., Å, and 1015 Å vacuum along the axis for the singlelayer (1L) slab structure. However, less information is given about the choice of the origin of the unit cell (atomic positions) and the parameter. Unfortunately, according to Bronsema et al. bronsema:zaac:540 () an inconsistency exists in literature regarding the choice of the atomic positions and the corresponding parameter. In fact, the latter determines the interatomic distances and thus the SMoS layer thickness. For this reason, the band structure of bulk and 1L has been calculated by LDA+SOC and +SOC for some of the crystal structures summarized in Tab. LABEL:tab:lattice. The +SOC results are depicted in Fig. 12.
For bulk , small variations of the parameter (compare red and blue lines in Fig. 12) and the lattice constants and (compare the red and green lines in Fig. 12) have only little influence on the bandstructure, since there is anyway a strong interlayer interaction that leads to a splitting of the valence and conduction bands. The situation is quite different for the singlelayer: as illustrated in Fig. 12(b). The significant role of the internal parameter that is defining the interatomic distances (MoS bond length and SMoS bond angle) is revealed. With increasing from 0.621 to 0.629, the MoS bond length is reduced from 2.42 Å to 2.35 Å, respectively. This favors hybridization between the Mo and S states that comprise the highest occupied and lowest unoccupied bands and therefore the band dispersion (band width) increases. As a consequence, the CBE at and the VBE at are pushed to higher energies and the CBE at becomes the CBM giving rise to a direct fundamental energy gap. It is strongly emphasized that an improper choice of atomic positions and corresponding parameter can fortuitously yield a direct fundamental band gap in 1L and may partly explain the inconsistency among band structuresKomsa2012 (); Ramasubramaniam2012 (); Cheiwchanchamnangij2012 (); Shi2013 () reported so far.
Another source of discrepancies between singlelayer calculations might be related to the relaxation of the atomic positions. In Fig. 13(a), the band structures of 1L calculated using the experimental bulk unit cell lattice constant (=3.16 Å) without and with LDArelaxed atomic positions in the singlelayer are compared. Due to the overbinding in DFTLDA, the MoS bond length gets reduced with relaxation, which again strengthens the Mo–S hybridization resulting in an increase of the band dispersion along . This results in a raise of the VBM at accompanied by an increase of the conduction band valley energy and the stabilization of the CBM at . Consequently, +SOC yields the correct direct gap band structure for 1L, if the atomic forces are minimized on the LDA level. As can be seen from Fig. 13(a), the direct gap at is reduced in the position relaxed case by 90 meV. While the CBM energy at is not affected, further reduction of the direct gap at by 40 meV is obtained by using the lattice constant of the optB86bVdW fully relaxed bulk unit cell (=3.164 Å) and LDArelaxed atomic positions in the singlelayer [not shown in Fig. 13(a)].
A final test on the level of for the influence of the structural details on the band structure of 1L was performed with a fully optB86bVdW optimized singlelayer structure, i. e., the inplane lattice constant = 3.162 Å and LDArelaxed atomic positions. Since the VdW interactions are not relevant in the singlelayer, the obtained structure is very close to the experimental bulk one. Thus the +SOC band structure resembles that one calculated without any atomic position relaxation [indirect gap, not shown in Fig. 13(a)]. From this analysis we conclude, that the location of the valence and conduction band extrema at and are very sensitive to the relaxation of the atomic positions and if the atomic positions in the singlelayer are relaxed within DFTLDA, the CBM at is stabilized with respect the CBE at .
The finding of an direct gap 1L on the level seemed to be controversial to the results obtained by Shi et al. Shi2013 (), who performed an analogous comparison of and calculations for 1L as presented in Fig. 13(b) and concluded that singleshot is insufficient in describing 1L. As shown here, omitting the relaxation of the atoms in the singlelayer structure constructed from the experimental bulk structure leads indeed to an incorrect description of 1L on the level (indirect band gap). The calculation cures this problem, but further increases the direct gap at . The tendency of to overestimate semiconductor band gaps is known shishkin:prl:99 () and thus one must assume that the direct gap of 1L is too large.
4.3 Performance of different methodologies
After the discussion of the relation between structural and electronic properties, we focus on how the results depend on the XC functionals. The results of this analysis for bulk, double, and singlelayer are summarized in Table 6 and Fig. 14. Note, the fully optimized structures using the van der Waals functional as previously described were used in these calculations.
(eV)  

Bulk  
LDA  1.64  1.86  1.07  1.40  0.83  0.24  2.08 
LDA Peelaers2012 ()  (1.80)  (0.81)  
PBEsol  1.65  1.87  1.10  1.42  0.87  0.23  2.11 
PBE Yun2012 ()  (0.87)  
PBE Peelaers2012 ()  (1.58)  (0.86)  
MBJ  1.62  1.82  1.21  1.56  1.15  0.06  2.27 
HSEsol  2.10  2.36  1.58  1.96  1.44  0.14  2.84 
HSE Peelaers2012 ()  (2.16)  (1.48)  
2.08  2.32  1.63  1.69  1.24  0.39  2.53  
Jiang2012 ()  (2.07)  (1.23)  
Komsa2012 ()  (2.00)  (1.30)  
sc  2.17  2.41  1.59  2.02  1.44  0.16  2.88 
sc Cheiwchanchamnangij2012 ()  (2.099)  (2.337)  (1.287)  
EXPT. Mak2010 ()  (1.78)  (1.29)  
EXPT. Frey1998 ()  (1.95)  (1.20)  
Single layer  
LDA  1.62  1.77  1.61  1.92  1.91  0.30  2.74 
PBEsol  1.65  1.79  1.69  1.89  1.93  0.24  2.78 
PBE Yun2012 ()  (1.75)  
PBE Ramasubramaniam2012 ()  (1.60)  
HSEsol  2.09  2.28  2.23  2.45  2.59  0.36  3.63 
HSE  2.06  2.25  2.13  2.51  2.58  0.45  3.60 
HSE Ramasubramaniam2012 ()  (2.05)  
2.45  2.60  2.61  2.59  2.74  0.14  3.60  
Ramasubramaniam2012 ()  (2.82)  
Komsa2012 ()  (2.97)  (3.26)  
sc  2.72  2.87  2.90  2.98  3.16  0.26  4.29 
sc Shi2013 () (w/o SOC)  (2.78)  
sc Cheiwchanchamnangij2012 ()  (2.759)  (2.905)  
EXPT. Mak2010 ()  (1.88)  (2.05)  (1.6) 
On all levels of theory the band structure of bulk shown in Fig. 14(a) and (b) corresponds to an indirect semiconductor with the VBM located at and the CBM at .
When comparing the numbers given in Table 6, LDA and PBEsol underestimate the transition (indirect gap of 0.85 eV) are found to be in agreement with previous calculations Yun2012 (); Ramasubramaniam2012 (). Compared to LDA, the inclusion of local exchange as provided by the MBJ potential mainly affects the and energies resulting in a larger indirect gap of 1.15 eV. However, the band dispersion along towards is reduced resulting in almost energetically balanced CBEs at and , i. e., is only 60 meV. Therefore the difference between indirect and the direct gap is less pronounced than in LDA. Further improvement is achieved by taking into account nonlocal exchange using the HSEsol functional that increases the indirect as well as the direct gap. The difference in HSEsol is roughly 140 meV, significantly less compared to LDA ( 240 meV). Taking into account the dielectric screening in the approach strongly enhances the difference to 390 meV. This results in an indirect for bulk of 1.24 eV, in good agreement with experiment (1.21.3 eV) as well as previous calculations. Jiang2012 (); Komsa2012 () The outstanding agreement of both, the indirect (1.24 eV) and direct (2.08eV) bulk gaps, with the values obtained by an allelectron code Jiang2012 () verify the accuracy of the present PAW results.
A higher level of accuracy is reached by performing sc calculations. Compared to the band structure, the difference is again significantly reduced (to 160 meV), since is pushed up in energy almost back to the HSEsol position. The fundamental indirect sc gap is 1.39 eV and slightly overestimated compared to experiment, which is due neglecting the attractive electronhole interactions via Vertex correction in . shishkin:prl:99 () In the – region, the sc band structure resembles the HSEsol, whereas we observe remarkable differences in the –– range. The and sc results are consistent to previous calculations Cheiwchanchamnangij2012 () given in Tab. 6 within the uncertainties originating from computational aspects.
Singlelayer is described as a semiconductor with a direct gap at on all levels of theory beyond standard DFTLDA, provided that the crystal structure is fully relaxed as stressed in Sec. 2. Standard DFT (LDA and PBEsol) severely underestimates the direct gap of the singlelayer structure. Besides that, LDA wrongly sets the VBM at at slightly higher energy than the VBE at . It is important to emphasize that the underestimated energy difference is observed if the singlelayer slab is constructed from the optB86bVdW relaxed bulk structure (inplane lattice constant =3.164 Å), but not in case of the experimental bulk lattice constant =3.16 Å. Therefrom we conclude that the relative positions of the and energies are strongly dependent on the inplane lattice constant. Hence is imperative the investigation of strain effects on the 1L band structure presented later in this section.
Employing the HSEsol functional to 1L shifts the conduction bands almost uniformly upwards in energy compared to DFTLDA resulting in a rather constant , i. e., 360 meV and 300 meV, in HSEsol and LDA, respectively. Band dispersions in the valence band region are increased within the HSEsol description and the VBM splitting at is enhanced to 200 meV compared to the LDA value of 150 meV. The VBM at is stabilized by roughly 200 meV compared to in HSEsol calculations. Concerning approaches, the subtle changes of the band dispersions between and result in a significant change of the energy difference: stabilizes the CBM at by 130 meV, whereas sc enhances this energy difference by a factor of two. Note that the energy differences strongly depend on the total number of bands (NBANDS) included in the calculations (see Fig. 11) and the convergence with NBANDS itself is influenced by the amount of vacuum included in the single layer cell (20 Å in the present case). This means that using a larger amount of vacuum requires an increase of the NBANDS parameter as well. For this reason, the comparison between the present results and previously reported values, as summarized in Tab. 6, is difficult. The values listed in Tab. 6 refer to calculations with NBANDS=512. Increasing NBANDS from 512 to 1920 reduces the direct gap at by roughly 80 meV, the indirect gap by 60 meV, but the energy difference increases by 40 meV.
Analogous to bulk , including nonlocal exchange by HSEsol increases the gap (2.09 eV) considerably. The calculated quasiparticle gap amounts to 2.45 eV, which is smaller by 0.30.5 eV to reported values. Ramasubramaniam2012 (); Komsa2012 (); Liang2013 () This difference is attributed to structural and computational details: A calculation performed with the experimental crystal structure and a reduced mesh of 882 yields 2.86 eV. Liang et al.reported a direct band gap of 1L of 2.75 eV, Liang2013 () which was obtained by calculations taking into account a Coulomb interaction truncation to avoid spurious interlayer interaction between the periodically repeated monolayers, but using the generalized plasmonpole model (GPP) for the dynamical screening and omitting SOC. The issues of the Coulomb interaction truncation, point sampling, and vacuum layer thickness were also addressed by Hüser et al.Huser2013 (), who argued that the band gap values converged with respect to point sampling and slab distance are rougly 0.4 eV too small compared to the free standing monolayer (including Coulomb truncation). Once again this reflects the difficulty to achieve accurate results and explains the plethora of band gap data in literature.
Compared to the singleshot result for the direct gap, the sc further increases the gap to 2.72 eV in agreement with previously reported values as listed in Table 6. The slow convergence of the 1L band gaps with the NBANDS parameter was only recently stressed Qiu2013 () and might explain the larger values previously reported.
At this point it is important to recall that quasiparticle gaps are singleparticle gaps. Their overestimation by roughly 0.5 eV compared to some experiment is explained by the missing electronhole interactions (excitonic effects), which are strong in 2D materials due to confinement and lead to the formation of bound electronhole pairs. These bound excitons reduce the direct band gap by their binding energy and define the optical gap, which is experimentally accessed by optical measurements such as photoluminescence or photoconductivity. Excitonic effects are addressed in the next section.
To conclude, it is evident based on the above discussion that the different levels of accuracy and/or complexity applied in the methods substantially alter the results. The nonlocal exchange and dynamical screening are inevitable for an accurate description of the electronic properties of . Based on our calculations of singlelayer , a reasonable estimate for the direct band gap is 2.40.2eV and the spinorbit splitting of the valence band edge at is 150160 meV. The energy difference between the two valence band extrema in 1L is much smaller than in bulk and very sensitive to the inplane lattice constant. As put forward by Kuc and HeineKuc2015 (), the estimations of the weak spinorbit splitting of the conduction band edge is strongly dependent on the XC functional used in the DFT calculation and a better description by methods beyond ground state DFT is required. From our calculations, we deduce a conduction band splitting of 10 meV, which however falls within the estimated uncertainty range.
A final remark concerns the starting point of the calculations. One should keep in mind that the result can be influenced by the wave functions (orbitals) used for calculating and as pointed out in Ref. Caruso2012 (). Thus, in the complexity of the method one can go a step forward by applying a full selfconsistent procedure in the Green’s function and the screened interaction , consisting in the iterative solution of the Dyson equation. However, the extraordinarily cumbersome calculations required for this procedure restrict the application of this approach to small systems, such as binary molecules as N or CO Caruso2012 (). Up to now this scheme has not been applied to singlelayer and we think that its implementation for layered materials is still far.
4.4 Strain effects in singlelayer
The ideal scenario of freestanding 2D layers as considered in most theoretical simulations is hardly fulfilled in reality. In the course of experiments or device fabrication with 2D materials, it is important to consider strain resulting e.g. from the lattice constant mismatch between the substrate and the 2D layer. Equally important in this context is the interaction of the 2D material with the substrate as shown in Ref. Jin2013 (). Therein, the ARPES scans of exfoliated singlelayer compared to those of chemical vapor deposition grown singlelayer on silicon revealed that the presence of substrate alone is sufficient to modify the band structure. In particular, the substrate interactions are responsible for the pronounced flattening of the VBM at of on silicon.
In addition, recent experiments have demonstrated that application of tensile strain changes the gap from direct to indirect Conley2013 (). In particular, the flake deposited on a flexible substrate which is subsequently deformed in a controllable manner, experiences uniaxial tensile strain up to 2.2 %. The photoluminescence spectra of these samples show a clear transformation of the band character, and an associated reduction of the integrated intensity of the optical signal.
The sensitivity of TMDs band structure on the lattice constant opens the possibility to modify the band gap and thus the optical properties in a controlled way by external strain. This issue has been theoretically addressed either through LDA/GGA calculations Ellis2011 (); Scalise2012 (); Scalise2014 (); Dong2014 (), or the method Shi2013 (). The effect of hydrostatic pressure on the vibrational, electronic, and optical properties of bulk, multi, and single layer was investigated by Nayak et al.Nayak2014 (); Nayak2015 () by combining various experiments (high resolution transmission electron microscopy, electrical resistance measurements, laser Raman spectroscopy, synchrotron Xray diffraction experiments under highpressure) with DFT calculations. Interestingly, while the direct bulk band gap decreases with increasing pressure, the direct band gap of 1L increases by 11.7% up to 12 GPa before it is reduced. Thus the pressure induced electronic transition from the semiconducting to a semimetallic state occurs at much larger pressures in the latter. Nayak2015 ()
Being aware of the importance of substrate interactions, we investigated the strain effects on the electronic properties of 1L within the approach and the model of freestanding 2D layers. Biaxial tensile strain has been realized by increasing the inplane lattice constant of 1L. The band structures of the strained materials were calculated with relaxed atomic positions and are shown in Fig. 15. The direct gap and interband transitions as a function of strain deduced from these band structures are collected in Tab. 7.
strain (%)  +SOC  sc+SOC  

(eV)  (eV)  (eV)  
0  2.51  2.71  0.08  2.76  2.85  0.39 
1  2.37  2.44  0.24  2.61  2.58  0.55 
2  2.23  2.18  0.41  2.47  2.31  0.72 
3  2.11  1.93  0.55  2.34  2.06  0.86 
4  1.98  1.69  0.69  2.21  1.82  0.99 
With increasing inplane lattice constant (biaxial tensile strain) the bond distances within the plane of the MoSMo sheets are changed, but also in the perpendicular direction through the relaxation along . According to Tab. 5, the valence states at (VBM) are mainly composed of S and Mo orbitals, whereas the valence states at have predominantly Mo character. Concerning the CBM at , the states are mainly Mo orbitals and the conduction band states around have predominantly Mo character. By changing the SMoS bond lengths and angles due to tensile strain, the overlap of the Mo with the S is reduced, whereas the coupling between the Mo and S is increased.Guzman2014 () As a consequence, the energy raises with respect to and the energy decreases compared to