Stacking-dependent exciton multiplicity in WSe$_2$ bilayers

Twisted layers of atomically thin two-dimensional materials realize a broad range of novel quantum materials with engineered optical and transport phenomena arising from spin and valley degrees of freedom and strong electron correlations in hybridized interlayer bands. Here, we report experimental and theoretical studies of WSe$_2$ homobilayers obtained in two stable configurations of 2H ($60^\circ$ twist) and 3R ($0^\circ$ twist) stackings by controlled chemical vapor synthesis of high-quality large-area crystals. Using optical absorption and photoluminescence spectroscopy at cryogenic temperatures, we uncover marked differences in the optical characteristics of 2H and 3R bilayer WSe$_2$ which we explain on the basis of beyond-DFT theoretical calculations. Our results highlight the role of layer stacking for the spectral multiplicity of momentum-direct intralayer exciton transitions in absorption, and relate the multiplicity of phonon sidebands in the photoluminescence to momentum-indirect excitons with different spin valley and layer character. Our comprehensive study generalizes to other layered homobilayer and heterobilayer semiconductor systems and highlights the role of crystal symmetry and stacking for interlayer hybrid states.


I. INTRODUCTION
The optical properties of transition metal dichalcogenide (TMD) semiconductors are governed by excitons in different spin, valley and layer configurations [1,2]. Among possible realizations of TMD systems, heterobilayers and homobilayers stand out as hosts of excitons with layer-indirect character. Initial studies of homobilayers were limited to natural 2H layer stacking with antiparallel alignment or 60 • twist angle [3][4][5], extracted by exfoliation from native crystals and structurally different from 3R stacking with parallel alignment or 0 • twist. Recently, variations of the twist angle in WSe 2 bilayers (BLs) away from 2H and 3R stackings revealed novel phenomena ranging from correlated electronic phases [6] to moiré exciton physics [7,8] with angle-controlled exciton valley coherence and dynamics [8][9][10], Coulomb correlations in effectively flat moiré exciton bands [11] or optical nonlinearities [12].
A particularly attractive feature of interlayer excitons in related homobilayer systems is provided by the permanent dipole moment of layer-indirect excitons [13][14][15] which promotes exciton-exciton dipolar interactions, allows employing electric fields to tune the optical properties via the Stark effect [10,[16][17][18][19][20][21] or implementing exciton trapping and routing [22]. The electrostatic dipole moment depends on the degree of exciton layer delocalization, which in turn is sensitive to the interlayer coupling and thus to the stacking order, as was shown recently by optical absorption for 2H and 3R MoS 2 BLs for transitions involving momentum direct KK excitons [23,24]. In contrast to absorption probing zero-momentum exciton transitions, the photoluminescence (PL) of BLs is dominated by momentum-indirect excitons [25,26] via luminescence phonon sidebands [27][28][29] which exhibit shift linear energy shifts in perpendicular electric fields according to the Stark effect [10,16,20,21]. The most recent observation of two different slopes in the energy dispersion of phonon sideband emission from 2H WSe 2 BL in electric field [20,21] indicates the presence of two degrees of electron-hole layer separation, with respective dipole moments attributed to excitons in distinct QK and QΓ reservoirs [20,21]. Experiments employing strain tuning [30] and magnetic fields [28], on the other hand, suggest that the PL sidebands stem exclusively from QK excitons, indicating shortcomings in the present understanding of the underlying lowest-energy exciton reservoirs with finite center-of-mass momenta.
To provide comprehensive insight into the nature of exciton reservoirs in BL WSe 2 with different spin, momentum and layer character, we performed experimental and theoretical studies of 2H and 3R stackings with contrasting spectroscopic responses. To this end we synthesized WSe 2 crystals by chemical vapor deposition (CVD), yielding 2H and 3R as two stable limits of relative layer orientation. As opposed to BL crystals of MoS 2 [23,24,31] and WS 2 [32,33] synthesized by CVD, studies of WSe 2 BLs in 3R stacking have remained elusive due to the challenge of perfect layer alignment. The only realization of nominal zero-angle twist so far has been obtained from exfoliation stacking [10], which can only approximate the ideal 3R layer order inherent to CVD growth. Moreover, the spectral features of BLs aligned near zero twist can be compromised by marginal-angle . c and d, Room-temperature Raman spectra of ML and BL WSe2 in 2H and 3R stacking (the spectra were recorded with excitation at 532 nm and normalized to the peak intensity of the Si phonon mode at 520 cm −1 ; note that the intensity of higher-order Raman modes associated with multi-phonon processes were scaled for increased visibility). The in-plane and out-of-plane phonon modes, E and A 1 , are degenerate in the ML but split by 2 cm −2 both in 2H and 3R homobilayer WSe2. e, Normalized photoluminescence (black) and differential reflectivity (red) spectra at 3.2 K, recorded in the ML region of hBN-encapsulated BL WSe2 crystals. The ML exciton (X) and biexciton (XX) peaks exhibit characteristic linear and superlinear scaling of the integrated PL intensity as a function of the excitation power, as evident from data and respective power-law fits IP L ∝ P α shown by solid lines in f.
In the following, we present a comprehensive study of excitons in 2H and 3R WSe 2 BLs, performed with cryogenic optical spectroscopy on CVD-synthesized highquality crystals. With a set of complementary spectroscopy techniques, we identify contrasting responses of 2H and 3R BL configurations and observe rich exciton multiplicity in the 3R case. We relate the differences between the optical spectra of 2H and 3R BL stackings to the nature of excitons in different spin, valley and layer configurations. Our interpretation is based on firstprinciples calculations of the band structure, optical absorption and exciton g-factors. Theory analysis shows that in contrast to 2H stacking, the top and bottom layers in 3R stacking differ in their band structure and exhibit different optical bandgaps. This insight resolves the recent puzzle of two distinct intralayer KK exciton transitions observed in absorption on nearly aligned WSe 2 BLs [10]. Moreover, our theoretical analysis highlights the multiplicity and degree of layer delocalization FOR interlayer excitons in both 2H in 3R stackings, thus providing an intuitive explanation for different electrostatic dipole moments of QK excitons without the requirement of involving QΓ states [20,21].

II. RESULTS AND DISCUSSION
WSe 2 BL crystals in 2H and 3R stacking were obtained from CVD synthesis. Atomic force micrographs in Fig. 1a and b show respective WSe 2 crystals with ML and BL regions in 2H and 3R configurations with 60°a nd 0°rotational alignment of the two WSe 2 monolayers (MLs), as can be deduced directly from the edges of CVD flakes with the same edge termination. The correspond-ing atomic registries are illustrated in the bottom panels of Fig. 1a and b, showing side (left) and top (right) views and the BL unit cell by dashed boxes.
Initial characterization of the two distinct crystal configurations was performed with Raman spectroscopy at ambient conditions. A typical Raman spectrum from a ML region is shown in Fig. 1c. It features degenerate inplane E and out-of-plane A 1 first-order Raman modes around 249 cm −1 , the double resonance 2LA mode at 257 cm −1 , and a series of multi-phonon modes at 358, 372, and 395 cm −1 , all of which are consistent with ML features [43][44][45][46][47][48]. The degeneracy of the E and A 1 modes, characteristic of MLs [47,48], is lifted in the Raman spectra of both 2H and 3R BLs shown in Fig. 1d. We determined by peak decomposition that in both BL stackings the respective E 1 2g and A 1g modes exhibit a ∼ 2 cm −1 frequency splitting due to interactions between the layers. We also observe a red-shift of the E /E 1 2g and a blue shift of the A 1 /A 1g modes for both stackings [49]. The peak at 308 cm −1 , labelled as A 2 , is absent in the ML limit but very pronounced in the Raman spectra of both 2H and 3R stackings as a hallmark of BL regions [43,[46][47][48]. The spectra of multi-phonon modes at frequencies above 350 cm −1 are similar in ML and BLs.
Having identified 2H and 3R BL crystals by AFM and Raman spectroscopy, we performed cryogenic micro-PL and differential reflection (DR) spectroscopy. To this end, the CVD-grown crystals were encapsulated in hexagonal boron nitride (hBN) by standard lift-off and exfoliation and transferred onto a Si/SiO 2 substrate. Typical PL and DR spectra of ML regions, recorded at 3.2 K, are shown in Fig. 1e. Both spectra feature the resonance of the fundamental exciton (X) around 1.71 eV with a spectrally narrow full-width at half-maximum (FWHM) linewidth of 4 meV, exceeding the transformlimited linewdith of 1 meV in best hBN-encapsulated MLs obtained from native crystals [50][51][52][53] yet slightly narrower than 5.2 meV reported for CVD-grown MoS 2 ML in hBN [53]. The second pronounced peak in the PL spectrum at 1.69 eV corresponds to the biexciton (XX) emission with a characteristic 20 meV red-shift [54][55][56][57][58] and superlinear dependence on the excitation power shown in Fig. 1f. All main features in PL and DR as well as the absence of trion signatures indicate vanishing residual doping and the overall high quality of our CVD-grown crystals.
The cryogenic PL and DR of BL regions with 2H and 3R stacking are shown in Fig. 2a and b with pronounced differences in the spectral characteristics. The 2H spectra in Fig. 2a agree with the results of previous studies on exfoliated BLs from bulk WSe 2 [10, 28-30, 59, 60], with the neutral exciton resonance X at around 1.72 eV associated with momentum direct transitions between conduction and valence band states at K, as well as a series of PL peaks between 1.55 and 1.61 eV we label as P 1 to P 4 with increasing red-shift. Consistent with the indirect band gap of WSe 2 BL between the conduction band minima at the six inequivalent Q-points of the first Bril- louin zone and the valence band maximum at K [48,61], the lowest energy exciton states form as QK and Q K with finite center-of-mass momentum and characteristic phonon sideband luminescence into the peaks P 1 through P 4 [28,29]. As on ML regions, the high-quality of our CVD-synthesized crystals is reflected by spectrally narrow features of BLs in both stackings. The characteristic features of 2H WSe 2 BL are contrasted by the spectra from 3R stacking. First, the peak multiplicity increases both in PL and DR spectra of Fig. 2b. We observe two distinct exciton resonances X 1 and X 2 around 1.70 and 1.68 eV with an energy splitting of 17 meV and spot-to-spot variations between 16 and 20 meV, which is very different from the single transition observed throughout the 2H stacked sample. As a second important observation, we find in the region of phonon sidebands between 1.50 and 1.55 eV three additional peaks emerging in between P 1 through P 4 , which we label as P 1 , P 2 and P 3 . Moreover, we observe a redshift of 16 meV from the exciton resonance X in 2H to X 1 in 3R, as well as a larger red-shift to the first peak P 1 , increasing from ∼ 100 meV in 2H to ∼ 150 meV in 3R as a consequence of stronger interlayer hybridization in the latter case.
To understand the increased multiplicity of exciton transitions as well as their energy shifts and splittings in 3R stacking, we carried out band structure calculations including excitonic effects via the GW+BSE approach. Our calculations, performed for WSe 2 BLs in vacuum without hBN-encapsulation, yield the band structure shown in Fig. 3a and b for 2H and 3R stacking, respectively. Due to the different atomic registries for top and bottom layers in 3R stacking, the band structure for these two layers is distinct, with different transition energies at the K-point. This is in stark contrast to 2H stacking, where both layers, and their associate intralayer states, are energetically degenerate.
For direct comparison of our theory with PL and DR experiments, we calculated the absorbance of momentum-direct exciton transitions for both stackings. For the 2H stacking we found only one transition at about 1.75 eV, which corresponds to the A-exciton or the X transition in the spectra of Fig. 2a. Note that the feature at 1.85 eV with very low oscillator strength is the interlayer KK exciton with a much lower oscillator strength here than in BL MoSe 2 and MoS 2 due to the larger spin orbit coupling in WSe 2 [9,24,62,63]. As in the experiment, we find in our calculated absorption for the 3R stacking two transitions around 1.75 eV which are separated by about 20 meV. From our calculations we infer that the two transitions correspond to A-excitons in the top and bottom layer. This insight provides a microscopic explanation for the splitting we observe for the 3R stacked sample in Fig. 2a. Note that a similar lifting of the intralayer exciton degeneracy was anticipated for 3R MoS 2 BLs but not clearly discernible in the experimental data [24]. The interlayer exciton feature is completely absent in 3R stacking as hybridization of electronic states at the K-point is forbidden by symmetry [9]. Lower-energy momentum-indirect transitions that dominate BL emission [64] are not captured by our absorbance calculation due to the vanishingly small oscillator strength.
To complement our theory analysis of excitons in 2H and 3R WSe 2 BLs, we determined the exciton gfactors from magneto-luminescence measurements and correlated them with theoretically calculated values in different spin, valley and layer configurations [28,65]. Out-of-plane magnetic fields break the time-reversal symmetry and lift the valley degeneracy between σ + and σ − polarized PL peaks with energies E + and E − to induce a valley Zeeman splitting ∆ = E + − E − = gµ B B proportional to the exciton g-factor, the Bohr magneton µ B , and the magnetic field B [66][67][68][69]. Polarization-resolved magneto-luminescence spectra of 2H and 3R WSe 2 BLs, recorded at B = −8 T with σ + and σ − circular detection under linearly polarized (π) excitation, are shown in Fig. 4a and b. The exciton g-factors follow from linear fits to the data shown by solid lines in Fig. 4c and d and grouped in color according to the g-factor values.
For 2H (Fig. 4c), we found for the fundamental exciton peak X the characteristic g-factor of −4.3, as well as g-factors of momentum-indirect QK excitons (i. e.   Coulomb-correlated states of the conduction band electron in the Q-valley and the valence band vacancy in the K-valley) in the range from −9.3 to −10.9 (peaks P 1 through P 4 ), consistent with previous findings for exfoliated WSe 2 BLs in 2H stacking [28]. In 3R stacking, our high-quality sample with spectrally narrow PL emission allows us to quantify two distinct g-factors of −4.3 and −4.8 for the momentum-direct KK exciton doublet X 1 and X 2 (Fig. 4d). For the emission peaks at lower energy between 1.50 and 1.55 eV, we determined g-factors of about −11 as in 2H (peaks P 1 through P 4 ) and an additional group of peaks (peaks P 1 through P 3 ) with a g-factor close to −14. According to our understanding of the band structure and exciton g-factors in 2H WSe 2 BL, the peaks P 1 through P 4 stem from phonon-assisted recombination of QK excitons, without contribution of QΓ states with small g-factors between 0 and 4 [28].
To interpret the origin of additional peaks in 3R stacking with distinct magneto-luminescence features as well as their difference to the 2H configuration, we calculated the exciton g-factors from first principles [28,65,[70][71][72]. The results for both stackings are presented in Table I and pictorially also in Fig. 4e and f. In the case of 2H, the g-factor of X corresponds to spin-like (↑↑) KK transition, whereas momentum-indirect exciton peaks (P 1 to P 4 ) relate to spin-like intralayer QK and interlayer Q K transitions with g = −9.51. For transitions in 3R, our results also identify spin-like KK excitons X 1 and X 2 with two different g-factors of −4.18 and −3.63. This difference in the theoretical values arises from different g-factors of the conduction bands at K in the top and bottom layers L 1 and L 2 .
Momentum-indirect excitons with peaks P 1 to P 4 in 3R with similar g-factors to the corresponding peaks in 2H stem from spin-like intralayer and interlayer QK transitions with g = −8.86 and g = −9.78, respectively. Additional 3R peaks with relatively large g-factors (P 1 to P 3 ), on the other hand, originate from spin-unlike intralayer and interlayer Q K transitions with g-factors of −13.80 and −13.76. We note that in 2H, spin-unlike intralayer Q K and interlayer QK transitions with g = −13.71 from Table I are not observed within the signal-to-noise ratio of our experiment. These states are degenerate with the two respective spin-like interlayer Q K and intralayer QK transitions, yet their exciton population seems to be vanishingly small, presumably due to spin-conserving re-TABLE I. Theoretical exciton g-factors in 2H and 3R WSe2 bilayers. We distinguish between intra-and interlayer excitons and show for momentum-direct KK transitions only the g-factors of intralayer excitons. For momentum-indirect QK and Q K transitions, we select spin-valley configurations of intra-and interlayer excitons with lowest energies. The transition column (L cb → L vb ) indicates the layers with conduction and valence band electrons involved in phononassisted recombination.
Exciton g-factors laxation [10]. However, once the degeneracy of the states is removed by the lack of inversion symmetry in 3R stacking, our analysis shows that spin-unlike Q K excitons light up in the PL of peaks P 1 to P 3 .

III. CONCLUSIONS
Our extensive experimental and theoretical studies of BL WSe 2 identify stacking-dependent optical response of 2H and 3R crystals. We find the marked differences to arise from distinct multiplicity of intralayer and interlayer excitons in the two different stable realizations of BL stackings. Using high-quality samples of both stacking configurations by CVD crystal synthesis and hBN-encapsulation, we obtain spectrally narrow resonances which allow us to identify in great detail the characteristics of momentum-direct excitons in absorption and momentum-indirect excitons in emission of 2H and 3R BLs. For both stackings we find that the lowestenergy PL is dominated by momentum-indirect QK and Q K excitons with intralayer and interlayer character. The different degrees of layer delocalization give rise to two different electrostatic dipoles, providing an explanation for the recent observations of two distinct slopes in the first-order Stark effect [20,21]. Moreover, our results demonstrate that exciton state multiplicity, transition energies and oscillator strengths sensitively depend on the BL stacking order, and highlight how layerdegenerate excitons in inversion symmetric 2H BLs split into exciton doublets in 3R stacking with broken inversion symmetry and different atomic registries for the top and bottom layers. Our findings are relevant for the un-derstanding of other TMD homobilayer systems as well as semiconductor van der Waals heterobilayers with varying degrees of symmetry, twist and interlayer hybridization.

IV. METHODS
Synthesis of homobilayer WSe 2 : Homobilayer WSe 2 crystals were synthesized on thermally oxidized silicon substrates (with a SiO 2 thickness of 285 nm) using WO 2 and Se powders (99.99%, Sigma-Aldrich) as precursors and the vapor phase chalcogenization method to obtain high-quality crystalline samples [73,74] with addition of NaCl [75] for higher yield. A three-zone furnace CVD system (Carbolite Gero) equipped with a 1-inch quartz tube was used for the growth. An alumina boat containing mixed powder of WO 2 (40 mg) and NaCl (5 mg) was placed at the center of the first zone and a SiO 2 /Si substrate was located face-down above the WO 2 powder. Another crucible boat containing Se was placed 25 cm upstream from the center of the first zone. After the tube was evacuated to ∼ 10 mTorr several times to remove air and moisture, the reaction chamber pressure was increased to ambient pressure through 500 sccm argon gas flow. Then, the furnace was heated with a ramping rate of 50 • C/min to the target growth temperature and kept there for 5 min with 140/10 sccm Ar/H 2 gas flow before cooling down. In general, the yield of 3R stacked crystals was highest at 960 • C, while 2H stacking was obtained at higher temperatures.
Sample fabrication: A PDMS/PC stamp was used to sequentially pick up the exfoliated hBN layers (NIMS) and CVD-grown homobilayer WSe 2 crystals using the dry transfer method [76,77]. Poly-(Bisphenol A-carbonate) pellets (Sigma Aldrich) were dissolved in chloroform with a weight ratio of 8. The mixture was stirred at 500 rpm using a magneton bar at room temperature overnight. The well-dissolved PC film was mounted on a PDMS block on a glass slide. First, the top hBN layer with a thickness of 163 nm for 2H (161 nm for 3R) was picked up with the stamp, followed by 2H or 3R WSe 2 BL and the bottom hBN layer with a thickness of 67 nm for 2H (66 nm for 3R). The pick-up temperatures for the hBN flakes and WSe 2 BLs were around 50 • C and 140 • C, respectively. The entire stack was released at a temperature of 180 • C onto a Si/SiO 2 target substrate, then soaked in chloroform solution for 20 min to remove PC residues, cleaned by acetone and isopropanol and annealed at 200 • C under ultrahigh vacuum for 15 hours.
Optical spectroscopy: Raman spectra were recorded at room temperature under ambient conditions using a Raman system (T64000, Horiba) with 40 µW laser excitation at 532 nm and a spectral resolution of 0.3 cm −1 with a 2400 grooves/mm grating. Low-temperature PL and DR measurements were performed with a micro-spectroscopy set-up assembled around a closed-cycle magneto-cryostat (attocube systems, attoDRY1000) equipped with nanopositioners (attocube systems, ANP101 series), a low-temperature apochromatic microobjective (attocube systems, LT-APO/VISIR/0.82) and a bi-directional solenoid for magnetic fields of up to 9 T in Faraday configuration. A halogen lamp (HL 2000, Ocean Optics) or pulsed supercontinuum laser (NKT, SuperK Extreme) was used for reflectivity measurements. Differential reflectivity spectra were obtained by normalizing the reflected light intensity from the hBN encapsulated homobilayer region to that from a bare hBN region. For PL measurements, the sample was excited at powers ranging from 10 to 50 µW either by a continuous-wave laser diode at 670 nm or pulsed supercontinuum laser with a 10 nm bandwidth around 670 nm. The spectra were detected with a CCD cooled by liquid nitrogen (Roper Scientific, Spec 10:100BR/LN) dispersed with a monochromator (Roper Scientific, Acton SP2500).

Theoretical calculations:
The atomic structures, the quasi-particle band structures and optical spectra have been obtained using the VASP package [78,79]. Core electrons have been treated by the plane-augmented wave scheme [80,81]. A lattice parameter value of 3.32 Å has been set for all calculation runs. A grid of 15 × 15 × 1 k-points has been used, in conjunction with a vacuum height of 21.9 Å, for all the calculation cells. The geometry's optimization process has been performed at the PBE-D3 level [82] in order to include van der Waals interaction between layers. All the atoms were allowed to relax with a force convergence criterion below 0.005 eV/Å. Heyd-Scuseria-Ernzerhof (HSE) hybrid functional [83][84][85] has been used as approximation of the exchangecorrelation electronic term, including SOC, to determine eigenvalues and wave functions as input for g-factor calculations and the full-frequency-dependent GW calculations [86] performed at the G 0 W 0 level. An energy cutoff of 400 eV and a gaussian smearing of 0.05 eV width have been chosen for partial occupancies, when a tight electronic minimization tolerance of 10 −8 eV was set to determine with a good precision the corresponding derivative of the orbitals with respect to k needed in quasi-particle band structure calculations. The total number of states included in the GW procedure is set to 1280, in conjunction with an energy cutoff of 100 eV for the response function, after a careful check of the direct band gap convergence (smaller than 0.1 eV as a function of k-points sampling). Band structures have been obtained after a Wannier interpolation procedure performed by the WAN-NIER90 program [81]. All optical excitonic transitions have been calculated by solving the Bethe-Salpeter equation [87,88], using the twelve highest valence bands and the sixteen lowest conduction bands to obtain eigenvalues and oscillator strengths on all systems. From these calculations, we report the absorbance values by using the imaginary part of the complex dielectric function constructed with a broadening of 12 meV.

V. ACKNOWLEDGEMENTS
We thank P. Altpeter and C. Obermayer for assistance in the clean room and the group of T. Weitz for providing access to their Raman spectrometer. Note: During the preparation of our manuscript we became aware of a related work on stacking-dependent optical properties in bilayer WSe 2 [89].