Landau levels and Shubnikov-de Haas oscillations in monolayer transition metal dichalcogenide semiconductors

New J. Phys.17(2015)103006 doi:10.1088/1367-2630/17/10/103006

PAPER

Landau levels and Shubnikov – de Haas oscillations in monolayer transition metal dichalcogenide semiconductors

Andor Kormányos¹, Péter Rakyta^2,3and Guido Burkard¹

1 Department of Physics, University of Konstanz, D-78464 Konstanz, Germany

2 Department of Theoretical Physics, Budapest University of Technology and Economics, H-1111 Budafoki út. 8, Hungary

3 MTA-BME Condensed Matter Research Group, Budapest University of Technology and Economics, H-1111 Budafoki út. 8, Hungary E-mail:andor.kormanyos@uni-konstanz.deandguido.burkard@uni-konstanz.de

Keywords:transition metal dichalcogenides, Landau levels, Shubnikov–de Haas oscillations, 2-dimensional systems

Abstract

We study the Landau level

(LL)

spectrum using a multi-band

k p·

theory in monolayer transition metal dichalcogenide semiconductors. We

find that in a wide magneticfield range the LL can be

characterized by a harmonic oscillator spectrum and a linear-in-magnetic

field term which describes

the valley degeneracy breaking. The effect of the non-parabolicity of the band-dispersion on the LL spectrum is also discussed. Motivated by recent magnetotransport experiments, we use the self- consistent Born approximation and the Kubo formalism to calculate the Shubnikov–de Haas

oscillations of the longitudinal conductivity. We investigate how the doping level, the spin-splitting of the bands and the broken valley degeneracy of the LLs affect the magnetoconductance oscillations. We consider monolayer MoS

2

and WSe

2

as concrete examples and compare the results of numerical calculations and an analytical formula which is valid in the semiclassical regime. Finally, we briefly analyze the recent experimental results

(Cuiet al

2015

Nat. Nanotechnol.10534)

using the theoretical approach we have developed.

1. Introduction

Atomically thin transition metal dichalcogenides semiconductors(TMDCs)[1–3]are recognized as a material system which, due to itsfinite band gap, may have a complementary functionality to graphene, the best known member of the family of atomically thin materials. The experimental evidence that TMDCs become direct band gap materials in the monolayer limit[4]and that the valley degree of freedom[5]can be directly addressed by optical means[6–9]have spurred a feverish research activity into the optical properties of these materials[10–

13]. Equally influential has proved to be the fabrication of transistors based on monolayer MoS2[14]which motivated a lot of subsequent research to understand the transport properties of these materials. Achieving good Ohmic contact to monolayer TMDCs is still challenging and this complicates the investigation of intrinsic properties through transport measurements. Nevertheless, significant progress has been made recently in reducing the contact resistance by e.g., using local gating techniques[15], phase engineering[16], making use of monolayer graphene as electrical contact[17–19], or selective etching procedure[20].

Our main interest here is to study magnetotransport properties of monolayer TMDCs. Unfortunately, the relatively strong disorder in monolayer TMDC samples have to-date hindered the observation of the quantum Hall effect. Nevertheless, the classical Hall conductance has been measured in a number of experiments [15,18,21–23]and was used to determine the charge densityneand to extract the Hall mobilityμH. In addition, three recent works have reported very promising progress in the efforts to uncover magneticfield induced quantum effects in monolayer TMDCs. Firstly, in[24]the weak-localization effect was observed in monolayer MoS2. Secondly, it was shown that in boron-nitride encapsulated mono- and few-layer MoS2[18]and in few layer WSe2[20]it was possible to measure the Shubnikov–de Haas(SdH)oscillations of the longitudinal resistance. Both of these developments are very significant and can provide complementary informations: the

OPEN ACCESS

RECEIVED

11 June 2015

REVISED

24 July 2015

ACCEPTED FOR PUBLICATION

4 September 2015

PUBLISHED

2 October 2015

Content from this work may be used under the terms of theCreative Commons Attribution 3.0 licence.

Any further distribution of this work must maintain attribution to the author(s)and the title of the work, journal citation and DOI.

© 2015 IOP Publishing Ltd and Deutsche Physikalische Gesellschaft

Konstanzer Online-Publikations-System (KOPS)

weak localization corrections about the coherence length and spin relaxation processes[25,26], whereas SdH oscillations about the cross-sectional area of the Fermi surface and the effective mass of the carriers.

Here wefirst briefly review the most important steps to calculate the Landau level(LL)spectrum in monolayer semiconductor TMDCs in perpendicular magneticfield using a multi-bandk p· model[3]. We show that for magneticfields ofB20 Ta simple approximation can be applied to capture all the salient features of the LL spectrum. Motivated by recent experiments in MoS2[18]and WSe2[20], we use the LL spectrum and the self-consistent Born approximation(SCBA)to calculate the SdH oscillations of the longitudinal conductances_xx.We discuss how the intrinsic spin–orbit coupling and the valley degeneracy breaking(VDB)of the magneticfield affect the magnetoconductance oscillations. We also point out the different scenarios that can occur depending on the doping level.

