Moir\'e-pattern fluctuations and electron-phason coupling in twisted bilayer graphene

TL;DR

This paper investigates how moiré-pattern fluctuations and electron-phason interactions in twisted bilayer graphene influence electronic properties, revealing mechanisms that affect Landau level degeneracy and resistivity, especially near magic angles.

Contribution

It introduces a detailed analysis of electron-phason coupling in twisted bilayer graphene and its effects on electronic degeneracy and resistivity, highlighting the role of soliton network softness.

Findings

Electron-phason coupling lifts Dirac cone degeneracy.

Resistivity shows linear temperature dependence influenced by twist angle.

Soliton network stiffness decreases near magic angles.

Abstract

In twisted bilayer graphene, long-wavelength lattice fluctuations on the scale of the moir\'e period are dominated by phason modes, i.e., acoustic branches of the incommensurate lattice resulting from coherent superpositions of optical phonons. In the limit of small twist angles, these modes describe the sliding motion of stacking domain walls separating regions of partial commensuration. The resulting soliton network is a soft elastic manifold, whose reduced rigidity arises from the competition between intralayer (elastic) and interlayer (adhesion) forces governing lattice relaxation. Shear deformations of the beating pattern dominate the electron-phason coupling to the leading order in $t_{\perp}/t$, the ratio between interlayer and intralayer hopping parameters. This coupling lifts the layer degeneracy of the Dirac cones at the corners of the Moir\'e Brillouin zone, which could explain the observed 4-fold (instead of 8-fold) Landau level degeneracy. Electron-phason scattering gives rise to a linear-in-temperature contribution to the resistivity that increases with decreasing twist angle due to the reduction of the stiffness of the soliton network. This contribution alone, however, seems to be insufficient to explain the huge enhancement of the resistivity of the normal state close to the magic angle.

Tables

Table 1. Table 1: Classification of electronic valley- and layer-diagonal operators according to irreducible representations (irreps) of D 6 subscript 𝐷 6 D_{6} and parity ( ± 1 plus-or-minus 1 \pm 1 , even/odd) under timer-reversal symmetry.

Electronic operators ($𝒯$)	Irrep	$ℰ$	${𝒞}_{2z}$	$2{𝒞}_{3z}$	$2{𝒞}_{6z}$	$3{𝒞}_{2x}$	$3{𝒞}_{2y}$
${\hat{\mu }}_z{\hat{\tau }}_z{\hat{\sigma }}_z(-)$	$A_1$	1	1	1	1	1	1
${\hat{\mu }}_z(+)$, ${\hat{\tau }}_z{\hat{\sigma }}_z(-)$	$A_2$	1	1	1	1	-1	-1
${\hat{\tau }}_z(-)$, ${\hat{\mu }}_z{\hat{\sigma }}_z(+)$	$B_1$	1	-1	1	- 1	1	-1
${\hat{\sigma }}_z(+)$, ${\hat{\mu }}_z{\hat{\tau }}_z(-)$	$B_2$	1	-1	1	-1	-1	1
$\left(\begin{array}{c}{\hat{\tau }}_z{\hat{\sigma }}_x\\ {\hat{\sigma }}_y\end{array}\right)(-)$,  $\left(\begin{array}{c}-{\hat{\mu }}_z{\hat{\sigma }}_y\\ {\hat{\mu }}_z{\hat{\tau }}_z{\hat{\sigma }}_x\end{array}\right)(-)$	$E_1$	2	-2	-1	1	0	0
$\left(\begin{array}{c}{\hat{\sigma }}_x\\ {\hat{\tau }}_z{\hat{\sigma }}_y\end{array}\right)(+)$ , $\left(\begin{array}{c}-{\hat{\mu }}_z{\hat{\tau }}_z{\hat{\sigma }}_y\\ {\hat{\mu }}_z{\hat{\sigma }}_x\end{array}\right)(+)$	$E_2$	2	2	-1	-1	0	0

Equations

$$\rho \left(\mathbf{r}\right)\approx \sum _{\left\{\mathbf{G}\right\}}{\rho }_{\mathbf{G}}{e}^{i\mathbf{G}\cdot \mathbf{r}}=\sum _{\left\{\mathbf{G}\right\}}\left|{\rho }_{\mathbf{G}}\right|e^{-i{\varphi }_{\mathbf{G}}+i\mathbf{G}\cdot \mathbf{r}},$$

$${\mathbf{R}}_{1,2}={\left[1-R^{-1}\left(\theta \right)\right]}^{-1}{\mathbf{a}}_{1,2}.$$

$${\varphi }_{{\mathbf{G}}_i}=\tilde{\mathbf{u}}\cdot {\mathbf{G}}_i+\frac{\gamma }{3},$$

$$\tilde{\mathbf{u}}$$

$$\frac{\lambda +\mu }{2}\nabla \left(\nabla \cdot \mathbf{u}\right)+\frac{\mu }{2}{\nabla }^2\mathbf{u}=\frac{\partial }{\partial \mathbf{u}}{\mathcal{V}}_{\text{ad}}\left[\mathbf{r},\mathbf{u}\left(\mathbf{r}\right)\right].$$

$${\mathcal{V}}_{\text{ad}}\left[\mathbf{u}\left(\mathbf{r}\right)\right]=V\sum _{i=1}^3\left\{\frac{1}{2}+\cos \left[{\mathbf{b}}_i\cdot \left(\mathbf{u}+{\mathbf{\Delta }}_0\right)\right]\right\},$$

$$\ell =\sqrt{\frac{\mu }{2\frac{{\partial }^2{\mathcal{V}}_{\text{ad}}}{\partial u^2}{|}_{\text{AB}}}}=\frac{a}{\pi }\sqrt{\frac{\mu }{2V}}\approx 3.2\text{nm}.$$

$$u_3\left(x\right)\approx \frac{2a}{\pi }\arctan \left(e^{\frac{x-x_3}{\ell }}\right),$$

$${\rho }_{\text{1D}}\left(\mathbf{r}\right)$$

$$\approx \sum _{n=0}^{\infty }\left|{\rho }_n\right|\cos \left[n\left({\hat{\mathbf{\varrho }}}_3\cdot {\mathbf{G}}_3\right)\left(x-x_3\right)\right],$$

$$x_i\approx {\hat{\mathbf{\varrho }}}_i\cdot \tilde{\mathbf{u}}-\frac{\gamma L}{6\pi }.$$

$${\partial }_x{\tilde{u}}_x+{\partial }_y{\tilde{u}}_y\sim A_1,$$

$${\partial }_x{\tilde{u}}_y-{\partial }_y{\tilde{u}}_x\sim A_2,$$

$$\displaystyle\left[\begin{array}[]{c}\partial_{x}\tilde{u}_{x}-\partial_{y}\tilde{u}_{y}\\ -\partial_{x}\tilde{u}_{y}-\partial_{y}\tilde{u}_{x}\end{array}\right]\sim E_{2},$$

$$\mathcal{F}\left[\tilde{\mathbf{u}}\left(\mathbf{r}\right)\right]=\frac{1}{2}\int d\mathbf{r}$$

$$\tilde{\lambda }\left(\theta \right)\approx \frac{\left(1+\nu \right)\ell V\sin \left(\frac{\theta }{2}\right)}{a\left(\nu -1\right)},$$

$$\tilde{\mu }\left(\theta \right)\approx \frac{\left(3-\nu \right)\ell V\sin \left(\frac{\theta }{2}\right)}{a\left(1-\nu \right)},$$

$$\gamma \left(\theta \right)\approx \frac{4\ell V\sin \left(\frac{\theta }{2}\right)}{a\left(1-\nu \right)},$$

$${\omega }_{\mathbf{q}}^{(L)}=\sqrt{\frac{\tilde{\lambda }+2\tilde{\mu }}{\tilde{\rho }}}\left|\mathbf{q}\right|\approx \sqrt{\frac{5-3\nu }{4\left(1-\nu \right)}}c\left|\mathbf{q}\right|,$$

$${\omega }_{\mathbf{q}}^{(T)}=\sqrt{\frac{\tilde{\mu }+\gamma }{\tilde{\rho }}}\left|\mathbf{q}\right|\approx \sqrt{\frac{7-\nu }{4\left(1-\nu \right)}}c\left|\mathbf{q}\right|,$$

$$\tilde{\rho }\left(\theta \right)\approx \rho \sqrt{\frac{8V}{{\pi }^2\mu }}\sin \left(\frac{\theta }{2}\right).$$

$$\displaystyle\hat{\mathcal{H}}^{\left(\tau\right)}=\left(\begin{array}[]{cc}\hat{\mathcal{H}}_{D}^{\left(\tau,t\right)}&\hat{T}^{\left(\tau\right)}\left(\mathbf{r}\right)\\ \left[\hat{T}^{\left(\tau\right)}\left(\mathbf{r}\right)\right]^{\dagger}&\hat{\mathcal{H}}_{D}^{\left(\tau,b\right)}\end{array}\right).$$

$${\hat{\mathcal{H}}}_D^{\left(\tau ,\mu \right)}=\hslash v_F{\mathbf{\Sigma }}^{\left(\tau ,\mu \right)}\cdot \left(\hat{\mathbf{k}}-{\mathbf{K}}_{\tau ,\mu }\right),$$

$${\mathbf{\Sigma }}^{\left(\tau ,\mu \right)}=\left(\tau e^{\frac{i\mu \theta {\hat{\ell }}_z}{2}}{\hat{\sigma }}_x{e}^{-\frac{i\mu \theta {\hat{\ell }}_z}{2}},e^{\frac{i\mu \theta {\hat{\ell }}_z}{2}}{\hat{\sigma }}_y{e}^{-\frac{i\mu \theta {\hat{\ell }}_z}{2}}\right),$$

$${\hat{T}}^{\left(\tau \right)}\left(\mathbf{r}\right)=t_{\perp }\left[{\hat{T}}_0^{\left(\tau \right)}+{\hat{T}}_1^{\left(\tau \right)}{e}^{-i\tau {\mathbf{G}}_1\cdot \left(\mathbf{r}-\tilde{\mathbf{u}}\right)}+{\hat{T}}_2^{\left(\tau \right)}{e}^{i\tau {\mathbf{G}}_2\cdot \left(\mathbf{r}-\tilde{\mathbf{u}}\right)}\right],$$

$$\displaystyle\hat{T}_{0}^{\left(\tau\right)}=\left(\begin{array}[]{cc}1&1\\ 1&1\end{array}\right),$$

