Moir\'e band model and band gaps of graphene on hexagonal boron nitride

TL;DR

This paper develops a moiré band model for graphene on hexagonal boron nitride, analyzing how symmetry and strain influence band gaps and pseudomagnetic fields, with implications for electronic properties.

Contribution

The paper introduces a comprehensive moiré band model for G/BN that incorporates symmetry and ab-initio data, elucidating the origins of primary and secondary Dirac point gaps.

Findings

Moiré patterns induce primary Dirac point gaps via inversion symmetry breaking.

Secondary gaps arise from specific pseudospin component mixing.

Pseudomagnetic fields up to 40 Tesla are mainly due to virtual hopping.

Abstract

Nearly aligned graphene on hexagonal boron nitride (G/BN) can be accurately modeled by a Dirac Hamiltonian perturbed by smoothly varying moir\'e pattern pseudospin fields. Here, we present the moir\'e-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\'e patterns, leading to primary Dirac point gaps when the moir\'e strains give rise to a finite average mass, and to secondary gaps when the moir\'e pseudospin components are mixed appropriately. The pseudomagnetic strain fields which can reach values of up to $\sim$40 Tesla near symmetry points in the moir\'e cell stem almost entirely from virtual hopping and dominate over the contributions arising from bond length distortions due to the moir\'e strains.

Tables

Table 1. Table 1: Two-way mapping between the parametrizations in Ref. Jung  et al. ( 2014 ) and Wallbank  et al. ( 2013 ) are explicitly presented here. These are equivalent to the more compact expressions in Eq. ( 8 )-( 10 ). In the left columns, we are expressing u i subscript 𝑢 𝑖 u_{i} and u ~ i subscript ~ 𝑢 𝑖 \tilde{u}_{i} in terms of the corresponding C μ subscript 𝐶 𝜇 C_{\mu} and ϕ μ subscript italic-ϕ 𝜇 \phi_{\mu} , and vice versa in the right columns. The approximate expression for small twist angles involving ε ~ ~ 𝜀 \tilde{\varepsilon} can be made exact by restoring ε / ε ~ ≃ χ θ ≡ cos ⁡ ( φ ) similar-to-or-equals 𝜀 ~ 𝜀 subscript 𝜒 𝜃 𝜑 \varepsilon/\widetilde{\varepsilon}\simeq\chi_{\theta}\equiv\cos(\varphi) , where φ 𝜑 \varphi is the rotation angle of the moiré pattern used in Eq. ( 28 ). The angle ϕ x ​ y subscript italic-ϕ 𝑥 𝑦 \phi_{xy} defined in this work is related to ϕ A ​ B subscript italic-ϕ 𝐴 𝐵 \phi_{AB} in Ref.  Jung  et al. ( 2014 ) through ϕ x ​ y = π / 6 − ϕ A ​ B subscript italic-ϕ 𝑥 𝑦 𝜋 6 subscript italic-ϕ 𝐴 𝐵 \phi_{xy}=\pi/6-\phi_{AB} for consistency of H A ​ B subscript 𝐻 𝐴 𝐵 H_{AB} term in Eq. (38) of Ref.  Jung  et al. ( 2014 ) with the H x ​ y subscript 𝐻 𝑥 𝑦 H_{xy} term defined in Eq. ( 6 ).

$u_0\upsilon G$	$C_0\cos \left({\varphi }_0\right)$	$C_0$	$\upsilon G\sqrt{u_0^2+{\tilde{u}}_0^2}$
${\tilde{u}}_0\upsilon G$	$-C_0\sin \left({\varphi }_0\right)$	${\varphi }_0$	$\arg \left[u_0-i{\tilde{u}}_0\right]$
$u_3\upsilon G$	$-C_z\sin \left({\varphi }_z\right)$	$C_z$	$\upsilon G\sqrt{u_3^2+{\tilde{u}}_3^2}$
${\tilde{u}}_3\upsilon G$	$C_z\cos \left({\varphi }_z\right)$	${\varphi }_z$	$\arg \left[{\tilde{u}}_3-i{u}_3\right]$
$u_1\upsilon G$	$-{\tilde{\epsilon }}^{-1}\epsilon C_{xy}\sin \left({\varphi }_{xy}\right)$	$C_{xy}$	$\left|{\epsilon }^{-1}\tilde{\epsilon }\right|\upsilon G\sqrt{u_1^2+{\tilde{u}}_1^2}$
${\tilde{u}}_1\upsilon G$	${\tilde{\epsilon }}^{-1}\epsilon C_{xy}\cos \left({\varphi }_{xy}\right)$	${\varphi }_{xy}$	$\arg \left[\epsilon \left({\tilde{u}}_1-i{u}_1\right)\right]$

Table 2. Table 2: Ab-initio moiré parameters in a perfectly aligned rigid G/BN ( z 0 subscript 𝑧 0 z_{0} = 3.35 Å ) presented in Refs.  Jung  et al. ( 2014 , 2015 ) . Each column represents different sets of parametrizations related to each other through 2 ​ π / 3 2 𝜋 3 2\pi/3 rotation as defined in Eqs. ( 11 - 12 ) which results in different local stacking at the moiré pattern center. The moiré Hamiltonian parameters C μ subscript 𝐶 𝜇 C_{\mu} , u i ​ υ ​ G subscript 𝑢 𝑖 𝜐 𝐺 u_{i}\upsilon G and u ~ i ​ υ ​ G subscript ~ 𝑢 𝑖 𝜐 𝐺 \tilde{u}_{i}\upsilon G are presented in meV units. We found that inversion symmetry is maximized (minimized) when BA (AA) stacking is chosen as the coordinate center.

Rigid G/BN	AA	AB	BA
$C_0$	10.13	10.13	10.13
${\varphi }_0$	$-{93.47}^{\circ }$	${26.53}^{\circ }$	${146.53}^{\circ }$
$C_z$	9.01	9.01	9.01
${\varphi }_z$	$-{171.57}^{\circ }$	$-{51.57}^{\circ }$	${68.43}^{\circ }$
$C_{xy}$	11.34	11.34	11.34
${\varphi }_{xy}$	${10.40}^{\circ }$	${130.40}^{\circ }$	$-{109.60}^{\circ }$
$u_0\upsilon G$	$-$0.613	9.06	$-$8.45
${\tilde{u}}_0\upsilon G$	10.11	$-$4.53	$-$5.58
$u_1\upsilon G$	2.05	8.64	$-$10.68
${\tilde{u}}_1\upsilon G$	$-$11.15	7.36	3.80
$u_3\upsilon G$	1.32	7.06	$-$8.38
${\tilde{u}}_3\upsilon G$	$-$8.91	5.60	3.31

Table 3. Table 3: Hamiltonian parameters based on the ab-initio calculations for rigid and relaxed graphene in Ref. Jung  et al. ( 2014 , 2015 ) , and the parameters obtained within the simplified relaxation scheme used in the present work. In a relaxed G/BN, there are two distinct effects which lead to the modifications in the moiré Hamiltonian parameters: (1) Change in the G/BN couplings due to local stacking modifications due to strains, (2) Strain-induced modification of electronic structure and pseudomagnetic fields in graphene. We denote with the label α 𝛼 \alpha the case in which only (1) is included, while β 𝛽 \beta refers to the case where both effects are present. All parameters are presented in meV except for the phases provided in degrees, and we label with “XY” the results obtained when only in-plane relaxation is allowed while we keep a constant interlayer separation distance at z 0 = 3.4 ​ Å subscript 𝑧 0 3.4 italic-Å z_{0}=3.4\,\,\AA . We have taken the BA stacking point as our coordinate reference that provides the parameter sets where the inversion asymmetric terms are smallest.

	$m_0$	$u_0\upsilon G$	$u_1\upsilon G$	$u_3\upsilon G$	${\tilde{u}}_0\upsilon G$	${\tilde{u}}_1\upsilon G$	${\tilde{u}}_3\upsilon G$	$C_0$	${\varphi }_0$	$C_{xy}$	${\varphi }_{xy}$	$C_z$	${\varphi }_z$
