Atomistic theory of twist-angle dependent intralayer and interlayer exciton properties in twisted bilayer materials

Introduction
The stacking and twisting two (or more) layers of two-dimensional (2D) materials results in the emergence of a large-scale moiré pattern which gives rise to novel properties. Recently, such moiré materials have attracted tremendous interest as a platform to explore exotic optical properties1,2. For example, intralayer excitons (where the constituent electron and hole reside in the same layer) exhibit topological properties3 and a splitting of the lowest-lying bright exciton at small twist angles3,4,5,6. On the other hand, interlayer excitons (where the constituent electron and hole reside in different layers) show remarkably long lifetimes up to a few hundred of nanoseconds7,8,9,10,11, which makes them promising candidates for realizing exotic quantum phenomena, such as superfluidity12,13,14, supersolidity15 and Bose-Einstein condensation16, and for designing efficient excitonic devices7,17.
To gain a microscopic understanding of the properties of excitons in moiré materials and how these properties depend on the twist angle, strain, electric and magnetic fields and the composition of the individual layers, theoretical calculations can play an important role. To date, most theoretical studies of excitons in moiré materials are based on the effective mass approximation3,18. These studies have provided many important insights, but it is challenging to incorporate important factors, such as atomic relaxations or multi-valley effects. In principle, such effects can be straightforwardly captured using first-principles techniques, such as the ab initio density functional theory (DFT) combined with many-body perturbation theory (GW) and the Bethe–Salpeter equation (BSE) approach. However, application of this approach to twisted bilayer systems is numerically extremely challenging because of its unfavorable scaling with system size of this approach19. As a consequence, additional simplifications are usually introduced20. For example, computed monolayer spectra are interpolated to approximate the optical properties of the twisted bilayer21,22 or the single-particle wavefunctions of the moiré system are approximated in terms of monolayer wavefunctions23. However, a full solution of the BSE for a twisted bilayer system has not yet been achieved at small twist angles. The intriguing discovery of intralayer charge transfer excitons by Naik and coworkers23 emphasizes the importance of accurate atomistic models for capturing excitons in moiré materials, as these excitons cannot be captured using continuum models based on effective mass approximations.
In this paper, we introduce an atomistic approach to solve the BSE which exploits the localization of Wannier functions. Using this approach, we investigate the dependence of intralayer exciton properties on the twist angle in WS2/WSe2, a prototypical transition metal dichalcogenide heterobilayer. In agreement with experiments4,24,25, we find that the low-energy intralayer exciton peak splits into three peaks at small twist angles, and we observe good agreement with the available experimental data. We also analyze the wavefunction of the low-energy intralayer excitons and find that the character of the second lowest exciton depends sensitively on twist angle. We also observe interlayer excitons with an excitation energy of approximately 1.4 eV.
Bethe–Salpeter Equation in Wannier basis
To study properties of excitons in twisted bilayer WS2/WSe2, we solve the Bethe–Salpeter equation (BSE) for the interacting electron-hole Green’s function26,27. For zero-momentum excitons of relevance to optical absorption, the BSE Hamiltonian is given by
where ϵc(v)k denotes the quasiparticle energy of an electron in conduction (valence) band c (v) with crystal momentum k with wavefunction ψc(v)k. Also, D and X are the direct and exchange interaction terms27,28 given by
where x = (r, s) is a composite index representing both spatial and spin degrees of freedom, and (W({bf{r}},{{bf{r}}}^{{prime} })) and (V({bf{r}},{{bf{r}}}^{{prime} })) are the screened and bare Coulomb interactions, respectively. The evaluation of these integrals for large systems is computationally challenging, in particular for standard implementations that employ a plane-wave expansion of the single-particle wavefunctions and the Coulomb interactions29.
To overcome this problem, we expand the single-particle wavefunctions in a basis of Wannier functions30 and then exploit the localization of these basis functions to efficiently construct the BSE Hamiltonian. In particular, we first perform an ab initio DFT calculation for the twisted bilayer system and then construct Wannier functions using the Wannier90 code31. We include five d-like Wannier functions for metal atoms and three p-like Wannier functions for chalcogen atoms (note that spin-orbit coupling is included after the Wannier function generation, see Supplementary Information (SI) (See SI for the details of atomic structures, theoretical formulations, further technical details, benchmarking against existing studies, convergence tests, electronic wavefunctions of a 7 nm moiré unit cell, optical conductivity of isolated WSe2 and WS2/WSe2 heterobilayer, and a range of intralayer exciton wavefunctions), Section C for details). The DFT Bloch states ψmk can be expanded in terms of the Wannier function ϕn as
where R denotes a lattice vector, ({C}_{nm}^{{bf{k}}}) are the expansion coefficients, tn represents the position of ϕn in the home unit cell at R = 0 and N is the number of k points used to sample the moiré Brillouin zone. More details of the Wannier function generation can be found in the Methods section. We provide evidence for the completeness of our Wannier function basis by comparing the band structures obtained from DFT with the one obtained from diagonalizing the DFT Hamiltonian in the Wannier function basis, see Section C of the SI. The Kohn–Sham band gaps of monolayer WS2 and WSe2 are corrected using the GW method through rigid shifts32.
In the Wannier function basis, the BSE Hamiltonian can be approximated by using the localized and orthogonal nature of Wannier functions, along with the translational invariance of Coulomb interactions
A detailed derivation of this expression is provided in Section B of the SI along with a discussion of all approximations that have been made.
To determine the screened interaction W in a twisted bilayer system, we consider each layer as a polarizable sheet and solve the resulting electrostatic problem33,34,35,36,37,38, see Section C of the SI. When the electron and the hole reside on the same layer, the screened interaction is given by
where rs = 2π(α1 + α2) is a characteristic length scale determined by the polarizabilities αi of the two layers35, ϵbg is the background dielectric constant due to the presence of the substrate, and H0 and Y0 are the zeroth-order Struve and Bessel functions of the second kind, respectively.
When the electron and hole reside on different layers, the screened interaction is given by the same expression as above, but rs is replaced by rs + d with d being the average interlayer distance. We use ({alpha }_{{{rm{WS}}}_{2}}=6.03) Å ({alpha }_{{{rm{WSe}}}_{2}}=7.18) Å and d = 7 Å to represent the polarizabilities of the WS2 and WSe2 (calculated from the inverse dielectric tensor using DFT in conjunction with random-phase-approximation35) and the average interlayer separation of the heterobilayer, respectively.
The optical conductivity of the twisted bilayer WS2/WSe2 is given by
where ωS, and ({A}_{cv{bf{k}}}^{S}) denote the eigenvalues and eigenvectors of the BSE Hamiltonian, respectively, and the matrix elements of the momentum operator are expressed as (langle v{bf{k}}| {bf{p}}| c{bf{k}}rangle ={sum }_{{n}_{1}{n}_{2}}{({C}_{{n}_{2}v}^{{bf{k}}})}^{* }{C}_{{n}_{1}c}^{{bf{k}}}({nabla }_{{bf{k}}}{H}_{{n}_{1}{n}_{2}}^{{bf{k}}}+i({{bf{t}}}_{{n}_{1}}-{{bf{t}}}_{{n}_{2}}){H}_{{n}_{1}{n}_{2}}^{{bf{k}}})) with Hk being the Hamiltonian at k in the Wannier function basis. The effect of spin-orbit coupling is included as a perturbation39. More details can be found in the SI, Section C.
The exciton wavefunction for a fixed hole position ({{bf{r}}}_{h}^{0}) is obtained from
with ({n}_{3}in {{bf{r}}}_{h}^{0}) denoting the set of orbitals centered at ({{bf{r}}}_{h}^{0}).
To assess the accuracy of our approach, we apply it to monolayer WS2 and WSe2. We find good agreement with experimental results and plane-wave BSE calculations at a significantly reduced computational cost, see Section D of the SI. We also note that similar approaches have recently been used to compute optical properties of monolayers of MoS240,41. Below, we apply our approach to study low-energy intralayer and interlayer excitons in a WS2/WSe2 heterobilayer.
Results
Intralayer excitons
Figure 1 compares the computed optical conductivity from low-energy intralayer excitons in twisted WS2/WSe2 to the measured reflection contrast spectrum4 for a range of twist angles. Low-energy intralayer excitons are localized in the WSe2 layer and this allows us to consider only the relaxed WSe2 layer in the BSE calculation, see Methods for details. As the twist angle is reduced, the computed spectra exhibit several significant changes: (i) at twist angles below ~2∘, the spectrum exhibits three peaks (labeled I, II, and III) in the low-energy region, (ii) the spectral weight contained in the additional peaks II and III increases at small twist angles relative to the spectral weight in peak I, and (iii) the energy separation between peaks I and II and between peaks I and III shrinks as the twist angle is reduced, as shown in Fig. 2. These findings are in qualitative agreement with the experimental observations. In particular, the peak separations at small twist angles agree within a few meV. However, our approach underestimates the spectral weight contained in peak II.

