Electronic spectrum of twisted bilayer graphene

We study the electronic properties of twisted bilayer graphene in the tight-binding approximation. The interlayer hopping amplitude is modeled by a function which depends not only on the distance between two carbon atoms, but also on the positions of neighboring atoms as well. Using the Lanczos algorithm for the numerical evaluation of eigenvalues of large sparse matrices, we calculate the bilayer single-electron spectrum for commensurate twist angles in the range $1^{\circ }\lesssim \theta \lesssim {30}^{\circ }$. We show that at certain angles $\theta$ greater than ${\theta }_c\approx 1.{89}^{\circ }$ the electronic spectrum acquires a finite gap, whose value could be as large as 80 meV. However, in an infinitely large and perfectly clean sample the gap as a function of $\theta$ behaves nonmonotonously, demonstrating exponentially large jumps for very small variations of $\theta$. This sensitivity to the angle makes it impossible to predict the gap value for a given sample, since in experiment $\theta$ is always known with certain error. To establish the connection with experiments, we demonstrate that for a system of finite size $\tilde{L}$ the gap becomes a smooth function of the twist angle. If the sample is infinite, but disorder is present, we expect that the electron mean-free path plays the same role as $\tilde{L}$. In the regime of small angles $\theta <{\theta }_c$, the system is a metal with a well-defined Fermi surface which is reduced to Fermi points for some values of $\theta$. The density of states in the metallic phase varies smoothly with $\theta$.

Figure 1

(a) Structure of the $AB$-stacked graphene bilayer. A twisted graphene bilayer is obtained by rotating the top layer by an angle $\theta$ around the axis connecting sites $A1$ and $B2$; quantity $t$ is the in-plane nearest-neighbor hopping, and ${\gamma }_1, {\gamma }_3$, and ${\gamma }_4$ are out-of-plane hopping amplitudes of the $AB$-stacked bilayer. These $\gamma \mathrm{s}$ are used to fix the fitting parameters of the function $t_{\perp }(\mathbf{r};{\mathbf{r}}^{\prime })$ (see the text). In this paper we use ${\gamma }_1=0.4\mathrm{eV}, {\gamma }_3=0.254\mathrm{eV}$, and ${\gamma }_4=0.051\mathrm{eV}$. (b) The large hexagons show the Brillouin zones of individual layers: the dashed hexagon (red) corresponds to the bottom layer, the dot-dashed hexagon (blue) corresponds to the top layer for the twist angle $\theta =21.{787}^{\circ }$ ($m_0=1, r=1$). The first Brillouin zone of the bilayer is shown by the central (green) thick solid hexagon. The next several Brillouin zones of the tBLG are depicted by the six surrounding (black) thin solid hexagons. The electronic spectra presented in Fig. 2 are calculated along the path specified by the triangle $\mathbf{\Gamma }{\mathbf{K}}_1{\mathbf{K}}_2$ (black). For the twisted bilayer, the Dirac point ${\mathbf{K}}^{\prime }$ (${\mathbf{K}}_{\theta }^{\prime }$) is equivalent to the point ${\mathbf{K}}_{\theta }$ ($\mathbf{K}$) if $r\ne 3n$. When $r=3n, {\mathbf{K}}_{\theta }\sim \mathbf{K}$ and ${\mathbf{K}}_{\theta }^{\prime }\sim {\mathbf{K}}^{\prime }$ (see the text). The tBLG Dirac points ${\mathbf{K}}_{1,2}$ are doubly degenerate: each of them is equivalent to one of two Dirac points of each graphene layer. For the particular case of the $(1,1)$ superstructure, ${\mathbf{K}}_1\sim \mathbf{K}\sim {\mathbf{K}}_{\theta }^{\prime }$ and ${\mathbf{K}}_2\sim {\mathbf{K}}^{\prime }\sim {\mathbf{K}}_{\theta }$.

Figure 2

(a)–(c) The spectra of the twisted bilayer graphene calculated for three different twist angles $\theta$ along the path $\mathbf{\Gamma }{\mathbf{K}}_1{\mathbf{K}}_2$ shown in Fig. 1. The spectrum shown in panel (a) corresponds to $\theta \approx 21.{787}^{\circ }$. It demonstrates a significant gap. The detailed behavior of this spectrum near the Dirac point is shown in panel (d) [the dispersion curves shown in panel (d) and (e) are calculated along the line parallel to vector $\mathrm{\Delta }\mathbf{K}={\mathbf{K}}_2-{\mathbf{K}}_1]$. Two almost-degenerate bands approaching the Fermi level $\mu$ from above and two almost-degenerate bands approaching it from below are clearly seen. For a much smaller angle $\theta \approx 4.{408}^{\circ }$, panel (b), the gap is much smaller, but still present; see panel (e). Panels (c) and (f) correspond to $\theta \approx 1.{89}^{\circ }$. The spectrum is gapless and three bands cross the Fermi energy forming the Fermi surface. The low-energy dispersion shown in panel (f) is calculated along the line passing through the Dirac point ${\mathbf{K}}_2$ perpendicular to the vector $\mathrm{\Delta }\mathbf{K}$ [the dot-dashed line in Fig. 8]. In panels (d)–(f) the Dirac point corresponds to $\mathbf{\delta }\mathbf{k}=0$.

