Pseudomagnetic fields in fully relaxed twisted bilayer and trilayer graphene

Twisted bilayer graphene shows a very complex electronic phase diagram for small twist angles arising from the existence of narrow electronic bands for certain angles referred to as 'the magic' angles [1–7]. The width and shape of the electronic bands near these magic angles are highly sensitive to many perturbations, such as strains [8, 9], interactions [10], or the effect of the substrate [11–15]. These effects typically increase the width of the narrow middle bands, (whose width varies between 4 and 40 meV), and can also give rise to small gaps between the valence and conduction bands.

The interaction between the two twisted layers in tBLG leads to lattice relaxation [16–25], which shrinks the highly unstable AA regions while enlarging the more energetically favorable AB regions effectively reducing the twist angle in the AB regions while increasing the local twist angle in the AA-stacking zones. Lattice relaxation has been shown to play a pivotal role in the piezoelectric, ferroelectric and electronic properties of twisted TMDs [26–39] as well as in the charge density wave phases of twisted NbSe₂ [40, 41]. In graphitic systems [42], atomic relaxation has been shown to determine the out-of-plane polarization of twin boundaries [43] while in supermoiré systems such as the non-symmetric twisted trilayer graphene, lattice relaxation leads to a clear separation between the flat band and the highly dispersive Dirac cone [44–60]. For sufficiently small twist angles, $\theta \lesssim 0.5^\circ$, a domain pattern, made of AB and BA regions separated by domain walls emerges, leading to a conducting network of topologically protected 1D channels [61–66]. Near the first magic angle, where the Fermi velocity vanishes [67, 68], $\theta \approx 1^\circ$, the in-plane atomic displacements lead to an effective magnetic field [16, 69, 70] with opposing sign in the two layers and two valleys which heavily renormalizes the intervalley hopping.

In this manuscript, we analytically obtain the lattice relaxation-induced pseudomagnetic field [71, 72] profile for twist angles near the first magic angle as well as the scalar potential derived from the in-plane and out-of-plane lattice relaxation. Then, we added to our continuum tBLG and tTLG Hamiltonian an additional term accounting for the energy difference between the interlayer p_z orbital dimer and non-dimer pairs [73]. We find that the in-plane deformations not only shift the position of the tBLG magic angle towards lower twist angles but also form an entire range of twist angles $\theta\in[0.8-0.85]^{\circ}$ where both the lowest conduction and highest valence bands are very flat [74]. Similarly, in a tTLG made of three AAA-stacked layers with the middle layer twisted at an angle θ relative to the top and bottom graphene membranes (see figure 2(b)), the collapse of the Fermi velocity towards lower angles and the formation of a magic angle plateau was found even more pronounced. Finally, lattice relaxation hugely altered the electron–hole asymmetry both in twisted bilayer and trilayer graphene while it preserved the topological protection of the Dirac cones.

Lattice relaxation originates from the competition between van der Waals forces which tend to enlarge the $AB/BA$ regions and the graphene in-plane stiffness which counteracts such tendency. We start our analysis considering the superlattice geometry presented in figure 2. Firstly we orient the monolayer graphene such that their lattice vectors are taken as $\mathbf{a}_{1} = a(1,0),\mathbf{a}_{2} = a (\frac{1}{2},\frac{\sqrt{3}}{2})$ where a = 2.46 Å, which prescribes the following two reciprocal lattice vectors of the monolayer graphene membrane, $\mathbf{g}_{1} = \frac{4\pi}{\sqrt{3}a} (0,-1),\mathbf{g}_{2} = \frac{4\pi}{\sqrt{3}a} (\frac{\sqrt{3}}{2},-\frac{1}{2})$. Considering an AA-stacked bilayer with lattice vectors $\mathbf{a}_{1},\mathbf{a}_{2}$ we first apply a $\theta/2$ counter-clockwise rotation of the top layer around the center of the hexagonal lattice. Then we rotate the bottom layer by an angle $-\theta/2$ relative to the same rotation axis, shifting the two originally overlayed carbon atoms located at r₀ relative to the rotation center by an amount $\Delta\mathbf{r} = \mathbf{r}_{t}-\mathbf{r}_{b} = \theta \hat{z}\times \mathbf{r}_{0}$. This new geometrical arrangement results in an approximately periodic pattern (see figure 2(a)) characterized by a new periodic Moiré lengthscale defined as,

$$\begin{align} L_{M} = \frac{a}{2\sin{\left(\theta/2\right)}}. \end{align} \tag{ 1 }$$

Given our previously defined monolayer orientation, the Moiré reciprocal superlattice vectors are obtained from the difference in wavevector between the rotated graphene reciprocal lattice vectors,

$$\begin{align} \mathbf{G}_{i} = R\left(-\theta/2\right)\mathbf{g}_{i}-R\left(\theta/2\right)\mathbf{g}_{i} = \mathbf{g}^{b}_{i}-\mathbf{g}^{t}_{i} = -\theta\hat{\mathbf{z}}\times \mathbf{g}_{i}. \end{align} \tag{ 2 }$$

Therefore, $\mathbf{G}_{1} = \frac{-4\pi}{\sqrt{3}L_{M}} (1,0)$ and $\mathbf{G}_{2} = \frac{-4\pi}{\sqrt{3}L_{M}} (\frac{1}{2},\frac{\sqrt{3}}{2})$. We model the adhesion van der Waals potential $W(\vec{r})$ between the displaced layers by Fourier expanding the potential in Moiré reciprocal lattice vectors as presented in [75, 76]. Considering the lowest and more dominant harmonics, we obtain the following approximate expression for the adhesion potential which satisfies all the symmetries of the Moiré supercell [16, 17],

$$\begin{align} W\left(\vec{r}\right)\approx V_{\mathrm{vdW}} \sum_{i = 1, 2, 3} \cos \left[ \vec{g}_i \left( \vec{r}_t - \vec{r}_b \right) \right], \end{align} \tag{ 3 }$$

where the three $\vec{g}_i$ are three trigonally symmetric reciprocal lattice vectors of each graphene layer and $\vec{r}_t,\vec{r}_b$ label the spatial coordinates within the xy-plane of a point in the top and bottom layer respectively. The amplitude $V_{\mathrm{vdW}}$ is set to $V_{\mathrm{vdW}} = 0.9$ meV Å⁻², following the electron-phason study presented in [23]. Due to the interplay between the adhesion and the elastic forces, the top and bottom graphene membrane experience an atomic displacement field $\vec{u}^{t}$ and $\vec{u}^{b}$ respectively. To balance the force that the top layer does onto the bottom one and vice-versa, we assume that lattice relaxation induces at each point of the two graphene membranes a perfectly antisymmetric displacement field such that $\vec{u}_{t} = -\vec{u}_{b} = \vec{u}$. Therefore, each atomic location of the top and bottom layer after twisting and relaxing becomes

$$\begin{align} \vec{r}_t & = R_{\theta/2} \left( \vec{r} \right)+\vec{u} \nonumber \\ \vec{r}_b & = R_{- \theta/2} \left( \vec{r} \right)-\vec{u}, \end{align} \tag{ 4 }$$

where $R_{\theta/2}$ is the rotation operator by a twist angle of $\theta/2$. Equation (3) can be rewritten in terms of the supercell reciprocal lattice vectors as,

$$\begin{align} \vec{g}_i \left[ R_{\theta/2} \left( \vec{r} \right) - R_{- \theta/2} \left( \vec{r} \right) \right] \approx \vec{G}_i \vec{r}, \end{align} \tag{ 5 }$$

which leads to

$$\begin{align} W \left( \vec{r} \right) &\approx V_{\mathrm{vdW}} \sum_{i = 1,2,3} \cos \left( \vec{G}_i \vec{r} + 2 \vec{g}_i \vec{u} \right). \end{align} \tag{ 6 }$$

Assuming that $V_{\mathrm{vdW}}\gt 0$, the function $W(\vec{r})$ has a maximum ($W_{\mathrm{max}}] = 3V_{\mathrm{vdW}}$) at the center of the real space unit cell (AA region), and two minima, $W_{\mathrm{min}} = -(3V_{\mathrm{vdW}})/2$ located at the corners of the Moiré unit cell [77] (AB and BA regions). Note that in equation (6), the reciprocal lattice vectors $\vec{G}_i$ and $\vec{g}_i$ are orthogonal with magnitude $|\vec{G}_i| = (4\pi)/(\sqrt{3}L_{M})$ and $|\vec{g}_i| = (4\pi )/(\sqrt{3}a)$ respectively, where L_M and a are the lattice constants of the Moiré supercell and of each isolated monolayer. The van der Waals forces on the atoms experienced at the top and bottom layers can be obtained by taking the gradient of the potential in equation (6) with respect to $\pm\vec{u}$. Assuming that $\vec{G}_{i}\vec{r}\gt\gt 2\vec{g}_{i}\vec{u}$, the force of the top and bottom layer become,

$$\begin{align} \vec{F}_{t,b}\left(\vec{r}\right)& = \mp V_{\mathrm{vdW}}\sum \vec{g}_i \sin\left(\vec{G}_i\vec{r}\right). \end{align} \tag{ 7 }$$

These forces define a transverse vector field normal to the van der Waals potential which induces a static atomic displacement field towards their equilibrium position. By balancing the force in equation (7) with that of a frozen transverse acoustic phonon whose elastic energy is only given by the Lamé coefficient µ = 9.25 eV Å⁻² [78], we obtain an equilibrium condition for the atomic displacement at each point of the two graphene membranes,

$$\begin{align} u_{t} = -\vec{u}_{b} &\approx \frac{V_{\mathrm{vdW}}}{\mu G^2} \sum_{i = 1,2,3} \vec{g}_i \sin \left( \vec{G}_i \vec{r} \right). \end{align} \tag{ 8 }$$

From the above atomic displacement fields, we predict a vanishing deformation potential [79] ($V = \Upsilon(u_{xx}+u_{yy})$ where $\Upsilon = 4.5$ eV [80] is the work function of graphene) and a twist-angle independent pseudomagnetic field [81, 82] given by $\vec{B} = \vec{\triangledown}\times\vec{A}$ where $\vec{A}(\vec{r}) = \pm\frac{3}{4}\frac{\gamma_{0}|\beta|}{v_F}(u_{xx}-u_{yy},-u_{xy}-u_{yx})$,

$$\begin{align} B \left(\vec{r}\right) = \mp \sqrt{3} \pi \frac{\gamma_{0} \beta}{v_F} \frac{V_{\mathrm{vdW}}}{\mu a} \sum_{i = 1,2,3} \sin \left( \vec{G}_i \vec{r} \right). \end{align} \tag{ 9 }$$

In the above expression, $\gamma_{0}\approx 2.7$ eV [83] is magnitude of the hopping parameter between nearest neighboring carbon-atom p_z orbitals and β is the Grünesein parameter [84] $\beta\equiv-\frac{\partial {\mathrm{ln}}(\gamma_{0})}{\partial {\mathrm{ln}}(a_{c})}\approx 3$, with $a_{c} = a/\sqrt{3}\approx 1.42$ Å being the distance between nearest neighboring atoms. The parameter $v_F = \frac{\gamma_{0}\sqrt{3}a}{2}\approx 5.59$ eV Å is the monolayer graphene Fermi velocity. From equation (9), the characteristic magnetic length of electrons located at the magnetic hot spots in the AB and BA is obtained as,

$$\begin{align} \ell_M = a\sqrt{\frac{1 }{2 \pi} \frac{1}{| \beta |}\frac{\mu}{V_{\mathrm{vdW}}}}\approx 5.6~\textrm{nm} \end{align} \tag{ 10 }$$

for $V_{\mathrm{vdW}}/\mu \approx 10^{-4}$. This pseudomagnetic field averages to zero over the unit cell, and reaches $\approx \pm 49$ T at the AB and BA corners (see inset in figure 1). To improve this simple formulation of the in-plane relaxation in tBLG, we plug the obtained $(u_{x},u_{y})$ from equation (8) into the adhesion term in equation (6) and solve the Euler–Lagrange equations via the Fourier transform method (see appendix A) to recalculate a new atomic displacement field $(\delta U_{x},\delta U_{y})$. The solutions of the resulting strain field $(\delta U_{x},\delta U_{y})$ have the form

