Ab initio theory of moiré superlattice bands in layered two-dimensional materials

When atomically thin two-dimensional (2D) materials are layered, they often form incommensurate noncrystalline structures that exhibit long-period moiré patterns when examined by scanning probes. In this paper, we present an approach that uses information obtained from ab initio calculations performed on short-period crystalline structures to derive effective Hamiltonians that are able to efficiently describe the influence of the moiré pattern superlattices on electronic properties. We apply our approach to the cases of graphene on graphene (G/G) and graphene on hexagonal boron nitride (G/BN), deriving explicit effective Hamiltonians that have the periodicity of the moiré pattern and can be used to calculate electronic properties of interest for arbitrary twist angles and lattice constants.

Figure 1

(Left) Schematic representation of two commensurate honeycomb layers with bottom layer sites indicated by light grey circles and top layer sites indicated by dark grey circles. The unit cell of the bilayer contains four sites, $A$ and $B$ for bottom layer and $A^{\prime }$ and $B^{\prime }$ for the top layer. The shaded region represents the primitive cell area $A_0$ used for the Fourier integrals described in the text. (Middle) The relative displacement between the honeycombs is specified by the displacement vector $\vec{d}$. We choose $\vec{d}=0$ for AA stacking in which the two honeycombs have no lateral displacement. For $\vec{d}=(0,a/\sqrt{3})$, we have the AB stacking where the top layer $A^{\prime }$ site is directly above the bottom layer $B$ site. The bilayer lattice is a periodic function of $\vec{d}$ and the primitive cell for this periodic dependence is shaded grey in the left figure. It is convenient to use the rectangular $a\times \sqrt{3}a$ area enclosed by a dotted line in the left figure to illustrate the dependence of the bilayer Bloch bands on $\vec{d}$. (Right) This panel specifically indicates the points within the rectangular area at which the high symmetry $AA$, $BA$, and $AB$ stacking arrangements occur and is helpful for the interpretation of later figures.

Figure 2

Energy and Dirac point gap landscape as a function of stacking. When the variations are smooth the whole map can be accurately interpolated from the values at the three symmetric stacking points at AA, AB, and BA using Eq. (21) expanded in the first shell of $\vec{G}$ vectors or first harmonic approximation. (Left) Total energy per unit cell relative to AA stacking as a function of displacement $\vec{d}$. These results are for G/G with constant in-plane lattice constant and vertical separation $c=3.35\text{Å}$. The highest energy configuration corresponds to AA stacking and the lowest to AB or BA stacking. Gap refers to the separation between conduction and valence bands at the Dirac point, which vanishes at AB and BA points. (Right) The same plots for G/BN with both sheets constrained to have the in-plane self-consistent LDA lattice constant of graphene and vertical separation $c=3.35\text{Å}$. The lowest-energy stacking configuration corresponds to BA stacking with one of the two carbon atoms in the graphene unit cell sitting on top of boron. The highest-energy stacking configuration corresponds to the AA arrangement in which the two carbon atoms in the unit cell sit on top of B and N. The AB configuration, which has C on top of N, has an intermediate energy. The scale of the dependence of total energy on $\vec{d}$ is similar in the graphene and hBN cases. Note that the Dirac point gap in the G/hBN case does not vanish at any value of $\vec{d}$, but that the typical gap scale is larger in the G/G case. This later property reflects stronger interlayer coupling.

Figure 3

