Effects of strain on band structure and effective masses in MoS
2
H. Peelaers
∗
and C. G. Van de Walle
Materials Department, University of California, Santa Barbara, California 93106-5050
(Dated: November 27, 2012)
We use hybrid density functional theory to explore the band structure and effective masses of
MoS
2
, and the effects of strain on the electronic properties. Strain allows engineering the magni-
tude as well as the nature (direct versus indirect) of the band gap. Deformation potentials that
quantify these changes are reported. The calculations also allow us to investigate the transition in
band structure from bulk to monolayer, and the nature and degeneracy of conduction-band valleys.
Investigations of strain effects on effective masses reveal that small un iaxial stresses can lead t o large
changes in the hole effective mass.
PACS numbers: 71.20.Nr,71.70.Fk,73.22.-f
2
(a)
(b)
(c)
M
K
L
H
A
Γ
Λ
Σ
T
Q
∆
P
FIG. 1. (Color online) (a) Side view of the bulk unit cell of MoS
2
. (b) A monolayer of MoS
2
. (c) The hexagonal Brillouin
zone. High symmetry points and lines are indicated. The conduction- band valleys at Λ
min
and K are schematically depicted,
with ellipsoids representing constant-energy surfaces.
Molybdenum disulfide (MoS
2
) is a semiconducting mater ial, widely used as a dry lubricant because of its structural
similarity to graphite. It consists of stacked hexagonal S-Mo-S layers (Fig. 1). These layers, conventionally referred
to as monolayers of MoS
2
, are weakly bound by van de r Waals forces. Similarly to the production of graphene, MoS
2
samples consisting of a single or a few monolayers can by produced by micromechanical exfoliation
1
—but, contrary to
graphene, MoS
2
actually has a band gap. Mono layer MoS
2
field-effect transistors have already been demonstrated,
2
but devices based on multilayers show great promise a s we ll.
3
Many details about the electronic prope rties of bulk MoS
2
are still lacking. In this work we report comprehensive
results fo r band structure, addressing direct and indirect band gaps a nd multiple conduction-band valleys, as well as
effective masses. We also investigate how str ain affects these properties, in the process also clarifying the differences
in band structure between bulk and monolayer. Strain effects on monolayers and bilayers have already been inves-
tigated,
4,5
but results for bulk MoS
2
are not yet available. Strain can result from externally applied stress, or arise
from pseudomorphic growth
6–8
or when a MoS
2
layer is clamped to a subs trate. We quantify the changes in the band
structure as a function o f strain in terms of deformation potentials. All of the quantities rep orted here ar e relevant
for further development of electronic applications of this material as well as for device modeling.
We investigate the two types of deformations that are relevant for MoS
2
: uniaxial and biaxial. Uniaxial strain
parallel to the c direction, which we denote by ε
zz
, directly affects the interlayer separation. In the limit of large
tensile strain, interlayer distances become large and interactions negligible, and we effectively reach the monolayer
limit. Thus the investigation of the effects of tensile uniaxial str ain provides insight into the trans ition betwee n
the bulk system and the monolayer . We impose a given value of the c lattice parameter, corresponding to a given
ε
zz
, and allow for a relaxation of the lattice parameters in the perp e ndic ular direction, as wo uld occur in a realistic
uniaxial stress geometry. For biaxial stress, we impose a value of the in-plane lattice parameters a (or equivalently
b), corresp onding to a strain ε
xx
= ε
yy
. In this case we also allow for lattice re laxation in the dir e ction parallel to c,
corres po nding to biaxial stre ss.
3
TABLE I. Band-gap energies (in eV) for bulk MoS
2
calculated using different methods.
LDA GGA HSE06 HSE06
15
GW
0
17
QSGW
16
Experiment
18
Indirect 0.81 0.86 1.48 1.48 1.23 1.29 1.29
Direct 1.80 1.58 2.16 2.33 2.07 2.10 1.95
All calculations are based on generalized Kohn-Sham theory with the HSE06 screened hybrid functional
9
and the
projector augmented-wave (PAW) pseudopotential-plane-wave method
10
as implemented in the VASP code.
11
For the
PAW pseudopotential fo r Mo we included the full n=4 shell (4s
2
, 4p
6
, and 4d
5
) plus 5s
1
as valence. For S, the n =3
shell is included as valence (3s
2
and 3p
4
). A 10×10×2 Monkhorst-Pack
12
k-point grid was used for all calc ulations
and a plane-wave basis set with an energy cutoff of 2 80 eV.
The hybrid functiona l approach is well suited to des c ribe both structural properties and band str uc tures,
13
in par-
ticular the band gaps of semiconductors, which are severely underestimated when using standa rd exchange-correlation
functionals such as the local-density approximation (LDA) or the generalized gradient approximation (GGA). Van
der Waa ls interactions, which govern the interlayer distance in MoS
2
, are not explicitly included in HSE06; in fact,
this is still an area of active research within density functional theory. We therefore fix the c lattice parameter (which
is overestimated by 5.6 % in HSE06) to its experimental value (12.29
˚
A, Ref. 14), a practice also applied in previous
computational studies.
15,16
Table I compares our calculated band gaps of bulk MoS
2
with prev iously published values. Bulk MoS
2
is an
indirect-band-gap material, with the valence-band maximum (VBM) located at the Γ point and the conduction-band
minimum (CBM) at a point on the Γ-K line; we denote this minimum as Λ
min
[see Fig. 1(c)]. We also examine the
direct K-K band gap. Both LDA and GGA underes tima te the band gaps. In the GW
0
method
17
the one-electron
Green’s function G is self-consistently updated, while the screened Coulomb interaction W is fixed at its initial
value. QSGW stands for the quasi-pa rticle self-consistent GW method utilized in Ref. 16. The GW methods slightly
underestimate the indirect band gap, but overestimate the direct gap. Our approach based on the HSE06 functional
slightly overestimates the gaps compared to experiment—though we note that the optica lly measured gaps probably
reflect excitonic contributions that are not included in our calculations. Overall, this comparison confirms that the
hybrid functional is the right choice for our study, since in addition to providing band structures of comparable
or better quality than those obtained with other approaches, it allows calc ulating forces (not available in the GW
approaches), needed to relax the system in respo nse to applied stresses.
Figure 2 shows the calculated band structure fo r different values of uniaxial stress along the c direction. The
valence-band maximum (VBM) of the bulk is chosen as the zero energy reference, and the band structures are aligned
using the Mo 4s electrons as reference states. Since we are interested in a comparison with the monolayer and in
transport within the layers, we focus on the in-plane part of the Brillouin zone.
As noted above, bulk MoS
2
has an indirect band gap Γ-Λ
min
. The monolayer, on the other hand, has a direct band
gap, located at the K point. As the layers are moved apart (corresponding to an increase in ε
zz
), a transition from
the indirect band gap to the direct gap occurs, but only when the interlayer distance is increased by almos t 50%. As
seen in Fig. 2 , the VBM at Γ moves down in energy, while the VB shift at the K point is much smaller. The location
of the VBM thus shifts from Γ to K. Simultaneously the Λ
min
CBM moves higher in energy compared to the K point,
shifting the over all minimum to the K point. The switch of location of the CBM oc c urs for much smaller strains than
the switch of the VBM. Therefore, with increasing tens ile strain the nature of the band gap switches from indirect
Γ-to-Λ
min
, via indirect Γ-to-K, to direct K-to-K, as illustrated in Fig. 3 (a).
The relative shift of the VB-1 and the CB+1 bands is opposite to that o f the highest VB and lowest CB, as these
bands originate fr om the two different S- Mo-S layers that form the bulk unit c ell [Fig . 1(a)], and thus have to merge
into one degenerate band for lar ge ε
zz
(large interlayer distance).
For compressive strains (ε
zz
<0), the VBM at Γ moves up and the CBM at K moves down; the band gap remains
indirect. We note that in unstrained bulk the CBM a t Λ
min
is only slightly higher in energy than the minimum at K,
and the energy difference between the Γ-K and Γ-Λ
min
band gaps remains smaller than 0.02 eV fo r strains les s than
5%. This near-degeneracy, along with the multiplicity o f these off-Γ CB valleys, is an important feature for device
applications since it affects the density of states. The CB valleys are schematically depicted in Fig. 1(c), where the
Λ
min
has a multiplicity of 6 (inner ellipses) and the valley at K a multiplicity of 2 (ellipses at edge, 6 · 1/3 = 2).
Our result for a near-degeneracy between the Λ
min
and K CB minima somewhat differs from a previous report
based o n HSE06 hybrid functional calculations,
15
in which the CB at K was found to be ∼0.3 eV higher in energy
than at Λ
min
. Possible reasons for the discrepancy could be the use of a different basis set (Gaussians versus plane
waves) or the use of pseudopotentials versus PAWs. This c an also be seen in Table I, where the direct band gap
(from K to K) differs in magnitude between our appr oaches. We note that our results are in good agreement with the
quasiparticle selfconsistent GW results of Ref. 16.
A similar study can be perfo rmed for biaxial stress, where now the c lattice parameter is allowed to relax in response
4
-2
-1
0
1
2
3
Γ Σ M T K Λ Γ
Energy (eV)
ε
zz
=-0.1
ε
zz
= 0
ε
zz
= 0.1
layer
FIG. 2. (Color online) Band structure of bulk MoS
2
under uniaxial stress along the c-axis, plotted along high-symmetry lines
in the in-plane Brillouin zone. The band structure of unstrained bulk and an unstrained monolayer are also shown. For clarity
we only show the highest two valence bands and lowest two conduction bands. The labels underneath the segments of the
Brillouin zone indicate the names of the high-symmetry lines [see Fig. 1(c)].
to an applied in-plane strain. The biaxial stress results are shown in Fig. 3(b); we used the data from Fig. 3(a) to
correct for the difference between the calculated and experimental values of the c lattice parameter. Fig. 3(b) shows
that for positive (tensile) strain the CBM at Λ
min
shifts up, a nd the CB at K shifts down. The VBM remains at Γ.
The band gap thus remains indirect for all biaxial stress co nditions, switching from Γ-to- Λ
min
to Γ-to-K at a very
small tensile strain.
The variation of band gaps with strain can be expressed in terms of deformation potentials, which are useful
quantities for device modeling. For small strains in the v icinity of the equilibrium lattice parameters the va riation of
the gaps is approximately linear , and hence can be expressed as
∆E
Γ-K
g
= D
Γ-K
zz
ε
zz
and ∆E
Γ-K
g
= D
Γ-K
⊥
(ε
xx
+ ε
yy
), (1)
where the second deformation potential is expressed with respect to the sum of the in-plane-strains, in analogy with
the definition of the deformation potentials D
2
and D
4
in a hexagonal system.
19
Similar expressions hold for the
Γ-to-Λ
min
and K- to-K gaps.
The calculated deformation potentials are lis ted in Table II. The range over w hich the linear approximation inherent
in Eq. (1 ) is valid is larger for biaxial strains than in the uniax ial case (see Fig. 3); for the latter, a change in c by 5%
changes the deformation potentials by as much as 1 eV.
Note that Eq. (1) is formulated in terms of strain, not stre ss, and hence the deformation potentials are calculated by
varying only the c lattice parameter (in the case of D
zz
), with the a and b lattice parameters fixed to their equilibrium
values. To obtain the variations under uniaxial or biaxial s tress conditions that were depicted in Fig. 3, strains along
directions parallel as well as perpendicular to c need to be taken into account. Under uniaxial stress, the ratio between
ε
xx
and ε
zz
is given by Poisson’s ratio ν:
ε
xx
= −νε
zz
with ν = c
13
/(c
11
+ c
12
) , (2)
where the c
ij
are the elastic constants of MoS
2
. Under bia xial stress, we have:
ε
zz
= −(2c
13
/c
33
)ε
xx
. (3)