Continuum models for twisted bilayer graphene: the effects of lattice deformation and hopping parameter

TL;DR

This paper investigates how lattice deformation and modified hopping parameters affect the electronic properties of twisted bilayer graphene, revealing complex impacts on flat bands and implications for superconductivity models.

Contribution

It introduces a comprehensive continuum model incorporating layer deformation and advanced hopping parameters, highlighting their significant effects on flat band characteristics and wave functions.

Findings

Flat bands are robust but vary in nature with model details.

Higher interlayer coupling strengths lead to additional Dirac points.

Complex models require more parameters and full in-layer dispersion.

Abstract

We analyze a description of twisted graphene bilayers, that incorporates deformation of the layers due to the nature modern interlayer potentials, and a modification of the hopping parameters between layers in the light of the classic Slonczewski-Weiss-McClure parametrisation. We shall show that flat bands result in all cases, but that their nature can be rather different. We will show how to construct a more general reduction to a continuum model, and show that even though such a model can be constructed, its complexity increases, requiring more coupling parameters to be included, and the full in-layer dispersion to be taken into account. We conclude that the combination of all these effects will have a large impact on the wave functions of the flat bands, and that changes in the detail of the underlying models can lead to significant changes. A robust conclusion is that the natural strength of the interlayer couplings is higher than usually assumed, which causes additional Dirac points to appear for the standard magic angles. This gives rise to a degeneracy at the $\Gamma$ point. Since the appearance of a gap at the $\Gamma$ point is crucial for the construction of the Wannier states which are used in the standard descriptions of superconductivity, such an approach not be robust.

Tables

Table 1. Table 1: Lattice harmonic expansion of the deformation: the values of 𝒖 𝒒 subscript 𝒖 𝒒 \boldsymbol{u}_{\boldsymbol{q}} (in units of Å ) for the points labeled as in Ref.  [ 50 ] .

$n_1,n_2$	LCBOP+KC	AIREBO-M+KC	AIREBO-M+ILP	N&K
(1,0)	(0.00042,0.04972)	(0.00141,0.07689)	(0.00129,0.07302)	(0.,0.02660)
(2,0)	(0.00006,0.00323)	(0.00025,0.01307)	(0.00026,0.01078)	(0.,0.00270)
(2,1)	(-0.0019,0.00347)	(-0.00442,0.00809)	(-0.0051,0.00928)	(-0.00100,0.0017)
(3,0)	(0.00001,0.00015)	(0.00008,0.00272)	(0.00007,0.00191)	(0.,0.00036)
(3,1)	(0.00001,0.00001)	(-0.00051,0.00182)	(-0.00075,0.00269)	(-0.00002,0.00035)
(3,2)	(-0.00005,0.00005)	(-0.00132,0.00153)	(-0.00189,0.00216)	(-0.00028,0.00020)
(4,0)	(0.00001,-0.00002)	(0.00003,0.00064)	(0.00002,0.00039)	–
(4,1)	(0.00004,-0.00014)	(0.,0.00016)	(-0.00013,0.00071)	–
(4,2)	(0.00003,-0.00006)	(-0.0002,0.00038)	(-0.00042,0.00078)	–
(4,3)	(0.00009,-0.0001)	(-0.00015,0.00015)	(-0.00055,0.00053)	–

Table 2. Table 2: Table of parameters used in Eqs. ( III – 11 ) for our two many-body screened hopping models.

parameters	screened-1[35]	screened-2
$V_0$	1.06191 eV	1.06191 eV
${\alpha }_1$	0.476	1.0
${\alpha }_2$	0.295 Å-1	0.295 Å-1
${\alpha }_3$	1.411	1.411
${\beta }_1$	6.811	6.811
${\beta }_2$	0.01	0.01
${\beta }_3$	19.176	20.5

Table 3. Table 4: The results for g i ​ u 𝒒 𝑔 𝑖 subscript 𝑢 𝒒 \frac{g}{i}u_{\boldsymbol{q}} using fourth order perturbation theory for q 𝑞 q vectors in one sixth of space.

(1,0)	$-\alpha ∕128\left(31{\alpha }^3-36{\alpha }^2-16\alpha +64\right)(0,1)$
(2,0)	${\alpha }^2\left(\left(3\beta +1579\right){\alpha }^2+392\alpha -1568\right)∕25088(0,1)$
(2,1)	${\alpha }^2\left(\left(3\beta -1277\right){\alpha }^2-1176\alpha +1568\right)∕75264(\sqrt{3},-3)$
(3,0)	$-{\alpha }^3\left(\alpha \left(3\beta -59\right)+392\right)∕37632(0,1)$
(3,1)	${\alpha }^3∕2458624(\sqrt{3}\left(392\left(19-5\beta \right)+\alpha \left(135{\beta }^2+276\beta -6193\right)\right),\alpha \left(81{\beta }^2+216\beta +30965\right)-392\left(3\beta +95\right))$
(3,2)	${\alpha }^3∕1229312(\sqrt{3}\left(392\left(19+2\beta \right)+\alpha \left(-54{\beta }^2-123\beta -6193\right)\right),\alpha \left(-81{\beta }^2-153\beta +12386\right)+392\left(3\beta -38\right))$
(4,0)	$-1∕512{\alpha }^4(0,1)$
(4,1)	${\alpha }^4∕8479744(7\sqrt{3}\left(655-3\beta \left(\beta +96\right)\right),-3\beta \left(3\beta +340\right)-32095)$
(4,2)	${\alpha }^4\left(3\beta -185\right)∕150528(-\sqrt{3},3)$
(4,3)	${\alpha }^4∕8479744(3\sqrt{3}\left(\beta \left(5\beta +506\right)+4585\right),3\beta \left(9\beta +838\right)-22925)$

Equations

$${\mathbf{r}}_{\parallel \mathbf{i}\sigma }$$

$${\mathbf{u}}^{\sigma }(\mathbf{r})$$

$${\mathbf{u}}_{R_{2\pi j/3}\mathbf{q}}=R_{2\pi j/3}{\mathbf{u}}_{\mathbf{q}}.$$

$$w_{AA}({\mathbf{r}}_{li})=\frac{1}{a^2}{\delta }_{⟨r_{li}{r}_{\bar{l}j}⟩}\left[3{\left({\mathbf{r}}_{li}-{\mathbf{r}}_{\bar{l}j}\right)}^2+\sum _k{\left({\mathbf{r}}_{lik+}-{\mathbf{r}}_{\bar{l}jk+}\right)}^2\right].$$

$$w_{AB}({\mathbf{r}}_{li})$$

$$w=w_{AA}-w_{AB}$$

$$H^{(l)}=t\sum _{⟨ij⟩}{c}_i^{(l)†}{c}_j^{(l)}+t^{\prime }\sum _{⟨⟨ij⟩⟩}{c}_i^{(l)†}{c}_j^{(l)}.$$

$$t(r)=0.4\exp (-a(r-r_0)),$$

