Tuning valleys and wave functions of van der Waals heterostructures by varying the number of layers: A first-principles study

In van der Waals heterostructures of two-dimensional transition-metal dichalcogenides (2D TMDCs) electron and hole states are spatially localized in different layers forming long-lived interlayer excitons. Here, we have investigated, from first principles, the influence of additional electron or hole layers on the electronic properties of a MoS2/WSe2 heterobilayer (HBL), which is a direct band gap material. Additional layers modify the interlayer hybridization, mostly affecting the quasiparticle energy and real-space extend of hole states at the G and electron states at the Q valleys. For a sufficient number of additional layers, the band edges move from K to Q or G, respectively. Adding electron layers to the HBL leads to more delocalized Q states, while G states do not extend much beyond the HBL, even when more hole layers are added. These results suggest a simple and yet powerful way to tune band edges and the real-space extend of the electron and hole wave function in TMDC heterostructures, strongly affecting the lifetime and dynamics of interlayer excitons.

For instance, Tongay et al. [22] demonstrated that the interlayer coupling of TMDC vdW HS can be tuned by vacuum thermal annealing, where light emission spectrum gradually evolves from single layers, contributing separately to the novel coupled spectrum as function of the interlayer distance.
Up to date, the majority of research focuses on vdW HS made of just two different TMDC layers. However, band structure engineering by systematically extending the number of layers of the heterostructure has been much less explored. [37,38] Prominent questions to be addressed are: (i) Can we exploit interlayer coupling by mixing different number of layers to tune the multi-valley, electronic structure near the band edges? (ii) Will additional electron or hole layers delocalize or localize the interlayer excitons in the vdW HS? (iii) For how many layers of the multi-stack, would the band gap be direct or indirect?
Here, in order to address these questions, we investigate the influence of additional electron or hole layers (nL, with n = 1-4) on the electronic properties of a MoS2/WSe2 heterobilayer (HBL).
Our results show that adding additional WSe2 layers to the HBL (forming 1L-MoS2/(n+1)L-WSe2 HS) will mostly affect the valence (hole) states, while adding MoS2 layers to the HBL (forming (n+1)L-MoS2/1L-WSe2 HS) will affect the conduction (electron) states. We find that the heterostructures of up to 1L-MoS2/4L-WSe2 or 3L-MoS2/1L-WSe2 stay direct gap semiconductors (at the K point), while more layers in both HS would open new transition channels at the Q and/or G points. We, therefore, provide the energy differences between k valleys in the conduction (Q-K) and the valence (K-G) bands, values which can be obtained experimentally, e.g., from carrier dynamics measurements. [21] The calculated partial charge densities (PCD) for the electron and hole states, analyzed for each relevant k point of the band edges, are either delocalized throughout all the layers or localized in the heterolayer, depending on the stoichiometry of vdW HS. Here, we apply a simple two-step approach of analyzing electron and hole states in large incommensurate HS using single unit cells, where we minimize the effect of the tensile strain or compression, resulting from the differences in the lattice vectors of constituent layers.