$$\displaystyle\hat{T}_{1}^{\left(\tau\right)}=\left(\begin{array}[]{cc}1&e^{i\tau\frac{2\pi}{3}}\\ e^{-i\tau\frac{2\pi}{3}}&1\end{array}\right),$$

$$\displaystyle\hat{T}_{2}^{\left(\tau\right)}=\left(\begin{array}[]{cc}1&e^{-i\tau\frac{2\pi}{3}}\\ e^{i\tau\frac{2\pi}{3}}&1\end{array}\right).$$

$$\delta {\hat{T}}^{\left(\tau \right)}\left(\mathbf{r}\right)\approx$$

$${\hat{H}}_{\text{e-ph}}=-\sum _{\tau ,{\mathbf{k}}^{\prime },\mathbf{k},\mathbf{q}}\tau t_{\perp }\{\sqrt{\frac{\hslash }{2A\tilde{\rho }{\omega }_{\mathbf{q}}^{\left(L\right)}}}{\hat{\text{A}}}_{\mathbf{q}}^{\left(L\right)}\left[\frac{\mathbf{q}\cdot {\mathbf{G}}_1}{\left|\mathbf{q}\right|}{\left({\Psi }_{{\mathbf{k}}^{\prime }-\tau {\mathbf{G}}_1}^{\left(\tau ,t\right)}\right)}^{†}{\hat{T}}_1^{\left(\tau \right)}{\Psi }_{\mathbf{k}}^{\left(\tau ,b\right)}-\frac{\mathbf{q}\cdot {\mathbf{G}}_2}{\left|\mathbf{q}\right|}{\left({\Psi }_{{\mathbf{k}}^{\prime }+\tau {\mathbf{G}}_2}^{\left(\tau ,t\right)}\right)}^{†}{\hat{T}}_2^{\left(\tau \right)}{\Psi }_{\mathbf{k}}^{\left(\tau ,b\right)}\right]$$

$$+\sqrt{\frac{\hslash }{2A\tilde{\rho }{\omega }_{\mathbf{q}}^{\left(T\right)}}}{\hat{\text{A}}}_{\mathbf{q}}^{\left(T\right)}\left[\frac{{\left(\mathbf{q}\times {\mathbf{G}}_1\right)}_z}{\left|\mathbf{q}\right|}{\left({\Psi }_{{\mathbf{k}}^{\prime }-\tau {\mathbf{G}}_1}^{\left(\tau ,t\right)}\right)}^{†}{\hat{T}}_1^{\left(\tau \right)}{\Psi }_{\mathbf{k}}^{\left(\tau ,b\right)}-\frac{{\left(\mathbf{q}\times {\mathbf{G}}_2\right)}_z}{\left|\mathbf{q}\right|}{\left({\Psi }_{{\mathbf{k}}^{\prime }+\tau {\mathbf{G}}_2}^{\left(\tau ,t\right)}\right)}^{†}{\hat{T}}_2^{\left(\tau \right)}{\Psi }_{\mathbf{k}}^{\left(\tau ,b\right)}\right]\}{\delta }_{{\mathbf{k}}^{\prime },\mathbf{k}+\mathbf{q}}+\text{h.c.},$$

$${\hat{\mathcal{H}}}^{\left(\tau ,\mu \right)}=\hslash v_F^{\ast }\left(\tau {\hat{\sigma }}_x,{\hat{\sigma }}_y\right)\cdot \mathbf{p}+{\hat{\mathcal{H}}}_{\text{e-ph}}^{\left(\tau ,\mu \right)},$$

$${\hat{\mathcal{H}}}_{\text{e-ph}}^{\left(\tau ,t/b\right)}=$$

$$\pm g_{E_2}^{(a)}\left[\tau \left({\partial }_x{\tilde{u}}_x-{\partial }_y{\tilde{u}}_y\right){\hat{\sigma }}_y+\left({\partial }_x{\tilde{u}}_y+{\partial }_y{\tilde{u}}_x\right){\hat{\sigma }}_x\right],$$

$$⟨{\zeta }^{\prime },\tau ,{\mathbf{k}}^{\prime }|{\hat{H}}_{\text{ep}}|\zeta ,\tau ,\mathbf{k}⟩=\sqrt{\frac{\hslash }{2A\tilde{\rho }{\omega }_{\mathbf{q}}^{\left(T\right)}}}{\varpi }_{T,\tau }^{\zeta ,{\zeta }^{\prime }}\left(\mathbf{q},\mathbf{k},{\mathbf{k}}^{\prime }\right){\delta }_{{\mathbf{k}}^{\prime },\mathbf{k}+\mathbf{q}}.$$

Full text

I Introduction

Two graphene layers rotated with respect to each other an incommensurate angle $\theta$ form a quasi-periodic structure known as moiré pattern.portu At a magic angle $\theta\sim 1^{\textrm{o}}$, quantum interference of electrons in the associated superlattice potential gives rise to narrow low-energy bands,macdonald setting the stage for strongly correlated phenomena.jarillo1 ; jarillo2 In addition to what appears to be Mott insulating states at half-filling of the lowest flat bandsjarillo1 and the onset of superconductivity under dopingjarillo2 or hydrostatic pressure,columbia new many-body insulating states have been reported for different filling factors, some of them with apparent topological character.stanford ; efetov ; UCSB

Qualitative differences in phenomenology from one device to another suggest that we should envision the moiré pattern not as a rigid potential landscape but as a spatially inhomogeneous and (as I will argue here) most likely fluctuating one. Structural inhomogeneities around the magic angle are revealed by variations in the electronic densities of full-filled (single-particle) insulating states, along with evidences from quantum interference measurements of the formation of insulating islands within the superconducting state.jarillo2 ; columbia For smaller twist angles, spatial variations of the moiré period have been directly visualized by transmission electron TEM and scanning tunnel microscopiesLeRoy as well as near-field optical techniques,Basov showing the formation of regions of partial commensuration separated by stacking domain walls.domain_walls_exp ; domain_walls_th These structural differences can induce/favor different symmetry-broken states. The interaction with the encapsulating boron nitride in transport devices can break the sublattice symmetry (usually referred as $\mathcal{C}_{2}\mathcal{T}$ symmetry in the literaturesymmetry1 ; symmetry2 ), opening a gap in the low-energy bands, which acquire nonzero valley-Chern numbers. This state could serve as a precursor for the formation of a (quantum) anomalous Hall ferromagnet at odd fillings as the bands become spin-valley polarized due to electron correlations.stanford ; Zaletel ; Zhang_etal Mean-field calculations suggest that this symmetry can also be spontaneously broken,Xie_MacDonald which has been invoked to explain some features in the magnetotransport of the most homogeneous samples at one-quarter filling of the hole band.efetov The electronic spectrum is also sensitive to the dielectric environment due to electrostatic effects.Paco_PNAS In open samples, suitable for local probes, tunneling spectroscopy has unveiled the emergence of charge ordering with broken 3-fold symmetry as the Fermi level varies.STM1 ; STM2 ; STM3

Another striking observation is the remarkably large, linear in temperature ($T$) resistance of the normal state.jarillo3 ; columbia_phonons Here there appear to be subtle difference between samples too: While in the devices of Cao et al.jarillo3 a linear-$T$ resistivity is only apparent for fractional fillings (more prominently at half-filling), suggesting a connection with the correlated state at lower temperatures, the devices of Polshyn et al.columbia_phonons show qualitatively the same behavior (albeit non-monotonic) for a broader range of carrier densities, with increasing values of the resistivity as the twist angle decreases. This latter behavior has been attributed to electron-phonon scattering.phonons_transport1 ; phonons_transport2 ; phonons_transport3 These studies focus on the original acoustic modes of graphene, but as noticed in Ref. phonons_Bernevig, , relative displacements of the two layers have a stronger impact on the moiré interference pattern and, potentially, on the electronic structure. Here I analyze the case with the account of interlayer interaction forces, which play a fundamental role in the energetics of these modes.

The structure of the article is as follows. I will discuss first in Sec. II the emergence of new collective modes, phasons, associated with the broken translational symmetry of the incommensurate lattice, or more accurately, the invariance of the system under relative translations of the two layers. These modes define two acoustic branches in the spectrum of small oscillations, which is determined by the competition of the two relevant length scales in the problem: the moiré period, and a length scale related to the curvature of the adhesion potential that define the characteristic width of stacking textures (solitons) connecting degenerate adhesion-energy minima. Already for angles $\theta\lesssim 3^{\textrm{o}}$, phonon softening is very pronounced, marking the instability of the system towards the formation of a soliton network. Phason dynamics is then governed by the effective elasticity of these objects. In Sec. III, I will analyze how phasons couple to electrons. Transverse modes (i.e., shear and rotational deformations of the interference pattern) have a stronger impact on the low-energy spectrum and, in particular, can explain the reduced Landau level degeneracy observed in magnetotransport.jarillo1 ; jarillo2 ; columbia Although the coupling considered here is intrinsically weak (proportional to the ratio between interlayer and intralayer hopping parameters), scattering by large (due to the reduced rigidity of stacking solitons) phason fluctuations limit the electron mobilities, leading to a linear-$T$ resistivity that increases with decreasing twist angle. However, this contribution seems to be insufficient to explain the dramatic increase of the resistivity around the magic angle reported in the experiments. I will finally discuss other possible scenarios in Sec. IV. Details of the calculations are saved for the appendices.

II Phasons in moiré-patterned bilayer graphene

Let me first introduce the geometry of the moiré superlattice. In the absence of lattice relaxation, atomic positions in the bottom and top layers are spanned by primitive vectors $\mathbf{a}_{1,2}$ and $\mathbf{a}_{1,2}^{\prime}=R(\theta)\,\mathbf{a}_{1,2}$ of a triangular Bravais lattice, where $R(\theta)$ is a SO(2) matrix describing the relative rotation. For concreteness, I am going to consider the situation depicted in Fig. 1, in which carbon atoms sit on top of each other (AA stacking) prior to the relative rotation. This choice is arbitrary but ultimately inconsequential in the limit of small (generically incommensurate) twist angles, for in that case the system explores all possible stacking configurations. For future reference, lateral positions will be referred to the rotation axis; the coordinate system is such that the $x$ and $y$ axes lie along the two-fold symmetry axes highlighted in Fig. 1 a. Positive twist angles will correspond to anti-clockwise rotations of the top layer.