$$V_1({\mathbf{r}}_1^{(1)},{\mathbf{r}}_2^{(2)})$$

$$\times (1-\tanh (\xi )),$$

$$\xi$$

$$f(x)$$

$$r_0=3.460\text{A˚}.$$

$$|\bar{\mathbf{k}}p⟩$$

$$+\sum _{\mathbf{m},s_2}{c}_{\mathbf{m},s_2}^{(p)2}(\mathbf{k}){|0⟩}_1\otimes {|\mathbf{k}+{\mathbf{G}}_{\mathbf{m}},s_2⟩}_2.$$

$$\text{innerproduct}\mathbf{r}{\mathbf{k},s_l}_l=\frac{1}{\sqrt{N}}\sum _{{\mathbf{r}}_{sl}}{e}^{i\mathbf{k}\cdot \mathbf{r}}\delta (\mathbf{r}-{\mathbf{r}}_{sl}),$$

$$H=\left(\begin{array}{cc}{H}_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right),$$

$$H_{12}$$

$$V_{12}({\mathbf{r}}_1,{\mathbf{r}}_2)\approx V_{12}\left(({\mathbf{r}}_1-{\mathbf{r}}_2)/2\right).$$

$$H_{12}$$

$$=U_{s^{\prime }s}(\mathbf{K},\mathbf{\kappa }).$$

$$U_{s^{\prime }s}(\mathbf{K},\mathbf{\kappa })\approx U_{s^{\prime }s}(\kappa )=U_{s^{\prime }s}({\mathbf{K}}_{\text{av}},\mathbf{\kappa }),$$

$${\tilde{\mathbf{t}}}_{ND}(\mathbf{k})=t\left|2e^{\frac{ia{k}_y}{2}}\cos \left(\frac{1}{2}\sqrt{3}a{k}_x\right)+e^{-ia{k}_y}\right|.$$

$${\mathbf{a}}_1=R_{-\theta /2}(1,0),{\mathbf{a}}_2=R_{-\theta /2}(1/2,\sqrt{3}/2).$$

$${\tilde{\mathbf{a}}}_i=R_{\theta }{\mathbf{a}}_i.$$

$${\mathbf{b}}_1=m{\mathbf{a}}_1+n{\mathbf{a}}_2,{\mathbf{b}}_2=(n+m){\mathbf{a}}_1-m{\mathbf{a}}_2,$$

$$\theta ={\cos }^{-1}\left(\frac{m^2+n^2+4mn}{2m^2+2n^2+2mn}\right).$$

$${\mathbf{b}}_1=n{\tilde{\mathbf{a}}}_1+m{\tilde{\mathbf{a}}}_2,{\mathbf{b}}_2=(m+n){\tilde{\mathbf{a}}}_1-n{\tilde{\mathbf{a}}}_2,$$

$$T{\mathbf{r}}_{kl}^{(0)+}={\mathbf{r}}_{k,-l}^{(0)-},$$

$$r_{kl}^{(0)+}\to r_{-k+3(m+n)/2,l+3(m-n)/2}^{(0)-}.$$

$${\bar{\mathbf{a}}}_i=({\mathbf{a}}_i+{\tilde{\mathbf{a}}}_i)/2=(R^{-1}+R)/{\tilde{\mathbf{a}}}_i/2=\cos (\theta /2)R^{1/2}{\mathbf{a}}_i,$$

$${\mathbf{r}}_{kl}^{\pm }={\mathbf{r}}_{kl}^{(0)\pm }+{\mathbf{u}}^{\pm }({\mathbf{r}}_{kl}^{(0)\pm })={\mathbf{r}}_{kl}^{(0)}\pm \frac{1}{2}\delta r_{kl}^{(0)}+{\mathbf{u}}^{\pm }({\mathbf{r}}_{kl}^{(0)\pm }),$$

$${\mathbf{G}}_i=\frac{1}{m-n}(1-R(\theta )){\mathbf{g}}_i$$

Full text

I Introduction

The discovery of strong interactions and superconductivity in twisted graphene bilayers has been one of the main achievement in two-dimensional materials in the past year; it has been chosen as the Physics World breakthrough of the year 2018 [1, 2, 3], see also Ref. [4]. This field has grown so rapidly that it now carries its own dedicated label, “twistronics”. Twisted graphene layers show a rich phenomenology, likely due to the interplay of a complex electronic structure and the effects of electron interactions. The core ideas build on previous work on the behavior of graphene superlattices on a BN substrate, see for example Refs. [5, 6, 7, 8, 9, 10]. In all of these cases we have a periodic, long wavelength, Moiré modulation, but for graphene on BN the mismatch in lattice spacing between the different materials in the layers limits the maximum wavelength, and thus the diversity of electronic structures for the accessible modulations [11, 12, 13]. On the other hand, the two graphene layers in a twisted bilayer have the same spacing and the periodicity of the Moiré structure has no limit, and diverges at small twist angles [14, 15, 16, 17, 18, 19], L_{M}=d/\bigl{(}2\sin(\theta/2)\bigr{)}, where $d\approx 2.42\,\text{\AA}$ is the lattice unit of graphene. For sufficiently small angles almost flat bands arise near the charge neutrality point [20, 16, 17]. The effects of the intrinsically small interaction effects in graphene are expected to be enhanced for special ‘magic’ angles where the width of the low energy bands is smallest. Novel magnetic phases become possible when the lowest band is half filled [21]. Layer dependent strains can also lead to Moiré structures and narrow bands [22, 23].

When we (almost) align two graphene layers, we have two minimum energy options as shown in Fig. 1. We can either replicate the two layers with only a change in the height ($AA$ alignment), or we can translate one of the layers over a single nearest neigbor distance, which gives $AB$ alignment. In areas with $AB$ alignment half of the carbon atoms in one layer align with those of the other one, but the other half aligns with the midpoints of the hexagons in the other layer. This situation has a lower energy than that with $AA$ alignment. If we consider a twisted bilayer, where both layers are perfectly hexagonal but rotated by an angle relative to a common axis, we find areas with both alignments that are of the same size. At a small cost, the graphene layers can warp, both in and out of plane, to enlarge the beneficial effect of the $AB$ alignment. Doing a fully microscopic calculation (which in this case would require a computationally extremely expensive Green’s function Monte Carlo analysis, since density functional theory calculations struggle to describe bilayer graphene [24], see also[25, 26]) is out of the question for the more than $10,000$ carbon atoms that are contained in a single unit cell, so we need to fall back to simpler models. A few DFT studies are available in the literature[27, 25, 26], although it does not (yet) seem feasible to carry out calculations at the size required to deal with small twist angles.

We can use elegant and simple continuum models when we have no deformation [14, 17], or we can use semi-microscopic atomistic models, such as classical force models for the interatomic forces, both within each layer and between different layers, combined with tight-binding methods for the electronic structure. As we shall discuss below, this latter approach, which relies implicitly on many-body interactions, is likely to give the most realistic description.