a Optical conductivity (left panel) and wavefunctions of the two lowest-energy intralayer excitons (right panel) of twisted bilayer WS2/WSe2 at a range of twist angles. Gray arrows indicate the energies of the excitons and blue dots indicate fixed hole positions at the center of the home unit cell. The optical conductivities are normalized such that the height of the main peak is unity. Additionally, the main peak is shifted to match the main peak of the monolayer spectrum. b Experimental reflection contrast spectrum from ref. 4. When experimental data for specific twist angles near 0∘ is not available, we compare to the corresponding twist angles near 60∘ (indicated by AB) since intralayer excitons at 60∘ + θ have similar properties as those at θ4.

Energy separation between peaks I and II (black dots) and peaks I and III (black squares) in the optical conductivity (shown in Fig. 1) of twisted bilayer WSe2/WS2 as function of twist angle.
Figure 1a also shows the exciton wavefunctions of the two lowest-energy intralayer excitons for a fixed hole position at the center of the home unit cell. For all twist angles, the lowest-energy intralayer exciton has a hydrogenic 1s character. Near a twist angle of zero degree, however, the distribution of the electron around the hole becomes elliptical. In contrast, the character of the exciton which produces the second peak undergoes several qualitative changes as the twist angle is reduced. For the largest twist angle (e.g., 16∘), the second peak is produced by a 2p hydrogenic state. At intermediate twist angles (5.7∘ and 4.7∘), the peak is instead generated by an exciton with mixed 2s and 2p character. At intermediate twist angles (e.g., 2.2∘), the wavefunction associated with the second peak looks similar to the lowest-energy 1s state. We interpret this as a folded finite-momentum state. Finally, at very small twist angles (e.g., 0.4∘), the state producing the second peak again has a 2p-like character. The full evolution of the spectrum of all low-energy excitons is shown in the SI, Section H.
This evolution of the character of the intralayer exciton as function of twist angle is a consequence of the interplay between electron-hole interactions and the moiré potential: at large twist angles, the moiré potential is negligibly small and the exciton properties are determined by electron-hole interactions which only depend on the relative distance between electron and hole. At small twist angles, the moiré potential is large and determines the region where the exciton localizes.
In order to demonstrate the importance of the moiré potential, the electron distribution for different positions of the hole is examined at a small twist angle, see Fig. 3. When the hole is fixed in certain regions of the moiré unit cell (e.g. at the corner of the unit cell), the exciton wavefunction is small in magnitude for all values of the electron coordinate. In general, we find that the electron “follows” the hole (i.e. the exciton has a Wannier-type character), but the shape of the distribution changes. In contrast to Naik and coworkers23, we do not observe the exciton corresponding to peak III to have significant intralayer charge-transfer character. This is likely a consequence of the smaller unit cell used in our calculations which gives rise to a weaker moiré potential, see SI, Section F for a comparison of the Kohn–Sham wavefunctions, where we show that the electron and hole wavefunctions always have overlap at relevant stackings.