Rigid, Ref. Jung et al. (2014)   $z_0=3.35\text{Å}$	0	$-$8.46	$-$10.69	$-$8.38	$-$5.58	3.81	3.32	10.13	146.53∘	11.34	$-{109.60}^{\circ }$	9.01	68.43∘
Rigid, Ref. Jung et al. (2015)   $z_0=3.40\text{Å}$	0	$-$7.24	$-$9.15	$-$7.23	$-$4.85	3.26	2.79	8.71	146.18∘	9.71	$-$109.60∘	7.75	68.90∘
Relaxed, Ref. DaSilva et al. (2015)	3.74	$-$8.41	$-$6.68	$-$4.70	$-$3.40	3.05	3.13	9.07	157.99∘	7.34	$-$114.53∘	5.64	56.34∘
Relaxed $\alpha$	3.62	$-$9.48	$-$7.49	$-$5.94	$-$4.33	3.85	3.41	10.42	155.45∘	8.43	$-$117.20∘	6.85	60.14∘
Relaxed $\beta$	3.62	$-$4.93	$-$12.57	$-$5.94	$-$5.01	3.09	3.41	7.03	134.54∘	12.94	$-$103.81∘	6.85	60.14∘
Relaxed-XY $\alpha$	1.55	$-$7.73	$-$8.50	$-$6.72	$-$4.60	3.36	2.89	9.00	149.24∘	9.14	$-$111.57∘	7.32	66.73∘
Relaxed-XY $\beta$	1.55	$-$3.22	$-$13.59	$-$6.72	$-$5.28	2.60	2.89	6.18	121.38∘	13.84	$-$100.83∘	7.32	66.73∘

Equations

$$\vec{g}_{m}=\,\hat{R}_{2\pi(m-1)/6}\,\bigg{(}0,\frac{4\pi}{3a}\bigg{)},\hskip 5.0ptm\in\{1,2,...,6\}\\$$

$$\begin{array}{rl}{\vec{G}}_m& ={\vec{g}}_m^{BN}-{\hat{R}}_{\theta }{\vec{g}}_m=[(1+\epsilon )-{\hat{R}}_{\theta }]{\vec{g}}_m\\ & \approx \epsilon {\vec{g}}_m-\theta (\hat{z}\times {\vec{g}}_m).\end{array}$$

$$H=\hslash \upsilon \vec{p}\cdot \vec{\sigma }{\tau }_0+H_0(\vec{r}){\sigma }_0{\tau }_0+H_z(\vec{r}){\sigma }_3{\tau }_3+{\vec{H}}_{xy}(\vec{r})\cdot \vec{\sigma }{\tau }_3$$

$$V(\vec{r})$$

$$m(\vec{r})$$

$$\vec{A}(\vec{r})$$

$$(f_1(\vec{r}),f_2(\vec{r}))=\sum _m(1,i(-1{)}^{m-1})\exp (i{\vec{G}}_m\cdot \vec{r}),$$

$$\upsilon G(u_0-i{\tilde{u}}_0)=C_0{e}^{i{\varphi }_0},$$

$$\upsilon G({\tilde{u}}_3-i{u}_3)=C_z{e}^{i{\varphi }_z},$$

$$\upsilon G({\tilde{u}}_1-i{u}_1)=\cos (\phi )C_{xy}{e}^{i{\varphi }_{xy}}\simeq \frac{\epsilon }{\tilde{\epsilon }}{C}_{xy}{e}^{i{\varphi }_{xy}}.$$

$$\left(\begin{array}{ccc}{u}_0& u_1& u_3\\ -{\tilde{u}}_0& {\tilde{u}}_1& {\tilde{u}}_3\end{array}\right)$$

$${\varphi }_{\mu }$$

$$m_0=\frac{1}{A_M}{\int }_{A_M}d\vec{r}{\tilde{H}}_z({\vec{d}}_0(\vec{r})),$$

$$H_{\mu ,{\vec{G}}_j}=\frac{1}{A_M}{\int }_{A_M}d\vec{r}\exp (-i{\vec{G}}_j\cdot \vec{r}){\tilde{H}}_{\mu }({\vec{d}}_0(\vec{r})),$$

$$\Phi (\vec{r})=2{\tilde{\epsilon }}^{-2}{C}_R(\tilde{\epsilon })\Re e[e^{i{\varphi }_R(\tilde{\epsilon })}f(\vec{r})],$$

$$\vec{u}(\vec{r})=u_x\hat{x}+u_y\hat{y}=[1+k_R(\hat{z}\times 1)]\vec{\nabla }\Phi (\vec{r})$$

$$\vec{d}(\vec{r})={\vec{d}}_0(\vec{r})+\vec{u}(\vec{r})+h(\vec{r})\hat{z},$$

$$\delta H_0=\gamma {}^{\prime }(u_{xx}+u_{yy}),\gamma {}^{\prime }=\frac{a}{2}\frac{\partial {ϵ}_p^{\prime }}{\partial a},$$

$$\delta\vec{H}_{xy}=\gamma\big{[}(u_{xx}-u_{yy})\hat{x}-2u_{xy}\hat{y}\big{]},\,\hskip 7.0pt\gamma=\frac{3a}{4}\frac{\partial t}{\partial a}.$$

$$\upsilon G(\delta u_0-i\delta {\tilde{u}}_0)=-\gamma {}^{\prime }{C}_R(\tilde{\epsilon })g^2{e}^{i{\varphi }_R},$$

$$\upsilon G(\delta u_1+i\delta {\tilde{u}}_1)=\gamma C_R(\tilde{\epsilon })g^{2}F(\phi )e^{i{\varphi }_R},$$

$$F(\phi )=\cos (3\phi )+k_R\sin (3\phi ).\text{vspace}-3pt$$

$${\Delta }_p$$

$$E_{\zeta \kappa ,m}\approx \frac{\upsilon G\zeta -\Re e\left[\exp \left(i\frac{2\pi m}{3}\right){\Delta }_{\zeta \kappa }\right]}{\sqrt{3}},$$

$${\Delta }_{\zeta \kappa }$$

$$\displaystyle\Delta_{s,\zeta\kappa}=\big{|}\Delta_{\zeta\kappa}\sin{\phi_{\zeta\kappa}}\big{|}$$

$${\Delta }_{s,\zeta \kappa }$$

$$\cos (\phi )\equiv {\chi }_{\theta }=\frac{1+\epsilon -\cos (\theta )}{(1+\epsilon {)}^2-2(1+\epsilon )\cos (\theta )+1}\approx \frac{\epsilon }{\tilde{\epsilon }}.$$

$$E[\vec{u},u_{ij}]={\int }_{A_M}d\vec{r}[U_e(u_{ij}(\vec{r}))+U_p(\vec{u}(\vec{r}))],$$

$$U_{e}(u_{ij}(\vec{r}))=\frac{\lambda_{g}}{2}\big{[}\sum_{i}u_{ii}(\vec{r})\big{]}^{2}+\mu_{g}\sum_{ij}u_{ij}(\vec{r})^{2},$$

$$U_p(\vec{u}(\vec{r}))=\frac{1}{A_g}[U_0+U_1{\tilde{f}}_1({\vec{d}}_0+\vec{u}(\vec{r}))+U_2{\tilde{f}}_2({\vec{d}}_0+\vec{u}(\vec{r}))],$$

$$\begin{array}{rl}({\tilde{f}}_1(\vec{d}),{\tilde{f}}_2(\vec{d}))& =\sum _n(1,i(-1{)}^{n-1})\exp (i{\vec{g}}_n\cdot \vec{d}),\end{array}$$

Full text

I I. Introduction

