Lattice relaxation and energy band modulation in twisted bilayer graphenes

TL;DR

This paper develops a continuum theory to analyze lattice relaxation in twisted bilayer graphene and demonstrates how relaxation significantly alters electronic properties, including band gaps and Fermi velocity enhancements.

Contribution

It introduces a new effective continuum model for lattice relaxation in TBG and reveals its impact on electronic band structure, especially at small twist angles.

Findings

Relaxed lattice reduces AA-stacking area and forms a triangular domain structure.

Energy gaps up to 20meV open at superlattice subband edges due to relaxation.

Lattice relaxation enhances the Fermi velocity compared to non-relaxed models.

Abstract

We theoretically study the lattice relaxation in the twisted bilayer graphene (TBG) and its effect on the electronic band structure. We develop an effective continuum theory to describe the lattice relaxation in general TBGs and obtain the optimized structure to minimize the total energy. At small rotation angles $< 2^{\circ}$, in particular, we find that the relaxed lattice drastically reduces the area of AA-stacking region, and form a triangular domain structure with alternating AB and BA stacking regions. We then investigate the effect of the domain formation on the electronic band structure. The most notable change from the non-relaxed model is that an energy gap up to 20meV opens at the superlattice subband edges on the electron and hole sides. We also find that the lattice relaxation significantly enhances the Fermi velocity, which was strongly suppressed in the non-relaxed model.

Equations

$$\mathbf{\delta }(\text{r})=\text{r}-{\text{r}}_0=(1-R^{-1})\text{r}.$$

$${\mathbf{L}}_i^{\mathrm{M}}=(1-R^{-1}{)}^{-1}{\mathbf{a}}_i(i=1,2).$$

$$L_{\mathrm{M}}=\frac{a}{2\sin (\theta /2)}.$$

$${\mathbf{G}}_i^{\mathrm{M}}=(1-R){\text{a}}_i^{\ast }={\text{a}}_i^{\ast }-{\tilde{\text{a}}}_i^{\ast }.(i=1,2),$$

$$L=|m{\mathbf{a}}_1+n{\mathbf{a}}_2|=a\sqrt{m^2+n^2+mn}=\frac{|m-n|a}{2\sin (\theta /2)},$$

$$\cos \theta =\frac{1}{2}\frac{m^2+n^2+4mn}{m^2+n^2+mn}.$$

$$\delta (x)=x-\frac{N-1}{N}x=\frac{a}{L_{\mathrm{M}}}x.$$

$$\delta (x)={\delta }_0(x)+u_2(x)-u_1(x),$$

$$U_E=\int \frac{1}{2}\kappa \left[{\left(\frac{\partial u_1}{\partial x}\right)}^2+{\left(\frac{\partial u_2}{\partial x}\right)}^2\right]\mathrm{d}x,$$

$$V[\delta ]=-2V_0\cos a^{\ast }\delta ,$$

$$U_B=\int V\left[\delta (x)\right]\mathrm{d}x.$$

$$\kappa \frac{{\partial }^2{u}_{+}}{\partial x^2}=0,$$

$$\kappa \frac{{\partial }^2{u}_{-}}{\partial x^2}-4a^{\ast }{V}_0\sin (G_{\mathrm{M}}x+a^{\ast }{u}_{-})=0.$$

$$u_{-}(x)=\sum _n{u}_n{e}^{in{G}_{\mathrm{M}}x}$$

$$\sin (G_{\mathrm{M}}x+a^{\ast }{u}_{-}(x))=\sum _n{f}_n{e}^{in{G}_{\mathrm{M}}x}.$$

$$u_n=-\frac{4a^{\ast }{V}_0}{\kappa (n{G}_{\mathrm{M}}{)}^2}{f}_n.$$

$$\frac{u_n}{a}=-\frac{2{\eta }^2}{\pi n^2}{f}_n,$$

$$\eta =\sqrt{\frac{V_0}{\kappa }}\frac{L_{\mathrm{M}}}{a}.$$

$$n\text{lower}-1.29167pt<\text{lower}3.01385pt\sim \sqrt{\frac{2}{\pi }}\eta .$$

$$\kappa \frac{{\partial }^2{u}_{-}}{\partial x^{{}^{\prime }2}}+4a^{\ast }{V}_0\sin (a^{\ast }{u}_{-})=0,$$

$$u_{-}(x^{\prime })=\frac{2a}{\pi }\arctan \left[\exp \left(4\pi \sqrt{\frac{V_0}{\kappa }}\frac{x^{\prime }}{a}\right)\right]-\frac{a}{2},$$

$$w_{\mathrm{d}}\approx \frac{a}{4}\sqrt{\frac{\kappa }{V_0}}.$$

$$\frac{w_{\mathrm{d}}}{L_{\mathrm{M}}}=\frac{1}{4\eta },$$

$$\mathbf{\delta }(\text{r})={\mathbf{\delta }}_0(\text{r})+{\text{u}}^{(2)}(\text{r})-{\text{u}}^{(1)}(\text{r}),$$