At the same time we need to ask ourselves what is the “best” tight-binding description for the electronic structure: For a single layer of graphene the standard approach is to use a nearest-neighbor hopping, and maybe a next nearest neighbor one, to describe the spectra. That approach work very well, even for systems with deformed lattices (typically Moiré supercells). The structure of classical potential models that describe the atomic positions of the atoms in a 2D layer is well understood, and most modern potential models describe the structure of graphene both near and far from equilibrium very well.

The description of both the binding of a bilayer, and the electronic hopping between the layers is much more challenging. The most realistic potential models contain complex many-body interactions, that are necessary to describe the complexities of intra- and inter-layer binding. It is also reasonably well established that one must include many-body effects in the hopping parameters for both graphite and graphene. The key signature of the problems with a two-body description is the difference between nearest-neighbor hopping parameters for different positions in an $AB$-aligned the lattice. As originally described for graphite in the Slonczewski-Weiss-McClure (SWM) model [28, 29, 30] the hopping parameter $\gamma_{1}$, between vertically displaced carbon atoms in $B_{1}A_{2}$ alignment, which has a value of about $0.4\,\text{eV}$ in graphene [31, 32, 33], differs strongly from the two hopping parameters for slightly larger distances: $\gamma_{4}=0.04-0.15\,\text{eV}$ for the hopping near vertical alignment ($B_{1}B_{2}$, etc.), and $\gamma_{3}=0.3\,\text{eV}$ for midpoint aligned carbon atoms (usually labelled $A_{1}B_{2}$). In Fig. 2 we show how $\gamma_{4}$ ($B_{1}B_{2}$) occurs next to $\gamma_{3}$, and that both have the same hopping distanc. Nevertheless $\gamma_{4}$ is much smaller than $\gamma_{3}$ in graphite, which is not captured by the standard distance-dependent two-center Koster-Slater hopping. As discussed in a recent review [34], for bilayer graphene there is a spread in the values found and used. The consensus is that the value of $\gamma_{4}$ is still substantially smaller than $\gamma_{3}$, see also Ref. [18]. A useful form of a model where the screening is dominated by in-layer nearest-neighbor atoms is given in Ref. [35], see also Refs. [36, 37]. This is very similar to the case of the interatomic potentials, which also require a many-body screening largely dominated by nearest-neighbors. Clearly both the graphene lattice deformation and the many-body effects in the hopping will play an important role in describing the band structure obtained in a tight-binding model.

Once we have determined the atomic positions and the hopping parameters for the tight-binding model, we need to deal with the large dimensionality which arises from the size of the unit cell, which leads to a large number of bands. Especially for small angles and thus long Moiré wavelengths, the matrices become extremely large. However, these matrices are very sparse and can be dealt with sparse matrix methods such as ARPACK [38]. Even using those methods numerical calculations are still time consuming. Thus, especially if we want to study many-body physics, we would like to reduce the full tight-binding model to a more efficient low-energy effective model. The one usually used is discussed in Refs. [14, 17], but only works for the simplest lattice and hopping parameters. Since we will use a more complex tight-binding model than normally considered and lattice deformation on top of that, we need to more be careful in making this reduction. We shall investigate this in detail, using an approach that incorporates and generalises the ideas of Ref. [39].

In this work we shall study in a holistic way both the effects of lattice deformation and the change in hopping due the change in alignment, which should be contrasted to related work in Refs. [40, 41]. Our calculation of deformation bears some similarity to the work by van Wijk et al [42, 43], but the authors of those references mainly study a single layer on either bulk graphite or hBN. There are a few other papers that take a related approach [44, 45, 46, 47, 48, 49] to lattice deformation, often in a slightly different context. We start out by selecting a few modern potentials for graphene, and will analyze in detail the deformation of the bilayer systems. This will be validated by comparison to experimental results for strain solitons in bilayers, and will also be compared to the results of a simplified technique originally developed by Nam and Koshino [50]. [We shall show in the Appendix that we can get a simple analytical series expansion for this model with minor modifications.] We then analyze the tight-binding model based on these data, and show that the lowest energy bands remain flat in the presence of a lattice deformation. Then we analyse a general way to extract a low-energy model from such data, and discuss potential issues there. In this work we concentrate on the study of lattice relaxation and electronic structure for a twisted sample at a fixed twist angle, $\theta\approx 1.05^{\circ}$. For this angle, the electronic properties depend on the choice of parameters used. In this respect, our analysis is rather different from those which select a given parametrization and modify the angle in order to obtain the narrowest band [51, 40, 26]. Note, finally, that experiments determine the twist angles mostly from measuring the electron density required to fill the bands in the Moiré superlattice. The angles observed in this way need not coincide with the theoretically defined “magic angles” where the Fermi velocity at the $K$ and $K^{\prime}$ points in the superlattice Brillouin Zone vanishes[17]. Also, it may make sense to use other definitions of the magic angles, such as those which lead to the narrowest bands, or to the largest gaps between the lowest states and the next ones.

Our goal is to quantify the uncertainty that exists in the basic description that is used as a starting point in most calculations of novel features in bilayer graphene. We shall not directly draw conclusions which approach is best; this should ideally be resolved by further experimental measurements of the Moiré structure and the local density of states, as should be accessible to STM measurements. We will, however, compare the lattice relaxation results to the measurements from Ref. [52], and show that we can obtain results that are rather comparable to the interface solitons seen in free-standing bilayers.

Finally, we discuss the robustness of results that rely on a particular Wannier function to describe superconductivity [53, 39, 51, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63]. We shall argue that the electron-assisted hopping model of Ref. [64] looks like the most robust way to obtain superconductivity, independent of the unknown details of the model.

II Classical atomistic simulations

In this section we shall investigate the deformation of free-standing graphene bilayers using atomistic potential models. We will employ a small number of well-established potential models, and for calculational simplicity we restrict our attention to those implemented in the LAMMPS package [65].

For a single layer graphene, we shall use AIREBO-M [66] form of the AIREBO potential [67], as well as the LCBOP-I potential [68], all of which work well for graphene. The reason we shall not use the AIREBO is its small equilibrium C-C spacing of $1.397\,\text{\AA}$, unlike the standard value of $1.420\,\text{\AA}$ recovered for the AIREBO-M and LCBOP potentials. The nature of the interlayer interaction is a subtle question; the long-range and many-body nature of these potentials is discussed in Refs. [68, 69, 42, 70, 71]. Most potential models are modifications of models first used for the interaction of graphene and HBN, and there is some indication that that this leads to a small underestimate of the corrugation of the graphene layers [42]. In this work we shall only use the Kolmogorov-Crespi (KC) potential [69] and the interlayer potential (ILP) [70, 71]. Note that in the “overlay” implementations of the ILP and KC potentials used in LAMMPS, the long range part of the AIREBO is switched off, effectively turning these potentials in re-parameterized REBO potentials [72].

