An extended plane wave framework for the electronic structure calculations of twisted bilayer material systems

The lack of periodicity poses a major difficulty for the theoretical studies of the electronic properties of twisted 2D material systems, where an abundance of fascinating correlation phenomena and physics are hosted [1, 2, 3, 4]. A common approach is approximating them by commensurate supercells, such that the periodic boundary condition is restored and the Bloch theorem can be applied. However, in many cases the size of an acceptable commensurate cell is overly large and generally out of the capability of most computational source. A good example is the magic angle ($\sim 1.05^{\circ}$) twisted bilayer graphene (TBG). Yet the smallest supercell to approximate this TBG system requires more than 10,000 atoms, which puts serious challenge for DFT calculations using PW basis [5]. Meanwhile, this approximation also prohibits a more systematic error analysis.

Alternatively one could switch to tight binding (TB) basis, which relieves the computational cost to some extent. We further note that recently there have emerged improved TB models and corresponding efficient numerical algorithms, which treat the incommensurate system from rigorous mathematical foundations [6, 7, 8]. Even though, one still needs to worry about parameters in the hamiltonian, constructing localized orbitals, and convergence issues. In many scenarios one would like to have similar footings for the PW basis, which is in nature almost free from previous issues, to better facilitate the theoretical studies of incommensurate material systems.

Our previous work has partially fulfilled that goal by introducing a general PW framework for the quantum eigenvalue problem of the incommensurate systems, where major theoretical aspects have been discussed [9]. However, the number of PWs needed to discretize the Hamiltonian is huge, due to the fact that it is in principle a higher dimensional problem. Effective schemes to reduce the dimensions of the incommensurate problems are therefore needed. With that in mind, layer-splitting methods have been developed to simulate the time evolution of quantum states in the incommensurate potentials, where the states are evolving sequentially by each composed periodic potentials rather than the full incommensurate one. This reduces the original problem to several low dimensional sub-problems [10].

On the other hand, when we were trying to understand the extended-localized transitions in the incommensurate systems, we observed that the low-lying eigenstates in the higher dimensional reciprocal space were distributed in a very directional manner [11]. This gave us hints of better cutoff scheme. Later on, we solidified this idea by rigorous numerical analysis and efficient algorithms [12]. By splitting the cutoffs into a traverse one and a quadratically increasing one, we showed in the numerical experiments that a much smaller PWs basis set as well as faster convergence can be achieved. Yet this scheme has not been applied to more practical calculations in the 2D material systems.

Besides that, previous work commonly neglect the interlayer dimension (denoted as z) for simplicity [9, 12, 10]. But in the TGB, the interlayer coupling is crucial to give rise to the flat bands, and consequently the correlated phenomena such as the superconductivity and Mott insulator [1, 2]. However when trying to put back this dimension under the full PW basis, one needs a very long $z$ cell to separate the spurious image interactions from the periodic boundary condition. As a result, this leads to very large number of PWs from $z$ dimension. And the problem is further amplified by the incommensurate nature since we already have a higher dimensional problem in the $x$, $y$ directions. A natural way to improve is to replace the PW by localized functions in $z$ direction, though many formula in the framework have to be adjusted accordingly.

Based on the above discussions, in this paper we extend our PW framework in the following aspects: (1) a tensor-producted basis with PWs in the incommensurate dimensions and localized basis in z direction, with which we reformulated the eigenvalue problem; (2) the application of our newly developed sampling and cutoff techniques in more practical systems; (3) a quasi-band structure formulated under the small twisted angles and weak interlayer coupling limits. In principle, with (1) and (2) we can treat the full electronic structure calculations of twisted 2D material systems with moderate computational cost affordable by most modern computers. (3) is also equivalent to the mini Brillouin zone (BZ) description used by most continuum model Hamiltonian calculations [5, 13]. And even though ergodicity is more fundamental, with (3) we can better organize our calculations as well as understand the results. Also it enables us to compare with those band structure results from commensurate cell approximations or model hamiltonians.

As for numerical examples, the electronic structures of the linear TBG system with twist angle $\sim 1.05^{\circ}$ are studied. We have reproduced the famous flat bands with key features in quantitative agreement with those from experiments and other theoretical models [4, 5, 14, 15]. Moreover, our calculations have much less computational cost compared to the commensurate cell approximations. Also it is more extendable compared to the traditional model hamiltonians and tight binding calculations, since in this work we utilize at most two parameters, and in principle it can be made fully ab initio. Finally, this framework can readily accommodate nonlinear terms like Hartree energy and exchange-correlation energy and will laid the foundations for more effective and accurate DFT calculations for the incommensurate 2D material systems, which will be addressed in our future work.