Graphene is a single-atom thick sheet of carbon atoms arranged in a honeycomb lattice Novoselov et al. (2004, 2005); Geim and Novoselov (2007); Castro Neto et al. (2009); Das Sarma et al. (2011). In the past few years, hexagonal boron nitride (hBN) which also consists of van der Waals coupled bipartite honeycomb lattice layers, has emerged as a miracle substrate for graphene Dean et al. (2010); Geim and Grigorieva (2013). While graphene is a semimetal with a linear band crossing at the neutrality point, hBN is an insulator with a large bandgap of $\sim 5.8$ eV Zunger et al. (1976); Watanabe et al. (2004); Pacilé et al. (2008) due to its lack of inversion symmetry. Recent experiments have made it clear that graphene is very flat with reduced density of puddles when it is placed on hBN substrate Xue et al. (2011), allowing high carrier mobilities without sacrificing mechanical stability Dean et al. (2010). This drastic improvement in sample quality opened the door to the observation of new physics, including the discovery of new graphene fractional quantum Hall states Du et al. (2009); Bolotin et al. (2009), Fermi velocity renormalization Elias et al. (2011), and anomalously large magneto-drag Gorbachev et al. (2012). However, the influences of interlayer coupling between graphene and hBN become much stronger and are readily observed when both lattices have similar orientations. The interlayer coupling between carbon atoms in graphene and boron and nitrogen atoms in boron nitride is very large ($300-450{\rm meV}$) Jung et al. (2014), allowing for the possibility of strong substrate-induced distortions of the isolated graphene electronic structure when placed on a hBN substrate. Ab-initio theory has predicted that commensurate graphene on hBN (G/BN) would inherit a $50~{}{\rm meV}$ bandgap from the substrate Giovannetti et al. (2007).

Due to $\sim 1.7\%$ lattice mismatch between graphene and hBN lattices Pacilé et al. (2008); Jung et al. (2014); DaSilva et al. (2015), moiré supperlattices whose periodicity depends on the twist angle, were observed in the scanning tunneling microscopy Xue et al. (2011) and atomic force microscopy Woods et al. (2014). It has been shown that these results do not rely on the twist angle taking on discrete values giving commensurate superlattices, but hold for any continuous value Bistritzer and MacDonald (2011). Collectively, these results imply that G/BN should not have a band gap. However, the experimental observation of sizeable band gaps Amet et al. (2013); Hunt et al. (2013); Woods et al. (2014); Wang et al. (2016) has led to theories where the nonzero average mass generation introduced by the partial commensuration of the G/BN layers Jung et al. (2015); San-Jose et al. (2014), and electron-electron interaction effects Song et al. (2013); Bokdam et al. (2014) play a relevant role. Other manifestations of the moiré pattern effects in G/BN include the Hofstadter butterfly Dean et al. (2013); Ponomarenko et al. (2013), topological valley current Sui et al. (2015); Shimazaki et al. (2015), tunable Van Hove singularities in the low-energy regime Li et al. (2010), and the emergence of secondary Dirac cones (sDC) at the edge of the moiré Brilluoin Zone (mBZ) Ponomarenko et al. (2013); Wang et al. (2016).

In this work we present a theory of the moiré band model of G/BN for arbitrary twist angles under a framework that combines symmetry considerations and microscopic ab initio models for the moiré patterns Wallbank et al. (2013); Jung et al. (2014). Electronic structure theories of nearly aligned G/BN are most simply modeled through the continuum Hamiltonian of graphene subject to moiré patterns that vary slowly on an atomic scale. In this case we can formulate effective low energy theories in which the Hamiltonian has the periodicity of the moiré pattern and use the simplifications of Bloch’s theorem to obtain the moiré band Bistritzer and MacDonald (2011), thus bypassing the need to diagonalize large supercell approximants of incommensurable crystals often done for studying twisted bilayer graphene Lopes dos Santos et al. (2007); Shallcross et al. (2008, 2010); Moon and Koshino (2012). By establishing a unified framework that uses symmetry considerations Wallbank et al. (2013) and microscopic ab initio moiré band models Jung et al. (2014, 2015) we provide realistic estimates of the first harmonics parameters in G/BN and analytical expressions for the behavior of the band gap near the primary and secondary Dirac points. Our analysis allows us to understand the dependence of the moiré band Hamiltonian parameters on twist as well as the atomic lattice configuration through the study on rigid and relaxed structures, which in turn allows us to distinguish the strain fields resulting from virtual hopping and the effects of bond distortions due to relaxation strains. For the latter we further distinguish the changes in the moiré pattern due to modifications in stacking registry from the modifications in the intrinsic band structure of graphene due to the bond length distorting strains. Using approaches established in Refs. Jung et al. (2014, 2015) we examine how the twist and lattice relaxation can influence the band gap at the primary and secondary Dirac cones identifying the relevance of inversion asymmetric terms for opening up a band gap.

The paper is structured as follows. Section II briefly introduces and compares the moiré band presented in earlier literature. Subsequently, in section III we explain the relaxation model and the procedure to obtain the modified Hamiltonian parameters at various twist angles. In section IV we present the analysis on the resulting moiré band and the band gaps for the primary and secondary Dirac points, before we present the summary in section V.

II II. Moiré band Model

The Dirac electrons of monolayer graphene deposited on a BN substrate at nearly perfect alignment experience periodic moiré pattern perturbations whose length scale depends on the lattice constant difference between crystals and the twist angle. Using the BN substrate as a fixed reference the crystal lattice constant difference is represented through $\varepsilon=(a_{G}-a_{BN})/a_{BN}\approx-1.7$ and the twist angle is represented by $\theta$. A generic approach to analytically describe the moiré pattern perturbations in van der Waals crystals is based upon the realization that the interlayer coupling is smoothly varying over the moiré unit cell since the interlayer distance is substantially larger than the interatomic distanceBistritzer and MacDonald (2011). Consequently, for a heterostructure with a small lattice mismatch at a nearly perfect alignment ($|\varepsilon|,|\theta|\ll 1$), the influence of the substrate is effectively captured by the long wavelength components of the moiré pattern as defined by first harmonics $\vec{G}_{m}$ Wallbank et al. (2013); Jung et al. (2014, 2015) which are related to graphene and BN reciprocal lattice vectors $\vec{g}_{m}$ and $\vec{g}^{BN}_{m}$ through the following relations, see Fig. 1:

[TABLE]

[TABLE]

We denote the carbon-carbon distance on a graphene sheet as $a\approx 1.42\mbox{\AA}$, and $\hat{R}_{\theta}$ is a rotation by an angle $\theta$. In the G/BN systems the reciprocal lattice vector magnitude for the first harmonics is approximately $G=|\vec{G}_{m}|\approx(4\pi/3a)\sqrt{\varepsilon^{2}+\theta^{2}}$. They define the mBZ, which is $\lesssim 5\%$ of graphene-BZ’s size when $\theta\lesssim 2^{\circ}$ considered in our calculations, see Fig. 2. The corresponding moiré period $l_{M}=\sqrt{3}a/\sqrt{\varepsilon^{2}+\theta^{2}}$ ranges from $\sim$6 nm up to a maximum of $\sim$14 nm which is attained in a perfectly aligned G/BN. The moiré band Hamiltonian, constructed by adding the moiré pattern perturbation on top of the pristine graphene Hamiltonian,

[TABLE]