A different interlayer potential makes a difference in the results reported below. A detailed comparison between a large variety of choices in Table 1 of Ref. [73], who derive a rather different form for the potential. In some of those more emphasis is placed on the vertical corrugation of bilayers (which is indeed important for the magnitude of interaction, and even though included in our work, may be slightly underestimated due to the nature of the potentials used). Others concentrate on strained graphene bilayers. The work by Jain et al [46] employs a potential that is specifically designed for the out of layer deformation, but may be less well suited to the details of the in-layer deformation. Nevertheless, this reference also contains an interesting discussion of the lattice deformation. Even though in Ref. [39] the importance of the corrugation is strongly emphasized, we shall argue that the in plane deformation of the lattice actually dominates when we take into account the subtleties of interlayer hopping in $AB$ stacking–rather than the pure two-body form used in that reference. Also, we expect vertical corrugation to be suppressed when the two layers are encapsulated within BN, as is the case in most experiments. In all cases we expect the formation of $AB$ and $BA$ aligned regions separated by domain walls (“interface solitons”). This problem is also discussed in Ref. [74] using an analytic description of domain wall formation, but for rectangular domains.

We have performed simulations for a variety of supercell sizes, but will concentrate here on the case of a superlattice with periodicity $32\boldsymbol{a}_{1}+31\boldsymbol{a}_{2}$, with an angle between the two graphene lattices at the “canonical value” of $1.05^{\circ}$, where we can also compare directly to the semi-analytical work by Nam and Koshino [50]. This last approach is discussed in detail, in a simplified version that is susceptible to analytic solution, in the appendix.

We relax the lattice using a single supercell, with the dimensions chosen to contain a graphene bi-layer lattice without deformation. We then relax, using a conjugate-gradient minimization, first the positions within flat layers, followed by a full relaxation of the carbon atoms. We have checked that these results do not depend on the method or specific order of relaxation used.

A useful way to analyze the in-plane deformation of the relaxed layers is to expand the new positions in terms of a lattice harmonics,

[TABLE]

Here $R_{\theta}$ is a 2D rotation over an angle $\theta$, $H_{r}$ is the first sextant of the reciprocal lattice, i.e., the yellow domain in Fig. 18, and $\sigma=\pm$ denotes the top (bottom) graphene layer. The vector $\boldsymbol{r}^{(0)}$ denotes the undeformed graphene position, and the parallel symbol means we only look at the in-plane component. Due to three-fold symmetry we find we only need to specify a fraction of the coefficients,

[TABLE]

In order to compare the size of the $AA$ and $AB$ aligned domains, we construct a measure of alignment, by combining measures for $AA$ and $AB$ alignment. We first define the measure of $AA$ alignment by the function

[TABLE]

Here $l$ labels the layer, $\bar{l}$ denotes the opposite layer, $\langle r_{li}r_{{\bar{l}}j}\rangle$ denotes the atom $j$ closest to atom $i$ but in the other layer, and $\boldsymbol{r}_{\bar{l}jk\sigma}$ denote positions displaced over a single lattice spacing from $\boldsymbol{r}_{li}$ in the same layer, $\boldsymbol{r}_{lik\sigma}=\boldsymbol{r}_{li}+\sigma\delta^{(k)}_{li}$, where $\delta^{(k)}_{li}$, $k=1,2,3$ denotes the three nearest neighbors $l$ of atom $i$. In a similar way we define the quality of $AB$ alignment as

[TABLE]

See Fig. 3 for a graphical representation of these terms. The factors of 3 in front of the terms involving the central atoms ensure that we use six atoms in every expression; they also weigh the central atom more heavily, when they are aligned. The value of $a$ is the graphene nearest-neighbor spacing.

We then use

[TABLE]

as a measure of alignment. We shall combine data from both layers in a single plot. Note that $w$ is extremal for perfect alignment, negative for $AB$ and positive for $AA$ alignment. See Fig. 3 for a graphical explanation of each of the terms.

Clearly most of the results with a sensible in-layer potential (AIREBO-M and LCBOP-I) fall into groups that largely only depend on the interlayer potential: there are small differences, but they are much smaller than the effect of the interlayer potential. Also, the deformation of the Nam and Koshino analytic result is surprisingly small compared to what we find with modern potentials, with a pattern that appears to be somewhat different as well in the structure of the lattice harmonics, see Table 1. The best way to gauge the quality of these results is to look at the width of the strain solitons between the $AB$ and $BA$ regions. According to Ref. [52], this should be in the order of $6\,\text{nm}$ for a shear boundary, which we believe applies here.

In order to see whether we can reproduce such results, we need to look at larger domains (the ones studied in Ref. [52] vary in size, but typical sizes seem to be at least of the order of $50\,\text{nm}$). As we can see in Figs. 6 and 8 the size of the interface soliton saturates, and we obtain values of the width that are in reasonable agreement with Ref. [52]; a bit narrow for the AIREBO-M+ILC calculations, but rather similar to experiment for the LCBOP+KC ones. The latter case also shows some slow growth with cell size, suggesting the results agrees even better with large-size results from Ref. [52].

As shown in Fig. 9, the pattern of atomic positions for the the strain soliton looks very similar to that presented in Ref. [52]. Of course our analysis is based on only on atomic positions, unlike the results in the paper cited, which are obtained either experimentally using an indirect measure of position, or described with substantial modelling of the probe from the position data. Nevertheless, the similarities are striking.

It is well-known from various simulations cited earlier that vertical corrugation of the graphene layers is important. We would expect a slight underestimate of the corrugation for our current choice of potentials. From potential models fitted specifically to reproduce deformation data [46] we would expect a corrugation of about $d_{AA}=0.360\,\text{nm}$ and $d_{AB}=0.335\,\text{nm}$. The values we find are $d_{AB}=0.335\,\text{nm}$ and $d_{AA}=0.351\,\text{nm}$ for the LCBOP+KC calculations, and $d_{AB}=0.336\,\text{nm}$ and $d_{AA}=0.356\,\text{nm}$ for the AIREBO-M+ILC one. This may show a small underestimate of the vertical corrugation, especially for the LCBOP+KC potentials. Since the $AA$ regions are very small, there is little sensitivity of the binding energy to the $AA$ distance, and thus small changes in the binding can have large effects on this distance without changing the in-lattice deformation and the energy balance appreciably.

III tight binding

Having determined the atomic positions, we need to turn our attention to the electronic degrees of freedom, which we describe using a tight-binding approach. We assume that the lattice unit of the Moiré superlattice is much larger than the graphene unit cell. Without deformation, the Moiré unit cell can then be divided into regions with $AA,AB$, and $BA$ stacking, which each occupy a similar fraction of the unit cell. tight-binding calculations suggest that the wavefunctions which describe the approximately flat bands near the neutrality point are then localized within the $AA$ regions[20]. Since this relies on many approximations, this deserves a detailed investigation.

We start from a tight-binding model for a single layer graphene given by

[TABLE]

