Twist-tailoring Coulomb correlations in van der Waals homobilayers

The recent discovery of artificial phase transitions induced by stacking monolayer materials at magic twist angles represents a paradigm shift for solid state physics. Twist-induced changes of the single-particle band structure have been studied extensively, yet a precise understanding of the underlying Coulomb correlations has remained challenging. Here we reveal in experiment and theory, how the twist angle alone affects the Coulomb-induced internal structure and mutual interactions of excitons. In homobilayers of WSe2, we trace the internal 1s–2p resonance of excitons with phase-locked mid-infrared pulses as a function of the twist angle. Remarkably, the exciton binding energy is renormalized by up to a factor of two, their lifetime exhibits an enhancement by more than an order of magnitude, and the exciton-exciton interaction is widely tunable. Our work opens the possibility of tailoring quasiparticles in search of unexplored phases of matter in a broad range of van der Waals heterostructures.

For some of the samples in Supplementary Fig. 1, the twist angle could be determined in an alternative way. For example, the structure with a twist angle of = 60° is a pristine, as-exfoliated BL. The remaining samples were produced starting from an extremely large exfoliated monolayer. Then, only roughly half of this monolayer was transferred onto the substrate. This was achieved by rapidly retracting the transfer stage from the substrate during the transfer process, effectively splitting the monolayer into two pieces. Thus, the part of the monolayer, which remained on the transfer substrate, was perfectly aligned with the monolayer on the diamond substrate. By a precise rotation of the diamond substrate, the twist angle could then be controlled in situ without resorting to polarization-resolved second-harmonic generation.
Supplementary Figure 2 | Determination of the twist angle by second harmonic generation. Component of the second-harmonic intensity I2,|| polarized parallel with respect to the linear polarization of the laser excitation. The characteristic six-fold patterns indicate the armchair directions of the bottom (blue dots) and top (orange dots) monolayers. Additionally, the second-harmonic intensity recorded on the BL region (black dots) is shown. The strong destructive interference attests to an alignment close to the natural 2H configuration whereas constructive interference implies an alignment close to the 3R configuration, from which twist angles of  = 53° (a),  = 59° (b), and  = 21° (c) can be deduced.

Supplementary Note 2
The experimental details of the ultrafast near-infrared (NIR) pump -mid-infrared (MIR) probe spectroscopy scheme are summarized in the Methods section of the manuscript. Supplementary Fig. 3 displays the NIR pump spectrum and the electric field of the MIR probe transient with its spectral amplitude and phase.
Supplementary Figure 3 | Ultrafast NIR pump-MIR probe spectroscopy. a, Electric field waveform of the MIR probe pulses as a function of the electro-optic delay time teos. Inset: Spectral intensity of the NIR pump pulses as a function of the photon energy. b, Spectral amplitude and phase of the MIR probe pulses as a function of frequency.

Supplementary Note 3
The twist-dependent delocalization of the excitons also manifests as an angle-dependent dynamical shift of the 1s-2p resonance. The dielectric responses of two BLs with  = 60° and  = 0° are shown in Supplementary Fig. 4a,b, for two representative pump-probe delay times (tpp = 0 ps and 10 ps). For tpp = 0 ps, a maximum in 1 at an energy of ħres = 93 meV arises for both samples. After 10 ps, clear deviations between the dielectric responses of the samples are discernible. For  = 60°, the resonance energy slightly blue shifts by ~7 meV, whereas for  = 0° a redshift of ~30 meV occurs. Extracting the 1s-2p resonance energies of the dielectric function for a larger set of pump-probe delay times and twist angles yields the complete temporal evolution of the internal transitions, as depicted in Supplementary  Fig. 4c. The resonance energy increases/decreases as a function of the pump-probe delay time for a twist angle of  = 60°/ = 0° (gray/orange spheres). We attribute this blue/red shift to the formation of K- excitons -subjected to different degrees of hybridization -after initial photoexcitation at the K valley.
For  = 60°, the resulting K- excitons feature a higher binding energy; the 1s-2p transition energy therefore increases as compared to the initially created K- excitons owing to the larger electron effective mass in the  valley 2,3 . However, for  = 0° the 1s-2p transition energy of the K- excitons shows a clear red shift as compared to transition energy observed in  = 60°. Furthermore, for  = 59°, the evolution of the 1s-2p transition energy is similar to  = 60°, where the hybrid K- exciton transition energy slightly increases due to the reduced interlayer hybridization. For intermediate twist angles of  = 53°, 30°, and 5°, the changes to the intraexcitonic resonance are less pronounced, but can be well explained by our microscopic theory. the fits to those data with a phenomenological model. The red/blue dashed line indicates the 1s-2p resonance of the K-K and K- exciton, respectively. c, 1s-2p resonance energy ħres as a function of the pump-probe delay time tpp for different twist angles . The samples were kept at a temperature of 5 K. The error bars represent the 95% confidence interval of the fitting procedure.