$$\begin{align} \delta U_{x}\left(\vec{r}\right) & = \frac{\sum^{3,\infty\infty,\infty}_{i = 1,l = 0,m = 0,n = 0}\left(\left(2\mu+\lambda\right)G_{l,m,n,iy}^{2}+\mu G_{l,m,n,ix}^{2}\right)A^{\pm l,\pm m,\pm n,i}_{x}\sin{\left(\left(\pm l\vec{G}_{1}\pm m\vec{G}_{2}\pm n\vec{G}_{3}+\vec{G}_{i}\right)\vec{r}\right)}}{\mu\left(2\mu+\lambda\right)G_{l,m,n,i}^{4}} \notag\\ & \quad -\frac{\sum^{3,\infty\infty,\infty}_{i^{^{\prime}} = 1,l^{^{\prime}} = 0,m^{^{\prime}} = 0,n^{^{\prime}} = 0}\left(\mu+\lambda\right)G_{l^{^{\prime}},m^{^{\prime}},n^{^{\prime}},i^{^{\prime}}y}G_{l^{^{\prime}},m^{^{\prime}},n^{^{\prime}},i^{^{\prime}}x}A^{\pm l^{^{\prime}},\pm m^{^{\prime}},\pm n^{^{\prime}},i^{^{\prime}}}_{y}\sin{\left(\pm l^{^{\prime}}\vec{G}_{1}\vec{r}\pm m^{^{\prime}}\vec{G}_{2}\vec{r}\pm n^{^{\prime}}\vec{G}_{3}\vec{r}+\vec{G}_{i^{^{\prime}}}\vec{r}\right)}}{\mu\left(2\mu+\lambda\right)G_{l^{^{\prime}},m^{^{\prime}},n^{^{\prime}},i^{^{\prime}}}^{4}}, \nonumber \end{align} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. (Top) Parameters used for the calculation of the electronic bandstructure in twisted bilayer and trilayer graphene. (Bottom) Parameters used for the calculation of the elastic deformations in the two membranes.

ωAA	79.7 meV	ωAB	97.5 meV
γ0	2.6 eV	γ2	−17 meV
Δ	25 meV	vF	5.59 eV Å
µ	9.57 eV Å−2	λ	3.25 eV Å−2
ϒ	4.5 eV	$V_{\mathrm{VdW}}$	0.9 meV Å−2
β	3	acc	1.42 Å
a	2.46 Å	d	3.35 Å
$\Delta h$	0.24 Å

with $\delta U_{y}(\vec{r})$ obtained by replacing G_x by G_y and G_y by G_x into equation (11),

$$\begin{align} & \delta U_{y}\left(\vec{r}\right)\notag\\ & \quad = \delta U_{x}\left(\vec{r}\right)\left[\left(G_{x},G_{y},A_{x},A_{y}\right)\longrightarrow \left(G_{y},G_{x},A_{y},A_{x}\right)\right]. \end{align} \tag{ 12 }$$

In equation (11), the terms $A^{\pm l,\pm m,\pm n,i}_{x}$ and $A^{\pm l,\pm m,\pm n,i}_{y}$ are defined as in equations (A14) and (A15), of appendix A, while µ = 9.57 eV Å⁻²,λ = 3.24 eV Å⁻² [78] are the graphene's Lamé and shear modulus coefficients respectively. Finally, the terms $G_{l,m,n,i}$ are defined as $G_{l,m,n,i} = l\vec{G}_{1}+m\vec{G}_{2}+ n\vec{G}_{3}+\vec{G}_{i}$ where $\vec{G}_{1} = \frac{4\pi}{\sqrt{3}L_{M}}(1,0)$, $\vec{G}_{2} = R(2\pi/3)\vec{G}_{1}$ and $\vec{G}_{3} = R(2\pi/3)\vec{G}_{2}$ . As shown in the comparison presented in figure 6 between our formalism and the numerical method developed in [16], the description provided in equations (11) and (12), correctly catches both the pseudomagnetic and pseudoscalar potential down to angles $\theta\gt 0.5^{\circ}$ .

While equations (11) and (12) accurately model the in-plane relaxation of two coupled membranes, a realistic model of atomic relaxation in tBLG requires taking into account the out-of-plane deformation [20, 85–88] of its constituent membranes. Given its limited coupling with the planar strain field (see figure 1), we will account for the effect of out-of-plane atomic displacements by simply plugging into the z-component of the strain field tensor a height profile of the form

$$\begin{align} h\left(\vec{r}\right) = \left(d+\frac{3h_{0}}{2}\right)+h_{0}\sum^{i = 3}_{i = 1}\cos{\left(\vec{G}_{i}\boldsymbol{\cdot}\vec{r}\right)}, \end{align} \tag{ 13 }$$

where $h_{0} = \frac{2\Delta h}{9}\sim 0.06$ Å [19] is a quarter of the difference in the interlayer separation between the AA and AB-stacked regions, and d = 3.35 Å is the interlayer distance in Bernal bilayer graphene. As shown in appendix A, this additional term adds to the solutions of the Euler–Lagrange equations $(\delta U_{x},\delta U_{y})$ a radially symmetric term $(\delta U^{h}_{x},\delta U^{h}_{x})\propto|G|h^{2}_{0}\sin{(\vec{G}_{i}+\vec{G}_{j})}$ which adds ∼0.1–0.2 meVs to the pseudoscalar potential $V(\vec{r})\propto (u_{xx}+u_{yy})$, becoming the main contributor to the pseudoscalar potential only at angles $\theta\gt 1^{\circ}$ (see figure 1 and figure 3).

Zoom In Zoom Out Reset image size

To compute the electronic bandstructure of tBLG and tTLG with atomic relaxation, we will tune the commonly used Bistritzer and MacDonald (BM) model [67] to accommodate all the effects derived from atomic relaxation. In its simplest form, the BM model consists of two rotated Dirac Hamiltonians connected by a unique interlayer hopping parameter [89] ω = 110 meV which couples any two Dirac cones in opposite layer separated an amount $\Delta K = k_{\theta} = 2K\sin(\theta/2)$ in momentum space [90]. Written in the four-basis $\Psi^{t}_{A},\Psi^{t}_{B},\Psi^{b}_{A},\Psi^{b}_{B}$, the BM Hamiltonian has the form,

$$\begin{align} \mathcal{H}_{\mathrm{BM}} & = \begin{pmatrix} R_{\theta/2}\left[-v_{F}i\vec{\sigma}\vec{\triangledown}\right]& \mathcal{T}\\ \mathcal{T}^{\dagger} & R_{-\theta/2}\left[-v_{F}i\vec{\sigma}\vec{\triangledown}\right] \\ \nonumber \end{pmatrix},\notag\\ \mathcal{T} & = \sum^{j = 2}_{j = 0}\omega\begin{pmatrix} 1& e^{ij\phi}\\ e^{-ij\phi}&1 \\ \end{pmatrix}e^{-i\Delta\vec{\mathbf{K}}_{j}\vec{r}}, \end{align} \tag{ 14 }$$

where $\phi = -2\pi/3$, $R_{\theta/2}$ is a rotation operator twisting the Dirac Hamiltonian an angle $\theta/2$, and $\Delta \vec{\mathbf{K}}_{0} = k_{\theta}(0,-1),\Delta \vec{\mathbf{K}}_{1} = k_{\theta}(\frac{\sqrt{3}}{2},1/2),$ and $\Delta \vec{\mathbf{K}}_{2} = k_{\theta}(-\frac{\sqrt{3}}{2},1/2)$ are the momentum differences between neighboring Dirac cones. Due to the different interlayer distance between the AA and AB-stacking patches, the interlayer hopping parameter in the AA regions (ω_AA ) and the AB or BA regions (ω_AB ) will be slightly different. To account for this difference, we will use the parameterization presented in [91] for tBLG, which obtained ω_AA and ω_AB by Fourier transforming the Slater–Koster interlayer tunneling amplitude of an AA and AB stacked-bilayer around the K-point. Therefore, the interlayer tunneling becomes

$$\begin{align} \mathcal{T} = \sum^{j = 2}_{j = 0}\begin{pmatrix} \omega_{AA}& \omega_{AB}e^{ij\phi}\\ \omega_{AB}e^{-ij\phi}& \omega_{AA} \\ \end{pmatrix}e^{-i\Delta\vec{\mathbf{K}}_{j}\vec{r}}. \end{align} \tag{ 15 }$$

On top of the interlayer hopping renormalization due to the difference in interlayer distance, it is necessary to consider to effect of in-plane relaxation in ω_AA and ω_AB . To account for such an effect, we start from the Fourier expansion of the interlayer tunneling term in equation (14), and, as in [71] we replace $\vec{r}_{t/b}$ with $\vec{r}_{t/b}\longrightarrow \vec{r}_{t/b}+\vec{u}_{t/b}$. Plugging the expression derived in equation (8) and Taylor expanding the complex exponential in equation (14), we can obtain a first-order correction to the interlayer tunneling terms ω_AA and ω_AB (see appendix B for a more detailed discussion) whose order of magnitude can be estimated as,

$$\begin{align} \frac{\delta \omega_{AA/AB}}{\omega_{AA/AB}} \approx i\frac{2V_{vdW}}{\sqrt{3}\mu}\frac{L_{M}}{a}\approx i\frac{0.042}{\theta}. \end{align} \tag{ 16 }$$

Since $\delta w_{AA}\ll\omega_{AA}$ and $\delta w_{AB}\ll\omega_{AB}$, we will neglect the effect of in-plane relaxation in the interlayer hopping amplitudes.

To better capture the effect of atomic relaxation within the continuum model, we need to account for the energy difference in each carbon $A/B$ site due to its surrounding atomic registry. In an AB bilayer graphene structure, it is well known that the $A^{b}-B^{t}$ dimer pairs are at higher energy compared to their non-dimer $A^{t}-B^{b}$ counterparts. Therefore, in systems such as tBLG, this energy difference labeled by the parameter $\Delta = 25$ meV [92] is expected to create a difference in chemical potential between the AB/BA stacked regions and the more energetically unfavorable AA regions. As shown in [73], this additional term can be introduced in the tBLG continuum model by projecting the crystal potential of one layer into the other one, leading to an additional term $\mathcal{E}^{t}$ for its top layer,

$$\begin{align} \mathcal{E}^{t} = \begin{pmatrix} \mathcal{E}^{tA} & 0\\ 0 & \mathcal{E}^{tB} \end{pmatrix} = \frac{2\Delta}{9}\sum^{i = 2}_{i = 0}\begin{pmatrix} \cos{\left(\mathbf{G}_{i}\mathbf{r}\right)}\left(1+\cos{\left(\mathbf{g_{i}}{\tau_\mathbf{b}}\right)}\right)-\sin{\left(\mathbf{G}_{i}\mathbf{r}\right)}\sin{\left(\mathbf{g_{i}}{\tau_\mathbf{b}}\right)} & 0\\ 0 & \cos{\left(\mathbf{G}_{i}\mathbf{r}\right)}\left(1+\cos{\left(\mathbf{g_{i}}{\tau_\mathbf{b}}\right)}\right)+\sin{\left(\mathbf{G}_{i}\mathbf{r}\right)}\sin{\left(\mathbf{g_{i}}{\tau_\mathbf{b}}\right)} \\ \end{pmatrix}, \end{align} \tag{ 17 }$$

where ${\tau}_{B} = (0,a/\sqrt{3})$ and whose bottom layer counterpart has the form,

$$\begin{align} \mathcal{E}^{b} = \begin{pmatrix} \mathcal{E}^{tB} & 0\\ 0 & \mathcal{E}^{tA} \end{pmatrix}. \end{align} \tag{ 18 }$$