$$\begin{gathered} U_E=\sum _{l=1}^2\int \frac{1}{2}\{(\lambda +\mu )(u_{xx}^{(l)}+u_{yy}^{(l)}{)}^2 \\ +\mu \left[(u_{xx}^{(l)}-u_{yy}^{(l)}{)}^2+4(u_{xy}^{(l)}{)}^2\right]\}{d}^2\text{r} \end{gathered}$$

$$V[\mathbf{\delta }]=\sum _{j=1}^{3}2V_0\cos [{\text{a}}_j^{\ast }\cdot \mathbf{\delta }],$$

$$U_B$$

$$V[\mathbf{\delta }(\text{r})]=\sum _{j=1}^{3}2V_0\cos [{\text{G}}_j^{\mathrm{M}}\cdot \text{r}+{\text{a}}_j^{\ast }({\text{u}}^{(2)}-{\text{u}}^{(1)})],$$

$$\begin{gathered} \frac{1}{2}(\lambda +\mu )\left(\frac{{\partial }^2{u}_x^{-}}{\partial x^2}+\frac{{\partial }^2{u}_y^{-}}{\partial x\partial y}\right)+\mu \left(\frac{{\partial }^2{u}_x^{-}}{\partial x^2}+\frac{{\partial }^2{u}_x^{-}}{\partial y^2}\right) \\ +\sum _{j=1}^{3}2V_0\sin [{\text{G}}_j^{\mathrm{M}}\cdot \text{r}+{\text{a}}_j^{\ast }\cdot {\text{u}}^{-}]g_j^x=0, \end{gathered}$$

$$\begin{gathered} \frac{1}{2}(\lambda +\mu )\left(\frac{{\partial }^2{u}_y^{-}}{\partial y^2}+\frac{{\partial }^2{u}_x^{-}}{\partial x\partial y}\right)+\mu \left(\frac{{\partial }^2{u}_y^{-}}{\partial x^2}+\frac{{\partial }^2{u}_y^{-}}{\partial y^2}\right) \\ +\sum _{j=1}^{3}2V_0\sin [{\text{G}}_j^{\mathrm{M}}\cdot \text{r}+{\text{a}}_j^{\ast }\cdot {\text{u}}^{-}]g_j^y=0, \end{gathered}$$

$${\text{u}}^{-}(\text{r})=\sum _{\text{q}}{\text{u}}_{\text{q}}^{-}{e}^{i\text{q}\cdot \text{r}},$$

Full text

I Introduction

Twisted bilayer graphene (TBG) is a two-dimensional material where two graphene layers are relatively rotated by an arbitrary angle. Such a rotational stacking structure has been widely observed in epitaxially-grown multilayer graphenes,Berger et al. (2006); Hass et al. (2007, 2008); Li et al. (2009); Miller et al. (2010); Luican et al. (2011); De Heer et al. (2010) and also it can be fabricated by manually aligning single layers using the transfer technique. Wang et al. (2013); Tan et al. (2016) The electronic properties of TBG has also been intensively studied in theory, where it was shown that the energy spectrum sensitively depends on its rotation angle $\theta$. Lopes dos Santos et al. (2007); Mele (2010); Trambly de Laissardière et al. (2010); Shallcross et al. (2010); Morell et al. (2010); Bistritzer and MacDonald (2011a); Kindermann and First (2011); Xian et al. (2011); Lopes dos Santos et al. (2012); Moon and Koshino (2013, 2012); Bistritzer and MacDonald (2011b) In a small $\theta$, in particular, the interference between the incommensurate lattice structures gives rise to a moiré pattern with a long spacial period, and it significantly influences the low-energy spectrum. For TBG less than a few degree, in particular, the Dirac cones of the two layers are strongly hybridized by the moiré interlayer interaction, where the linear dispersion is distorted into superlattice subbands. Lopes dos Santos et al. (2007); Mele (2010); Trambly de Laissardière et al. (2010); Shallcross et al. (2010); Morell et al. (2010); Bistritzer and MacDonald (2011a); Kindermann and First (2011); Xian et al. (2011); Lopes dos Santos et al. (2012); Moon and Koshino (2013, 2012); Bistritzer and MacDonald (2011b). The characteristic features of superlattice band structure of TBG were actually observed in recent experiments. Kim et al. (2016); Cao et al. (2016); Kim et al. (2017)

Most of the band calculations for TBG assumes that the two graphene layers are rigid and simply stacked without changing the original honeycomb lattices. In a real system, however, the lattice structure spontaneously relaxes to achieve an energetically favorable structurePopov et al. (2011); Uchida et al. (2014); van Wijk et al. (2015); Dai et al. (2016); Jain et al. (2016), and it should influence the electronic spectrum. If we consider a TBG with a small rotation angle as in Fig. 1, for instance, we notice that the lattice structure locally resembles the regular stacking such as AA, AB or BA, depending on the position. Here AA represents the perfect overlapping of hexagons, while AB and BA are shifted configurations in which A(B) sublattice is right above B(A). Since the interlayer binding energy is the lowest in AB and BA and the highest in AA stacking Popov et al. (2011); Lebedeva et al. (2011); Gould et al. (2013), the TBG spontaneously deforms so as to maximize the AB/BA areas while minimize AA area. In fact, such an AB/BA domain structure was experimentally observed in multilayer graphenes grown by chemical vapor depositionBrown et al. (2012); Lin et al. (2013); Alden et al. (2013), and also captured in the theoretical calculations, Popov et al. (2011); van Wijk et al. (2015); Dai et al. (2016); Jain et al. (2016) while its implication on the electronic band structure is still unclear. Similar lattice relaxation was also found in another moiré superlattice of graphene on hexagonal boron-nitride, Dean et al. (2013); Ponomarenko et al. (2013); Hunt et al. (2013); Yu et al. (2014); Chen et al. (2014); Moon and Koshino (2014) where the sublattice inequality in hBN results in a hexagonal domain pattern. Woods et al. (2014); San-Jose et al. (2014a, b); Jung et al. (2015)

