Valley Zeeman effect and Landau levels in Two-Dimensional Transition Metal Dichalcogenides

We provide a unified description of the valley Zeeman effect (g-factor) and Landau levels in monolayer transition metal dichalcogenides using the Luttinger-Kohn approximation with spin-orbit coupling in the Hamiltonian. The Landau level indices are symmetric in the K and K' valleys and the Landau level energy spacings are linear in the magnetic field strength. Besides the valley, atomic and spin contributions to the g-factor, we find an additional cross term due to coupling between the valley and atomic terms. We develop an approach to compute the Berry curvature and g-factors from density functional theory including spin-orbit coupling, without the need for unoccupied states. Self-energy corrections within the static GW approximation are added in the calculation of the g-factors. Our first principles results give spin- and valley-split Landau levels, in agreement with recent optical experiments. The exciton g-factors deduced in this work are also in good agreement with experiment for the bright and dark excitons in monolayer WSe2, as well as the lowest-energy bright excitons in MoSe2-WSe2 heterobilayers with different twist angles.

Many experiments have measured the exciton Laudé g-factor, g X0 , defined as ∆E Z = g X0 µ B B z , where ∆E Z is the difference between the Zeeman exciton shifts in the two distinct K and K valleys in the presence of an out-of-plane magnetic field B z .These experiments all point to g-factors of about −4 in ML TMDs [6,11,12] and a g-factor of −9.4 for the spinforbidden dark exciton in ML WSe 2 [13][14][15][16].However, the origin of the g-factor is not fully understood [4,17].It is widely believed that the g-factor consists of three terms, namely the valley, spin and atomic orbital terms [17].Quantitative predictions of these terms can be summarized into two approaches [4].The first approach is based on a phenomenological model [18], where the atomic orbital term is computed by assuming that the valence band maximum (VBM) and conduction band minimum (CBM) are composed entirely of d orbitals with magnetic quantum numbers ±2 and 0, respectively [19].The valley term is defined to be inversely proportional to the effective mass, a formula that was derived for a two-band model [8,20], wherein by definition, the valley terms for the conduction and valence bands are equal.The second approach is based on Peierls substitution into an effective multi-band k • p Hamiltonian [21].Contrary to the phenomelogical model, the atomic term vanishes and the consideration of remote bands is believed to be important [22].Both approaches have given g-factors in reasonable agreement with experiment.Neither approach includes in its theoretical formulation the formation of LLs in the presence of the external magnetic field.Instead, LLs have been derived separately in a massive Dirac fermion model [23][24][25], obtained by Landau-Peierls substitution into a two-band Hamiltonian matrix.This model results in LL indices that are asymmetric between the K and K valleys [23][24][25], and these results have been widely used to interpret both transport and optical measurements in the literature [7,[26][27][28].To include both the LLs and Zeeman effect, some authors have assumed that the LLs computed in this model have to be shifted by the Zeeman term [28].In contrast, Niu et al. have suggested that the computed LLs do not need to be further shifted by the Zeeman term.This idea led to the interpretation that the LLs are actually symmetric with respect to the asymmetric Zeeman-shifted bands, under the assumption that only the valley part of the orbital magnetic moment contributes to the Zeeman effect [24].This confusion arises from the fact that the Zeeman effect and LLs have not been derived within the same theoretical framework.Optical probes of the valley-dependent LLs and Zeeman effect reveal fully spin-and valley-split LLs [7], whereas transport measurements suggest a LL degeneracy of two [26], consistent with the proposal in Ref. [24].
In this work, we present a unified description of the valley Zeeman effect and LLs in two-dimensional (2D) TMDs using the Luttinger-Kohn approximation [29] to treat the uniform external magnetic field as a perturbation to the periodic crystal potential.The valley bands split into LLs with a constant valley-dependent Zeeman shift.In contrast to previous predictions, the LL indices obtained from this unified description are symmetric with respect to the K and K valleys and the energy spacing between LLs is linear in B z .The orbital magnetic moment describing the Zeeman shift is shown to be equivalent to the Berrycurvature-like expression derived using a semiclassical approach [30].Deriving analytically separate components of the orbital magnetic moment using a tight-binding basis, we provide clear definitions of the valley and atomic terms of the g-factor and furthermore, uncover an additional cross term in the g-factor, which arises from a coupling between the phase winding of Bloch states and the parent atomic orbital.First principles calculations of the single band g-factors are performed using density functional theory (DFT), and within the static GW approximation, including spin-orbit coupling (SOC).These g-factors, together with our predictions for the LL spacings, result in spin-and valley-split LLs, consistent with the optical experiment by Mak et al. [7].We also deduce the exciton g-factors, which are in good agreement with experiment.
We start with the Hamiltonian of an electron in a local periodic potential with SOC, in the presence of an external magnetic field B = ∇ × A = (0, 0, B z ) : where m e is the electron mass, V (r) is the local periodic potential, µ B is the Bohr magneton, σ refers to Pauli matrices, and g s ≈ 2 is the free electron g-factor.Because s z is a good quantum number at K and K in TMDs, the spin diagonal components of Eq. ( 1) give the solutions for spin up and down states.Henceforth, the spin Zeeman term H spin is dropped from our equations, and added on again in the final result for the energy levels.
We first consider the energy levels of the Hamiltonian of ML TMD without the magnetic field.Using the k • p expansion around K, the non-degenerate energy levels E n (K + q) are given by where Here π = p + 4mec 2 (σ × ∇V ) and q is the crystal momentum.u nK are the periodic part of the Bloch eigenstates of the Hamiltonian without any magnetic field.Luttinger and Kohn [29] have shown that if the magnetic field is treated as a perturbation, then, dropping H spin , the energy levels of the Hamiltonian in Eq. ( 1) can be obtained by solving where we have replaced q in E n (K + q) by the operator q = p/ + eA/ .The operator E n (K + q) required in Eq. ( 4) can then be shown to be where the effective mass m * is defined as: with the matrix elements, Π nm =< u nK |π|u mK >, and the orbital magnetic moment m nK given by Here, we have neglected terms involving Π z nm as these are small due to the depolarization effect in 2D materials [31].The numerator from Eq. ( 3) results in a term of the form . Using [q x , qy ] = −ieB z / , Eq. ( 5) can be obtained from the three-fold rotational symmetry in TMDs, which gives Π x nm = ±iΠ y nm .
Substituting Eq. ( 5) in Eq. ( 4), we obtain the equation for a free electron in the presence of B = (0, 0, B z ), with energy levels shifted by a constant given by −m nK • B. Thus, solving Eq. ( 4) and adding on H spin gives LLs with energies where the cyclotron frequency ω c = eB z /m * , the LL indices N = 0, 1, 2, ..., and the singleband g-factor is defined as g nK µ B = m z nK .Notably, our result shows that the LL indices are symmetric in K and K , in contrast to the results obtained from the massive Dirac fermion model where the LL indices are valley-dependent [23,24].The only valley-dependent term in Eq. ( 8) is −g nK µ B B z .
In contrast to previous work, our derivation results directly in an energy expression that includes both the LLs and Zeeman effects.Our final result is qualitatively similar to Niu et al.'s interpretation that the LLs are symmetric with respect to the Zeeman-shifted bands [24].However, in contrast to their result [24], we obtain spin-and valley-split LLs, as seen from the computed LLs of ML WSe 2 in Fig. 1a (details below).
We now focus on interpretation and first principles calculations of the g-factors.Using the result from perturbation theory that we find that m nK in Eq. ( 7) can be written as which has the same form as the semiclassical formula for m nK derived in Ref. [30] in the absence of SOC.We emphasize that Eq. ( 10) is obtained at non-degenerate band extrema with SOC, which is different from the total orbital magnetic moment [32].
To obtain physical insight into the origin of the orbital magnetic moment, we use a tight-binding basis where the Einstein summation convention is used.We then have Substitution into Eq.( 10) gives three terms [33], which we call valley (V), atomic (A) and cross (X) terms: where the matrix elements are defined as In the above, k is evaluated at K and K = −K.If the system has time reversal symmetry, m nK .The expressions for valley and atomic terms, depend, to a large extent, on the phase winding of the Bloch state and the parent atomic orbital angular momentum, respectively.Thus, they can be regarded as analytical definitions for the valley and atomic terms used in the literature [18].The cross term, on the other hand, has not been discussed before, suggesting that the current model of the Landé g-factor is incomplete.
A similar result can be obtained for the Berry curvature.There, the atomic term vanishes and we have: where r l l =< β l |r|β l >.Similar expressions can be found in Ref. [34] which refers to Ω (X) nk as a dipole term.
In most of the tight-binding (TB) literature, u nk is approximated as C l nk .Thus, the second term in the expression for ∂ k u nk (Eq.( 12)) is ignored, and only the valley terms are present in the computed g-factors and Berry curvatures.A direct evaluation of Eq. ( 13) is challenging, and TB results for the Berry curvature give widely varying results for the same TMD material (Table I).We therefore proceed with a first principles evaluation of the g-factors, using Eq. ( 10).
We implement our evaluation of the orbital magnetic moment using the density functional theory (DFT) code, QuantumESPRESSO [35], and the BerkeleyGW [36] package.Our u nk are obtained as the periodic part of the Kohn-Sham eigenstates of the DFT Hamiltonian.
The periodic potential V in Eq. ( 1) is taken to be the effective potential felt by an electron in the TMD material.For our DFT calculations, we use the Perdew-Berke-Ernzerhof (PBE) [37] parametrization of the exchange-correlation (xc) functional, and neglect any explicit dependence of the xc functional on the current density [38].Similar approximations have been made for the computation of nuclear magnetic resonance chemical shifts, and good agreement with experiment was obtained [39].We use optimized norm-conserving pseudopotentials [40] with an energy cutoff of 60 Ry, and our ground state charge density is obtained using a 21×21×1 k-grid.To investigate the effects of self-interactions, we extend the results of Eq. ( 10) to a non-local GW Hamiltonian in the static limit.We use the quasiparticle GW eigenvalues for E nk and replace the local Hamiltonian H, by H GW = H DF T − V xc + Σ, where Σ is the self-energy computed from the DFT PBE starting point (without SOC), and V xc represents the exchange-correlation potential present in the DFT calculation.In this way, we obtain a correction to the g-factor from DFT without SOC, and then add this correc- tion to the g-factor computed using two-component spinor wavefunctions in DFT, thereby including SOC effects.For the GW calculations, we use a cutoff of 35 Ry for the dielectric matrix and a non-uniform sampling [41] of the Brillouin Zone, starting with a 12 × 12 k-grid.
Our g-factors are unchanged when the k-grid is increased to 18 × 18.
We avoid the sum over unoccupied states, traditionally used in the Kubo formula for the Berry curvature [42], by directly evaluating ∂ k u nk as where e iθ =< u nk |u nk+dk > /| < u nk |u nk+dk > |.The term e iθ eliminates the random phase factor between u nk and u nk+dk , allowing us to use the parallel transport gauge.We have checked the numerical convergence of both the Berry curvature and the g-factor with respect to the magnitude of dk, which we have taken to be 10 −5 times the reciprocal lattice constant.Our self-consistent cycle is converged using an energy criterion of 10 −8 Ry.The Berry curvature thus computed for ML MoS 2 agrees very well with the literature [42], while the g-factors and Berry curvatures we compute using TB and DFT are quite different from one another (Table I).
We compute the single-band g-factors for ML WSe 2 at K and K , using the experimental bulk lattice constant.Using our computed effective masses of 0.32, 0.44 and −0.39 m e for the c1, c2 and v1 bands, respectively, we predict from first principles the energy levels at the valleys including both LL and Zeeman effects (Eq.( 8)), and our results are shown in Fig. 1.Both the LL diagram in Fig. 1a and the energies of c1 as a function of B z in Fig. 1b-c are in excellent agreement with the experimental results in Ref. [7] (the definitions of K and K are reversed here).In particular, our results are consistent with their conclusion [7] that the LLs are spin-and valley-polarized.In most of the literature [6,18,21], the exciton g-factors are theoretically deduced from the single-particle g-factors ignoring electron-hole (e-h) interactions.Using this assumption, we obtain GW g-factors of −4.50 and −10.32 for the excitons corresponding, respectively, to the lowest energy spin-allowed and spin-forbidden transitions in ML WSe 2 .These g-factors, as well as the g-factors obtained from DFT, are in good agreement with experiment (Table II).We also compute the excitons corresponding to the lowest energy spin-allowed transitions in an MoSe 2 -WSe 2 bilayer heterostructure, considering both AA and AA stackings, which correspond to twist angles of close to 60 • and 0 • in Ref. [45].Here, we see that GW gives better quantitative agreement with experiment than DFT (Table II).We note that hybridization between the d-orbitals of Mo and W changes the computed g-factors from those estimated separately from MoSe 2 and WSe 2 ML.The opposite signs for the exciton g-factors in the AA and AA stackings agree with experiment, and originate from the optical selection rules for circularly polarized light in the two heterostructures.Based on the sign conventions used in Ref. [45], we predict that these experimental setups should give g-factors of ∼ +4 and not −4 for ML WSe 2 .
The effects of e-h interactions have been discussed by some authors [47].Our consideration of e-h interactions will be published separately [48].Briefly, we consider the effect of frequency-dependence in the kernel used in the Bethe-Salpeter equation (BSE), and, using

FIG. 1 .
FIG. 1. Energy levels predicted at K and K in ML WSe 2 in the presence of a uniform external magnetic field.(a) Horizontal dashed lines denote the energy levels (LLs) obtained using a magnetic field of 9 T. Solid curves represent the energy bands in the absence of a magnetic field, and dashed curves represent the Zeeman-shifted bands (ignoring LL effects).GW g-factors are used for this plot; a plot using DFT g-factors looks similar.(b-c) Energy of band c1 as a function of magnetic field strength at (b) K and (c) K .Solid lines and dashed lines are computed using GW and DFT g-factors, respectively.Blue denotes spin down and red denotes spin up.

TABLE I .
Single-band g-factors and Berry curvatures (Ω in Å2 ) of ML MoS 2 computed using TB

TABLE II .
Exciton g-factors of ML WSe 2 , as well as AA -stacked and AA-stacked MoSe 2 -WSe 2 bilayers.The numbers in brackets include the effects of a frequency-dependent BSE kernel.X0 and D0 refer to the lowest-energy bright (spin-allowed) and dark (spin-forbidden) excitions, respectively.