2. LLs in monolayer TMDCs

Electronic states in theKand-Kvalleys are related by time reversal symmetry in monolayer TMDCs and hence in the presence of a magneticfield their degeneracy should be lifted.(Note that in the case of graphene the inversion symmetry, which is present there but not in monolayer TMDCs, ensures that in the non-interacting limit the LLs remain degenerate in theKand-Kvalleys.)Recently several works have calculated the LL spectrum of monolayer TMDCs using the tight-binding(TB)method[27–29]and found that the magneticfield can indeed lift the degeneracy of the LLs in different valleys. However, due to the relatively large number of atomic orbitals that is needed to capture the zero magneticfield band structure, for certain problems, such as the SdH oscillations of longitudinal conductance, the TB methodology does not offer a convenient starting point.

On the other hand, a simplified two-bandk p· model was used to predict unconventional quantum Hall effect [30]and to discuss valley polarization[31]and magneto-optical properties[32]. This model, however, did not capture the VDB and was therefore in contradiction with the TB results and the considerations based on symmetry arguments.

Wefirst show that the VDB in perpendicular magneticfield can be described by starting from a more general, seven-bandsk p· model[3]. To this end we introduce an extended two-band continuum model which can be easily compared to previous works[30–33]. We then show a relatively simple approximation for the LL energies which will prove to be useful for the calculation of the SdH oscillations in section3.

2.1. LLs from an extended two-band model

Our starting point to discuss the magneticfield effects in monolayer TMDCs is a seven-bandk p· model (fourteen-band, if the spin degree is also taken into account), we refer the reader to[3]for details. In order to take into account the effects of a perpendicular magneticfield, one may use the Kohn–Luttinger prescription, i.e., we replace the wavenumbersq=(q_x,q_y)appearing in the seven-band model with operators:

q q ^{1 e}A,

ˆ i _

 = + whereA^T =(0,B x_z , 0)is the vector potential in Landau gauge ande>0is the magnitude of the electron charge. Note that due to this replacementqˆ₊=qˆ_x+iqˆ_yandqˆ_-=qˆ_x-iqˆ_ybecome non-commuting operators:[ ˆq_-,qˆ ]₊ = ^2eB_^z,where∣Bz∣is the strength of the magneticfield and[ ]... denotes the commutator. Working with a seven-band model is not very convenient and therefore one may want to obtain an effective model that involves fewer bands. This can be done using Löwdin-partitioning to project out those degrees of freedom from the seven-band Hamiltonian that are far from the Fermi energy. Sinceqˆ₊andqˆ_-are non-commuting operators, it is important to keep their order when one performs the Löwdin-partitioning. To illustrate this point wefirst consider a two-band model(four-band including spin)which involves the valence and the conduction bands(VB and CB). We will follow the notation used in[3]. Onefinds that the low-energy effective Hamiltonian in a perpendicular magneticfield is given by

H_eff^t,s=H₀+H_so^t,s+H_{k p}^t_·^,s, ( )1 wheres=1(s=−1)denotes spin()and

H m

q q q q

g B s

2 2

1

2 _{z z} 2

0 2

e e B

ˆ ˆ ˆ ˆ

 ( )

m

= +

+ - - + +

is the free electron term(g_e »2is theg-factor andμBis the Bohr magneton). Furthermore,

H s

s 0

0 3

s z

so, vb z

cb

t ( )

= D t

D

t ⎛

⎝⎜ ⎞

⎠⎟

describes the spin–orbit coupling in VB and CB(szis a spin Pauli matrix)andτ= ±1 for the±Kvalleys. The k p· HamiltonianH_{k p}^t_·^,sreads

H_{k p}^t_·^,s=H_D^t,s+H_as^t,s+H₃^t_w^,s+H_cub^t,s, ( )4 where

H q

q , 5a

s s

s

D, vb ,

, cb

ˆ

ˆ ( )

*

e t g

t g e

t = t t

t t -

+

⎛

⎝⎜⎜ ⎞

⎠⎟⎟

H q q

q q b

0

0 , 5

s s

s

as, ,

,

ˆ ˆ

ˆ ˆ ( )

a

= b

t t t t

t t t

+ -

- +

⎛

⎝⎜⎜ ⎞

⎠⎟⎟

H

q q

c 0

0

, 5

ws s

s

3, ,

2

, 2

( ) ( )

ˆ ˆ

* ( ) k

= k

t t t

t t

+

-

⎛

⎝

⎜⎜⎜

⎞

⎠

⎟⎟⎟

H

q q q q q q

q q q q q q d

2

0

0 . 5

s s s

s s

cub,1, ,

1

, 2

, 1

,

( ) ( )

2

ˆ ˆ ˆ ˆ ˆ ˆ

ˆ ˆ ˆ ˆ ˆ ˆ

( )

( ) ( )

( ) * ( ) *

t h h

h h

= - +

+

t t t t t

t t t t

t t t t

t t t t

+ - - - - +

+ + - - + +

⎛

⎝

⎜⎜

⎞