Upon adding to the tBLG Hamiltonian the relaxation-induced gauge field ($\mathcal{A}_{x}(\vec{r}),\mathcal{A}_{y}(\vec{r})$), the dimerization energy shifts $\mathcal{E}(\vec{r})$ as well as the two strain-induced pseudoscalar pseudopotentials ($V_{z}(\vec{r}) = \Upsilon[\partial_{x}\delta U^{h}_{x}(\vec{r})+\partial_{y}\delta U^{h}_{y}(\vec{r})]$ and $V_{||}(\vec{r}) = \Upsilon[\partial_{x}\delta U_{x}(\vec{r})+\partial_{y}\delta U_{y}(\vec{r}))$]) the fully relaxed tBLG graphene Hamiltonian written on the four-state basis, $\Psi^{t}_{A},\Psi^{t}_{B},\Psi^{b}_{A},\Psi^{b}_{B}$, becomes

$$\begin{align} H_{tBLG} = \begin{pmatrix} \vec{\sigma}_{\theta/2}\left(-v_{F}i\vec{\triangledown}+\mathcal{A}\left(\vec{r}\right)\right)+\mathcal{E}^{t}\left(\vec{r}\right)+V_{||}\left(\vec{r}\right)+V_{z}\left(\vec{r}\right) & \mathcal{T}\left(\vec{r}\right)\\ \mathcal{T}^{\dagger}\left(\vec{r}\right) & \vec{\sigma}_{-\theta/2}\left(-v_{F}i\vec{\triangledown}-\mathcal{A}\left(\vec{r}\right)\right)+\mathcal{E}^{b}\left(\vec{r}\right)-V_{||}\left(\vec{r}\right)-V_{z}\left(\vec{r}\right) \end{pmatrix}. \end{align} \tag{ 19 }$$

Note that in equation (19), the pseudogauge field $\mathcal{A}(\vec{r})$ and the two pseudoscalar potentials V_z and $V_{||}$, have different signs in each layer due to the opposite direction of the atomic displacement field in each graphene membrane. As the twist angle increases, the unit cell becomes smaller suppressing the expansion of the AB and BA regions. This permits us to safely approximate $\vec{u}$ by the expression in equation (8) and B by the formula in equation (9). Due to the large momentum mismatch between nearest-neighboring Dirac cones in the high-twist angle limit, the energy difference between $\mathbf{K}^{\mathrm{top}}$ and $\mathbf{K}^{\mathrm{bottom}}$ becomes larger. As pointed out in [90], this energetic separation permits to safely truncate the BM Hamiltonian into an $8\times 8$ matrix which only considers the resonant tunneling between a top-layer Dirac cone and its three nearest-neighboring K-point states of the bottom layer. Ignoring both the pseudoscalar and the dimerization potential, we account for atomic relaxation by adding a pseudogauge field coupling only the bottom layer K-point states (see the inset in figure 5(a)). The resulting Hamiltonian written in the eight-state basis made of 4 Dirac cones ($\Psi^{A_{1}},\Psi^{B_{1}},\Psi^{B_{2}},\Psi^{B_{3}}$) × 2 sublattice sites ($A,B$) per cone reads as,

$$\begin{align} & H_{8\times 8} \notag\\ & \quad = \begin{pmatrix} v_{F}\mathbf{k}\vec{\sigma} & \mathcal{T}^{1} & \mathcal{T}^{2} & \mathcal{T}^{3} \\ \mathcal{T}^{1\dagger} & v_{F}\left(\mathbf{k}-\mathbf{q}_{1}\right)\vec{\sigma} & \mathcal{A}^{2}\vec{\sigma} & \mathcal{A}^{1}\vec{\sigma} \\ \mathcal{T}^{2\dagger} & \mathcal{A}^{2}\vec{\sigma} & v_{F}\left(\mathbf{k}-\mathbf{q}_{2}\right)\vec{\sigma} & \mathcal{A}^{3}\vec{\sigma} \\ \mathcal{T}^{3\dagger} & \mathcal{A}^{1}\vec{\sigma} & \mathcal{A}^{3}\vec{\sigma} & v_{F}\left(\mathbf{k}-\mathbf{q}_{3}\right)\vec{\sigma} \end{pmatrix}, \end{align} \tag{ 20 }$$

where $q_{1} = k_{\theta}(0,1)$, $q_{2} = k_{\theta}(-\sqrt{3}/2,-1/2)$ and $q_{3} = k_{\theta}(\sqrt{3}/2,-1/2)$ are the momentum differences between the A¹-B¹,A¹-B² and A¹-B³ $K-$point states. The terms $\mathcal{A}^{j}$ and $\mathcal{T}^{j}$ in equation (20) are defined as,

$$\begin{align} \mathcal{T}^{j} & = \omega_{0}\sigma_{0}+\omega_{1}\left[\cos{\left(\frac{2\pi}{3}\left(j-1\right)\right)}\sigma_{x} \right. \notag\\ & \quad \left. + \sin{\left(\frac{2\pi}{3}\left(j-1\right)\right)} \sigma_{y}\right], \notag \\ \mathcal{A}^{j}\vec{\sigma} & = v_{F}\begin{pmatrix} 0 & Ae^{-i\pi/2}e^{-ij2\pi/3} \\ Ae^{i\pi/2}e^{ij2\pi/3} & 0 \\ \end{pmatrix}, \end{align} \tag{ 21 }$$

where A is the magnitude of the gauge field. The C₃ invariance of the Hamiltonian in equation (20) (commonly termed 'the tripod Hamiltonian') imposes that

$$\begin{align} \hat{C}^{\dagger}_{3}\hat{H}\left(R_{2\pi/3}\left[\vec{r}\right]\right)\hat{C}_{3} = \hat{H}\left[\mathbf{A}\right], \end{align} \tag{ 22 }$$

where $R_{2\pi/3}$ is the rotation operator by 120^∘ and $\hat{C}_{3}\equiv e^{\frac{i2\pi}{3}\sigma_{z}}$. Using the perturbative approach presented in appendix C, we formulated an effective low-energy model of tBLG [93] projecting the high-energy bands in equation (20) into the lowest-energy states. The commonly used formula for the renormalized Fermi velocity $v^{*} = v_{F}\frac{(1-3\alpha^{2})}{(1+6\alpha^{2})}, \alpha = \omega/(v_{F}\Delta K)$ for the isotropic model of tBLG [67] ($\omega_{1} = \omega_{0} = \omega$) becomes

$$\begin{align} v^{*}\approx v_{F}\frac{\left(1-\frac{3\omega_{1}^{2}}{\left(v_{F}A+v_{F}\Delta K\right)^{2}}\right)}{\left(1+\frac{6\omega_{1}^{2}}{\left(v_{F}A+v_{F}\Delta K\right)^{2}}\right)}, \end{align} \tag{ 23 }$$

revealing that the effect of atomic relaxation is no other than an effective increase in ΔK. This implies that atomic reconstruction effectively increases the twist angle an amount $\delta \theta = A/K$. Within the linear response theory framework previously presented, the gauge field magnitude is $A = (\frac{3}{4}\frac{\gamma_{0}|\beta |}{v_F}\frac{V_{\mathrm{vdW}}}{\mu\theta} )$ which means an effective increase in twist angle by $\delta \theta = (\frac{9\pi}{16}\frac{\gamma_{0}|\beta |}{v_F}\frac{V_{\mathrm{vdW}}a}{\mu\theta} )$ as clearly shown in see figure 5(a).

The continuum formalism previously developed for twisted bilayer graphene can be generalized to the twisted trilayer structure presented in figure 4(b). In tTLG, the top and bottom layer graphene dispersions are described by two Dirac Hamiltonians rotated by an angle $-\theta/2$. Conversely, the middle layer term consists of a Dirac cone rotated an amount $\theta/2$ in the opposite direction to its top and bottom layer counterpart. While the interlayer tunneling terms connecting neighboring layers remain the same as in tBLG, atomic relaxation leads to significant differences compared to tBLG. Firstly, due to the $z\longleftrightarrow-z$ symmetry of the twisted trilayer structure, the twisted middle layer does not experience any out-of-plane atomic relaxation while the outer ones do. Therefore, the scalar potential $V_{z}(\vec{r})$ originated from the out-of-plane superlattice corrugation is only included at the very top and bottom layer. However, given the enhanced in-plane atomic distortion of the middle layer, the scalar potential induced by the adhesion energies $V_{||}(\vec{r})$ as well as the pseudogauge field $\mathbf{\mathcal{A}}(\vec{r})$, is twice as large in the middle layer as compared to the top and bottom ones. Finally, we account for the energy difference between the dimer and non-dimer pairs formed between the top and bottom graphene layers. Given that both the top and bottom layers lie perfectly AA-aligned to each other, we introduce into the tTLG Hamiltonian a hybridization term connecting the overlaying carbon atoms. This term usually labeled as $|\gamma_{2}|/2 = 8.5$ meV [92] in the Slonczewski–Weiss–McClure terminology, connects the very top and bottom layer via the AA and BB top and bottom layer sublattices,

$$\begin{align} \mathcal{G} = \begin{pmatrix} \gamma_{2}/2 & 0\\ 0 & \gamma_{2}/2 \end{pmatrix}. \end{align} \tag{ 24 }$$

As noted in [94, 95], a negative γ₂ fits better the electronic bandstructure of an ABA trilayer while a positive γ₂ describes more accurately the electronic structure of an ABC (rhombohedral) trilayer graphene. Since our trilayer system hybridizes two AA-stacked graphene monolayers, a negative sign will be used for our calculations. The full tTLG Hamiltonian can therefore be written in a six-state basis $\Psi^{t}_{A},\Psi^{t}_{B},\Psi^{m}_{A},\Psi^{m}_{B},\Psi^{b}_{A},\Psi^{b}_{B}$ as,

$$\begin{align} H_{\mathrm{tTLG}} = \begin{pmatrix} \vec{\sigma}_{-\frac{\theta}{2}}\left(-v_{F}i\vec{{\triangledown}}-\mathcal{A}\right)+\mathcal{E}^{b} -V_{||}-V_{z} & \mathcal{T} & \mathcal{G} \\ \mathcal{T}^{\dagger} & \vec{\sigma}_{\frac{\theta}{2}}\left(-v_{F}i\vec{\triangledown}+2\mathcal{A}\right)+2\mathcal{E}^{t}+2V_{||} & \mathcal{T}^{\dagger} \\ \mathcal{G} & \mathcal{T} & \vec{\sigma}_{-\frac{\theta}{2}}\left(-v_{F}i\vec{\triangledown}-\mathcal{A}\right)+\mathcal{E}^{b}-V_{||}-V_{z} \end{pmatrix}. \end{align} \tag{ 25 }$$

Zoom In Zoom Out Reset image size

Atomic relaxation in tBLG was found to significantly alter its band structure. Firstly, the pseudoscalar potentials were responsible for a slight energy shift between the κ and $\kappa^{^{\prime}}$ points while the pseudomagnetic field widened the Γ-point bandwidth (see figure 5(c)). Furthermore, the dimer-non dimer term Δ shifted the middle bands towards higher energy leading to a greater electron–hole asymmetry appreciable in the density of states (see figure 4(d)).

Zoom In Zoom Out Reset image size

Due to an increase in the local twist angle in the AA regions, we found that a smaller twist angle of the order of 0.8^∘–0.9^∘ was necessary to flatten the middle band of tBLG. Moreover, between 0.8^∘–0.9^∘ nearly constant pseudomagnetic fields of approximately 30 T were found very localized in the vicinity of the AA-regions. These pseudomagnetic fields heavily suppress the Dirac dispersion since the electron wavefunctions of the middle band are also confined around the AA patches where the pseudomagnetic field maxima are located. This additional confining force leads to a very stable flat-band plateau for twist angles ranging between 0.8^∘–0.85^∘. Likewise, the magic angle of fully relaxed tTLG also shifted towards lower angles and formed a magic-angle plateau between 1^∘–1.1^∘ due to the same confining mechanism as in tBLG. Nonetheless, due to the enhanced pseudomagnetic fields in the middle layer AA regions, the flat-band plateau becomes even more pronounced than in tBLG.

We have analyzed in detail the terms that modify the standard model used to study twisted graphene bilayers and symmetrically twisted graphene trilayers.