The moiré superlattice is defined by a beating pattern resulting from the periodicity of the individual triangular lattices, i.e., the Fourier components of the atomic density $\rho(\mathbf{r})$ are peaked at vectors $\mathbf{G}$ of the incommensurate reciprocal lattice,

[TABLE]

where in the last expression I have separated the Fourier components in modulus and phase (note that $\phi_{-\mathbf{G}}=-\phi_{\mathbf{G}}$, so $\rho(\mathbf{r})$ is real). The set $\{\mathbf{G}\}$ corresponds to the lattice spanned by primitive vectors $\mathbf{G}_{1,2}=\mathbf{b}_{1,2}-\mathbf{b}_{1,2}^{\prime}$, where $\mathbf{b}_{1,2}$ and $\mathbf{b}_{1,2}^{\prime}=R(\theta)\,\mathbf{b}_{1,2}$ are primitive vectors of the bottom and top reciprocal lattices, respectively (see Fig. 1 b). The moiré superlattice is just the dual to $\{\mathbf{G}\}$, spanned by primitive vectorsfoot1

[TABLE]

It is also convenient to introduce the function $\bm{\Delta}(\mathbf{r})\equiv\mathbf{r}-R^{-1}(\theta)\,\mathbf{r}$, which measures the distance between a given point $\mathbf{r}$ in the top layer with respect to its original position in the bottom layer before the twist (assuming yet no lattice relaxation). Note that at the maxima (local AA stacking) of the interference pattern $\bm{\Delta}$ coincides with a lattice vector of the graphene Bravais lattice, $\bm{\Delta}(n\,\mathbf{R}_{1}+m\,\mathbf{R}_{2})=n\,\mathbf{a}_{1}+m\,\mathbf{a}_{2}$, with $n$, $m$ integers. The beating pattern maxima are separated by a distance $L_{M}=|\mathbf{R}_{1,2}|=\sqrt{3}a/(2\sin\theta/2)$, defining the moiré superlattice period.

So far, interlayer forces have not appeared in our discussion. Lattice relaxation of structures with a nominal twist angle ($|\rho_{\mathbf{G}}|\neq 0$) could be generically described by a Landau-like free-energy expansion in powers of the beating-pattern density which, as inferred from the previous discussion, gives us an idea of the degree of overlap in the lateral position of the atoms in both layers (and therefore of the interlayer coupling). Minimization of this functional would lead to constraints of the form $\sum_{i=1...n}\phi_{\mathbf{G}_{i}}=\gamma$, with $\sum_{i=1...n}\mathbf{G}_{i}=0$, coming from terms in the $n$th power of $\rho(\mathbf{r})$. The reader should note that such relaxed structures are only metastable states that, in real devices, are likely to be stabilized by the unavoidable tensions generated during the fabrication process. In fact, it is likely that the samples present patches with different twist angles, i.e., different period of the moiré patterning. Here we are concerned only about long-wavelength fluctuations (smooth on the scale of the moiré period $L_{M}$) around a local minimum.

The previous constraints defining the configuration of minimum energy leave a certain number $N$ of the phases $\phi_{\mathbf{G}_{i}}$ unspecified; in other words, variations of $N$ phases generate physically distinct but energetically equivalent quasi-periodic structures and, thus, represent soft modes of the system. In this case, the number of soft modes is $N=2$, for $\{\mathbf{G}\}$ is a two-dimensional Bravais lattice. Consider, for example, a beating pattern truncated to the first six Fourier harmonics (first star), which set the six-fold symmetry and the single independent length scale of the moiré superlattice: $\rho(\mathbf{r})\approx\rho_{0}+\rho_{1}\sum_{i=1,2,3}\cos(\mathbf{G}_{i}\cdot\mathbf{r}-\phi_{\mathbf{G}_{i}})$, with $\mathbf{G}_{3}\equiv-\mathbf{G}_{1}-\mathbf{G}_{2}$. A generic cubic term in the phenomenological free energy fix one of the three phases. A convenient parametrization is then of the form

[TABLE]

where $\mathbf{\tilde{u}}$ is a two-dimensional vector describing the phason modes of the incommensurate lattice.Levine_etal Note that changes in $\mathbf{\tilde{u}}$ translate rigidly the beating pattern, while an arbitrary change in $\gamma$ (which is not a soft mode) distorts the pattern. In a totally incommensurate or floating state, i.e., when lattice relaxation is ignored, $\mathbf{\tilde{u}}$ can be identified straightforwardly with a relative displacement between layers, $\mathbf{u}$. Note that in that case $\bm{\Delta}(\mathbf{r})\rightarrow\mathbf{r}-R^{-1}(\theta)(\mathbf{r}-\mathbf{u})=\bm{\Delta}(\mathbf{r})+R^{-1}(\theta)\,\mathbf{u}$, therefore, the maxima of the beating pattern are translated a distancefoot2

[TABLE]

More generically, phason modes correspond to coherent superpositions of optical phonons with momenta separated by a superlattice vector $\mathbf{G}$. When lattice relaxation is taken into account, these coherent superpositions describe more complex atomic rearrangements than a simple relative displacement of the two layers. The coordinate $\mathbf{\tilde{u}}$ can be identified then with the sliding motion of domain walls separating regions of partial commensuration (alternating AB and energetically equivalent BA Bernal stackings, where atoms of different sublattices lie on top of each other).

II.1 Soliton network

The problem of lattice relaxation can be addressed in the framework of a two-dimensional version of the Frenkel-Kontorova model.rusos ; Mele ; Koshino Neglecting entropic terms in the free energy, the problem reduces to solve the following equation for $\mathbf{u}(\mathbf{r})$, understood now as a field in the continuum describing smooth (in the scale of the interatomic distance $a\approx 1.42$ Å) relative displacements of the layers:

[TABLE]

Here $\mu\approx 3\lambda\approx 9$ eV/Å2 are the Lamé coefficients of graphene;Lame I am disregarding corrugations provided that the bending energy is inconsequential on length scales longer than $\sqrt{\kappa/(\lambda+2\mu)}\sim a$, where $\kappa\approx 1$ eV is the bending rigidity.

The right-hand side of Eq. (5) describes variations in configurational space (stackings) of the adhesion potential (here with units of energy density) between layers. The notation emphasizes the different periodicities in real and configurational spaces: while the dependence on $\mathbf{r}$ is modulated on the scale of the moiré pattern, i.e., it admits a Fourier expansion in $\{\mathbf{G}\}$, the dependence on stacking configurations changes on the atomic unit cell. A first-star expansion of the latter compatible with six-fold symmetry reduces to

[TABLE]

where $\bm{\Delta}_{0}$ is the displacement of the top layer with respect to the bottom layer in Bernal (AB) stacking, which sets the reference in energies. The phenomenological parameter $V\approx 90$ meV/nm2 measures the energy difference between Bernal and AA stackings.Carr_etal The adhesion-energy landscape deduced from this model is shown in Fig. 2 a. In this approximation, the modulation in real space follows from Eq. (6) just by noting that when layers are rotated, the separation between unit cells depends on the position $\mathbf{r}$ via the substitution $\bm{\Delta}_{0}\rightarrow\bm{\Delta}(\mathbf{r})$, and $\mathbf{b}_{i}\cdot\mathbf{\Delta}\left(\mathbf{r}\right)=\mathbf{G}_{i}\cdot\mathbf{r}$ from the previous definitions.

In addition to the periodicity of the moiré pattern, there is another length scale encrypted in Eq. (5) related to the curvature of the adhesion potential,

[TABLE]

This length scale characterizes the spatial extension of stacking textures connecting equivalent minima (AB and BA stackings). The competition between their cost in elastic energy and the adhesion energy of large incommensurate areas governs the degree of lattice relaxation. When $L_{M}$ and $\ell$ are comparable, lattice relaxation is negligible and the beating pattern is well approximated by two floating layers. However, when $L_{M}$ is larger than $\ell$, it is energetically cheaper for the system to form regions of partial commensuration separated by domain walls of characteristic width $\ell$, where the cost in elastic and adhesion energies is concentrated.

The competition between these two length scales is reflected in the the spectrum of small oscillations around a local minimum, see Appendix A. The harmonic expansion of the adhesion potential around a floating state only couples modes with momentum mismatch in the first star, originally separated by a frequency of the order of $\omega_{M}=4\pi c/(\sqrt{3}L_{M})$, where $c=\sqrt{\mu/\rho_{0}}\approx 14$ Km/s is the sound velocity of transverse phonons in graphene. When $L_{M}\sim\ell$, the strength of these harmonics is negligible with respect to $\omega_{M}$, and the spectrum resembles the acoustic phonons of single-layer graphene folded into the moiré Brillouin zone. When $L_{M}$ is ostensibly larger than $\ell$, this simple model predicts a strong softening of the two acoustic branches, indicating that the floating state is unstable and lattice relaxation is not longer negligible. Sharper stacking textures are formed,Koshino giving rise to more harmonics in the expansion of the adhesion potential that have to be included in a new calculation of the spectrum of oscillations. Already in this regime, the simple identification in Eq. (4) breaks down and the energetics of phasons are no longer described by the elasticity of individual layers.

The first step is to determine the stacking texture. In order to construct approximated long-wavelength density profiles, let me consider first one-dimensional solutions of Eq. (5) of the form $\mathbf{u}(\mathbf{r})=u(\varrho)\,\mathbf{\hat{u}}$, where the spatial dependence (along unit vector $\bm{\hat{\varrho}}$) is not necessarily collinear with lattice relaxation. A simple inspection of Eq. (6) (plugging the one-dimensional ansatz and projecting over $\mathbf{\hat{u}}$ and $\mathbf{\hat{z}}\times\mathbf{\hat{u}}$) shows that only tensile ($\mathbf{\hat{u}}\parallel\bm{\hat{\varrho}}$) or shear ($\mathbf{\hat{u}}\perp\bm{\hat{\varrho}}$) solutions exist in this approximation.foot3 I am going to focus on the latter, which are less energetic; this assumption agrees with numerical calculationsKoshino ; Carr_etal (note also that the softening is more pronounced for the transverse oscillation mode) and is obviously favored by the orientation of the moiré superlattice with respect to the graphene lattice in the limit of small twist angles. As implied by the energy landscape in Figs. 2 a, the relaxation is more likely to occur along an armchair direction. Let me focus first on the case $\mathbf{\hat{u}}_{3}\approx\mathbf{\hat{y}}$, $\bm{\hat{\varrho}}_{3}\approx\mathbf{\hat{x}}$ (solitons marked in red in the figure); the other two solutions with the same energy follow from rotating the axes $120^{\textrm{o}}$. The notation emphasizes that $\mathbf{G}_{i}\approx-\frac{2\pi}{L}\bm{\hat{\varrho}}_{i}$ for small twist angles, with $L=|\mathbf{\hat{x}}\cdot\mathbf{R}_{1,2}|\approx\sqrt{3}L_{M}/2$.

