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

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 siliconOhkawa ; Sham , but the ability to exploit valley polarizations in these materials has been limited NatRevMat . 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 NatRevMat ; RMP2018 . The spin-valley locking effect in these MLs Xiao2012 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 Xuan2020 ; Rohlfing2020 ; PRB2020 ; PRB2020 ; NC2020 ; Marie2020 ; 2018Ensslin ; LL2019 ; Mak ; 2020PRL ; MakPRL ; jump ; NanoLett ; NatMater ; LL2017 ; 2018PRBR is also significantly larger than in traditional semiconductors Roth ; Yafet ; Janak ; Fang ; Lak ; ExpPRB ; Shayegan .

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 Xuan2020 ; Rohlfing2020 ; PRB2020 ; Marie2020 than that in silicon, where only the spin magnetic moment dominates Roth ; Yafet . 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 NatRevMat . 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 NatRevMat ; RMP2018 . It has been found that carrier doping can dramatically enhance the g-factors in TMD monolayers 2020PRL ; MakPRL ; jump ; NanoLett ; NatMater ; LL2017 ; 2018PRBR , 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 2020PRL ; MakPRL ; jump ; NanoLett ; NatMater ; LL2017 ; 2018PRBR ; Janak ; ExpPRB , 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-1970sJanak ; Ando . 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^{*}_{\text{enh}}$ for which the discontinuous change in g-factor at $B_{c}$ results in a LL alignment and valley-filling instability for $B\gtrsim B_{c}$. Such a phenomenon has not been observed or predicted for silicon and other conventional valleytronic materials. Our computed interaction-enhanced g-factors for hole-doped ML WSe₂ agree well with experiment 2020PRL ; jump ; NatMater and can be tuned by dielectric screening. The predicted valley-filling instability for $B\gtrsim B_{c}$ provides theoretical insights into recent experimental observations jump of a pronounced Landau level-filling instability at a critical magnetic field which closely matches our predicted values. The associated alignment of LLs is of interestjump to investigate quantum phase transitions in these doped TMDs Braz ; Donk ; Miserev ; Roch ; Nature1999 ; Science2000 . The recent observation of fractional quantum Hall states associated with non-abelian anyons in ML WSe₂ FQHWSe2 highlights the potential of creating pseudo-spinors from aligned LLs for topological quantum computing applications QCRMP .

RESULTS

Interaction-enhanced g-factors

For a many-electron system described by a static mean-field Hamiltonian $H\ket{n\mathbf{k}}=E_{n\mathbf{k}}\ket{n\mathbf{k}}$, it has been shown that an out-of-plane magnetic field $B$ results in the following expression for the LLs at the $K$ valley Xuan2020 :