Since we have allowed for deformation, we in principle have $t\rightarrow t_{ij}$. Since we shall concentrate on the intralayer coupling, and the changes in $t_{ij}$ are actually very small, we take $t_{ij}=t=\gamma_{0}=-2.7\,\text{eV}$111Please note that for some reason the value used in Ref. [39] is about $10\%$ smaller. and for simplicity we shall use $t^{\prime}=0$ (we have checked this makes no appreciable difference to our results). The fact that we use the same value of $t$ independent of lattice deformation is important: it means that the in-plane wave functions, which only depend on the in-plane hopping parameters, are the same as those of the undeformed lattice; this simplifies the calculations, and is not a real restriction since bond-stretching is extremely small, as explained above.

We use three sets of interlayer hopping parameters; first of all a Koster-Slater exponential parametrisation

[TABLE]

where we use $r_{0}$ as the flat-layer average distance as defined in Eq. 12, which means that we cannot use this parametrisation for the deformed lattices, since the $AB$ couplings become too strong due to the shorter inter-layer distance in the $AB$ regions. In principle, we could replace $r_{0}$ by the $AB$ distance, but since it is not clear that this makes sense, we will not do so, but only apply this parametrisation for flat layers at an interlayer distance $r_{0}$. (Note, however, that Ref. [39] appears to have carried out this program).

We use two sets of environmentally dependent (many-body) hopping parameters, both based on the work in [35], who have designed a many-body screening function that is completely saturated by nearest neighbors only. The form we use is (with $\boldsymbol{r}=\boldsymbol{r}^{(2)}_{2}-\boldsymbol{r}^{(1)}_{1}$)

[TABLE]

We choose two sets of parameters; one, called “screened-1”, is essentially the parameter set from Ref. [35] (with minor modifications); in the other one, “screened-2”, a few parameters have been modified to even more closely represents the parameters in the bilayer SWM parametrisation as reported in Ref. [34]. The parameters for these two potentials are given in Table 2, and we study the behavior of the resulting $\gamma_{i}$ as a function of distance in Fig. 10. We see that our “screened-2” potential only has a weak dependence on interlayer spacing, and gives $\gamma_{1}\approx 0.4\text{eV}$, $\gamma_{3}\approx 0.3\text{eV}$ and $\gamma_{4}\approx 0.1-0.2\text{eV}$, in agreement with the values quoted in [34]. The original screened potential is essentially identical for $\gamma_{1}$, has a slightly smaller and more variable $\gamma_{3}$, and has a value of $\gamma_{4}$ more appropriate to graphite. The Koster-Slater coupling has a great sensitivity to the interlayer spacing, which is especially problematic for $\gamma_{1}$, and follows the relation $\gamma_{3}=\gamma_{4}$, where $\gamma_{3}$ is rather small. Also, it has a 6-fold symmetry for the couplings near $\gamma_{3}$, whereas only a threefold symmetry is present.

We use two sets of deformation parameters, LCBOP+Kolmogorov-Crispi (LKC) and AIREBO-M+ILP (AILP). We also study the effect of fixing the separation, keeping the in plane deformation. For such a flat layer, as might be more appropriate when graphene bilayers are each mounted on HBN, we fix the separation of the layers at an average value of

[TABLE]

For the case studied here (with a unit vector of $32\boldsymbol{a}_{1}+31\boldsymbol{a}_{2}$, and a Moiré angle of $1.05^{\circ}$) the distance in the deformed lattices is given in Table 3.

In Fig. 11 we analyze the effect on the spectrum from both the deformation and interlayer coupling. Again, we only show results for a superlattice twist angle $\theta=1.05^{\circ}$, where the length of the superlattice unit vector is $L_{M}=|\boldsymbol{b}_{1}|\approx 134.2\,\text{\AA}$, and the unit cell contains $11908$ carbon atoms.

All of these results are for a regular bilayer without deformation. So what is the effect of deformation?

As we can see in Fig. 11a, for an undeformed flat lattice and the Koster-Slater hopping parameters, we indeed get flat bands. There also is a secondary Dirac point, so we have probably gone a little bit beyond the first magic angle, which for this interaction is slightly larger. Both of the environment-dependent potentials are a bit too long-range for flat layers, and leads to a larger spitting, but still of the order of $40\,\mathrm{meV}$. Adding lattice deformation leads to much more complicated spectra; secondary Dirac points appear in many places, and culminate in the complicated spectra seen in Fig. 11j-m. These still have a high density of states near the Fermi energy, mostly in a range of $\pm 5,\text{meV}$, so are likely to be susceptible to superconducting instabilities. None of these show a gap between the “flat bands” and the remaining states at the $\Gamma$ point. This will be investigated further below, but it seems unavoidable with the strength of the SWM parameters required, unless we look at a larger twist angle: Whenever we have a second Dirac point for the flat-band calculation, we find that bands touch at the $\Gamma$ point.

The tight-binding models shown here, even though for the canonical angle, show that this is not the magic angle as determined by the band structure. Since our results should at least be close to those by Koshino et al [39], who find clear flat bands and a gap, we first investigate what effects reducing the interlayer coupling and reducing the in-layer Fermi velocity have (these authors use a 10% reduced Fermi velocity). Note that experimental STM data seem more consistent with the larger bandwidth, and the larger Fermi velocity[76, 77] ,(see also related capacitance measurements in[78]).

As can be seen in Fig. 12, we see that the most important effect is the scaling of the Fermi velocity, as introduced in Ref. [39], which removes the secondary Dirac point, and opens a gap at the $\Gamma$ point. A reduction in only the interlayer hopping has almost the same effect, but the secondary Dirac point still remains. This also means that no gap opens at the $\Gamma$ point, where a degeneracy remains. Finally, making both changes has an effect that seems very similar to the bands studied elsewhere. From the discussion in this paper, it should become clear that at an angle of $1.05^{\circ}$ this is not the behavior seen; note that only a model with a gap at $\Gamma$ gives the possibility to project on the 2-band Wannier states.

Clearly we could have reached a similar result by choosing a larger alignment angle; again, keep in mind that all calculations have been done at the same angle, but note that the combination of interlayer couplings and Fermi velocities means that in many cases our chosen twist angle is smaller than the first magic angle for those parameters.

We believe that the in-plane deformation is crucial; the out of plane deformation is likely to be suppressed by the encapsulation of the graphene layer by BN.

IV Continuum projection

Most of the work on studying bilayer graphene has been done using the continuum model, using a $\boldsymbol{k}\cdot\boldsymbol{p}$ model expanded around the a point halfway between the nearest layer, i.e., graphene, Dirac points [17]. In most cases a simple symmetric model is used; the main exception is the work of Koshino et al [39], where the effect of the rippling of the graphene layers is used to modify the coupling strength in the $\boldsymbol{k}\cdot\boldsymbol{p}$ model, but with a simple two-body Koster-Slater interlayer hopping only. As explained in the previous section, we probably under-estimate the rippling, but our results also include the effects of lattice deformation and the many-body screening in the hopping, which have a much stronger effect.