Results and Discussion
First, electronic of bulk and 1L MoS2 and WSe2 were obtained (see Figure S1 in Supporting Information (SI)). Good agreement with previous works supports the choice of our simulation methods. [39,40] While 1L-MoS2 is a direct band gap semiconductor (at K), 1L WSe2 is an indirect gap material, [41] with the corresponding direct band gap at K being only 50 meV higher in energy. Both 1L systems exhibit spin splitting in the valence band at the K point of 197 meV and 606 meV for MoS2 and WSe2, respectively, due to strong SOC. It is well known that HSE06 overestimates the spin split for TMDCs. [42] Next, we have calculated band structures of vdW HS, HBL/nL-MoS2 and HBL/nL-WSe2 (with n = 1-4; see Figure 1), using the fully relaxed structures (method I) and our proposed 2-step approach (method II) for comparison (see Computational Details in Section 4 for detailed explanations on employed methods). In method II, we perform two different calculations to analyze the electron and hole states separately. By using different lattice constants for the electron and hole calculations, our approach removes the problematic tensile strain or compression, which arises from the incommensurability of the layers. This approach is valid for our systems, because they are type-II HS. Because with method I both " " and " $ stackings result in very similar electronic structures (see Figure S2 to Figure S5 in SI) and energy differences, we have only used the " " systems for method II. Here, we first shortly discuss the results from approach I and then more detailed results are shown for the approach II. The stacking polytype of the HBL is " " or " $ and the one of the nL and its interface with HBL is " " . (b and c) Top and side views of the considered stacking polytypes ( " " , also known as 2H, and " $ , also known as 3R) in HBLs with 0° and 60° twist angles.
While stacking results in very similar band structures, more significant differences are observed for the multi-stacks: 1) the band gap increases (decreases) for HBL/nL-MoS2 (HBL/nL-WSe2); 2) when n is varied for HBL/nL-MoS2 (HBL/nL-WSe2), the conduction (valence) bands are    The most significant difference between approach I and II is how fast the Q valley in the CB (G valley in the VB) decreases (increases) in energy. In other words, the energy differences between the Q and K valleys in the CB and the G and K valleys in the VB, indicating direct or indirect band gaps, differ for both approaches (see Figure 5). Approach I would indicate that both types of HS with at least 6 layers in the multi-stack are still direct-gap materials. However, approach II shows that these changes happen much faster and already 4L-MoS2/1L-WSe2 and 1L-MoS2/5L-WSe2 transform to indirect-gap materials. Additionally, we have calculated the PCD of electron and hole states (cf. Figure 3 and 4). We  [21] In present work, we avoid discussing the absolute band gap values. Our 2-step approach can also be used to analyze band gap values when the band edges are calculated with respect to the vacuum energy. From approach I, we have noticed that in HS with nL-MoS2, the band gaps increase, while in HS based on nL-WSe2, the band gaps decreased with n (cf. Figure S2 to S6).
This subject, however, requires more detailed investigations and is beyond the scope of this work.

Conclusion
We have investigated the influence of additional electron or hole layers (nL, with n = 1-4) on the electronic properties of a MoS2/WSe2 heterobilayer (HBL). We showed that this approach can tune the valleys and the wave functions of the multi-stack in different ways, depending on the stoichiometry of the layers. We have studied two types of systems, namely (n+1)L-MoS2/1L-WSe2 and 1L-MoS2/(n+1)L-WSe2. Furthermore, we compared a very common approach, where the heterostructure is modelled with small, but strongly strained, unit cells with our simple 2-step approach, where two calculations are performed to analyze the electron and hole states separately. By using different lattice constants for the electron and hole calculations, our approach removes the problematic tensile strain or compression, which arises from the incommensurability of the layers. It can be used, because the HS are type-II materials, where the electrons reside predominantly in the MoS2 layers and holes predominantly in WSe2 layers.
We have analyzed the energy differences between the band edges at the K valley and other band edge valleys, i.e., the Q valley in the conduction and the G valley in the valence band. These energy differences can be obtained experimentally, e.g., from carrier dynamics measurements.
We show that increasing the number of additional layers, the Q-K and G-K energy differences decrease and increase in energy, respectively. Our calculations suggest that the materials stay Our study offers a simple experimental approach to understand and control excitons in TMDC heterostructures.