In the asymptotic limit $L_{M}\gg\ell$, we can neglect the modulation on the scale of the moiré (first argument in the adhesion potential) and write $\mu\,u^{\prime\prime}=(V\pi/a)\sin(2\pi u/a)$, whose general solution is a train of domain walls separated a certain fixed length. The latter should be determined from energetic considerations,book but in our heuristic construction we can just identify this length scale with $L$. In this limit, the domain-wall profile is well described by a sine-Gordon soliton,

[TABLE]

as numerical calculations confirm.Koshino The corresponding stacking texture is represented in Fig. 2 c. Associated with the domain wall, there is a tension $\sigma=a\sqrt{2\mu V}/\pi$ characterizing its energy cost. This energy does not depend on the soliton center, $x_{3}$: while the separation between domain walls is set by the moiré periodicity, rigid displacements of the soliton solution does not cost energy.

The density wave associated with the one-dimensional train of domain walls is of the form

[TABLE]

where the Fourier components $|\rho_{n}|$ can be estimated from the approximate profile in Eq. (8). The two-dimensional beating pattern can be approximated then by the superposition of three density waves like Eq. (9) with director vectors rotated $120^{\textrm{o}}$. At rigid soliton crossings, the system explores all possible stacking configurations. Adhesion forces will relax the structure, but some areas will remain pinned to saddle points of the potential, including AA stackings, introducing a large free-energy cost. Therefore, the configuration of minimum energy is such that the number of soliton crossings is minimized, and consequently, from the three phases related to the soliton centers, $x_{1,2,3}$, only two are really independent. In the parametrization of Eq. (3), the soliton centers are related to the phason field as

[TABLE]

Minimizing the number of soliton crossings corresponds to the condition $x_{0}+x_{1}+x_{2}\approx 0$ (modulo $L/3$), or equivalently, $\gamma$ must be an integer of $2\pi$.foot4 The resulting beating pattern can be envisioned as a triangular lattice of AA-stacked regions connected by sine-Gordon domain walls, as illustrated in Fig. 2 b.

II.2 Generalized elasticity

Provided that $\ell$ is smooth on the interatomic distance, the free-energy cost of long-wavelength fluctuations of the beating pattern can be described by a phenomenological expansion in terms of derivatives of the phason field, $\partial_{i}\tilde{u}_{j}$, constrained only by the symmetries of the continuum model. Since the superlattice vectors imposes a preferential orientation of the beating pattern (i.e., only translations and not rotations are softfoot_soft ), introducing a symmetric strain tensor is not sensible for this problem. Instead, $\partial_{i}\tilde{u}_{j}$ can be arranged in irreducible representations of $D_{6}$ (see Fig. 1),

[TABLE]

representing compressional, rotational (or tilting), and shear deformations of the soliton network, respectively. A generic harmonic expansion reads then

[TABLE]

where the cost of compressional and shear deformations is expressed in terms of new Lamé coefficients and $\gamma$ is the tilt modulus accounting for the cost of rotations with respect to the moiré superlattice. The elastic constants can be estimated from the approximated profile in Eq. (8) as (see Appendix A)

[TABLE]

where $\nu\approx 0.2$ is graphene’s Poisson ratio.

The soliton network appears to have, within the limitations of the present calculation, a negative Poisson’s ratio ($\tilde{\lambda}<0$).note_bulk This is not so surprising for a two-dimensional incommensurate structure,Halperin_etal and manifests the tendency of the resulting beating pattern to preserve the six-fold symmetry of the moiré superlattice. Here I should emphasize that a relative compression/expansion of the layers (like the uniaxial hetero-strain considered in Ref. fu_strain, ) introduces a shear deformation of the beating pattern, and viceversa, relative shear between layers (which is energetically cheaper) introduces a longitudinal distortion. The latter can be interpreted as a modulation of the moiré period or a nonuniform twist angle, which is systematically observed in topographic images acquired by scanning tunneling microscopy.STM1 ; STM2 ; STM3 The lower energy for longitudinal deformations is ultimately ascribed to the fact that the soliton network connecting AA stacked regions must be understood as a system of strings under tension, whose energy scales linearly (instead of quadratically, like in a system of springs) with length. The tension of the domain walls manifest the metastability of bilayers with a nominal twist angle. These structures can be further stabilized by entropic terms in the free energy, not included in the purely mechanical model discussed here. Two-dimensional soliton configurations carry a lot of entropy,Villain ; book which is also associated with the fact that their energy scales linearly with length. Consequently, thermal fluctuations can contribute to renormalize the elastic constants. Pokrovsky_Tapalov ; Halperin_etal The existence of patches or domains with a nominal twist angle can also be interpreted as the fact that the transition from commensurate to incommensurate structures is first order and dominated by entropy.

Longitudinal and transverse phason modes are represented in Fig. 2 d and e, respectively. Their frequencies deduced from the harmonic expansion in Eq. (12) reads

[TABLE]

where I have introduced in the mass density of the soliton network (see Appendix A):

[TABLE]

The inertia of the soliton system is reduced because the formation of sharper (on the scale of $L_{M}$) stacking textures implies that, effectively, a smaller fraction of the atoms within the moiré cell takes part of the sliding motion of one layer with respect to the other. This compensates the reductions of the elastic constants, so the dispersion of the phason modes do not change much with respect to the acoustic phonons of individual graphene layers, in agreement with numerical calculations.phonons_2013 This does not imply, as we have seen, that interlayer adhesion forces are negligible, or that the elasticity of individual layers can describe the energetics of phason modes.

III Electron-phason coupling

In order to evaluate the effect of phason fluctuations on the electronic spectrum and transport properties, I am going to consider the continuum model usually discussed in the literature:portu ; macdonald

[TABLE]

The Hamiltonian is written in block form, each of them acting on a spinor wave function in a given valley sector (labelled by $\tau=\pm 1$) from top and bottom layers (upper and lower blocks, respectively) describing electronic states around points $K_{\tau,\mu}$ in Fig. 1. The block-diagonal terms are Dirac Hamiltonians of the form

[TABLE]

where $\hat{\mathbf{k}}=-i\bm{\partial}$ is the crystalline-momentum operator in real-space representation, $\hbar v_{F}=3ta/2$, with $t\approx 2.8$ eV being the intralayer hopping parameter, and $\mathbf{\Sigma}^{\left(\tau,\mu\right)}$ is a vector of Pauli matrices defined in the spinor (sublattice) space. Since crystalline momentum is expressed in a common frame of reference (defined by $\mathcal{C}_{2x}$ and $\mathcal{C}_{2y}$ axes in Fig. 1) the Pauli matrices have to be properly rotated,

[TABLE]

where $\hat{\ell}_{z}=\tau\hat{\sigma}_{z}/2$ is the generator of spinor rotations (see Appendix B) and $\mu=\pm 1$ for top/bottom layer blocks. Interlayer tunneling processes are described by the off-diagonal blocks, given by

[TABLE]

where $t_{\perp}\approx 110$ meV parametrizes the strength of the interlayer hopping; matrices $\hat{T}_{i}^{(\tau)}$ contain the pertinent phases acquired by the wave function when electrons hop between different sublattices:

[TABLE]

In addition to time-reversal and superlattice-translation symmetries, along with valley conservation imposed by construction, the continuum model is also invariant under $D_{6}$ point-group operationssymmetry1 ; symmetry2 and infinitesimal rigid translations of one layer with respect to the other: a uniform phason field $\mathbf{\tilde{u}}$ parametrizing the center of the beating pattern can be absorbed in a unitary rotation of the wave function (specifically, it can be absorbed as a phase in the Bloch states defined in different copies of the moiré Brillouin zone), so the mini-band spectrum remains invariant. Note also that the model is derived assuming no lattice relaxation (the main steps are reproduced in Appendix B following Ref. Koshino2, ). Its effect can be incorporated by means of $i$) strain fields within the Dirac-Hamiltonian blocks and $ii$) more harmonics in $\hat{T}(\mathbf{r})$ along with smaller amplitudes for interlayer hoppings between the same sublattice, reflecting the shrinking of AA stacked areas. Although these terms have an important effect on the electronic spectrum,Koshino ; Koshino3 ; Guinea_Walet I am going to neglect them in my estimation of the electron-phason coupling; the observation is that, regardless of the changes in the Hamiltonian, the invariance under relative translations of the two layers imposes the parametric dependence on $\mathbf{\tilde{u}}$ already contained in Eq. (21). Thus, we can obtain the leading contribution in $t_{\perp}$ by expanding Eq. (21) (and its hermitic conjugate) up to linear order in the smooth (on the scale of the moiré period) nonuniform phason field, just like in the case of a charge-density wave,cdw

[TABLE]

Introducing Fourier series for the phason field, decomposing the Fourier components in longitudinal and transverse components, and promoting the latter to boson operators in conventional fashion, provided that $\bm{\pi}=\tilde{\rho}\,\dot{\bm{\tilde{u}}}$ is the canonical momentum density conjugate to soliton displacements, I arrive at the following general expression for the electron-phason coupling in second quantization,

[TABLE]

where $A$ is the area of the system, $\hat{\textrm{A}}_{\mathbf{q}}^{(\nu)}=a_{\mathbf{q}}^{\nu}+(a_{-\mathbf{q}}^{\nu})^{\dagger}$, and $\Psi_{\mathbf{k}}^{(\tau,\mu)}=(c_{A,\mathbf{k}}^{\tau,\mu},c_{B,\mathbf{k}}^{\tau,\mu})^{T}$; operators $(c_{\alpha,\mathbf{k}}^{\tau,\mu})^{\dagger}$/$c_{\alpha,\mathbf{k}}^{\tau,\mu}$ create/annihilate Bloch states with crystalline momentum $\mathbf{k}$ around valley $\tau$ in sublattice $\alpha$ of layer $\mu$, while $(a_{\mathbf{q}}^{\nu})^{\dagger}$/$a_{\mathbf{q}}^{\nu}$ creates/annihilates longitudinal ($\nu=L$) and transverse ($\nu=T$) phasons with momentum $\mathbf{q}$.