In this work, we present a theoretical study on the lattice relaxation in TBG and its effect on the electronic band structure. First, we develop a minimum phenomenological model to describe the AB/BA domain formation in TBGs with general angles. Using the elastic theory and a simple interlayer adhesion potential, we express the total energy as a functional of the lattice deformation $\mathbf{u}(\mathbf{r})$, and optimize it by solving the Euler-Lagrange equation. In decreasing the rotation angle $\theta$, we actually find that $\mathbf{u}(\mathbf{r})$ increases and eventually forms a sharp domain structure with AB and BA regions clearly separated into a triangular pattern.

We then investigate the effect of the domain formation on the electronic spectrum of TBG. The lattice distortion is expected to affect the band structure by modifying the moiré pattern, and also by adding the strain-induced vector potential to graphene’s Dirac electron. Suzuura and Ando (2002); Pereira and Neto (2009); Guinea et al. (2010) Here we calculate the band structure of relaxed TBGs by the tight-binding model and compared it with the non-relaxed case. A significant difference is observed in $\theta<2^{\circ}$, where the domain structure becomes pronounced. The most notable effect is that an energy gap up to 20meV appears at the superlattice subband edge between the lowest and the second subband, while it was hardly found in the non-relaxed model. We also find that the lattice deformation significantly enhances the Fermi velocity, which was strongly suppressed in the non-relaxed model.de Laissardière et al. (2012); Bistritzer and MacDonald (2011c) The associated lattice distortion induces the pseudo magnetic field more than 30T, which alternates in space with the moiré superlattice period.

The paper is organized as follows. In Sec. II. we introduce the lattice geometry of TBG and the description of the moiré pattern. In Sec. III, we develop the effective continuum theory for the lattice relaxation. First we consider the simpler one-dimensional model as an intuitive example, and then we extend the model to two-dimensional TBG. In Sec. IV, we calculate the electronic band structure of relaxed TBGs and discuss the effect of the lattice deformation. A short conclusion is given in Sec. V.

II Geometry of twisted bilayer graphene

Let us consider a TBG lattice as illustrated in Fig. 1. Here we specify the stacking geometry by starting from the AA-stacking bilayer graphene, and rotate the layer 2 with angle $\theta$ with respect to the layer 1. In Fig. 1, we take $\theta=2.65^{\circ}$ as an example. We define the primitive lattice vectors of layer 1 as $\textbf{a}_{1}=a(1,0)$, $\textbf{a}_{2}=a(1/2,\sqrt{3}/2)$, where $a=0.246$ nm is graphene’s lattice constant. The primitive lattice vectors of layer 2 can be obtained by rotating those of layer 1 as $\tilde{\textbf{a}}_{i}=R(\theta)\textbf{a}_{i}\,(i=1,2)$ where $R(\theta)$ is rotation matrix. The reciprocal lattice vectors of layer 1 are given by $\mathbf{a}^{*}_{1}=2\pi/a(1,-1/\sqrt{3})$ and $\mathbf{a}^{*}_{2}=2\pi/a(0,2/\sqrt{3})$, and those of layer 2 by $\tilde{\mathbf{a}}^{*}_{i}=R(\theta)\mathbf{a}^{*}_{i}$ ($i=1,2$).

The Brillouin zone of layer 1 and layer 2 are shown in Fig. 2 by two large black, red dashed hexagons, respectively. In TBG, they are folded into reduced Brillouin zones shown by small blue hexagons. We label the corner points of the folded Brillouin zone by $\bar{K}$ and $\bar{K}^{\prime}$, the midpoint of each side by $\bar{M}$, and the zone center by $\bar{\Gamma}$.

When the rotation angle is small, the mismatch of the lattice periods of two rotated layers gives rise to the long-period moiré beating pattern, of which spatial period is estimated as follows. In the rotation from the AA stacking, an atom on layer 2 originally located at site $\textbf{r}_{0}$ (right above the layer 1’s atom) is moved to the new position $\textbf{r}=R(\theta)\textbf{r}_{0}$. Then we define the interlayer atomic shift $\bm{\delta}(\textbf{r})$ as the in-plane position of an layer 2’s atom at $\mathbf{r}$ measured from its counterpart on layer 1, i.e.,

[TABLE]

When $\bm{\delta}(\mathbf{r})$ coincides with a lattice vector of layer 1, then the position $\mathbf{r}$ (layer 2’s atom) is occupied also by an atom of layer 1, so that the local lattice structure at $\mathbf{r}$ takes AA arrangement as in the origin. Therefore, the primitive lattice vector of the moiré superlattice $\mathbf{L}_{i}^{\rm M}$ is obtained from the condition $\bm{\delta}(\mathbf{L}_{i}^{\rm M})=\mathbf{a}_{i}$, which leads to