⎠

⎟⎟

Here the operatorqˆ_^tis defined asqˆ_^t=qˆ_x  itqˆ_y.The material specific properties are encoded in the parametersevb,e_cb(band-edge energies in the absence of SOC),γ_τ,s(coupling between the VB and the CB)and

s,

a_t, βτ,s,κτ,s,h^{( )}_t¹_,s,h^{( )}_t²_,s,which describe the effects of virtual transitions between the VB(CB)and the other bands in the seven-band model. In general, the off-diagonal material parametersγs,τ,κs,τandh_s,^{( )1}_t,h^{( )}_s,²_tare complex numbers such that for the-Kvalley(τ=−1)they are the complex conjugates of theKvalley case (τ=1). In the absence of a magneticfield, the material parameters appearing in equations(5a)–(5d)can be obtained by, e.g.,fitting the eigenvalues ofH_eff^t,sto the band structure obtained from density functional theory (DFT)calculations. We refer to[3]for the details of thisfitting procedure and for tables containing the extracted parameters for monolayer semiconductor TMDCs. Here we only mention that such afitting procedure yields real numbers which depend on the spin indexsbut do not depend explicitly on the valley indexτ.(The

parametersh^{( )}_t¹_,sandh^{( )}_t²_,scannot be obtained separately fromfitting the DFT band structure, only their sum,h_t,s can be extracted. Fortunately, as we will see below, the effect ofH_cub,1^t,s is very small in the magneticfield range we are primarily interested in.)

We note that ak p· model, similar to ours, was recently used in[29,33]to calculate the LL spectrum.There are two differences between ourk p· Hamiltonian equations(4)and the model in[29,33]. Thefirst one is that higher order terms that would correspond to ourH₃^t_w^,sandH_cub,1^t,s were not considered in[29,33]. We keep these terms in order to see more clearly the magneticfield range where the approximation discussed in section2.2is valid. The second difference can be found in ourH_as^t,s(5b)and the corresponding Hamiltonian used in[29,33]. This difference can be traced back to the way the magneticfield is taken into account in the effective models that are obtained from multi-band Hamiltonians. In[29,33]first an effective zerofield two-band model was derived and then in a second step the Luttinger-prescription was performed in this effective model. Therefore the terms which are∼q²in the zerofield case become~q qˆ ˆ_{+ -}+q qˆ ˆ_{- +}after the Luttinger-prescription. In contrast, as mentioned above, we perform the Luttinger prescription in the multi-band Hamiltonian and obtain the effective two-band modelH_eff^t,s(1)in the second step. The two approaches may lead to different results because the operatorsqˆ₊,qˆ_-do not commute and this should be taken into account in the Löwdin-partitioning which yields the effective two-band model.

The spectrum ofH_eff^t,scan be calculated numerically using harmonic oscillator eigenfunctions as basis states.

TakingBz>0 for concreteness, one can see that the operatorsaanda^†_{defined as}q a,

l 2

b

ˆ_-= q a ,

l 2

b

ˆ₊= ^† where l_B=  ( ∣e B_z∣),satisfy the bosonic commutation relation[a a, ^†]=1.(ForBz<0 one has to define

q a,

l 2

b

ˆ₊= q a

l 2

b

ˆ_-= ^†). Therefore one can calculate the matrix elements ofH_eff^t,sin a large, butfinite harmonic oscillator basis and diagonalize the resulting matrix. For a large enough number of basis states the lowest eigenvalues ofH_eff^t,swill not depend on the exact number of the basis states. Such a LL calculation is shown in figure1for MoS2and infigure2for WSe2(we have used the material parameters given in[3]). One can see that the LLs in different valleys are not degenerate and that the magnitude of the VDB is different in the VB and CB and for the lower and higher-in-energy spin-split bands. While the results in the VB are qualitatively similar for MoS2and WSe2, considering the CB, for MoS2the valley splitting of the LLs is smaller in the higher spin-split band, whereas the opposite is true for WSe2. This is a consequence of the interplay of the Zeeman term in equation(2)and other, band-structure related terms which lead to VDB.(For MoS2the valley splitting in the higher spin-split CB(purple and cyan lines)is very small for the material parameter set used in these calculations and can only be noticed for large magneticfields.)One can also observe that in the CB the lowest LL is in valleyK, whereas in the VB it is in valley-K.

Further details of the VDB, including its dependence on the parameter set that can be extracted from DFT calculations, will be discussed in section2.2. Here we point out that these results qualitatively agree with the TB calculations of[27–29], i.e., the continuum approach can reproduce all important features of multi-band TB calculations. A more quantitative comparison between our results and the TB results[27–29]is difficult, partly because the details may depend on the way how the material parameters are extracted from the DFT band structure and also because in the TB calculations the Zeeman effect was often neglected.

The LL energies can also be obtained analytically in the approximation whereH₃^t_w^,sandH_cub,1^t,s are neglected.

We will not show these analytical results here because it turns out that an even simpler approximation yields a good agreement with the numerical calculations shown infigures1and2(see section2.2)and offers a suitable starting point to develop a theory for the SdH oscillations of the longitudinal conductivity.