that is periodic over one moiré length, and we take the Fermi velocity $\upsilon\sim 10^{6}$ m/s Elias et al. (2011). We represent the sublattice and valley pseudospins through $\sigma$ and $\tau$ Pauli matrices and use the conventions in Refs. Wallbank et al. (2013); Tomadin et al. (2014) where $H$ acts on the four-component state $(\Psi_{AK},\,\Psi_{BK},\,\Psi_{BK^{\prime}},\,-\Psi_{AK^{\prime}})^{T}$. We have explicitly separated the $\vec{G}=(0,0)$ Fourier components $V_{0}$ for the potential and $m_{0}$ for the average mass or sublattice potential difference DaSilva et al. (2015); Jung et al. (2015). The constant $V_{0}$ represents a global shift in the Dirac cone’s energy that plays no role in the physical properties and is set to zero in our present analysis. The average mass $m_{0}$ plays an important role in opening a gap at the primary Dirac cone and it can be increased by commemsuration strains Jung et al. (2015). The spatially periodic G/BN couplings give rise to three distinct effects which can be represented by three local terms in the sublattice pseudospins basis: (1) $V(\vec{r})$ reflects the periodic sublattice potential, (2) $m(\vec{r})$ describes the local mass which opens a local gap at the neutrality point, and (3) $\vec{A}(\vec{r})$ can be interpreted as an in-plane gauge field that arises due to asymmetric hopping between one carbon atom and its neighbors on the opposite honeycomb sublattice Vozmediano et al. (2010). The Pauli matrices $\sigma_{i}$ and $\tau_{i}$ act on the sublattice and valley degrees of freedom respectively. The form of Eq. (II) reflects the time-reversal symmetry present in the G/BN system. (In the representation we use $\tau_{3}$ and $\vec{\sigma}$ are odd under time reversal.) In this work, we limit our analysis to the primary valley $\tau_{3}=1$.

There are two common ways to construct the moiré patterns Hamiltonians in the literature. In Ref. Wallbank et al. (2013) the moiré pattern Hamiltonians are parametrized using symmetry considerations that are distinguished by inversion-symmetric and inversion-asymmetric coupling coefficients, which are denoted as $u_{i}$ and $\tilde{u}_{i}$ respectively. The other representation is found in the ab-initio studies of G/BN couplings, which extract the moiré couplings from Fourier analysis and were represented through complex numbers in a magnitude-phase representation $C_{\mu}$ and $\phi_{\mu}$ Jung et al. (2014). Here we show that the two parametrizations can be mathematically related to each other exactly in the following manner

[TABLE]

where $\cos(\varphi)\simeq\varepsilon/\tilde{\varepsilon}$ and $\sin(\varphi)\simeq\theta/\tilde{\varepsilon}$ can be expressed in terms of $\tilde{\varepsilon}=\sqrt{\varepsilon^{2}+\theta^{2}}$ in the small angle approximation, and where $\varphi$ is the rotation angle of the moiré pattern discussed in Appendix A. The real periodic functions $f_{1}$ and $f_{2}$ are defined in the following way

[TABLE]

and the complex valued function $f=(f_{1}+if_{2})/2$ is similar to the first nearest neighbor structure factor of graphene’s tight binding model Jung and MacDonald (2013) and has been used in Ref. Jung et al. (2014, 2015) to represent the triangular modulation of the moiré pattern in real space. In the above equations the Hamiltonian parameters $u_{i},\,\tilde{u}_{i}$ are inversely proportional to $\upsilon G$ and decrease when the twist angle is increased, while the parameters in the magnitude-phase representation, $(C_{\mu},\varphi_{\mu})$, capture the local stacking dependent interlayer coupling and are insensitive to twist angle. We particularly note that for a complete mapping of the $\vec{H}_{xy}=\vec{A}_{xy}$ term in a twisted G/BN as defined in Ref. Jung et al. (2014) we need to include all of the four parameters $(u_{1},\tilde{u}_{1},u_{2},\tilde{u}_{2})$ that define the gauge fields in the $xy$-plane. However, the contributions from $u_{2}$ and $\tilde{u}_{2}$ can be absorbed into the global phase and can be neglected in our subsequent analysis of the moiré band Wallbank et al. (2013). This is equivalent to dropping the $\sin(\varphi)$ term in $\vec{A}(\vec{r})$ of Eq. (6). Their dependence on the twist angle can be seen clearly from our analytical mapping in Eqs. (6) and (10) where both terms are scaled according to $\chi_{\theta}=\cos(\varphi)\approx{\varepsilon}/{\tilde{\varepsilon}}$, which reflects the amount of rotation in the moiré first harmonic wave vectors under twist, see Fig. 2. The unification of both representations can be summed up by the following equations relating $\{u_{i},\tilde{u}_{i}\}$ and $\{C_{\mu},\phi_{\mu}\}$:

[TABLE]

[TABLE]

It becomes transparent now that $C_{\mu}$ quantifies the strength of the moiré pattern pseudospins and the phase $\phi_{\mu}$ determines the ratio between the inversion symmetric $u_{i}$ and asymmetric $\tilde{u}_{i}$ parameters. The explicit equations relating them are presented in Appendix A, together with the specific parameter values that we have used to model the band structures.

It should be noted that there are three different but equivalent parameter sets to describe the same bands which depend on the choice of stacking configuration at the origin $\vec{r}=0$ Wallbank et al. (2015). Such a possibility can be attributed to the presence of three centers on the AA, AB and BA stacking points around which $c_{3}$ symmetry is respected, see Appendix A for more discussions. The different sets of moiré patterns are related to each other through a rotation of $\pm 2\pi/3$ and is equivalent to a translation that changes the reference frame origin to a different local symmetry point $H(\vec{r})\rightarrow H(\vec{r}\pm\frac{4\pi}{3G^{2}}\vec{G}_{1})$

[TABLE]

while in the magnitude-phase representation the parameter rotations are achieved by shifting the phases by $\pm 2\pi/3$. Seemingly diffent solutions just represent changes in the reference stacking point that permute the mapping between boron, nitrogen and empty sites under inversion. This equivalence between parameter sets is illustrated further in Appendix A where we represent each moiré pseudospin term in real space for different stacking.

III III. Ab-initio Moiré patterns in rigid and relaxed G/BN

Initial attempts to model the moiré patterns of graphene on hexagonal boron nitride have used simplifying assumptions based on experimental observations or physical intuition to restrict the large parameter space Kindermann et al. (2012); Ortix et al. (2012); Yankowitz et al. (2012); Wallbank et al. (2013); Mucha-Kruczyński et al. (2016). In this work we use as reference the model Hamiltonian parameters defined through ab initio calculations that resolves this uncertainty that incorporate information for all possible stacking configurations of the local crystal Hamiltonian going beyond the two center approximation for the interatomic hopping Jung et al. (2014). In the presence of a small lattice mismatch $\varepsilon$ or twist angle $\theta$ the crystals are in general incommensurable but the local Hamiltonian $H(\vec{r})\equiv\widetilde{H}(\vec{d}(\vec{r}))$ at a given point $\vec{r}$ can be captured with short period commensurate geometry calculations containing few atoms. The electronic structure of incommensurable G/BN with rigid lattices leads to vanishingly small gaps Jung et al. (2014); Ortix et al. (2012); Ponomarenko et al. (2013); Sachs et al. (2011); DaSilva et al. (2015), in contrast to the $\sim 50$ meV single particle gaps when the lattice constants of G and BN are perfectly aligned Giovannetti et al. (2007). As we have shown in the previous section, the ab-initio model can be mapped naturally into the moiré Hamiltonian through Eq. (II). The parameters of the moiré Hamiltonian implied by these ab initio calculations are presented in Appendix A.

The global average of the mass term $H_{z}(\vec{r})$ in the moiré unit cell of area $A_{M}$ is given by

[TABLE]

and the first harmonics contributions to the Hamiltonian $H_{\mu,\vec{G}_{j}}$ that use the first shell of $\vec{G}$ vectors can be obtained through the corresponding Fourier transforms

[TABLE]

where the index $\mu\in\{0,x,y,z\}$ label the sublattice pseudospins, and the integrands are obtained from stacking dependent ab-initio calculations. Despite of the sublattice asymmetry introduced by the BN layer on the graphene sheet, the average mass term $m_{0}$ of rigid graphene remains approximately zero in the absence of strains Jung et al. (2014) and only a small gap is expected to open due to higher order perturbation terms Jung et al. (2015); Kindermann et al. (2012).

