Valley-Filling Instability and Critical Magnetic Field for Interaction-Enhanced Zeeman Response in Doped WSe$_2$ Monolayers

Carrier-doped transition metal dichalcogenide (TMD) monolayers are of great interest in valleytronics due to the large Zeeman response (g-factors) in these spin-valley-locked materials, arising from many-body interactions. We develop an \textit{ab initio} approach based on many-body perturbation theory to compute the interaction-enhanced g-factors in carrier-doped materials. We show that the g-factors of doped WSe$_2$ monolayers are enhanced by screened exchange interactions resulting from magnetic-field-induced changes in band occupancies. Our interaction-enhanced g-factors $g^*$ agree well with experiment. Unlike traditional valleytronic materials such as silicon, the enhancement in g-factor vanishes beyond a critical magnetic field $B_c$ achievable in standard laboratories. We identify ranges of $g^*$ for which this change in g-factor at $B_c$ leads to a valley-filling instability and Landau level alignment, which is important for the study of quantum phase transitions in doped TMDs. We further demonstrate how to tune the g-factors and optimize the valley-polarization for the valley Hall effect.

Valleytronics, the control and manipulation of the valley degree of freedom (valley pseudospin), is being actively considered as the next new paradigm for information processing.
The field of valleytronics dates back to investigations on traditional semiconductors such as silicon 1,2 , but the ability to exploit valley polarizations in these materials has been limited 3 .
A major impetus for the renaissance of valleytronics is the recent discovery that H-phase transition metal dichalcogenide (TMD) semiconductor monolayers (MLs) are excellent candidates for valleytronics applications 3,4 . The spin-valley locking effect in these MLs 5 leads to long lifetimes for spin-and valley-polarization, while individual valleys can be probed and controlled using circularly-polarized light, paving the way to use the valley pseudospin for information processing. The valley Zeeman response in TMD MLs 6-20 is also significantly larger than in traditional semiconductors [21][22][23][24][25][26][27] .
When an external magnetic field is applied normal to the TMD ML, the energies of the valleys shift in equal magnitude and opposite directions. This Zeeman effect is quantified by the orbital and spin magnetic moments, which contribute to the Landé g-factors. In TMD MLs, the intrinsic Landé g-factors are about 6 times larger [6][7][8]10 than that in silicon, where only the spin magnetic moment dominates 21,22 . The larger g-factors in TMD MLs allow for greater control in tuning the energetics of the valley pseudospins, and results in a larger valley-polarized current, which is important for observations of the valley Hall effect 3 .
Besides the Zeeman effect, an external magnetic field also results in a quantization of states to form Landau levels (LLs).
Much of the current research on valleytronics seeks to understand how to manipulate the valley pseudospins in TMDs 3,4 . It has been found that carrier doping can dramatically enhance the g-factors in TMD monolayers [14][15][16][17][18][19][20] , opening up the possibility to tune the valley pseudospin by gating in a magnetic field. This enhancement in g-factors has been attributed to many-body interactions [14][15][16][17][18][19][20]23,26 , but a quantitative understanding is lacking. To interpret the g-factor enhancement in doped TMDs, experimentalists have typically relied on existing theoretical literature dating back to the 1960s-1970s 23,28 . However, there are two shortcomings of these theoretical approaches. Firstly, they are not ab initio methods and cannot provide quantitative predictions. Secondly, these studies all focused on silicon or III-V semiconductors, which are very different in nature from the TMD monolayers. It is important to question if the experimental observations on g-factors in doped TMDs serve only to validate in a new material what was already known for silicon, or if it is possible to observe novel phenomena not known before in conventional valleytronics materials.
In this work, we develop an ab initio approach based on many body perturbation theory to compute the interaction-enhanced Landé g-factors in carrier-doped systems. We predict that the larger intrinsic g-factors in TMD MLs enable the observation of a critical magnetic field B c above which the interaction-induced enhancement in the g-factors vanishes in doped TMDs. We identify ranges of the enhanced g-factor g * enh for which the discontinuous change in g-factor at B c results in a LL alignment and valley-filling instability for B B c . Such a phenomenon has not been observed or predicted for silicon and other conventional val- highlights the potential of creating pseudo-spinors from aligned LLs for topological quantum computing applications 36 .