A similar coupling was considered in Ref. umklapp, to estimate the contribution from phonon umklapp scattering in graphene on boron nitride. In fact, Eq. (24) captures both normal and umklapp scattering processes when the reconstruction of the electronic spectrum is taken into account and the resulting mini-bands are represented in a reduced-zone scheme, with crystal momentum $\mathbf{k}$ restricted to the first moiré Brillouin zone. For low carrier concentration, umklapp processes are mediated by phonons with momenta of the order of $\mathbf{G}_{i}$ and, therefore, their contribution to the resistivity scales as $\varrho\propto e^{-\hbar\omega_{\mathbf{G}_{i}}/k_{B}T}/T$. In what follows, I am going to restrict the analysis to quasi-elastic scattering processes within the Fermi circles around the mini-Dirac points of the spectrum.

III.1 Low-energy Hamiltonian

Before tacking the calculation of the resistivity, let me discuss first the impact of phason fluctuations on the low-energy sector of the electronic spectrum. In a commensurate approximant, valleys $K_{\tau,\mu}$ are folded onto corners $\tilde{\kappa}_{\eta}$ of the moiré Brillouin zone; Figure 1 illustrates the case for type-I (in the nomenclature of Ref. symmetry2, , sublattice-exchange odd in Ref. Mele0, ) commensurate structures, which are dense in the limit of small angles.portu2 A $\mathbf{k}\cdot\mathbf{p}$ expansion around these points reads

[TABLE]

where $\mathbf{p}$ is the crystalline momentum around $\tilde{\kappa}_{\eta}$. This Hamiltonian acts on a new spinor basis for each valley sector $\tau$, each entry corresponding to envelope wave functions mostly localized on a given sublattice of layer $\mu$; the hybridization with other sublattices/layers is measured by the parameter $\alpha\equiv t_{\perp}/(\hbar v_{F}|\tilde{\kappa}_{\eta}|)$. The new Fermi velocity can be estimated in perturbation theorymacdonald as $v_{F}^{*}/v_{F}=(1-3\alpha^{2})/(1+6\alpha^{2})$; $\alpha^{2}=1/3$ defines the first magic angle in this approximation.

The second term in Eq. (25) represents the electron-phason coupling projected onto the lowest-energy bands. The most general phenomenological expansion allowed by symmetry reads

[TABLE]

where the upper/lower sign applies to top/bottom layer sectors. Note that the electron-phason coupling can only depend on derivatives of the phason field due to the invariance of the electronic spectrum with respect to infinitesimal translations of the moiré pattern. The combinations in Eq. (11) can be paired with electronic operators transforming under the same irreducible representation of $D_{6}$ to form invariants under the point group and time reversal operations; details can be found in Appendix B. The coupling constants $g_{i}$ are phenomenological parameters with units of energy.

The first and third terms in Eq. (III.1) resemble the scalar and vector electron-phonon couplings in graphene and are expected to be subleading in layer hybridization, $g_{A_{1}},g_{E_{2}}^{(s)}\sim\mathcal{O}(\alpha^{2})$; consequently, the electron-phason coupling is dominated by transverse modes. The other two couplings act with opposite signs on valleys $K_{\tau,\pm}$ coming from different layers. In particular, tilting the soliton network ($\bm{\nabla}\times\mathbf{\tilde{u}}\neq 0$) with respect to the preferential direction imposed by the moiré superlattice lifts the degeneracy of these points, as illustrated in Fig. 3 a; the coupling parametrized by $g_{E_{s}}^{(a)}$ includes the effect of distortions of the beating pattern that break the 3-fold rotational symmetry, displacing the positions of the Dirac crossings in $\mathbf{k}$-space. These couplings reproduce the effect of relative strain between layers in the band structurefu_strain and, as noted before, could explain the reduced Landau level degeneracy reported in magnetotransport.

Figure 3 b shows the numerical evaluation (blue dots) of the matrix elements of the coupling with transverse phason modes in the second line of Eq. (24):

[TABLE]

Here $|\zeta,\tau,\mathbf{k}\rangle$ represents electronic states from valley $\tau$ in minin-band $\zeta$ with crystalline momentum $\mathbf{k}$ restricted to the first moiré Brillouin zone, which are obtained by diagonalizing the Hamiltonian in Eq. (18); the resulting matrix is truncated to a finite number of momentum values in each layer coming from different copies of the moiré zone, ranging from $7$ for the largest twist angle ($\theta=5^{\textrm{o}}$) to $61$ for the smallest ($\theta=1.5^{\textrm{o}}$). The results are normalized by the incident momentum (measured with respect to the corresponding mini-Dirac point) along the $x$ axis corresponding to a filling of $n=10^{11}$ cm*-2* in the lowest-energy electron band and fitted (red curves) to the phenomenological expression derived from Eq. (III.1),

[TABLE]

from which I estimate the coupling constants shown in Fig. 3 c for different twist angles. The calculation is restricted to intra-band ($\zeta=\zeta^{\prime}$), quasi-elastic processes, for which $|\mathbf{q}|=|\mathbf{p}^{\prime}-\mathbf{p}|=2|\mathbf{p}|\sin(\frac{\theta_{\mathbf{p}^{\prime}}-\theta_{\mathbf{p}}}{2})$. Note that in these last expressions momenta are measured with respect to $\tilde{\kappa}_{\eta}$, with $\eta=\mu\times\tau$ as prescribed by the folding scheme in Fig. (1), and $\theta_{\mathbf{k}}$ parametrizes the direction of momentum $\mathbf{k}$, namely, $\mathbf{k}=|\mathbf{k}|(\cos\theta_{\mathbf{k}},\sin\theta_{\mathbf{k}})$. In particular, the calculations of Fig. 3 corresponds to $\tau=+1$ and $\mu=+1$; the values of the couplings do not change appreciably when the calculation is performed for different incident momenta or band/valley numbers within the low-energy Dirac cones. The accuracy of the fitting curves is very good down to angles of the order of $\theta=1.2^{\textrm{o}}$.

The layer-symmetric coupling $g_{E_{2}}^{(s)}$ is only appreciable for the smallest angles, as expected. The layer-asymmetric couplings are non-monotonic with twist angle. This behavior is reproduced by the $\mathbf{k}\cdot\mathbf{p}$ perturbative expansion presented in Appendix B, which predicts a stronger coupling with tilting deformations of the beating pattern; second-order perturbation theory gives

[TABLE]

The coupling grows first as the layer hybridization of the electronic wave function $\alpha/(1+6\alpha^{2})$ increases, but around $\theta\lesssim 3^{\textrm{o}}$ it starts to decrease following the same trend as the Fermi velocity.

III.2 Phason-limited electronic transport

Electron scattering by strong phason fluctuations is expected to contribute to the $T$-dependent resistivity of twisted bilayer graphene at small incommensurate angles. I am going to consider a semiclassical treatment in the framework of Boltzmann transport theory, which implicitly assumes that $k_{F}\ell\gg 1$, where $k_{F}$ is the momentum of carriers within the Fermi surface and $\ell$ is the phason-limited mean free path. The electron-phason coupling enters explicitly in the collision integral via a scattering probability rate, which can be computed from Fermi’s golden rule as $\mathcal{W}_{i}^{f}=2\pi\hbar^{-1}|\langle f|\hat{H}_{\textrm{ep}}|i\rangle|^{2}\delta(E_{f}-E_{i})$; the initial state reads $|i\rangle=|\zeta,\tau,\mathbf{k}\rangle\otimes|n_{\mathbf{q}}^{(\nu)}\rangle$, where $|n_{\mathbf{q}}^{(\nu)}\rangle$ represents a state with $n_{\mathbf{q}}^{(\nu)}$ phasons in branch $\nu$. Note that, although the electron-phason coupling projected onto the low-energy bands results from coherent superpositions of electronic states in the two layers, at this point we are neglecting inter-band coherences in the calculation of the resistivity, which could alter the results as the system approaches the neutrality point.coherence1 ; coherence2 I am also going to assume that the phason ensemble thermalize much faster than electrons, so $n_{\mathbf{q}}^{(\nu)}$ reduces to a equilibrium Bose-Einstein distribution function. Phason emission/absorption processes scatterer the initial state into $|f\rangle=|\zeta^{\prime},\tau,\mathbf{k}^{\prime}\rangle\otimes|n_{\mathbf{q}}^{(\nu)}\pm 1\rangle$ at a rate

[TABLE]

In order to obtain analytical formulas, I am going to restrict the analysis to intra-band, quasi-elastic processes within the Dirac cones, dominated by the $g_{A_{2}}$ coupling with transverse phason modes according to the estimates presented in the previous subsection. The resistivity can be obtained from a variational methodZiman applied to the linearized Boltzmann equation describing the evolution of deviations of the electronic distribution function from equilibrium. The calculation is analogous to the case of the phonon-limited resistivity in graphenephonons_bilayer and I am not going to reproduce the details here. The final result reads

[TABLE]

where $I(x)$ is a dimensionless function defined by

[TABLE]

The temperature scale $T_{\textrm{BG}}$, akin to the Bloch-Grüneisen temperature in the problem of electron-phonon scattering, is related to the maximum momentum transfer ($2k_{F}$) between electronic states in a quasi-elastic scattering event,

[TABLE]

At temperatures much lower than $T_{\textrm{BG}}$, the resistivity is dominated by small-angle scattering events and scales as $\varrho\sim(T/T_{\textrm{BG}})^{4}$. In the high-temperature regime, $T\gg T_{\textrm{BG}}$, the resistivity grows linearly with $T$ as

[TABLE]

Given the dispersion relation in Eq. (14b) and the relation between the Fermi momentum and the carrier concentration in the low-energy Dirac model, $k_{F}=\sqrt{\pi n/2}$, the crossover takes place at temperatures of the order of $T_{\textrm{BG}}\approx 12\sqrt{n}$ K, with $n$ measured in units of $10^{11}$ cm*-2*. Figure 4 (black curve) shows the dependence on temperature for fixed values of the carrier density ($n=10^{11}$ cm*-2*) and the dimensionless electron-phason coupling $\bar{g}_{A_{2}}\equiv ag_{A_{2}}/\hbar v_{F}^{*}=0.1$.

IV Discussion