2.2. Approximation of the LLs spectrum

In zero magneticfield, the trigonal warping term equation(5c)and the third order term equation(5d)are important in order to understand the results of recent angle resolved photoelectron spectroscopy measurements and in order to obtain a goodfit to the DFT band structure, respectively[3]. However, as we will show for the calculation of LLs the termsH₃^t_w^,sandH_cub,1^t,s are less important. To see this one can perform another Löwdin- partitioning onH_eff^t,sto obtain effective singe-band Hamiltonians for the VB and the CB separately. Keeping only lowest order terms inBzonefinds that these single-band Hamiltonians correspond to a harmonic oscillator Hamiltonian(with different effective masses in the VB and CB and for the spin-split bands)and a term which describes a linear-in-Bzsplitting of the energies of the LLs in the two valleys. Therefore the LL spectrum can be approximated by

E n 1 g B s g B a

2 1 2

1

2 , 6

n s s s

z s

,vb, z vb,

vb ,

e B vl,vb B ( )

( ) ( )

e w m m t

= + + + +

t t t _⎜⎛ _⎟

⎝ ⎞

⎠

Figure 1.Numerically calculated LL spectrum of MoS2.(a)Thefirst few LL in the higher spin-split VB. Red lines: theKvalley(τ=1), blue lines: the-Kvalley(τ=−1). The inset shows the LLs in the lower spin-split VB.(b)Thefirst few LL in the CB. LLs both in lower spin-split band and in the higher spin-split band are shown. Red and purple lines: theKvalley, blue and cyan: the-Kvalley.

Figure 2.Numerically calculated LL spectrum of WSe₂.(a)Thefirst few LL in the higher spin-split VB. Red lines: theKvalley(τ=1), blue lines: the-Kvalley(τ=−1). The inset shows the LLs in the lower spin-split VB.(b)Thefirst few LL in the CB. LLs both in lower spin-split band and in the higher spin-split band are shown. Red and purple lines: theKvalley, blue and cyan:the-Kvalley.

E n 1 g B s g B b 2

1 2

1

2 . 6

n s s s

z s

,cb, z cb,

cb ,

e B vl,cb B ( )

( ) ( )

e w m m t

= + + + +

t t t _⎜⎛ _⎟

⎝ ⎞

⎠

Here, the following notations are introduced:n=0,1,2, ...is an integer denoting the LL index,

s s

vb cb, z

vb cb vb cb

( ) ( ) ( )

e^t =e + Dt are the band edge energies in the VB(CB)for a given spin-split bandsand

s eB

vb cb m

, z

s vb cb,

( ) ( )

( ) ( )

w^t = _t are cyclotron frequencies. In terms of the parameters appearing in equations(2)–(4), forτ=1 the effective massesm_{vb cb}^{( )s} _{( )}that enter the expression of the cyclotron frequencies are given by[3]

m m E a

2 2 , 7

s s s

2 vb1,

2 e

2 bg

∣ ∣ ( )

( ) ( )

 

a g

=⎛ + -

⎝⎜⎜ ⎞

⎠⎟⎟

m m E b

2 _s 2 s s , 7

2 cb1,

2 e

2 bg

∣ ∣ ( )

( ) ( )

 =⎛  +_b + g

⎝⎜⎜ ⎞

⎠⎟⎟

whereE_bg^{( )s} =e_cb^1,s -e^1,_vb^s.The corresponding expressions forτ=−1 can be easily found from the requirement electronic states that are connected by time reversal symmetry have the same effective mass. This means that bands corresponding to the same value of the productτshave the same effective mass. The third term in equations(6a)and(6b)comes from the free-electron term(2). The VDB is described by the last term in equations(6a),(6b)and the valleyg-factors are given by

g m

E a

4 , 8

s s s

vl,vb e 2

2 bg

∣ ∣ ( )

( )

 ^a ( )

= ⎛ + g

⎝⎜⎜ ⎞

⎠⎟⎟

g m

E b

4 . 8

s

s s

vl,cb e 2

2 bg

∣ ∣ ( )

( )

 ( )

g b

= ⎛ -

⎝⎜⎜ ⎞

⎠⎟⎟

As one can see from(8a)–(8b),gvl(s)depends on the(virtual)inter-band transition matrix elementsa_s,_βsandγ. Due to the intrinsic spin–orbit coupling, the magnitude of these matrix elements is spin-dependent[3]. Note, thatgvlis different in the VB and the CB. This is in agreement with numerical calculations based on multi-band TB models[27,29]. For the CB, the details of the derivation that leads to(6b)can be found in[34], for the VB the derivation of(6a)is analogous and therefore it will not be detailed here. We note that in variance to[34], we do not define separately an out-of-plane sping-factor and a spin independent valleyg-factor, these twog-factors are merged ingvl(s). The response to magneticfield also depends on the free electron Zeeman term. The spin-indexs to be used in the evaluation of the Zeeman term in equations(6a)–(6b)follows the spin-polarization of the given spin-split band. For MoS2, the spin polarizationssof each band are shown infigure5, other MoX2(X={S, Se, Te})monolayer TMDCs have the same polarization. For monolayer WX2TMDCs the spin polarization in the VB is the same as for the MoX2, but in the CB the polarization of the lower(higher)spin-split band is the opposite [3]. We are mainly interested in how the magneticfield breaks the degeneracy of those electronic states which are connected by time reversal in the absence of the magneticfield. Using equations(6a)–(6b), the valley splitting

