Ab-initio study of the energy competition between
$\Gamma$ and K valleys in bilayer transition metal dichalcogenides

Herein we consider twist-free, homo- and hetero- bilayer systems with chemical formula $\textrm{MX}_{2}$ where $\textrm{M}=\textrm{Mo},\textrm{W}$ and $\textrm{X}=\textrm{S},\textrm{Se},\textrm{Te}$. For heterobilayers, we consider only combinations containing common chalcogens in the parent layers to avoid possible inaccuracies due to large lattice constant mismatches. While many bilayer stacking arrangements exist, in this work we consider only the two high symmetry stacking orders denoted as 2H and AA and displayed in Fig. 1. The bilayer AA stacking occurs when the $\hat{x},\hat{y}$ coordinates of the metal (chalcogen) atom(s) in the top layer are the same as the metal (chalcogen) atom(s) in the bottom layer Devakul et al. (2021). The 2H stacking is obtained by a $180^{\circ}$ rotation about the $\hat{z}$-axis of the top layer relative to the AA stacking case, and is predicted to be the most stable of the bilayer configurations Li et al. (2023); He et al. (2014). In this case, the metal (chalcogen) atom(s) are directly above the chalcogen (metal) atom(s). These arrangements can not be made equal by relative translations of the layers and therefore do not occur at different positions in the same moiré pattern.

We employ the Vienna ab-initio Simulation Package (VASP)Kresse and Hafner (1993); Kresse and Furthmüller (1996, 1996) to perform structural relaxation for each 2H system under the following procedure. Firstly, full relaxation in which both the volume and shape of the unit cell may vary to minimize the total energy and force is performed on the bulk structures using three functionals; the Perdew-Burke-Ernzerhof exchange-correlation functional with spin orbit coupling, both with (PBE-SO-D3) and without (PBE-SO) the van der Waals correction, and the local density approximation LDA+SO. Additionally, we include the relaxation via PBE alone (as performed in [Jain et al., 2013]) for comparison. Secondly, we construct each ’free-standing’ bilayer system from the bulk lattice constants found from the LDA+SO relaxation. The total length of the c-axis is set to 35Å(Fig. 1), to isolate the bilayer and limit the unphysical interaction between periodic unit cells. This structure is then relaxed with fixed volume using LDA-SO in VASP, allowing atomic positions to change. From these structures, the bands are computed as a function of a variety of tuning parameters. For AA systems, we only employ the LDA-SO functional, for reasons explained in Sec. IV.2

In order to comprehensively explore the energy competition between $\Gamma$ and K valleys, we perform Density Functional Theory (DFT) calculations using the full-potential Linearized Augmented Plane Wave (LAPW) method as provided by the WIEN2k package Blaha et al. (2019); Schwarz and Blaha (2003a). In addition to the standard Local Density Approximation (LDA) functional, we employ the modified Becke-Johnson (mBJ) functional to investigate the effect of electronic correlation Becke and Johnson (2006); Tran and Blaha (2009). It has been shown that the mBJ functional gives very accurate bandgaps in many transition metal oxidesZhu and Schwingenschlögl (2012); Wahila et al. (2019) and semiconductors,Borlido et al. (2019) including VO₂Paez et al. (2020); Singh et al. (2022) and monolayer TMDs,Li et al. (2013) with much less computational time than hybrid functional, or $GW$ methods. Furthermore, the Local mBJ functional has been developed Rauch et al. (2020) for accurate prediction of band gaps in systems with vacuum space, and works best for our study of the free standing bilayer TMDs illustrated in Fig. 1. The energy convergence was chosen to be 0.1 mRy while the charge convergence was set to 0.001 $e^{-}$. In all cases, we include spin-orbit interactions and employ an additional $p_{1/2}$ radial basis function, called the Relativistic Local Orbital (RLO) provided by WIEN2k, for the metal atoms to improve the basis functions, aiding in convergence. For the band structure, we use an in plane momentum $\vec{k}$-mesh of $12\times 12\times 1$, and trace a $\Gamma\rightarrow\textrm{K}\rightarrow\textrm{M}\rightarrow\Gamma$ path in the Brillouin-zone.

