Moiré band model and band gaps of graphene on hexagonal boron nitride

Nearly aligned graphene on hexagonal boron nitride (G/BN) can be accurately modeled by a Dirac Hamiltonian perturbed by smoothly varying moiré pattern pseudospin fields. Here, we present the moiré-band model of G/BN for arbitrary small twist angles under a framework that combines symmetry considerations with input from ab initio calculations. Our analysis of the band gaps at the primary and secondary Dirac points highlights the role of inversion symmetry breaking contributions of the moiré patterns, leading to primary Dirac point gaps when the moiré strains give rise to a finite average mass, and to secondary gaps when the moiré pseudospin components are mixed appropriately. The pseudomagnetic strain fields, which can reach values of up to $\sim 40$ T near symmetry points in the moiré cell stem almost entirely from virtual hopping and dominate over the contributions arising from bond length distortions due to the moiré strains.

Figure 1

The first harmonics reciprocal space wave vectors and the moiré Brillouin zone (mBZ) in a perfectly aligned G/BN (${\vec{\mathrm{G}}}_m, \theta =0$) and at a finite twist angle (${\vec{\mathrm{G}}}_m^{\prime }, \theta >0$). The relative size of the mBZ with respect to the graphene BZ is defined by the factor $\tilde{\varepsilon }=\sqrt{{\varepsilon }^2+{\theta }^2}$, and is rotated by $\phi$ under twist, such that ${\chi }_{\theta }=cos\left(\phi \right)\approx \varepsilon /\tilde{\varepsilon }$. The Dirac electrons in moiré band are coupled by interger multiples of these wave vectors which generally lead to band splitting at the mBZ edges. The primary valley is located at the mBZ center ($\mathrm{\Gamma }$), while $K$ and $K^{\prime }$ refer to the secondary valleys at which sDC might be found.

Figure 2

(Top) Plot of the relaxation coefficient $C_R$ which defines the strain magnitude as a function of twist angle $\theta$ and a schematic representation of different symmetric stacking arrangements. It is shown in Appendix pp1 that $C_R\propto {\tilde{\varepsilon }}^{-2}$. (Bottom) In-plane strain ${\vec{u}}_{\vec{r}}$ and interlayer spacing $h\left(\vec{r}\right)$ at two different twist angles. In-plane strain causes a reconfiguration in the positions of graphene unit cells, which prefer the BA stacking (red), where one carbon atom sits on top of boron atom, and one carbon is in the middle of BN's hexagon. In this case the strain field points away from the BA stacking point since $a_{\mathrm{G}}<a_{\mathrm{BN}}$. When G/BN layers are at zero twist angle, in-plane relaxation is maximized and allows graphene unit cells to move a distance of up to $\sim 0.085\AA$. The out-of-plane relaxation brings the interlayer distance around BA to a minimum, while it is maximized around AA (black), consistent with BA (AA) as the most (least) energetically favorable configuration. Similar plots are obtained for twisted G/BN ($\theta =2^{\circ }$) at the same energy scale. Under this twist, the maximum displacement attainable by a carbon site decreases to $\sim 0.04\AA$. The twist introduces an additional contribution $k_R(\widehat{z}\times \nabla \mathrm{\Phi }\left(\vec{r}\right))$ to the strain field, which leads to a nonzero curl.

Figure 3

Moiré pattern Hamiltonian parameters for rigid and relaxed G/BN plotted against twist angle $\theta$ using the BA stacking as the reference origin. For the relaxed cases, we distinguish the $\alpha$ case that includes only stacking registry modification effects from the $\beta$ case that also includes the strain induced electronic structure change in graphene. The top figures represent the inversion symmetric terms $u_i$, while the bottom figures represent the inversion asymmetric terms ${\tilde{u}}_i$. The relaxation gives rise to nonzero average mass $m_0$, which we plot in the right bottom panel with a black solid line with a value of $\sim 3.5$ meV near zero twist angle. The changes in the moiré pattern Hamiltonian parameters as a function of twist due to full relaxation are greatest for site potential ($u_0, {\tilde{u}}_0$) and in-plane pseudomagnetic fields ($u_1, {\tilde{u}}_1$).

Figure 4

(Top) (a)–(c) Moiré band diagrams for $\theta =0^{\circ }$ G/BN obtained for different relaxation schemes using the parameters listed in Appendix pp1. The k-path in the mBZ follows $\mathrm{\Gamma }\to K^{\prime }\to M\to \mathrm{\Gamma }\to K\to M\to \mathrm{\Gamma }$. In the rigid case, we expect to find a sDC on $K^{\prime }$ at the valence band when graphene is hole doped. When relaxation is accounted for including only stacking registry change, the $\alpha$ scheme described in the main text, we find an indirect sDC gap due to a change in the band shape near the $M$ point. In the more complete $\beta$ scheme where the strain induced band structure change in graphene is accounted for the pseudospin fields arising due to stacking registry modification is partially canceled and the direct gap at the sDC is restored. Bottom panel: (d)-(f) Moiré band diagram for rigid (red) and relaxed $\beta$ (black) G/BN calculated at various twist angles and the corresponding density of states (DOS). The interlayer coupling in G/BN gives rise to gapped sDC around the $K^{\prime }$ point in the moiré Brillouin zone. This finding is consistent with the experimental observation of a resistivity peak and the reversal of the Hall resistivity sign at the hole side reported in Ref. [29]. The relaxation effects introduces strongest changes in the band structure at small twist angles whereas they are closely similar to the rigid band structure when the twist angle is increased. (a) Rigid, (b) Relaxed $\alpha$, (c) Relaxed $\beta$, (d) $\theta =0^{\circ }$, (e) $\theta =1^{\circ }$, and (f) $\theta =2^{\circ }$.