Interaction-enhanced g-factors
For a many-electron system described by a static mean-field Hamiltonian H |nk = E nk |nk , it has been shown that an out-of-plane magnetic field B results in the following expression for the LLs at the K valley 6 : where n is the corresponding band index, E nK is the mean-field single-particle energy at K, m * is the valley effective mass, µ B is the Bohr magneton and N = 0, 1, 2, ... is the LL index. The total intrinsic single band g-factor consists of the orbital and spin contribution, where g s is taken to be 2.0 and s z,nK is the spin quantum number. The orbital g-factor is defined using the orbital magnetic moment g orb nK µ B = m z nK 6 (see Methods).
However, the above static mean-field description does not account for the energydependent electron self-energy Σ(E) that is necessary for a many-body description of the quasiparticle (QP) energies. Σ(E) represents the change in the energy of the bare particle due to the interaction of the particle with itself via the interacting many-body system. The change in self-energy as the QP energy shifts with B leads to an effective renormalized g-factor, g * nK , defined as (E QP,1 Thus we have, Such a renormalization effect is missing in previous first principles calculations of g-factors in TMDs 6-8 . The electron self-energy in this work is computed within the GW approximation 37 (see Methods), which uses the first order term in the perturbative expansion of Σ in terms of the screened Couloumb interaction W . Table 1 shows the renormalized and intrinsic g-factors computed for undoped monolayer WSe 2 . We see that the magnitudes of the renormalized g-factors are reduced by ∼ 20% compared to the intrinsic GW g-factors, because ∂Σ(E) ∂E is in general negative 38 . The exciton g-factors deduced using the renormalized g-factors are in good agreement with experiment [39][40][41][42][43] .
In contrast to the undoped system, doped systems have partially occupied bands. Thus, if the band occupancies are also changing in response to the magnetic field, there is an additional term in dΣ(E) dE : where f is the Fermi-Dirac distribution function. The second term in Eq (4) can be simplified to give (see Methods): where n is the band index of the frontier doped band, and E F and k F are respectively the Fermi energy and Fermi wave vector. This term comes from the screened exchange contribution to the self-energy.W nk F (defined in Methods) is an effective quasi-2D screened Coulomb potential which can be evaluated completely from first principles. SinceW nk F is positive, the second term of Eq. (4) leads to an enhancement effect for the g-factor. For an ideal 2D fermion gas, this second term reduces to the term dΣ(E)/dE derived in Ref. 23 (see also Ref. 44) where the first term in Eq. (4) is ignored. We note that our numerical results forW vk F as a function of hole density differ from those computed for an ideal 2D fermion gas 1,23 , although the corresponding numerical results for the bare Coulomb potential match well with the ideal 2D case (see Supplementary Figure 1). This observation implies that an ab initio non-local description of the dielectric function of the quasi-2D system is important for a quantitative prediction of the renormalized g-factors.
The g-factors in this work are all computed for the valence band at K, and henceforth, the subscript vK is omitted. The magnitudes of our computed valence band g-factors |g * | are plotted as red squares in Fig. 1b for ML WSe 2 with different hole densities. Our predicted gfactors agree well with those deduced from multiple experiments 10,14,16,18 on hole-doped ML WSe 2 . The renormalized g-factor for the undoped system is labeled g * 0 . Due to interactioninduced enhancement, the g-factor increases significantly once hole carriers are introduced.
This enhancement reduces as the hole density is increased as expected from the density dependence of the many-body Coulomb interactions (see Supplementary Figure 1).

Critical magnetic field
A subtle but important point is that Eq. (5) applies only when the band occupancies are changing with B, and the Fermi level E F is fixed. In electrostatic gating experiments, the carrier concentration is fixed rather than the absolute Fermi level. However, the Zeeman shifts in K and K are equal in magnitude but opposite in sign (Fig. 2a, B < B c ), and the density of states for the quadratic bands in 2D is independent of energy. So for B small enough that both valleys have carriers (mixed polarized regime; Fig. 1a), E F is fixed while the band occupancies change and both terms in Eq. (4) will apply, leading to the interaction-enhanced g-factors, which we label as g * enh . However, above a critical magnetic field B c ∼ |E F |/(|g * enh |µ B ), only one valley has carriers (Fig. 2a, B > B c ) (see Supplementary Note for a more precise expression for B c ). As B increases beyond B c , a constant hole density is maintained when E F shifts with the bands without changing the band occupancies. Thus, for B > B c , only the first term in Eq. (4) applies, similar to the undoped case, leading to an abrupt drop in g * at B = B c (Fig. 2b), with a corresponding piecewise-linear Zeeman split E Z (Fig. 2c). For a hole density of 5.2 × 10 12 cm −2 , B c ≈ 36 T, and we have g * ≈ −11.0 for B < B c and g * ≈ −6.6 for B > B c .
This abrupt drop in g * at a critical magnetic field has never been reported or predicted before in traditional valleytronic materials such as silicon. Indeed, B c is inversely related to |g * enh |, and it is the large intrinsic g-factors and hence large renormalized g-factors for TMDs that allow for B c to be small enough to be reached in standard laboratories. For the same hole density of 5.2 × 10 12 cm −2 , we predict B c in silicon to be ∼ 400 T. The larger intrinsic g-factors for TMD MLs arise from the large orbital g-factors, which consist of a valley term, an orbital term, and a cross term that involves coupling between the phase-winding of the Bloch states and the parent atomic orbitals. 6 Since B c is the value of B characterizing the onset of the fully polarized regime, B c can be deduced using optical measurements of the exciton and polaron energies for K and K 16 .
In  which g * is enhanced by interactions (Fig. 3a). However, the magnitude of g * enh decreases as the hole concentration increases (Fig. 1). These competing effects imply that for any given B field B 0 , there is an optimal hole concentration ρ 0 which maximizes the Zeeman split E Z (Fig. 3b). This optimal hole concentration ρ 0 corresponds to the hole concentration for which B c = B 0 (Fig. 3), and yields a maximum concentration of valley-polarized carriers.
These predictions are useful for realizations of the valley Hall effect and other applications where a high concentration of valley-polarized carriers is desired.

LL Alignment and Valley Filling Instability
The abrupt change in g * at B = B c also has other interesting implications. As B increases beyond B c , the decrease in |g * | results in a decrease in magnitude of the slopes of the LL fan diagrams (Fig. 4a), leading to a crossing between the energies of 0K and N K (purple circle, N = 5 in Fig. 4a). If N K has carriers, this LL alignment results in a valley-filling instability, where the hole population is transferred back and forth between the two valleys for small changes in B.
LLs are filled with holes from higher energy to lower energy. In Fig. 4c, the blue and pink shading indicate the filling of the LLs with holes. If N K is fully occupied, the plot is shaded from N K + ωc 2 down to N K − ωc 2 . At B = B c , an integer number of LLs (5 in Fig. 4c) are completely filled with holes while 0K is completely empty. As B decreases slightly below B c , the LL degeneracy decreases, and the LL with the next lower energy (0K in Fig. 4c) is required to contain the holes, leading to a symmetric zigzag fine structure about E F . As B increases slightly above B c , the LL degeneracy increases. As long as 5K > 0K , only the K valley is filled with holes (B = B 1 in Fig. 4b), and the blue shaded area represents the constant hole density. But when B is slightly larger than B X (Fig. 4d) where the energies of 0K and 5K cross, holes will start to fill the K valley again (B = B 2 in Fig. 4b), until the B field is large enough that the LLs 0K to 4K can contain all the holes (yellow circle, Fig. 4d) and the system becomes fully polarized again. This represents a valley-filling instability, where K is depleted of holes from B = B c to B = B X , and filled again up to B X + ∆B. In practice, when holes begin to fill 0K , the mixed polarized regime is reached and g * becomes enhanced, leading to a change in the slope of the fan diagram that is expected to result in a LL alignment not just for B = B X but also for B up to B X + ∆B.  Table 1).
The alignment of LLs is of interest to investigate quantum phase transitions in these doped TMDs 16,[29][30][31][32][33][34] . Given that the LLs are expected to align for B between B X and B X + ∆B, it is interesting to predict how large ∆B can be and how sensitive ∆B is to fluctuations in g * enh . Not all values of g * enh will result in the instability (see Supplementary  Figure 2 and Supplementary Note). In particular, if the energies of 0K and N K cross at B X , the valley-filling instability only occurs if N K is occupied with holes. Fig. 4e plots the ranges of g * enh for which a valley-filling instability will occur, as well as the corresponding ∆B values. We see that an optimal range of |g * enh | for LL alignment is 10.4 < |g * enh | < 11.2. Here, ∆B is quite large (∼ 8T for |g * enh | ∼ 10.4) and is also fairly robust to changes in g * enh . The corresponding values of B c and B X fall within 30 to 40 T (Supplementary Table 1 g * enh can be tuned such that 0K matches exactly with N K for some N . However, once g * enh deviates slightly from this value, due to fluctuations in the hole density or dielectric environment (see Fig. 5), the LLs are no longer aligned. Our predictions above enable the alignment of LLs while allowing for some fluctuations in g * enh .

Tunability of interaction-enhanced g-factor
We further note that in addition to electrostatic gating which changes the carrier concentration and thus g * enh (Fig. 1), g * enh can also be tuned by dielectric screening (Fig. 5). The tunability of g * enh with the background dielectric constant can be understood from the fact that g * enh is related to the effective quasi-2D screened Coulomb potential at the Fermi surface (Eq. 5). This tunability of g * enh provides a handle to control the valley-polarized current, B c and ∆B.

DISCUSSION
In summary, our ab initio calculations show that many-body interactions in doped TMD MLs enhance the g-factors compared to the undoped MLs, up to a critical magnetic field

Calculation of intrinsic g-factor
The orbital component of the intrinsic g-factor g orb nK is defined as g orb nK µ B = m z nK 6 : We use the PBE exchange-correlation functional 46  The intrinsic single-band g-factor reduces by only 0.2 when an energy cutoff of 2Ry with 200 empty bands is used.

Calculation of renormalized g-factor
The renormalized g-factor g * is computed from the intrinsic g-factor g I and dΣ(E) dE using Eq. 3. As the dependence of E QP nk on B is no longer linear, we generalize Eq. (2) to the case where E QP nk refers to the quasiparticle energies in the presence of a B-field, and the applied B represents a small increment in B.
We approximate g I using the value in the undoped system. The intrinsic single band g-factors from DFT calculations do not change when ML WSe 2 is doped with holes.
where we define the quasi-2D screened Coulomb potential: Here, Ω is the cell volume, N q is the number of q-points and v q is the Coulomb potential with the slab Coulomb truncation scheme applied 50 .W mq is an effective quasi-2D screened Coulomb potential defined in valley K and L is the height of the supercell for ML WSe 2 .
The second term in Eq. (4) is then given by (see Supplemental Note): We compute g * vK by evaluatingW vk F ab initio using the random phase approximation for the dielectric matrix. Due to the partial occupancies, we calculate the dielectric matrix using a dense reciprocal space sampling of 120 × 120, a 2Ry G-vector cut off and 29 bands. g * is unchanged when we use instead 4Ry and 299 bands, and reduces by ∼ 3% when a 240 × 240 k-mesh is used. Care is taken to include the effect of spin-orbit splitting at the valleys. For the effective mass, we use our DFT value of m * = −0.48m e , which agrees well with the experimentally deduced value for hole-doped WSe 2 51 . If electronic screening is ignored, the effective quasi-2D bare Coulomb potentialV mq is defined by: Our first principles results forV vk F (Supplementary Figure 1) agrees with the analyticallyderived 2D Coulomb potential. At low doping densities,V vk F is very large, which would change the sign of g * compared to the intrinsic g-factor, indicating that screening is important for a meaningful description of g * .

Background dielectric constant
A uniform background dielectric constant ( med ) can be simply added to the dielectric function of the system to obtain the total dielectric function: (r, r , ω) = WSe 2 (r, r , ω) + med −1. In our first principles calculation, the dielectric function is expanded in a plane wave basis: (r, r , ω) = qGG e i(q+G)·r GG (q, ω)e −i(q+G )·r . Thus we approximate the effect of screening by a dielectric medium by modifying the static dielectric matrix as follows:  Here, g * enh = −8.0 and no valley-filling instability occurs. The screened exchange term in our ab initio GW self-energy can be written as: where L is the cell height, S is the cell area, Ω is the cell volume, N q is the number of q points.
In the above expression, the inverse dielectric function is written in a reciprocal space andW mq (E) is defined as ∂Σ(E) ∂f ∂f ∂E can be simplified as: where the delta function picks up the integrand at the Fermi surface and n is the band index of the frontier doped band. We set E = E F , which is equivalent to taking the approximation 0), where E K is the energy at the band extremum in K. This approximation is valid because E K − E F is on the order of meVs, and the dielectric function is fairly constant in this energy range. 3