In the following we discuss effects of strains resulting from partial commensuration in a G/BN heterojunction. For the study about the effect of strains we follow closely work previously presented in Ref. Jung et al. (2015). When we allow the lattice structure to relax in a G/BN heterojunction the lattices undergo a partial commensuration expanding and compressing regions with different local stacking depending on their energy landscape. It was shown that such strains become relevant in the limit of long moiré lengths because the elastic energy resisting deformation decreases proportionally to $\varepsilon^{2}$, or equivalently with the inverse square of the moiré length. We use the Born-von Karman plate theory to capture the relaxation of the atoms and we use interlayer potential energies that depend on the local displacement between graphene and BN unit cells. This type of continuum approximation is justified for moiré lengths on the order of $\sim$10 nm that is between one to two orders of magnitude larger than the interatomic distance. The following Lamé parameters: $\lambda_{g}=3.25\text{ eV}\,\mbox{\AA}^{-2}$ and $\mu_{g}=9.57\text{ eV}\,\mbox{\AA}^{-2}$ for graphene and $\lambda_{BN}=3.5\,\,\text{eV}\mbox{\AA}^{-2}$ and $\mu_{BN}=7.8\,\,\text{eV}\mbox{\AA}^{-2}$ for BN characterize the elastic properties of the layers. In the analysis of strains presented here we focus on the restricted relaxation scheme where only the graphene sheet is allowed to deform while we keep the underlying rigid BN substrate rigid. A more general analysis shows that the hexagonal boron nitride sheet in contact with the graphene layer relaxes by an approximately equal magnitude but in the opposite sense when compared with the strains produced in the graphene sheet.

We obtain the strain vector fields by taking the gradient on a scalar function $\Phi(\tilde{\varepsilon};\vec{r})$ which respects the periodicity of the moiré pattern that we define using the magnitude-phase representation

[TABLE]

where the coefficient $C_{R}$ quantifies the degree of relaxation, and $\phi_{R}$ sets the direction of displacements around the symmetry points. The magnitude of $C_{R}$ depends on the strength of the van der Waals interaction and the elastic constants of the layer, and its magnitude decreases with growing twist angle. At an angle of $\theta\sim 2^{\circ}$, the coefficient $C_{R}$ is reduced by a factor of $\sim 6$ compared to the strain for zero twist angle, see Fig. 2. On the other hand, the phase $\phi_{R}$ exhibits little dependence on $\tilde{\varepsilon}$, and can be assumed to maintain a constant value for different twist angles. Further details for the relationship between $C_{R}$ and the adhesion potentials are presented in Appendix B.

The in-plane displacement vectors $\vec{u}(\vec{r})$ for a general twist angle can be readily obtained from the scalar function in the following manner

[TABLE]

where the factor k_{R}=\big{(}2+\lambda_{g}/\mu_{g}\big{)}(\theta/\varepsilon) grows with misalignment. The first term $\vec{\nabla}\Phi(\vec{r})$ contributes to translation of the local unit cell leading to local compressions and expansions, while the latter term $k_{R}(\hat{z}\times\vec{\nabla}\Phi(\vec{r}))$ is a non-zero curl term that distorts the carbon unit cell without changing its local area. Since ab-initio calculations on commensurate lattices show that G/BN prefers the BA stacking, the displacement vectors of graphene lattices should be pointing away from the position of BA stacking to partiatlly compensate for the lattice mismatch due to $a_{G}<a_{BN}$. The local displacement at certain positions can be as large as $\sim 5\%$ of the carbon-carbon distance, e. g. $\sim 0.07\,\AA$ when $\theta\sim 0$ while the magnitude of the local displacement $\vec{u}(\vec{r})$ decreases with twist, as we illustrate in Fig. 2. The cumulative dispacement within the moiré cell will be smaller than one lattice constant per moiré cell required for global commensuration.

To describe the out-of-plane stacking-dependent interlayer distance of G/BN we assume that the spatial profile of the out-of-plane relaxation is purely determined by $\vec{d}(\vec{r})$, and which follows the moiré periodicity, instead of explicitly solving the out-of-plan height variations based on the interlayer coupling potentials. Typically the height variations range between $3.28\,\,\mbox{\AA}$ and $3.49\,\,\mbox{\AA}$, see Fig. 2, with an average interlayer distance of $z_{0}=3.4\,\,\mbox{\AA}$ within the local density approximation (LDA) Jung et al. (2014, 2015). Fig. 2 shows that the preference over BA configuration in G/BN layers is also reflected in the smaller interlayer separation around BA stacking points.

The local relaxation leads to significant changes to the effective G/BN couplings experienced by the Dirac electrons. We capture this effect by modifying the local sliding vector $\vec{d}$ between the graphene and BN sites

[TABLE]

where $h(\vec{r})=\widetilde{h}(\vec{d})$ is the local interlayer distance for each position. The effects of strains alter the Hamiltonian in such a way that $\widetilde{H}(\vec{d}_{0})\rightarrow\widetilde{H}(\vec{d})$, which in turn modifies the spatial averages of the different pseudospin terms in Eq. (13) and (14) that capture the moiré pattern effects. We discuss in Appendix C the analytical expressions for $m_{0}$, $u_{i}$ and $\tilde{u}_{i}$ in the limit of $|\vec{u}|\ll a$.

The changes in the moiré patterns captured by Eq. (14) result from the local modification in the effective couplings of graphene to the BN substrate due to the reconfiguration of the relative positions between graphene and BN sites. Additionally, we have to bear in mind the modification of the electronic structure in graphene due to bond distortions that can be related to variations in the local site energy $\varepsilon^{\prime}_{p}$ and inter-site hopping parameter $t$ Vozmediano et al. (2010); Suzuura and Ando (2002). The pseudospin terms in the Hamiltonian $H_{0}$ and $\vec{H}_{xy}$ are consequently changed due to these additional contributions by

[TABLE]

[TABLE]

where $\gamma\,^{\prime}\approx 4.0\text{ eV}$ quantifies the rate of change of the site potential, and $\gamma\approx-4.5\,\text{ eV}$ captures the gauge fields produced by hopping asymmetry in the presence of bond distortions Wallbank et al. (2015); Mucha-Kruczyński and Fal’ko (2012); Ferone et al. (2011).

It is possible to gain further insight on how the strains modify the Hamiltonian assuming that they depend linearly on the symmetrized strain tensor u_{ij}\equiv\frac{1}{2}\big{(}{\partial u_{i}}/{\partial x_{j}}+{\partial u_{j}}/{\partial x_{i}}\big{)}+{\partial h}/{\partial x_{i}}\,{\partial h}/{\partial x_{j}}, where the potential and pseudomagnetic field terms respectively arise from the dependence of the local site energy $\varepsilon^{\prime}_{p}$ and the hopping parameters $t$ to the neighboring carbon atoms. With an in-plane strain profile which respects the lattice symmetry and the moiré periodicity, we are able to map the strain-induced potential and gauge field into Eqs. (4-6). The resulting changes in the first harmonics parameters are given by

[TABLE]

[TABLE]

where $g={4\pi}/{3a}$ is the magnitude of the reciprocal lattice vector of graphene, while $F(\varphi)$ is a dimensionless function that respects the three-fold rotational symmetry,

[TABLE]

The presence of $k_{R}$-dependent term only in Eq. (22) indicates that the non-zero curl contributions modify the gauge field terms but do not influence the electrostatic potential terms. As a consequence the site potential $\delta H_{0}$ contribution decreases together with $C_{R}(\tilde{\varepsilon})$ whereas the gauge field contribution $\delta\vec{H}_{xy}$ has a more complex behavior since $F(\varphi)$ changes sign around $\theta_{c}\approx 1.3^{\circ}$ as observed in the variation of the moiré parameters $u_{1}$ and $\tilde{u}_{1}$ with respect to twist angle, see Fig. 3. This critical angle is related with the elastic constants of graphene.