The interlayer separation, $h$, defined henceforth as the distance measured purely along the $\hat{z}$-axis between metal atoms in each layer (shown in Fig. 1) has a strong influence on the bilayer tunneling and as such is a key physical parameter Xu et al. (2017) that strongly influences the energy competition between the K and $\Gamma$ valleys. We therefore vary this parameter for values in the neighborhood of the relaxed interlayer separation, $h_{r}$, and examine the valley competition for each material in the 2H configuration. The total length along $\hat{z}$ (vacuum + c lattice constant) is kept fixed while the interlayer separation $h$ is varied. The magnitude of the a and b lattice constants are fixed throughout all calculations. This approach keeps the total system + vacuum unit cell volume fixed throughout the calculations, ensuring a consistent total energy.

In addition to interlayer separation, which can be adjusted by applying pressure, applied gate voltages are an experimentally accessible tuning parameter that influences the K-$\Gamma$ competition. Liu et al. (2012); Ramasubramaniam et al. (2011); Huang and Kaxiras (2016); Santos and Kaxiras (2013); Sharma et al. (2014); Zheng et al. (2016); Xu et al. (2023). To finalize our discussion on key factors contributing to TMD band engineering, we perform electronic structure calculations for these materials under various electric fields.

There are several element-related trends that can be seen in our results. To make the analysis concrete, we investigate the role of the element type in the 2H homobilayers only (Fig. 3). Firstly, we consider the trend of varying the metal atom in the $\textrm{MX}_{2}$ with the chalcogenide atom X fixed. It can be seen clearly that in this case, $\textrm{WX}_{2}$ always has a larger value of $\Delta E_{K-\Gamma}$ than $\textrm{MoX}_{2}$ does. This can be attributed to the fact that the W atom has a larger atomic number than the Mo atom, and therefore larger spin-orbit coupling. To demonstrate this point, we define $\Delta E_{SO}$ as the energy difference between the top two valence states at the K point as shown in Fig. 2, which can serve as a good estimation of the spin-orbit coupling. Fig. 5 (top) plots $\Delta E_{SO}$ as a function of the absolute interlayer separation $h$ for $\textrm{MoX}_{2}$ and $\textrm{WX}_{2}$ with $\textrm{X}=\textrm{S},\textrm{Se}$ and Te. Clearly, $\Delta E_{SO}$ depends strongly on the metal atom but weakly on the height $h$ and the chalcogenide atom X, which is expected since the spin orbit coupling is an intrinsic property of the metal atom.

Now we study the trend of varying the chalcogen atom in $\textrm{MX}_{2}$ with the metal atom M fixed. We see a monotonic trend where $\Delta E_{K-\Gamma}$ is generally increased when going from $\textrm{S}\to\textrm{Se}\to\textrm{Te}$. We note that because chalcogenide atoms X exist between metal atoms, as shown in Fig. 1, the channels going through the $p$ orbitals of X make the major contribution to the interlayer coupling. Generally speaking, as the atomic number of X increases, orbitals on the outermost shell becomes more de-localized, which leads to a stronger wavefunction overlap with $d$ orbitals of the metal atom and consequently a larger bilayer splitting at the $\Gamma$ point. We define $\Delta E_{BL}$ as the energy difference between the top two valence states at $\Gamma$ point as shown in Fig. 2, which can be a good estimation of the interlayer coupling. $\Delta E_{BL}$ as a function of the absolute interlayer separation $h$ for $\textrm{MoX}_{2}$ and $\textrm{WX}_{2}$ with $\textrm{X}=\textrm{S},\textrm{Se}$ and Te is presented in Fig. 5 (bottom). It is clearly observed that bilayer splitting depends strongly on interlayer separation, and that the equilibrium separation of each material monotonically increases as $\textrm{S}\to\textrm{Se}\to\textrm{Te}$, indicating that stronger bilayer splitting heavily influences the relaxation.

The structural relaxation is performed in VASP, following the procedure outlined in Sec. II.1. The resulting structural data is summarized in Table 1. For clarity and discussion, we additionally include a plot of the bulk a and c lattice constants in Fig. 6. There are a number conclusions that may be drawn from the relaxation procedure. Inclusion of the vacuum after LDA-SO bulk relaxation can be seen to increase the metal-to-metal atom distance in bilayers, which will reduce interlayer tunneling. The reduction of $h$ in bulk systems is expected since the van der Waals attraction between layers is expected to be stronger when every layer has two neighboring layers. If this consideration is correct, then the encapsulating layers present in most experimental systems may also play a role in determining the equilibrium layer separation. Our DFT calculations confirm this assertion as we see a boost in the K valley while comparing bulk to bilayer results. As we employ both VASP and WIEN2k, we may also comment that the LDA+SO band structures of VASP’s pseudo-potential method and the all-electron method of WIEN2k are nearly identical if the same structure is provided. ¹¹1Furthermore, we find that the total energy of the bilayers predicted by LDA+SO using WIEN2k is in general at the global minimum when the equilibrium interlayer separation is used. The energy generally grows as we consider the shifting of $\pm 0.2$ Å, with a sharper growth with positive shifting in interlayer separation.

