Probing the Role of Interlayer Coupling and Coulomb Interactions on Electronic Structure in Few-Layer MoSe2 Nanostructures

Despite the weak nature of interlayer forces in transition metal dichalcogenide (TMD) materials, their properties are highly dependent on the number of layers in the few-layer two-dimensional (2D) limit. Here, we present a combined scanning tunneling microscopy/spectroscopy and GW theoretical study of the electronic structure of high quality single- and few-layer MoSe2 grown on bilayer graphene. We find that the electronic (quasiparticle) bandgap, a fundamental parameter for transport and optical phenomena, decreases by nearly one electronvolt when going from one layer to three due to interlayer coupling and screening effects. Our results paint a clear picture of the evolution of the electronic wave function hybridization in the valleys of both the valence and conduction bands as the number of layers is changed. This demonstrates the importance of layer number and electron–electron interactions on van der Waals heterostructures and helps to clarify how their electronic properties might be tuned in future 2D nanodevices.

Of these, three stacking configurations, referred to as AB-stacked, have an inversion center and two configurations, referred to as AA-stacked, do not. There are two configurations (AA1 and AB1) where a Mo atom in the top layer is above a Se atom in the bottom layer. There is one configuration (AB2) where a Se atom in the top layer is above the hexagon center of the bottom layer, and two configurations (AA3 and AB3) where a Se atom in the top layer is directly above a Se atom in the bottom layer. To determine the stacking configuration of our samples, we performed fully relativistic density functional theory (DFT) calculations using the local density approximation (LDA) and scalar relativistic DFT calculations using both LDA and the van der Waals density functional with Cooper exchange (vdW-DF-c09x) 1, 2 . We find that AA1 and AB1 are the most stable stacking configurations.
AA1 and AB1 also minimize the interlayer distance, suggesting that steric repulsion effects determine the stability of MoSe 2 stacking configurations, as is also the case in MoS 2 3 . The other stacking configurations, AA3, AB2, and AB3, are between 0.03 and 0.09 eV/(unit cell) higher in energy than AB1 and AA1. The Kohn-Sham band gaps and total energies of the different MoSe 2 stacking configurations are reported in Table 1.
Since AA1 has no inversion symmetry, spin-orbit coupling splits the valence bands near the K point into four bands that are roughly equally spaced (Fig. S1f). In the bandstructure of AB1 (Fig S1g), however, time reversal and inversion symmetries ensure that the four highest valence bands become 3 two doubly degenerate bands, which are separated by 0.23 eV. Our STM data ( Fig S2b) and previously published ARPES measurements of our bilayer MoSe 2 sample 4 show two valence bands at the K point separated by 0.2 eV. This strongly suggests that our sample has AB1 stacking. AB1 is also the bulk stacking configuration of 2H-MoSe 2 . Hence, all our GW calculations are done on AB1-stacked MoSe 2 .  Bandstructure of AB1-stacked bilayer MoSe 2 from fully-relativistic LDA calculation.

Theoretical dI/dV Calculation:
We calculated the LDOS and used it to simulate the dI/dV curve within the Tersoff-Hamann formalism 5 . For each MoSe 2 structure, we used a spherically localized tip 4Å above the top Se atom (center-to-center). The tip position was averaged in a plane parallel to the MoSe 2 surface to reduce numerical noise.

Electronic Structure Calculations:
We performed mean-field density functional theory (DFT) calculations in the local density approximation (LDA) using the Quantum Espresso code 6 . The calculations were done in a supercell arrangement with a plane-wave cutoff basis using norm-conserving pseudopotentials with a 125 Ry 6 wave function cutoff. We included the Mo semicore 4d, 4p and 4s states as valence states for our DFT and GW calculations. The distance between repeated supercells in the out-of-plane direction is 50 Å.
We fully relaxed the few-layer MoSe 2 structures including the bilayer graphene substrate, and included spin-orbit interactions as a perturbation as in Ref. ( 7 ). We calculated the substrate screening due to the bilayer graphene following the same procedure as in Ref. ( 8 ).
As mentioned in the main text, we found the GW quasiparticle bandgap to converge very slowly with respect to the number of k-points in the Brillouin zone. This slow convergence originates from the large lattice constant in the out-of-plane direction, which leads to a sharp variation in the dielectric matrix ߝ ۵۵ ᇲ ሺ‫ܙ‬ሻ for small wave vectors ‫.ܙ‬ A Monkhorst-Pack grid of at least a 60x60x1 kpoints is necessary to converge the quasiparticle bandgap to within 100 meV. In order to address this computational challenge, we employed a non-uniform sampling scheme for the Brillouin zone. We first set up a regular 6x6x1 Monkhorst-Pack k-grid, and we then overlaid an extra set of 10 q-points around Γ to accurately capture the long wave length behavior of the screening for few-layer MoSe 2 . This nonuniform approach yields the same gaps as those obtained on regular 60x60x1 k-grids, but at a small fraction of the computational cost.
We estimate our quasiparticle band structures to be converged to within 100 meV for our GW calculations. The error arises mainly from the cutoff of the dielectric matrix (20 Ry) and number of bands (5000) included in the expression for the self-energy operator. We estimate the error as roughly 150 meV for our GW calculations that include substrate screening. The extra error of ~ 50 meV is estimated from neglecting hybridization between MoSe 2 and BLG, as well as from the discretized frequency sampling employed to calculate the frequency-dependent dielectric response from BLG.