E_{cb vb}^{( )i}_{( )} g_{eff,cb vb}^{( )i} _{( ) B}B_z

d = m of these states can be characterized by an effectiveg-factor

g_{eff,cb vb}^{( )i} _{( )} =(g s_e t+g_{vl,cb vb}^{( )s} _{( )}),wherei=1(2)denotes the higher-in-energy(lower-in-energy)spin-split band.

In the VB the upper index(1)[(2)]is equivalent to(),but in the CB the relation depends on the specific material being considered because the polarization is different for MoX2and WX2materials. Takingfirst the MoX2monolayers onefinds that(see alsofigure5)

g_eff,vb^{( )1 = - +}

(

g_e g_vl,vb^

)

g_eff,vb^{( )2 =}

(

g_e⁺g_vl,vb^

)

, ⁽9a⁾ g_eff,cb^{( )1 =}

(

g_e ⁺g_vl,cb^

)

g_eff,cb^{( )2 = - +}

(

g_e g_vl,cb^

)

. ⁽9b⁾ For WX2monolayersgeff,vbⁱ can also be calculated by(9a), whereas in the CB

g_eff,cb^{( )1 = - +}

(

g_e g_vl,cb^

)

g_eff,cb^{( )2 =}

(

g_e ⁺g_vl,cb^

)

. ⁽10⁾ As an example the numerical values of the variousg-factors defined above are given in table1for MoS2and in table2for WSe2. One can see thatgvl,cb(s) (vb)can be comparable in magnitude toge. This explains why the valley splitting is very small for MoS2in the case of the upper spin-split band in the CB(seefigure1), whereas the opposite is true for WSe2(figure2).

As one can see from equations(8a)and(8b),gvl(s)depends explicitly on the band-gapEbg(s)of a given spins. In addition, the parametersγ,αsandβsimplicitly also depend onEbg(s)due to thefitting procedure that is used to obtain them from DFT band structure calculations[3]. It is known thatE_bg^{( )s} is underestimated in DFT calculations and its exact value at the moment is not known for most monolayer TMDCs. Therefore in[3]we have obtained two sets of thek p· band structure parameters, thefirst one usingEbg(s)from DFT and the second

one usingEbg(s)extracted from GW calculations. The calculations shown infigures1and2were obtained with the former parameter set. As shown in table1, the calculatedg-factors depend quite significantly on the choice of the parameter set. While there is an uncertainty regarding the magnitude ofg_vl^{( )s},we expect that theg-factors obtained by using the DFT and the GW parameter sets will bracket the actual experimental values. On the other hand, the effective masses are probably captured quite well by DFT calculations and therefore thefirst term in equations(6a)–(6b)is less affected by the uncertainties of the band structure parameters. The calculations in figures1and2correspond to the‘DFT’parameter set in tables1and2.

In order to see the accuracy of the approximation introduced in equation(7a)–(7b), infigure3we compare the LL spectrum obtained in this approximation and calculated numerically using the Hamiltonian(1). As one can see the approximation is very good both in the VB and in the CB up to magneticfields20T.For larger magneticfields and large LL indices(n>7)deviations start to appear between the full quantum results and the approximation. The deviations are stronger in the VB which we attribute to the larger trigonal warping[3]of the band structure in the VB. To our knowledge the effects of the non-parabolicity of the band-dispersion on the LL spectrum has not been discussed before for monolayer TMDCs.

Given the noticeable uncertainty regarding the exact values of the effectiveg-factors, one may ask which features of the LL spectrum are affected or remain qualitatively the same. Looking at tables1and2, one can see that in some cases only the magnitude of an effectiveg-factor changes, in other cases both the magnitude and the sign. Firstly, we consider a case which illustrates possible effects of the uncertainty in the magnitude of an effectiveg-factor. Infigure4we show the LLs in the lower spin-split CB in MoS2for the two differentg_eff,cb^{( )2} given in table1. One can see that infigure4(a)the VDB is small, except for the lowest LL, which is clearly separated from the other LLs. If one assumes that the LLs acquire afinite broadening then all LLs would appear as doubly

Table 1.Valleyg-factors in MOS2. in thefirst row theg-factors are obtained with the help of DFT band gap, in the second row the g-factors are calculated with a band gap taken from theGWcalculations.

E_bg^ E_bg^ g_vl,vb^ g_vl,vb^ g_eff,vb^{( )1} g_eff,vb^{( )2} g_vl,cb^ g_vl,cb^ g_eff,cb^{( )1} g_eff,cb^{( )2}

DFT 1.66 eV^a 1.838 eV^a 0.98 0.96 -1.02 2.96 -2.11 -2.05 -0.05 -4.11