It is clear from the step pattern in Fig. 6 that the magnitude of the lattice constants is determined mainly by the type of chalcogen in the structure. This is particularly true for the methods of PBE-SO-D3 and LDA-SO in the case of the c lattice constant. With comparison to the range of literature reported bulk values Puotinen and Newnham (1961); Coehoorn et al. (1987); Wakabayashi et al. (1975); Schutte et al. (1987); Molina-Sánchez and Wirtz (2011); Podberezskaya et al. (2001), we conclude that PBE and PBE-SO both tend to overestimate the lattice constants alone, and the inclusion of van der Waals is necessary. As expected Molina-Sánchez and Wirtz (2011), even in the absence of the van der Waals force LDA-SO agrees with PBE-SO-D3 quite well for layered systems such as the group VI TMDs, and is in some cases more accurate with respect to experiments. Hence, we employ only LDA-SO for the AA stacking. Our bilayer electronic structure calculations are carried out using constants from the LDA-SO approach, as explained in II.1

The energy competition between $\Gamma$ and K valleys is influenced by the strengths of interlayer tunneling, spin orbit coupling, and electronic correlations. Generally speaking, interlayer tunneling produces a more sizable bilayer splitting at the $\Gamma$ point pushing the top valence state up, with negligible effect at the K point. As a result, a smaller distance between the bilayers favors the valence band maxima to be at $\Gamma$ point due to the larger bilayer splitting. For heterobilayers, the symmetry between layers is broken and the bilayer splitting at the K is not negligible, which explains why the heterobilayers tend to favor the K valley when compared to the homobilayers of their parent atoms. On the other hand, the energy splitting due to spin-orbit coupling is very small near the $\Gamma$ point but much more pronounced near the K point. Consequently, stronger spin-orbit coupling favors the top valence band maxima to be at K point. Since varying the distance between bilayers can modulate the interlayer tunneling strength, without influencing the spin orbit coupling strongly (see Fig. 5), it is an ideal experimental parameter to engineer the $\Gamma$-K competition.

Because the highest energy valence bands are composed primarily of $d$-orbitals from the transition metal atoms, electronic correlations also influence the $\Gamma$-K energy competition. Given that the wave function is more extended for states near the $\Gamma$ point and more localized for states near the K point, correlations tend to push states near the K point to higher relative energy. This trend is observed in all of our DFT calculations. One may observe that in Molybdenum-based structures, the inclusion of electronic correlation produces a larger upward push, due to the 4$d$ orbital characteristic as opposed to the 5$d$ of the stable Tungsten-based structures.

With respect to the stacking order, we note that the AA stacked materials possess a much larger metal to metal distance when fully relaxed compared to their 2H counterparts. As a result, we expect that the AA bilayer splitting is smaller than in 2H, contributing to a lower $\Gamma$ valley. In Table 2, we see that all of the equilibrium AA stacked materials are K-valleys for LmBJ, due primarily to this weaker interlayer tunneling. As in the 2H cases, the correlation effect of the LmBJ potential is to push the K-valley up, but the magnitude of the increase is lesser in AA materials than in 2H.

Pressure engineering has been widely used as one of the ways to tune an interlayer distance in bilayer and multilayer TMDs Bhattacharyya and Singh (2012b); Fan et al. (2015); Dou et al. (2014); Su et al. (2014); Xia et al. (2021); Nayak et al. (2014, 2015). Decrease in the interlayer distance has been achieved experimentally by applying hydrostatic pressure through a DAC (Diamond anvil cell) Nayak et al. (2014, 2015); Dou et al. (2014); Xia et al. (2021). Pressure used to decrease the interlayer separation can be used to modify the band-gap, and in the semiconducting regime the change of interlayer distance estimated from optical probes can be up to 10% within 10 GPa. Nayak et al. (2015, 2014); Xia et al. (2021). While we do not attempt to make a quantitative study of the pressure-induced effects, we would like to see if our results are reasonably consistent with experiments. For this purpose, we employ the formula used in Ref. [Bhattacharyya and Singh, 2012c] that approximates pressure as