Dirac-point $\pi$-band Wannier-representation Hamiltonian matrix elements as a function of sliding vector $\vec{d}$ for graphene/graphene. (Top) Left to right real and then imaginary parts of the $A{A}^{\prime }$ and then $A{B}^{\prime }$ interlayer matrix elements as a function of position $\vec{d}$ in the rectangular cell of Fig. 1. The $B{A}^{\prime }$ matrix element is closely related to the $A{B}^{\prime }$ matrix element as shown in Eq. (20) and the $B{B}^{\prime }$ matrix element is identical to the $A{A}^{\prime }$ matrix element. Interlayer coupling matrix elements have a typical magnitude $\sim$300 meV. We show later that the dependence of these four complex numbers on $\vec{d}$ is accurately described by a single real number. The color scales show energies in units of eV. (Bottom) Intralayer Wannier-representation Hamiltonian matrix elements. Left to right the real parts of the $AA$ and $BB$ matrix elements followed by the real then imaginary parts of the $AB$ matrix element. Typical matrix element values are $\sim$5 meV. For graphene on graphene, the spatial variation of intralayer matrix elements has a negligible influence on electronic properties. The color scale shows energy in units of meV.

Figure 4

Real part of the Fourier transform of $H_{A{A}^{\prime }}(\vec{K}:\vec{d})$ evaluated at $\vec{K}=(4\pi /3a,0)$. In Fourier space, interlayer coupling is strong only for three reciprocal lattice vectors, $\vec{G}=0$ and the two nonzero reciprocal lattice vectors for which $|\vec{K}+\vec{G}|=|\vec{K}|$. The imaginary part of $H_{A{A}^{\prime }}(\vec{K}:\vec{d})$ vanishes. At these values of $\vec{G}$, $H_{A{A}^{\prime }}(\vec{K}:\vec{G})$ is real with identical values of $0.113\pm 0.001$ eV.

Figure 5

(Left) Representation of the first shell of $\vec{G}$ reciprocal lattice vectors with their corresponding numeral labels used in the main text. The three circles in red correspond to the ${\vec{G}}_0,{\vec{G}}_{\pm }$ vectors with large interlayer tunneling coefficients. The three vectors ${\vec{K}}_j=\vec{K}+{\vec{G}}_j$ $j=0,\pm$ have the same magnitude. (Right) First shell moiré reciprocal lattice vectors $\tilde{G}=-\theta \widehat{z}\times \vec{G}$ for graphene on graphene result differ from the honeycomb reciprocal lattice vectors by a clockwise 90${}^{\circ }$ rotation and a reduction in size by a factor proportional to $\theta$. The solid black arrows represent the ${\vec{Q}}_j$ vectors that connect a $\vec{q}$ vector in the bottom layer to one in the top layer and the red hexagon encloses the moiré pattern Brillouin zone.

Figure 6

Real and imaginary parts of the Fourier transforms of the intralayer AA site diagonal Hamiltonian matrix element. Because of the symmetries of the honeycomb lattice the $\vec{G}$-dependent site potentials satisfy $H_{AA,\vec{G}}=H_{BB,\vec{G}}^{*}=H_{A^{\prime }{A}^{\prime },\vec{G}}^{*}=H_{B^{\prime }{B}^{\prime },\vec{G}}$. These contributions to the twisted layer Hamiltonian for G/G are small and can often be neglected.

Figure 7

Real and imaginary parts of the Fourier transform of the intralayer, intersubband Hamiltonian matrix element $H_{AB,\vec{G}}$. It follows from symmetry that $H_{AB,\vec{G}}=H_{A^{\prime }{B}^{\prime },\vec{G}}$. These matrix elements are also small.

Figure 8

Comparison of LDA G/G band structures (solid lines) and the one-parameter moiré band model (blue dashed lines), which retains only the first shell Fourier-expansion of interlayer coupling. Results are shown for commensurate G/G with AA, AB, and ${\vec{\tau }}_{\mathrm{br}}=(0,a/2\sqrt{3})$ bridge stacking. The electronic structure for BA stacking is identical to AB stacking. We find excellent agreement between the direct and approximate calculations, demonstrating the accuracy of the first shell approximation for interlayer coupling. The intralayer tight-binding Hamiltonian uses the models in Refs. [42, 43] with the experimental lattice constants of $a=2.46\text{Å}$, whereas the interlayer coupling is given by the first shell approximation as parametrized in Eq. (20). The small differences can be attributed mainly to the approximations involved in the first shell approximation for describing the interlayer coupling.