[TABLE]

The lattice constant $L_{\rm M}=|\mathbf{L}_{1}^{\rm M}|=|\mathbf{L}_{2}^{\rm M}|$ is given by

[TABLE]

The corresponding moiré reciprocal lattice vectors satisfying $\mathbf{G}^{\rm M}_{i}\cdot\mathbf{L}_{j}^{\rm M}=2\pi\delta_{ij}$ are written as

[TABLE]

where we used $R^{\dagger}=R^{-1}$.

In general TBGs, the lattice structure is not exactly periodic in the atomic level, since the moiré interference pattern is not generally commensurate with the lattice period. However, the superlattice becomes rigorously periodic at some special $\theta$, where vector $m\textbf{a}_{1}+n\textbf{a}_{2}$ meets vector $n\textbf{a}^{\prime}_{1}+m\textbf{a}^{\prime}_{2}$ with certain integers $m$ and $n$.Mele (2010) The exact superlattice period is then given by

[TABLE]

which is $|m-n|$ times as big as the moiré period $L_{\rm M}$. The rotation angle $\theta$ is equal to the angle between two lattice vectors $m\mathbf{a}_{1}+n\mathbf{a}_{2}$ and $n\mathbf{a}_{1}+m\mathbf{a}_{2}$, and it is explicitly given by

[TABLE]

In Table 1, we present $(m,n)$, the rotation angle $\theta$, the size of the moiré unit cell $L_{\rm M}$, and the dimensionless parameter $\eta$ (introduced in Sec. III) for several TBGs considered in this paper.

III Optimized lattice structure

III.1 1-D atomic chain

To describe the lattice relaxation in the continuum theory, we start with a one-dimensional (1D) model Popov et al. (2011) as a simple and intuitive example. The extension to TBG is straightforward as we will see in the next section. Here we consider a 1D moiré superlattice as shown in Fig. 3(a), which is composed of two atomic chains 1 and 2 having slightly different lattice periods, $a=L_{\rm M}/N$ and $a^{\prime}=L_{\rm M}/(N-1)$, respectively, with a large integer $N$. The common period of the whole system is given by $L_{\rm M}$. Inside a supercell, there are $N$ sites and $N-1$ sites in chain 1 and 2, respectively. This model can be viewed as an interacting two-chain version of Frenkel-Kontorova model.

In TBG, the atoms in different layers can be associated by rotation. Likewise, the atoms of chain 1 and 2 are associated by expansion. We can then define the interchain atomic shift $\delta(x)$ as the relative position of the site on chain 2 located at $x$ measured from the position of its counterpart on chain 1, or

[TABLE]

This corresponds to Eq. (1) in TBG. Obviously we have $\delta(nL_{\rm M})=na$ for integer $n$, i.e. the atoms on different chains are vertical aligned at $x=nL_{\rm M}$.

Now we introduce the attractive interaction between the atoms of chain 1 and 2, while allowing the atoms move only in parallel to chain. We expect that the atoms move their positions to reduce the interchain binding energy $U_{B}$. As a result, the system tends to increase the vertically-aligned region, and then a domain structure should be formed as schematically illustrated in Fig. 3(b). At the same time, however, such deformation increases the elastic energy $U_{E}$, so the optimized state can be obtained by minimizing the total energy $U=U_{E}+U_{B}$.

Now we define $u_{1}(x)$ and $u_{2}(x)$ as displacement of atomic positions on layer 1 and layer 2, respectively, measured from the non-relaxed state. The interchain atomic shift in presence of deformation is then

[TABLE]

where $\delta_{0}(x)=(a/L_{\rm M})x$ is that in absence of the deformation, Eq. (7). Following the standard elastic theory, we assume that the elastic energy is expressed as

[TABLE]

where $\kappa$ is elastic constant to characterize the stiffness of the lattice.

If $a$ is very close to $a^{\prime}$, the moiré superperiod $L_{\rm M}$ is much greater than the lattice constant $a$, and then the local lattice structure resembles commensurate chains with the identical lattice period $a$, which are relatively shifted by some specific $\delta$ [Fig. 3(c)]. Let define $V[\delta]$ as the interchain binding energy per unit length of the commensurate chains. Here we assume an attractive interaction described by the sinusoidal function for $V[\delta]$,

[TABLE]

where $V_{0}>0$ and $a^{*}=2\pi/a$. Obviously $V[\delta]$ is periodic with period $a$, because the sliding by the lattice spacing $a$ is equivalent to no sliding. It takes minimum at vertically aligned arrangement, $\delta=na$ ($n$: integer), and maximum at the half shift $\delta=(n+1/2)a$. Now in the incommensurate chains in which $a$ and $a^{\prime}$ are slightly different, the interchain atomic shift $\delta$ is not a constant but slowly varying as a function of $x$. Therefore, the interchain binding energy of the incommensurate chains as a whole is written as

[TABLE]