A linear-$T$ resistivity is expected in general if electron transport is limited by scattering off boson fluctuations above a certain temperature scale, the latter defined by phase-space constraints in these scattering events. In particular, for the coupling with either phasons or layer-symmetric phonons (corresponding to acoustic in-phase vibrations of both layers, for which the crossover temperature is of the same order as $T_{\textrm{BG}}$ defined above), the slope of the resistivity with temperature can be written as

[TABLE]

where $R_{0}=\pi\hbar/e^{2}$ is the quantum of resistance, $\mathcal{K}_{i}$ is the suitable elastic modulus ($\tilde{\mu}+\gamma$ for transverse phasons, $\lambda+2\mu$ and $\mu$ for longitudinal and transverse phonons, respectively), and the dimensionless parameter $\bar{g}_{i}$ measures the strength of the coupling within the mini-Dirac bands in units of the bandwidth parameter $\hbar v_{F}^{*}/a$. In the case of phasons, we have $\bar{g}_{A_{2}}\sim t_{\perp}/t$, where the exact fraction depends on the grade of layer hybridization in the wave function and, therefore, is expected to be sensitive to the effect of lattice relaxation on the electronic spectrum, which is not taken into account in the estimates of Sec. III. Regarding layer-symmetric phonons, since the vector potential couples directly to the electron velocity operator, interference of the electronic wave function leads to the same cancelation as in the Fermi velocity, giving $\bar{g}_{v}\approx\beta\equiv-\partial\ln t/\partial\ln a\sim 2-3$. In the case of the scalar coupling, the same cancelation takes place due to the reconstruction of the electronic spectrum, which enters through the electrostatic screening of the deformation potential;Ziman in a Thomas-Fermi treatment,phonons_bilayer the dimensionless coupling reduces to

[TABLE]

where $D\approx 20-30$ eV and $k_{e}=1/4\pi\epsilon_{0}$ are the bare deformation potentialAndo and Coulomb constants, respectively.

The conclusions of this analysis are the following: $i$) The effective couplings (normalized by the bandwidth) within the reconstructed Dirac cones should not depend strongly on the twist angle with the exception of, maybe, the coupling with phason modes due to its sensitivity to layer hybridization. $ii$) Despite the weaker coupling with phasons, these modes can dominate the resistivity at small twist angles due to the reduced stiffness of the stacking domain-wall system. The inset of Fig. 4 shows the dependence of the slope of the resistivity as a function of the twist angle prescribed by the scaling of the elastic constants in Eqs. (13) (i.e., neglecting changes in the effective couplings). The resistivity increases as the twist angle decreases due to the reduction of the rigidity of the soliton network, but the slope is at least two orders of magnitude smaller than those reported in the experiments.phonons_transport1 ; phonons_transport2 This points to a different mechanism, possibly related to strong electron correlations around the magic angle. Recently, González and Stauber have argued that perfect nesting for fillings close to the emergent van Hove singularity at the band edge gives rise to a marginal Fermi liquid scaling of the quasiparticle lifetime.Gonzalez_Stauber1 Partially related to this observation, a recent model for the linear-$T$ resistivity in the normal state of cuprate superconductors invoking umklapp scatteringcuprates starts from the assumption of a Fermi surface reconstruction to maximize commensurate nesting; this could be induced by a spin-wave instability that, also in the present case,phonons_Fu ; Gonzalez_Stauber2 can compete with the superconducting order. Other indications are the fact that the low-temperature crossover is systematically smaller than $T_{\textrm{BG}}$, as pointed out in Ref. phonons_transport2, , and the absence of saturation of the resistivity around the Fermi temperature (only when higher-energy bands start to be populated the resistivity drops), which may indicate that electronic quasiparticles are not well defined.

The conclusions presented here are based on Boltzmann transport theory, which does not include mesoscopic effects ascribed to the intrinsically disordered nature of the devices. In particular, in the comparison with the contribution from conventional phonons, I have neglected the effect of disorder on the beating pattern, which could pin the soliton network and open a gap in the phason spectrum. The reduction of the effective rigidity of the soliton network implies that phason fluctuations are enhanced and, consequently, the moiré pattern is also more sensitive to perturbations induced by the substrate. The simplest perturbation to the harmonic theory in Eq. (12) is a weak (in the sense of LarkinLarkin ) disorder potential of the form

[TABLE]

Here $\mathbf{f}\left(\mathbf{r}\right)$ represent forces acting independently on the stacking solitons. This potential breaks the translational invariance of the incommensurate lattice and pin the soliton system. Positional order in the moiré superlattice is lost at distances of the order of

[TABLE]

where $\sigma_{f}^{2}\equiv\overline{f^{2}}$ represents the dispersion in the distribution of forces. In the case of encapsulated samples, for example, we can estimate this parameter as $\sigma_{f}\sim V_{\textrm{BN}}/L_{\textrm{BN}}$, where $V_{\textrm{BN}}\approx 0.1$ eV/nm2 characterizes the adhesion energy between graphene and boron nitride,Katsnelson_numerics and $L_{\textrm{BN}}\approx 14$ nm is the characteristic size of the moiré due to the lattice mismatch, about $1.8\%$. According to Eqs. (13) and (38), around the magic angle $\theta\sim 1^{\textrm{o}}$ positional order in the corresponding moiré superlattice is lost at distances $L_{c}\sim 26$ nm, about twice the size of the moiré period; $L_{c}$ collapses to $L_{M}$ at angles $\theta\sim 0.7^{\textrm{o}}$.

Equation (38) must be interpreted as a collective pinning lengthLarkin2 ; Lee below which the soliton network responds elastically or, more accurately, the harmonic expansion in Eq. (12) holds. What happens beyond that length it is difficult to say due to the evident shortcomings of the model in Eq. (37), which, for example, neglects the quasi-periodicity of the beating pattern and, possibly, of the interaction with boron nitride if the latter is aligned with the sample. The difficulty arises from the fact that realistic disorder potentials vary on length scales much shorter than the moiré period. These considerations will be taken into account in a future study. Nevertheless, it is worth emphasizing that, although these are only tentative estimates, it is precisely the low stiffness of the soliton network defining the beating pattern at small twist angle what makes the system so sensitive to structural disorder, explaining the widespread presence of twist angle disorder in the samples.

Finally, another possible mechanism for the softening of the soliton network is the presence of strong nematic fluctuations of electronic origin, as suggested by tunneling spectroscopy.STM1 ; STM2 ; STM3 This observation is particularly relevant for the nematicity around the charge neutrality point, which seems to be more pronounced at stacking domain walls.STM2 ; STM3 The phenomenological Hamiltonian presented in Eq. (III.1) provides a good description of the coupling with phason fluctuations close to neutrality. The symmetry analysis presented in Appendix B can be used to construct phenomenological theories describing the coupling of nematic order parameters with lattice degrees of freedom. The interplay between electrons, phasons, and nematic fluctuations is also of potential relevance to understand the $T$-dependent resistivity observed in twisted bilayer graphene. On the experimental front, it would be interesting to have a systematic study of the anisotropy in transport as a function of carrier density and in the absence and presence of applied tensions (although this could be affected by other mesoscopic effects). It would be also very useful to compare the resistivity of devices under different hydrostatic pressure (which have been shown to display a similar phenomenology than magic angle bilayerscolumbia ) in order to determine the dominant role of phasons in transport, since the application of pressure can pin the soliton network and suppress this contribution to the resistivity.

In conclusion, long-wavelength fluctuations of a moiré beating pattern in the limit of small twist angles are dominated by phason modes. Their contribution to resistivity grows linearly with $T$, with increasing slope as the twist angle decreases due to the reduction of the stiffness of the soliton network. This contribution alone, however, seems to be insufficient to explain the fast growth of the resistivity when the magic angle is approached, pointing to a different mechanism that might involve the presence of strong nematic fluctuations or a Fermi surface reconstruction linked to the correlated phenomena at lower temperatures.

Acknowledgements.

I would like to thank Francisco Guinea, Abhay Pasupathy, Matthew Yankowitz, and Cory Dean for valuable discussions, and Ricardo Zarzuela for his careful reading of an early version of this manuscript. This work has been supported by the NSF MRSEC program DMR-1420634.

Note: While finishing this manuscript, I came across a recent preprint,preprint which presents a full calculation of the spectrum of oscillation with the account of lattice relaxation, following the same recipe as in Appendix A. This work highlights the role of stacking domain walls, reaching the same conclusions as in Sec. II, in particular, the scaling of the elastic constants of the moiré superlattice in Eqs. (13).

Appendix A Mechanical model

Lattice relaxation is described by an elastic free energy of the form

[TABLE]

The first terms accounts for in-plane elastic distortions of graphene layers,

[TABLE]

where $\mathbf{u}^{(\mu)}$ describes displacements of unit cells in layer $\mu$ with respect to their equilibrium positions in the absence of interlayer forces. It is convenient to introduce the relative displacement, $\mathbf{u}=\mathbf{u}^{(t)}-\mathbf{u}^{(b)}$, and total displacement, $\mathbf{v}=\mathbf{u}^{(t)}+\mathbf{u}^{(b)}$. In the simplest approximation, the adhesion energy, second term in Eq. (39), is a functional of the former field only,

[TABLE]

where $\mathcal{V}_{\textrm{ad}}$ is the adhesion potential introduced in the main text. Equation (5) comes from minimizing the free energy of variations with respect to $\mathbf{u}$. Dynamical equations are derived from the total Lagrangian $L=K-F$, where the kinetic energy reads

[TABLE]

Here $\rho=7.6\times 10^{-7}$ kg/m2 is the mass density of individual graphene layers.

The first field $\mathbf{v}$ represents in-phase displacements of both layers and can be identified with the original acoustic phonons.phonons_bilayer From this point on, I am going to focus on the dynamics of relative displacements, $\mathbf{u}$. Let me consider deviations from a metastable configuration, $\delta\mathbf{u}(t,\mathbf{r})=\mathbf{u}(t,\mathbf{r})-\mathbf{u}^{(0)}(\mathbf{r})$, where $\mathbf{u}^{(0)}(\mathbf{r})$ is a solution of Eq. (5). Plugging this ansatz into the previous equations gives

[TABLE]

where the first and second terms are just Eq. (40) evaluated with the metastable solution and the corresponding deviation, while the third term mixes both,

[TABLE]

The last term comes from the adhesion energy; expanding up to quadratic order gives

[TABLE]

Hereafter repeated latin indices are summed up. Integration by parts (dropping boundary terms) leads to

[TABLE]