• We present an analytic scheme, based on linear response, to study the in-plane lattice relaxation, based on two parameters, the Lamé elastic coefficient µ, and the difference in the adhesion potential between the AA and AB regions of a graphene bilayer, $V_{\mathrm{VdW}}$. These two parameters lead to a dimensionless number $V_{\mathrm{vdW}}/ \mu \sim 10^{-4}$ which quantifies the angle dependence of the in-plane deformations and the strength of the resulting pseudomagnetic fields acting on the electrons. It is worth noting that this dimensionless parameter also characterized the competition between the elastic and the adhesion forces in other situations such as in the shape of bubbles [96] or in mechanically deformed nanoribbons [97]. Furthermore, our proposed method can also be generalized to more complicated graphitic systems [42] by adjusting the magnitude of the van der Waals adhesion potential at each interface. Within the linear response approximation, the magnitude of the pseudomagnetic fields is roughly independent of the twist angle B ≈ 49 T. For angles similar or lower to the first magic angle, $\theta \lesssim 1^\circ$, we computed the in-plane atomic displacements by going beyond linear response theory. We have also formulated an effective two-band Hamiltonian capable of quantifying the shift in magic angle due to atomic relaxation.
• The in-plane atomic displacements in the central layer of a symmetric twisted trilayer graphene [4] is increased by a factor of 2 with respect to the relaxation in a twisted bilayer, as the central layer experiences forces of the same sign from the two surface layers. In contrast to tBLG, atomic relaxation in twisted trilayer graphene significantly stabilizes the flat-band plateau around 1^∘ while having two very energetically separated regimes, a middle flat-band similar to the one in tBLG and a highly dispersive Dirac cone above charge neutrality.
• We formulated a simple model based on a few harmonics to describe the pseudoscalar potential due to the change in the interlayer distance. Together with the pseudoscalar field produced by the adhesion forces, the energetic degeneracy between the κ and $\kappa^{^{\prime}}$ points in the Moiré Brillouin zone was lifted by a few meV (see figure 4 and figure 7).
• Finally, we analyzed the effect of a sublattice-dependent potential [73] $\Delta\sim 25$ meV which is known to play a role in the band structure of Bernal bilayers. We found that this term energetically shifted the flat middle band towards higher energies, producing a visibly particle-hole asymmetric density of states profile.

We thank Zhen Zhan and Pierre A Pantaleón for fruitful discussions. F Guinea and A Ceferino acknowledge the support from the Graphene Flagship, Core 3, Grant Number 881603 and from Grants NMAT2D (Comunidad de Madrid, Spain), SprQuMat and (MAD2D-CM)-MRR MATERIALES AVANZADOS-IMDEA-NC. We also acknowledge support from NOVMOMAT, Grant PID2022-142162NB-I00 funded by MCIN/AEI/ 10.13039/501100011033 and, by 'ERDF A way of making Europe'. IMDEA Nanociencia acknowledges support from the 'Severo Ochoa' Programme for Centres of Excellence in R&D (Grant No. SEV-2016-0686).

All data that support the findings of this study are included within the article (and any supplementary files).

To obtain the relaxation-induced atomic displacement field in tBLG, we start with an energy functional describing the sum of the elastic and the adhesion energy,

$$\begin{align} U & = E_{\mathrm{el}}+E_{\mathrm{ad}} = \sum_{n = t,b}\int \left(\left[\frac{\kappa}{2}\left(\frac{\partial h^{\left(n\right)}}{\partial x^{2}} +\frac{\partial h^{\left(n\right)} }{\partial y^{2}}\right)^{2} \right.\right. \notag\\ & \quad +\frac{\left(\mu+\lambda\right)}{2}\left( u^{\left(n\right)}_{xx}+ u^{\left(n\right)}_{yy}\right)^{2} +\frac{\mu}{2}\left(\left(u^{\left(n\right)}_{xx}-u^{\left(n\right)}_{yy}\right)^{2} \right. \notag\\ & \quad \left.\left. +\left(u^{\left(n\right)}_{yx}+u^{\left(n\right)}_{xy}\right)^{2}\right) \vphantom{\left(\frac{\partial h^{\left(n\right)}}{\partial x^{2}} +\frac{\partial h^{\left(n\right)} }{\partial y^{2}}\right)^{2}}\right] +\left[\sum^{i = 3}_{i = 1}V_{\mathrm{vdW}}\cos \left(\mathbf{G}_{i}\mathbf{r} \right. \right. \notag\\ & \quad \left.\left.\left. +\mathbf{g}_{i}\cdot\left(\mathbf{u}^{t}-\mathbf{u}^{b}\right)\right) \vphantom{\sum^{i = 3}_{i = 1}}\right]\right) \mathbf{dR} \nonumber \end{align} \tag{ A1 }$$

where $\mathbf{G} = -\theta\hat{z}\times\mathbf{g}$. Due to a net zero force in each point of the two membranes, we require $u^{top} = -u^{bottom}$ and $h^{top} = -h^{bottom}$. Therefore, we obtain the following equation for the total energy as a function of $u\equiv u^{top}$

$$\begin{align} U & = 2E_{\mathrm{el}}+E_{\mathrm{ad}} = \int \left(2\left[\frac{\kappa}{2}\left(\frac{\partial h}{\partial x^{2}}+\frac{\partial h}{\partial y^{2}}\right)^{2} \right.\right. \notag\\ & \quad +\frac{\left(\mu+\lambda\right)}{2}\left( u_{xx}+ u_{yy}\right)^{2} +\frac{\mu}{2}\left(\left(u_{xx} -u_{yy}\right)^{2} \right. \notag\\ & \quad\left. \left. +\left(u_{yx}+u_{xy}\right)^{2}\right) \vphantom{\left(\frac{\partial h^{\left(1\right)}}{\partial x^{2}}+\frac{\partial h^{\left(1\right)} }{\partial y^{2}}\right)^{2}}\right] +\left[\sum^{i = 3}_{i = 1}V_{\mathrm{vdW}}\cos\left(\mathbf{G}_{i}\mathbf{r}\right. \right. \notag\\ & \quad \left.\left. \left. +2\mathbf{g}_{i}\cdot\mathbf{u}\right) \vphantom{\sum^{i = 3}_{i = 1}}\right]\right) \mathbf{dR}, \end{align} \tag{ A2 }$$

where

$$\begin{align} u_{\alpha\beta} = \frac{1}{2}\left(\partial_{\beta}u_{\alpha}+\partial_{\alpha}u_{\beta}+\partial_{\alpha}h\partial_{\beta}h\right), \end{align} \tag{ A3 }$$

and

$$\begin{align} h\left(\vec{\mathbf{r}}\right) = \left(d+\frac{3h_{0}}{2}\right)+h_{0}\sum^{i = 3}_{i = 1}\cos{\left(\vec{G}_{i}\boldsymbol{\cdot}\vec{r}\right)}. \end{align} \tag{ A4 }$$

To minimize the energy, we use the generalized Euler–Lagrange equations for the field $\vec{u} = (u_{x},u_{y})$,

$$\begin{align} \frac{\partial L}{\partial u_{x}}-\frac{\partial}{\partial x}\frac{\partial L}{\partial \left(\partial_{x}u_{x}\right)}-\frac{\partial}{\partial y}\frac{\partial L}{\partial \left(\partial_{y}u_{x}\right)} & = 0, \notag \\ \frac{\partial L}{\partial u_{y}}-\frac{\partial}{\partial x}\frac{\partial L}{\partial \left(\partial_{x}u_{y}\right)}-\frac{\partial}{\partial y}\frac{\partial L}{\partial \left(\partial_{y}u_{y}\right)} & = 0, \nonumber \end{align} \tag{ A5 }$$

resulting in the following two equations,

$$\begin{align} & \left(\mu+\lambda\right)\frac{\partial_{x}}{2}\left(2\partial_{x}u_{x}+\left(\partial_{x}h\right)^{2}+2\partial_{y}u_{y}+\left(\partial_{y}h\right)^{2}\right) +\mu\frac{\partial_{x}}{2}\left(2\partial_{x}u_{x}+\left(\partial_{x}h\right)^{2}-2\partial_{y}u_{y}-\left(\partial_{y}h\right)^{2}\right) \notag\\ & \quad +\mu\partial_{y}\left(\partial_{x}u_{y}+\partial_{y}u_{x}+\partial_{y}h\partial_{x}h\right) = -\sum^{i = 3}_{i = 1}V_{vdW}\sin{\left(\vec{G}_{i}\vec{r}+2\vec{g}_{i}\vec{u}\right)}g_{ix}, \left(\mu+\lambda\right)\frac{\partial_{y}}{2}\notag\\ & \quad \times \left(2\partial_{x}u_{x}+\left(\partial_{x}h\right)^{2}+2\partial_{y}u_{y}+\left(\partial_{y}h\right)^{2}\right) +\mu\frac{\partial_{y}}{2}\left(-2\partial_{x}u_{x}-\left(\partial_{x}h\right)^{2}+2\partial_{y}u_{y}+\left(\partial_{y}h\right)^{2}\right)\notag\\ & \quad +\mu\partial_{x}\left(\partial_{x}u_{y}+\partial_{y}u_{x}+\partial_{y}h\partial_{x}h\right) = -\sum^{i = 3}_{i = 1}V_{vdW}\sin{\left(\vec{G}_{i}\vec{r}+2\vec{g}_{i}\vec{u}\right)}g_{iy}. \end{align} \tag{ A6 }$$

By simply Fourier transforming equation (A6), we obtain

$$\begin{align} \begin{pmatrix} \delta \mathbf{U}_{x}\left(\vec{q}\right) \\ \delta \mathbf{U}_{y}\left(\vec{q}\right) \\ \end{pmatrix} & = \frac{\begin{pmatrix} \left(2\mu+\lambda\right)q^{2}_{y}+\mu q^{2}_{x} & -\left(\mu+\lambda\right) q_{x}q_{y}\\ -\left(\mu+\lambda\right)q_{x}q_{y} & \left(2\mu+\lambda\right)q^{2}_{x}+\mu q^{2}_{y} \end{pmatrix}}{\left[\left(2\mu+\lambda\right)q^{2}_{x}+\mu q^{2}_{y}\right]\left[\left(2\mu+\lambda\right)q^{2}_{y}+\mu q^{2}_{x}\right]-\left(\left(\mu+\lambda\right)\left(q_{x}q_{y}\right)\right)^{2}} \notag \\ & \quad \times \sqrt{2\pi}\sum_{i,j}\begin{pmatrix} \\ & V_{\mathrm{vdW}}g_{ix}\frac{\left(\delta\left(\vec{q}-\vec{G}_{i}\right)-\delta\left(\vec{q}+\vec{G}_{i}\right)\right)}{2i}+\left[\left(2\mu+\lambda\right)h^{2}_{0}G^{2}_{ix}G_{jx}+\left(\mu+\lambda\right) h^{2}_{0}G_{iy}G_{ix}G_{jy}+h^{2}_{0}\mu G^{2}_{iy}G_{jx}\right] \\ & \times \frac{\left(\delta\left(\vec{q}-\vec{G}_{i}-\vec{G}_{j}\right)-\delta\left(\vec{q}-\vec{G}_{i}+\vec{G}_{j}\right)+\delta\left(\vec{q}+\vec{G}_{i}-\vec{G}_{j}\right)-\delta\left(\vec{q}+\vec{G}_{i}+\vec{G}_{j}\right)\right)}{4i} \\ & \\ & V_{\mathrm{vdW}}g_{iy}\frac{\left(\delta\left(\vec{q}-\vec{G}_{i}\right)-\delta\left(\vec{q}+\vec{G}_{i}\right)\right)}{2i}+\left[(2\mu+\lambda)h^{2}_{0}G^{2}_{iy}G_{jy}+(\mu+\lambda) h^{2}_{0}G_{ix}G_{iy}G_{jx}+h^{2}_{0}\mu G^{2}_{ix}G_{jy}\right] \\ & \times \frac{\Big(\delta(\vec{q}-\vec{G}_{i}-\vec{G}_{j})-\delta(\vec{q}-\vec{G}_{i}+\vec{G}_{j})+\delta(\vec{q}+\vec{G}_{i}-\vec{G}_{j})-\delta(\vec{q}+\vec{G}_{i}+\vec{G}_{j})\Big)}{4i} \end{pmatrix},\nonumber \end{align} \tag{ A7 }$$