By using Eqs. (7), (8) and (10), we have $V\left[\delta(x)\right]=-2V_{0}\cos[G_{\rm M}x+a^{*}(u_{1}-u_{2})]$, where $G_{\rm M}=2\pi/L_{\rm M}$ is the reciprocal vector for the moiré superlattice, and we used the relation $a^{*}\delta_{0}(x)=G_{\rm M}x$.

The total energy $U=U_{B}+U_{E}$ is the functional of $u_{1}(x)$ and $u_{2}(x)$. Here we define the coordinates $u_{\pm}=u_{1}\pm u_{2}$, and rewrite $U$ as a functional of $u_{\pm}$. The optimized state to minimize the total energy is obtained by solving the Euler-Lagrange equations,

[TABLE]

In the following, we assume that the lattice deformation keeps the original superlattice period, i.e., $u_{\pm}(x)$ is periodic in $x$ with period $L_{\rm M}$. Then $u_{+}(x)={\rm const.}$ is the only solution of the first equation. To solve the second, we apply the Fourier transformation,

[TABLE]

Eq. (13) is then reduced to

[TABLE]

Eqs.(14), (15) and (16) are a set of self-consistent equations to be solved.

By scaling the displacement $u_{i}(x)$ by $a$, Eq. (16) can be written in a dimensionless form,

[TABLE]

where $\eta$ is a dimensionless parameter defined by

[TABLE]

Roughly speaking, the parameter $\eta$ characterizes how many harmonics are relevant in the displacement $u_{-}(x)$. Since $f_{n}$ is of the order of 1, the condition that $u_{n}$ is comparable to $a$ is given by 2\eta^{2}/(\pi n^{2})\,\lower-1.29167pt\hbox{>}\kern-7.5pt\lower 3.01385pt\hbox{\sim}\,1, or

[TABLE]

When $\eta$ is small such that $(2/\pi)^{1/2}\eta\ll 1$, only the first harmonic term is relevant so $u_{-}(x)$ is well approximated by a single sinusoidal function. This situation occurs in stiff lattice (large $\kappa$), weak interchain interaction (small $V_{0}$) or small moiré period (small $L_{\rm M}$). When $\eta$ is large, on the contrary, the large number of harmonics are significant so that $u_{-}(x)$ becomes a sharp function with respect to the moiré period $L_{\rm M}$. This condition corresponds to soft lattice, strong interaction, and large moiré period.

The self-consistent equations Eqs.(14), (15) and (16) can be solved by numerical iteration with higher harmonics appropriately truncated. In Fig. 4(a), we plot the interlayer atomic shift $\delta(x)$ for the optimized state at some different $\eta$’s. The line of $\eta=0$ represents non-relaxed case $\delta_{0}(x)=ax/L$, and the relative shift from this line represents the displacement $u_{-}(x)$. The actual displacement on each chain is given by $u_{1}=-u_{-}/2$ and $u_{2}=u_{-}/2$. At $\eta=0.3$, $u_{-}$ is small compared to the atomic spacing $a$, and it contains only low frequency Fourier components. In increasing $\eta$, the $u_{-}$ becomes larger and at the same time higher harmonics become more relevant. In $\eta=3.$, we clearly see a step-like structure consisting of two plateau regions of $\delta=0$ and $a$, which are nothing but domains where the atoms are locked to the vertically-aligned positions. Fig. 4(b) presents the corresponding plot for the local binding energy $V[\delta(x)]$. We see that the original cosine function at $\eta=0$ is gradually deformed so as to expand the plateau regions. In $\eta=3$, the system achieves the minimum energy almost everywhere, except for a thin domain boundary in the middle.

Actually, the sharp domain boundary observed in large $\eta$ is well approximated by an analytical soliton solution. If we concentrate on a small region near the domain boundary centered at $x=L_{\rm M}/2$, Eq. (13) is reduced to

[TABLE]

where $x^{\prime}=x-L_{\rm M}/2$, and the term $G_{\rm M}x$ is approximated by $G_{\rm M}(L_{\rm M}/2)=\pi$, assuming the domain boundary is much narrower compared to $L_{\rm M}$. This has an exact solution Popov et al. (2011)

[TABLE]

which is found to nicely agree with the numerically-obtained $u_{-}$ near the boundary. Therefore, the typical width of the domain boundary is characterized by

[TABLE]

Using Eqs. (18) and (22), we have

[TABLE]

so the parameter $\eta$ characterizes the ratio of the domain wall width to the moiré unit cell.

III.2 Twisted bilayer graphene

The above formulation for 1D moiré superlattice can be extend to TBG system in a straightforward manner. Let us consider a TBG with a long-period moiré pattern as illustrated in Fig.  1, and introduce the lattice deformation which is specified by the displacement vector $\textbf{u}^{(l)}(\textbf{r})$ for layer $l=1,2$. The interlayer atomic shift under the deformation is then given by

[TABLE]

where $\bm{\delta}_{0}(\textbf{r})$ is one without lattice deformation, which is defined by Eq. (1). Here we neglect the out-of plane component of the displacement vectors and concentrate on the in-plane motion, in order to describe the domain formation within the simplest framework. The expected effect of the out-of-plane corrugation will be discussed later. The elastic energy is expressed by Suzuura and Ando (2002); San-Jose et al. (2014b)

