Accurate tight-binding models for the $\pi$ bands of bilayer graphene

We derive an ab initio $\pi$-band tight-binding model for $AB$ stacked bilayer graphene based on maximally localized Wannier wave functions centered on the carbon sites, finding that both intralayer and interlayer hopping is longer in range than assumed in commonly used phenomenological tight-binding models. Starting from this full tight-binding model, we derive two effective models that are intended to provide a convenient starting point for theories of $\pi$-band electronic properties by achieving accuracy over the full width of the $\pi$ bands, and especially at the Dirac points, in models with a relatively small number of hopping parameters. The simplified models are then compared with phenomenological Slonczewski-Weiss-McClure–type tight-binding models in an effort to clarify confusion that exists in the literature concerning tight-binding model parameter signs.

Figure 1

Definition of the unit cell for Bernal or $AB$ stacked bilayer graphene. The $A^{\prime }$ site of the top layer sits right above a $B$ site of the lower layer. The $A$ and $B^{\prime }$ sites do not have a corresponding vertical neighbor. The solid lines in the right panel represent the large-amplitude hopping processes in literature tight-binding models, whereas the dotted lines indicate weaker processes that are responsible for particle-hole symmetry breaking and trigonal distortions. In the main text we quantify the role they play in bilayer graphene electronic structure.

Figure 2

Band structure of $AB$ stacked bilayer graphene obtained through Wannier interpolation of first principles LDA results for $6\times 6$ and $30\times 30$ $k$-point sampling densities. The $6\times 6$ interpolation is accurate except near the $M$ point, where the maximum deviation is $\sim 0.2$ eV, and the $30\times 30$ interpolation is accurate throughout the Brillouin zone. The difference near the $M$ point between the $6\times 6$ and $30\times 30$ Wannier interpolation is larger than in the single-layer graphene case.

Figure 3

Intralayer hopping terms in Bernal bilayer graphene evaluated starting from $A$ to $A$ and $B$ sites in the bottom layer with $30\times 30$ $k$-point sampling density. The intralayer Hamiltonian matrix elements are extremely similar to those of isolated graphene layers obtained in Ref. [16] in spite of the interlayer coupling effect. The difference between isolated layer and bilayer intralayer hopping is not visible on the scale of this figure.

Figure 4

Neighbor classification relative to a central reference point labeled 0. Neighbors of a given site can be grouped into symmetry equivalent subsets with identical hopping parameters. The neighbor groups are ordered by their distance from the reference point. $m\text{th}$ nearest neighbors of the reference point are labeled in this figure by the integer $m$. Left panel: Neighbor classifications for intralayer $AB$ pairs and for interlayer $B{A}^{\prime }$ and $B{B}^{\prime }$ pairs. Right panel: Neighbor classification for interlayer $A{A}^{\prime }$, $A{B}^{\prime }$ pairs. The structure factors (see text) in $A{A}^{\prime }$ matrix elements are identical to those in $AB$ elements. The structure factors in $A{B}^{\prime }$ acquire a mirror reflection in the $y$ direction that leads to a complex conjugation in the matrix elements. The integer labels are derived by grouping neighbors according to in-plane separation so that some $A{A}^{\prime }$ and $A{B}^{\prime }$ neighbor pairs share labels with pairs in the left panel.

Figure 5

Upper panel: Hopping parameters connecting site $A$ of the bottom layer to the $A^{\prime }$ and $B{}^{\prime }$ sites in the top layer. The distances on the horizontal axis correspond to in-plane projections of displacements. Lower panel: Hopping parameters connecting site $B$ of the bottom layer to the top layer sites. Note that two different values occur for some values of $d$ ($m=5,10,15,17$) because of the occurrence of neighbor pairs which share the same projected displacements that are not symmetry equivalent. These cases are marked with an $*$ symbol in Fig. 4.

Figure 6

Evolution of the $\vec{k}·\vec{p}$ model coefficients $C_{\alpha \beta 1}$ and $C_{\alpha \beta 2}$ with neighbor number truncation. These coefficients are used in construction of the effective tight-binding model in Eqs. (20, 21, 22). The comparison of results obtained by $6\times 6$ and $30\times 30$ calculations points to the need for denser sampling to capture distortions of the local $p_z$ orbitals by interlayer coupling.

Figure 7

Comparison of electronic structure near the Dirac point for the full tight-binding (FTB) model, the simplified tight-binding model explained in the text (F1G0), and the effective continuum model (C) defined by Eqs. (14, 15, 16) with the second order term in $k$ neglected. Deviations increase as we move away from the Dirac point, but the relative error is of the order of a few % for $q\sim 1/a$ and F1G0 is superior to the continuum approximation. The dotted square represents the region over which the deviation from the Dirac point $\left|q\right|<1/a$.

Figure 8

Comparison of the full tight-binding bands (FTB) as presented in Sec. 2, and the simplified models presented in Sec. 3 consisting of the continuum model (C), the simplified tight-binding models with a single $f_1$ structure factor (F1G0), and with structure factors up to $f_2$ and $g_2$ (F2G2). We can observe that the simplest models C, F1G0 have substantial deviations away from the Dirac point, whereas F2G2 provides a good compromise between simplicity and accuracy.

Figure 9

Two dimensional color scale plot of the LUMO band dispersion over the full Brillouin zone (top left) and near the ${\mathbf{k}}_D=(4\pi /3a,0)$ Dirac point. The points at which the conduction and valence bands are degenerate are represented by black dots which lie at the vertices of a triangle. Right panel: Particle-hole symmetry breaking illustrated by plotting the sum of the LUMO (lowest conduction) and HOMO (highest valence) band energies [$E_c\left(\mathbf{k}\right)+E_v\left(\mathbf{k}\right)$]. The main contribution comes from the $A{A}^{\prime }$ and $B{B}^{\prime }$ hopping processes as explained in the main text.