GW 2.8 eV^b 2.978 eV^b 2.57 2.38 0.57 4.38 -0.52 -0.6 1.4 -2.52

aAdapted from[3].

bAdapted from[35].

Table 2.Valleyg-factors in WSe₂. In thefirst row theg-factors are obtained with the help of DFT band gap, in the second row theg-factors are calculated with a band gap taken from theGWcalculations.

E_bg^ E_bg^ g_vl,vb^ g_vl,vb^ g_eff,vb^{( )1} g_eff,vb^{( )2} g_vl,cb^ g_vl,cb^ g_eff,cb^{( )1} g_eff,cb^{( )2}

DFT 1.337 eV^a 1.766eV^a -0.38 -0.23 -2.38 1.77 -2.71 -2.81 -4.71 -0.81

GW 2.457 eV^b 2.886eV^b 2.55 1.9 0.55 3.9 -0.67 0.13 -2.67 2.13

aAdapted from[3].

bAdapted from[36].

Figure 3.Comparison of the LL spectrum in MoS₂obtained from the two-band model and from the single band model.(a)The numerically calculated LLs using(1)for theτ=1,s=−1 in the VB(squares)and the approximation(6a) (solid lines)for LL indices n=0 ... 9.₍b)The same as in(a)but for the for theτ=1,s=−1 band in the CB(squares)and the approximation(6b) (solid lines).

degenerate except the lowest one in, e.g. an STM measurement. In contrast, the LLs are infigure4(b)are more evenly spaced and may appear as non-degenerate even if they are broadened.

Secondly, in some cases also the sign ofgeffchanges depending on which parameter set is used. Forgeff>0 the LLs in theKvalley have higher energy than the LLs in the-Kvalley, while for negativegeffthe opposite is true. We note that in[37]equations(6a)–(6b)were used to understand the VDB in the excitonic transitions in MoSe2. Theexciton valley g-factor gvl,excwas obtained by considering the energy difference between the lowermost LL in the CB and the uppermost LL in the VB in each valley:

gex,vl BBz En 0,cb1, En En En . 11

0,vb1,

0,cb1,

0,vb1,

( ) ( )

^{( )}

m = ^t₌^{= } - ^t₌^{= } - ^t₌^{=- }- ^t₌^{=- }

Using equations(7a)–(8b), one can easily show that in this approximation the exciton valley g-factor is independent of the band gap and it can be expressed in terms of the effective masses in the CB and VB[37,38]:

g m

m

m m

4 2 _s . 12

ex,vl s

e cb

e vb

( )

= - -

⎛

⎝

⎜⎜

⎞

⎠

⎟⎟

Therefore, albeit the effectiveg-factors in the CB and VB separately are affected by uncertainties, the excitong- factor, in principle, can be calculated more precisely so long the effective masses are captured accurately by DFT calculations. The comparison of DFT results and ARPES measurements[3]suggest that the DFT effective masses in the VB match the experimental results quite well. At the moment, however, it is unclear how accurate are the DFT effective masses in the CB.

Finally, we make the following brief comments.

(i) In the gapped-graphene approximation, i.e., if one neglects the free electron term and the terms~a b_s, _sin equations(7a)–(7b)and in(8a)–(8b)then the lowest LL in the CB and the highest one in the VB will be non- degenerate, but for all other LLs the valley degeneracy would not be lifted[31]due to a cancellation effect between thefirst and last terms in equations(6a)and(6b).

(ii) By measuring the valleyg-factors and the effective masses one can deduce theDiracnessof the spectrum [48], i.e., the relative importance of the off-diagonal and diagonal terms inH_D^t,s(5a)andH_as^t,s(5b),

respectively.

3. SdH oscillations of longitudinal conductivity

As we will show, the results of the section2.2provide a convenient starting point for the calculation of the SdH oscillations of the magnetoconductance. Our main motivation to consider this problem comes from the recent experimental observation of SdH oscillations in monolayer[18]and few-layer[18,20]samples. Regarding previous theoretical works on magnetotransport in TMDCs, quantum corrections to the low-field magneto- conductance were studied in[25,26]. A different approach, namely, the Adams–Holstein cyclotron-orbit migration theory[39], was used in[40]to calculate the longitudinal magnetoconductanceσxx. This theory is applicable if the cyclotron frequency is much larger than the average scattering rate1 t¯_sc.By using the effective mass obtained from DFT calculations[3]and taking the measured values of the zerofield electron mobility

e n e m

2e sc

cb

m = ¯^t and the electron densitynegiven in[18]for monolayer MoS2, a rough estimate fort¯sccan be

Figure 4.Comparison of the LL spectrum in in the lower-in-energy spin-split CB of MoS₂obtained with(a)g_{eff, cb}² =−4.11 and(b) g_eff,cb^{( )2} = -2.52.LLs in different valleys are denoted by different colours.

