Construction of low-energy symmetric Hamiltonians and Hubbard parameters for twisted multilayer systems using ab initio input

A computationally efficient workflow for obtaining the low-energy symmetric tight-binding Hamiltonians for twisted multilayer systems is presented in this work. We apply this scheme to twisted bilayer graphene at the first magic angle. As the initial step, the full-energy tight-binding Hamiltonian is generated by the Slater-Koster model with parameters fitted to ab initio data at larger angles. Then, the low-energy symmetric four-band and 12-band Hamiltonians are constructed using the maximum-localization procedure subjected to crystal- and time-reversal-symmetry constraints. Finally, we compute extended Hubbard parameters for both models within the constrained random phase approximation for screening, which again respect the symmetries. Our workflow, exemplified in this work on twisted bilayer graphene, is straightforwardly transferable to other twisted multilayer materials.

Figure 1

The figure gives a workflow for creating TB Hamiltonians for TBG. The first line describes the process of creating the full energy range SK Hamiltonian using a $21.{79}^{\circ }$-TBG DFT calculation with subsequent Wannierization (see Sec. 2a) followed by an SK parameter fit. The resulting SK band structure for magic-angle TBG is further Wannierized (second line) to give the low-energy TB Hamiltonians (third line).

Figure 2

Frozen windows (between dashed horizontal lines) for Wannierization of the DFT band structures for single-layer and twisted bilayer graphene at $21.{79}^{\circ }$, respectively. The resulting Wannier-interpolated band structures (solid lines) and original DFT ones (dots) are also shown. The derived TB model for SLG was used in additional studies of Appendix pp4 and is placed here for a visual comparison.

Figure 3

Data for the out-of plane (left) and in-plane (right) hopping amplitudes for $21.{79}^{\circ }$-TBG (empty circles) computed with the ab initio TB scheme. Its SK parametric fit (red) together with the SK curve generated with commonly used parameters given in Sec. 2d [11, 15, 16, 17, 18] (orange). The result of a similar fit for $13.{17}^{\circ }$-TBG for $t_{\perp }$ is given by filled circles on the left-hand side. The good agreement between the SK fit at $21.{79}^{\circ }$ and the $13.{17}^{\circ }$ data supports our assumption that SK parameters are transferable between different angles.

Figure 4

The corrugation structure of the magic-angle TBG together with all Wickoff positions used in this work as centers of trial orbitals. One of the $6g$ sites is taken as (0.45,0.1) in lattice coordinates; the others are related by symmetry operations. The model of the out-of-plane corrugations seen on the plot is the one described by Eq. (20) [11, 34], which gives comparable sizes of AA and AB domains of larger and smaller interlayer distance, respectively, encoded by the color. All distances are in Å.

Figure 5

Wannier basis functions for four-band TB model. The color shows the weight of the $p_z$ component relative to its maximal value specified at the top right of each subplot. The given Wannier functions transform according to the atomiclike character shown at the bottom left of each plot.

Figure 6

The same as Fig. 5, but for 12-band model. The difference is that the centers of the site-symmetry irreducible representations and the center of mass of the orbitals do not generally match. In each plot, only the centers corresponding to the plotted orbital are shown.

Figure 7

Left: Averaged absolute error in a band energy per band and $\mathbf{k}$ point computed with a low-energy Hamiltonian as a function of nearest-neighbor cutoff $R_{\mathrm{cut}}$. Note that the 12-band model shows better convergence of flat bands (FB, green circles) than the four-band one. Right: Relative values of extended Hubbard parameters against the distance between centers of masses of orbitals. Fit to a soft Coulomb potential truncated by a Gaussian is given by the solid lines [Eq. (34)]. Fitted parameters are given as numbers in meV/Å units. The inset shows the relative scale of the on-site intra-orbital Hubbard matrix elements for each orbital with value for $s$ orbitals given by $U_0$. The exchange parameters correspond to the (ex.) label.

Figure 8

Comparison between the original band structure of the SK approach and its Wannier interpolation with the 12-band model.

Figure 9

Band structure computed with a 12-band low-energy Hamiltonian without introducing a crystal-symmetry constraint, in comparison with the original SK band structure.

Figure 10

The same as Fig. 6, but taken from the nonsymmetric Wannierization.

Figure 11

Top layer (left) and bottom layer (right) vertical coordinate of atoms in the vicinity of the long diagonal of the TBG UC in three different crystal structures. The model structure is given by Eq. (20), the structure relaxed using classical forces (lammps code) was given by the authors of Ref. [19], while the DFT-relaxed structure (vasp code) was taken from Ref. [35].

Figure 12

Electronic bands computed with various in-plane and out-of-plane TB model parametrizations for $1.{08}^{\circ }$-TBG. The first two rows are obtained with the SK parametrization for $t_{\perp }$ developed in this work (Sec. 2e), while the standard SK parameters for $t_{\perp }$ were used in the lower two rows. The standard SK parameters for $t_{\parallel }$ were used in the second and the fourth row, while the first and the third row correspond to using the ab initio TB parameters for SLG for an in-plane Hamiltonian of TBG with, however, reversed sign of the second-nearest-neighbor hopping amplitude. Column names are explained in the caption of Fig. 11, and all vertical energy axes are in eV. Finally, values for the standard SK parameters can be found in Fig. 3 of the main text and in references of its caption.

Figure 13

Valley- and layer-resolved band structure of $1.{79}^{\circ }$-TBG ignoring the out-of-plane coupling in the tight-binding model. The ${\mathcal{V}}_{\mathrm{tot}}$ is a difference between the valley characters ($\mathcal{V}$) of two layers. Note that odd and even bands are shifted vertically with respect to each other for a better observation of the valley character, which is normalized to span the $[-1,1]$ range.

Figure 14

The same as Fig. 13, but with finite-$t_{\perp }$ parametrization used in our work.