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

We study the Landau level spectrum using a multi-band $\mathbf{k}\cdot\mathbf{p}$ theory in monolayer transition metal dichalcogenide semiconductors. We find that in a wide magnetic field range the Landau levels 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 Landau level 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 Landau levels affect the magnetoconductance oscillations. Motivated by recent experiments 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 (Reference [18]) using the theoretical approach we have developed.


Introduction
Atomically thin transition metal dichalcogenides semiconductors (TMDCs) [1,2,3] are recognized as a material system which, due to its finite 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,7,8,9] have spurred a feverish research activity into the optical properties of these materials [10,11,12,13]. Equally influential has proved to be the fabrication of transistors based on monolayer MoS 2 [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,18,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,22,23] and was used to determine the charge density n e and to extract the Hall mobility µ H . In addition, three recent works have reported very promising progress in the efforts to uncover magnetic field induced quantum effects in monolayer TMDCs. Firstly, in Reference [24] the weak-localization effect was observed in monolayer MoS 2 . Secondly, it was shown that in boron-nitride encapsulated mono-and few layer MoS 2 [18] and in few layer WSe 2 [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 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 we first briefly review the most important steps to calculate the LL spectrum in monolayer semiconductor TMDCs in perpendicular magnetic field using a multi-band k · p model [3]. We show that for magnetic fields of B 20 T a simple approximation can be applied to capture all the salient features of the LL spectrum. Motivated by recent experiments in MoS 2 [18] and WSe 2 [20], we use the LL spectrum and the self-consistent Born approximation (SCBA) to calculate the SdH oscillations of the longitudinal conductance σ xx . We discuss how the intrinsic spin-orbit coupling and the valley degeneracy breaking (VDB) of the magnetic field affect the magnetoconductance oscillations. We also point out the different scenarios that can occur depending on the doping level.

Landau levels in monolayer TMDCs
Electronic states in the K and −K valleys are related by time reversal symmetry in monolayer TMDCs and hence in the presence of a magnetic field 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 the K and −K valleys.) Recently several works have calculated the Landau level (LL) spectrum of monolayer TMDCs using the tight-binding (TB) method [27,28,29] and found that the magnetic field 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 magnetic field 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-band k·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 valley degeneracy breaking and was therefore in contradiction with the TB results and the considerations based on symmetry arguments.
We first show that the VDB in perpendicular magnetic field can be described by starting from a more general, seven-bands k · p model [3]. To this end we introduce an extended two-band continuum model which can be easily compared to previous works [30,31,32,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 Section 3.

LLs from an extended two-band model
Our starting point to discuss the magnetic field effects in monolayer TMDCs is a sevenband k · p model (fourteen-band, if the spin degree is also taken into account), we refer the reader to Refs. [3] for details. In order to take into account the effects of a perpendicular magnetic field, one may use the Kohn-Luttinger prescription, i.e., we replace the wavenumbers q = (q x , q y ) appearing in the seven-band model with operators: is the vector potential in Landau gauge and e > 0 is the magnitude of the electron charge. Note that due to this replacement q + =q x +iq y andq − =q x −iq y become non-commuting operators: [q − ,q + ] = 2eBz , where |B z | is the strength of the magnetic field 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 we first 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 Reference [3]. One finds that the low-energy effective Hamiltonian in a perpendicular magnetic field is given by where s = 1 (s = −1) denotes spin ↑ (↓) and is the free electron term (g e ≈ 2 is the g-factor and µ B is the Bohr magneton). Furthermore, where Here the operatorq τ ± is defined asq τ ± =q x ± iτq y . The material specific properties are encoded in the parameters ε vb , ε cb (band-edge energies in the absence of SOC), γ τ,s (coupling between the VB and the CB) and α τ,s , β τ,s , κ τ,s , η τ,s , which describe the effects of virtual transitions between the VB (CB) and the other bands in the sevenband model. In general, the off-diagonal material parameters γ s,τ , κ s,τ and η (1) s,τ , η (2) s,τ are complex numbers such that for the −K valley (τ = −1) they are the complex conjugates of the K valley case (τ = 1). In the absence of a magnetic field, the material parameters appearing in Eqs. (5a) -(5d) can be obtained by, e.g., fitting the eigenvalues of H τ,s eff to the band structure obtained from density functional theory (DFT) calculations. We refer to Reference [3] for the details of this fitting procedure and for tables containing the extracted parameters for monolayer semiconductor TMDCs. Here we only mention that such a fitting procedure yields real numbers which depend on the spin index s but do not depend explicitly on the valley index τ . (The parameters η (1) τ,s and η (2) τ,s cannot be obtained separately from fitting the DFT band structure, only their sum, η τ,s can be extracted. Fortunately, as we will see below, the effect of H τ,s cub,1 is very small in the magnetic field range we are primarily interested in. ) We note that a k · p model, similar to ours, was recently used in References [29,33] to calculate the LL spectrum. There are two differences between our k · p Hamiltonian Eqs. (4) and the model in References [29,33]. The first one is that higher order terms that would correspond to our H τ,s 3w and H τ,s cub,1 were not considered in References [29,33]. We keep these terms in order to see more clearly the magnetic field range where the approximation discussed in Section 2.2 is valid. The second difference can be found in our H τ,s as (5b) and the corresponding Hamiltonian used in [29,33]. This difference can be traced back to the way the magnetic field is taken into account in the effective models that are obtained from multi-band Hamiltonians. In References [29,33] first an effective zero field two-band model was derived and then in a second step the Luttingerprescription was performed in this effective model. Therefore the terms which are ∼ q 2 in the zero field case become ∼q +q− +q −q+ after the Luttinger-prescription. In contrast, as mentioned above, we perform the Luttinger prescription in the multiband Hamiltonian and obtain the effective two-band model H τ,s eff (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 of H τ,s eff can be calculated numerically using harmonic oscillator eigenfunctions as basis states. Taking B z > 0 for concreteness, one can see that the operators a and a † defined asq − = √ 2 . Therefore one can calculate the matrix elements of H τ,s eff in a large, but finite harmonic oscillator basis and diagonalize the resulting matrix. For a large enough number of basis states the lowest eigenvalues of H τ,s eff will not depend on the exact number of the basis states. Such a LL calculation is shown in Figure 1 for MoS 2 and in Figure 2 for WSe 2 (we have used the material parameters given in Reference [3]). One can see that the LLs in different valleys are not degenerate and that the magnitude of the valley degeneracy breaking 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 MoS 2 and WSe 2 , considering the CB, for MoS 2 the valley splitting of the LLs is smaller in the higher spin-split band, whereas the opposite is true for WSe 2 . This is a consequence of the interplay of the Zeeman term in Eq. (2) and other, band-structure related terms which lead to VDB. (For MoS 2 the 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 magnetic fields.) One can also observe that in the CB the lowest LL is in valley K, 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 Section 2.2. Here we point out that these results qualitatively agree with the TB calculations of Reference [27,28,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,28,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 where H τ,s 3w and H τ,s cub,1 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 in Figures 1 and 2 (see Section 2.2) and offers a suitable starting point to develop a theory for the SdH oscillations of the longitudinal conductivity.

Approximation of the LLs spectrum
In zero magnetic field, the trigonal warping term Eq. (5c) and the third order term Eq. (5d) are important in order to understand the results of recent angle resolved photoelectron spectroscopy measurements and in order to obtain a good fit to the DFT band structure, respectively [3]. However, as we will show for the calculation of LLs the terms H τ,s 3w and H τ,s cub,1 are less important. To see this one can perform another Löwdinpartitioning on H τ,s eff to obtain effective singe-band Hamiltonians for the VB and the CB separately. Keeping only lowest order terms in B z one finds 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-B z splitting of the energies of the LLs in the two valleys. Therefore the LL spectrum can be approximated by Here, the following notations are introduced: n = 0, 1, 2, . . . is an integer denoting the LL index, ε τ,s vb(cb) = ε vb(cb) + τ ∆ vb(cb) s z are the band edge energies in the VB (CB) for a given spin-split band s and ω  (1,s) vb 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 τ s have the same effective mass. The third term in Eqs (6a) and (6b) comes from the free-electron term (2). The VDB is described by the last term in Eqs. (6a), (6b) and the valley g-factors are given by As one can see from (8a)-(8b), g vl depends on the (virtual) inter-band transition matrix elements α s , β s and γ. Due to the intrinsic spin-orbit coupling, the magnitude of these matrix elements is spin-dependent [3]. Note, that g vl is different in the VB and the CB. This is in agreement with numerical calculations based on multi-band tight-binding 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 Reference [34], we do not define separately an out-of-plane spin g-factor and a spin independent valley g-factor, these two g-factors are merged in g (s) vl . The response to magnetic field also depends on the free electron Zeeman term. The spin-index s to be used in the evaluation of the Zeeman term in Eqs. (6a)-(6b) follows the spin-polarization of the given spin-split band. For MoS 2 , the spin polarizations s of each band are shown in Figure 5, other MoX 2 (X={S, Se, Te}) monolayer TMDCs have the same polarization. For monolayer WX 2 TMDCs the spin polarization in the VB is the same as for the MoX 2 , but in the CB the polarization of the lower (higher) spin-split band is the opposite [3]. We are mainly interested in how the magnetic field breaks the degeneracy of those electronic states which are connected by time reversal in the absence of the magnetic field. Using Eqs. (6a)-(6b), the valley splitting δE ef f,cb (vb) µ B B z of these states can be characterized by an effective g-factor g vl,cb (vb) ), where i = 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 polarisation is different for MoX 2 and WX 2 materials. Taking first the MoX 2 monolayers one finds that (see also Figure 5) ef f,vb can also be calculated by (9a), whereas in the CB g As an example the numerical values of the various g-factors defined above are given in Table 1 for MoS 2 and in Table 2 for WSe 2 . One can see that g vl,cb (vb) can be comparable in magnitude to g e . This explains why the valley splitting is very small for MoS 2 in the case of the upper spin-split band in the CB (see Figure 1), whereas the opposite is true for WSe 2 ( Figure 2).
As one can see from Eqs. (8a) and (8b), g vl depends explicitly on the band-gap E (s) bg of a given spin s. In addition, the parameters γ, α s and β s implicitly also depend on E (s) bg due to the fitting procedure that is used to obtain them from DFT band structure calculations [3]. It is known that E (s) bg is underestimated in DFT calculations and its exact value at the moment is not known for most monolayer TMDCs. Therefore in Reference [3] we have obtained two sets of the k · p band structure parameters, the first one using E (s) bg from DFT and the second one using E (s) bg extracted from GW calculations. The calculations shown in Figures 1 and 2 were obtained with the former parameter set. As shown in Table 1, the calculated g-factors depend quite significantly on the choice of the parameter set. While there is an uncertainty regarding the magnitude of g (s) vl , we expect that the g-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 the first term in Eqs. (6a)-(6b) is less affected by the uncertainties of the band structure parameters. The calculations in Figures 1 and 2 correspond to the "DFT" parameter set in Tables 1 and 2.
In order to see the accuracy of the approximation introduced in Eq. (7a)-(7b), in Figure 3 we 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 magnetic fields 20T. For larger magnetic fields 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 Table 1. Valley g-factors in MoS 2 . In the first row the g-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 the GW calculations.
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 effective g-factors, one may ask which features of the LL spectrum are affected or remain qualitatively the same. Looking at Tables 1 and 2, one can see that in some cases only the magnitude of an effective g-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 effective g-factor. In Figure 4  for the two different g (2) ef f,cb given in Table 1. One can see that in Figure 4(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 a finite broadening then all LLs would appear as doubly degenerate except the lowest one in, e.g. an STM measurement. In contrast, the LLs are in Figure 4(b) are more evenly spaced and may appear as non-degenerate even if they are broadened.
Secondly, in some cases also the sign of g ef f changes depending on which parameter set is used. For g ef f > 0 the LLs in the K valley have higher energy than the LLs in the −K valley, while for negative g ef f the opposite is true. We note that in Reference [37] Eqs (6a)-(6b) were used to understand the VDB in the excitonic transitions in MoSe 2 . The exciton valley g-factor g vl,exc was obtained by considering the energy difference between the lowermost LL in the CB and the uppermost LL in the VB in each valley: Using Eqs. (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]: Therefore, albeit the effective g-factors in the CB and VB separately are affected by uncertainties, the exciton g-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 ∼ α s , β s in Eqs. (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 the first and last terms in Eqs.(6a) and (6b).
ii) By measuring the valley g-factors and the effective masses one can deduce the Diracness of the spectrum [48], i.e., the relative importance of the off-diagonal and diagonal terms in H τ,s D (5a) and H τ,s as (5b), respectively.

Shubnikov-de Haas oscillations of longitudinal conductivity
As we will show, the results of the Section 2.2 provide 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 References [25,26]. A different approach, namely, the Adams-Holstein cyclotron-orbit migration theory [39], was used in Reference [40] to calculate the longitudinal magnetoconductance σ xx . This theory is applicable if the cyclotron frequency is much larger than the average scattering rate 1/τ sc . By using the effective mass obtained from DFT calculations [3] and taking the measured values of the zero field electron mobility µ e = e 2 neτsc m cb and the electron density n e given in Reference [18] for monolayer MoS 2 , a rough estimate forτ sc can be obtained. This shows that for magnetic fields B 15 T the samples are in the limit of ω cbτsc 1 and therefore the Adams-Holstein approach cannot be used to describe σ xx . Therefore we will extend the approach of Ando [42] to calculate σ xx in 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 in Figure 5 (a) and (b), respectively. Considering first the CB, for electron densities n e ∼ 10 13 /cm 2 measured in Reference [18] both the upper and lower spin-split bands would be occupied. In contrast, due to the much larger spinsplitting, for hole doped samples E F would 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 MoS 2 , e.g., in MoTe 2 or WSe 2 . For strong doping other extrema in the VB and CB, such as the Γ and Q points 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. The first one is that one can neglect inter-valley scattering and also intra-valley scattering between the spinsplit bands. Clearly, this is a simplified model whose validity needs to be checked against experiments. One can argue that in the VB (see Figure 5(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 WSe 2 [41] indeed seems to show a strong supression of inter-valley scattering. In the CB, for the case shown in Figure 5(b), the inter-valley scattering is not forbidden by spin selection rules. Even if E F was 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 MoS 2 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 Eq. (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 References [18,24], some evidence for the presence of short range scatterers in monolayer MoS 2 has 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 Section 2.2, for not too large magnetic fields 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-magnetic field 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 Reference [42], we treat impurity scattering in the self-consistent Born approximation (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 ′ ) = λ sc δ(r − r ′ ), the self-energy Σ τ,s R = Σ τ,s r + iΣ τ,s i in a given band (τ, s) does not depend on the LL index n. It is given by the implicit equation where E τ,s n is given by Eqs (6a)-(6b). The term λ sc /2πl 2 B on the right-hand side of Eq.(13) can be rewritten as λsc is the scattering rate calculated in the Born-approximation in zero magnetic field. As in Section 2.2, the upper index i = 1(2) refers to the higher(lower)-in-energy spin-split band in a given valley (see also Figure 5). Using the Kubo-formalism the conductivity coming from a single valley and band σ τ,s xx is calculated as where f (E) is the Fermi function and Here G τ,s R (n, E) and G τ,s A (n, E) are the retarded and advanced Greens-functions, respectively. Vertex corrections are neglected in this approximation. Since we neglect inter-valley and intra-valley inter-band scattering, the disorder-averaged Greensfunction G τ,s R,A (n, E) = [E − E τ,s n − Σ τ,s R,A ] −1 is diagonal in the indices τ, s and in the LL representation it is also diagonal in the LL index n. The total conductivity is then given by σ xx = τ,s σ τ,s xx where the summation runs over occupied subbands for a given total electron (hole) density n e (n h ). In general, one has to determine Σ τ,s r + iΣ τ,s i by soving Eq. (13) numerically. The Greens-functions G τ,s R,A can be then calculated and σ τ,s xx follows from Eq. (15). It can be seen from Eq. (14) that at zero temperature Σ τ,s (E) and σ τ,s xx (E) has to be evaluated at E = E F . In the semiclassical limit, when there are many occupied LLs below E F , i.e., ω (i) c ≪ E F , one can derive an analytical result for σ τ,s xx , see References [42,43] for the details of this calculation. Here we only give the final form of σ xx and compare it to the results of numerical calculations.
As mentioned above, the situation depicted in Figure 5(a), i.e., when there is only one occupied subband in each of the valleys is probably most relevant for p-doped samples. One finds that in this case the longitudinal conductance is Here σ 0 = e 2 τ (1) n h 2 is the zero field conductivity per single valley and band, n h is the total charge density and we assumed Σ τ,s r ≪ Σ τ,s i ≪ E F . The amplitudes A 1,2 and B are given by where k B is the Boltzmann constant and T is the temperature. One can see that Eqs. (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 A 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 effective g-factors affects the amplitude of the oscillations but not their phase. The term proportional to µ B B z /E F in Eq. (16) is usually much smaller than the first 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 Eq. (16) is only accurate if ω vb ≪ E F . However, in semiconductors, especially at relatively low doping, one can reach magnetic field values where the cyclotron energy becomes comparable to E F . In this case the numerically calculated σ xx may differ from Eq. (16) ‡. It is known that, e.g., WSe 2 can be relatively easily gated into the VB, and a decent Hall mobility was recently demonstrated in fewlayer samples in Reference [15]. As a concrete example we take the following values [15]: n h = −4 * 10 12 /cm 2 and Hall mobility µ H = 700cm 2 /V s. By taking m sc ≈ 1, there are around six occupied LLs. One can expect that for B z 15 T the LL spectrum is well described by Eq(6a), however, since the number of LLs is relatively low, there might be deviations between the analytically and numerically calculated σ xx . In Figure 6(a) we show a comparison between the analytical result Eq. (16) and the numerically calculated longitudinal conductance at zero temperature. The effective g-factors g (1) ef f,vb used in these calculations are given in Table 2.
One can see that for larger magnetic fields the amplitude of the oscillations is not captured very precisely by Eq. (16) but the overall agreement with the numerical results is good. Next, in Figures 6(b) and (c) we compare the oscillatory parts σ xx,osc of the longitudinal conductivity obtained from numerical calculations and from Eq (16) using two different g ef f,vb values. In the case of the numerical calculations σ xx,osc was obtained by subtracting the smooth function 2/(1 + (ω (1) vb τ (1) sc ) 2 ) from σ xx . According to ‡ From a theoretical point of view, in strong magnetic fields one should also calculate vertex correlations to σ xx , but this is not considered here.  Figure 6. Comparison of the numerically (symbols) and analytically (solid line) calculated zero temperature conductivity for WSe 2 for the situation depicted in Figure  5. a) total conductivity, g ef f,vb = 0.55. b) and c) comparison of the oscillatory parts of σ xx . In b) g (1) ef f,vb = 0.55, while in c) g (1) ef f,vb = −2.38 (see Table 2). The figures correspond to a magnetic field range of about 5.7 − 15.7 T.
Eq. (16), the valley-splitting of the LLs and the different effective g-factors should only affect the amplitude of the oscillations. While the amplitude of the oscillations indeed depends on g (1) ef f,vb , one can see that the agreement is better for g (1) ef f,vb = 0.55 than for g (1) ef f,vb = −2.38. For the latter case the position of the conductance minimuma start to differ for large magnetic field, whereas the maximuma in σ xx agree in both figures. These calculations illustrate that Eq. (16) may not agree with the numerical results when there are only a few LLs below E F .
We now turn to the case shown in Figure 5(b) when both spin-split subbands are populated. The total conductance is given by the sum of the conductances coming from the two spin-split subbands : σ xx = σ (1) xx + σ (2) xx . Since the effective masses in the spin-split bands are, in general, different, the associated scattering times τ (1) sc and τ (2) sc calculated in the Born-approximation are also different. We defineτ sc = τ (1) sc + τ (2) sc and σ 0 = e 2τ sc 2π 2 E F , and obtain for the magnetoconductance Here In Eq. (18) we have neglected terms which are ∼ µ B B z /E F . The result shown in Eq. (18) is similar to the multiple subband occupation problem in 2DEG [45,46,47]. The valley splitting affects the amplitude of the oscillations (see Eq. (19b)), whereas the intrinsic SOC can influence the amplitude of the oscillations [see Eq. (19a)] and it also leads to a phase difference [Eq. (18)] between the oscillations coming from the two spin-split subbands.
The situation depicted in Figure 5(b) is easily attained, e.g., in the CB of monolayer MoS 2 , where DFT calculations predict that the spin-splitting is 2∆ cb = 3 meV and therefore both subbands can be populated for relatively low densities. Our choice of the parameters for the numerical calculations shown below is motivated by the recent experiment of Cui et al. [18], where SdH oscillations in mono-and few layer MoS 2 samples have been measured. We use n e = 10 13 /cm 2 and mobility µ H = 1000cm 2 /Vs. The effective masses are chosen as m cb = 0.43m e and the spin-splitting in the CB is 2∆ cb = 3 meV [3]. Using these parameters we find E F = 28.43 meV. Since the effective masses are rather similar, the scattering times calculated from µ H are close to each other: τ sc ≈ 2.6 × 10 −13 s, i.e., they are almost twice as long as in the case of WSe 2 . The oscillations in σ xx should become discernible for B z 7 T, and at B z = 10T there are ten LLs in both the lower and the upper spin-split CB in each valley. We will focus on the oscillatory part σ xx,osc = σ (1) xx,osc + σ (2) xx,osc of the conductance, since this contains information about the spin and valley splittings. As in the previous example, we first calculate σ xx and σ (2) xx for the two different sets of g-factor values given in Table 1 as a function ofω cτsc which was introduced as a dimensionless scale of the magnetic field. Hereω c = eBz m cb with m = m (1) cb m (2) cb andτ sc = τ (1) sc τ (2) sc . All calculations are at zero temperature. One can observe that due to the CB spin splitting 2∆ cb the oscillations of σ (1) xx and σ (2) xx will not be in-phase for larger magnetic field. This effect is expected to be even more important for TMDCs having larger 2∆ cb than MoS 2 and leads to more complex oscillatory features in the total conductance σ xx than in the previous example of p-doped WSe 2 where only one band in each valley contributed to the conductance. One can also observe that in Figure 7(b) additional peaks with smaller amplitude appear in σ (2) xx,osc for larger magnetic fields, while there are no such peaks in σ (2) xx,osc in Figure 7(a). The origin of this behaviour can be traced back to the different valley-splitting patterns shown in Figure 4. The valley splitting of the LLs in Figure 4(a) is small (except for the lowest LL), while in Figure 4(b) all LLs belonging to different valleys are well-separated for larger fields and this leads to the appearance of the additional, smaller amplitude peaks in σ (2) xx,osc in Figure 7(b). The comparison between the numerically calculated total oscillatory part σ xx,osc = σ (1) xx,osc + σ (2) xx,osc and the corresponding analytical result given in Eq. (18) is shown in Figures 7(c) and (d). The agreement between the two approaches (red squares) and σ (2) osc (blue circles) using g (1) ef f,cb = −0.05 and g (2) ef f,cb = −4.11. b) numerically calculated σ (1) osc (red squares) and σ (2) osc (blue circles) using g (1) ef f,cb = 1.44 and g (2) ef f,cb = −2.55. c) The total oscillatory conductance σ osc = σ is qualitatively good forω cτ sc 1. However, for larger magnetic fields whereω cτ sc 1 the amplitude of the oscillations start to differ. In this regime the oscillatory behaviour in σ xx,osc can be quite complex, influenced by both the valley splitting and also by the intrinsic SOC splitting of the bands.
We have tried to analyze the experimental results by Cui et al. [18] using the theoretical approach outlined above. To this end we have first calculated σ xx,exp (B z ) by inverting the experimentally obtained resistance matrix and normalized it by the zerofield conductance σ xx,exp (0). To simplify the ananlysis, we assumed that the effective masses are the same in the two spin-split CB: m   sc . We then fitted σ xx,exp (B z )/σ xx (0) by the function f 0 (B z ) = C + A/(1 + (µ q B z ) 2 ), where the amplitudes A, C and the quantum mobility µ q are fit parameters. This function, according to Eq. (18), should give the smooth part of the conductance. The fit was performed in the magnetic field range [4T − 15T]: for smaller fields the weaklocalization corrections might be important which are not considered in this work, while in larger magnetic field the semiclassical approximation may not be accurate. We have found that σ xx,exp can be approximated quite well by f 0 (B z ). The most important parameter that can be extracted from the fit is the quantum scattering time τ sc,q , which is obtained from τ sc,q = m cb µq e . We find that it is roughly 3.5 times shorter than the transport scattering time τ sc,tr that follows from the measured Hall mobility µ H = 1000 cm 2 /Vs. The ratio τ sc,tr /τ sc,q depends to some extend on the fitting range that is used, but typically it is τ sc,tr /τ sc,q > 2. This difference may be explained by the fact that small-angle scattering is unimportant for τ sc,tr but it can affect τ sc,q . We note that Cui et. al. [18] has also found that the τ sc,tr is larger than τ sc,q , but they have used the amplitude of the longitudinal resistance oscillations in the magnetic field range 10 − 25T to extract τ sc,q and obtained τ sc,tr /τ sc,q ≈ 1.5. The significantly shorter τ sc,q makes it difficult to analyze the magnetic oscillations in a quantitative way using Eq. (18). Namely, it implies that oscillations should be discernible for B z 15T, i.e., for magnetic fields where only a few LLs are occupied and the semiclassical approximation may not be accurate. Using f 0 (B z ) we then extracted  [18]. The measured conductance oscillations σ osc xx,exp (squares) and the fitting of the function f osc (B z ) (see Eq. (20)) using n e = 1.1 * 10 13 /cm 2 , 2∆ cb = 3 meV (blue) n e = 1.31 * 10 13 /cm 2 , 2∆ cb = 3 meV (black), and n e = 1.31 * 10 13 /cm 2 , 2∆ cb = 5 meV (purple). the oscillatory part σ xx,osc (B z ) = σ xx,exp (B z ) − f 0 (B z ) of the conductance and fitted it with the function where D 1,2 are fitting parameters. As one can see in Figure 8, the fit can qualitatively reproduce the meauserements, but the complex oscillations between 15 − 22T are not captured. We also note that a somewhat better fit can be obtained if we assume that the charge density is larger than what is deduced from the classical Hall measurements (see the black line in Figure 8) and if we choose the spin-splitting larger than the value obtained from DFT calculations (purple line). In all cases we find, however, that the fit parameters D 1 and D 2 differ quite significantly in their magnitude, which is difficult to interpret in the present theoretical framework. This might indicate that additional scattering channels, such as inter-valley scattering, would have to be taken into account for a quantitative theory.
Finally, we would briefly comment on the relevance of the other valleys in the band structure for the SdH oscillations. Regarding p-doped samples, the Γ point might, in principle, be important for MoS 2 . However, according to DFT calculations [3] and ARPES measurements [49] the effective mass at the Γ point is significantly larger than at ±K and therefore we do not expect that states at Γ would lead to additional SdH oscillations. Nevertheless, they can be important for the level broadening of the sates at ±K because scattering from ±K to Γ does not require a spin-flip [3]. In the case of other monolayer TMDCs the Γ valley is most likely too far away in energy from the top of the VB at ±K to influence the transport for realistic dopings [3]. The situation can be more complicated for n-doped samples, especially for WS 2 and WSe 2 . For these two materials the states in the six Q valleys are likely to be nearly degenerate with the states in the ±K valleys. Therefore the Q valleys might be relatively easily populated for finite n-doping and, in contrast to the Γ point, the effective mass is comparable to that in the ±K [3] valleys. Therefore they may contribute to the SdH oscillations. They would also affect the level broadening of the ±K valley states because scattering from K (−K) to three of the six Q valleys is not forbidden by spin selection rules [3]. Furthermore, we note that in the absence of a magnetic field the six Q valleys are pairwise connected by time reversal symmetry. Therefore, taking into account only the lowest-in-energy spin-split band in the Q valleys, the LLs belonging to the Q valleys will be three-fold degenerate: the magnetic field, similarly to the case of the K and −K points, would lift the six-fold valley-degeneracy. The effective valley g-factors, however, might be rather different from the ones in the ±K valleys. For n-doped monolayer MoX 2 materials the situation is probably less complicated because the Q valleys are higher in energy and are not as easily populated as for the WX 2 monolayers. For MoS 2 monolayers, therefore, one can neglect the Q valleys in first approximation.