[TABLE]

where $S_{\rm M}=(\sqrt{3}/2)L_{\rm M}^{2}$ is the area of moiré unit cell, $\lambda\approx 3.5$ eV/$\text{\AA}^{2}$ and $\mu\approx 7.8$ eV/$\text{\AA}^{2}$ are typical values graphene’s $\text{Lam}\acute{\text{e}}$ factor San-Jose et al. (2014a), and $u^{(l)}_{ij}=(\partial_{i}u_{j}^{(l)}+\partial_{j}u_{i}^{(l)})/2$ is strain tensor.

When the moiré superperiod $L_{\rm M}$ is much greater than the lattice constant $a$, the local lattice structure resembles non-rotated bilayer graphene relatively shifted by $\bm{\delta}$ depending on the position [Fig. 5]. We define as $V[\bm{\delta}]$ the interlayer binding energy per area of non-rotated bilayer graphene. In the simplest approximation, it can be written as a cosine function $\bm{\delta}(\textbf{r})$ as

[TABLE]

where $\textbf{a}^{*}_{3}=-\textbf{a}^{*}_{1}-\textbf{a}^{*}_{2}$. The function takes the maximum value $6V_{0}$ at AA stacking ($\bm{\delta}=0$) and the minimum value $-3V_{0}$ at AB and BA stacking. The difference between the binding energies of AA and AB/BA structure is $9V_{0}$ per area, and this amounts to $\Delta\epsilon=9V_{0}S_{G}/4$ per atom where $S_{G}$ is the area of graphene’s unit cell. In the following calculation, we use $\Delta\epsilon=0.0189$ (eV/atom) as a typical value Lebedeva et al. (2011); Popov et al. (2011). The potential profile of $V[\bm{\delta}]$ is presented in Fig. 5.

In TBG, $\bm{\delta}$ is not a constant but slowly varying as a function of the 2D position. Then the inter-layer binding energy of TBG as a whole is written as

[TABLE]

$V[\bm{\delta}]$ is periodic in $\bm{\delta}$ with the lattice period of graphene. By using Eqs. (24) and (26), we have

[TABLE]

where $\textbf{G}^{\rm M}_{3}=-\textbf{G}^{\rm M}_{1}-\textbf{G}^{\rm M}_{2}$ and we used the relation $\textbf{a}^{*}_{j}\cdot\bm{\delta}_{0}(\textbf{r})=\textbf{G}_{j}^{\rm M}\cdot\textbf{r}$.

The relaxed state can be obtained by minimizing total energy $U=U_{E}+U_{B}$ as a functional of $\textbf{u}^{(l)}(\textbf{r})$. We define $\textbf{u}^{\pm}=\textbf{u}^{(2)}\pm\textbf{u}^{(1)}$ and rewrite $U$ as a functional of $\textbf{u}^{\pm}$. In the following we consider only $\textbf{u}^{-}$ as we are interested in relative displacement between atoms on two layers. The Euler-Lagrange equations for $\textbf{u}^{-}$ read

[TABLE]

[TABLE]

We define the Fourier components $\textbf{u}^{-}_{\textbf{q}}$ and $f^{j}_{\textbf{q}}(j=1,2,3)$ as

[TABLE]

where $\textbf{q}=m\textbf{G}^{\rm M}_{1}+n\textbf{G}^{\rm M}_{2}$ are vectors of reciprocal superlattice. Euler-Lagrange equations (29) and (30) are rewritten in a matrix form as

[TABLE]

Eq. (31), (32) and (35) are a set of self-consistent equations. Following Eq. (18) in the 1D model, the number of the relavant harmonics in $\textbf{u}^{-}_{\textbf{q}}$ is roughly characterized by a dimensionless parameter

[TABLE]

In TBG, we have two elastic constants $\lambda$ and $\mu$, and it is ambiguous which should replace the position of $\kappa$ in Eq. (18) for 1D model. Here we adopt the simple sum $\lambda+\mu$ in Eq. (36). Figure 6 plots $\eta$ as a function of the rotation angle $\theta$. The approximation with the lowest harmonics (i.e., six q-points of $\pm\textbf{G}^{\rm M}_{1},\pm\textbf{G}^{\rm M}_{2},\pm\textbf{G}^{\rm M}_{3}$) is valid when $\eta\ll 1$, or $\theta>\sim 2^{\circ}$. The contribution of high frequency harmonics is not negligible when $\eta$ is of the order of 1.

We numerically solve the self-consistent equation for several TBGs by numerical iterations with sufficiently large cut-off in $q$-space. Figure 7 presents an example of calculated result for lattice relaxation in $\theta=1.05^{\circ}$, where the central panel plots the displacement vector $\textbf{u}^{-}(\textbf{r})$ as a function of position, and the left (right) panels show the local atomic structure near AA (AB) stacked point before and after the relaxation. The actual displacement on each layer is given by $\textbf{u}^{(1)}=-\textbf{u}^{-}/2$ and $\textbf{u}^{(2)}=\textbf{u}^{-}/2$. We actually observe that the $\textbf{u}^{-}$ rotates around the center of AA region, and as a result, the AA region is significantly shrunk while AB region is expanded.

