Exciton g-factors in monolayer and bilayer WSe2 from experiment and theory

The optical properties of monolayer and bilayer transition metal dichalcogenide semiconductors are governed by excitons in different spin and valley configurations, providing versatile aspects for van der Waals heterostructures and devices. Here, we present experimental and theoretical studies of exciton energy splittings in external magnetic field in neutral and charged WSe2 monolayer and bilayer crystals embedded in a field effect device for active doping control. We develop theoretical methods to calculate the exciton g-factors from first principles for all possible spin-valley configurations of excitons in monolayer and bilayer WSe2 including valley-indirect excitons. Our theoretical and experimental findings shed light on some of the characteristic photoluminescence peaks observed for monolayer and bilayer WSe2. In more general terms, the theoretical aspects of our work provide additional means for the characterization of single and few-layer transition metal dichalcogenides, as well as their heterostructures, in the presence of external magnetic fields.

M onolayer (ML) and bilayer (BL) transition metal dichalcogenides (TMDs) such as WSe 2 represent semiconductor building blocks for novel van der Waals heterostructures. By virtue of sizable light-matter coupling governed by excitons 1 , they exhibit versatile potential for applications in photonics and optoelectronics 2,3 , opto-valleytronics 4,5 , and polaritonics 6 . Most recently, the optical interface to TMDs has been instrumental for the observation of strongly correlated electron phenomena in twisted homobilayer and heterobilayer moiré systems [7][8][9] .
The key to further developments of van der Waals heterostructures for fundamental studies and practical devices using TMD MLs and BLs is the detailed understanding of their optical properties. While substantial understanding of zero-momentum excitons in ML and BL WSe 2 has been established 1 , some important aspects remain subject of debate 10 . This holds, in particular, for valley-dark excitons with finite center-of-mass momentum that escape direct optical probes by virtue of momentum mismatch with photons. In MLs, they complement the notion of intravalley spin-bright and spin-dark excitons 1 , and they entirely dominate the photoluminescence (PL) from the lowest-energy states in native homobilayers of WSe 2 (ref. 11 ).
Within the realm of optical spectroscopy techniques, magnetospectroscopy provides means for studying the exciton spin and valley degrees of freedom. Magneto-luminescence experiments on ML WSe 2 in the presence of out-of-plane and in-plane magnetic fields, for instance, have been used to quantify the valley Zeeman splitting of bright excitons [12][13][14][15][16][17] or to brighten spin-dark excitons [18][19][20] , respectively. To date, however, a rigorous assignment of exciton g-factors to intervalley excitons with finite momentum falls short mainly due to the lack of theoretical predictions 10 .
In this work, we develop theoretical methods to evaluate gfactors for excitons in different spin and valley configurations, and provide explicit values for WSe 2 ML and BL excitons composed from electron and hole states away from high symmetry points of the first Brillouin zone. Our calculations go beyond the existing tight-binding models by employing the density functional theory (DFT). We compare our theoretical results with experimentally determined g-factors of intravalley excitons and use them to interpret ambiguous peaks in the PL spectra of ML and BL WSe 2 attributed to intervalley excitons. The technique can be expanded to other materials like WS 2 with similarly complex spectra 21 and large g-factors of ML 17 and BL 22 excitons.