In order to avoid confusion, we shall denote the bilayer’s first Brillouin zone points by a bar in the following; unbarred quantities refer to the single layer graphene points. The technique is straightforward, if a little confusing at first. We refer to Fig. 13 for a graphical representation of the edge of the Brillouin zones of the graphene lattice. Since these are slightly twisted, the reciprocal space is also not perfectly aligned, and the $K$ points in the two layers, $K_{1}$ and $K_{2}$, no longer coincide. For small angles these points are relatively close together, and develop a continuum Hamiltonian around the point $\boldsymbol{K}_{\text{av}}$ halfway between the two $K$ points. On folding to the bi-layer graphene Brillouin zone these $K$ points map to inequivalent points $\bar{K}$ and $\bar{K}^{\prime}$ in the bilayer-superlattice Brillouin zone. The point $\boldsymbol{K}_{\text{av}}$ maps to $\bar{M}$. Since the Fermi velocity of graphene is rather large, we expect that only momenta near these two Dirac points play a role. To make that more precise, we expand the bilayer wave functions in products of the states of the graphene layers. The fact that the Moiré pattern is periodic means that only states that differ in momentum by the superlattice reciprocal lattice vectors mix. More precisely, we write for a single electron state of momentum $\boldsymbol{k}$ in the $p$th band:

[TABLE]

Here we choose for convenience $\boldsymbol{k}$ as the “unfolded” momentum corresponding to the momentum $\bar{\boldsymbol{k}}$ in the FBZ, i.e., the equivalent momentum nearest the two Dirac points, and $s_{l}$ is a sublattice index for each layer. The states $\ket{\boldsymbol{k},s_{l}}_{l}$ are the standard plane wave solutions (since we have not modified the in-plane hopping parameters, the positions used here are the undeformed lattice positions)

[TABLE]

where $\boldsymbol{r}_{sl}$ are the positions in sublattice $s$ in layer $l$. The expression (13) is exact, and only becomes approximate on restriction of the superlattice sums. Before we do that, we first look at at the representation of the tight-binding Hamiltonian in this basis.

Clearly we can write a block diagonal form

[TABLE]

where each block itself is a block of $2\times 2$ matrices in sublattice space, with the dimension determined by the number of vectors $\boldsymbol{G}_{M}$ included. Thus $(H_{11})_{\,\boldsymbol{m}s_{1},\boldsymbol{m}^{\prime}s^{\prime}_{1}}=\epsilon_{\boldsymbol{k}+\boldsymbol{G}_{\boldsymbol{m}},s_{1}}\delta_{\boldsymbol{m}\boldsymbol{m}^{\prime}}\delta_{s^{\prime}_{1}\bar{s}_{1}}$ and $(H_{22})_{\,\boldsymbol{n}s_{2},\boldsymbol{n}^{\prime}s^{\prime}_{2}}=\epsilon_{\boldsymbol{k}+\boldsymbol{G}_{\boldsymbol{n}},s_{2}}\delta_{\boldsymbol{n}\boldsymbol{n}^{\prime}}\delta_{s^{\prime}_{2}\bar{s}_{2}}$. For small $\boldsymbol{n}$ and $\boldsymbol{m}$ this is a slightly modified Dirac Hamiltonian (see below). The off-diagonal terms do allow coupling between different momenta due to the periodic Moiré, and the allowed couplings are of the form

[TABLE]

The momentum dependence of the tight-binding Hamiltonian originates from the imposition of periodic boundary conditions, and also from the fact that, even though short-ranged, the interlayer potential, $V_{12}(\boldsymbol{r}_{1},\boldsymbol{r}_{2})$, where $\boldsymbol{r}_{1}$ and $\boldsymbol{r}_{2}$ reside in different layers, is non local.

The standard continuum approximation makes the assumption that the interlayer potential is only significantly different from zero if $|\boldsymbol{r}_{1}-\boldsymbol{r}_{2}|\ll L_{M}$, where $L_{M}$ is the Moiré lattice unit. Then, the position dependence of the interlayer potential should be well approximated by

[TABLE]

When expressed in momentum space, this approximation neglects the dependence on the average momentum, $\boldsymbol{k}+(\boldsymbol{G}_{m}+\boldsymbol{G}_{n})/2$. Even though in most cases considered here this is a small effect, the gaps we observe are also very small, and we would like to take a more careful approach

We will make the approximation $H(\boldsymbol{k})=H(\boldsymbol{k}+(\boldsymbol{G}_{\boldsymbol{m}}+\boldsymbol{G}_{\boldsymbol{n}})/2)$. This is rigorously true only if $(\boldsymbol{G}_{\boldsymbol{m}}+\boldsymbol{G}_{\boldsymbol{n}})/2$ is a superlattice vector. Nevertheless, we write $\boldsymbol{K}=\boldsymbol{k}+(\boldsymbol{G}_{m}+\boldsymbol{G}_{n})/2$ and $\boldsymbol{\kappa}=(\boldsymbol{G}_{m}-\boldsymbol{G}_{n})/2$ and

[TABLE]

The basic idea of the continuum model [14, 17], see also [39], is that for low energy states, and thus momenta near the Dirac points, we can make the approximation that the dependence on the average momentum can be replaced by the momentum at the point $\boldsymbol{K}_{\text{av}}$ halfway between the two $K$ points. This would mean that for momenta near $\boldsymbol{K}_{\text{av}}$ we only consider the following quantity

[TABLE]

which is slightly more satisfying approach to the local-potential approximation. Since the interlayer coupling usually falls of quickly with momentum, $U$ is thus dominated by a few points on the triangular $\boldsymbol{G}_{\boldsymbol{m}}$ lattice [17]. Actually, for reasons not perfectly clear to us, it seems better to use $\boldsymbol{K}_{\text{av}}=\boldsymbol{K}_{1}$ for a low-order truncation to $U$–this preserves the three-fold symmetry normally imposed on the model. We find that even that is not the optimal approximation, as is shown below.

Let us first look at what these matrix elements (18) are for the problems studied previously; we study all of the cases shown in Fig. 11 in Fig. 14. We indeed find that for a Koster-Slater potential and a flat lattice the couplings are dominated by 3 wave vectors (which is the model underlying Refs. [17, 39]). We clearly see that in all cases the three nearest-neighbor vectors dominate but that the decay is slower both due to lattice deformation and the change of the interlayer hopping parameters. For the $AA$ coupling we always find a small asymmetry between the $\boldsymbol{G}=0$ coupling and the other two strong couplings (by a few percent), removing some of the symmetries of the model, which can be restored, see below. For an undeformed graphene lattice and the Koster-Slater hopping parameters (a), the parameters are essentially those quoted in Ref.  [17], after a small rescaling of the strength. The $\gamma_{3}$–$\gamma_{4}$ asymmetry in the remaining results clearly has a big impact. For an undeformed lattice (b/c), we see larger $AB$ than $AA$ couplings, but there is an indication that the coupling decays slightly more slowly, and some additional couplings may be thus be required in the continuum model. For the relaxed and deformed lattices, we find a more substantial difference between the $AA$ and $AB$ couplings, where the $AA$ coupling is smaller (by $15-30\%$) than the $AB$ one. It should come as no surprise that the AILP results, which have the smallest $AA$ regions, show the largest difference.