Summary
In summary, we have studied the LL spectrum of monolayer TMDCs in a k · p theory framework. We have shown that in a wide magnetic field range the effects of the trigonal warping in the band structure are not very important for the LL spectrum. Therefore the LL spectrum can be approximated by a harmonic oscillator spectrum and a linear-in-magntic field term which describes the VDB. This approximation and the assumption that only intra-valley intra-band scattering is relevant allowed us to extend previous theoretical work on SdH oscillations to the case of monolayer TMDCs. In the semiclassical limit, where analytical calculations are possible, it is found tha the VDB affects the amplitude of the SdH oscillations, whereas the spin-splitting of the bands leads to a phase difference in the oscillatory components. Since in actual experimental situations there might be only a few occupied LL below E F , we have also performed numerical calculations for the conductance oscillations and compared them to the analytical results. As it can be expected, if there are only a few LLs populated the amplitude of the SdH oscillations obtained in the semiclassical limit does not agree very well with the results of numerical calculations. This should be taken into account in the analysis of the experimental measurements. We used our theoretical results to analyze the measured SdH oscillations of Reference [18]. It is found that the quantum scattering time relevant for the SdH oscillations is significantly shorter than the transport scattering time that can be extracted from the Hall mobility. Finally, we briefly discussed the effect of other valleys in the band structure on the SdH oscillations.

Acknowledgement
A.K. thanks Xu Cui for discussions and for sending some of their experimental data prior to publication. This work was supported by Deutsche Forschungsgemeinschaft (DFG) through SFB767 and the European Union through the Marie Curie ITN S3NANO. P. R. would like to acknowledge the support of the Hungarian Scientific Research Funds (OTKA) K108676.