It is found that the intrinsic changes in the Hamiltonian due to out-of-plane strains are negligible with respect to the in-plane components, since the variation in the interlayer distance $\Delta h$ are orders of magnitude smaller than the moiré length. We can thus neglect the out-of-plane effects in its contribution to Eqs. (18) and (19).

Our analysis for $\delta H_{0}$ and the calculations of $\gamma\,^{\prime}$ have uncertainties related with the effect of electrostatic screenings by the carriers Gibertini et al. (2009), whereas the pseudomagnetic field effects are unaffected by electrostatic effects. We found that in G/BN structures with long moiré periods, the pseudomagnetic field realistically reaches 40 T near the AA and BA stacking points. These originate almost entirely from virtual strain fields produced by hopping of electrons to and fro from graphene to B and N sites Jung et al. (2014, 2015) while the contributions due to real strains have maximum values of $\sim$5 Te which amounts to  15% of the total pseudomagnetic field at selected points. The specific Hamiltonian parameter values for $\theta=0^{\circ}$ and how they are influenced by the lattice relaxation are discussed further in Appendix C. We present in Fig. 3 the behavior of the the Hamiltonian coefficients as a function of twist angle where we label with $\alpha$ the Hamiltonian parameters that incorporates only the modifications in the shape of the moiré pattern due to change in the local stacking, while the $\beta$ solutions include also the modifications of graphene’s Hamiltonian itself due to bond distortions resulting from the relaxation strains. When compared with the rigid lattice, we find that relaxations tend to enhance ($u_{0},\,\tilde{u}_{1}$ and $\tilde{u}_{3}$), and reduce ($\tilde{u}_{0},\,u_{1}$ and $u_{3}$) in the limit of small twist angles. The relaxation also allows the generation of a finite average mass term $m_{0}$ leading to a band gap on the order of $\Delta_{p}\approx 2|m_{0}|\sim 7$ meV within a graphene relaxation only scheme. The somewhat larger gaps than previous values in Ref. Jung et al. (2015) can be attributed to the simpler form of the strains assumed to calculate the elastic energy functional in the present calculations. Further enhancement of the band gap results when electron-electron interaction effects are included, making it possible to obtain magnitudes for the band gaps that are comparable to experiments Jung et al. (2015); Hunt et al. (2013). The additional effects of bond distortions in the Hamiltonian parameters included in the $\beta$-relaxation case through the first harmonic terms in Eqs. (18) and (19) do not affect the average mass term. The effects of $\beta$ relaxations in Fig. 3 tend to partially restore the parameters of the rigid Hamiltonian and this explains in part why the rigid model Jung et al. (2014, 2015) already gives a relatively accurate description of experimentally observed moiré band features Leconte et al. (2016); Yu et al. (2014). The modifications of the Hamiltonian parameters show most clearly near the mBZ edges and they introduce subtle changes in the behavior of the sDC features. The strain-induced gauge field terms change sign around $\theta\sim 1.3^{\circ}$ with a crossover in the values of $u_{1}$ and $\tilde{u}_{1}$ in the relaxed $\alpha$ and $\beta$ cases. As a result, for moderately large twist angles the complete relaxation effects indicates an overall weakening of the pseudomagnetic fields. We present additional discussions on the relationship between strain and pseudomagnetic fields in Appendix B.

IV IV. Band gaps at the primary and secondary Dirac points

The presence of band gaps in G/BN has been a subject of interest thanks to the prospect of tailoring a high mobility 2D semiconductor based on graphene. Although the LDA predicts the formation of a band gap of $\sim 50$ meV in lattice matched G/BN Giovannetti et al. (2007), the lattice constant mismatch between the graphene and hexagonal boron nitride should suppress the band gap near charge neutrality Ortix et al. (2012) because the spatial average of the local mass term distributed sinusoidally in real space cancels out to a value close to zero Jung et al. (2014). The experimental observation of a band gap at the primary Dirac point came as a surprise Amet et al. (2013); Hunt et al. (2013) and soon it was speculated that the band gap near charge neutrality at the primary Dirac point originates from Coulomb interaction effects Song et al. (2013); Bokdam et al. (2014). Another plausible scenario for the formation of the gap is the generation of an average mass term resulting from the moiré strains Jung et al. (2015) observed experimentally through atomic force microscopy Tang et al. (2013); Woods et al. (2014). These strains originate from stacking registry dependent total energy differences on the order of a few tens of meV Sachs et al. (2011); Jung et al. (2014) that lead to sizeable in-plane strains in the limit of long moiré patterns thanks to the quick decrease of the elastic energy with the increase of the moiré length Jung et al. (2015). The average of the band gap at the primary Dirac cone as a function of strain obtained neglecting higher order contributions in $G$ reads

[TABLE]

with $\widetilde{C}_{R}(\theta)=12\,C_{R}(\theta)\,C_{z}\left(\varepsilon+\theta^{2}(2+\lambda_{g}/\mu_{g})/\varepsilon\right)g^{2}$, where $g$ is the magnitude of the reciprocal lattice vector of graphene. The constant $\Delta_{0}\sim 4\,{\rm meV}$ accounts for the small gap that already develops with only out-of-plane relaxation, and $C_{z}$, $\phi_{z}$ are the moiré parameters for the mass field obtained when height variation is allowed. This result shows that the band gap at the primary Dirac point depends almost linearly to the coefficient $C_{R}(\theta)$ used to quantify the magnitude of the in-plane deformation. Further details of this derivation are presented in Appendix C.

We now turn our attention to the features of the moiré band near the secondary Dirac points observed in G/BN heterojunctions Park et al. (2008); Yankowitz et al. (2012); Ponomarenko et al. (2013). The band structure and the associated density of states (DOS) profile resulting from our model Hamiltonian in Fig. 4 confirm the presence of such features near the mBZ corner at energies of $\sim\pm\hbar\upsilon G/\sqrt{3}$. It was noted that the strong asymmetry between electron and holes commonly observed in experiments with prominent hole features are intrinsic to the band structure model and originate due to opposite chiral winding of the bands combined with the moiré pseudospin terms leading to destructive or constructive contributions of the secondary Dirac cone features DaSilva et al. (2015). A proper modeling of the Hamiltonian parameters for the moiré patterns is necessary to capture the correct band features near the mBZ corners including the sDC shapes and number.

The modifications in the Hamiltonian in the presence of lattice relaxation lead to mild modifications of the sDC features as shown in Fig. 4. Depending on how we account for the relaxation effects in the $\alpha$ and $\beta$ cases, the gap at the secondary Dirac point is an indirect one in the former due to the increase in energy of the valence bands near the $M$ point in the mBZ, whereas in the latter case a Dirac cone shape similar to the rigid band structure is restored. This suggests that the modifications in the electronic structure of graphene due to bond length changes as shown in the $\beta$ solutions can influence details in the optical transitions within the secondary Dirac cones in the terahertz range DaSilva et al. (2015). In the following we discuss how the strain and twist can influence the primary and sDC gaps. In Fig. 5 we represent the magnitude of the gaps as a function of a strain parameter $\eta$ and twist angle $\theta$. While both gaps shrink with growing twist angle, they show opposite behaviors under the influence of an in-plane strain because the physical origin of the gaps are different. In the case of the primary gap we find that $\Delta_{p}$ which is essentially proportional to the global mass $m_{0}$ grows monotonically with the magnitude of in-plane strain. On the other hand, the sDC feature at the mBZ edges is determined by the interplay between all of the first harmonics pseudospin components and the increase of in-plane strain results in reduction of its magnitude. Further insight about the behavior of the gaps at the sDC can be achieved by analysing the bands based on perturbation theory where a triply degenerate band crossing splits under the influence of moiré Hamiltonian parameters Wallbank et al. (2013, 2015); DaSilva et al. (2015). The energy splitting at the mBZ $K$-point can be neatly expressed in the following analytical form:

[TABLE]

where the term

[TABLE]

can be viewed as a sum of three complex numbers that define the three pseudospin components of the moiré pattern. The $m\in\{0,1,2\}$ indices represent the bands on the secondary valley $K$ ($\kappa=1$) or $K^{\prime}$ ($\kappa=-1$) of the mBZ, while $\zeta=1$ ($\zeta=-1$) labels the conduction (valence) band.

The way that the band triplet of Eq. (24) splits in the seondary valleys is schematically represented in Fig. 6. The bands can be described as singlet-doublet structure when only inversion-symmetric Hamiltonians with $\tilde{u}_{i}=0$ is used Wallbank et al. (2013), while in the limit where inversion-asymmetric terms are dominant the triplets are separated by approximately equal gaps Abergel et al. (2013). Defining the sDC gap $\Delta_{s,\zeta\kappa}$ as the energy difference at the Dirac point of the mBZ between bands that are nearer to the primary Dirac cone they can be written as

[TABLE]

where the phase factor $\phi_{\zeta\kappa}=\arg[-\zeta\Delta_{\zeta\kappa}^{3}]/3$ is defined in such a way that $\phi_{\zeta\kappa}\in(-\pi/3,\,\pi/3]$. The above equation for the gap $\Delta_{s,\zeta\kappa}$ is a generally valid expression independent of the choice on coordinate reference for the moiré patterns. The triplets exhibit a periodicity in phase of $2\pi/3$ and the form of $\Delta_{\zeta\kappa}$ is consistent with the possibility of having three different sets of parameters with equivalent electronic structures where the only difference is the choice of the reference origin for the moiré patterns as we discussed earlier in Eqs. (11,12). In the case we choose the Hamiltonian parameter set that already satisfies $-\pi/3<\arg[-\zeta\Delta_{\zeta\kappa}]\leq\pi/3$, the secondary gap can be expressed as the imaginary part of $\Delta_{s,\zeta\kappa}$ involving only inversion asymmetric Hamiltonian terms $\tilde{u}_{i}$

[TABLE]

We note that the expression for the sDC gap in Refs. Abergel et al. (2013); Wallbank et al. (2015) did not impose the restriction in the phase of $\Delta_{\zeta\kappa}$ required for a correct description of the results and whose parameter sets can have phases rotated by $\pm 2\pi/3$ depending on the choirce in the reference point.

V V. Summary

We have presented a moiré band Hamiltonian for G/BN within a framework that uses the symmetry of the moiré patterns and relies on input from ab-initio theories, unifying notation used in the literature for the moiré pattern models in the first harmonics approximation. Our G/BN model accounts for the average mass term that develops in the presence of commensuration strains when the lattice is allowed to relax. We have calculated the lattice relaxation by using a continuum elastic model of graphene for several misalignment angles based on simplifying assumptions where the strain fields are represented within a first harmonics expansion. Taking advantage of the formal simplicity where the strains can be expressed in terms of two parameters, effectively reducing to a magnitude, which is variable, and a constant phase term we have obtained numerically and analytically the twist angle dependence of the moiré pattern Hamiltonian parameters. The presence of strains influences the moiré pattern Hamiltonian parameters first by modifying the real-space distribution of local stacking profiles which in turn change the interlayer coupling maps, and then through intrinsic modifications in the electronic structure of graphene arising due to the bond distortions. These two contributions partially cancel each other resulting in a Hamltonian where the final electronic structure is not too different from the rigid model except for the presence of a band gap at the primary Dirac point. The virtual strain terms in the Hamiltonian due to second order hopping processes of the electrons that are already present for rigid lattices is found to dominate over the relatively small corrections in the Hamiltonian introduced by bond length distortions due to the moiré strains.

The evolution of the band gaps at the primary and secondary Dirac points were studied analytically and as numerically based on the strain models we have developed. Our analysis shows that an overall increase of the commensuration strains, and thus the average mass term, opens up further the primary gap whereas the magnitude of the secondary gap is reduced when strains are larger. Both gaps are found to progressively decrease with the increase of twist angle due to the impact the shortening of the moiré pattern period has on the electronic structure. The primary gap reduces due to the quick decrease of the strain magnitudes, whereas the reshaping of the bands near the mBZ corners that reduce the secondary gap is more strongly influenced by the decreases of the virtual strain terms rather than by the modifications in the Hamiltonian introduced by commensuration strains.

VI Acknowledgement

E.L. acknowledges helpful discussions with J.N.B. Rodrigues. This work is supported by Korean NRF through the Grant NRF 2016R1A2B4010105 (J.J.), EPSRC Grant EP/L013010/1 (M.M.K.), US DoE grant DE-FG02-ER45118 and Welch Foundation TBF1473 (A.M.D. and A.H.M.), and National Research Foundation of Singapore Grant NRF-NRFF2012-01 (S.A.).

Appendix A APPENDIX A: Analytical mapping of the ab-initio parameters in rigid G/BN

The ab-initio results for the moiré pattern in Ref. Jung et al. (2014) were calculated using input from short period commensurate structures with the minimum size of unit cell within self-consistent LDA using an equal lattice constant $2.439\,\,\mbox{\AA}$ for both graphene and hBN. The projection of the pseudospin Hamiltonian components $\widetilde{H}_{\mu}(\vec{d})$ calculated for various displacement vectors $\vec{d}$ are used as input to obtain the first harmonics approximation for the Fourier components of three pairs of parameters $\{C_{\mu},\phi_{\mu}\}$ which correspond to each pseudospin term. These can be related with the $u_{i},\,\tilde{u}_{i}$ parameters Wallbank et al. (2013) through the relations in Eqs. (8-10) whose explicit correspondence is presented in Table 1.

Here we show the analytical mapping between the two parametrization systems. We found that only $u_{1}\upsilon G$ and $\tilde{u}_{1}\upsilon G$ are changed with twist according to $\chi_{\theta}$, see Fig. 2, which can be analytically expressed in the following way for general twist angle

[TABLE]

Three equivalent parameter sets for the same solution that can be obtained using the transformations in Eq. (11-12) are presented in Table 2, and the corresponding plots of the moiré patterns in real space $H_{\mu}(\vec{r})$ are presented in Fig. 7.

The Hamiltonian parameters for rigid graphene presented in Table 2 are equivalent to those in Eq. (40) of Ref. Jung et al. (2014) where the pseudospin magnitudes have been defined as positive numbers and their phases have been changed accordingly. In addition, we have also redefined $\phi_{xy}$ such that it relates to $\phi_{AB}$ in Ref. Jung et al. (2014) through $\phi_{xy}=\pi/6-\phi_{AB}$ which allows to represent the $H_{AB}$ term in Eq. (38) of Ref. Jung et al. (2014) in the compact form given in Eq. (6) in the main text of this work.

Appendix B APPENDIX B: Graphene in-plane relaxation model and strain-induced pseudomagnetic field

In our analysis on the relaxed structure, we assume the simplest case in which only graphene layer relaxes under the influence of a rigid BN substrate. It was shown that the BN sheet in contact with graphene acquires strains that are almost equal in magnitude Jung et al. (2015) and are important to correctly account for the modified Hamiltonian coupling parameters. This additional effect can be included as a later correction and here we exclude it from the self-consistent calculation to give preference to the formal simplicity. The graphene layer is modeled using the Born-von Karman plate theory, in which the elastic properties of graphene are represented by the Lamé parameters $\lambda_{g}$ and $\mu_{g}$ Washizu (1984); Zakharchenko et al. (2009). We also approximate the relaxation by decoupling the in-plane components from the out-of-plane components, so that graphene in-plane relaxation allows the system to minimize its energy with respect to the rigid structure as given by the following functional

