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

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 properties^1,2. For example, intralayer excitons (where the constituent electron and hole reside in the same layer) exhibit topological properties³ and a splitting of the lowest-lying bright exciton at small twist angles^3,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 nanoseconds^7,8,9,10,11, which makes them promising candidates for realizing exotic quantum phenomena, such as superfluidity^12,13,14, supersolidity¹⁵ and Bose-Einstein condensation¹⁶, and for designing efficient excitonic devices^7,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 approximation^3,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 approach¹⁹. As a consequence, additional simplifications are usually introduced²⁰. For example, computed monolayer spectra are interpolated to approximate the optical properties of the twisted bilayer^21,22 or the single-particle wavefunctions of the moiré system are approximated in terms of monolayer wavefunctions²³. 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 coworkers²³ 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 WS₂/WSe₂, a prototypical transition metal dichalcogenide heterobilayer. In agreement with experiments^4,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.

To study properties of excitons in twisted bilayer WS₂/WSe₂, we solve the Bethe–Salpeter equation (BSE) for the interacting electron-hole Green’s function^26,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 terms^27,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 interactions²⁹.

To overcome this problem, we expand the single-particle wavefunctions in a basis of Wannier functions³⁰ 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 code³¹. 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 WSe₂ and WS₂/WSe₂ 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, t_n 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 WS₂ and WSe₂ are corrected using the GW method through rigid shifts³².

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 problem^{33,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 r_s = 2π(α₁ + α₂) is a characteristic length scale determined by the polarizabilities α_i of the two layers³⁵, ϵ_bg is the background dielectric constant due to the presence of the substrate, and H₀ and Y₀ 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 r_s is replaced by r_s + 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 WS₂ and WSe₂ (calculated from the inverse dielectric tensor using DFT in conjunction with random-phase-approximation³⁵) and the average interlayer separation of the heterobilayer, respectively.

The optical conductivity of the twisted bilayer WS₂/WSe₂ 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 H^k being the Hamiltonian at k in the Wannier function basis. The effect of spin-orbit coupling is included as a perturbation³⁹. 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 WS₂ and WSe₂. 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 MoS₂^40,41. Below, we apply our approach to study low-energy intralayer and interlayer excitons in a WS₂/WSe₂ heterobilayer.

Intralayer excitons

Figure 1 compares the computed optical conductivity from low-energy intralayer excitons in twisted WS₂/WSe₂ to the measured reflection contrast spectrum⁴ for a range of twist angles. Low-energy intralayer excitons are localized in the WSe₂ layer and this allows us to consider only the relaxed WSe₂ 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 WS₂/WSe₂ 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. ⁴. 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 θ⁴.

Full size image

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 WSe₂/WS₂ as function of twist angle.

Full size image

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 coworkers²³, 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 WS₂/WSe₂ for different hole positions indicated by a blue dot. The home unit cell is marked with black solid lines.

Full size image

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 B^W/S and B^Se/W stacking in WS₂/WSe₂. The high-symmetry stackings are labeled. c–f Electronic band structures of isolated WSe₂ for several twist angles obtained from ab initio density functional theory calculations that include spin-orbit coupling. The WSe₂ layer was isolated after performing structural relaxation of the WS₂/WSe₂ heterobilayer. The blue solid lines indicate the results of the relaxed structures, while the red dashed lines indicate the results of the flat monolayer WSe₂ with the same dimensions.

Full size image

Atomic relaxations have an important effect on the electronic band structure, see Fig. 4c–f. Here we compare the band structures of the WSe₂ layer (after relaxation and removal of the WS₂ 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 decreases⁴². 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 WS₂/WSe₂. 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 WS₂ layer and the hole is localized on the WSe₂ 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 excitons^24,43,44. The calculated interlayer exciton energy is in agreement with previous experiments, with measured energies ranging from 1.35 to 1.45 eV^4,10,24,45.

a Electronic band structure, b optical conductivity of the WS₂/WSe₂ 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 WSe₂ layer. The conduction spin-split bands in (a) give rise to two peaks in (b).

Full size image

The two peaks originate from transitions between spin-split partners of the CBM of WS₂ and VBM of WSe₂, 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.

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 WSe₂/WS₂. 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.

Structure generation

All the twisted WS₂/WSe₂ heterobilayer structures were generated using the TWISTER package⁴⁶ using lattice constants of 3.32 Å for monolayer WSe₂ and 3.18 Å for WS₂, 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 layers⁴⁷ and a Kolmogorov–Crespi potential for the interactions between the layers⁴⁸. All relaxations are performed using the LAMMPS package (https://lammps.org/)⁴⁹. We used the FIRE algorithm⁵⁰ to relax the atoms within a fixed simulation box with a force tolerance of 10⁻⁶ 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 stacking⁵¹. 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 basis⁵². We used norm-conserving Troullier-Martins pseudopotentials⁵³ and the local density approximation to describe exchange-correlation effects⁵⁴. We used a single-ζ plus polarization basis for the expansion of wavefunctions. For monolayer WSe₂, 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 method³⁰ as implemented in the Wannier90 code³¹. 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 (p_x, p_y, and p_z) 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 used^30,55. We also compared our exciton properties to those obtained from Wannier functions generated from the plane-wave DFT code Quantum ESPRESSO package^56,57. Again, we find good agreement between the two approaches for monolayer WS₂ and WSe₂, 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 mpi4py⁵⁸, numpy⁵⁹, scipy⁶⁰, cython⁶¹, and h5py⁶² libraries. The Kohn–Sham band gaps of monolayer WS₂ and WSe₂ are corrected using the GW method through rigid shifts³². 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 W₀ = U ⋅ W_r=a where a is the pristine unit cell lattice constant of the WS₂ or WSe₂, and the parameter U is chosen to be 1²⁸. 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 WSe₂ layer since it has a smaller band gap than the WS₂ layer. Moreover, the low-energy valence and conduction bands of the WSe₂ 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 WS₂ layer after the relaxation and carry out large-scale ab initio DFT and BSE calculations on the WSe₂ 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 WS₂/WSe₂ heterobilayer. We have incorporated spin-orbit coupling effects perturbatively. We also compute the optical conductivity of isolated WSe₂ and the WS₂/WSe₂ heterobilayer at a 5.7^∘ twist angle using the same k-grid for comparison, see SI, Section G.