where $E_{tot}(h)$ is the total energy (from LDA+SO) at interlayer separation $h$, $h_{r}$ is the relaxed interlayer separation, and $A$ is the unit cell area in the $ab$ plane. For the heterobilayer MoSe₂/WSe₂ that has been experimentally studiedXia et al. (2021), we find that in order to achieve 3% change of the height, the required pressure is predicted to be $\sim 1$ GPa. Although the order of magnitude seems to be reasonably within the experimental range, more calculations must be done to have a quantitative description of the pressure effects and more experimental data will be needed for comparison. Applying pressure that decreases the layer separation can only enable transitions from K to $\Gamma$ - and not the reverse transition.

While applying pressure is an experimentally feasible way to decrease the interlayer distance, the strain engineering offers new possibilities to increase the interlayer distance. By applying tensile or compressive strain to the bilayer TMD, interlayer distance can be either increased or decreased. Transition between direct and indirect bandgaps has been reported in experimental and theoretical studies.Kumar and Ahluwalia (2013); Pak et al. (2017)

The most convenient way of engineering bilayer electronic structure is to apply electrical gate fields. Liu et al. (2012); Ramasubramaniam et al. (2011); Huang and Kaxiras (2016); Santos and Kaxiras (2013); Sharma et al. (2014); Zheng et al. (2016). We employ the ”zig-zag” approach in WIEN2k Schwarz and Blaha (2003b); Blaha et al. (2020) to apply a simulated gate voltage across the layers of $\sim 0.2$V, yielding transverse electric fields with magnitudes of $\sim 0.01-0.1$ V/A in the samples. Examples of the resulting band structure change for the experimentally relevant Cai et al. (2023); Mak et al. (2018); Zeng et al. (2023); Jin et al. (2021); Rivera et al. (2016) cases of homobilayer $\textrm{MoTe}_{2}$, and heterobilayer of $\textrm{MoSe}_{2}/\textrm{WSe}_{2}$ are shown in Fig. 7, and Figs. 8, 9, respectively. We note that the electric field produces the stronger effect near K than $\Gamma$ valleys. That the K valley is influenced more than the $\Gamma$ valley by the external field can be explained by the same model presented in the study of the bilayer grapheneMin et al. (2007). In essence, since the external electric field produces an electric potential difference between layers, its effect can be described by a two-band model capturing the layer degrees of freedom in Equation 5 of Ref. [Min et al., 2007]. Firstly, in the two-band model, the electric potential difference appears in the diagonal terms, whose effect in the eigenenergies will be reduced in general in the presence of the off-diagonal terms due to the interlayer tunneling. As a result, the shift of the eigenenergies is expected to be larger at momentum where tunneling is small. This explains our observation of the external electrical field having stronger effect near K than $\Gamma$ valleys. We note that the evolution of the $\Delta E_{K-\Gamma}$ is further complicated by the screening effect due to the Hartree potential that is shown to be large and well-captured by DFT calculations. In general, the screening effect decreases the electric potential difference actually seen in the bilayer systems in a non-monotonic way.

Another intriguing observation is that for homobilayers (Fig. 7) direction of the field does not matter as both layers maintain the same on-site potential in each layer and have a mirror symmetry between layers, and the field makes a positive contribution to $\Delta E_{K-\Gamma}$. For heterobilayers (Figs. 8, 9), we observe that changing the field direction produces opposite effects. If the field points from $\textrm{MoSe}_{2}$ to $\textrm{WSe}_{2}$, the top valence band (which is dominated by d orbitals of W atom) is pushed away from the Fermi energy and the gap is enlarged resulting in a reduction in $\Delta E_{K-\Gamma}$. If the field points from $\textrm{WSe}_{2}$ to $\textrm{MoSe}_{2}$, the top valence band is pushed closer to the Fermi energy and the gap is reduced. This can also been explained by the two-band modelMin et al. (2007). Because in the heterobilayer the on-site energy is already different between layers, a diagonal term in the two band model breaking the layer symmetry is already present at zero electric field. Therefore, the external field can drive K closer to or further away from the VBM in heterobilayers, depending on whether the direction of the electric field aligns with the chemical potential difference. Accordingly, the gate voltage is a powerful tuning parameter in TMD systems.