This paper is organized as follow: in Sec. 2, the incommensurate quantum eigenvalue problem is reformulated using the PW$\otimes$localized-function basis set. In Sec. 3, the practical application of new cutoff techniques are discussed. In Sec. 4, we discuss the quasi-band structure that can be formulated under the small twisted angles and weak interlayer coupling limits. In Sec. 5, the linear TGB system with magic twisted angle is taken as numerical examples, and we show and discuss the results on the electronic structures and density distributions. In the last section, conclusions are given.

The linear eigenvalue problem of a twisted bilayer system can be written as:

where the interlayer dimension z has been explicitly included. x compactly represents the inplane dimensions $x$ and $y$. By definition, $V_{1}(\textbf{x},z)$ and $V_{2}(\textbf{x},z)$ are periodic potentials in x: $V_{j}(\textbf{x}+n\textbf{R}_{j},z)=V_{j}(\textbf{x},z)$ with $\textbf{R}_{j}$ the lattice constants for layer $j=1,2$. And incommensurate constraint is given by:

Under the full PW basis, a common practice to model slab systems (layers, surfaces and interfaces) is adding a vacuum layer (usually about 20 Å) in $z$ direction, which serves to screen the spurious interactions between image parts. This will make the number of PWs a few orders of magnitude larger compared to bulk calculations and may cause some convergence issues [16]. These problems are further amplified in the incommensurate systems since one already has a large number of PWs in x. Here, we adopt a tensor-product basis set with PWs in x and some localized functions in $z$:

where $A$ is some large area within the x dimensions in which all coupled PWs are normalized. With properly chosen localized functions, such as atomic-like orbitals or finite elements, a much small basis set as well as a comparable accuracy with PW code can be expected.

Now we reformulate the eigenvalue problem under the new basis. For the PW part, elements in the set $\{\textbf{k}_{i}\}$ are still differed by the sum of two sets reciprocal lattice vectors $\{\textbf{G}_{1},\textbf{G}_{2}\}$ and $\{\textbf{Q}_{1},\textbf{Q}_{2}\}$, as derived in [9]:

where $n,m,k,l$ are integers and $i,j$ should be viewed as compact notation of these integers. For simplicity of the presentation, we assume that the localized functions are orthogonal (though it is not necessary):

And one can easily verify that $\{\Omega_{in}\}$ are also orthogonal. The discretization of the Hamiltonian in (1) follows similar procedures in [9], with some complexity comes in due to the localized functions. The matrix elements of kinetic operator take the form:

The matrix elements of each potential components can be calculated as:

where in (7):

In the above equations, $\textbf{k}_{ij}$ is a compact notation for $\textbf{k}_{i}-\textbf{k}_{j}$, $A_{c}$ the unit cell area, $V^{at}_{l}$ the atomic potential and $S_{l}$ the corresponding structural factor. The integration of $z$ can either be performed directly in the real space (6) or in the reciprocal space (7). The matrix elements of other ingredients in DFT, including the Hartree term, exchange-correlation term and nonlocal part of the potential can be derived similarly. Since we focus on the linear system in this paper, they will be reserved for our future studies.

With Eqs. (5), (6) and (7), one can in principle construct the full Hamiltonian matrix and solve for the eigenpairs. The above efforts relieve a large part yet unnecessary computational cost from the inefficient PW representation of the localized wavefunctions in $z$ direction. To further save the computational cost, we will briefly introduce our recently developed cutoff scheme and apply it to practical calculations. Some insights on the $k$-point sampling will be given as well, which helps to form the quasi-band structure description later on.

As termed and discussed in [9, 11], the coupled PWs $\{\textbf{k}_{i}\}$ in (4) are ergodic, meaning that they densely and uniformly distributed in 2d reciprocal space for infinitely large enough cutoffs. Though one could argue that in principle single point sampling with large cutoff is enough, in many practical cases the multi $k$-point sampling is still preferred for faster convergence. This is best manifested in systems that are very close to periodic systems, for instance the small twisted angle bilayer systems to be discussed. First we briefly recap some basis on $k$-point sampling and ergodicity.

In Fig. 1 we compare the distribution of PWs in bilayer hexagonal lattice systems with small ($1.05^{\circ}$) and large ($17^{\circ}$) twisted angles. As a variant of Eq. (4), the PW sets $\{\textbf{k}^{i}_{0}\}$ generated by $\textbf{k}_{0}$ are explicitly given by:

where different $\textbf{k}_{0}$’s generate same PW sets with infinite cutoffs due to ergodicity. Here, the index $i$ is in fact $i_{(n,m,k,l)}$, which depends on $n,m,k,l$. For simplicity, here and hereafter, we omit $(n,m,k,l)$ and denote it as $i$. In plotting the coupled PWs, we restrict the integers in Eq. (8) within a threshold $n_{t}$, namely $|n,m,k,l|<n_{t}$. And we adopt a brute cutoff scheme:

And generally we have $n_{t}\cdot|\textbf{G}|=2\sim 3~{}k_{c}$ to achieve satisfying uniformity.

Under small twisted angles, the PWs tend to form clustering sites with a quasi lattice arrangement (not strictly periodic) on the scale of the generating lattice vectors. The size of each cluster is proportional to the threshold $n_{t}$. In such case, multi-$k$ point sampling is preferred since single $k$ point would require a way too large $n_{t}$, which results in a extremely large number of PWs to cover the reciprocal space. As the angle increases, the spacing between the points increases and the clusters start to grow in size. With further increase, the sites contact and overlap, resulting in a rather uniform distribution in Fig. 1(b). In this case, single $k$ point sampling would suffice. The above observation also key to draw the quasi-band structure picture.

Yet the cutoff scheme can be improved. If we switch to the higher dimensional representation, previous scheme basically introduces a cutoff in some directions of the higher dimensional reciprocal space, while leaving others unaffected. This is the place we can make improvements. Based on our observations in [11] and more detailed analysis in [12], the energies of PWs increase quadratically along some directions as in the normal situation, but they also remain almost constant along some other directions. This is better illustrated in a 2-component 1d incommensurate system, as shown in Fig. 2.

Along the diagonal direction, the energies of PWs change quadratically, while along the anti-diagonal direction, they only fluctuate around certain value. Given this nature, one can replace the brute cutoff scheme (green dash) by our newly proposed one (black box), which also truncates from other directions, thus significantly reduce the number of PWs. Detailed analysis and numerics can be found in [12].

To be more specific for the twisted bilayer system, we define a conjugate $\tilde{\textbf{k}}^{i}_{0}$ to each $\textbf{k}^{i}_{0}$:

Note the sign change compared to Eq. (8). New cutoff scheme requires the PWs in $\{\textbf{k}^{i}_{0}\}$ and their traverse conjugate satisfy:

Generally $k_{L}>k_{W}$ due to the fact that they converge as $\textbf{e}^{-ck_{W}}$ and $1/k_{L}$, respectively [12]. In practical calculations we choose $k_{L}=3\sim 5~{}k_{W}$.

The combination of new basis and new cutoffs produces significant reduction in computational cost, making the eigenvalue calculations (either model calculations in this work or the full DFT calculations) within the capability of most modern computers. We illustrate this point by previous example. In Fig. 3, we plot the PW distribution with the new cutoff scheme, where $k_{L}=20$ and $k_{W}=6$.

As shown in the figure, new cutoff scheme has much less PWs. For the $1.05^{\circ}$ twisted angle, the number of PWs reduces from 42000 to 8000 in this case. To make further estimations, we assume the systems to be bilayer graphene. Together with the localized functions in $z$, the dimension of the Hamiltonian matrix falls in the range of $15000\sim 80000$ (since for Carbon atom, $2\sim 10$ localized functions would be suffice). The scale of this matrix is comparable with small to medium systems containing few tens of atoms in the full PW discretization, in drastic contrast with the 11164 atom supercell in the commensurate approximation scheme [5].

Before moving to numerics, we show that under the small twisted angles and weak interlayer coupling limits, a quasi-band structure picture can be formulated, though the ergodicity generally forbids such interpretations and is more fundamental. Here we give a rudimentary analysis under our framework. Though in the following we use TGB lattice as example, the conclusion should be general to other lattice systems.

The quasi-band structure picture largely roots in our previous observation that for small enough twisted angles and not too large $n_{t}$, the PWs in $\{\textbf{k}^{i}_{0}\}$ distribute in a clustering way, as shown in Figs. 1(a) and 3(a). If we zoom into each cluster, they further constitutes a hexagonal lattice with the lattice vectors given by:

Now we can choose a point in the central unit cell of the cluster (denoted as $\textbf{k}_{0}$) to generate the coupled PW set. If (1) the PW sets from different generating points are independent, and more importantly, (2) the eigenvalue problem can be well described on these PWs, then we are in a good position to recall the band structure description. Statement (1) is obviously true given a finite of cutoffs. To make statement (2) true, we need the assumption that the interlayer interaction is weak, which generally holds in the layered 2d materials.

Here we demonstrate it in a heuristic way. For the ease of presentation, we define following notations:

where each subset contains the PWs only coupled by one of the periodic potentials. By properly regrouping the terms in Eq. (8), $\{\textbf{k}^{i}_{0}\}$ can be partitioned using both subsets:

Other generating points, such as $\textbf{k}_{1}$, $\textbf{k}_{2}$, are the points in the same cluster of $\textbf{k}_{0}$, as shown in Fig. 4(b). In addition, the full basis is denoted by $\{\textbf{k}^{i}_{0}\otimes\xi_{u},\textbf{k}^{i}_{0}\otimes\xi_{d}\}$, in which $\{\xi_{u}\}$ represents the localized basis for the upper layer and $\{\xi_{d}\}$ the bottom layer.

First assume there is no interaction between the twisted bilayer. Under our basis set, the hamiltonian is block-diagonal in the sense that:

where blocks $H_{\{\textbf{k}^{G}_{n}\xi_{u}\}}$ and $H_{\{\textbf{k}^{Q}_{m}\xi_{d}\}}$ are the periodic Hamiltonians by Bloch theorem. For now, the upper half block only contains the eigenstates in the top layer and the lower half contains the the eigenstates from the bottom. And these eigenstates can be labeled by the PWs in the central cluster such as $\textbf{k}_{0},\textbf{k}_{1},\textbf{k}_{2}$, since there is no off-diagonal term to connect them.

By turning on the interlayer coupling, PWs in the same cluster start to mix. The full hamiltonian matrix now takes the following form:

In the next section, we will explicitly calculate $T_{\textbf{k}_{i},\textbf{k}_{j}}$. But here we refrain ourselves to qualitative argument. First, the terms in $T_{\textbf{k}_{i},\textbf{k}_{j}}$ are generally small, since they rely on the tunneling between layers, and are $2\sim 3$ orders of magnitude smaller compared to diagonal terms (or intralayer hopping) in the 2d materials. Furthermore, to connect nearby $\textbf{k}_{i}$ states in the same cluster, a process with a forward jump of $\textbf{G}_{i}$ to another cluster plus a backward jump of $-\textbf{Q}_{i}$ to the neighbor of the original site, is needed, which is $O(T^{2})$ in scale. One can further imagine the coupling between more distant PWs is even more weak. Therefore, the mixing may only extend to few neighboring sites, as shown in Fig. 4(b). And if we cares more about the electronic structures near the central region of the cluster, the new eigenstates can be still well described by the same PW set.

The above analysis underpines the justifiability of quasi-band structure picture, though such picture may only be useful at some specific regions of the original BZ, usually near the Fermi level with some states stand out due to extra symmetries and topological properties. For example, it is of more significance if $\textbf{K}_{m}$ coincides with the K of monolayer’s BZ, otherwise one only gets highly folded band states and a well separated band gap, which is hard to get any further insights. The non-trivial case, which often named as flat bands in the literature, will be discussed in detail next.

In this section, we take the bilayer graphene lattice systems as numerical example. For simplicity, we choose the atomic potential to be the one-body, local pseudopotentials. Even within this linear model, the crucial flat band features can be quantitatively reproduced ¹¹1We also note that in MacDonald’s original paper [13], in which the flat band model is first proposed, they used a even more simplified Dirac cone model. We believe linear models better manifest the effect of incommensurate geometries on the emergent non-trivial electronic structures. On the other hand, to study the related correlation phenomena and physics, one definitely has to go beyond single-particle picture.. Moreover the extension to DFT calculations is straightforward and we reserve the related technical issues, for instance the self-consistent schemes and error analysis, for our future work.

In this paper, we focus on the practical calculations of twisted 2d material systems and introduce extensions of our PW framework in the following aspects: (1) a tensor-producted basis set with PWs in the incommensurate dimensions and localized functions in $z$ direction, (2) the practical application of our newly developed cutoff techniques, and (3) a quasi-band structure picture under the small twisted angles and weak interlayer coupling limits. With (1) and (2), we have remarkably reduced the dimensions of hamiltonian matrix, which makes the electronic structure calculations of twisted bilayer 2D material systems affordable to most modern computers. And (3) helps us better organize the calculations as well as understand results. We further use the linear TGB system with magic twisted angles ($\sim 1.05^{\circ}$) as numerical examples. We have reproduced the famous flat bands with key features in good quantitative with other theoretical and experimental results. In terms of efficiency, our framework has much less computational cost compared to the commensurate cell approximations. While it is also more extendable compared to the traditional model hamiltonians and tight binding calculations. Lastly, nonlinear terms like Hartree energy and exchange-correlation energy can be readily included in the framework thus more effective and accurate DFT calculations of incommensurate 2D material systems can be expected in the near future.

The authors would like to thank the valuable discussions with Prof. Huajie Chen, Dr. Ting Wang and Prof. Daniel Massatt. This work was partially supported by National Key R & D Program of China under grants 2019YFA0709600 and 2019YFA0709601. Y. Zhou’s work was also partially supported by the National Natural Science Foundation of China under grant 12004047.