Figure 5

The magnitudes of the band gaps for the primary Dirac cone ${\mathrm{\Delta }}_p$ (left) and secondary Dirac cone ${\mathrm{\Delta }}_s$ (right) for the valence band $K^{\prime }$ valley in the mBZ plotted against twist angle and in-plane strain $\eta C_R\left(\theta \right)$, where $\eta$ is a dimensionless multiplicative factor used to locally increase or reduce the relaxation. When $\eta =0$, only out-of-plane relaxation is present in the system and our model of relaxation corresponds to the white dashed line along the horizontal axis for $\eta =1$. As a general feature we find that the primary gap increases for larger strains with its value doubling when the strain is three times larger. However, the sDC gap decreases with increasing strains and a well defined sDC feature is only found within a certain window in twist angle. The shaded region on the right panel represents the parameter space where the direct sDC gap is obscured by the higher-energy band which dips below the sDC energy at the $M$ point when the twist angle is increased.

Figure 6

Schematic representations (a)–(c) of the band triplets at the secondary valleys in the mBZ used to analyze the magnitude of the sDC gaps of G/BN. We show three representative cases in the Hamiltonian parameter space for the behavior of the bands at $K^{\prime }$. The three arrow heads separated by a phase of $2\pi /3$ represent the energy levels that can be found projecting them on the $y$ axis following Eq. (24). The length of the arrows given by the magnitude of ${\mathrm{\Delta }}_{\zeta \kappa }$ given in Eq. (25) quantifies the interplay between the moiré patterns that split the bands at the sDC valleys with respect to the energy coordinate origin located at $\zeta \upsilon G/\sqrt{3}$. In the insets, we plot ${\mathrm{\Delta }}_{\zeta \kappa }$ and the phase ${\varphi }_{\zeta \kappa }$, which determines the size of the sDC gap, see Eq. (26). (a) In a system with only inversion symmetric couplings $u_i$, the band triplets separate into a singlet-doublet structure, with an energy difference of $\sqrt{3}\left|{\mathrm{\Delta }}_{\zeta \kappa }\right|/2$. A gap at the sDC appears if the energy level crossing of the doublet located closer to the charge neutrality is lifted. (b) In the presence of inversion asymmetric couplings ${\tilde{u}}_i$, we expect to find an energy band splitting of $\tilde{\mathrm{\Delta }}=\left|{\mathrm{\Delta }}_{\zeta \kappa }\right|/2$ that leads to a gapped sDC. (c) A general situation involves both nonzero $u_i$ and ${\tilde{u}}_i$ terms. This diagram represents the situation of G/BN valence sDC ($\zeta =-1$), which is found on $K^{\prime }$ ($\kappa =-1$). The analytical expression for the sDC gap given in Eq. (26) is most conveniently represented when the phase ${\varphi }_{\zeta \kappa }$ lies within $-\pi /3<\arg [-\zeta {\mathrm{\Delta }}_{\zeta \kappa }]\le \pi /3$, see Eq. (27). (a) $u_i\ne 0,{\tilde{u}}_i=0$, (b) $u_i=0,{\tilde{u}}_i\ne 0$, and (c) $u_i,{\tilde{u}}_i\ne 0$.

Figure 7

The spatial plots of the pseudospin components are obtained in rigid G/BN using the ab initio parameters [16]. We also illustrate here the possibility of using three different stacking points as our coordinate reference: AA (black), AB (blue), and BA (red). The threefold symmetry of the system is respected around these special points, so that a translation from one symmetry point to another is equivalent to a $2\pi /3$-rotation of the moiré pattern parameters, see Eq. (11), that does not change the moiré band. The figures also indicate that different degrees of inversion symmetry depending on the chosen coordinate reference. We found that in G/BN, inversion symmetric couplings are maximized (minimized) around BA(AA), see Table 2.

Figure 8

(a) Potential energy density profile $U_p\left(\vec{r}\right)$, which results from graphene interaction with the BN substrate, calculated at $z_0=3.4\AA$. The panel (b) represents the pseudomagnetic field profile in a relaxed G/BN at $\theta =0^{\circ }$ and (c) as a function of twist angle at special stacking points in rigid (dotted) and relaxed G/BN (solid). An increase in the pseudomagnetic field is expected at small twist angle, with fields of magnitude 40 T are attainable in the vicinity of AA and BA stacking points. The relatively long magnetic lengths compared with the moiré period preclude the formation of well defined Landau levels originated by pseudomagnetic fields but may lead to snake states for appropriate spatial distributions of positive and negative fields.