B. Critical magnetic field B c and condition for valley-filling instability
The condition for the valley-filling instability is given by ∆B > 0. In the following, we derive expressions for B c , B X , and ∆B. Let N K be the K-valley LL so that the energy of 0K lies in between those of N K and (N + 1)K when B < B c . This is equivalent to the condition that N m * < |g * enh | < N + 1 m * B c is defined as the minimum value of B at which 0K is just completely empty with holes. This implies that N K must also be completely occupied with holes at B = B c . All the hole population originally in 0K for B < B c has been transferred to N K at B = B c .
Thus, at B = B c , we have: So, where E F is the position of the Fermi level in the absence of the magnetic field. This condition leads to the expression: This expression for B c gives results very similar to B c ∼ |E F |/(|g * enh |µ B ) given in the main text. B X is obtained by finding the value of B at which N K = 0K , for B > B c : where g * 0 is the renormalized g-factor for B > B c , i.e. the renormalized g-factor of the undoped system, and g * enh refers to the interaction-enhanced g-factor. B X + ∆B is the value of B at which the K-valley LLs from index 0 to (N − 1) can just contain all the holes, and can be obtained as the value of B for the crossing point indicated by the yellow circle in Figure 4d of the main text. Thus, we have: This gives: The instability condition for g * enh can be obtained from ∆B > 0. Together with Eq. 16, we have This expression also implies a condition on N given by Using our prediction of |g * 0 | = 6.6 and m * = −0.48, the minimum N is 4, giving the minimum |g * enh | for a valley-filling instability to be 8.33. * phyqsy@nus.edu.sg 1 Stern, F. Polarizability of a two-dimensional electron gas. Phys. Rev. Lett. 18, 546 (1967).