Figure 9

Moiré band structure and density of states of two graphene layers for four different relative orientation angles. Our results are similar to those obtained in Refs. [22, 23]. We plot the band structure as a function of momentum along the straight lines in $k$-space connecting points A, B, C, and A in Fig. 5. The accompanying density-of-states plots demonstrate the complex influence of interlayer coupling, which is responsible for many van Hove singularities.

Figure 10

Schematic Brillouin zones of graphene in black lines, of hBN in blue, and the moiré Brillouin zone (MBZ) in red for incommensurate G/BN. The difference between the lattice constants has been exaggerated to aid visualization. Both the size and orientation of the MBZ change continuously with twist angle due to the lattice constant mismatch. The three ${\vec{Q}}_j$ vectors given in Eq. (16) connect the $K$ points of graphene to those of hBN. The moiré pattern reciprocal lattice vectors can be constructed by summing pairs of ${\vec{Q}}_j$ vectors.

Figure 11

Displacement vector $\vec{d}$ dependence of Wannier representation interlayer Hamiltonian matrix elements for G/BN. (Top) Matrix elements $A{A}^{\prime }$, $A{B}^{\prime }$ connecting boron with carbon, and matrix elements $B{A}^{\prime }$, $B{B}^{\prime }$ connecting nitrogen with carbon. The interlayer coupling matrix elements vary over a large range $\sim 600$ meV. (Bottom) Displacement vector $\vec{d}$ dependence of intralayer Wannier representation Hamiltonian matrix elements for G/BN. On-site energies in the graphene layers vary by $\sim$60 meV. In these plots the site energies are plotted relative to their spatial averages $C_{0ii}$. [See Eq. (31).] The BN layer $AB$ intersublattice terms vary over a range of $\sim$35 meV, whereas the graphene layer $A^{\prime }{B}^{\prime }$ terms vary by $\sim$15 meV.

Figure 12

G/BN interlayer moiré band model parameters obtained by evaluating Fourier expansion coefficients for the layer separation dependence of Dirac point Wannier representation Hamiltonian matrix elements. As in the G/G case, three Fourier coefficients dominate interlayer coupling.

Figure 13

Fourier expansion of the Wannier-representation matrix elements $H_{AA}(\vec{K}:\vec{d})$, $H_{BB}(\vec{K}:\vec{d})$, $H_{A^{\prime }{A}^{\prime }}(\vec{K}:\vec{d})$, and $H_{B^{\prime }{B}^{\prime }}(\vec{K}:\vec{d})$ that describe the interlayer displacement dependence of site-energies in the hBN and graphene layers. These numerical results demonstrate that the site energies are accurately approximated by the model that includes only the first shell of reciprocal lattice vectors. The parameters of this model are listed in the main text. We have chosen the average energy on the carbon sites as the zero of energy. With this choice, the elements corresponding to $\vec{G}=0$ are $H_{AA}(\vec{K}:\vec{G}=0)=3.332$ eV for boron, $H_{BB}(\vec{K}:\vec{G}=0)=-1.493$ eV for nitrogen, $H_{A^{\prime }{A}^{\prime }}(\vec{K}:\vec{G}=0)=0$ eV, and $H_{B^{\prime }{B}^{\prime }}(\vec{K}:\vec{G}=0)=0$ eV.

Figure 14

Fourier expansion of the Wannier-representation matrix element $H_{AB}(\vec{K}:\vec{d})$, which describes intersublattice tunneling within the hBN layer, and $H_{A^{\prime }{B}^{\prime }}(\vec{K}:\vec{d})$, which describes intersublattice tunneling within the graphene layer. These numerical results demonstrate that the local coordination dependence of interlayer hopping processes is accurately approximated by the model that includes only the first shell of reciprocal lattice vectors. The parameters of this model are listed in the main text.

Figure 15

