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

Atomically thin transition metal dichalcogenides semiconductors (TMDCs) [1–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–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 MoS₂ [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 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 [24] the weak-localization effect was observed in monolayer MoS₂. Secondly, it was shown that in boron-nitride encapsulated mono- and few-layer MoS₂ [18] and in few layer WSe₂ [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 Landau level (LL) spectrum in monolayer semiconductor TMDCs in perpendicular magnetic field using a multi-band ${\bf{k}}\cdot {\bf{p}}$ model [3]. We show that for magnetic fields of $B\lesssim 20\;{\rm{T}}$ a simple approximation can be applied to capture all the salient features of the LL spectrum. Motivated by recent experiments in MoS₂ [18] and WSe₂ [20], we use the LL spectrum and the self-consistent Born approximation (SCBA) to calculate the SdH oscillations of the longitudinal conductance ${\sigma }_{{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.

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 LL spectrum of monolayer TMDCs using the tight-binding (TB) method [27–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 ${\bf{k}}\cdot {\bf{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.

We first show that the VDB in perpendicular magnetic field can be described by starting from a more general, seven-bands ${\bf{k}}\cdot {\bf{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 section 3.

2.1. LLs from an extended two-band model

Our starting point to discuss the magnetic field effects in monolayer TMDCs is a seven-band ${\bf{k}}\cdot {\bf{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 magnetic field, one may use the Kohn–Luttinger prescription, i.e., we replace the wavenumbers ${\bf{q}}=({q}_{x},{q}_{y})$ appearing in the seven-band model with operators: ${\bf{q}}\to \hat{{\bf{q}}}=\frac{1}{{\rm{i}}}{\boldsymbol{\nabla }}+\frac{e}{{\hslash }}{\bf{A}},$ where ${{\bf{A}}}^{T}=(0,{B}_{z}x,0)$ is the vector potential in Landau gauge and $e\gt 0$ is the magnitude of the electron charge. Note that due to this replacement ${\hat{q}}_{+}={\hat{q}}_{x}+{\rm{i}}{\hat{q}}_{y}$ and ${\hat{q}}_{-}={\hat{q}}_{x}-{\rm{i}}{\hat{q}}_{y}$ become non-commuting operators: $[{\hat{q}}_{-},{\hat{q}}_{+}]=\frac{2{{eB}}_{z}}{{\hslash }},$ 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. Since ${\hat{q}}_{+}$ and ${\hat{q}}_{-}$ 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 [3]. One finds that the low-energy effective Hamiltonian in a perpendicular magnetic field is given by

$$\begin{eqnarray}&&{H}_{{\rm{eff}}}^{\tau ,s}={H}_{0}+{H}_{{\rm{so}}}^{\tau ,s}+{H}_{{\bf{k}}\cdot {\bf{p}}}^{\tau ,s},\end{eqnarray} \tag{ 1 }$$

where s = 1 (s = −1) denotes spin $\uparrow$ $(\downarrow )$ and

$$\begin{eqnarray}&&{H}_{0}=\displaystyle \frac{{{\hslash }}^{2}}{2{m}_{{\rm{e}}}}\displaystyle \frac{{\hat{q}}_{+}{\hat{q}}_{-}+{\hat{q}}_{-}{\hat{q}}_{+}}{2}+\displaystyle \frac{1}{2}{g}_{{\rm{e}}}{\mu }_{{\rm{B}}}{B}_{z}{s}_{z}\end{eqnarray} \tag{ 2 }$$

is the free electron term (${g}_{{\rm{e}}}\approx 2$ is the g-factor and μ_B is the Bohr magneton). Furthermore,

$$\begin{eqnarray}{H}_{{\rm{so}}}^{\tau ,s}=\left(\begin{array}{cc}\tau {{\rm{\Delta }}}_{{\rm{vb}}}{s}_{z} & 0\\ 0 & \tau {{\rm{\Delta }}}_{{\rm{cb}}}{s}_{z}\end{array}\right)\end{eqnarray} \tag{ 3 }$$

describes the spin–orbit coupling in VB and CB (s_z is a spin Pauli matrix) and τ = ±1 for the ±K valleys. The ${\bf{k}}\cdot {\bf{p}}$ Hamiltonian ${H}_{{\bf{k}}\cdot {\bf{p}}}^{\tau ,s}$ reads

$$\begin{eqnarray}&&{H}_{{\bf{k}}\cdot {\bf{p}}}^{\tau ,s}={H}_{{\rm{D}}}^{\tau ,s}+{H}_{{\rm{as}}}^{\tau ,s}+{H}_{3w}^{\tau ,s}+{H}_{{\rm{cub}}}^{\tau ,s},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}{H}_{{\rm{D}}}^{\tau ,s}=\left(\begin{array}{cc}{\varepsilon }_{{\rm{vb}}} & \tau \;{\gamma }_{\tau ,s}{\hat{q}}_{-}^{\tau }\\ \tau \;{\gamma }_{\tau ,s}^{*}{\hat{q}}_{+}^{\tau } & {\varepsilon }_{{\rm{cb}}}\end{array}\right),\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}{H}_{{\rm{as}}}^{\tau ,s}=\left(\begin{array}{cc}{\alpha }_{\tau ,s}{\hat{q}}_{+}^{\tau }{\hat{q}}_{-}^{\tau } & 0\\ 0 & {\beta }_{\tau ,s}^{}{\hat{q}}_{-}^{\tau }{\hat{q}}_{+}^{\tau }\end{array}\right),\end{eqnarray} \tag{ 5b }$$

$$\begin{eqnarray}{H}_{3w}^{\tau ,s}=\left(\begin{array}{cc}0 & {\kappa }_{\tau ,s}{({\hat{q}}_{+}^{\tau })}^{2}\\ {\kappa }_{\tau ,s}^{*}{({\hat{q}}_{-}^{\tau })}^{2} & 0\end{array}\right),\end{eqnarray} \tag{ 5c }$$

$$\begin{eqnarray}{H}_{\mathrm{cub},1}^{\tau ,s}=-\displaystyle \frac{\tau }{2}\left(\begin{array}{cc}0 & {\eta }_{\tau ,s}^{(1)}{\hat{q}}_{+}^{\tau }{\hat{q}}_{-}^{\tau }{\hat{q}}_{-}^{\tau }+{\eta }_{\tau ,s}^{(2)}{\hat{q}}_{-}^{\tau }{\hat{q}}_{-}^{\tau }{\hat{q}}_{+}^{\tau }\\ {({\eta }_{\tau ,s}^{(1)})}^{*}{\hat{q}}_{+}^{\tau }{\hat{q}}_{+}^{\tau }{\hat{q}}_{-}^{\tau }+{({\eta }_{\tau ,s}^{(2)})}^{*}{\hat{q}}_{-}^{\tau }{\hat{q}}_{+}^{\tau }{\hat{q}}_{+}^{\tau } & 0\end{array}\right).\end{eqnarray} \tag{ 5d }$$

Here the operator ${\hat{q}}_{\pm }^{\tau }$ is defined as ${\hat{q}}_{\pm }^{\tau }={\hat{q}}_{x}\pm {\rm{i}}\tau {\hat{q}}_{y}.$ The material specific properties are encoded in the parameters ${\varepsilon }_{{\rm{vb}}},$ ${\varepsilon }_{{\rm{cb}}}$ (band-edge energies in the absence of SOC), γ_{τ, s} (coupling between the VB and the CB) and ${\alpha }_{\tau ,s},$ β_{τ, s}, κ_{τ, s}, ${\eta }_{\tau ,s}^{(1)},$ ${\eta }_{\tau ,s}^{(2)},$ 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, τ} and ${\eta }_{s,\tau }^{(1)},$ ${\eta }_{s,\tau }^{(2)}$ 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 equations (5a)–(5d) can be obtained by, e.g., fitting the eigenvalues of ${H}_{{\rm{eff}}}^{\tau ,s}$ to the band structure obtained from density functional theory (DFT) calculations. We refer to [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 ${\eta }_{\tau ,s}^{(1)}$ and ${\eta }_{\tau ,s}^{(2)}$ cannot be obtained separately from fitting the DFT band structure, only their sum, ${\eta }_{\tau ,s}^{}$ can be extracted. Fortunately, as we will see below, the effect of ${H}_{\mathrm{cub},1}^{\tau ,{\text{}}s}$ is very small in the magnetic field range we are primarily interested in. )

We note that a ${\bf{k}}\cdot {\bf{p}}$ model, similar to ours, was recently used in [29, 33] to calculate the LL spectrum.There are two differences between our ${\bf{k}}\cdot {\bf{p}}$ Hamiltonian equations (4) and the model in [29, 33]. The first one is that higher order terms that would correspond to our ${H}_{3w}^{\tau ,s}$ and ${H}_{\mathrm{cub},1}^{\tau ,s}$ were not considered in [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}_{{\rm{as}}}^{\tau ,s}$ (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 [29, 33] first an effective zero field 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 zero field case become $\sim {\hat{q}}_{+}{\hat{q}}_{-}+{\hat{q}}_{-}{\hat{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 model ${H}_{{\rm{eff}}}^{\tau ,s}$ (1) in the second step. The two approaches may lead to different results because the operators ${\hat{q}}_{+},$ ${\hat{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}_{{\rm{eff}}}^{\tau ,s}$ 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}^{\dagger }$ defined as ${\hat{q}}_{-}=\frac{\sqrt{2}}{{l}_{b}}a,$ ${\hat{q}}_{+}=\frac{\sqrt{2}}{{l}_{{\rm{b}}}}{a}^{\dagger },$ where ${l}_{{\rm{B}}}=\sqrt{{\hslash }/(e| {B}_{z}| )},$ satisfy the bosonic commutation relation $[a,{a}^{\dagger }]=1.$ (For B_z < 0 one has to define ${\hat{q}}_{+}=\frac{\sqrt{2}}{{l}_{{\rm{b}}}}a,$ ${\hat{q}}_{-}=\frac{\sqrt{2}}{{l}_{{\rm{b}}}}{a}^{\dagger }$). Therefore one can calculate the matrix elements of ${H}_{{\rm{eff}}}^{\tau ,s}$ 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}_{{\rm{eff}}}^{\tau ,s}$ will not depend on the exact number of the basis states. Such a LL calculation is shown in figure 1 for MoS₂ and in figure 2 for WSe₂ (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 MoS₂ and WSe₂, considering the CB, for MoS₂ the valley splitting of the LLs is smaller in the higher spin-split band, whereas the opposite is true for WSe₂. 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 MoS₂ 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.$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 [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 where ${H}_{3w}^{\tau ,s}$ and ${H}_{\mathrm{cub},1}^{\tau ,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 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.

2.2. Approximation of the LLs spectrum

In zero magnetic field, 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 good fit to the DFT band structure, respectively [3]. However, as we will show for the calculation of LLs the terms ${H}_{3w}^{\tau ,s}$ and ${H}_{\mathrm{cub},1}^{\tau ,s}$ are less important. To see this one can perform another Löwdin-partitioning on ${H}_{{\rm{eff}}}^{\tau ,s}$ 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

$$\begin{eqnarray}&&{E}_{n,{\rm{vb}}}^{\tau ,s}={\varepsilon }_{{\rm{vb}}}^{\tau ,s}+{\hslash }{\omega }_{{\rm{vb}}}^{(\tau ,s)}\left(n+\displaystyle \frac{1}{2}\right)+\displaystyle \frac{1}{2}{g}_{{\rm{e}}}{\mu }_{{\rm{B}}}{B}_{z}s+\displaystyle \frac{1}{2}{g}_{\mathrm{vl},{\rm{vb}}}^{(s)}{\mu }_{{\rm{B}}}{B}_{z}\;\tau ,\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{E}_{n,{\rm{cb}}}^{\tau ,s}={\varepsilon }_{{\rm{cb}}}^{\tau ,s}+{\hslash }{\omega }_{{\rm{cb}}}^{(\tau ,s)}\left(n+\displaystyle \frac{1}{2}\right)+\displaystyle \frac{1}{2}{g}_{{\rm{e}}}{\mu }_{{\rm{B}}}{B}_{z}s+\displaystyle \frac{1}{2}{g}_{\mathrm{vl},{\rm{cb}}}^{(s)}{\mu }_{{\rm{B}}}{B}_{z}\;\tau .\end{eqnarray} \tag{ 6b }$$

Here, the following notations are introduced: n = 0,1,2, ...is an integer denoting the LL index, ${\varepsilon }_{{\rm{vb}}(\mathrm{cb})}^{\tau ,s}={\varepsilon }_{{\rm{vb}}(\mathrm{cb})}+\tau {{\rm{\Delta }}}_{{\rm{vb}}(\mathrm{cb})}{s}_{z}$ are the band edge energies in the VB (CB) for a given spin-split band s and ${\omega }_{\mathrm{vb}\;\ (\mathrm{cb})}^{(\tau ,s)}=\frac{{{eB}}_{z}}{{m}_{\mathrm{vb}\;\ (\mathrm{cb})}^{(\tau ,s)}}$ are cyclotron frequencies. In terms of the parameters appearing in equations (2)–(4), for τ = 1 the effective masses ${m}_{\mathrm{vb}\;\ (\mathrm{cb})}^{(s)}$ that enter the expression of the cyclotron frequencies are given by [3]

$$\begin{eqnarray}&&\displaystyle \frac{{{\hslash }}^{2}}{2{m}_{{\rm{vb}}}^{(1,s)}}=\left(\displaystyle \frac{{{\hslash }}^{2}}{2{m}_{{\rm{e}}}}+{\alpha }_{s}-\displaystyle \frac{{| \gamma | }^{2}}{{E}_{{\rm{bg}}}^{(s)}}\right),\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{\hslash }}^{2}}{2{m}_{{\rm{cb}}}^{(1,s)}}=\left(\displaystyle \frac{{{\hslash }}^{2}}{2{m}_{{\rm{e}}}}+{\beta }_{s}+\displaystyle \frac{{| \gamma | }^{2}}{{E}_{{\rm{bg}}}^{(s)}}\right),\end{eqnarray} \tag{ 7b }$$

where ${E}_{{\rm{bg}}}^{(s)}={\varepsilon }_{{\rm{cb}}}^{1,s}-{\varepsilon }_{{\rm{vb}}}^{1,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 τ s have 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 valley g-factors are given by

$$\begin{eqnarray}&&{g}_{\mathrm{vl},{\rm{vb}}}^{(s)}=4\displaystyle \frac{{m}_{{\rm{e}}}}{{{\hslash }}^{2}}\left({\alpha }_{s}+\displaystyle \frac{{| \gamma | }^{2}}{{E}_{\mathrm{bg}}^{(s)}}\right),\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&{g}_{\mathrm{vl},{\rm{cb}}}^{(s)}=4\displaystyle \frac{{m}_{{\rm{e}}}}{{{\hslash }}^{2}}\left(\displaystyle \frac{{| \gamma | }^{2}}{{E}_{\mathrm{bg}}^{(s)}}-{\beta }_{s}\right).\end{eqnarray} \tag{ 8b }$$

As one can see from (8a)–(8b), g_vl^(s) depends on the (virtual) inter-band transition matrix elements ${\alpha }_{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 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 spin g-factor and a spin independent valley g-factor, these two g-factors are merged in g_vl^(s). 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 equations (6a)–(6b) follows the spin-polarization of the given spin-split band. For MoS₂, the spin polarizations s of each band are shown in figure 5, other MoX₂ (X={S, Se, Te}) monolayer TMDCs have the same polarization. For monolayer WX₂ TMDCs the spin polarization in the VB is the same as for the MoX₂, 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 equations (6a)–(6b), the valley splitting $\delta {E}_{{\rm{cb}}(\mathrm{vb})}^{(i)}={g}_{\mathrm{eff},\mathrm{cb}\;\ (\mathrm{vb})}^{(i)}{\mu }_{{\rm{B}}}{B}_{z}$ of these states can be characterized by an effective g-factor ${g}_{\mathrm{eff},\mathrm{cb}\;(\mathrm{vb})}^{(i)}=({g}_{e}s\tau +{g}_{\mathrm{vl},\mathrm{cb}(\mathrm{vb})}^{(s)}),$ 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 $\downarrow$ $(\uparrow ),$ but in the CB the relation depends on the specific material being considered because the polarization is different for MoX₂ and WX₂ materials. Taking first the MoX₂ monolayers one finds that (see also figure 5)

$$\begin{eqnarray}&&{g}_{\mathrm{eff},{\rm{vb}}}^{(1)}=(-{g}_{{\rm{e}}}+{g}_{\mathrm{vl},{\rm{vb}}}^{\downarrow })\quad {g}_{\mathrm{eff},{\rm{vb}}}^{(2)}=({g}_{{\rm{e}}}+{g}_{\mathrm{vl},{\rm{vb}}}^{\uparrow }),\end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray}&&{g}_{\mathrm{eff},{\rm{cb}}}^{(1)}=({g}_{{\rm{e}}}+{g}_{\mathrm{vl},{\rm{cb}}}^{\uparrow })\quad {g}_{\mathrm{eff},{\rm{cb}}}^{(2)}=(-{g}_{{\rm{e}}}+{g}_{\mathrm{vl},{\rm{cb}}}^{\downarrow }).\end{eqnarray} \tag{ 9b }$$

For WX₂ monolayers gⁱ_eff,vb can also be calculated by (9a), whereas in the CB

$$\begin{eqnarray}&&{g}_{\mathrm{eff},{\rm{cb}}}^{(1)}=(-{g}_{{\rm{e}}}+{g}_{\mathrm{vl},{\rm{cb}}}^{\downarrow })\quad {g}_{\mathrm{eff},{\rm{cb}}}^{(2)}=({g}_{{\rm{e}}}+{g}_{\mathrm{vl},{\rm{cb}}}^{\uparrow }).\end{eqnarray} \tag{ 10 }$$

As an example the numerical values of the various g-factors defined above are given in table 1 for MoS₂ and in table 2 for WSe₂. One can see that g_vl,cb(vb)^(s) can be comparable in magnitude to g_e. This explains why the valley splitting is very small for MoS₂ in the case of the upper spin-split band in the CB (see figure 1), whereas the opposite is true for WSe₂ (figure 2).

Table 1.  Valley g-factors in MOS₂. 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.

	${E}_{{\rm{bg}}}^{\downarrow }$	${E}_{{\rm{bg}}}^{\uparrow }$	${g}_{\mathrm{vl},{\rm{vb}}}^{\downarrow }$	${g}_{\mathrm{vl},{\rm{vb}}}^{\uparrow }$	${g}_{\mathrm{eff},{\rm{vb}}}^{(1)}$	${g}_{\mathrm{eff},{\rm{vb}}}^{(2)}$	${g}_{\mathrm{vl},{\rm{cb}}}^{\downarrow }$	${g}_{\mathrm{vl},{\rm{cb}}}^{\uparrow }$	${g}_{\mathrm{eff},{\rm{cb}}}^{(1)}$	${g}_{\mathrm{eff},{\rm{cb}}}^{(2)}$
DFT	1.66 eVa	1.838 eVa	0.98	0.96	$-1.02$	2.96	$-2.11$	$-2.05$	$-0.05$	$-4.11$
GW	2.8 eVb	2.978 eVb	2.57	2.38	0.57	4.38	$-0.52$	$-0.6$	1.4	$-2.52$

^aAdapted from [3]. ^bAdapted from [35].

Table 2.  Valley g-factors in WSe₂. 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.

	${E}_{{\rm{bg}}}^{\downarrow }$	${E}_{{\rm{bg}}}^{\uparrow }$	${g}_{\mathrm{vl},{\rm{vb}}}^{\downarrow }$	${g}_{\mathrm{vl},{\rm{vb}}}^{\uparrow }$	${g}_{\mathrm{eff},{\rm{vb}}}^{(1)}$	${g}_{\mathrm{eff},{\rm{vb}}}^{(2)}$	${g}_{\mathrm{vl},{\rm{cb}}}^{\downarrow }$	${g}_{\mathrm{vl},{\rm{cb}}}^{\uparrow }$	${g}_{\mathrm{eff},{\rm{cb}}}^{(1)}$	${g}_{\mathrm{eff},{\rm{cb}}}^{(2)}$
DFT	1.337 eVa	1.766eVa	$-0.38$	$-0.23$	$-2.38$	1.77	$-2.71$	$-2.81$	$-4.71$	$-0.81$
GW	2.457 eVb	2.886eVb	2.55	1.9	0.55	3.9	$-0.67$	0.13	$-2.67$	2.13

^aAdapted from [3]. ^bAdapted from [36].

As one can see from equations (8a) and (8b), g_vl^(s) depends explicitly on the band-gap E_bg^(s) of a given spin s. In addition, the parameters γ, α_s and β_s implicitly also depend on E_bg^(s) due to the fitting procedure that is used to obtain them from DFT band structure calculations [3]. It is known that ${E}_{\mathrm{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 the ${\bf{k}}\cdot {\bf{p}}$ band structure parameters, the first one using E_bg^(s) from DFT and the second one using E_bg^(s) 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}_{\mathrm{vl}}^{(s)},$ 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 equations (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 equation (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 $\lesssim 20{\rm{T}}.$ 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 effects of the non-parabolicity of the band-dispersion on the LL spectrum has not been discussed before for monolayer TMDCs.

Zoom In Zoom Out Reset image size

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 we show the LLs in the lower spin-split CB in MoS₂ for the two different ${g}_{\mathrm{eff},{\rm{cb}}}^{(2)}$ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Secondly, in some cases also the sign of g_eff changes depending on which parameter set is used. For g_eff > 0 the LLs in the K valley have higher energy than the LLs in the $-K$ valley, while for negative g_eff the opposite is true. We note that in [37] equations (6a)–(6b) were used to understand the VDB in the excitonic transitions in MoSe₂. 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:

$$\begin{eqnarray}&&{g}_{\mathrm{ex},\mathrm{vl}}{\mu }_{{\rm{B}}}{B}_{z}=({E}_{n=0,{\rm{cb}}}^{\tau =1,\downarrow }-{E}_{n=0,{\rm{vb}}}^{\tau =1,\downarrow })-({E}_{n=0,{\rm{cb}}}^{\tau =-1,\uparrow }-{E}_{n=0,{\rm{vb}}}^{\tau =-1,\uparrow }).\end{eqnarray} \tag{ 11 }$$

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]:

$$\begin{eqnarray}&&{g}_{\mathrm{ex},\mathrm{vl}}=4-2\left(\displaystyle \frac{{m}_{{\rm{e}}}}{{m}_{\mathrm{cb}}^{s}}-\displaystyle \frac{{m}_{{\rm{e}}}}{| {m}_{\mathrm{vb}}^{s}| }\right).\end{eqnarray} \tag{ 12 }$$

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 $\sim {\alpha }_{s},{\beta }_{s}$ in 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 the first and last terms in equations (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}_{{\rm{D}}}^{\tau ,s}$ (5a) and ${H}_{{\rm{as}}}^{\tau ,s}$ (5b), respectively.

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 [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 rate $1/{\bar{\tau }}_{\mathrm{sc}}.$ By using the effective mass obtained from DFT calculations [3] and taking the measured values of the zero field electron mobility ${\mu }_{{\rm{e}}}=\frac{{e}^{2}{n}_{{\rm{e}}}{\bar{\tau }}_{\mathrm{sc}}}{{m}_{\mathrm{cb}}}$ and the electron density n_e given in [18] for monolayer MoS₂, a rough estimate for ${\bar{\tau }}_{\mathrm{sc}}$ can be obtained. This shows that for magnetic fields $B\lesssim 15\;{\rm{T}}$ the samples are in the limit of ${\omega }_{\mathrm{cb}}{\bar{\tau }}_{\mathrm{sc}}\lesssim 1$ and therefore the Adams–Holstein approach cannot be used to describe σ_xx. Therefore we will extend the approach of Ando [42] to calculate ${\sigma }_{{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 figures 5 (a) and (b), respectively. Considering first the CB, for electron densities ${n}_{{\rm{e}}}\sim {10}^{13}{\mathrm{cm}}^{-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 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₂, e.g., in MoTe₂ or WSe₂. 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 spin-split 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₂ [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₂ 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 MoS₂ 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 [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({\bf{r}})$ with short range correlations $\langle V({\bf{r}})V({{\bf{r}}}^{\prime })\rangle ={\lambda }_{\mathrm{sc}}\delta ({\bf{r}}-{{\bf{r}}}^{\prime }),$ the self-energy ${{\rm{\Sigma }}}_{{\rm{R}}}^{\tau ,s}={{\rm{\Sigma }}}_{r}^{\tau ,s}+{\rm{i}}{{\rm{\Sigma }}}_{i}^{\tau ,s}$ in a given band (τ, s) does not depend on the LL index n. It is given by the implicit equation

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{r}^{\tau ,s}+i{{\rm{\Sigma }}}_{i}^{\tau ,s}=\displaystyle \frac{{\lambda }_{{sc}}}{2\pi {l}_{{\rm{B}}}^{2}}\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{E-{E}_{n}^{\tau ,s}-({{\rm{\Sigma }}}_{r}^{\tau ,s}+{\rm{i}}{{\rm{\Sigma }}}_{{\rm{i}}}^{\tau ,s})},\end{eqnarray} \tag{ 13 }$$

where ${E}_{n}^{\tau ,s}$ is given by equations (6a)–(6b). The term ${\lambda }_{\mathrm{sc}}/2\pi {l}_{{\rm{B}}}^{2}$ on the right-hand side of equation (13) can be rewritten as $\frac{{\lambda }_{\mathrm{sc}}}{2\pi {l}_{{\rm{B}}}^{2}}=\frac{1}{2\pi }{\hslash }{\omega }_{c}^{(i)}\frac{{\hslash }}{{\tau }_{\mathrm{sc}}^{(i)}},$ where $1/{\tau }_{\mathrm{sc}}^{(i)}={\lambda }_{\mathrm{sc}}{m}^{(i)}/{{\hslash }}^{3}$ 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 ${\sigma }_{{xx}}^{\tau ,s}$ is calculated as

$$\begin{eqnarray}&&{\sigma }_{{xx}}^{\tau ,s}=\displaystyle \frac{{e}^{2}}{{\pi }^{2}{\hslash }}\int {\rm{d}}E\left(-\displaystyle \frac{\partial f(E)}{\partial E}\right){\sigma }_{{xx}}^{\tau ,s}(E),\end{eqnarray} \tag{ 14 }$$

where f(E) is the Fermi function and

$$\begin{eqnarray}{\sigma }_{{xx}}^{\tau ,s}(E) & = & {({\hslash }{\omega }_{c}^{(i)})}^{2}\displaystyle \sum _{n=0}^{\infty }(n+1){\rm{Re}}\left[{G}_{{\rm{A}}}^{\tau ,s}(n,E){G}_{{\rm{R}}}^{\tau ,s}(n+1,E)\right.-\left.{G}_{{\rm{A}}}^{\tau ,s}(n,E){G}_{{\rm{A}}}^{\tau ,s}(n+1,E)\right]].\end{eqnarray} \tag{ 15 }$$

Here ${G}_{{\rm{R}}}^{\tau ,s}(n,E)$ and ${G}_{{\rm{A}}}^{\tau ,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-function ${G}_{{\rm{R}},{\rm{A}}}^{\tau ,s}(n,E)={[E-{E}_{n}^{\tau ,s}-{{\rm{\Sigma }}}_{{\rm{R}},{\rm{A}}}^{\tau ,s}]}^{-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 ${\sigma }_{{xx}}={\displaystyle \sum }_{\tau ,s}{\sigma }_{{xx}}^{\tau ,s}$ where the summation runs over occupied subbands for a given total electron (hole) density n_e (n_h). In general, one has to determine ${{\rm{\Sigma }}}_{r}^{\tau ,s}+{\rm{i}}{{\rm{\Sigma }}}_{i}^{\tau ,s}$ by soving equation (13) numerically. The Green-functions ${G}_{R,A}^{\tau ,s}$ can be then calculated and ${\sigma }_{{xx}}^{\tau ,s}$ follows from equation (15). It can be seen from equation (14) that at zero temperature ${{\rm{\Sigma }}}^{\tau ,s}(E)$ and ${\sigma }_{{xx}}^{\tau ,s}(E)$ has to be evaluated at E = E_F. In the semiclassical limit, when there are many occupied LLs below E_F, i.e., ${\hslash }{\omega }_{c}^{(i)}\ll {E}_{{\rm{F}}},$ one can derive an analytical result for σ_xx^{τ, s}, see [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

$$\begin{eqnarray}{\sigma }_{{xx}}/{\sigma }_{0} & = & \displaystyle \frac{2}{1+{({\omega }_{{\rm{vb}}}^{(1)}{\tau }_{\mathrm{sc}}^{(1)})}^{2}}\left[1-\displaystyle \frac{4{({\omega }_{{\rm{vb}}}^{(1)}{\tau }_{\mathrm{sc}}^{(1)})}^{2}}{1+{({\omega }_{{\rm{vb}}}^{(1)}{\tau }_{\mathrm{sc}}^{(1)})}^{2}}{e}^{-\frac{\pi }{({\omega }_{{\rm{vb}}}^{(1)}{\tau }_{\mathrm{sc}}^{(1)})}}\mathrm{cos}\left(\displaystyle \frac{2\pi {E}_{{\rm{F}}}}{{\hslash }{\omega }_{{\rm{vb}}}^{(1)}}\right){{\mathcal{A}}}_{1}\;{\mathcal{B}}\right.\\ & & +\left.\displaystyle \frac{{g}_{\mathrm{eff}}^{(i)}}{2}\displaystyle \frac{{\mu }_{{\rm{B}}}{B}_{z}}{{E}_{{\rm{F}}}}\displaystyle \frac{4{({\omega }_{{\rm{vb}}}^{(1)}{\tau }_{{sc}}^{(1)})}^{2}}{1+{({\omega }_{{\rm{vb}}}^{(1)}{\tau }_{{sc}}^{(1)})}^{2}}{e}^{-\frac{\pi }{({\omega }_{{\rm{}}\;{\rm{vb}}}^{(1)}{\tau }_{\mathrm{sc}}^{(1)})}}\mathrm{sin}\left(\displaystyle \frac{2\pi {E}_{{\rm{F}}}}{{\hslash }{\omega }_{{\rm{vb}}}^{(1)}}\right){{\mathcal{A}}}_{2}\;{\mathcal{B}}\right].\end{eqnarray} \tag{ 16 }$$

Here ${\sigma }_{0}=\frac{{e}^{2}{\tau }_{\mathrm{sc}}^{(1)}}{2\pi {{\hslash }}^{2}}{E}_{{\rm{F}}}=\frac{{e}^{2}{\tau }_{\mathrm{sc}}^{(1)}}{{m}_{{\rm{vb}}}^{(1)}}\frac{{n}_{{\rm{h}}}}{2}$ is the zero field conductivity per single valley and band, n_h is the total charge density and we assumed ${{\rm{\Sigma }}}_{r}^{\tau ,s}\ll {{\rm{\Sigma }}}_{i}^{\tau ,s}\ll {E}_{F}.$ The amplitudes ${{\mathcal{A}}}_{\mathrm{1,2}}$ and ${{\mathcal{B}}}_{}$ are given by

$$\begin{eqnarray}&&{{\mathcal{A}}}_{1}=\mathrm{cos}\left(\displaystyle \frac{\pi }{2}{g}_{\mathrm{eff},{\rm{vb}}}^{(1)}\displaystyle \frac{{m}_{{\rm{vb}}}^{(1)}}{{m}_{{\rm{e}}}}\right),\quad \quad {{\mathcal{A}}}_{2}=\mathrm{sin}\left(\displaystyle \frac{\pi }{2}{g}_{\mathrm{eff},{\rm{vb}}}^{(1)}\displaystyle \frac{{m}_{{\rm{vb}}}^{(1)}}{{m}_{{\rm{e}}}}\right);\end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray}&&{{\mathcal{B}}}_{}=\displaystyle \frac{2{\pi }^{2}{k}_{{\rm{B}}}T/{\hslash }{\omega }_{{\rm{vb}}}^{(1)}}{\mathrm{sinh}\left(2{\pi }^{2}{k}_{{\rm{B}}}T/{\hslash }{\omega }_{{\rm{vb}}}^{(1)}\right)},\end{eqnarray} \tag{ 17b }$$

where k_B is the Boltzmann constant and T is 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 ${{\mathcal{A}}}_{\mathrm{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 ${\mu }_{{\rm{B}}}{B}_{z}/{E}_{{\rm{F}}}$ in equation (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 equation (16) is only accurate if ${\hslash }{\omega }_{{\rm{vb}}}^{(1)}\ll {E}_{{\rm{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⁴ equation (16). It is known that, e.g., WSe₂ can 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\times {10}^{12}\;{\mathrm{cm}}^{-2}$ and Hall mobility μ_H = 700 cm² V⁻¹ s⁻¹. By taking m_vb¹ = −0.36 m_e [3] the Fermi energy is ${E}_{{\rm{F}}}\approx -26.6\;\mathrm{meV}$and using that ${\mu }_{{\rm{H}}}=e{\tau }_{\mathrm{sc}}^{(1)}/{m}_{\mathrm{vb}}^{(1)}$ we obtain ${\tau }_{\mathrm{sc}}^{(1)}=1.4\times {10}^{-13}$ s. The amplitude of the oscillations should become discernible when ${\omega }_{\mathrm{vb}}^{(1)}{\tau }_{\mathrm{sc}}^{(1)}={\mu }_{{\rm{H}}}{B}_{z}\lesssim 1,$ i.e., for magnetic fields ${B}_{z}\gtrsim 10\;{\rm{T}}$, while at B_z = 14.28 T, which corresponds to ${\omega }_{\mathrm{vb}}^{(1)}{\tau }_{\mathrm{sc}}^{(1)}\approx 1,$ there are around six occupied LLs. One can expect that for ${B}_{z}\lesssim 15\;{\rm{T}}$ the LL spectrum is well described by equation (6a), however, since the number of LLs is relatively low, there might be deviations between the analytically and numerically calculated ${\sigma }_{{xx}}.$ In figure 6(a) we show a comparison between the analytical result equation (16) and the numerically calculated longitudinal conductance at zero temperature. The effective g-factors ${g}_{\mathrm{eff},{\rm{vb}}}^{(1)}$ used in these calculations are given in table 2.

Zoom In Zoom Out Reset image size

One can see that for larger magnetic fields the amplitude of the oscillations is not captured very precisely by equation (16) but the overall agreement with the numerical results is good. Next, in figures 6(b) and (c) we compare the oscillatory parts ${\sigma }_{{xx},\mathrm{osc}}$ of the longitudinal conductivity obtained from numerical calculations and from equation (16) using two different ${g}_{\mathrm{eff},{\rm{vb}}}$ values. In the case of the numerical calculations σ_xx,osc was obtained by subtracting the smooth function $2/(1+{({\omega }_{{\rm{vb}}}^{(1)}{\tau }_{\mathrm{sc}}^{(1)})}^{2})$ from σ_xx. According to equation (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}_{\mathrm{eff},{\rm{vb}}}^{(1)},$ one can see that the agreement is better for g_eff,vb¹ = 0.55 than for ${g}_{\mathrm{eff},{\rm{vb}}}^{(1)}=-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 equation (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 : ${\sigma }_{{xx}}={\sigma }_{{xx}}^{(1)}+{\sigma }_{{xx}}^{(2)}.$ Since the effective masses in the spin-split bands are, in general, different, the associated scattering times τ_sc⁽¹⁾ and ${\tau }_{\mathrm{sc}}^{(2)}$ calculated in the Born-approximation are also different. We define ${\tilde{\tau }}_{\mathrm{sc}}={\tau }_{\mathrm{sc}}^{(1)}+{\tau }_{\mathrm{sc}}^{(2)}$ and ${\tilde{\sigma }}_{0}=\frac{{e}^{2}{\tilde{\tau }}_{\mathrm{sc}}^{}}{2\pi {{\hslash }}^{2}}{E}_{{\rm{F}}},$ and obtain for the magnetoconductance

$$\begin{eqnarray}{\sigma }_{{xx}}/{\tilde{\sigma }}_{0} & = & 2{{\mathcal{C}}}^{(2)}\displaystyle \frac{1}{1+{({\omega }_{{\rm{cb}}}^{(2)}{\tau }_{\mathrm{sc}}^{(2)})}^{2}}\left[1-\displaystyle \frac{4{({\omega }_{{\rm{cb}}}^{(2)}{\tau }_{\mathrm{sc}}^{(2)})}^{2}}{1+{({\omega }_{{\rm{cb}}}^{(2)}{\tau }_{\mathrm{sc}}^{(2)})}^{2}}{e}^{-\frac{\pi }{({\omega }_{{\rm{cb}}}^{(2)}{\tau }_{\mathrm{sc}}^{(2)})}}\mathrm{cos}\left(\displaystyle \frac{2\pi {E}_{{\rm{F}}}}{{\hslash }{\omega }_{{\rm{cb}}}^{(2)}}\right){{\mathcal{A}}}_{1}^{(2)}\;{{\mathcal{B}}}^{(2)}\right]\\ & & +2{{\mathcal{C}}}^{(1)}\displaystyle \frac{1}{1+{({\omega }_{{\rm{cb}}}^{(1)}{\tau }_{\mathrm{sc}}^{(1)})}^{2}}\left[1-\displaystyle \frac{4{({\omega }_{{\rm{cb}}}^{(1)}{\tau }_{\mathrm{sc}}^{(1)})}^{2}}{1+{({\omega }_{{\rm{cb}}}^{(1)}{\tau }_{\mathrm{sc}}^{(1)})}^{2}}{e}^{-\frac{\pi }{({\omega }_{{\rm{cb}}}^{(1)}{\tau }_{\mathrm{sc}}^{(1)})}}\mathrm{cos}\left(\displaystyle \frac{2\pi ({E}_{{\rm{F}}}-2{{\rm{\Delta }}}_{\mathrm{cb}})}{{\hslash }{\omega }_{{\rm{cb}}}^{(1)}}\right){{\mathcal{A}}}_{1}^{(1)}\;{{\mathcal{B}}}^{(1)}\right].\end{eqnarray} \tag{ 18 }$$

Here

$$\begin{eqnarray}&&{{\mathcal{C}}}^{(1)}=\left(1-\displaystyle \frac{2{{\rm{\Delta }}}_{\mathrm{cb}}}{{E}_{{\rm{F}}}}\right)\displaystyle \frac{{\tau }_{\mathrm{sc}}^{(1)}}{{\tilde{\tau }}_{\mathrm{sc}}},\quad {{\mathcal{C}}}^{(2)}=\displaystyle \frac{{\tau }_{\mathrm{sc}}^{(2)}}{{\tilde{\tau }}_{\mathrm{sc}}},\end{eqnarray} \tag{ 19a }$$

$$\begin{eqnarray}&&{{\mathcal{A}}}_{1}^{(i)}=\mathrm{cos}\left(\displaystyle \frac{\pi }{2}{g}_{\mathrm{eff}}^{(i)}\displaystyle \frac{{m}_{{\rm{cb}}}^{(i)}}{{m}_{{\rm{e}}}}\right),\end{eqnarray} \tag{ 19b }$$

$$\begin{eqnarray}&&{{\mathcal{B}}}^{(i)}=\displaystyle \frac{2{\pi }^{2}{k}_{{\rm{B}}}T/{\hslash }{\omega }_{{\rm{cb}}}^{(i)}}{\mathrm{sinh}\left(2{\pi }^{2}{k}_{{\rm{B}}}T/{\hslash }{\omega }_{{\rm{cb}}}^{(i)}\right)}.\end{eqnarray} \tag{ 19c }$$

In equation (18) we have neglected terms which are ∼μ_BB_z/E_F. The result shown in equation (18) is similar to the multiple subband occupation problem in 2DEG [45–47]. The valley splitting affects the amplitude of the oscillations (see equation (19b)), whereas the intrinsic SOC can influence the amplitude of the oscillations (see equation (19a)) and it also leads to a phase difference (equation (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₂, where DFT calculations predict that the spin-splitting is $2{{\rm{\Delta }}}_{\mathrm{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₂ samples have been measured. We use n_e = 10¹³ cm⁻² and mobility μ_H = 1000 cm² V⁻¹s⁻¹. The effective masses are chosen as ${m}_{\mathrm{cb}}^{(1)}=0.46{m}_{{\rm{e}}},$ m_cb⁽²⁾ = 0.43 m_e and the spin-splitting in the CB is $2{{\rm{\Delta }}}_{\mathrm{cb}}=3\;\mathrm{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: ${\tau }_{\mathrm{sc}}^{(1)}\approx {\tau }_{\mathrm{sc}}^{(2)}\approx 2.6\times {10}^{-13}$ s, i.e., they are almost twice as long as in the case of WSe₂. The oscillations in σ_xx should become discernible for B_z ≳ 7 T, and at B_z = 10 T there are ten LLs in both the lower and the upper spin-split CB in each valley. We will focus on the oscillatory part ${\sigma }_{{xx},\mathrm{osc}}={\sigma }_{{xx},\mathrm{osc}}^{(1)}+{\sigma }_{{xx},\mathrm{osc}}^{(2)}$ of the conductance, since this contains information about the spin and valley splittings. As in the previous example, we first calculate σ_xx⁽ⁱ⁾ numerically using equations (13)–(15) and obtain ${\sigma }_{{xx},\mathrm{osc}}^{(i)}$ by subtracting the smooth function $2{{\mathcal{C}}}^{(i)}/[1+{({\omega }_{c}^{(i)}{\tau }_{\mathrm{sc}}^{(i)})}^{2}].$ We then compare these results to the oscillations that can be obtained from equation (18).

In figures 7(a) and (b) we show the numerically calculated ${\sigma }_{{xx}}^{(1)}$ and ${\sigma }_{{xx}}^{(2)}$ for the two different sets of g-factor values given in table 1 as a function of ${\bar{\omega }}_{c}{\bar{\tau }}_{\mathrm{sc}}$ which was introduced as a dimensionless scale of the magnetic field. Here $\bar{{\omega }_{{\rm{c}}}}=\frac{{{eB}}_{z}}{{\bar{m}}_{\mathrm{cb}}}$ with $\bar{m}=\sqrt{{m}_{\mathrm{cb}}^{(1)}{m}_{\mathrm{cb}}^{(2)}}$ and ${\bar{\tau }}_{\mathrm{sc}}=\sqrt{{\tau }_{\mathrm{sc}}^{(1)}{\tau }_{\mathrm{sc}}^{(2)}}.$ All calculations are at zero temperature. One can observe that due to the CB spin splitting $2{{\rm{\Delta }}}_{\mathrm{cb}}$ the oscillations of ${\sigma }_{{xx}}^{(1)}$ and ${\sigma }_{{xx}}^{(2)}$ 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₂ and leads to more complex oscillatory features in the total conductance ${\sigma }_{{xx}}^{}$ than in the previous example of p-doped WSe₂ 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 ${\sigma }_{{xx},\mathrm{osc}}^{(2)}$ for larger magnetic fields, while there are no such peaks in ${\sigma }_{{xx},\mathrm{osc}}^{(2)}$ 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 ${\sigma }_{{xx},\mathrm{osc}}^{(2)}$ in figure 7(b). The comparison between the numerically calculated total oscillatory part ${\sigma }_{{xx},\mathrm{osc}}^{}={\sigma }_{{xx},\mathrm{osc}}^{(1)}+{\sigma }_{{xx},\mathrm{osc}}^{(2)}$ and the corresponding analytical result given in equation (18) is shown in figures 7(c) and (d). The agreement between the two approaches is qualitatively good for ${\bar{\omega }}_{{\rm{c}}}{\bar{\tau }}_{\mathrm{sc}}^{}\lesssim 1.$ However, for larger magnetic fields where ${\bar{\omega }}_{{\rm{c}}}{\bar{\tau }}_{\mathrm{sc}}^{}\gtrsim 1$ the amplitude of the oscillations start to differ. In this regime the oscillatory behaviour in ${\sigma }_{{xx},\mathrm{osc}}^{}$ can be quite complex, influenced by both the valley splitting and also by the intrinsic SOC splitting of the bands.

Zoom In Zoom Out Reset image size

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 zero-field conductance σ_{xx, exp}(0). To simplify the ananlysis, we assumed that the effective masses are the same in the two spin-split CB: ${m}_{{\rm{cb}}}^{(1)}={m}_{{\rm{cb}}}^{(2)}=0.43{m}_{{\rm{e}}},$ and hence ${\tau }_{\mathrm{sc}}^{(1)}={\tau }_{\mathrm{sc}}^{(2)}.$ We then fitted ${\sigma }_{{xx},\mathrm{exp}}({B}_{z})/{\sigma }_{{xx}}(0)$ by the function ${f}_{0}({B}_{z})=C+A/(1+{({\mu }_{{\rm{q}}}{B}_{z})}^{2}),$ where the amplitudes A, C and the quantum mobility μ_q are fit parameters. This function, according to equation (18), should give the smooth part of the conductance. The fit was performed in the magnetic field range [4T–15T]: for smaller fields the weak-localization 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 ${\tau }_{\mathrm{sc},q}=\frac{{m}_{{\rm{cb}}}{\mu }_{{\rm{q}}}}{e}.$ We find that it is roughly 3.5 times shorter than the transport scattering time ${\tau }_{\mathrm{sc},\mathrm{tr}}$ that follows from the measured Hall mobility μ_H = 1000 cm² V⁻¹s⁻¹. The ratio τ_sc,tr/τ_sc,q depends to some extend on the fitting range that is used, but typically it is ${\tau }_{\mathrm{sc},\mathrm{tr}}/{\tau }_{\mathrm{sc},{\rm{q}}}\gt 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 ${\tau }_{\mathrm{sc},\mathrm{tr}}$ is larger than ${\tau }_{\mathrm{sc},{\rm{q}}},$ but they have used the amplitude of the longitudinal resistance oscillations in the magnetic field range 10–25 T to extract τ_sc,q and obtained ${\tau }_{\mathrm{sc},\mathrm{tr}}/{\tau }_{{\rm{sc}},{\rm{q}}}\approx 1.5.$ The significantly shorter ${\tau }_{{\rm{sc}},{\rm{q}}}$ makes it difficult to analyze the magnetic oscillations in a quantitative way using equation (18). Namely, it implies that oscillations should be discernible for ${B}_{z}\gtrsim 15$ T, i.e., for magnetic fields where only a few LLs are occupied and the semiclassical approximation may not be accurate. Using f₀(B_z) we then extracted the oscillatory part ${\sigma }_{{xx},\mathrm{osc}}({B}_{z})={\sigma }_{{xx},\mathrm{exp}}({B}_{z})-{f}_{0}({B}_{z})$ of the conductance and fitted it with the function

$$\begin{eqnarray}{f}_{\mathrm{osc}}({B}_{z}) & = & -A\displaystyle \frac{4{({\mu }_{{\rm{q}}}{B}_{z})}^{2}}{{(1+{({\mu }_{{\rm{q}}}{B}_{z})}^{2})}^{2}}\mathrm{exp}\left(\displaystyle \frac{-\pi }{{\mu }_{{\rm{q}}}{B}_{z}}\right)\\ & & \times \left[{D}_{1}\mathrm{cos}\left(\displaystyle \frac{2\pi {E}_{{\rm{F}}}}{{\hslash }{\omega }_{{\rm{cb}}}}\right)+{D}_{2}\mathrm{cos}\left(\displaystyle \frac{{E}_{{\rm{F}}}-2{{\rm{\Delta }}}_{\mathrm{cb}}}{{\hslash }{\omega }_{{\rm{cb}}}}\right)\right],\end{eqnarray} \tag{ 20 }$$

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 and 22 T 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₁ and D₂ 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.

Zoom In Zoom Out Reset image size

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₂. 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 $\pm 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₂ and WSe₂. 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₂ 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₂ monolayers. For MoS₂ monolayers, therefore, one can neglect the Q valleys in first approximation.

In summary, we have studied the LL spectrum of monolayer TMDCs in a ${\bf{k}}\cdot {\bf{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 populated LLs the amplitude of the SdH oscillations obtained in the semiclassical limit does not agree wery 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 [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.

AK 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. PR would like to acknowledge the support of the Hungarian Scientific Research Funds (OTKA) K108676.