Figure  8 summarizes the results for TBGs from $\theta=2.65^{\circ}$ down to $0.547^{\circ}$. Here the panels in (a) show the absolute value of the displacement vector $\textbf{u}^{-}(\textbf{r})$ as a function of position. The distribution of $\textbf{u}^{-}(\textbf{r})$ on two-dimensional place looks all similar among all the cases, where it takes the maximum on a ring-like region near the AA spot. On the other hand, its magnitude strongly depends on $\theta$, where the $|\textbf{u}^{-}|$ is much smaller than the atomic scale $a$ in $\eta\ll 1$, while it eventually becomes comparable when $\eta$ is of the order of 1. Figure 8(b) presents the corresponding plots for the local binding energy $V[\bm{\delta}(\textbf{r})]$. When $\textbf{u}^{-}(\textbf{r})$ is much smaller than $a$, as in $\theta=2.65^{\circ}$, the potential profile is approximately given by $V[\bm{\delta}_{0}(\textbf{r})]$, which is essentially a sum of three plain waves. In decreasing $\theta$, the spots of $AA$-regions shrink and AB and BA regions eventually dominate. The result looks especially dramatic in small angles less than $1^{\circ}$, where the relaxed lattices clearly exhibits a triangular domain pattern of AB and BA regions. Similar to Eq. (22) for the 1D model, the characteristic width of the domain boundary is given by

[TABLE]

Indeed it roughly agrees with the typical scale of the AB/BA domain wall in Fig. 8(b). It is also consistent with the experimental observation of the shear boundary, which estimates the averaged width about 6 nm.Alden et al. (2013)

IV Band structure

To calculate the energy band structures in the presence of the lattice strain, we use the tight-binding method. The Hamiltonian is written as

[TABLE]

where $\textbf{R}_{i}$ is the atomic coordinate, $|\textbf{R}_{i}\rangle$ is the wavefunction at site $i$, and $t(\textbf{R}_{i}-\textbf{R}_{j})$ is the transfer integral between atom $i$ and $j$. We adopt the Slater-Koster type formula for the transfer integral, Slater and Koster (1954)

[TABLE]

[TABLE]

where $\textbf{d}=\textbf{R}_{i}-\textbf{R}_{j}$ is the distance between two atoms. $\textbf{e}_{z}$ is the unit vector on $z$ axis. $V_{pp\pi}^{0}\approx-2.7$eV is the transfer integrals between nearest-neighbor atoms of monolayer graphene which are located at distance $a_{0}=a/\sqrt{3}\approx 0.142$nm. $V_{pp\sigma}^{0}\approx 0.48$eV is the transfer integral between two nearest-vertically aligned atoms. $d_{0}\approx 0.334$nm is the interlayer spacing. The decay length $r_{0}$ of transfer integral is chosen at $0.184a$ so that the next nearest intralayer coupling becomes $0.1V_{pp\pi}^{0}$. Uryu (2004); Trambly de Laissardière et al. (2010) At $d>\sqrt{3}a$, the transfer integral is very small and negligible.

Using the optimized structure obtained in the last section, we specify the lattice position of each single atom in the relaxed TBG, construct the tight-binding Hamiltonian, and calculate the energy bands. Figure 9 compares the electronic band structure of relaxed (black solid lines) and non-relaxed (red dashed lines) TBGs at several rotation angles. The horizontal axis are labeled by the symmetric points of the Brillouin zone for the moiré superlattice [Fig. 2], and it scales in proportion to $2\pi/L_{\rm M}$. At $\theta=2.65^{\circ}$, we only see a minor difference in accordance with the small change in the lattice structure in Fig. 8(a). A significant deviation is observed below $2^{\circ}$. The most notable change from the non-relaxed case is that a band gap opens between the lowest subband near the Dirac point and the first excited subband both in the electron side and the hole side. Fig. 10 (a) shows the size of gap versus rotation angle $\theta$. The gap is observed in TBGs of 1^{\circ}\,\lower-1.29167pt\hbox{<}\kern-7.5pt\lower 3.01385pt\hbox{\sim}\,\theta\,\lower-1.29167pt\hbox{<}\kern-7.5pt\lower 3.01385pt\hbox{\sim}\,1.5^{\circ}, and the maximum energy width is about 18meV.

The lattice strain also strongly modifies the band velocity. Figure 10 (b) plots the band velocity at the Dirac point as a function of $\theta$ for relaxed and non-relaxed cases. In both cases, the central band at the Dirac point is gradually flattened in decreasing $\theta$, and the Fermi velocity vanishes at a certain angle de Laissardière et al. (2012); Bistritzer and MacDonald (2011c). We find that the band flattening is a little slower in the relaxed case, i.e at the same angle, the band width is larger in the relaxed TBG than in non-relaxed counterpart, so that the critical angle for the vanishing velocity is shifted to the lower rotation angle in the relaxed TBG.

The lattice relaxation affects the electronic structure in two different ways, by a change of interlayer Hamiltonian associated with the modified moiré pattern, and also by a change of the intralayer Hamiltonian through distortion of the lattice. The latter is known to be described by the pseudo magnetic field in the effective mass Dirac Hamiltonian.Suzuura and Ando (2002); Pereira and Neto (2009); Guinea et al. (2010) The vector potential for the pseudo-magnetic field on layer $l(=1,2)$ is given by Suzuura and Ando (2002); Pereira and Neto (2009); Guinea et al. (2010)