Relative displacement $\vec{d}$-dependence of the matrix elements of the two-band low-energy effective model for graphene on a hBN substrate. In the effective model map the $H_0=(H_{AA}+H_{BB})/2$ term represents a sublattice independent potential and $H_z=(H_{AA}-H_{BB})/2$ acts as a mass term in the Dirac equation. The off-diagonal matrix-elements $H_{AB}$ accounts for changes in the bonding pattern within the graphene layer. Here, $A,B$ refer to the sublattice sites of graphene.

Figure 16

Mass and pseudospin field terms in the effective Hamiltonian as a function of displacement $\vec{d}$. The $H_z$ term is due to sublattice potential difference and vanishes along lines of this two-dimensional plot. The simultaneous presence of finite $H_x$ and $H_y$ and $H_z$ terms in the effective Hamiltonian implies that the Dirac-point gap [Eq. (39)] does not vanish at any relative displacement.

Figure 17

(Top) Fourier expansion coefficients for the effective model matrix element $H_0(\vec{K}:\vec{d})$. This term captures the variation of the site-averaged potential across the moiré pattern. The first shell of reciprocal lattice vectors dominates. (Bottom) Fourier expansion coefficients for the mass term in the effective model, $H_z(\vec{K}:\vec{d})=[H_{AA}(\vec{K}:\vec{d})-H_{BB}(\vec{K}:\vec{d})]/2$.

Figure 18

Fourier expansion coefficients for the off diagonal effective model matrix element $H_{AB}(\vec{K}:\vec{d})$. The first shell of reciprocal lattice vectors dominates.

Figure 19

Comparison of the LDA band structure (solid black), the four-band moiré band model (dashed blue lines) and the low energy two-band model (dashed red lines) with the first shell used for the superlattice potentials. The commensurate G/BN arrangements plotted are AA, AB, BA, and a bridge stacking with ${\vec{\tau }}_{\mathrm{br}}=(0,a/2\sqrt{3})$. Note that he electronic structures of AB and BA stacking are different in the G/B case. For the intermediate bridge stacking we see a substantial reduction of the band gap due to a shifting in the Dirac cone momentum space location caused by in-plane pseudospin terms. The intralayer Hamiltonian of graphene is approximated using the massless Dirac model with the LDA Fermi velocity [42] while the boron nitride Hamiltonian is modeled with the same Dirac model with a mass term compatible with the LDA gaps. The interlayer coupling is given by the first shell approximation using the parametrizations of Eq. (20), with tunneling from boron and carbon sites distinguished as in Eq. (30) and the parameter set in Eq. (31). In these plots, the energy origin of the represented bands has been adjusted so that zero is in the middle of the band gap.

Figure 20

Modulation of the local potential fluctuations $H_0\left(\vec{r}\right)$ in real space for different twist angles. These plots illustrate the rotation of the moire pattern when $\left|\varepsilon \right|\sim \theta$, and the property that the Moire periodicity $L_M\sim a/\sqrt{{\varepsilon }^2+{\theta }^2}$ becomes shorter with increasing twist angle. The other pseudospin components of the local Hamiltonian illustrated Fig. 15 produce similar spatial superlattice patterns.

Figure 21

(Left) Schematic illustration of potential variations and local band edges as a function of $d_y$ for fixed $d_x=0$. The approximate conditions for local electron/hole charging discussed in the main text are satisfied over the segments identified by bold black lines. The difference between the Dirac point gap $\tilde{\Delta }$ and the absolute value of the mass $\left|H_z\right|$ reflects the influence of the in-plane pseudospin fields terms. (Right) Sliding vector $\vec{d}$ dependent map of electron and hole puddles.

Figure 22

Changes in total energy, Dirac point gaps, average layer separation distance, and Fermi energy resulting from allowing self-consistent LDA relaxation in the out of plane $z$ axis for G/G.

Figure 23

Changes in total energy, Dirac point gaps, average layer separation distance measured from the minimum separation distance, and Fermi energy resulting from allowing self-consistent LDA relaxation in the out of plane $z$ axis for G/BN.