What we have not shown is the imaginary parts: normally one assumes that the coefficients in $U$ are real after removing a trivial phase-dependence. In our case they seem to develop small but significant imaginary parts.

We shall now apply the expansion of $U$ in two different approaches: Since, due to the large energy cost associated with moving up the Dirac cones, only momenta $\boldsymbol{K}$ near $\boldsymbol{K}_{\text{av}}$ will contribute, it is usually considered sufficient to replace the average momentum dependence by the central value, and expand the graphene dispersion to linear order about this same point. The second idea is based on the fact that we can do better at little cost: for the momenta that are relevant, a linear approximation of dependence of $U$ on $\boldsymbol{K}$ (expanded near $\boldsymbol{K}_{\text{av}}$), can easily be combines with the full in-layer dispersion. We truncate the matrix diagonalization to the $n$th hexagon, and we find that a projection with $n=3-5$ (depending on the range of $U$) is sufficient to reproduce the energy of the flat bands, which is similar to the truncation proposed in the literature; we usually use a few more hexagons to ensure convergence. Slightly more concerning is the effect of an expansion of the Dirac Hamiltonian about the $K^{(1,2)}$ points. In the most complete calculation we use the exact dispersion, described by the off diagonal element of the in-plane Hamiltonian

[TABLE]

In Figs. 15 and 16 we give two examples of calculations for two extreme cases; a complete set is shown in the supplementary material.

Let us look at the “standard case”, Fig. 15 first. we see a rapid convergence of the results with the range of $U$; the three dominant matrix elements are almost sufficient. We see a small symmetry breaking along the $\Gamma$–$M$ lines for the Bistritzer-McDonald calculation. This could have been avoided by replacing $K_{\text{av}}$ by $K_{1}$, and we would get the correct degenerate spectrum, but not the particle-hole asymmetry–i.e., the Fermi energy is incorrect. Interestingly enough, by using the linear $K$-dependence of $U$ and the full dispersion (and both are required) we get a perfect reproduction of the tight-binding spectrum.

The situation gets much more interesting for Fig. 16, which corresponds to Fig. 11j. Clearly even for this complicated case the full calculation converges to something close to the tight-binding results (we could have added probably one more hexagon of couplings, which would have converged). The standard approximation, based on just three harmonics, gives a rater poor approximation.

We now present some results for the continuum projection (a complete set can be found in the supplementary material). We selected two cases of most interest: the SWM model without deformation and the AILP+KC deformation, as probably the most reasonable cases to investigate.

As we can see in Fig. 15, the inclusion of a full dispersion and the the dependence of the intralayer matrix elements is required to get the degeneracy of the energies in the two valleys along the line $\Gamma$-$K$.

This changes in spectrum will clearly also have important consequences for the wave function–which in turn can be used to construct the Wannier functions. These are shown in the supplementary material.

V Conclusions.

We have presented a comprehensive analysis of the lattice relaxation in twisted graphene bilayers, and its effect on the electronic properties, due to the modulation of the interlayer hopping. Calculations have been carried out for a Moiré superlattice with lattice vector $\boldsymbol{L}_{M}=32\boldsymbol{a}_{1}+31\boldsymbol{a}_{2}$, where $\boldsymbol{a}_{1}$ and $\boldsymbol{a}_{2}$ are the unit vectors of the graphene lattice. The twist angle is $\theta\approx 1.05^{\circ}$. Note that our approach is complementary to other studies, where one selects the angle which gives the narrowest bands near the neutrality point, and keeps the parametrization used fixed[51, 40, 26].

The relaxation is calculated using classical interatomic force models, and the electronic states are determined using tight-binding models. We have compared different force models, and different dependencies of the interlayer electronic hopping parameters on atomic positions, and find rather similar results. The relaxed positions of the atoms are used as input for the calculation of the electronic structure, calculated using tight-binding models. Different parametrizations of the couplings are used: i) hoppings between orbitals in different layers which combine a form factor which reflects the symmetry of $p$ orbitals, and a simple exponential dependence on distance, and ii) hoppings that depend on the distance and the local environment of the two orbitals involved in the process. Models of type ii) reproduce the difference between the SWM parameters $\gamma_{3}$ and $\gamma_{4}$ needed to describe aligned bilayers and graphite. For a fixed twist angle, $\theta\approx 1.05^{\circ}$, the low energy bands show a significant dependence on both the range of the interaction and whether the hopping parameters depend solely on interatomic distances, or they also include other features of the environment. To some extent, the results can be interpreted as a parameter dependent shift of the “magic angle”, where the low energy bands are narrowest. When the choice of parameters is such that the magic angle is greater than $1.05^{\circ}$, we find new band crossings and Dirac points[79].

The low-energy electronic bands show a significant dependence on the amount of lattice relaxation and on the dependence of the interlayer hopping parameters on distance and local environment. The bandwidth of the lowest bands at neutrality is probably a bit larger than for the non-relaxed case, but still has a large density of states within a few meVs. However, a number of features, such as the number and location of additional band crossings (Dirac points) and saddle points (van Hove singularities) varies considerably as function of the model being used. The overlap, or lack thereof, of the lowest bands and neighboring bands is also quite sensitive the choice of parameters, within a range of physically sensible ones.

Finally, we have studied the connection between tight-binding and continuum $\boldsymbol{k}\cdot\boldsymbol{p}$ models. We find that the number of harmonics required in a continuum approximation is dependent on the strength of the lattice relaxation and details of the interlayer hopping, but that effective continuum models can be defined in all cases.

We have analyzed the minimal continuum models required to approximate the electronic bands obtained from tight-binding calculations defined at the atomic scale. The complexity of the continuum models depends significantly on the range of the hoppings, and on whether they depend significantly on the local environment. Isotropic couplings which do not decay too abruptly with distance are reasonably described with the standard model based on an expansion with three harmonics of the interlayer hoppings. A continuum description is possible for all tight-binding models considered, although more than three harmonics are required in some cases, especially when the hopping parameters depend on the local environment. Even with a large number of harmonics, such a continuum model can be an effective way to study a tight-binding model, especially when also adding residual interactions. This is of course dependent on the method for extracting the coupling parameters from a tight binding model. This calculation can be done quite simply, and only relies on the construction of a tight-binding Hamiltonian, not its diagonalization.