Results
Magneto-luminescence spectroscopy of charge-controlled ML and BL WSe 2 . In our experiments, we used a field-effect heterostructure based on an exfoliated WSe 2 crystal with extended ML and BL regions encapsulated in hexagonal boron nitride (hBN). The device layout is shown schematically in Fig. 1a (see the Methods section for details) and the first Brillouin zone of ML and BL WSe 2 with most relevant points in Fig. 1b. The sample was cooled down to 3.2 K, and the PL was probed as a function of voltage-controlled doping with laser excitation at 1.85 eV and powers below the regimes of neutral and charged biexcitons [23][24][25][26] . Magneto-luminescence experiments were performed in Faraday configuration with a bi-directional solenoid at magnetic fields of up to 9 T (see the Methods section for experimental details).
The evolution of the PL with the gate voltage is shown in Fig. 1c, d for representative spots of ML and BL regions, respectively. In Fig. 1c, the ML reaches the intrinsic limit at gate voltages <−5 V consistent with residual n-doping of the exfoliated crystal 27,28 . The neutral regime is characterized by the bright exciton PL (X 0 ) at 1.72 eV and a series of red-shifted peaks that we label as M 0 1 , M 0 2 , and M 0 3 . None of these peaks with respective red-shifts of 35, 60 and 75 meV from the bright exciton peak is to be attributed to the PL of dark excitons (D 0 ) with 42 meV red-shift 20,29,30 . In our sample, this feature is a rather weak shoulder at the low-energy side of M 0 1 . At positive gate voltages, the ML is charged with electrons and thus exhibits the characteristic signatures of a bright trion doublet (X À 1 and X À 2 ) split by the exchange energy of~6 meV (ref. 28 ), the dark trion (D − ) at 28 meV red-shift from X À 1 (refs. [31][32][33][34], and a hBN SiO 2  series of low-energy peaks dominated by the peak M À 1 at 44 meV red-shift 33,34 . The PL from the BL region in Fig. 1d is characterized by a multi-peak structure, >100 meV below X 0 . It exhibits the same limits of charge neutrality and electron doping as a function of the gate voltage, consistent with the charging behavior of the ML in Fig. 1c. The BL peaks, labeled by an increasing subscript number with decreasing peak energy as B 0 1 through B 0 3 and B À 1 through B À 3 in the neutral and negative regime, respectively, correspond to phonon sidebands of neutral and charged momentum-indirect excitons with a global red-shift of 22 meV at about −7 V (ref. 11 ) in Fig. 1d. According to the single-particle band structure of BL WSe 2 (refs. 35,36 ), the field-induced electron concentration is accommodated at the conduction band edge by the six inequivalent Q-valleys. However, the nature of the hole states that constitute the lowest-energy momentum-dark excitons as long-lived reservoirs of phonon-assisted PL remains ambiguous. The energetic proximity of the valance band edge states at K and Γ in BL WSe 2 (ref. 37 ) renders QK and QΓ excitons and trions (composed from electrons at Q and holes at K or Γ) nearly degenerate, which in turn complicates their energetic ordering 11 .
To examine the origin of the BL peaks and to shed light on the nature of ML peaks with ambiguous or partly controversial interpretation, we performed magneto-spectroscopy in the two well-defined limits of charge neutrality and negative doping. The external magnetic field B was applied along the z-axis perpendicular to the sample. It removes the valley degeneracy and splits the exciton reservoirs by their characteristic Zeeman energies proportional to the exciton g-factor in WSe 2 (refs. [12][13][14][15][16][17] ). The respective polarization-contrasting spectra recorded at −8 T under linearly polarized excitation (π) and circularly polarized detection (σ + and σ − ) for the neutral (negatively charged) ML and BL are shown in the top (bottom) panel of Fig. 2a At each magnetic field, we quantified the experimental Zeeman splitting for every PL peak as the energy difference Δ = E + − E − between the peak energies E + and E − recorded under σ + and σ − polarized detection. The left and right panels of Fig. 3a, b show Δ as a function of the magnetic field for all peaks of the neutral and negatively charged ML and BL, respectively. The set of data derived from magneto-PL measurements was complemented for X 0 , X À 1 , and X À 2 by performing magneto-reflectivity under circular excitation and detection. The corresponding experimental exciton g-factors, obtained from Δ = gμ B B as the slopes of best linear fits to the data in Fig. 3 scaled by the Bohr magneton μ B , are summarized in Table 1. The negative sign of the g-factors reflects the energy ordering of exciton states that exhibit higher (lower) energy for σ − (σ + ) polarized PL peaks at positive magnetic fields.
In ML WSe 2 , the g-factors of both neutral and negatively charged excitons with the corresponding PL peaks X 0 , D 0 , X À 1 , X À 2 , and D − have been established in previous experiments on a wide range of different samples [12][13][14][15][16][17][18][19][20][31][32][33][34] . Our results for the bright exciton and the trion doublet in Table 1 agree well with these reports if we discard the magneto-luminescence result for X À 1 that is compromised by both a vanishingly small PL intensity at high magnetic fields and the relatively broad linewidth of 6 meV in our sample. Due to this inhomogeneous broadening, we are unable to track the dispersion of the relatively weak spin-dark exciton peak D 0 , with g-factors ranging between 9.1 and 9.9 in previous reports 20,31,33,34 nor its chiral-phonon replicum with the same g-factor at 65 meV red-shift from X 0 (refs. 33,34,38 ). The signature of the latter is overwhelmed in our spectra by the peak M 0 2 with 60 meV red-shift and a g-factor of −12.9 ± 0.7 in agreement with values reported from samples with spectrally narrow PL 33,34 . The red-most peak M 0 3 features the same g-factor within the experimental error bars as M 0 1 , suggesting a joint reservoir as their origin. The negatively charged trion D − was reported to have the same g-factor as its neutral counterpart 31-34 , whereas we determine −12.2 ± 0.1. The agreement with previous reports is better for the peak M À 1 with a g-factor of −9.0 ± 0.1 that is supposed to be a phonon sideband of D − (refs. 33,34 ). The latter studies also reported an intense PL peak between M À 1 and D − with a remarkably small g-factor of −4.1 (ref. 33 ) and −3.4 (ref. 34 ). This peak of unidentified origin is missing in our spectra from the negative doping regime.
There are other peaks in ML WSe 2 without conclusive assignment, and in particular M 0 1 has received controversial interpretation as phonon-assisted PL from virtual trions 39 Fig. 2 Magneto-luminescence spectroscopy of charge-controlled monolayer and bilayer WSe 2 . a, b Photoluminescence spectra of monolayer and bilayer WSe 2 , respectively, in a perpendicular magnetic field of −8 T. The neutral and negatively charged regimes are shown in the top and bottom panels, respectively. The spectra were recorded with linearly polarized excitation (π) and circularly polarized detection (σ + and σ − ).
configuration 34 that we denote as K 0 L . Due to the lack of theory for the g-factors of excitons with finite center-of-mass momentum, the task of confronting the competing hypotheses with the characteristic valley Zeeman splittings of controversial ML peaks has remained elusive. The same shortcoming holds for both neutral and charged BL excitons with finite center-of-mass momentum. To shed additional light on the nature of PL peaks in both ML and BL WSe 2 , we calculate in the following the g-factors for excitons in different spin and valley configurations from DFT.
Ab initio calculations of exciton g-factors. We consider a crystal electron in a Bloch state ψ nk ðrÞ ¼ S À1=2 expðikrÞu nk ðrÞ with energy E nk , where n is the band number, k is the wave vector, u nk (r) is the periodic Bloch amplitude, and S is the normalization area. In the presence of a weak perturbation by a static magnetic field B, the first-order correction to the electron energy is proportional to B and given by 40 : where μ B = |e|ℏ/(2m 0 c) is the Bohr magneton, e and m 0 are the charge and mass of the free electron, ℏ is the Planck constant, and c is the speed of light. The expression in square brackets is usually called the effective magnetic moment 41,42 , which contains both spin and orbital contributions. In particular, the first term is proportional to the free electron Landé factor g 0 ≃ 2 and the spin angular momentum s = σ/2, where σ denotes the Pauli matrix. The second term, L n ðkÞ ¼ ψ nk ðrÞ L ψ nk ðrÞ , is the orbital angular momentum with the operator L = ℏ −1 [r × p]. To obtain its matrix elements, one can reduce the calculation to the interband matrix elements of the space coordinate operator r 14,41-43 : where m is the sum over all bands with energy E nk but the band of interest, and ξ nm ðkÞ ¼ i u nk ðrÞ h j∂=∂k u mk ðrÞ j i is the interband matrix element of the crystal coordinate operator.
In the following, we restrict our analysis to the orientation of the magnetic field along the z-axis and define the electron Zeeman splitting as the difference between the energy of the electron state with wave vector +k and spin projection +s along the z-axis and the state with −k and −s as: Thus the electron g-factor of Bloch electrons in the nth band can be written as: with +(−) for s = +1/2 (−1/2) corresponding to spin up (down) projections along z denoted as ↑ (↓), and the explicit expression for the z-component of the orbital angular momentum: where ξ ð ± Þ mn ¼ ðξ ðxÞ mn ± iξ ðyÞ mn Þ= ffiffi ffi 2 p . To calculate the contributions of the conduction (c) band electron with k c , s c and the hole (h) with k h , s h to the exciton gfactor, we neglect electron-hole Coulomb interactions 14,44 . In this case, the exciton Zeeman splitting simplifies to the sum of the Zeeman splittings of the electron and the hole. Using time reversal symmetry that relates the spin and wave vector of the hole to the corresponding spin and wave vector of the empty   Experimental exciton g-factors obtained from magneto-luminescence ( a complementary data from magneto-reflectivity) of neutral and negatively charged monolayer and bilayer WSe2. a X 0 : −4.3 ± 0.1; X À 1 : À4:7 ± 0:3; X À 2 : À6:5 ± 0:4. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18019-1 electron state in the valence (v) band (s h = −s v and k h =−k v ), we obtain the exciton g-factor as Finally, by reference to the valence band electron with k v = K or Γ with spin-up projection s v = +1/2, we discriminate spin-like (L) excitons (with s c = s v ) from spin-unlike (U) excitons (with s c = −s v ). Their respective exciton g-factors are given by: Using Eqs. (7) and (8), we calculate in the following the exciton g-factors from the orbital angular momenta L c (k c ) and L v (k v ) of conduction and valence bands obtained from Eq. (5) within DFT calculations on the Γ-centered k ! grid of 12 × 12 divisions with 300 (600) bands (see the Methods section for details of DFT calculations). In Fig. 4a, b, we show the convergence of the orbital angular momenta L n (k) within our ML and BL calculations as a function of the number of bands taken into account in the sum of Eq. (5). For the ML, Fig. 4a shows the results for the top-most valence band state v at K (blue solid line) and the highest valence band state v at Γ (gray solid line), as well as the two lowest conduction band states c and c + 1 at K and Q (red and black solid and dashed lines). As the BL bands are doubly degenerate, each k-point of the Brillouin zone has at least two bands with L n (k) = L n+1 (k) or L n (k) = L n−1 (k). For the BL in Fig. 4b, we consider the same k-points as for the ML and show the corresponding bands where the orbital angular momenta have the same sign as in the ML case of Fig. 4a.
For the orbital angular momenta of these states, convergence is observed above 275 and 550 bands in the case of ML and BL in Fig. 4a, b, respectively, with the factor of two difference related to the doubled number of atoms in BL calculations. We note that the values for the valence band states at Γ must vanish by symmetry arguments, whereas our numerical calculations yield ±0.01 for both ML and BL. This marginal discrepancy is due to a finite number of bands taken into account and can be used to estimate the precision of our numerical calculations. The corresponding bound on the absolute error of the exciton g-factors from DFT, given explicitly in Table 2 for selected exciton configurations, is thus in the order of ±0.05.
As evident from Fig. 4a, particular bands make decisive contributions to the g-factor. To discuss this behavior for the ML case in more detail, we consider the 24th and 26th bands that correspond to the highest valence band (v) and the second conduction band (c + 1), respectively, and give rise to largest mutual contributions in the g-factors. This is expected according to Eq. (5), where the orbital angular momentum is proportional to the product of the interband matrix elements, which in turn are largest for the fundamental A-exciton transition X 0 between the 24th and 26th bands. Similar arguments apply for the mutual contributions of the 23rd and 25th bands to the g-factor of Bexcitons. It is also instructive to note the different dependencies of the orbital momenta for the two lowest conduction bands (L c and L c+1 ) and the top valence band (L v ) on the number of bands included. In Fig. 4a, L c and L c+1 exhibit jumps at 23rd and 24th bands, respectively, and then increase only marginally. In contrast, L v in Fig. 4a increases nearly monotonously beginning from m~30 on. A closer inspection shows that for m > 26, the sign of the square bracket in Eq. (5) alternates with increasing m, and the terms of comparable absolute values therefore cancel each other for both L c and L c+1 . For L v , on the other hand, the absolute values of the positive terms systematically exceed the negative terms and thus L v continues to grow with increasing band number. Further analysis will be required to understand this behavior in more detail.
The DFT results for L n (k) within the first Brillouin zone are shown in Fig. 5. Since spin-orbit effects were included at the DFT level, it is instructive to show both spin-orbit split highest valence bands (v and v − 1) and lowest conduction bands (c and c + 1). With the matrix elements of the orbital angular momenta of the valence and conduction bands in Fig. 5, it is straight forward to  Fig. 4 Electron orbital angular momentum in monolayer and bilayer WSe 2 from DFT. a, b Electron orbital angular momentum from DFT calculations for the highest valence bands at K and Γ and the lowest conduction bands at K and Q in monolayer and bilayer WSe 2 , respectively.
Exciton g-factors for selected spin-valley configurations of excitons in monolayer and bilayer WSe2. Note that without further assumptions the sign of the g-factor is meaningful only for zero-momentum spin-like excitons with valley-contrasting dipolar selection rules.
calculate the g-factors of the lowest-energy ML excitons in various configurations. In Table 2, we list the g-factors obtained from our DFT results for excitons in different configurations of valleys (k c , k v ) and spins (s c , s v , with ↑ or ↓ projection along z).
In the top block of Table 2, we list excitons with the hole at K and the electron at K or K 0 in spin-like and spin-unlike configurations with short exciton notation for zero-momentum bright and dark neutral excitons X 0 and D 0 and their finitemomentum counterparts K 0 L and K 0 U . The block below shows the results for the spin-like and spin-unlike Q-excitons with the electron in Q and the hole in K, followed by two blocks without short exciton notation for momentum-indirect excitons composed from electrons in K or Q and holes in Γ. Note that the sign of the g-factor can be determined without further assumptions only for X 0 with established valley-contrasting dipolar selection rules. For Zeeman-split momentum-indirect excitons, on the other hand, the sign will depend on the symmetry of the actual phonons involved in phonon-assisted PL 34 . This is analogous to the case of spin-dark excitons D 0 with linearly polarized in-plane zero-phonon emission 30 contrasted by circularly polarized PL sidebands of the same reservoir mediated by chiral phonons 33,45 . In principle, not only the g-factor signs but also the absolute values of the respective PL peaks should be distinct due to different exciton-phonon coupling not accounted for in our model. However, we expect such higher-order corrections to be well below the resolution of our measurements.

Discussion
First, we discuss the results of our calculations for excitons in ML WSe 2 . The g-factor of −4.0 from our DFT model is in excellent agreement with the experimental value of −4.1 for X 0 (refs. [12][13][14][15][16][17]. The good agreement in the g-factor of spin-dark excitons with g ≃ 9.4 in experiment 20 and 10.1 in DFT provides further confidence in our model. According to our calculations, the states K 0 L and K 0 U , which are the momentum-indirect counterparts of X 0 and D 0 , exhibit different g-factors with large values of 13.6 and 19.6, respectively. The g-factors of Q-momentum excitons (9.2-14.5) are similar to those of D 0 and K 0 L , whereas excitons with the hole at Γ are predicted to have rather small g-factors (<5.8) except for the spin-unlike configuration with K 0 electron (9.8). As expected, the g-factors of intralayer excitons in BL WSe 2 are close to the values of the corresponding ML excitons 46 . In addition to intralayer excitons, the BL hosts interlayer counterparts (e.g., intralayer Q L and interlayer Q 0 L , intralayer Q U and interlayer Q 0 U , so on) that exhibit the same g-factors within our model, which neglects Coulomb corrections for intralayer and interlayer excitons.
By providing explicit g-factor values for momentum-indirect excitons, our DFT results complement the experimental observations in ML and BL WSe 2 . In the framework of neutral MLs, however, they do not resolve the ambiguity between the two competing explanations of the peak M 0 1 . The assignment of the peak as a phonon sideband of Q-momentum excitons 21 , on the one hand, is consistent with the g-factors of 9.2 and 14.5 for Q L and Q 0 U states in Table 2 (note that Q U and Q 0 L excitons, 250 meV above degenerate Q L and Q 0 U states 47 , are irrelevant in this context) and the structured peak M 0 1 in Fig. 2 with a g-factor of 11.5. On the other hand, the interpretation of the peak as direct PL emission by momentum-dark K 0 L excitons 34 is also consistent with the theoretical g-factor of 13.6 from DFT. Our DFT results also identify KΓ and QΓ with small g-factors as potential candidates to explain the bright PL peak between M À 1 and D − in the negatively charged regime of high-quality samples with narrow spectra 33,34 .
For the neutral BL, our results help to rule out QΓ excitons and suggest spin-unlike interlayer QK and intralayer Q 0 K exciton reservoirs rather than K 0 Γ as a joint origin of phonon sidebands B 0 1 , B 0 2 , and B 0 3 (ref. 11 ). Whereas a detailed assignment of the neutral BL peaks to the specific reservoirs and phonon sidebands is yet to be developed, the values of the exciton g-factors in the charged regime can be understood, as in the ML case, by regarding the additional electron in the charged complex simply as a spectator to the Zeeman effect of the neutral finitemomentum exciton reservoir.
In summary, our work provides exciton g-factors for neutral and charged ML and BL WSe 2 from both experiment and DFT. For ML WSe 2 , the g-factors obtained from first-principles calculations are in excellent quantitative agreement with previous reports and complement these studies by providing theoretical gfactors for momentum-indirect excitons in different configurations of spins and valleys. For BL WSe 2 , our work adds insight into the origin of PL peaks on the basis of theoretical g-factor values. In the broad context of research on layered semiconductors and their applications, the theoretical aspects of our work provide guidelines for magneto-optical studies of singlelayer TMDs, homobilayer or heterobilayer systems, and other realizations of TMD-based van der Waals heterostructures.
Note: During the submission of our manuscript, we became aware of three related works on the theory of exciton g-factors in TMD MLs and heterostructures from first principles [48][49][50] .

Methods
Experimental methods. The field-effect heterostructure consisted of an exfoliated WSe 2 crystal (HQ Graphene) with extended ML and BL regions encapsulated in hBN (NIMS). To control the charge doping, the crystal was contacted by a gold electrode deposited on a 50-nm-thick thermal silicon oxide layer of a p-doped silicon substrate. With the electrode grounded, a gate voltage applied to the highly doped silicon was used to control the doping level in ML and BL WSe 2 . The sample was mounted in a cryogenic confocal microscope and cooled down in a closedcycle magneto-cryostat (attocube systems, attoDRY1000) with a base temperature of 3.2 K. The PL was excited at 1.85 eV with a few μW power of a continuous-wave laser diode focused to the diffraction-limited confocal excitation and detection spot of a low-temperature apochromatic objective (attocube systems, LT-APO/VISIR/ 0.82), dispersed with a monochromator (Roper Scientific, Acton SP2500), and detected with a nitrogen-cooled CCD (Roper Scientific, Spec 10:100BR/LN). Magneto-luminescence experiments were performed in Faraday configuration with a bi-directional solenoid at magnetic fields of up to 9 T.
DFT calculations. DFT calculations were performed within the generalized gradient approximation with the PBEsol exchange-correlation functional 51 as implemented in the Vienna ab initio simulation package. The van der Waals interactions were considered with the DFT-D3 method with Becke-Johnson damping 52,53 ; the spin-orbit interaction was included at all stages. Elementary cells with a vacuum thickness of 30 Å were used in order to minimize interactions between periodic images. The atomic positions were relaxed with a cut-off energy of 400 eV until the change in the total energy was <10 −6 eV. The band structure of ML (BL) was calculated on the Γ-centered k ! grid of 12 × 12 divisions with 300 (600) bands.
Data availability