The last term is [math] just from Eq. (5). The first term represents the free energy of the equilibrium solution,

[TABLE]

while $U[\delta\mathbf{u}]$ describes the spectrum of harmonic oscillations,

[TABLE]

The Euler-Lagrange equations describing oscillations around the metastable state reads then

[TABLE]

By introducing Fourier series,

[TABLE]

the problem reduces to solve the following secular equation,

[TABLE]

where I have introduced

[TABLE]

A.1 Floating layers

For large twist angles such that $L_{M}$ and $\ell$ (defined in Eq. 7) are comparable, we can neglect lattice relaxation, $\mathbf{u}^{(0)}\approx 0$, and consider oscillations around two floating layers. The average free energy per moiré cell is just $\bar{F}_{0}=3V/2$. The adhesion potential reduces to

[TABLE]

In this approximation, all the harmonics $K_{ij}(\mathbf{G})$ are [math] except for the ones in the first star ($\alpha=1,2,3$)

[TABLE]

The spectra in Fig. 5 are the result of truncating the secular equation to include the first 36 harmonics, giving a $37\times 37$ matrix to diagonalize. The lowest-frequency branches are strongly softened (blue curves) when the matrix elements in Eq. (54) start to be comparable with $\omega_{M}=4\pi c/(\sqrt{3}L_{M})$.

A.2 Soliton network

Let me assume now that the lattice relaxes and the solution is roughly described by the superposition of the three density waves introduced in the main text,

[TABLE]

where $\mathbf{u}_{\alpha}\left(\varrho\right)\approx u\left(\varrho\right)\,\mathbf{\hat{u}}_{\alpha}$ is a soliton train along $\bm{\hat{\varrho}}_{\alpha}=\mathbf{\hat{z}}\times\mathbf{\hat{u}}_{\alpha}$, $u(\varrho)$ is given by Eq. (8), and the unit vectors $\mathbf{\hat{u}}_{\alpha}$ lie along the three armchair directions:

[TABLE]

To the leading order in $\ell/L_{M}$, we can neglect soliton interactions at the crossings and approximate

[TABLE]

For a sine-Gordon soliton, the tension reduces to

[TABLE]

As argued in the main text, the acoustic branches in this limit can be identified with traveling-wave modes of the soliton network parametrized by collective coordinates $\mathbf{\tilde{u}}$. Specifically, if $\mathbf{u}^{(0)}\left(\mathbf{r}-\mathbf{\tilde{u}}\right)$ is the soliton-network solution centered at $\mathbf{\tilde{u}}$, then a smoothly distorted profile can be approximated by the functional

[TABLE]

where the collective coordinate has been promoted to a field $\mathbf{\tilde{u}}(\mathbf{r})$. Spatial derivatives can be approximated as

[TABLE]

and therefore we can identify $\partial_{i}\delta u_{j}\approx-\partial_{k}u_{j}^{(0)}\partial_{i}\tilde{u}_{k}$. Plugging this result in the first line of Eq. (48) gives

[TABLE]

We can now estimate from this last result the coefficients in the free-energy expansion of Eq. (12). Since the spatial dependence of the phason field $\mathbf{\tilde{u}}(\mathbf{r})$ must be smooth on the scale of the moiré superlattice, the derivatives of $\mathbf{u}^{(0)}$ inside the integral can be approximated by their average over a moiré unit cell. We end up with

[TABLE]

where

[TABLE]

Note that the symmetry-adapted expansion in the main text adopts the same form, where the tensor of elastic coefficients reads

[TABLE]

In order to evaluate the elastic constants, we can proceed as before and plug the superposition of the three soliton waves into Eq. (63). Ignoring again the subleading contribution from soliton crossings, I find

[TABLE]

where the prefactor reads

[TABLE]

The tensors between parenthesis are just

[TABLE]

Comparing these expressions with Eq. (64), I arrive at the final formulas in Eqs. (13) of the main text.

The effective mass density of the soliton network can be derived in the same manner by approximating the time derivatives in Eq. (42) as $\dot{u}_{i}\approx-\partial_{j}u_{i}^{(0)}\dot{\tilde{u}}_{j}$. The mass tensor of a soliton reads then

[TABLE]

where the integral is extended over a moiré unit cell. Proceeding just as before, the final result reads $M_{ij}=M\delta_{ij}$, where $M$ is the inertia of a stacking domain wall,

[TABLE]

Dividing this quantity by the area of the moiré unit cell $A_{M}=\sqrt{3}L_{M}^{2}/2$ gives the mass density $\tilde{\rho}$ in Eq. (15). The dispersion relations in Eqs. (14) follow from the corresponding Euler-Lagrange equations, where the Fourier components of the phason field can be decomposed in longitudinal and transverse components as usual,

[TABLE]

These can be promoted to phason creation/annihilation operators in conventional fashion via the identification (here $\nu=L,T$ labels the branch)

[TABLE]

where phason operators satisfy the boson algebra $[a_{\mathbf{q}_{1}},(a_{\mathbf{q}_{2}})^{\dagger}]=\delta_{\mathbf{q}_{1},\mathbf{q}_{2}}$, following from the conjugacy relations $\{\tilde{u}_{i}(\mathbf{r}),\pi_{j}(\mathbf{r}^{\prime})\}=\delta_{ij}\,\delta(\mathbf{r}-\mathbf{r}^{\prime})$, where $\bm{\pi}=\tilde{\rho}\,\dot{\bm{\tilde{u}}}$.

Appendix B Electronic model

Let me start by writing a generic tight-binding Hamiltonian for two graphene layers floating on top of each other, $\hat{H}_{\textrm{tb}}=\hat{H}_{b}+\hat{H}_{t}+\hat{T}+\hat{T}^{\dagger}$. The first two terms represent intralayer hopping processes; up to first nearest neighbors, we have

[TABLE]

The Hamiltonian is written in a monoelectronic basis of Wannier $\pi$-orbitals localized on sites $A/B$ of unit cell $i$ of layer $\mu$. We can introduce Bloch states as

[TABLE]

with crystalline momenta $\mathbf{k}_{\mu}$ restricted to the first Brillouin zone of layer $\mu$ ($N$ is the number of unit cells on each layer that I assume the same). The previous Hamiltonian reduces to

[TABLE]

where the matrix in sublattice space reads

[TABLE]

The sum is extended over the three vectors connecting $A$ with nearest B sites in layer $\mu$; in my notation, $\bm{\delta}_{1,2}^{(\mu)}=\bm{\delta}_{3}^{(\mu)}-\mathbf{a}_{1,2}^{(\mu)}$, and $\bm{\delta}_{3}^{(\mu)}=\mathbf{R}_{i,B}^{(\mu)}-\mathbf{R}_{i,A}^{(\mu)}$ connects the two sites within the unit cell.

The terms $\hat{T}$ ($\hat{T}^{\dagger}$) describe interlayer tunneling processes,

[TABLE]

In a two-center, Slater-Koster-like approximation, $T_{\alpha\beta}^{ij}$ depends only on the relative distance between Wannier centers, so it must admit a Fourier expansion of the form

[TABLE]

where $A_{c}$ is graphene’s unit-cell area. Introducing Bloch states and plugging this last expression into Eq. (78), we can rewrite the latter as

[TABLE]

where the first two sums are on reciprocal vectors of the top and bottom lattices. In deriving this expression, I have made use of the identities

[TABLE]

where $\bm{\delta}_{\alpha}^{(\mu)}$ is the position of site $\alpha$ within a reference unit cell in layer $\mu$, such that $\bm{\delta}_{B}^{(\mu)}-\bm{\delta}_{A}^{(\mu)}=\bm{\delta}_{3}^{(\mu)}$.

Equation (80) describes scattering events satisfying the general umklapp condition

[TABLE]

where $n_{i}$, $n_{i}^{\prime}$ are integers. Here I have introduced the notation of the main text, namely, $\mathbf{b}_{i}$ and $\mathbf{b}_{i}^{\prime}=R(\theta)\mathbf{b}_{i}$ are the primitive vectors of the reciprocal lattice of the bottom and top layers, respectively; similarly, I write $\bm{\delta}_{\alpha}^{(b)}\equiv\bm{\delta}_{\alpha}$, $\bm{\delta}_{\alpha}^{(t)}=R(\theta)\bm{\delta}_{\alpha}+\mathbf{u}$ where, in addition to the relative rotation along a common hexagon center, I consider a relative displacement $\mathbf{u}$ of the top with respect to the bottom layer. Commensurate approximants to the electronic structure simplifies the condition in Eq. (82) such that $n_{i}=n_{i}^{\prime}$, i.e., $\mathbf{k}_{t}=\mathbf{k}_{b}+\mathbf{G}$, where $\mathbf{G}$ is a vector of the moiré superlattice; equation (80) simplifies then to

[TABLE]

The sum in $\left\{\mathbf{b}\right\}$ is also restricted in practice, provided that $T(\mathbf{q})$ is a rapidly decaying function of momentum; in particular, $T(\mathbf{q})$ is strongly suppressed for values $|\mathbf{q}|>1/d$, where $d$ is the separation between layers.portu ; macdonald

For small twist angles, the electronic spectrum will be dominated by the low-energy Dirac bands of decoupled graphene layers lying around the two inequivalent corners (valleys) of the respective Brillouin zones labelled by $K_{\tau,\mu}$ in Fig. 1; here $\tau=\pm 1$ labels the valleys, which is assumed to be a good quantum number. In a continuum description,portu ; macdonald band dispersion is simplified by expanding the phases in Eq. (77) for small momenta around $\mathbf{K}_{\tau,\mu}$,

[TABLE]

where the $x$, $y$ components are adapted to the high-symmetry (zig-zag and armchair, respectively) axes of individual graphene layers. As mentioned in the main text, it is convenient to expresses crystalline momenta in a common frame of reference (defined in this case by the high-symmetry axes highlighted in Fig. 1), so the spinor basis in the corresponding sublattice space of each layer has to be rotated, leading to the block-diagonal terms in Eq. (18).

For the tunneling terms, we can, in the same spirit, neglect small deviations from $\mathbf{K}_{\tau,\mu}$ in the argument of $T(\mathbf{q})$. If we include only one harmonic, $t_{\perp}\equiv T\left(|\mathbf{K}_{\tau,\mu}|\right)$, the sum in $\{\mathbf{b}\}$ has to be restricted to $\mathbf{b}=\mathbf{0},-\tau\,\mathbf{b}_{1},\tau\,\mathbf{b}_{2}$, corresponding to the three equivalent positions of a given valley within a single-layer Brillouin zone. Equation (83) simplifies to