Experimental Section
Computational Details: All simulations were performed using density functional theory (DFT) methods as implemented in the Vienna ab-initio simulation package (VASP). 39 The projectoraugmented wave (PAW) method was used to describe the interactions between electrons and nuclei. The exchange and correlation potentials were treated using the generalized-gradient approximation (GGA) 40 functional proposed by Perdew, Burke, and Ernzerhof (PBE) for full geometry optimization (lattice vectors and atomic positions). Plane-wave cutoff of 500 eV was used in all simulations. The D3 dispersion correction was included, as proposed by Grimme. 41 A vacuum layer of at least 15 Å was added along the out-of-plane direction to avoid spurious interactions with the next periodic images. All structures were relaxed until all the forces acting on each atom were less than 2´10 -2 eV Å -1 and the total energy change between two selfconsistent steps was less than 1´10 -4 eV. A G-centered 12´12´1 k-point mesh was used to sample the Brillouin zone during structural optimizations. It is well known that standard GGA DFT functionals underestimated the electronic band gaps, thus, electronic properties were calculated using the Heyd-Scuseria-Ernzerhof hybrid (HSE06) 42,43 functional, considering relativistic effects, such as spin-orbit coupling (SOC). Even though, we do not discuss the band gap sizes in the present work, we have also noticed strong deviations in the shape of bands in band structures between PBE and HSE06, thus, we used the latter for the electronic structure analysis. The PCD was, however, calculated using PBE with SOC and analyzed for exemplary structures, i.e., (HBL + 2L) 4L vdW HS. Structural illustrations and PCD plots were created using the 3D visualization package VESTA. 44 The systems studied in this work are shown in Figure 1. We consider a MoSe2/WSe2 HBL and systematically add one to four additional layers (nL) of either WSe2 or MoS2, creating either a 1L-MoS2/(n+1)L-WSe2 or a (n+1)L-MoS2/1L-WSe2 vdW HS. The (n+1)L multilayers were kept in their most stable stacking arrangement, " " (also known as 2H). The interface between MoS2 and WSe2 is realized either in the " " or " $ (also known as 3R) stacking fashion, which are among the most stable high-symmetry stacking arrangement observed in 0° and 60° twist angle moiré structures. 25 Other stacking polytypes, e.g., " $ or " " , are, of course, also possible, 21,25 but beyond the scope of the present work, because the electronic properties are affected much stronger by the number of layers than the stacking itself.
A lattice mismatch between the building blocks of HS results in an incommensurate moiré superstructure that, in the best case, would need to be approximated by supercells containing hundreds or thousands of atoms (see Figure S6). Such large systems are, however, too demanding for DFT calculations, especially when including SOC with hybrid functionals. The most commonly used approach is to relax such a HS in a small unit cell (referred to as method I in the text below), however, this results in strain or compression of the constituent layers, which is an unrealistic representation of the HS. As previously shown in different theoretical and experimental works, tensile strain and compression strongly affect electronic structures of TMDC layers. 45,46 Thus here, we propose a simple 2-step approach (referred to as method II) to separately analyze the electron and hole states of an incommensurate HS by performing two calculations within primitive unit cells. We suggest that, in order to analyze the electron states near the Fermi level, we use a lattice constant (a) of the relaxed (n+1)L-MoS2, while for the hole states, we use a of the relaxed (n+1)L-WSe2. This is possible, because of the type-II band alignment, which is discussed further below. We used the same approach in our recent work on the W-based HBLs. 47 In Table S1 in SI, we summarized the lattice parameters for (n+1)L-MoS2 and (n+1)L-WSe2, and we compare them to the lattice parameters of fully relaxed HS.  material, while 1L WSe2 is indirect (between K and Q points) band gap system, [2] with the direct gap at the K point being by 50 meV higher in energy. Figure S2. Band Structures of (n+1)L-MoS2/1L-WSe2 heterostructures with the " " stacking at the interface. Grey arrows indicate the fundamental band gaps, which are in all cases direct at the K point. Increasing number of layers in MoS2 stack affects mostly the conduction band (CB), with Q point lowering in energy, when increasing number of layers. The valence band (VB) is mostly affected at the G point, which is, however, much lower in energy than the K point. The difference in energy between the K and Q points in CB is given in red for each system. For convenience, the Fermi level was shifted to VBM at zero and indicated by horizontal dashed lines. Figure S3. Band Structures of 1L-MoS2/(n+1)L-WSe2 heterostructures with the " " stacking at the interface. Grey arrows indicate the fundamental band gaps, which are in all cases direct at the K point. Increasing number of layers in WSe2 stack affects mostly the valence band (VB), with the G point increasing in energy, when increasing number of layers. The difference in energy between the K and G points in VB is given in red for each system. For convenience, the Fermi level was shifted to VBM at zero and indicated by horizontal dashed lines. Figure S4. Band Structures of (n+1)L-MoS2/1L-WSe2 heterostructures with the " $ stacking at the interface. Grey arrows indicate the fundamental band gaps, which are in all cases direct at the K point. Increasing number of layers in MoS2 stack affects mostly the conduction band, with Q point lowering in energy with increasing number of layers. The valence band is mostly affected at the G point, which is, however, much lower in energy than the K point. The difference in energy between the K and Q points in the conduction band is given in red for each system. For convenience, the Fermi level was shifted to VBM at zero and indicated by horizontal dashed lines the K and G points in the valence band is given in red for each system. For convenience, the Fermi level was shifted to VBM at zero and indicated by horizontal dashed lines. Figure S6. Schematic representation of the moiré pattern of the MoS2/WSe2 heterobilayer with the 60° twist angle. Red, orange, and brown hexagons mark regions with high-symmetry stacking configurations. The 0° twist angle results in a similar representation with different high-symmetry stackings. Both require a few thousand of atoms in the supercell to minimize the strain, which is a result of the lattice mismatch between the individual layers.