Figure 3

Single-electron spectrum properties as functions of the twist angle $\theta$. In panel (a) the dependencies of the band gap $\mathrm{\Delta }$ (red solid curve) and band splitting ${\mathrm{\Delta }}_s$ (blue dashed curve) are shown for $r=1$ structures. In panel (b) the Fermi velocity $v_F$ for both $r=1$ structures (small blue dots), and $r\ne 1$ structures (larger red dots) is plotted. To extract the gap $\mathrm{\Delta }$, splitting ${\mathrm{\Delta }}_s$, and Fermi velocity $v_{\mathrm{F}}$, the numerically determined low-energy bands ${\overline{E}}_{\mathbf{k}}^{\left(\nu \right)}$ were fitted by Eq. (27).

Figure 4

The spectra of the twisted bilayer graphene calculated for structures $(9,4)$ (a) and $(2,1)$ (b). The spectrum for $(9,4)$ is calculated along the line connecting the points ${\mathbf{G}}_1(9,4)$ and ${\mathbf{G}}_2(9,4)$. The spectrum for $(2,1)$ is calculated in the folded (four times) reciprocal cell of the superlattice along the line connecting the points ${\mathbf{G}}_1(2,1)/4\approx {\mathbf{G}}_1(9,4)$ and ${\mathbf{G}}_2(2,1)/4\approx {\mathbf{G}}_2(9,4)$.

Figure 5

(a) Finite temperature ($T=0.01t$) density of states ${\rho }_T\left(E\right)$ calculated for the superstructure $(2,1)$ ($\theta \approx 13.{174}^{\circ }$), the superstructure $(3,2)$ ($\theta \approx 16.{426}^{\circ }$), and the superstructure $(7,3)$ ($\theta \approx 11.{635}^{\circ }$). The peaks in the range $0.25<|E/t|<0.5$ are Van Hove singularities due to the overlapping of two tBLG Dirac cones. The inset shows the zero-temperature density of states $\rho$ calculated for energies close to the Fermi level. The structures $(2,1)$ and $(3,2)$ are gapped; the gap for the structure $(3,2)$ is much lower than that for the $(2,1)$. The structure $(7,3)$ is gapless and has a finite density of states at the Fermi level. (b)–(c) Typical tBLG spectra close to the Dirac point for $r\ne 3n$ (b) and $r=3n$ (c) structures. The value of the band splitting parameter ${\mathrm{\Delta }}_s$ is shown by double arrows.

Figure 6

The band splitting ${\mathrm{\Delta }}_s$ as a function of the twist angle $\theta$ calculated for all superstructures with $N<2000$ (gray solid line). Different dashed curves connect the points corresponding to superstructures with the same value of $r$. We can see that, if $r$ is fixed, the band splitting decreases monotonously as $m_0$ grows (when $m_0$ grows, the angle $\theta$ decreases). However, when $r$ is not restricted, the splitting ${\mathrm{\Delta }}_s$ can change exponentially for weak variations of $\theta$.

Figure 7

Dependence of the band gap $\mathrm{\Delta }$ on the deviation $\delta \theta$ of the twist angle from the value $\theta (1,1)\cong 21.{787}^{\circ }$, calculated for finite sample of rhombic shape containing $151\times 151$ unit cells in each graphene layer. The total number of atoms in the sample $N_{\text{atoms}}=91204$. The inset shows the energy distribution of the first 200 electron levels close to the zero energy, calculated for $\delta \theta =0$.

Figure 8

Fermi surfaces of the superstructures $(17,1) [\theta \cong 1.{89}^{\circ }$, (a)] and $(18,1) [\theta \cong 1.{79}^{\circ }$, (b)] calculated at half filling. Different colors correspond to different bands intersecting the Fermi level $\mu$ at half filling. The first Brillouin zone (hexagon) and the reciprocal supercell (rhombus) are also shown. The dot-dashed line in (a) shows the way along which the spectrum presented in Fig. 2 is calculated.

Figure 9

Fermi surfaces of the superstructures $(17,1) [\theta \cong 1.{89}^{\circ }$, (a)], $(35,2) [\theta \cong 1.{84}^{\circ }$, (b)], and $(18,1) [\theta \cong 1.{79}^{\circ }$, (c)] calculated at half filling. The Fermi surfaces for the structures $(17,1)$ and $(18,1)$ are calculated by band folding of the original Fermi surfaces. Different colors correspond to different bands intersecting the Fermi level $\mu$ at half filling. Panels (a) and (c) show the same Fermi surfaces as Figs. 8 and 8, but in the folded Brillouin zone. (d) The low-energy density of states calculated for three superstructures with similar twist angles $\theta <{\theta }_c$. The density of states is calculated at finite temperature $T/t={10}^{-5}$ by numerical integration over the momentum; see Eq. (31). Note that in units of $t$ the energy window where the weight is enhanced is very narrow.