Wavefunctions of the low-energy intralayer excitons in 1.3∘ twisted WS2/WSe2 for different hole positions indicated by a blue dot. The home unit cell is marked with black solid lines.
To understand the origin of the twist-angle dependent exciton properties in more detail, we analyze the atomic structure of the twisted bilayer and its electronic band structure, see Fig. 4. We find that atomic relaxations play an important role at twist angles smaller than 2∘. For example, we find that the size of the AA stacking regions shrinks by 4% at a twist angle of 5.7∘ after atomic relaxations. However, at a twist angle of 1.3∘, the size of these regions shrinks by 50%, see Fig. 4a, b. See Methods for additional details on the estimation of the AA stacking area.

a, b Atomic relaxations lead to the shrinking of regions with unfavorable AA stacking and the growth of regions with low-energy BW/S and BSe/W stacking in WS2/WSe2. The high-symmetry stackings are labeled. c–f Electronic band structures of isolated WSe2 for several twist angles obtained from ab initio density functional theory calculations that include spin-orbit coupling. The WSe2 layer was isolated after performing structural relaxation of the WS2/WSe2 heterobilayer. The blue solid lines indicate the results of the relaxed structures, while the red dashed lines indicate the results of the flat monolayer WSe2 with the same dimensions.
Atomic relaxations have an important effect on the electronic band structure, see Fig. 4c–f. Here we compare the band structures of the WSe2 layer (after relaxation and removal of the WS2 layer) with those of a flat unrelaxed monolayer for several twist angles. At twist angles less than 2∘, the band structure of the relaxed system deviates significantly from that of the unrelaxed system. Importantly, relaxations lead to a flattening of the highest valence band and also of the lowest conduction band and a reduction of the band gap as the twist angle decreases42. In particular, the widths of the highest valence band and the lowest conduction band are reduced by ~30% at a twist angle of 1. 3∘ compared to the results for a flat unrelaxed monolayer.
Interlayer excitons
We also study interlayer excitons of twisted WS2/WSe2. For this, we solve the BSE calculations for the fully relaxed bilayer. Figure 5b shows the optical conductivity of the twisted bilayer at a twist angle of 5. 7∘. Near 1.45 eV, two peaks with very small intensities are observed, see inset. These peaks arise from interlayer excitons in which the electron is localized on the WS2 layer and the hole is localized on the WSe2 layer, see Fig. 5c, d. The spatial separation of the electron and the hole results in a very small transition dipole moment and oscillator strength for these excitons24,43,44. The calculated interlayer exciton energy is in agreement with previous experiments, with measured energies ranging from 1.35 to 1.45 eV4,10,24,45.