where $n$ is the corresponding band index, $E_{n\mathbf{K}}$ is the mean-field single-particle energy at $K$, $m^{*}$ is the valley effective mass, $\mu_{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, $g^{I}_{n\mathbf{K}}=g^{\text{orb}}_{n\mathbf{K}}-g_{s}s_{z,n\mathbf{K}}$, where $g_{s}$ is taken to be $2.0$ and $s_{z,n\mathbf{K}}$ is the spin quantum number. The orbital g-factor is defined using the orbital magnetic moment $g^{\text{orb}}_{n\mathbf{K}}\mu_{B}=m^{z}_{n\mathbf{K}}$ Xuan2020 (see Methods).

However, the above static mean-field description does not account for the energy-dependent electron self-energy $\Sigma(E)$ that is necessary for a many-body description of the quasiparticle (QP) energies. $\Sigma(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^{*}_{n\mathbf{K}}$, defined as $(E_{n\mathbf{k}}^{\text{QP},1}-E_{n\mathbf{k}}^{\text{QP}})=g^{*}_{n\mathbf{K}}\mu_{B}B$, where

Thus we have,

Such a renormalization effect is missing in previous first principles calculations of g-factors in TMDs Xuan2020 ; Rohlfing2020 ; PRB2020 .

The electron self-energy in this work is computed within the $GW$ approximation Hedin (see Methods), which uses the first order term in the perturbative expansion of $\Sigma$ in terms of the screened Couloumb interaction $W$. Table 1 shows the renormalized and intrinsic g-factors computed for undoped monolayer WSe₂. We see that the magnitudes of the renormalized g-factors are reduced by $\sim 20$% compared to the intrinsic $GW$ g-factors, because $\frac{\partial\Sigma(E)}{\partial E}$ is in general negative Louie . The exciton g-factors deduced using the renormalized g-factors are in good agreement with experiment 2d2015 ; NC2019 ; NatPhys2015 ; 2d2019 ; Liu2019 .

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 $\frac{d\Sigma(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. $\bar{W}_{nk_{F}}$ (defined in Methods) is an effective quasi-2D screened Coulomb potential which can be evaluated completely from first principles. Since $\bar{W}_{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\Sigma(E)/dE$ derived in Ref. Janak, (see also Ref. Quinn, ) where the first term in Eq. (4) is ignored. We note that our numerical results for $\bar{W}_{vk_{F}}$ as a function of hole density differ from those computed for an ideal 2D fermion gasstern ; Janak , 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 $v\mathbf{K}$ is omitted. The magnitudes of our computed valence band g-factors $|g^{*}|$ are plotted as red squares in Fig. 1b for ML WSe₂ with different hole densities. Our predicted g-factors agree well with those deduced from multiple experiments jump ; 2020PRL ; Marie2020 ; NatMater on hole-doped ML WSe₂. The renormalized g-factor for the undoped system is labeled $g^{*}_{0}$. Due to interaction-induced 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^{\prime}$ 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^{*}_{\text{enh}}$. However, above a critical magnetic field $B_{c}\sim|E_{F}|/(|g^{*}_{\text{enh}}|\mu_{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\times 10^{12}\text{cm}^{-2}$, $B_{c}\approx 36$ T, and we have $g^{*}\approx-11.0$ for $B<B_{c}$ and $g^{*}\approx-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^{*}_{\text{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\times 10^{12}\text{cm}^{-2}$, we predict $B_{c}$ in silicon to be $\sim 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. Xuan2020

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^{\prime}$jump . In Fig. 3a, we plot these values of $B_{c}$ (blue circles) and compare them with our predicted values (red squares). The predicted dependence of $B_{c}$ on the hole density agrees well with experiment.

How can one maximize the concentration of valley-(and spin-)polarized carriers in the TMD ML? As the hole concentration increases, $B_{c}$ increases, giving a larger range of $B$ for which $g^{*}$ is enhanced by interactions (Fig. 3a). However, the magnitude of $g^{*}_{\text{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 $\rho_{0}$ which maximizes the Zeeman split $E_{Z}$ (Fig. 3b). This optimal hole concentration $\rho_{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^{\prime}$ and $NK$ (purple circle, $N=5$ in Fig. 4a). If $NK$ 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 $NK$ is fully occupied, the plot is shaded from $\epsilon_{NK}+\frac{\hbar\omega_{c}}{2}$ down to $\epsilon_{NK}-\frac{\hbar\omega_{c}}{2}$. At $B=B_{c}$, an integer number of LLs ($5$ in Fig. 4c) are completely filled with holes while $0K^{\prime}$ is completely empty. As $B$ decreases slightly below $B_{c}$, the LL degeneracy decreases, and the LL with the next lower energy ($0K^{\prime}$ 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 $\epsilon_{5K}>\epsilon_{0K^{\prime}}$, 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^{\prime}$ and $5K$ cross, holes will start to fill the $K^{\prime}$ 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^{\prime}$ is depleted of holes from $B=B_{c}$ to $B=B_{X}$, and filled again up to $B_{X}+\Delta B$. In practice, when holes begin to fill $0K^{\prime}$, 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}+\Delta B$.

Our predictions provide important theoretical insights into a recent experiment on doped ML WSe₂jump , where optical absorption plots showed a pronounced signature of the peak positions changing from one inter-LL transition to another over a small range of $B$ close to the onset of the fully-polarized regime in the experiment. This is consistent with the highest occupied LL in the $K^{\prime}$ valley being emptied and partially filled with holes at $B\sim B_{c}$ in our predictions. The authors of Ref. jump, attributed this observation to the oscillatory g-factors predicted for traditional semiconductors such as silicon Ando . However, in this theory, the changes in the g-factors are directly related to the position of the Fermi level relative to the LLs, and the g-factors therefore have an “oscillatory” Ando dependence on $B$ rather than a pronounced change at one particular value of $B$ as seen in the experiment. Furthermore, such a pronounced instability was not observed in experiments Fang ; Lak ; ExpPRB ; Shayegan on the g-factors in doped silicon and other traditional valleytronic materials for which these oscillatory g-factors were predicted. Thus, this pronounced instability observed in doped ML WSe₂jump is in fact a manifestation of the valley-filling instabilities that are predicted here to emerge specifically for doped TMDs. Our conclusion is further supported by the fact that the measured values of $B_{c}$ and $B_{X}$ are respectively $32$ and $38$ T for $|g^{*}_{\text{enh}}|\sim 11$jump , close to our predicted values of $31$ and $37$ T for the same $|g^{*}_{\text{enh}}|$ (see also Supplementary Table 1).

The alignment of LLs is of interest to investigate quantum phase transitions in these doped TMDs jump ; Braz ; Donk ; Miserev ; Roch ; Nature1999 ; Science2000 . Given that the LLs are expected to align for $B$ between $B_{X}$ and $B_{X}+\Delta B$, it is interesting to predict how large $\Delta B$ can be and how sensitive $\Delta B$ is to fluctuations in $g^{*}_{\text{enh}}$. Not all values of $g^{*}_{\text{enh}}$ will result in the instability (see Supplementary Figure 2 and Supplementary Note). In particular, if the energies of $0K^{\prime}$ and $NK$ cross at $B_{X}$, the valley-filling instability only occurs if $NK$ is occupied with holes. Fig. 4e plots the ranges of $g^{*}_{\text{enh}}$ for which a valley-filling instability will occur, as well as the corresponding $\Delta B$ values. We see that an optimal range of $|g^{*}_{\text{enh}}|$ for LL alignment is $10.4<|g^{*}_{\text{enh}}|<11.2$. Here, $\Delta B$ is quite large ($\sim 8$T for $|g^{*}_{\text{enh}}|\sim 10.4$) and is also fairly robust to changes in $g^{*}_{\text{enh}}$. The corresponding values of $B_{c}$ and $B_{X}$ fall within $30$ to $40$ T (Supplementary Table 1), well within the reach of experiments.

The alignment of LLs in different valleys can in principle be achieved for $B<B_{c}$ if $g^{*}_{\text{enh}}$ can be tuned such that $0K^{\prime}$ matches exactly with $NK$ for some $N$. However, once $g^{*}_{\text{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^{*}_{\text{enh}}$.

Tunability of interaction-enhanced g-factor

We further note that in addition to electrostatic gating which changes the carrier concentration and thus $g^{*}_{\text{enh}}$ (Fig. 1), $g^{*}_{\text{enh}}$ can also be tuned by dielectric screening (Fig. 5). The tunability of $g^{*}_{\text{enh}}$ with the background dielectric constant can be understood from the fact that $g^{*}_{\text{enh}}$ is related to the effective quasi-2D screened Coulomb potential at the Fermi surface (Eq. 5). This tunability of $g^{*}_{\text{enh}}$ provides a handle to control the valley-polarized current, $B_{c}$ and $\Delta 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 $B_{c}$ above which the g-factors revert to those in the undoped systems. Such a phenomenon has not been predicted or observed in silicon and other traditional valleytronic materials, because the corresponding $B_{c}$ would be much larger due to the smaller g-factors in these materials.

The enhancement in g-factors arises from the effect of a magnetic-field-induced change in occupancies on the screened exchange interactions. This effect is only present in the mixed-polarized regime ($B<B_{c}$). As the carrier concentration increases, $g^{*}$ decreases and $B_{c}$ increases, so that for any value of the magnetic field $B_{0}$, the valley-polarization is maximized when the carrier concentration is such that $B_{c}=B_{0}$. This prediction has implications for maximizing the valley- and spin-polarized current for the valley Hall effect.

The computed interaction-enhanced g-factors agree well with experiment for different doping concentrations. We further identify the values of $g^{*}_{\text{enh}}$ and corresponding ranges of $B$ that lead to a valley-filling instability and expected LL alignment, which are of interest for the investigation of quantum phase transitions in doped TMDs Braz ; Donk ; Miserev ; Roch ; Nature1999 ; Science2000 . Recent observations of fractional quantum Hall states associated with non-abelian anyons in ML WSe₂ FQHWSe2 suggest that the creation of pseudo-spinors from a linear combination of valley-aligned LLs can be useful for topological quantum computing applications QCRMP .

METHODS

Calculation of intrinsic g-factor

The orbital component of the intrinsic g-factor $g^{\text{orb}}_{n\mathbf{K}}$ is defined as $g^{\text{orb}}_{n\mathbf{K}}\mu_{B}=m^{z}_{n\mathbf{K}}$ Xuan2020 :

We use the PBE exchange-correlation functional PBE1996 for the DFT mean-field calculations QE and the details follow those in Ref. Xuan2020, . For $GW$ calculations of the intrinsic g-factors, we use a non-uniform sampling method sub of the Brillouin Zone starting with a $12\times 12$ k-grid as implemented in the BerkeleyGW code BGW . The energy-dependence of the dielectric function is treated within the generalized plasmon pole (GPP) model Louie . An energy cutoff of $35$ Ry with $4000$ empty bands is used for the reciprocal space expansion of the dielectric matrix. The intrinsic single-band g-factor reduces by only $0.2$ when an energy cutoff of $2$Ry 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 $\frac{d\Sigma(E)}{dE}$ using Eq. 3. As the dependence of $E_{n\mathbf{k}}^{\text{QP}}$ on $B$ is no longer linear, we generalize Eq. (2) to the case where $E_{n\mathbf{k}}^{\text{QP}}$ 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₂ is doped with holes.

The first term of Eq. 4 can be obtained directly from the BerkeleyGW outputBGW .

$\Sigma=iGW$ can be partitioned Louie into the dynamical non-local screened-exchange (SEX) and Coulomb-hole (COH) interaction terms $\Sigma=\Sigma^{\text{SEX}}+\Sigma^{\text{COH}}$. Only the screened exchange term depends on the occupancies $f$ and contributes to the second term of Eq. 4. The screened exchange energy $\Sigma^{\text{SEX}}$ in our ab initio plane-wave calculation can be written as (see Supplemental Note):

where we define the quasi-2D screened Coulomb potential:

Here, $\Omega$ is the cell volume, $N_{\mathbf{q}}$ is the number of q-points and $v_{\mathbf{q}}$ is the Coulomb potential with the slab Coulomb truncation scheme applied trunc . $\bar{W}_{m\mathbf{q}}$ is an effective quasi-2D screened Coulomb potential defined in valley $K$ and $L$ is the height of the supercell for ML WSe₂. The second term in Eq. (4) is then given by (see Supplemental Note):

We compute $g^{*}_{v\mathbf{K}}$ by evaluating $\bar{W}_{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\times 120$, a $2$Ry $G$-vector cut off and $29$ bands. $g^{*}$ is unchanged when we use instead $4$Ry and $299$ bands, and reduces by $\sim 3\%$ when a $240\times 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₂ LL2016 . If electronic screening is ignored, the effective quasi-2D bare Coulomb potential $\bar{V}_{mq}$ is defined by:

Our first principles results for $\bar{V}_{vk_{F}}$ (Supplementary Figure 1) agrees with the analytically-derived 2D Coulomb potential. At low doping densities, $\bar{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 ($\epsilon^{\text{med}}$) can be simply added to the dielectric function of the system to obtain the total dielectric function: $\epsilon(\mathbf{r},\mathbf{r}^{\prime},\omega)=\epsilon^{\text{WSe}_{2}}(\mathbf{r},\mathbf{r}^{\prime},\omega)+\epsilon^{\text{med}}-1$. In our first principles calculation, the dielectric function is expanded in a plane wave basis: $\epsilon(\mathbf{r},\mathbf{r}^{\prime},\omega)=\sum_{\mathbf{qGG^{\prime}}}e^{i(\mathbf{q}+\mathbf{G})\cdot\mathbf{r}}\epsilon_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega)e^{-i(\mathbf{q}+\mathbf{G}^{\prime})\cdot\mathbf{r}^{\prime}}$. Thus we approximate the effect of screening by a dielectric medium by modifying the static dielectric matrix as follows:

ACKNOWLEDGEMENTS

This work is supported by the NUS Provost’s Office, the Ministry of Education (MOE 2017-T2-2-139) and the National Research Foundation (NRF), Singapore, under the NRF medium-sized centre programme. Calculations were performed on the computational cluster in the Centre for Advanced 2D Materials and the National Supercomputing Centre, Singapore.

Supplementary Information for: Valley-Filling Instability and Critical Magnetic Field for Interaction-Enhanced Zeeman Response in Doped WSe₂ Monolayers