[TABLE]

where $\gamma_{0}=t(a_{0})$ is the nearest neighbor transfer energy of intrinsic graphene, $v=(\sqrt{3}/2)a\gamma_{0}$ is the band velocity of the Dirac cone, and

[TABLE]

In the present model Eq.  (39), we have $\beta=a_{0}/r_{0}\approx 3.14$. The pseudo magnetic field is given by $B^{(l)}_{\rm eff}=[\nabla\times\mathbf{A}^{(l)}]_{z}$.

Figure 8(c) shows the distribution of or the pseudo magnetic field on the layer 1 for several TBG’s. The field direction is opposite between layer 1 and 2, because $\textbf{u}^{(1)}=-\textbf{u}^{(2)}$. We observe a triangular pattern with positive and negative field domains, which are centered at the AB and BA stacking regions, respectively. The field amplitude is huge, but it does not necessarily results in a strong effect on the electronic structure, since it is rapidly oscillating in space with nano-meter scale. The pseudo-magnetic field enters in the Hamiltonian as a form of $evA$ with the pseudo vector potential $A$. When the field spatially modulates with the wave length $L_{\rm M}$, the associated matrix element opens a band gap at the energy $E\sim\hbar v/L_{\rm M}$ measured from the Dirac point. Therefore, the effect of the pseudo field significantly affects the band structure when $evA>\sim\hbar v/L_{\rm M}$, while otherwise it is just perturbative. Now the scale of $evA$ is roughly estimated as

[TABLE]

where the strain tensor $u_{ij}$ is estimated about $u/L_{\rm M}=2(u/a)\sin(\theta/2)$, considering the displacement field $\mathbf{u}$ is modulating with the moire wavelength $\sim 1/L_{\rm M}$. The typical scale of $u/a$ can be read from Fig. 8(a). For $\theta=2.65^{\circ}$, for example, $u/a\sim$ 0.03 gives $evA\sim$ 10meV, and it is much smaller than $\hbar v/L_{\rm M}\sim$ 180meV. So the effect of $evA$ is perturbative, and this is consistent with a small change in the band structure observed in Fig.  9 (a). For $\theta=1.05^{\circ}$, on the other hand, $evA\sim$ 30meV is comparable with $\hbar v/L_{\rm M}\sim$ 50meV, so the pseudo field plays a significant role in the modification of the low-energy bands.

V Conclusion

We have developped the effective theory to calculate the spontaneous relaxation in TBG, and studied the atomic and electronic structures. In rotation angle larger than $2^{\circ}$, the lattice is hardly deformed and so the effect on the electronic structure is minimal, while in smaller rotation angle below $2^{\circ}$, the lattice is significantly modified to form AB/BA triangular domain structure. The electronic band structure is then strongly modified where a band gap up to 20 meV opens above and below the lowest band. The lattice deformation also significantly relaxes the band flattening observed in non-relaxed case, and it lowers the critical angle at which the Fermi velocity vanishes.

Actually a recent experiment observed an insulating gap about 50 meV at the superlattice subband edges in TBG with $\theta\approx 1.8^{\circ}$. Cao et al. (2016) This seems qualitatively consistent with the present result, although 1.8∘ is out of the gap-opening range in our model calculation, and also 50 meV is a bit too large compared to the typical gap width obtained here. As we see in the present work, however, the lattice relxation and the electronic structure sensitively depend on the parameter $\eta$, and it might be possible that the real system could have a greater interlayer interaction, allowing a greater gap and a wider range in the rotation angle for gap opening. It is also conceivable that the gap could be enhanced when the Fermi energy is right at the superlattice gap position, while the doping effect is not considered in the present study. We leave the further quantitative arguments for a future problem.

The present model takes account of only the in-plane components of the lattice distortion, as it is aimed to describe the domain formation within the simplest theoretical framework. Inclusion of the out-of-plane motion is known to give rise to a corrugation in the perpendicular direction, Uchida et al. (2014); van Wijk et al. (2015); Dai et al. (2016); Jain et al. (2016) where the interlayer spacing modulates by 10%. In a small angle TBG less than 2∘, in particular, the detailed computational study that the the interlayer spacing is largest only near AA spot while it is almost flat otherwise. van Wijk et al. (2015); Dai et al. (2016) The corrugation is small even at the AB/BA domain boundary, presumably because it is a shear boundary with no tensile strain, and also the optimized interlayer spacing does not strongly depend on the stacking structure around there. Therefore, we expect that the corrugation effect on the electronic structure exclusively comes from AA spots, where the interlayer distance change should reduce the interlayer electronic coupling by a few 10%. Since the system is dominated by AB/BA regions, the change of the electronic band structure is expected to be minor compared to the change caused by AB/BA domain formation itself.

ACKNOWLEDGMENTS

The authors thank Pilkyung Moon for helpful discussions. This work was supported by JSPS KAKENHI Grants No. JP25107001, No. JP25107005 and JP15K21722.