[TABLE]

where I have used the relation in Eq. (4) along with the identity in footnote foot1, , so that $\mathbf{b}\cdot(R^{T}\mathbf{u})=-\mathbf{G}\cdot\mathbf{\tilde{u}}$. Equation (85) transformed back to real space corresponds to Eq. (21), with matrices $\hat{T}_{i}$ given in Eqs. (22).

B.1 Symmetry-adapted electronic operators

In order to construct phenomenological couplings within the low-energy sector of the spectrum, it is convenient to introduce a basis of electronic operators adapted to the irreducible representations (irreps) of the point-group symmetry $D_{6}$ of the continuum model. The internal Hilbert space of the continuum model is spanned by spin, sublattice, layer, and valley degrees of freedom. I am going to ignore the spin since relativistic corrections are weak. Note that neither sublattice nor layer are good quantum numbers, since inter-layer hopping terms mix them, but we can still refer to them to label the transformation properties of the Bloch wave function around $\tilde{\kappa}_{\eta}$ points. Valley $\tau=\pm 1$ and layer $\mu=\pm 1$ indices label four Dirac crossings reminiscent of the Dirac points of two decoupled graphene layers, while sublattice indices span the associated subspace of Dirac doublets at each crossing. Electronic operators can be written in a basis of Pauli matrices acting on each subspace: $\hat{\sigma}_{i}$ for sublattice, $\hat{\tau}_{i}$ for valley, and $\hat{\mu}_{i}$ for layer. Since we are interested in the coupling to long-wavelength phason fluctuations on the scale of the moiré period, I am going to restrict the analysis to diagonal operators in valley and layer (or $\tilde{\kappa}_{\eta}$-point) numbers with the pertinent identifications $\hat{\tau}_{z}\rightarrow\tau$, $\hat{\mu}_{z}\rightarrow\mu$. I find convenient, however, to keep the matrix notation here, because the transformation rules of electronic operators are more easily identified from the Pauli-matrix algebra.

For example, elementary rotations along the 6-fold principal axis are implemented by unitary operators

[TABLE]

where $\hat{\ell}_{z}=\frac{1}{2}\hat{\tau}_{z}\hat{\sigma}_{z}$. Two-fold rotations along $x$ and $y$ axes are given by

[TABLE]

The rest of point-group operations follow from matrix multiplication. Finally, time-reversal symmetry is implemented by the anti-unitary operator

[TABLE]

where $\mathcal{K}$ denotes complex conjugation.

All the possible combinations of valley- and layer-diagonal operators can be classified according to the irreps of $D_{6}$ and time-reversal symmetry. The result is summarized in Tab. 1, along with the characters of $D_{6}$. It is worth emphasizing that there are two ways to lift the Dirac degeneracies without breaking time-reversal symmetry: by breaking both $\mathcal{C}_{2z}$ and $\mathcal{C}_{2y}$ symmetries ($\hat{\mu}_{z}\hat{\sigma}_{z}\sim B_{1}$), i.e., by removing the physical equivalence of both sublattices and layers (a staggered potential of opposite sign in each layer), or by breaking $\mathcal{C}_{2z}$ and $\mathcal{C}_{2x}$ symmetries ($\hat{\sigma}_{z}\sim B_{2}$), i.e., by removing only the physical equivalence between sublattices (the same staggered potential in both layers). Only the latter leads to bands with nonzero valley-Chern number. The competition between these two mass terms give rise to the phase diagram discussed in Ref. Zaletel, .

The phenomenological electron-phason Hamiltonian in Eq. (III.1) consists of all the possible invariants formed from combinations of the phason-field derivatives in Eqs. (11) with the corresponding electronic operators transforming under the same irrep. Here I should note that the phason field transforms as $\mathbf{\tilde{u}}=(\tilde{u}_{x},\tilde{u}_{y})\sim E_{1}$, while the relative displacement between layers transforms as $(-u_{y},u_{x})\sim E_{1}$ since, recall, $\mathcal{C}_{2x,y}$ rotations exchange the layers. The relation between these two transformation rules can be understood from the twist (second, dominant term) in Eq. (4) or, with the account of lattice relaxation, the fact that solitons involve relative shear between layers.

B.2 Perturbative calculation of the electron-phason matrix element

The Hamiltonian of the continuum model can be diagonalized in a basis of Bloch states by restricting the values of crystalline momentum to the first Brillouin zone of the moiré reciprocal lattice and introducing new (band) quantum numbers associated with different copies separated by momenta in $\{\mathbf{G}\}$. In this process, we can absorb a uniform phason field as a phase in the new electronic basis, so the spectrum remains invariant under translations of the moiré pattern. It is useful to consider the positions of the valleys folded into the first moiré Brillouin zone,

[TABLE]

This folding scheme applies to type-Isymmetry2 or sublattice-exchange oddMele0 commensurate approximants, but once the first sum in Eq. (83) is restricted to one harmonic, this choice is inconsequential.portu2

Next, we consider the dispersion of low-energy electronic states with momentum $\mathbf{p}$ around these points, in the previous notation, $\mathbf{k}_{t}=\mathbf{K}_{\tau,t}+\mathbf{p}$ (defining the top layer sector, $\mu=+1$), and $\mathbf{k}_{b}=\mathbf{K}_{\tau,b}+\mathbf{p}$ (corresponding to the bottom layer sector, $\mu=-1$). For each moiré Brillouin zone in a given layer I am going to truncate the number of copies in the opposite layer to three, so each (decoupled) sector labelled by valley and layer numbers is described by a $8\times 8$ matrix Hamiltonian of the form

[TABLE]

The first term reads just (in block form)

[TABLE]

for the top-layer sector, and similarly for the bottom layer,

[TABLE]

From this point on, I am going to neglect the rotation of the spinor basis, which restores the electron-hole symmetry of the spectrum (the error scales with $\theta^{2}$); while these matrices read the same now, note that they are expressed in a different basis. The second term in Eq. (90) reads

[TABLE]

for the top-layer sector, and

[TABLE]

for the bottom-layer sector. The low-energy subspace is defined by the zero-energy eigenstates of $\hat{\mathcal{U}}^{\left(\tau,\mu\right)}$, specifically,

[TABLE]

These states form a basis for the $D_{6}$ irreps introduced in the previous subsection. The rest of eigenstates are separated by energies $E_{i}=\pm\hbar v_{F}|\tilde{\kappa}|$, $\pm\hbar v_{F}\sqrt{1+6\alpha^{2}}\,|\tilde{\kappa}|$ (the last ones are 2-fold degenerate in this approximation).

Next, I am going to produce a $\mathbf{k}\cdot\mathbf{p}$ expansion of the Hamiltonian by projecting-out the high-energy states in a Löwdin perturbative scheme.Lowdin Let me introduce first the following projection operators:

[TABLE]

where $i$ labels the high-energy eigenstates. Up to second order in perturbation theory, where any term in the Hamiltonian of Eq. (90) (generically written as $\hat{V}^{(\tau,\mu)}$) is treated as a perturbation to $\hat{\mathcal{U}}^{\left(\tau,\mu\right)}$, we have

[TABLE]

First-order perturbation theory in $\hat{\mathcal{H}}_{\mathbf{p}}^{\left(\tau,\mu\right)}$ leads to the first term in Eq. (25) with a new Fermi velocity $v_{F}^{*}$ reduced by a factor $(1-3\alpha^{2})/(1+6\alpha^{2})$, the result first obtained by Bistritzer and MacDonald in Ref. macdonald, . For the electron-phason coupling, we should consider first its expression in the $8\times 8$ Hilbert space; from Eq. (23), we have

[TABLE]

for the top-layer sector, and

[TABLE]

for the bottom layer. First-order perturbation theory gives [math], as expected, since the invariance under translations of the soliton network implies that there must be momentum transfer between electronic states.foot_robustness Treating both $\hat{V}_{\textrm{e-ph}}^{\left(\tau,\mu\right)}$ and $\hat{\mathcal{H}}_{\mathbf{p}}^{\left(\tau,\mu\right)}$ up to second order in perturbation theory gives

[TABLE]

where $\mathbf{p}$ must be interpreted now as the average momentum during a electron-phason scattering event, $\mathbf{p}\equiv(\mathbf{k}^{\prime}+\mathbf{k})/2$, with $\mathbf{k}$, $\mathbf{k}^{\prime}$ labelling the initial and final states, respectively.

Equation (B.2) may look odd at first glance, for it is difficult to recognize the symmetry-allowed couplings in the phenomenological expansion of Eq. (III.1) from this expression. Note, however, that Eq. (B.2) is compatible with time-reversal and $D_{6}$ point-group symmetries and, in fact, the combination of valley, sublattice and layer indices is such that when projected over a band state, i.e., and eigenstate of the first term in Eq. (25),

[TABLE]

where $\zeta=\pm 1$ labels the electron/hole bands, Eq. (B.2) produces the expected matrix element of $g_{A_{2}}$ coupling. Note first that for intra-band processes we have

[TABLE]

When the phason field is decomposed in longitudinal and transverse components, the former give rise to scalar products of the form $i\left(\mathbf{k}+\mathbf{k}^{\prime}\right)\cdot\mathbf{\hat{q}}$, while for transverse modes we have $i\left(\mathbf{k}+\mathbf{k}^{\prime}\right)\cdot\left(\mathbf{\hat{z}}\times\mathbf{\hat{q}}\right)$; here $\mathbf{\hat{q}}$ is a unit vector along the transferred momentum, $\mathbf{q}=\mathbf{k}^{\prime}-\mathbf{k}$. For quasi-elastic processes ($|\mathbf{k}^{\prime}|=|\mathbf{k}|$) illustrated in Fig. 6, we have

[TABLE]

Collecting all the pieces, I arrive at the following expression for the matrix element of the coupling with transverse phasons,

[TABLE]

where I have regrouped some factors in $v_{F}^{*}/v_{F}$. This is indeed the matrix element of the coupling with $(\bm{\nabla}\times\mathbf{\tilde{u}})_{z}$ since for quasi-elastic processes $|\mathbf{q}|=2|\mathbf{k}|\sin(\frac{\theta_{\mathbf{k}^{\prime}}-\theta_{\mathbf{k}}}{2})$. We end up then with the identification of the phenomenological parameter in Eq. (29).