which, upon inverse Fourier transform becomes

$$\begin{align} \delta U_{x}\left(\mathbf{x}\right) & = \sum_{i = 1-3}\frac{\left(V_{\mathrm{vdW}}\left(\left(2\mu+\lambda\right)G^{2}_{iy}+\mu G^{2}_{ix}\right)g_{ix}+\left(-V_{vdW}\left(\mu+\lambda\right) G_{ix}G_{iy}\right)g_{iy}\right)\sin{\left(\vec{G}_{i}\vec{r}\right)}}{\mu\left(2\mu+\lambda\right)G^{4}}+\sum_{i = 1-6,j = 1-3}\sin{\left[\left(\vec{G}_{i}+\vec{G}_{j}\right)\vec{r}\right]} \notag\\ & \quad \times \left[\frac{1}{2}\frac{\left(\left[\left(2\mu+\lambda\right)h^{2}_{0}G^{2}_{ix}G_{jx}+\left(\mu+\lambda\right)h^{2}_{0}G_{iy}G_{ix}G_{jy}+h^{2}_{0}\mu G^{2}_{iy}G_{jx}\right]\right)\left(\left(2\mu+\lambda\right)\left(G_{iy}+G_{jy}\right)^{2}+\mu \left(G_{ix}+G_{jx}\right)^{2}\right)}{\mu\left(2\mu+\lambda\right)\left(\left(G_{ix}+G_{jx}\right)^{2}+\left(G_{iy}+G_{jy}\right)^{2}\right)^{2}} \right.\nonumber\\ & \quad \left. -\frac{1}{2}\frac{\left(\left[\left(2\mu+\lambda\right)h^{2}_{0}G^{2}_{iy}G_{jy}+\left(\mu+\lambda\right)h^{2}_{0}G_{ix}G_{iy}G_{jx}+h^{2}_{0}\mu G^{2}_{ix}G_{jy}\right]\right)\left(\mu+\lambda\right)(G_{iy}+G_{jy})(G_{ix}+G_{jx})}{\mu(2\mu+\lambda)((G_{ix}+G_{jx})^{2}+(G_{iy}+G_{jy})^{2})^{2}} \vphantom{\frac{1}{2}\frac{\left(\left[\left(2\mu+\lambda\right)h^{2}_{0}G^{2}_{ix}G_{jx}+\left(\mu+\lambda\right)h^{2}_{0}G_{iy}G_{ix}G_{jy}+h^{2}_{0}\mu G^{2}_{iy}G_{jx}\right]\right)\left(\left(2\mu+\lambda\right)\left(G_{iy}+G_{jy}\right)^{2}+\mu \left(G_{ix}+G_{jx}\right)^{2}\right)}{\mu\left(2\mu+\lambda\right)\left(\left(G_{ix}+G_{jx}\right)^{2}+\left(G_{iy}+G_{jy}\right)^{2}\right)^{2}}}\right], \nonumber \end{align} \tag{ A8 }$$

with $\delta U_{y}(\mathbf{x})$ defined as,

$$\begin{align} & \delta U_{y}\left(\mathbf{x}\right) \notag\\ & \quad = \delta U_{x}\left(\mathbf{x}\right)\left[\left(G_{x},G_{y},A_{x},A_{y}\right)\longrightarrow \left(G_{y},G_{x},A_{y},A_{x}\right)\right]. \end{align} \tag{ A9 }$$

Since the $h^{2}_{0}$ term is much smaller than the $V_{\mathrm{vdW}}$ term, we can simplify equation (A8) as

$$\begin{align} \vec{u}\approx\frac{V_{\mathrm{vdW}}}{\mu G^2} \sum_{i = 1,2,3} \vec{g}_i \sin \left( \vec{G}_i \vec{r} \right). \end{align} \tag{ A10 }$$

Zoom In Zoom Out Reset image size

While equations (A8) and (A9) correctly captures the membrane deformations at very large angles, we can improve this model by plugging equation (A10) into the adhesion term and apply the Jacobi–Anger identity,

$$\begin{align} \cos{\left(z\sin{\left(\theta \right)}\right)} & = J_0\left(z\right)+2\sum _{n = 1}^{\infty } J_{2n}\left(z\right) \cos \left(\text{2n}\theta\right), \notag\\ \sin{\left(z\sin\left(\theta \right)\right)} & = 2 \sum _{n = 1}^{\infty } J_{2n-1}\left(z\right) \sin \left(\text{(2n-1)}\theta\right). \end{align} \tag{ A11 }$$

The adhesion potential then becomes,

$$\begin{align} \sum^{i = 3}_{i = 1}V_{\mathrm{vdW}}\sin{\left(\vec{G}_{i}\vec{r}+2\vec{g}_{i}\vec{u}\right)} & = V_{\mathrm{vdW}}\sum^{i = 3}_{i = 1} \sin{\left(\vec{G}_{i}\vec{r}\right)}\left[\sum^{\infty}_{l = 0}\left(J^{1}_{l}\cos{\left(2l\vec{G}_{1}\vec{r}\right)}\right)\sum^{\infty}_{m = 0}\left(J^{2}_{m}\cos{\left(2m\vec{G}_{2}\vec{r}\right)}\right)\sum^{\infty}_{n = 0}\left(J^{3}_{n}\cos{\left(2n\vec{G}_{3}\vec{r}\right)}\right)\right.\notag\\ & \quad -\left(\sum^{\infty}_{l = 0}J^{1}_{l}\cos{\left(2l\vec{G}_{1}\vec{r}\right)}\right)\left(\sum^{\infty}_{m = 1} J^{2}_{2m-1}\sin{\left(\left(2m-1\right)\vec{G}_{2}\vec{r}\right)}\right)\left(\sum^{\infty}_{n = 1} J^{3}_{2n-1}\sin{\left(\left(2n-1\right)\vec{G}_{3}\vec{r}\right)}\right) \nonumber \\ & \quad -\left(\sum^{\infty}_{l = 1} J^{1}_{2l-1}\sin{\left(\left(2l-1\right)\vec{G}_{1}\vec{r}\right)}\right)\left(\sum^{\infty}_{m = 0} J^{2}_{m}\cos{\left(2m\vec{G}_{2}\vec{r}\right)}\right)\left(\sum^{\infty}_{n = 1} J^{3}_{2n-1}\sin{\left(\left(2n-1\right)\vec{G}_{3}\vec{r}\right)}\right) \nonumber \\ & \quad \left. -\left(\sum^{\infty}_{l = 1} J^{1}_{2l-1}\sin{((2l-1)\vec{G}_{1}\vec{r})}\right)\left(\sum^{\infty}_{m = 1} J^{2}_{2m-1}\sin{((2m-1)\vec{G}_{2}\vec{r})}\right)\left(\sum^{\infty}_{n = 0} J^{3}_{n}\cos{(2n\vec{G}_{3}\vec{r})}\right)\right] \nonumber \\ & \quad + \sum^{i = 3}_{i = 1}V_{\mathrm{vdW}}\cos{(\vec{G}_{i}\vec{r})}\left[\left(\sum^{\infty}_{l = 1} J^{1}_{2l-1}\sin{((2l-1)\vec{G}_{1}\vec{r})}\right)\left(\sum^{\infty}_{m = 0} J^{2}_{m}\cos{(2m\vec{G}_{2}\vec{r})}\right)\left(\sum^{\infty}_{n = 0}J^{3}_{n}\cos{(2n\vec{G}_{3}\vec{r})}\right) \right.\nonumber \\ & \quad +\left(\sum^{\infty}_{l = 0} J^{1}_{l}\cos{(2l\vec{G}_{1}\vec{r})}\right)\left(\sum^{\infty}_{m = 1} J^{2}_{2m-1}\sin{((2m-1)\vec{G}_{2}\vec{r})}\right)\left(\sum^{\infty}_{n = 0}J^{3}_{n}\cos{(2n\vec{G}_{3}\vec{r})}\right) \nonumber \\ & \quad +\left(\sum^{\infty}_{l = 0}J^{1}_{l}\cos{(2l\vec{G}_{1}\vec{r})}\right)\left(\sum^{\infty}_{m = 0}J^{2}_{m}\cos{(2m\vec{G}_{2}\vec{r})}\right)\left(\sum^{\infty}_{n = 1}J^{3}_{2n-1}\sin{((2n-1)\vec{G}_{3}\vec{r})}\right) \nonumber \\ & \quad \left. -\left(\sum^{\infty}_{l = 1} J^{1}_{2l-1}\sin{((2l-1)\vec{G}_{1}\vec{r})}\right)\left(\sum^{\infty}_{m = 1} J^{2}_{2m-1}\sin{((2m-1)\vec{G}_{2}\vec{r})}\right) \left(\sum^{\infty}_{n = 1} J^{3}_{2n-1}\sin{((2n-1)\vec{G}_{3}\vec{r})}\right) \right], \end{align} \tag{ A12 }$$

where $J^{1}_{l},J^{2}_{m},J^{3}_{n}$, are the prefactors proportional to Bessel functions of each harmonic in equation (A11) with arguments $(\frac{2V_{\mathrm{vdW}}}{\mu G^{2}}\vec{g}_{1}\cdot \vec{g}_{i} ), (\frac{2V_{\mathrm{vdW}}}{\mu G^{2}}\vec{g}_{2}\cdot \vec{g}_{i} )$ and $(\frac{2V_{\mathrm{vdW}}}{\mu G^{2}}\vec{g}_{3}\cdot \vec{g}_{i} )$ respectively. The above expression can be simplified as,

$$\begin{align} & \sum^{i = 3}_{i = 1}V_{\mathrm{vdW}}\sin{\left(\vec{G}_{i}\vec{r}+2\vec{g}_{i}\vec{u}\right)} = \sum^{3,\infty\infty,\infty}_{i = 1,l = 0,m = 0,n = 0}A^{\pm l,\pm m,\pm n,i}\sin{\left(\vec{G}_{i}\vec{r}\pm l\vec{G}_{1}\vec{r}\pm m\vec{G}_{2}\vec{r}\pm n\vec{G}_{3}\vec{r}\right)}, \end{align} \tag{ A13 }$$

where

$$\begin{align} \vert A^{\pm l,\pm m,\pm n,i} \vert & = V_{\mathrm{vdW}}J_{|l|}\left(\frac{2V_{\mathrm{vdW}}}{\mu G^{2}}\vec{g}_{1}\cdot \vec{g}_{i}\right)J_{|m|}\left(\frac{2V_{\mathrm{vdW}}}{\mu G^{2}}\vec{g}_{2}\cdot \vec{g}_{i}\right)J_{|n|}\left(\frac{2V_{\mathrm{vdW}}}{\mu G^{2}}\vec{g}_{3}\cdot \vec{g}_{i}\right), \notag \\ A^{\pm l,\pm m,\pm n,i} & = \vert A^{l,m,n,i}\vert Sign\left(\pm l,\pm m,\pm n\right), \end{align} \tag{ A14 }$$

and $Sign(l,m,n)$ is defined as,

$$\begin{align} & \bullet Sign\left(l,m,n\right) = 1~\textit{for no odd indices within (l,m,n)}. \notag\\ & \bullet Sign\left(l,m,n\right) = sgn\left(p\right)~\textit{for } 1 \textit{ odd index p within (l,m,n)}. \notag\\ & \bullet Sign\left(l,m,n\right) = sgn\left(p\right)\times sgn\left(q\right)~\textit{for } 2 \textit{ odd indices p and q within (l,m,n) }.\notag\\ & \bullet Sign\left(l,m,n\right) = sgn\left(l\right)\times sgn\left(m\right)\times sgn\left(n\right)~\textit{for } 3 \textit{ odd indices within (l,m,n)}. \end{align} \tag{ A15 }$$