Supplementary Note 4
To investigate the density-dependent renormalization of the hybrid exciton binding energy, we recorded the dielectric response function of WSe2 BLs with twist angles of  = 60°, 5°, and 0° at different pump fluences using a fixed pump-probe delay time of tpp = 5.1 ps ( Supplementary Fig. 5). For all samples and pump fluences, we observe a clear maximum in the real part of the optical conductivity 1 and a corresponding dispersive feature in the real part of the dielectric function 1, arising from the resonant response of the intraexcitonic 1s-2p transition. The characteristic spectral features of the pristine BL and the sample with a twist angle of  = 5° remain at an energy of ~100 meV for all pump fluences. For  = 60°, no significant pump-probe signal was observed for the lowest pump fluence of  = 7 µJcm -². Employing our phenomenological model to fit the experimental spectra (see Methods), we quantify the renormalization of the 1s-2p resonance energy ħres (Fig. 4c) as a function of the pump fluence While ħres of the BL with 2H stacking and the one with  = 5° are hardly affected and remain nearly unchanged at ~100 meV, the pump fluence has a dramatic impact for  = 0°. Here the 1s-2p resonance energy drops from ħres = 83 meV to 63 meV as the pump fluence is increased from 7 µJcm -² to 27 µJcm -². These observations are in line with our microscopic understanding of the twist angle-dependent hybridization of excitonic states as discussed in the main text.
Here, we present the theoretical approach used to determine eigenenergies and wavefunctions of the energetically lowest lying excitonic states in the BL. Since K- excitons are known to be the energetically lowest states in the WSe2 BL(ref. 2 ), we start with a BL Hamilton operator restricted to electrons and holes at these high symmetry points. Here, the electronic contribution reads with annihilation (creation) operators of conduction k ( ) and valence k ( ) band electrons with the momentum k and in the layer = 1,2. Here, the eigenenergies of the layer-localized basis functions are given by the effective mass approximation k λ =  λ ± ℏ /(2 λ ) with effective masses λ of the  valley electrons and K valley holes. Effective masses, bandgaps and energetic separations of valleys entering 0 λ are taken from DFT calculations 3 . The second term of Eq. 1 contains the Coulomb interaction between electrons and holes with the matrix elements q . The latter are derived in a semiclassical approach assuming anisotropic homogeneous slabs to account for the dielectric environment created by the transition metal dichalcogenide (TMD) BL system 3 . Here, we use ab initio parameters 4 for the dielectric tensor of WSe2 and diam = 5.8 (ref. 5 ). The last term of Eq. 1 reflects the overlap of electronic wavefunctions of the two layers giving rise to interlayer hybridization. We focus on the hybridization of electrons at the  valley, which is known to be orders of magnitude stronger than at the K valley 6 . The interlayer hopping is momentum conserving, which in the valley-local momentum coordinates requires a momentum transfer ℎ q = ℎ q,Λ Λ (ref. 7 This resembles the well-known anti-crossing behavior of hybridizing particles. Furthermore, for all configurations, the lowest hybrid exciton is mostly composed of intralayer states. Therefore, the energy reduction µQ hyb = µQ intra − µQ hyb, of the intralayer exciton due to delocalization (measure for the degree of hybridization) is given by µQ hyb = − µQ + Δ µQ + 4 µ , with the detuning of inter-and intralayer states Δ µQ = µQ inter − µQ intra > 0. We find that hyb > 0 and hyb / Δ < 0. Therefore, the red shift reduces monotonously, meaning that the hybridization becomes weaker with increasing detuning of the mixing states. This explains the different hybridization of states at  = 0° and 60° stacking, since the layer inverted spin-orbit coupling leads to a large separation of intra-and interlayer states at 60°. In contrast, the mixing electronic bands are degenerate at 0°. For ≠ 0° or 60°, the dispersion of intraand interlayer excitons becomes misaligned in momentum, leading to a Q-dependent detuning ( Supplementary Fig. 6). This already leads to an upshift of the hybrid dispersion (violet lines), since the interlayer state (orange line) hybridizing with the minimum of the intralayer dispersion (blue line) is moving upwards in energy and the relevant detuning increases. There are two more effects leading to a reduced exciton hybridization in twisted BLs: (i) As a result of stacking dependent interlayer distances, the electronic hopping integral ℎ decreases for twisted samples 7,9 (cf. section 6). (ii) The excitonic hopping integral µ (dΛ) ∝ µ (dΛ) = ∑ µ inter ( ) * µ intra ( + dΛ) decreases with increasing momentum mismatch |dΛ| = |Λ − Λ 2 |~.
This reflects the decreasing in-plane overlap of intra-and interlayer exciton wavefunctions. The phenomenon (ii) explains why the 1s-2p transition energy increases when hybridization decreases. The large 2p state orbitals correspond to narrower momentum space wavefunctions compared to the 1s ground state (see Fig. 3c of the main manuscript). Therefore, the in-plane overlap µ decreases faster with the twist angle for the excited states than for the ground states, i.e. | 1 / | < | 2 / |. Hence, the 1s-2p transition energy increases, since the 2p state blue-shifts with a larger rate than the 1s state.
Higher-order corrections to the energy levels of excitons are expected to have a negligible effect on our current theoretical analysis, for the following reasons: (i) The expected splitting between the p+ and p-states lies in the range of 10 -25 meV for K-K excitons in monolayer TMDs in vacuum as derived in the literature [10][11][12] . Recently, this splitting was also confirmed experimentally 13 and an energy difference of ~14 meV was reported, in line with the theoretical predictions. Owing to the enhanced dielectric screening experienced by excitons in a bilayer TMD on a diamond substrate, we expect that the corresponding Berry curvature-induced splitting should be even smaller. Specifically, exciton binding energies on the order of 200 meV suggest energy differences between 2p+ and 2p-states within the range of only 5 -13 meV.
(ii) The considerations so far have neglected the fact that the electrons of the energetically favourable excitons in WSe2 reside in the -valley, where the Berry curvature is significantly smaller 14 than in the K-valley. Hence, we expect that the energetic difference between the 2p+ and 2p-states of the hybrid excitons in WSe2 homobilayers should lie well below 10 meV.
On the experimental side, we cannot currently resolve such higher-order corrections. First, the midinfrared probe pulses are linearly polarized for the extraction of the pump-induced dielectric response of the homobilayers. Since linearly polarized light is a superposition of right-and left-circularly polarized light, we simultaneously drive 1s-2p+ and 1s-2p-transitions. With this approach, we effectively obtain an averaged 1s-2p transition energy. Second, the experimental line width of the internal excitonic transitions amounts to tens of meV, which obscures the inhomogeneous broadening induced by the p+/p-splitting. Therefore, even without Berry curvature-induced corrections to the excitonic energy levels, our effective-mass model already effectively captures the essential mechanism responsible for the renormalization of the exciton binding energy with twist angle: the hybridization of intra-and interlayer excitons.
Calculating the energy levels of hybrid excitons in WSe2 homobilayers at twist angles close to  = 30° is not feasible at this moment because the interlayer hopping cannot be extracted from the single-particle bandstructure in this case. Consequently, we have decided to limit our theoretical consideration to small misalignments from 0° and 60° for the following reasons: The calculation of interlayer-intervalley excitons in twisted TMDs based solely on DFT is extremely challenging and has so far not been demonstrated yet. Instead, we use a combination of DFT and an effective-mass Wannier model to predict intervalley exciton energies in twisted homobilayers. In perfectly aligned bilayers, the initially degenerate bands of the constituent monolayers are split depending on the strength of the interlayer hopping. For twisted bilayers, this does, not hold because zone folding, for example, gives rise to a mixing of multiple bands. Consequently the band splitting is no longer a good measure for the interlayer hopping. Therefore, we use the DFT bandstructure of perfectly aligned bilayers and assume that the hopping term extracted from the band splitting to a first approximation only changes slowly with twist angle in the vicinity of a given high symmetry point and can therefore also be used for slightly misaligned samples. This approximation is, however, limited to small twist angles. At large twist angles such as 30°, the hybridization between different high-symmetry points such as  of one layer and  (which lies between  and ') of the other layer will strongly affect the interlayer hopping.
The electronic band structures of the WSe2 BLs were simulated on the basis of density functional theory in the Perdew-Burke-Ernzerhof approximation, as implemented in the Quantum ESPRESSO package 15 .
Integrations in reciprocal space were performed on a Monkhorst-Pack grid of 12x12x1 points in the Brillouin zone. The core electrons were replaced by normconserving pseudopotentials from the Pseudo Dojo library 16 , where the s and p semi-core electrons of tungsten and the semicore d electrons of selenium were included in the set of valence electrons. We used a cutoff energy of 90 Ry (1224 eV) for the employed planewave basis set, ensuring well converged ground state electron densities and wavefunctions. Using these parameters, we optimized the atomic positions and in-plane lattice constants of twisted WSe2 BLs with different relative twist angles, such as 60° (corresponding to the 2H-stacking order in bulk WSe2) and 46.8°. We used the resulting atomic geometries to derive an average equilibrium distance between the two tungsten sublayers. We included semi-empirical van der Waals corrections from the PBE+D3 (ref. 17 ) scheme for these calculations, which yield excellent predictions of the in-and out-of-plane lattice constants of layered systems 18,19 . Based on the optimized geometry of 2H-stacked BL WSe2 with an optimized interlayer distance of 6.5 Å, we subsequently simulated the evolution of the electronic bandstructure (neglecting spin-orbit-coupling effects) with the interlayer interaction by rigidly shifting the upper layer in out-of-plane direction, thus varying the interlayer distance in the range between 6.5 Å and 7.25 Å. The interlayer hopping parameter ℎ introduced in section 5 is extracted from the BL bandstructure without spin-orbit interaction. In particular, the interlayer hybridization gives rise to a k-dependent splitting of the electronic bands Δ k , allowing to determine the hopping integral ℎ = Δ /2 (ref. 6 ). We find e.g. ℎ ≈ 230 meV for the interlayer distance of 6.5 Å ( = 60°) and an almost linear dependence on interlayer distance, giving rise to a reduced hopping of about ℎ ≈ 200 meV at 6.75 Å ( = 46.8°) as the result of a reduced orbital overlap. To obtain a continuous twist angle dependence of the hopping parameter, the interlayer distances obtained for the computationally feasible angles are interpolated linearly for the intermediate angles.
To obtain the full bandstructure for 2H-stacked WSe2, we combined spin-orbit coupling and quasiparticle energy corrections from the G0W0 approximation. Here, we first used the YAMBO code 20 to calculate the quasiparticle corrections with SOI effects on a regular 9x9 k-point grid. The dynamic dielectric function was represented on a discrete 9×9 q-point grid of transferred momenta in reciprocal space and by the Godby-Needs plasmon-pole approximation in frequency space. Local field effects and exchange interactions were included up to a cutoff energy of 300 eV and 650 eV, respectively. We truncated the Coulomb interaction in the non-periodic direction through a Wigner-Seitz cutoff scheme 18 . The correct asymptotic behavior of the dielectric screening for transferred momenta close to the  point was imposed through an isotropic model dielectric function that was used previously to great success 18,21 . While it is sufficient to capture the correct quasiparticle band dispersion, it causes an overestimation of the electronic band gap sizes. We thus performed additional GW calculations (neglecting spin-orbit interactions) with the BerkeleyGW code 22 , where we used the recently proposed Nonuniform Neck Subsampling method 23 to derive an additional scissor shift correction for the electronic band gaps. The electronic band structures were then derived from the corrected electronic bands through Wannier interpolation.