obtained. This shows that for magneticfieldsB15 Tthe samples are in the limit ofw t_{cb sc}¯ 1and therefore the Adams–Holstein approach cannot be used to describeσxx. Therefore we will extend the approach of Ando [42]to calculates_xxin monolayer TMDCs because it can offer a more direct comparison to existing

experimental results.

Before presenting the detailed theory of SdH oscillations we qualitatively discuss the role of the doping and the assumptions that we will use. The most likely scenarios in the VB and the CB are shown infigures5(a)and (b), respectively. Consideringfirst the CB, for electron densitiesn_e~10 cm^{13 -2}measured in[18]both the upper and lower spin-split bands would be occupied. In contrast, due to the much larger spin-splitting, for hole doped samplesEFwould typically intersect only the upper spin-split VB. Such a situation may also occur for n-doped samples in those monolayer TMDCs where the spin-splitting in the CB is much larger than in MoS2, e.g., in MoTe2or WSe2. For strong doping other extrema in the VB and CB, such as theΓandQpoints may also play a role, this will be briefly discussed at the end of this section.

We will have two main assumptions in the following. Thefirst one is that one can neglect inter-valley scattering and also intra-valley scattering between the spin-split bands. Clearly, this is a simplified model whose validity needs to be checked against experiments. One can argue that in the VB(seefigure5(a))in the absence of magnetic impurities the inter-valley scattering should be strongly suppressed because it would also require a simultaneous spin-flip. A recent scanning-tunneling experiment in monolayer WSe2[41]indeed seems to show a strong supression of inter-valley scattering. In the CB, for the case shown infigure5(b), the inter-valley scattering is not forbidden by spin selection rules. Even ifEFwas smaller, such that only one of the spin-split bands is populated in a given valley, the inter-valley scattering would not be completely suppressed because the bands are broadened by disorder which can be comparable to the spin-splitting 2Δcb(2Δcb=3 meV for MoS2

and 20–30 meV for other monolayer TMDCs.)On the other hand, the intra-valley scattering between the spin- split bands in the CB should be absent due to the specific form of the intrinsic SOC, see equation(3). We note that strictly speaking any type of perturbation which breaks the mirror symmetry of the lattice, such as a substrate or certain type of point defects(e.g., sulphur vacancies)would(locally)lead to a Rashba type SOC and hence induce intra-valley coupling between the spin-split bands. It is not known how effective is this

mechanism, in the present study we neglect it. The second assumption is that we only consider the effect of short range scatterers. This assumption is widely used in the interpretation of SdH oscillations as it facilitates to obtain analytical results[42]. We note that according to[18,24], some evidence for the presence of short range

scatterers in monolayer MoS2has indeed been recently found.While short-range scatterers can, in general, cause inter-valley scattering, on the merit of its simplicity as a minimal model we only take into account intra-valley intra-band scattering.

Using these assumptions it is straightforward to extend the theory of Ando[42]to the SdH oscillations of monolayer TMDCs. Namely, as it has been shown in section2.2, for not too large magneticfields the LLs in a given band can be described by a formula which is the same as for a simple parabolic band except that it contains a term which describes a linear-in-magneticfield valley-splitting. Then, because of the assumption that one can neglect inter-valley and intra-valley inter-band scattering, the total conductance will be the sum of the

conductances of individual bands with valley and spin indicesτ,s. This simple model allows us to focus on the effects of intrinsic SOC and valley splitting on the SdH oscillations, which is our main interest here.

Following[42], we treat impurity scattering in the SCBA and use the Kubo-formalism to calculate the longitudinal conductivityσxx(for a recent discussion see, e.g.,[43,44]). Assuming a random disorder potential V( )r with short range correlationsáV( )r V( )r¢ ñ =l d_sc (r- ¢r),the self-energyS_R^t,s = S_r^t,s + Si _i^t,sin a given

band(τ,s)does not depend on the LL indexn. It is given by the implicit equation

Figure 5.Schematics of the dispersion in the VB and in the CB around theKand-Kpoints of the band structure. The spin-split bands are denoted by red and blue lines, different colours indicate different spin-polarization. The arrows show the spin-polarization for MoS₂. For typical values of doping, the Fermi-levelE_F(denoted by a dashed line)would intersect only the upper spin-split band in the VB or both spin-split bands in the CB. The index(1)and(2)denote the upper and lower spin-split band.

i 2 l E E 1

i

, 13

rs

is sc

n n s

rs s

, ,

B2

0 , ,

i,

( )

^{( )}

l

å

S + S = p

- - S + S

t t

t t t

=

¥

whereE_n^t,sis given by equations(6a)–(6b). The termlsc 2plB2on the right-hand side of equation(13)can be

rewritten as ,

l ci

2 1

2 ⁱ

sc B2

sc

( )

w ^( )

l =

p p t where1 t_sc( )ⁱ =lscm( )ⁱ 3is the scattering rate calculated in the Born- approximation in zero magneticfield. As in section2.2, the upper indexi=1(2)refers to the higher(lower)-in- energy spin-split band in a given valley(see alsofigure5).

Using the Kubo-formalism the conductivity coming from a single valley and bands_xx^t,sis calculated as

e E f E

E E

d , 14