In defining $A_{x},A_{y}$ as,

$$\begin{align} A^{l,m,n,i}_{x} & = A^{\pm l,\pm m,\pm n,i}g_{x}, \notag\\ A^{l,m,n,i}_{y} & = A^{\pm l,\pm m,\pm n,i}g_{y}, \end{align} \tag{ A16 }$$

we can write the full solution for $(\delta U_{x},\delta U_{y})$ as,

$$\begin{align} \begin{pmatrix} \delta \textbf{U}_{x}(\vec{q}) \\ \delta \textbf{U}_{y}(\vec{q}) \\ \end{pmatrix} & = \frac{\begin{pmatrix} (2\mu+\lambda)q^{2}_{y}+\mu q^{2}_{x} & -(\mu+\lambda) q_{x}q_{y}\\ -(\mu+\lambda)q_{x}q_{y} & (2\mu+\lambda)q^{2}_{x}+\mu q^{2}_{y} \end{pmatrix} }{\left[(2\mu+\lambda)q^{2}_{x}+\mu q^{2}_{y}\right]\left[(2\mu+\lambda)q^{2}_{y}+\mu q^{2}_{x}\right]-((\mu+\lambda)(q_{x}q_{y}))^{2}} \notag \\ & \quad \times \sqrt{2\pi}\begin{pmatrix} & \sum^{3,\infty\infty,\infty}_{i = 1,l = 0,m = 0,n = 0}A^{l,m,n,i}_{x}\frac{(\delta(\vec{q}-(\pm l\vec{G}_{1}\pm m\vec{G}_{2}\pm n\vec{G}_{3}+\vec{G}_{i}))-\delta(\vec{q}+(\pm l\vec{G}_{1}\pm m\vec{G}_{2}\pm n\vec{G}_{3}+\vec{G}_{i}))}{2i} \\ & \\ & \sum^{3,\infty\infty,\infty}_{i = 1,l = 0,m = 0,n = 0}A^{l^{^{\prime}},m^{^{\prime}},n^{^{\prime}},i^{^{\prime}}}_{y}\frac{(\delta(\vec{q}-(\pm l^{^{\prime}}\vec{G}_{1}\pm m^{^{\prime}}\vec{G}_{2}\pm n^{^{\prime}}\vec{G}_{3}+\vec{G}_{i^{^{\prime}}}))-\delta(\vec{q}+(\pm l^{^{\prime}}\vec{G}_{1}\pm m^{^{\prime}}\vec{G}_{2}\pm n^{^{\prime}}\vec{G}_{3}+\vec{G}_{i^{^{\prime}}})))}{2i} \end{pmatrix}, \end{align} \tag{ A17 }$$

leading to

$$\begin{align} \delta U_{x}\left(\vec{x}\right) & = \frac{\sum^{3,\infty\infty,\infty}_{i = 1,l = 0,m = 0,n = 0}\left(\left(2\mu+\lambda\right)G_{l,m,n,iy}^{2}+\mu G_{l,m,n,ix}^{2}\right)A^{\pm l,\pm m,\pm n,i}_{x}\sin{\left(\left(\pm l\vec{G}_{1}\pm m\vec{G}_{2}\pm n\vec{G}_{3}+\vec{G}_{i}\right)\vec{r}\right)}}{\mu\left(2\mu+\lambda\right)G_{l,m,n,i}^{4}} \notag\\ & \quad -\sum^{3,\infty\infty,\infty}_{i^{^{\prime}} = 1,l^{^{\prime}} = 0,m^{^{\prime}} = 0,n^{^{\prime}} = 0}\frac{\left(\mu+\lambda\right)G_{l^{^{\prime}},m^{^{\prime}},n^{^{\prime}},i^{^{\prime}}y}G_{l^{^{\prime}},m^{^{\prime}},n^{^{\prime}},i^{^{\prime}}x}A^{\pm l^{^{\prime}},\pm m^{^{\prime}},\pm n^{^{\prime}},i^{^{\prime}}}_{y}\sin{\left(\pm l^{^{\prime}}\vec{G}_{1}\vec{r}\pm m^{^{\prime}}\vec{G}_{2}\vec{r}\pm n^{^{\prime}}\vec{G}_{3}\vec{r}+\vec{G}_{i^{^{\prime}}}\vec{r}\right)}}{\mu\left(2\mu+\lambda\right)G_{l^{^{\prime}},m^{^{\prime}},n^{^{\prime}},i^{^{\prime}}}^{4}}, \nonumber\\ \delta U_{y}\left(\mathbf{x}\right) & = \delta U_{x}\left(\mathbf{x}\right)\left[\left(G_{x},G_{y},A_{x},A_{y}\right)\longrightarrow \left(G_{y},G_{x},A_{y},A_{x}\right)\right]. \nonumber \end{align} \tag{ A18 }$$

In this appendix, we review how the gauge field $\mathcal{A}$ is derived from the strain field $(\delta U_{x},\delta U_{y})$ following the analysis presented in [24, 69] and we discuss the effect of atomic relaxation in the interlayer hoppings. We start from the usual Slater–Koster form of the nearest-neighboring p_z orbital hopping,

$$\begin{align} t\left(d\right) = V_{\sigma}\left[\frac{\mathbf{d}\cdot\hat{\mathbf{z}}}{d}\right]+V_{\pi}\left[1-\left(\frac{\mathbf{d}\cdot\hat{\mathbf{z}}}{d}\right)^{2}\right], \end{align} \tag{ B1 }$$

where $V_{\sigma}(R) = V^{0}_{\sigma}e^{-(\frac{|R|-a_{cc}}{r_{0}})}$, $V_{\pi}(R) = V^{0}_{\pi}e^{-(\frac{|R|-d}{r_{0}})}$, d is the interlayer distance and $a_{cc} = a/\sqrt{3}$ is the distance between two neighboring carbon atoms in a graphene membrane. Using the monolayer graphene lattice vectors $\vec{a}_{1} = a (\frac{1}{2},\frac{\sqrt{3}}{2} ),\vec{a}_{2} = a (-\frac{1}{2},\frac{\sqrt{3}}{2} )$ and following the analysis presented in [24], we evaluate down to first order in u the change in the three nearest-neighboring intralayer p_z orbital hopping,

$$\begin{align} \delta t_{1} & = -t\beta u_{yy}+O\left(u^{2}\right), \notag\\ \delta t_{2} & = -t\beta \left(\frac{3}{4}u_{xx}+\frac{\sqrt{3}}{4}u_{xy}+\frac{\sqrt{3}}{4}u_{yx} +\frac{1}{4}u_{yy}\right) \notag\\ & \quad + O\left(u^{2}\right), \nonumber \\ \delta t_{3} & = -t\beta\left(\frac{3}{4}u_{xx}-\frac{\sqrt{3}}{4}u_{xy}-\frac{\sqrt{3}}{4}u_{yx} +\frac{1}{4}u_{yy}\right) \notag\\ & \quad + O\left(u^{2}\right), \nonumber \end{align} \tag{ B2 }$$

where $\beta = \frac{a_{cc}}{r_{0}}$ and $t = -\gamma_{0} = -2.6$ eV (see [83]). The tight-binding Hamiltonian of monolayer graphene is well known to read as

$$\begin{align} H_{G} = \begin{pmatrix} 0 & g\left(k\right)\\ g\left(k\right) & 0 \end{pmatrix}, \end{align} \tag{ B3 }$$

where $g(k) = t+\sum^{i = 2}_{i = 1} te^{i\mathbf{k}\mathbf{a}_{i}}$. Expanding H_G around $\vec{K} = (\frac{4\pi}{3a},0 )$ and taking into account the intralayer hopping corrections $\delta t_{1},\delta t_{2}$ and $\delta t_{3}$, in g(k), we obtain the following modified graphene Hamiltonian

$$\begin{align} H_{G} = \begin{pmatrix} 0 & v\left(k_{x}-ik_{y}\right)+\left(A_{x}-iA_{y}\right)\\ v\left(k_{x}+ik_{y}\right)+\left(A_{x}+iA_{y}\right) & 0 \end{pmatrix}, \end{align} \tag{ B4 }$$

where $v\equiv\frac{-\sqrt{3}}{2}ta$ and

$$\begin{align} A_{x} & = \frac{3\beta t}{4}\left(\frac{\partial u_{x}}{\partial x}-\frac{\partial u_{y}}{\partial y}\right)+O\left(u^{2}\right), \nonumber\\ A_{y} & = -\frac{3\beta t}{4}\left(\frac{\partial u_{x}}{\partial y}+\frac{\partial u_{y}}{\partial x}\right)+O\left(u^{2}\right). \end{align} \tag{ B5 }$$

Following the above tight-binding approach, we review the effect that relaxation has in the interlayer tunneling terms following the expansion performed in [18, 24, 91]. Starting from the Bloch state

$$\begin{align} \vert \mathbf{k},X, l\rangle = \frac{1}{\sqrt{N}}\sum_{R}e^{i\mathbf{k}\mathbf{R}}\vert \mathbf{R}+\mathbf{u}\rangle, \end{align} \tag{ B6 }$$

we firstly define the 2D Fourier transform of the interlayer hopping $T(\vec{r})$ where $\vec{r}$ is the distance between two carbon atoms,

$$\begin{align} t\left(p\right) = \frac{1}{S_{0}d_{0}}\int d^{3}r T\left(\mathbf{r}\right)e^{-i\mathbf{p}\mathbf{r}}. \end{align} \tag{ B7 }$$

We therefore obtain the following interlayer tunneling matrix element,

$$\begin{align} & \langle \mathbf{k}^{^{\prime}},X^{^{\prime}},2\vert U \vert \mathbf{k},X,1\rangle \notag\\ & \quad = \frac{-1}{N}\sum_{R\in l^{1}}\sum_{R^{^{\prime}}\in l^{2}}e^{i\mathbf{k}\mathbf{R}-i\mathbf{k}^{^{\prime}}\mathbf{R}^{^{\prime}}}T\left(\mathbf{R}^{^{\prime}}+\mathbf{u}^{\left(2\right)}-\mathbf{R}-\mathbf{u}^{\left(1\right)}\right)\notag\\ & \quad = \frac{-S_{0}d_{0}}{N\left(2\pi\right)^{3}}\sum_{R\in l^{1}}\sum_{R^{^{\prime}}\in l^{2}}\int d^{3}p t\left(\mathbf{p}\right)e^{i\left(\mathbf{k}-\mathbf{p}\right)\mathbf{R}-i\left(\mathbf{k}^{^{\prime}}-\mathbf{p}\right)\mathbf{R}^{^{\prime}}}\notag\\ & \qquad \times e^{i\mathbf{p}\left(\mathbf{u_{X^{^{\prime}}}}^{2}-\mathbf{u_{X}}^{1}\right)}, \nonumber \end{align} \tag{ B8 }$$

where S₀ and d₀ are the unit cell and interlayer distance respectively and $\mathbf{u}^{1(2)}$ the atomic displacement fields of the top and bottom layer respectively. We then expand the exponential function $e^{i\mathbf{p}\mathbf{u}^{(l)}}$ and obtain

$$\begin{align} & \langle \mathbf{k}^{^{\prime}}, X^{^{\prime}}, 2\vert U \vert \mathbf{k}, X, 1\rangle \notag\\ & \quad = \sum_{\mathbf{g,g^{^{\prime}}}}\sum_{n_{1},n_{2},\ldots}\sum_{n^{^{\prime}}_{1},n^{^{\prime}}_{2},\ldots} \Gamma^{n^{^{\prime}}_{1},n^{^{\prime}}_{2}\ldots}_{n_{1},n_{2},\ldots}\left(\mathbf{Q}\right)e^{i\left(-\mathbf{g}{\tau}^{1}_{X}+\mathbf{g^{^{\prime}}}{\tau}^{2}_{X^{^{\prime}}}\right)}\notag\\ & \qquad \times \delta_{\mathbf{k}+\mathbf{g}+n_{1}\mathbf{q}_{1}+n_{2}\mathbf{q}_{2}+\cdots,\mathbf{k}^{^{\prime}}+\mathbf{g^{^{\prime}}}-n^{^{\prime}}_{1}\mathbf{q}_{1} -n^{^{\prime}}_{2}\mathbf{q}_{2}+\cdots}, \end{align} \tag{ B9 }$$

where $\mathbf{g},\mathbf{g^{^{\prime}}}$ are the reciprocal lattice vectors of layer of layer 1 and 2 respectively, q are the Fourier modes of $\mathbf{u}^{(l)}$ and $\mathbf{Q} = \mathbf{k}+\mathbf{g}+n_{1}\mathbf{q}_{1}+n_{2}\mathbf{q}_{2}+\cdots$ are the resonant momentum transfer between the k and k' state. Finally, the term Γ is defined as,

$$\begin{align} \Gamma^{n^{^{\prime}}_{1},n^{^{\prime}}_{2}\ldots}_{n_{1},n_{2},\ldots}& = -t_{||}\left(\mathbf{Q},d_{0}\right)\left[\frac{-i\mathbf{Q}\mathbf{u}^{1}_{X\mathbf{q}_{1}}}{n_{1}!}\right]^{n_{1}}\left[\frac{-i\mathbf{Q}\mathbf{u}^{1}_{X\mathbf{q}_{2}}}{n_{2}!}\right]^{n_{2}}\cdots \notag\\ & \quad \times\left[\frac{-i\mathbf{Q}\mathbf{u}^{2}_{X^{^{\prime}}\mathbf{q}_{1}}}{n^{^{\prime}}_{1}!}\right]^{n^{^{\prime}}_{1}}\left[\frac{-i\mathbf{Q}\mathbf{u}^{2}_{X^{^{\prime}}\mathbf{q}_{2}}}{n^{^{\prime}}_{2}!}\right]^{n^{^{\prime}}_{2}}\cdots, \end{align} \tag{ B10 }$$

where $t_{||}(Q,d_{0})$ is the Q Fourier mode of T(r) at a fixed interlayer distance d₀. Equation (B10) clearly shows the very quick decay of the corrections in the interlayer tunneling terms with an increasing number of harmonics.

Zoom In Zoom Out Reset image size

To find the renormalized Fermi velocity at sufficiently large twist angles θ, we will first formulate a low-energy 2-band model of tBLG around the K-point accounting for atomic relaxation. We start by truncating the BM Hamiltonian into a simpler $8\times 8$ matrix Hamiltonian which only considers the resonant tunneling between the bottom layer Dirac cone and its three nearest top layer counterparts. This approximate model was demonstrated in [67] to provide a very accurate prediction for the first magic angle in tBLG. As previously discussed in section 3, lattice relaxation adds a three-fold symmetric pseudogauge field $\mathcal{A}$ connecting Dirac cones separated by integer multiples of the Moiré reciprocal lattice vectors $\vec{G}_{1} = \frac{4\pi}{\sqrt{3}L_{M}}(1,0),\vec{G}_{2} = \frac{4\pi}{\sqrt{3}L_{M}}(-1/2,\sqrt{3}/2)$ and $\vec{G}_{3} = \frac{4\pi}{\sqrt{3}L_{M}}(-1/2,-\sqrt{3}/2)$. Therefore, the tripod Hamiltonian of fully relaxed tBLG written in the basis $(\Psi^{A1},\Psi^{B1},\Psi^{B2},\Psi^{B3})$ designating the four Dirac cones depicted in figure 5(a) becomes

$$\begin{align} H\Psi & = \begin{pmatrix} v_{F}\mathbf{k}\vec{\sigma} & \mathcal{T}^{1} & \mathcal{T}^{2} & \mathcal{T}^{3} \\ \mathcal{T}^{1\dagger} & v_{F}\left(\mathbf{k}-\mathbf{q}_{1}\right)\vec{\sigma} & \mathcal{A}^{2}\vec{\sigma} & \mathcal{A}^{1}\vec{\sigma} \\ \mathcal{T}^{2\dagger} & \mathcal{A}^{2}\vec{\sigma} & v_{F}\left(\mathbf{k}-\mathbf{q}_{2}\right)\vec{\sigma} & \mathcal{A}^{3}\vec{\sigma} \\ \mathcal{T}^{3\dagger} & \mathcal{A}^{1}\vec{\sigma} & \mathcal{A}^{3}\vec{\sigma} & v_{F}\left(\mathbf{k}-\mathbf{q}_{3}\right)\vec{\sigma} \end{pmatrix} \notag\\ & \quad \times \begin{pmatrix} \Psi^{A1} \\ \Psi^{B1} \\ \Psi^{B2} \\ \Psi^{B3}\end{pmatrix} = E\begin{pmatrix} \Psi^{A1} \\ \Psi^{B1} \\ \Psi^{B2} \\ \Psi^{B3} \end{pmatrix}, \end{align} \tag{ C1 }$$

where the interlayer tunneling matrices $\mathcal{T}^{j}$ and the gauge field $\mathcal{A}^{j}$ were previously defined in equation (21) of section 3 while the vectors $\vec{q}_{1},\vec{q}_{2},\vec{q}_{3}$ are $\Delta K(0,1),\Delta K(-\sqrt{3}/2,-1/2)$ and $\Delta K(\sqrt{3}/2,-1/2)$ respectively. The $8\times 8$ Hamiltonian in equation (C1) is then subdivided into two square quadrants Q_l and Q_h ,

$$\begin{align} Q_{l} & = v_{F}\mathbf{k}\vec{\sigma},\notag\\ & \notag\\[-10pt] Q_{h} & = \begin{pmatrix} v_{F}\left(\mathbf{k}-\mathbf{q}_{1}\right)\vec{\sigma} & \mathcal{A}^{2}\vec{\sigma} & \mathcal{A}^{1}\vec{\sigma} \\ \mathcal{A}^{2}\vec{\sigma} & v_{F}\left(\mathbf{k}-\mathbf{q}_{2}\right)\vec{\sigma} & \mathcal{A}^{3}\vec{\sigma} \\ \mathcal{A}^{1}\vec{\sigma} & \mathcal{A}^{3}\vec{\sigma} & v_{F}\left(\mathbf{k}-\mathbf{q}_{3}\right)\vec{\sigma} \end{pmatrix}, \end{align} \tag{ C2 }$$

which are connected by the interlayer tunneling terms $\mathcal{T}^{j}$. The matrices Q_l and Q_h can be identified as the low and the high-energy states of the tripod Hamiltonian. Considering the diagonal $v_{F}\mathbf{k}\vec{\sigma}$ terms in Q_h as a perturbation, we perform a unitary transformation $\hat{U}$ in Q_h which diagonalizes the Q_h matrix,

$$\begin{align} \tilde{\mathbf{Q}}_{h} & = \hat{U}^{-1}Q_{h}\hat{U} = \begin{pmatrix} \left(v_{F}A+v_{F}\Delta K\right)\mathbb{I}_{2} & 0 & 0 \\ 0 & -\left(v_{F}A+v_{F}\Delta K\right)\mathbb{I}_{2} & 0 \\ 0 & 0 & \left(-2v_{F}A+v_{F}\Delta K\right)\sigma_{z} \end{pmatrix}, \text{where} \notag\\ U & = \begin{pmatrix} -\frac{1}{2} e^{-\frac{i\pi}{4}} & \frac{e^{-\frac{i\pi}{4}}}{2 \sqrt{3}} & \frac{1}{2} e^{-\frac{i\pi}{4}} & -\frac{e^{-\frac{i\pi}{4}}}{2 \sqrt{3}} & \frac{e^{-\frac{i\pi}{4}}}{\sqrt{6}} & -\frac{e^{-\frac{i \pi}{4}}}{\sqrt{6}} \\ 0 & \frac{e^{-\frac{i\pi}{12}}}{\sqrt{3}} & 0 & \frac{e^{-\frac{i\pi}{12}}}{\sqrt{3}} & -\frac{e^{-\frac{i\pi}{12}}}{\sqrt{6}} & -\frac{e^{-\frac{i\pi}{12}}}{\sqrt{6}} \\ 0 & \frac{e^{\frac{i\pi }{12}}}{\sqrt{3}} & 0 & -\frac{e^{\frac{i \pi }{12}}}{\sqrt{3}} & -\frac{e^{\frac{i \pi }{12}}}{\sqrt{6}} & \frac{e^{\frac{i \pi }{12}}}{\sqrt{6}} \\ \frac{1}{2} e^{-\frac{5i\pi}{12}} & \frac{e^{-\frac{5i\pi}{12}}}{2 \sqrt{3}} & \frac{1}{2} e^{-\frac{5i\pi}{12}} & \frac{e^{-\frac{5i\pi}{12}}}{2\sqrt{3}} & \frac{e^{-\frac{5i\pi}{12}}}{\sqrt{6}} & \frac{e^{-\frac{5i\pi}{12}}}{\sqrt{6}} \\ \frac{1}{2} e^{\frac{5i\pi }{12}} & \frac{e^{\frac{5i\pi }{12}}}{2 \sqrt{3}} & -\frac{1}{2} e^{\frac{5i\pi }{12}} & -\frac{e^{\frac{5i\pi }{12}}}{2 \sqrt{3}} & \frac{e^{\frac{5i\pi }{12}}}{\sqrt{6}} & -\frac{e^{\frac{5i\pi }{12}}}{\sqrt{6}} \end{pmatrix}.\nonumber \end{align} \tag{ C3 }$$

Regarding Q_l , we projected this $2\times2$ Hamiltonian in the orthogonal basis defined by its two eigenvectors $\Psi_{c} = \frac{1}{\sqrt{2}}\begin{pmatrix} e^{-i\theta_{k}/2} \\ e^{i\theta_{k}/2} \end{pmatrix}$ and $\Psi_{v} = \frac{1}{\sqrt{2}}\begin{pmatrix} e^{-i\theta_{k}/2} \\ -e^{i\theta_{k}/2} \end{pmatrix}$, where $\theta_{k}\equiv\arctan(k_{y}/k_{x})$. Hence, the diagonal matrix emanating from this projection has the form,

$$\begin{align} \tilde{\mathbf{Q}}_{l} & = O^{-1}_{\theta_{k}}Q_{l}O_{\theta_{k}} = \begin{pmatrix} v_{F}k & 0 \\ 0 & -v_{F}k \end{pmatrix} \text{where}\notag\\ O_{\theta_{k}} & = \frac{1}{\sqrt{2}}\begin{pmatrix} e^{-i\theta_{k}/2} & e^{-i\theta_{k}/2} \\ e^{i\theta_{k}/2} & -e^{i\theta_{k}/2} \end{pmatrix}. \end{align} \tag{ C4 }$$

For simplicity, we will fix the momentum orientation θ_k in the $\kappa-\kappa^{^{\prime}}$ direction ($\theta_{k} = \pi/2$), neglecting the trigonal symmetry around the Dirac point. From the matrices U and $O_{\theta_{k}}$, we construct a projection operator $\mathcal{P}$ defined as

$$\begin{align} \mathcal{P}\equiv\begin{pmatrix} O_{\theta_{k}} & \tilde{0} \\ \tilde{0} & U\end{pmatrix}, \end{align} \tag{ C5 }$$

where $\tilde{0}$ is a $2\times 6$ zero matrix. Considering the limit of $|\vec{k}|\longrightarrow 0$ in $\tilde{\mathbf{Q}}_{l}$, we split H_trip into two parts, H₀ and H' labeling the bare and the perturbative part of the tripod Hamiltonian. Within H₀ we have the diagonal matrices $\tilde{\mathbf{Q}}_{l}$ in the limit of $\mathbf{k}\longrightarrow0$ as well as $\tilde{\mathbf{Q}}_{h}$ and the projected interlayer tunneling terms which become $\mathcal{T}^{j}\longrightarrow O_{\pi/2}^{-1}\mathcal{T}^{j}U$,

$$\begin{align} H_{0} = \begin{pmatrix} H^{t}_{0} & \boldsymbol{0} \\ \boldsymbol{0}^{T} & H^{b}_{0} \end{pmatrix}, \end{align} \tag{ C6 }$$

where 0 is a $6\times 2$ zero matrix independent of θ_k , $H^{b}_{0} = (-2A+v_{F}\Delta K)\sigma_{z}$ and

$$\begin{align} H^{t} = \left( \begin{array}{cccccc} 0 & 0 & \frac{\sqrt{3}\omega_{0}-3\omega_{1}}{2 \sqrt{2}} & \frac{1}{2} \sqrt{\frac{3}{2}} \left(\sqrt{3} \omega_{0}+\omega_{1}\right) & -\frac{i \left(3\omega_{0}-\sqrt{3}\omega_{1}\right)}{2 \sqrt{2}} & \frac{1}{2} i \sqrt{\frac{3}{2}} \left(\omega_{0}+\sqrt{3}\omega_{1}\right) \\ 0 & 0 & -\frac{i \left(3\omega_{0}+\sqrt{3}\omega_{1}\right)}{2 \sqrt{2}} & \frac{1}{2} i \sqrt{\frac{3}{2}} \left(\omega_{0}-\sqrt{3}\omega_{1}\right) & \frac{\sqrt{3}\omega_{0}+3\omega_{1}}{2 \sqrt{2}} & \frac{1}{4} \left(3 \sqrt{2}\omega_{0}-\sqrt{6}\omega_{1}\right) \\ \frac{\sqrt{3}\omega_{0}-3\omega_{1}}{2 \sqrt{2}} & \frac{i\left(3\omega_{0}+\sqrt{3}\omega_{1}\right)}{2 \sqrt{2}} & v_{F}A+v_{F}\Delta K & 0 & 0 & 0 \\ \frac{3\omega_{0}+\sqrt{3}\omega_{1}}{2 \sqrt{2}} & -\frac{i \left(\sqrt{3}\omega_{0}-3\omega_{1}\right)}{2 \sqrt{2}} & 0 & v_{F}A+v_{F}\Delta K & 0 & 0 \\ \frac{i \left(3\omega_{0}-\sqrt{3}\omega_{1}\right)}{2 \sqrt{2}} & \frac{\sqrt{3}\omega_{0}+3\omega_{0}}{2 \sqrt{2}} & 0 & 0 & -\left(v_{F}A+v_{F}\Delta K\right) & 0 \\ -\frac{i \left(\sqrt{3}\omega_{0}+3\omega_{1}\right)}{2 \sqrt{2}} & \frac{3 \omega_{0}-\sqrt{3}\omega_{1}}{2 \sqrt{2}} & 0 & 0 & 0 & -\left(v_{F}A+v_{F}\Delta K\right) \\ \end{array} \right). \end{align} \tag{ C7 }$$

Since there is no coupling between the zero-energy states in $H^{t}_{0}$ and the high-energy levels in $H^{b}_{0}$, we will neglect $H^{b}_{0}$ in the unperturbed Hamiltonian such that $H_{0} = H^{t}_{0}$.

The perturbative term H' consist of the diagonal $v_{F}\mathbf{k}\vec{\sigma}$ terms projected in the orthogonal basis $\mathcal{P}$ defined in equation (C5). Akin to H₀, the coupling terms between the $v_{F}A+v_{F}\Delta K$ and the $(-2v_{F}A+v_{F}\Delta K)$ states are neglected as they would contribute quadratically to the momentum dependence of the low-energy dispersion of $H^{t}_{0}$. Therefore, the perturbative H' term becomes

$$\begin{align} H^{^{\prime}} = v_{F}\left( \begin{array}{cccccc} k & 0 & 0 & 0 & 0 & 0 \\ 0 & -k & 0 & 0 & 0 & 0 \\ 0 & 0 & -\frac{k}{4} & \frac{\sqrt{3}k}{4} & -\frac{1}{4} i \sqrt{3}k & -\frac{1}{4}\left(ik\right) \\ 0 & 0 & \frac{\sqrt{3}k}{4} & \frac{k}{4} & -\frac{1}{4}\left(ik\right) & \frac{1}{4}i\sqrt{3}k \\ 0 & 0 & \frac{1}{4}i\sqrt{3}k & \frac{i k}{4} & \frac{k}{4} & -\frac{1}{4} \left(\sqrt{3}k\right) \\ 0 & 0 & \frac{i k}{4} & -\frac{1}{4} i \sqrt{3} k & -\frac{1}{4} \left(\sqrt{3} k\right) & -\frac{k}{4} \\ \end{array} \right). \end{align} \tag{ C8 }$$

The doubly degenerate zero-energy eigenstates of H₀ were obtained by diagonalizing H₀,

$$\begin{align} \begin{aligned} \Psi^{c} & = \frac{1}{\sqrt{\Big(2+\frac{2(v_{F}A+v_{F}\Delta K)^{2}}{3(\omega^{2}_{0}+\omega^{2}_{1})}\Big)}} \begin{pmatrix} \frac{i (v_{F}A+v_{F}\Delta K) \left(\omega^{2}_{0}+2 \sqrt{3}\omega_{0}\omega_{1}+3\omega_{1}^2\right)}{\sqrt{2}(\sqrt{3}\omega_{0}+3\omega_{1})(\omega^{2}_{0}+\omega_{1}^2} \\ \frac{(v_{F}A+v_{F}\Delta K) \left(3\omega_{0}-\sqrt{3}\omega_{1}\right)}{3 \sqrt{2} \left(\omega^{2}_{0}+\omega^{2}_{1}\right)}\\ -\frac{i \left(\omega^{2}_{0}-\omega^{2}_{1}\right)}{\omega^{2}_{0}+\omega^{2}_{1}}\\ -\frac{2 i \omega_{0}\omega_{1}}{\omega^{2}_{0}+\omega^{2}_{1}} \\ 0 \\ 1 \\ \end{pmatrix},\\ \Psi^{v} & = \frac{1}{\sqrt{\Big(2+\frac{2(v_{F}A+v_{F}\Delta K)^{2}}{3(\omega^{2}_{0}+\omega^{2}_{1})}\Big)}}\begin{pmatrix} -\frac{i (v_{F}A+v_{F}\Delta K) \left(3\omega_{0}-\sqrt{3}\omega_{1}\right)}{3 \sqrt{2} \left(\omega^{2}_{0}+\omega^{2}_{1}\right)} \\ \frac{(v_{F}A+v_{F}\Delta K) \left(\sqrt{3}\omega_{0}+3\omega_{1}\right)}{3 \sqrt{2} \left(\omega^{2}_{0}+\omega^{2}_{1}\right)} \\ -\frac{2 i\omega_{0}\omega_{1}}{\omega^{2}_{0}+\omega^{1}_{0}} \\ \frac{i \left(\omega^{2}_{0}-\omega^{2}_{1}\right)}{\omega^{2}_{0}+\omega^{2}_{1}} \\ 1 \\ 0 \\ \end{pmatrix}. \end{aligned} \end{align} \tag{ C9 }$$

This new set of eigenstates permits to construct a 2-band model projecting H' into the orthogonal basis ($\Psi^{c},\Psi^{v}$),

$$\begin{align} H_{\mathrm{eff}} & = \begin{pmatrix} \langle \Psi^{c} |H^{^{\prime}}| \Psi^{c} \rangle & \langle \Psi^{c} |H^{^{\prime}}| \Psi^{v}\rangle \\ \langle \Psi^{v} |H^{^{\prime}}| \Psi^{c} \rangle & \langle \Psi^{v} |H^{^{\prime}}| \Psi^{v} \rangle \end{pmatrix} \notag\\ & = v_{F}\begin{pmatrix} \frac{3 k\left(-\omega^{4}_{0}+10\omega_{0}^2 \omega_{1}^2+8\sqrt{3}\omega_{0}\omega_{1}^{3}+3\omega_{1}^{4}\right)\left(\left(v_{F}A+v_{F}\Delta K\right)^{2}-3\omega_{1}^2\right)}{2 \left(\sqrt{3}\omega_{0}+3\omega_{1}\right)^{2} \left(\omega_{0}^{2}+\omega_{1}^{2}\right)\left(\left(v_{F}A+v_{F}\Delta K\right)^{2}+3 \left(\omega^{2}_{0}+\omega^{2}_{1}\right)\right)} & \frac{k\left(-3\omega_{0}^{3}-5\sqrt{3}\omega_{0}^{2}\omega_{1}-3 \omega_{0}\omega^{2}_{1}+3\sqrt{3}\omega_{1}^3\right) \left(\left(v_{F}A+v_{F}\Delta K\right)^{2}-3\omega^{2}_{1}\right)}{2\left(\sqrt{3}\omega_{0}+3\omega_{1}\right)\left(\omega_{0}^2+\omega_{1}^2\right)\left(\left(v_{F}A+v_{F}\Delta K\right)^2+3\left(\omega^{2}_{0}+\omega^{2}_{1}\right)\right)} \\ \frac{k\left(-3\omega_{0}^{3}-5\sqrt{3}\omega_{0}^{2}\omega_{1}-3\omega_{0} \omega_{1}^2+3\sqrt{3}\omega_{1}^3\right)\left(\left(v_{F}A+v_{F}\Delta K\right)^{2}-3 \omega_{1}^2\right)}{2\left(\sqrt{3} \omega_{0}+3\omega_{1}\right)\left(\omega_{0}^2+\omega_{1}^2\right) \left(\left(v_{F}A+v_{F}\Delta K\right)^{2}+3\left(\omega_{0}^{2}+\omega_{1}^2\right)\right)} & \frac{k (\omega_{0}^2-2 \sqrt{3}\omega_{0}\omega_{1}-\omega_{1}^2)((v_{F}A+v_{F}\Delta K)^2-3\omega_{1}^2)}{2(\omega_{0}^2+\omega_{1}^2)((v_{F}A+v_{F}\Delta K)^{2}+3 (\omega^{2}_{0}+\omega^{2}_{1}))}. \end{pmatrix}. \nonumber \end{align} \tag{ C10 }$$

Diagonalizing $H_{\mathrm{eff}}$, the Fermi velocity was obtained for three different situations. Firstly, when $\omega_{1} = \omega_{0}$, secondly for $\omega_{1}\neq\omega_{0}$ and thirdly for $\omega_{0} = 0$,

$$\begin{align} v_{\omega_{0} = \omega_{1}} & = v_{F}\frac{(1-3\alpha^{^{\prime} 2})}{(1+6\alpha^{^{\prime} 2})}, \notag \\ v_{\omega_{0}\neq\omega_{1}} & = v_{F}\frac{3(-(v_{F}A+v_{F}\Delta K)^{2}+3\omega^{2}_{1})\sqrt{(\omega^{2}_{0}+\omega^{2}_{1})^{2}(\omega^{4}_{0}+4\sqrt{3}\omega^{3}_{0}\omega_{1}+18\omega^{2}_{0}\omega^{2}_{1}+12\sqrt{3}\omega_{0}\omega^{3}_{1}+9\omega^{4}_{1})}}{(\sqrt{3}\omega_{0}+3\omega_{1})(\omega^{2}_{0}+\omega^{2}_{1})((v_{F}A+v_{F}\Delta K)^{2}+3(\omega^{2}_{0}+\omega^{2}_{1})}, \nonumber\\ v_{\omega_{0} = 0} & = v_{F}\frac{(1-3\alpha^{^{\prime} 2})}{(1+3\alpha^{^{\prime} 2})}. \nonumber \end{align} \tag{ C11 }$$

where $\alpha^{^{\prime}} = \omega_{1}/(v_{F}A+v_{F}\Delta K)$.We noted that in the limit of no relaxation $(A\longrightarrow 0)$, the expressions for $v_{\mathrm{chiral}}$ and $v_{\omega_{1} = \omega_{0}}$ coincided with the ones reported in [67, 71] respectively. Furthermore, from equations (C7) and (C11) it is transparent that the effect of atomic relaxation is no other than a simple replacement of $v_{F}\Delta K$ by $v_{F}\Delta K+v_{F}A$, effectively increasing the twist angle in the AA-sites.