a Electronic band structure, b optical conductivity of the WS2/WSe2 twisted heterobilayer with a twist angle of 5.7∘. c, d Exciton envelope associated with the interlayer exciton (shown in the inset) with a hole fixed at the center of the home unit cell in the WSe2 layer. The conduction spin-split bands in (a) give rise to two peaks in (b).
The two peaks originate from transitions between spin-split partners of the CBM of WS2 and VBM of WSe2, see Fig. 5a. The magnitude of the spin splitting is approximately 30 meV. From the electronic band structure of the twisted bilayer, it can be observed that the bandgap is indirect with the valence band maximum located at Γ and the conduction band minimum at the K point of the moiré Brillouin zone. This suggests that interlayer excitons with finite momenta can have lower energies than zero-momentum excitons. We leave the exploration of the Bethe–Salpeter equation for finite-momentum excitons to future work.
Discussion
We have developed an efficient approach that exploits the localization of Wannier functions to solve the Bethe–Salpeter equation of twisted bilayer materials. We have applied this approach to study the properties of intra- and interlayer excitons in twisted WSe2/WS2. In agreement with experimental measurements, we find three peaks in the optical conductivity of intralayer excitons at small twist angles. Moreover, the energy separation between the peaks is well reproduced. The excitation energy of interlayer excitons is about 250 meV smaller than that of intralayer excitons, but their oscillator strengths are very weak. In the future, our approach can be straightforwardly applied to other multilayers of two-dimensional materials but also to other systems whose optical properties are difficult to describe with standard ab initio techniques, such as defects in solids or disordered materials.
Methods
Structure generation
All the twisted WS2/WSe2 heterobilayer structures were generated using the TWISTER package46 using lattice constants of 3.32 Å for monolayer WSe2 and 3.18 Å for WS2, see SI, Section A, for additional details.
Atomic relaxations
We performed atomic relaxations using accurate classical interatomic potentials fitted to density functional theory calculations. We used the Stillinger–Weber potential to represent interactions within individual layers47 and a Kolmogorov–Crespi potential for the interactions between the layers48. All relaxations are performed using the LAMMPS package (https://lammps.org/)49. We used the FIRE algorithm50 to relax the atoms within a fixed simulation box with a force tolerance of 10−6 eV/Å for any atom along any direction. To determine the size of the AA region in the relaxed and unrelaxed structures, we define the order parameter u as the minimum in-plane translation needed to transform any stacking configuration within the moiré pattern to AA stacking51. We used ∣u∣ < 0.5 Å as criterion to define the AA region.
Electronic structure calculations
The electronic structure calculations were performed on both unrelaxed and relaxed structures using the SIESTA package that uses localized atomic orbitals as the basis52. We used norm-conserving Troullier-Martins pseudopotentials53 and the local density approximation to describe exchange-correlation effects54. We used a single-ζ plus polarization basis for the expansion of wavefunctions. For monolayer WSe2, we also carried out calculations using a double-ζ plus polarization basis, but did not observe any significant differences in the resulting exciton properties, see SI, Section E for details. For all calculations we used Γ point sampling to obtain the charge density and a plane-wave energy cutoff of 100 Rydberg. We used an energy shift of 0.02 Ry. A large vacuum spacing of 20 Å was used in the out-of-plane direction.
Wannier function generation
We used the relevant d and p orbitals as starting guess for the generation of Wannier functions via the one-shot projection method30 as implemented in the Wannier90 code31. For example, we project the Kohn–Sham wavefunctions onto atom-centered d (({d}_{xy},{d}_{yz},{d}_{zx},{d}_{{x}^{2}-{y}^{2}},{d}_{{z}^{2}})) orbitals for every W atom, and p (px, py, and pz) orbitals for every S and Se atom and then orthogonalize. We have compared the resulting exciton properties to those obtained using maximally localized Wannier functions and found no significant differences, see SI, Section E. A disentanglement procedure was used30,55. We also compared our exciton properties to those obtained from Wannier functions generated from the plane-wave DFT code Quantum ESPRESSO package56,57. Again, we find good agreement between the two approaches for monolayer WS2 and WSe2, see SI, Section D.
Exciton calculations
We have developed a PyMEX programme, a Python package for Moiré EXcitons (The updated package will be freely available on GitHub (https://github.com/imaitygit/PyMEX)), to solve the Bethe–Salpeter Equation (BSE). The software package uses the mpi4py58, numpy59, scipy60, cython61, and h5py62 libraries. The Kohn–Sham band gaps of monolayer WS2 and WSe2 are corrected using the GW method through rigid shifts32. We note that by including only a few valence and conduction bands while constructing the BSE Hamiltonian, we do not reach absolute convergence of the lowest lying exciton for small twist angle moiré unit cell. However, we carefully check the relative convergence of the lowest few excitons with respect to the number of bands. See SI, Section E for more details of these calculations. Throughout the paper, the onsite interaction for the Keldysh potential is regularized as W0 = U ⋅ Wr=a where a is the pristine unit cell lattice constant of the WS2 or WSe2, and the parameter U is chosen to be 128. Furthermore, in all the results presented in the main text, we do not include bare Coulomb interactions. The impact of including the bare Coulomb interaction is discussed in the SI.
Low-energy intralayer excitons are localized in the WSe2 layer since it has a smaller band gap than the WS2 layer. Moreover, the low-energy valence and conduction bands of the WSe2 in the twisted bilayer are derived from the K and (K^{prime}) valleys of the flat monolayer. The wavefunctions associated with these valleys are predominantly composed of d-orbitals of W atoms which are only weakly affected by interlayer hybridization (in contrast to states derived from the flat monolayer Γ valley)11,23,63. To study intralayer excitons in the twisted heterobilayer, we remove the WS2 layer after the relaxation and carry out large-scale ab initio DFT and BSE calculations on the WSe2 layer with the same atomic structure as in a twisted bilayer system, see Sec. A of the SI for additional details. For the intralayer excitons at twist angles of 16.1∘ and 5.7∘, we use a 15 × 15 × 1 grid and a 9 × 9 × 1 grid, respectively. A 5 × 5 × 1 grid is used for a twist angle of 4. 7∘, and a 3 × 3 × 1 grid is used for all other twisted bilayer calculations when setting up the BSE Hamiltonian. Results of convergence tests are shown in Section E of the SI.
For the interlayer excitons, we have used a 3 × 3 × 1 k-grid, 12 valence, and 28 conduction bands to construct the BSE Hamiltonian for 5.7∘ twisted WS2/WSe2 heterobilayer. We have incorporated spin-orbit coupling effects perturbatively. We also compute the optical conductivity of isolated WSe2 and the WS2/WSe2 heterobilayer at a 5.7∘ twist angle using the same k-grid for comparison, see SI, Section G.
Responses