We have compared results from various models, both for the interatomic forces and for the electronic hopping parameters, using the same twist angle, $\theta=1.05^{\circ}$. This choice is motivated by the fact that the value of the twist angle is the magnitude most accessible experimentally. It is yet unclear how precisely the experimentally studied twist angles correspond to the theoretical definition of magic angles. The dependence found here of the electronic properties on the choice of parameters suggests that the observed tendency towards broken symmetry phases must be quite robust. The appearance of superconductivity and insulating behavior in twisted graphene bilayers is likely to arise from rather general properties of the models.

Acknowledgements.

FG was supported by the European Commission under the Graphene Flagship, contract CNECTICT-604391; NRW is supported by UK STFC under grant ST/P004423/1.

Appendix A Analytical model for lattice deformation

Here we derive an analytic expression for the elastic deformation of bilayer graphene, based on the work by Nam and Koshino [50].

We assume that the lattice vectors of the two unperturbed graphene lattices, which are rotated by a relative angle $\theta$, for each layer are given by (from now we use the graphene lattice spacing, $1.42\sqrt{3}\,\text{\AA}$, as a length unit)

[TABLE]

and for the second layer we have

[TABLE]

The lattice vectors of the super cell are

[TABLE]

and the angle between the two layers can be expressed as

[TABLE]

We can also express this in terms of $\tilde{\boldsymbol{a}}_{i}$, where $m$ and $n$ change roles:

[TABLE]

In the remainder we shall always implicitly assume that the angle $\theta$ is small (normally we will only consider the case $m=n+1$ where $\theta=\sin^{-1}\left(\frac{\sqrt{3}(2m-1)}{6m^{2}-6m+2}\right)\approx 1/(\sqrt{3}m)$). We will denote $R_{\theta/2}$ as $R$.

There is a symmetry between the layers, as can be seen in Fig. 17. We label the layers by $+$ (top) and $-$ (bottom). It is easy to show that with the lattice positions given by $\boldsymbol{r}^{(0)+}_{kl}=k\boldsymbol{a}_{1}+l\boldsymbol{a}_{2}$ and $\boldsymbol{r}^{(0)-}_{kl}=k\tilde{\boldsymbol{a}}_{1}+l\tilde{\boldsymbol{a}}_{2}$ we have an additional symmetry under reflection $T_{x}$ in the $x$-axis,

[TABLE]

We have a similar symmetry for reflections in the line connecting $\boldsymbol{b}_{1}$ to $\boldsymbol{b}_{2}$, Without writing down the detailed form of the transformation matrices, we see that this maps

[TABLE]

We now assume that the two lattices will deform in a similar way, respecting the reflection symmetry. If we define an average lattice by the vectors

[TABLE]

We can write for the lattice vectors in the two lattices, labeled as $\pm$,

[TABLE]

with $\boldsymbol{r}_{0}=k\bar{\boldsymbol{a}}_{1}+l\bar{\boldsymbol{a}}_{2}$. If we assume $\delta r^{(0)}_{kl}$ is small222That is not true over the whole supercell; the misalignment is one lattice spacing at the far corner of the supercell. Fortunately, that is where $r$ is large as well., then we can make the approximation that $\boldsymbol{u}^{\pm}(\boldsymbol{r}^{(0)\pm}_{kl})=\pm\boldsymbol{u}(\boldsymbol{r}^{(0)}_{kl})$, and we can simplify this expression. We use $kl$ as labels to show that their range is either $(m,n)$ or $(n,m)$, depending on the layer. As we can see from Fig. 17, this makes most sense in half the Brillouin zone; we can, however, work with the hexagonal Brillouin zone where this approach works well everywhere.

We define the three reciprocal lattice vectors $\boldsymbol{g}_{i}$ to $\boldsymbol{a}_{j}$, and similar for $\tilde{\boldsymbol{a}}_{j}$. We then define the superlattice reciprocals,

[TABLE]

It is straightforward to see that $\boldsymbol{G}_{i}\cdot(\boldsymbol{g}_{i}+\tilde{\boldsymbol{g}}_{i})=0$. [Note the slightly awkward labeling: $\boldsymbol{G}_{1}$ and $\boldsymbol{G}_{3}$ are the dual vectors to $\boldsymbol{b}_{1}$ and $\boldsymbol{b}_{2}$.]

We now minimize the combination of the misalignment of the lattices and the elastic energy as done by Nam and Koshino, with a minor change in the vectors used in the misalignment energy, assuming that we can write the continuum approximation (notice that here there is an important difference with Koshino, who have no reference to the mean displacements, but work in one of the two sub-lattices, so the meaning of $\boldsymbol{r}$ is very different, and their final results lacks the layer symmetry found below)

[TABLE]

where $\boldsymbol{\delta}(\boldsymbol{r})$ is a field in the average lattice, with $\boldsymbol{\delta}(\boldsymbol{r})\cdot\boldsymbol{r}=0$. Since the $\boldsymbol{\delta}$ is the vector from the top to the bottom lattice, we would like to align this displacement with the favourable positions for the top lattice, but then we would like to align $-\boldsymbol{\delta}$ with the bottom lattice. Thus we see we need to minimise the potential

[TABLE]

We find that, using the average $\bar{\boldsymbol{g}}_{j}=(\boldsymbol{g}_{j}+\tilde{\boldsymbol{g}}_{j})/2$,

[TABLE]

Thus,

[TABLE]

We can now follow Nam and Koshino, and the standard continuum elastic energy to the energy derived here. This lead to the requirement to solve the coupled equations, where $\boldsymbol{q}_{\perp}=(q_{y},-q_{x})$:

[TABLE]

If we make the simplest approximation for the sine, neglecting completely the contribution from $\boldsymbol{u}$, we find that

[TABLE]

and thus, since $\boldsymbol{G}_{j}$ and $\bar{\boldsymbol{g}}_{j}$ are orthogonal, we find that ($(1)$ for first order)

[TABLE]

Since the two vectors $(G^{j}_{y},-G^{j}_{x})$ and $\bar{\boldsymbol{g}}_{j}$ are parallel, this can be written as

[TABLE]

Thus we find that the dimensionless quantity

[TABLE]

The expansion parameter $\alpha=\frac{4V_{0}g^{2}}{\mu G^{2}}$ grows with the size of the unit cell, showing that for very small angles a perturbative approach must fail.

With the help of a simple mathematica code it is now straightforward to find the higher order terms, which involves expanding Eq. (37) to higher order in $V_{0}/\mu$. Results following the notation by Nam and Koshino are given in Table. 4. Our results are a universal (lattice-size independent) expression when we scale $u_{\boldsymbol{q}}$ as $\frac{g}{i}u_{\boldsymbol{q}}$, and express the values in terms of the parameters

[TABLE]

Here we use $g^{2}=\frac{8\pi^{2}}{3a^{2}}(1+\cos\theta)$ and $G^{2}=\frac{64\pi^{2}}{3a^{2}}\sin^{2}\theta/2$.

When using this for finite discrete lattices, we shall use $\boldsymbol{r}^{(0)\pm}$ as the argument of $\boldsymbol{u}$, which restores the broken reflection symmetry.