[TABLE]

and the elastic energy density $U_{e}$ can be fully expressed in terms of the symmetric tensor in-plane components within linear approximation $u_{ij}=\frac{1}{2}[\partial_{i}u_{j}+\partial_{j}u_{i}]$, where the potential energy density $U_{p}$ is purely a function of the displacement vector $\vec{u}(\vec{r})$,

[TABLE]

[TABLE]

where $A_{g}=3\sqrt{3}a_{0}^{2}/2$ is graphene’s unit cell area. The $\widetilde{f}_{1}(\vec{d})$ and $\widetilde{f}_{2}(\vec{d})$ functions are similar to the first harmonics functions $f_{1}(\vec{r})$ and $f_{2}(\vec{r})$ but now expressed in terms of the displacement vector $\vec{d}(\vec{r})=\vec{d_{0}}(\vec{r})+\vec{u}(\vec{r})$

[TABLE]

where $\vec{g}_{n}$ are the first shell reciprocal lattice vectors of graphene, similar to Eq. (2). In the absence of strain, $\widetilde{f}_{1}(\vec{d_{0}})$ and $\widetilde{f}_{2}(\vec{d_{0}})$ are equivalent to $f_{1}(\vec{r})$ and $f_{2}(\vec{r})$ respectively, see Eq. (7). It should also be noted that in Eq. (31), we have offset the average energy to zero, and the values for $U_{1}$ and $U_{2}$ can be easily deduced from ab-initio calculations of the potential energy on AA, AB and BA stacking points Jung et al. (2015) through the following relations:

[TABLE]

[TABLE]

[TABLE]

Optimising the functional in Eq. (29) with respect to $\vec{u}$ and its derivative $\partial_{i}u_{j}$ requires solving

[TABLE]

The potential energy profile $U(\vec{u}(\vec{r}))$ as expressed in Eq. (31), describing the interlayer coupling at zero twist is illustrated in Fig. 8(a). Treating the problem within the first harmonics, the term on each side of the equation can be expressed as follows:

[TABLE]

[TABLE]

in which we denote the asymmetric part of the strain tensor as $u_{xy}^{A}=\frac{1}{2}[\partial_{x}u_{y}-\partial_{y}u_{x}]$. We solve Eq. (36) using the following ansatz which respects the moiré periodicity of the system

[TABLE]

where the scalar function $\Phi(\vec{r})$ defined in Eq. (15) can also be written in terms of $C_{1}$ and $C_{2}$ parameters that satisfy

[TABLE]

whose magnitude-phase expression is given by

[TABLE]

The solutions that relate the strain parameters and interlayer potentials are

[TABLE]

where $\rho=(A_{g}g^{2})^{-1}\approx 2.19\times 10^{-2}$ is a term that depends on the lattice constant of graphene. The analytic solutions illustrate the way the twist angle influences the magnitudes $C_{1},\,C_{2}\propto\tilde{\varepsilon}^{-4}$ and the shape distortions encoded in the factor $k_{R}$. Within our estimates, $C_{1}$ and $C_{2}$ are always related by a constant phase $\phi_{R}\approx-171^{\circ}$, in keeping with the value of $\phi_{R}=-51^{\circ}$ in Ref. Jung et al. (2015), where we added a $-2\pi/3$ phase for AA$\rightarrow$BA stacking origin representation change.

We should keep in mind that even in the absence of any relaxation the pseudomagnetic fields are already present in G/BN due to the influence of the substrate on the hopping asymmetryJung et al. (2014, 2015), leading to finite values in $u_{1}$ and $\tilde{u}_{1}$, see Fig. 3. This is expected to give rise to a constant pseudomagnetic field profile despite the changing moiré periodicity, see the dashed lines in Fig. 8(c) where a pseudomagnetic field of 35 T is found on the AA stacking point. The pseudomagnetic fields can be obtained by calculating the curl on the vector potential and results in

[TABLE]

where the relaxation strains for different twist angles can introduce changes in the value of the parameters $u_{1}$ and $\tilde{u}_{1}$ or $C_{R}$ and $\phi_{R}$ upon relaxation. The pseudomagnetic field profile at zero twist in a relaxed structure is plotted in Fig. 8, together with its behavior under twist on some special stacking points. For small twist angles, relaxation enhances the pseudomagnetic fields resulting in $u_{1}$ values that are significantly larger than $\tilde{u}_{1}$ with local maxima at AA and AB stacking points where magnitudes of up to 40 Teslas are expected. This enchancement due to relaxation is present for a very small window in twist angle before it decreases and even counters the underlying pseudomagnetic field in the rigid structure when the twist goes beyond the critical angle $\theta_{c}$ of $1.3^{\circ}$.

Appendix C APPENDIX C: Modifications in the moiré pattern Hamiltonian parameters and the primary gap in the presence of strains

As we have shown in the main text, the linearity in the displacement field with respect to position in rigid G/BN lead to a practically zero spatial average for each pseudospin term $\widetilde{H}_{\mu}(\vec{d}_{0}(\vec{r}))$, thus resulting in a zero global mass, see Eq. (13). Accordingly, the band features in such a system can be fully described by the first harmonics contributions to the Hamiltonian. This is in line with the experimental observation of a nearly vanishing gap on the primary valley of G/BN despite the sublattice asymmetry of BN Ponomarenko et al. (2013). However, structural relaxation modifies the displacement field $\vec{u}(\vec{r})$, resulting in a global mass term $m_{0}$ in order to describe the band features properly.

The changes in the interaction between the carbon atoms with the underlying BN substrates due to the additional displacement $\vec{u}(\vec{r})$ lead to modifications in the first harmonics functions $\widetilde{f}_{1}(\vec{d})$ and $\widetilde{f}_{2}(\vec{d})$. The calculation of the effective moiré couplings under relaxation requires numerical computations of the Fourier components of $H_{\mu,\vec{G}_{j}}$, see Eq. (14) but we can obtain the following analytical approximation in the small strain limit where $|\vec{u}|\ll a$

[TABLE]

[TABLE]

where $\varepsilon^{\prime}=\varepsilon+k_{R}\theta$, and $k_{R}$ was defined in Eq. (42). This approximation allows to decompose the contributions into two different relaxation modes quantified by the coefficients $C_{1}$ and $C_{2}$. The finite average results from the non-zero cancellation of the spatial average of $\widetilde{f}_{1}(\vec{d})$ ($\widetilde{f}_{2}(\vec{d})$) due to the presence of a constant $6\varepsilon^{\prime}C_{1}g^{2}$ ($6\varepsilon^{\prime}C_{2}g^{2}$). Considering the transformations in Eq. (45) and (46) on $H_{z}(\vec{r})$, the mass $m_{0}$ given in Eq. (23) can also be expressed as

[TABLE]

which is expected to grow proportionally to the interlayer coupling strength and the strain magnitude in graphene. The constant $\Delta_{0}$ is the finite average gap that opens when the system is allowed to relax on the $z$ axis only. Experimental observation of the primary gap can thus be understood in the light of Eq. (47), which makes the contributions of out-of-plane and in-plane deformation to the resulting gap more transparent.

Taking into account the additional terms introduced by the potential energy and the gauge field, see Eqs. (20) and (21), the effective Hamiltonian parameters $u_{i}^{\prime}$ and $\tilde{u}_{i}^{\prime}$ can be expressed in terms of the parameters $u_{i}$ and $\tilde{u}_{i}$ in a system with mere corrugation using the following expansion:

[TABLE]

The changes in the average mass and the Hamiltonian parameter modifications have signatures in observable features of the bands such as the primary gap $\Delta_{p}$ and the sDC gap $\Delta_{s}$. The explicit parameter values for the moiré Hamiltonians corresponding to rigid and relaxed configurations for zero twists angle obtained within different approximations are listed in Table 3.