xxs

xxs

, 2

2

( ) ,

( ) ( )



ò

s =p -¶ s

¶

t ⎛ t

⎝⎜ ⎞

⎠⎟ wheref(E)is the Fermi function and

E n 1 Re G n E G, n 1,E G n E G, n 1,E . 15

xxs

ci n

s s s s

, 2

0

A,

R,

A,

A,

( )

( )  ^{( )}

å

( ) [ ( ) ( ) ( ) ( )] ( )

s^t = w + ^{t t} + - ^{t t} +

=

¥

HereG_R^t,s(n E, )andG_A^t,s(n E, )are the retarded and advanced Green-functions, respectively. Vertex corrections are neglected in this approximation. Since we neglect inter-valley and intra-valley inter-band scattering, the disorder-averaged Green-functionG_R,A^t,s(n E, )=[E-E_n^t,s- S_R,A^t,s ]^-1is diagonal in the indicesτ,sand in the LL representation it is also diagonal in the LL indexn. The total conductivity is then given by xx s xxs

,

å

,

s = _t s^t where the summation runs over occupied subbands for a given total electron(hole)densityne(nh). In general, one has to determineS_r^t,s + Si _i^t,sby soving equation(13)numerically. The Green-functionsG_{R A}^t_,^,scan be then calculated ands_xx^t,sfollows from equation(15). It can be seen from equation(14)that at zero temperature

s E

, ( )

S^t ands_xx^t,s( )E has to be evaluated atE=EF. In the semiclassical limit, when there are many occupied LLs belowEF, i.e.,w_c( )ⁱ EF,one can derive an analytical result forσxxτ,s

, see[42,43]for the details of this calculation. Here we only give thefinal form ofσxxand compare it to the results of numerical calculations.

As mentioned above, the situation depicted infigure5(a), i.e., when there is only one occupied subband in each of the valleys is probably most relevant for p-doped samples. Onefinds that in this case the longitudinal conductance is

e E

g B

E e E

2 1

1 4

1

cos 2

2

4 1

sin 2

. 16

xx

i

z sc

sc 0

vb1 sc1 2

vb1 sc1 2 vb1

sc1 2

F vb

1 1

eff B F

vb

1 1 2

vb1 1 2

F vb1 2

vb1 sc1

vb 1 sc1

( )

( )

( )

( )

( )

( )

( )

_{( )}

( ) ( )

( ) ( )

( ) ( ) ( )

( ) ( ) ( )

( ) ( ) ( )

( ) ( )

( ) ( )





 

 

s s w t

w t w t

p w

m w t

w t

p w

= + -

+

+ +

-

-

p w t

p w t

⎡

⎣

⎢⎢

⎛

⎝⎜ ⎞

⎠⎟

⎛

⎝⎜ ⎞

⎠⎟

⎤

⎦

⎥⎥

Here ^e E ^e

m n

0 2 F 2

2sc1 2

2sc1

vb 1

( ) ( ) h

 ( )

s = ^t =

p

t is the zerofield conductivity per single valley and band,nhis the total charge density and we assumedS_r^t,s S_i^t,s E_F.The amplitudes_1,2andare given by

g m

m g m

m a

cos 2 , sin

2 ; 17

1 eff,vb

1 vb1 e

2 eff,vb

1 vb1 e

( )

( ) ( )

( ) ( )

 = ⎛p  = p

⎝⎜ ⎞

⎠⎟ ⎛

⎝⎜ ⎞

⎠⎟ k T

k T 2 b

sinh 2

, 17

2 B vb1

2 B vb1

( )

^{( )}

( ) ( )



 p w

p w

=

wherekBis the Boltzmann constant andTis the temperature. One can see that equations(16)–(17b)are very similar to the well known expression derived by Ando[42]for a two-dimensional electron gas(2DEG). The valley-splitting, which leads to the appearance of the amplitudes_1,2,plays an analogous role to the Zeeman spin-splitting in 2DEG. Therefore, under the assumption we made above, the uncertainty regarding the value of the effectiveg-factors affects the amplitude of the oscillations but not their phase. The term proportional to

B E_z

B F

m in equation(16)is usually much smaller than thefirst term. Thus, it can be neglected in the calculation of the total conductance, but may be important if one is interested only in the oscillatory part ofσxx, see below.

We emphasize that equation(16)is only accurate ifw_vb^{( )1} E_F.However, in semiconductors, especially at relatively low doping, one can reach magneticfield values where the cyclotron energy becomes comparable to EF. In this case the numerically calculatedσxxmay differ from⁴equation(16). It is known that, e.g., WSe2can be relatively easily gated into the VB, and a decent Hall mobility was recently demonstrated in few-layer samples in [15]. As a concrete example we take the following values[15]:n_h = - ´4 10 cm^{12 -2}and Hall mobility μH=700 cm²V⁻¹s⁻¹. By takingmvb1 =−0.36me[3]the Fermi energy isE_F» -26.6 meVand using that

4From a theoretical point of view, in strong magneticfields one should also calculate vertex correlations toσxx, but this is not considered here.