Is magic-angle twisted bilayer graphene near a quantum critical point?

TL;DR

This study uses quantum Monte Carlo simulations on an extended Hubbard model to identify a quantum critical point related to symmetry breaking, which may shed light on phenomena in magic-angle twisted bilayer graphene.

Contribution

It provides the first unbiased numerical evidence of a quantum critical point in a model relevant to twisted bilayer graphene, with precise critical exponents and phase transition details.

Findings

Identification of a quantum phase transition from Dirac semi-metal to VBS insulator.

Critical exponents consistent with 3D N=4 Gross-Neveu chiral XY universality.

Discovery of two distinct VBS phases and their transition properties.

Abstract

An extended Hubbard model on a honeycomb lattice with two orbitals per site at charge neutrality is investigated with unbiased large-scale quantum Monte Carlo simulations. The Fermi velocity of the Dirac fermions is renormalized as the cluster charge interaction increases, until a mass term emerges and a quantum phase transition from Dirac semi-metal to valence bond solid (VBS) insulator is established. The quantum critical point is discovered to belong to 3D $N=4$ Gross-Neveu chiral XY universality with the critical exponents obtained at high precision. Further enhancement of the interaction drives the system into two different VBS phases, the properties and transition between them are also revealed. Since the model is related to magic-angle twisted bilayer graphene, our results may have relevance towards the symmetry breaking order at the charge neutrality point of the material, and associate the wide range of universal strange metal behavior around it with quantum critical fluctuations.

Tables

Table 1. Table S1: Comparison of critcal exponents for 3D N = 4 𝑁 4 N=4 Gross-Nevue chiral XY universality class calculated by different methods.

$N=4$	$\eta$	$\nu$
Monte Carlo, this work	0.80(2)	1.01(3)
$4-ϵ$, four loop, $P_{\left[2/2\right]}$ Zerf et al. (2017b)	0.929	1.130
$4-ϵ$, four loop, $P_{\left[3/1\right]}$ Zerf et al. (2017b)	0.911	1.130
functional RG (LPA’) Classen et al. (2017)	0.946	1.082
$4-ϵ$, two loop Rosenstein et al. (1993)	0.82	0.97
Monte Carlo Li et al. (2017)	0.80(4)	1.11(3)
large-$N$ Hands et al. (1993)	0.932	1.135

Equations

$$H_0=-t\sum _{⟨ij⟩l\sigma }\left(c_{il\sigma }^{†}{c}_{jl\sigma }+h.c.\right)-t_2\sum _{⟨ij{⟩}^{\prime }l\sigma }\left(i^{2l-1}{c}_{il\sigma }^{†}{c}_{jl\sigma }+h.c.\right)$$

$$H_U=U\sum _{\text{varhexagon}}(Q_{\text{varhexagon}}-4{)}^2,$$

$$C_B(\mathbf{k})=\frac{1}{L^4}\sum _{i,j}{e}^{i\mathbf{k}\cdot ({\mathbf{r}}_i-{\mathbf{r}}_j)}⟨B_{i,\delta }{B}_{j,\delta }⟩$$

$$H_{\text{eff}}(\mathbf{k})=-t\left(\tilde{\mathbf{k}}\cdot \tilde{\mathbf{s}}{\tau }^2+m{s}^0{\tau }^3\right),$$

$$H_0=-t\sum _{⟨ij⟩l\sigma }\left(c_{il\sigma }^{†}{c}_{jl\sigma }+h.c.\right)-t_2\sum _{⟨ij{⟩}^{\prime }l\sigma }\left(i^{2l-1}{c}_{il\sigma }^{†}{c}_{jl\sigma }+h.c.\right).$$

$$⟨\hat{O}⟩=\frac{⟨{\Psi }_0|\hat{O}|{\Psi }_0⟩}{⟨{\Psi }_0|{\Psi }_0⟩}=\underset{\Theta \to \infty }{\lim }\frac{⟨{\Psi }_T|e^{-\frac{\Theta }{2}\mathbf{H}}\hat{O}{e}^{-\frac{\Theta }{2}\mathbf{H}}|{\Psi }_T⟩}{⟨{\Psi }_T|e^{-\Theta \mathbf{H}}|{\Psi }_T⟩}.$$

$$⟨{\Psi }_T|e^{-\Theta H}|{\Psi }_T⟩=⟨{\Psi }_T|{\left(e^{-\Delta \tau H_U}{e}^{-\Delta \tau H_0}\right)}^{L_{\tau }}|{\Psi }_T⟩+\mathcal{O}(\Delta {\tau }^2)$$

$$e^{-\Delta \tau U(Q_{\text{varhexagon}}-4{)}^2}=\frac{1}{4}\sum _{\{s_{\text{varhexagon}}\}}\gamma (s_{\text{varhexagon}})e^{\alpha \eta (s_{\text{varhexagon}})\left(Q_{\text{varhexagon}}-4\right)}$$

$$⟨{\Psi }_T|e^{-\Theta H}|{\Psi }_T⟩=\sum _{\{s_{\text{varhexagon},\tau }\}}\left[\left(\prod _{\tau }\prod _{\text{varhexagon}}\gamma (s_{\text{varhexagon},\tau })e^{-4\alpha \eta (s_{\text{varhexagon},\tau })}\right)\det \left[P^{†}B(\Theta ,0)P\right]\right]$$

$$B(\tau +1,\tau )=e^{V[\{s_{\text{varhexagon},\tau }\}]}{e}^{-{\Delta }_{\tau }K}$$

$$G(\tau ,\tau )=1-R(\tau ){\left(L(\tau )R(\tau )\right)}^{-1}L(\tau )$$

$$\mathcal{U}=i{\sigma }_0{\tau }_2\mathcal{K},$$

$$h(k)=-\tilde{t}\left(\begin{array}[]{cccccc}0&(1+m)e^{-ik_{y}}&0&e^{-i\frac{\sqrt{3}}{2}k_{x}+\frac{i}{2}k_{y}}&0&e^{i\frac{\sqrt{3}}{2}k_{x}+\frac{i}{2}k_{y}}\\ (1+m)e^{ik_{y}}&0&e^{i\frac{\sqrt{3}}{2}k_{x}-\frac{i}{2}k_{y}}&0&e^{-i\frac{\sqrt{3}}{2}k_{x}-\frac{i}{2}k_{y}}&0\\ 0&e^{-i\frac{\sqrt{3}}{2}k_{x}+\frac{i}{2}k_{y}}&0&(1+m)e^{i\frac{\sqrt{3}}{2}k_{x}+\frac{i}{2}k_{y}}&0&e^{-ik_{y}}\\ e^{i\frac{\sqrt{3}}{2}k_{x}-\frac{i}{2}k_{y}}&0&(1+m)e^{-i\frac{\sqrt{3}}{2}k_{x}-\frac{i}{2}k_{y}}&0&e^{ik_{y}}&0\\ 0&e^{i\frac{\sqrt{3}}{2}k_{x}+\frac{i}{2}k_{y}}&0&e^{-ik_{y}}&0&(1+m)e^{-i\frac{\sqrt{3}}{2}k_{x}+\frac{i}{2}k_{y}}\\ e^{-i\frac{\sqrt{3}}{2}k_{x}-\frac{i}{2}k_{y}}&0&e^{ik_{y}}&0&(1+m)e^{i\frac{\sqrt{3}}{2}k_{x}-\frac{i}{2}k_{y}}&0\end{array}\right)$$

$$U=\left(\begin{array}[]{ccc}-1&-1&1\\ 0&1&1\\ 1&0&1\end{array}\right)\otimes\tau^{0}.$$

$$H_{\text{eff}}(\mathbf{k})=-t\left(\tilde{\mathbf{k}}\cdot \tilde{\mathbf{s}}{\tau }^2+m{s}^0{\tau }^3\right),$$

Full text

Introduction —  Twisted bilayer graphene (TBG) forms Morié patterns in real space with the size of the Morié unit cell tuned by the twisting angle. The Fermi velocity of the Dirac fermions of monolayer graphene is renormalized in TBG. At some magic angle, the Fermi velocity vanishes Bistritzer and MacDonald (2011); Suárez Morell et al. (2010); Lopes dos Santos et al. (2012); Fang and Kaxiras (2016); Trambly de Laissardière et al. (2012); Tarnopolsky et al. (2018), such that flat bands are formed and the system consequently become susceptible towards many instabilities. In 2018, the gate tunable magic-angle TBG is realized in laboratory Cao et al. (2018a, b), and interesting phenomena encompassing the correlated insulating phase Cao et al. (2018a), unconventional superconductivity Cao et al. (2018b); Yankowitz et al. (2019) and strange metal behavior Cao et al. (2019) are quickly discovered. Those results hint that, unlike its monolayer cousin, the gate tunable magic-angle TBG is a strongly correlated system in nature and share many common features of the phase diagram of doped cuprates, consequently spur the interests of theoretical and experimental communities on Morié physics Chen et al. (2018); Xu and Balents (2018); Yuan and Fu (2018); Po et al. (2018, 2018); Liu et al. (2018a); Dodaro et al. (2018); Huang et al. (2018); Roy and Juričić (2019); Guo et al. (2018); Xu et al. (2018); Zhang et al. (2018); Koshino et al. (2018); Wu et al. (2018); Tang et al. (2018); Cao et al. (2019); Yankowitz et al. (2019); Zou et al. (2018); Kang and Vafek (2018); Thomson et al. (2018); Venderbos and Fernandes (2018); Zhu et al. (2018a); Zhang (2018); You and Vishwanath (2018); Lian et al. (2018); Song et al. (2018); Liu et al. (2018b); Xie and MacDonald (2018); Zhu et al. (2018b); Shen et al. (2019); Liu et al. (2019); Lu et al. (2019).

Compared with cuprates, magic-angle TBG also acquires unique properties and two of them are related with the modeling of the material. First, although there are huge number of electrons in one unit cell which fill thousands of energy bands, various band calculations show that there exist an isolated band branch with four bands around charge neutrality point Bistritzer and MacDonald (2011); Suárez Morell et al. (2010); Lopes dos Santos et al. (2012); Fang and Kaxiras (2016); Trambly de Laissardière et al. (2012); Koshino et al. (2018); Kang and Vafek (2018). The four bands are made up of the spin and valley degrees of freedom of untwisted graphene. Second, the charge center forms a triangular lattice, but symmetry obstacles force one to define the effective model on a honeycomb lattice if the different band degeneracy at $\Gamma$ and $\mathbf{K}$ of the BZ were to be respected Yuan and Fu (2018); Po et al. (2018). Meanwhile, there are also obstacles from deriving a single valley tight-binding model due to chirality or mirror symmetry Zou et al. (2018). Putting all these factors together, a two orbital (counts the two valleys of the untwisted graphene) spinful lattice model on a honeycomb lattice with cluster charge interaction (considering the charge center form triangle lattice) is a good starting point to describe the system Po et al. (2018); Koshino et al. (2018); Xu et al. (2018); Kang and Vafek (2019).

But such a model is still a strongly correlated one and cannot be solved analytically. In light of the situation, we performed unbiased sign-problem-free quantum Monte Carlo (QMC) simulation to investigate such a system at charge neutrality and map out its precise phase diagram. By gradually increasing the interaction strength, the phase diagram exhibits – in a consecutive manner – a Dirac semi-metal (DSM), a plaquette valence bond solid (pVBS) and a columnar valence bond solid (cVBS) phases at weak, intermediate and strong interaction regions. The quantum phase transitions between these phases are revealed with scrutiny and we found the DSM-pVBS transition is continuous, belonging to the 3D $N=4$ Gross-Neveu chiral XY universality class; the pVBS-cVBS transition on the other hand is first order but bestowed with a sign change in the mass term of fermion bilinear, and implies that a quantum pseudo-spin Hall effect can be generated between the zigzag domains of those two insulators. The experimental relevance of our discoveries in quantum criticality and phase transitions towards to on-going investigations of TBGs is also discussed.

Model and Method —  We study an extended Hubbard model with two orbitals of spinful fermions on a honeycomb lattice. The model contains two parts $H=H_{0}+H_{U}$, where

[TABLE]

is the tight-binding part introduced in Ref. Koshino et al. (2018) and serves as a minimal model to describe of the low energy band structure of magic-angle TBG with Dirac points at charge neutrality and band splitting along $\Gamma$-M direction. Here $c^{\dagger}_{il\sigma}$ ($c_{il\sigma}$) is the creation (annihilation) operator of electron at site $i$, orbital $l=1,2$ with spin $\sigma=\uparrow,\downarrow$. Throughout this Letter, we take the nearest neighbor hopping $t$ as the energy unit. The fifth neighbor hopping ($it_{2}$ for $l=1$ and $-it_{2}$ for $l=2$) is purely imaginary and breaks orbital degeneracy along $\Gamma$-M direction. As $t_{2}/t$ is small in the material, we focus on the range of $t_{2}/t<0.6$.

For the Coulomb interaction term $H_{U}$, as the Wannier orbitals are quite extended in TBG, onsite, first, second and third neighbor repulsions are all important Koshino et al. (2018); Po et al. (2018); Kang and Vafek (2019). To capture such non-local interactions, a cluster charge Hubbard term which maintain the average filling of each elemental hexagon on the honeycomb lattice to be 4 is the genuine choice, therefore we write down

[TABLE]

where the cluster charge $Q_{\varhexagon}\equiv\sum_{i\in\varhexagon}\frac{n_{i}}{3}$ with $n_{i}=\sum_{l\sigma}c^{\dagger}_{il\sigma}c_{il\sigma}$ summing over all the six sites of the elemental hexagon. If we expand Eq. (2), the onsite, first, second and third neighbor interaction strength are $\frac{2}{3}U$, $\frac{4}{9}U$, $\frac{2}{9}U$ and $\frac{2}{9}U$, with ratio 3:2:1:1. As different range of interactions favor different kinds of order, it may require larger interaction to open a fermion gap compared to a local Hubbard model.

At any finite $U/t$, the model $H=H_{0}+H_{U}$ is nonperturbative in nature, but we found it actually immune from sign-problem due to an antiunitary symmetry Wu and Zhang (2005) at charge neutrality point, and is readily exposed to large-scale projection QMC (PQMC) simulations Blankenbecler et al. (1981); Hirsch (1985); Assaad and Evertz (2008). PQMC speaks out the ground state phase diagram, correlation functions (to determine the pattern of symmetry-breaking) and dynamical information (single-particle and collective excitation gaps above the ground state), and has been employed in several our previous studies Meng et al. (2010); Xu et al. (2017); He et al. (2018); Xu et al. (2018). The symmetry analysis and numeric implementation of PQMC are discussed in SM sup , we only mention here that the projection length is set to $\Theta=2L$ and the simulations are performed with linear system size upto $L=24$, amounts to $N_{e}=4\times L^{2}=2304$ interacting electrons on TBG model.

Phase diagram and quantum criticality —  Our phase diagram, expanded by axes $U/t$ and $t_{2}/t$, is shown in Fig. 1, the DSM, pVBS and cVBS phases are in place. The two VBS are gapped insulators. In the following parts, $t_{2}/t=0$ if it is not specified. It is interesting to notice that different from the local Hubbard Sorella and Tosatti (1992); Meng et al. (2010), $t-J$ Lang et al. (2013) and extended cluster charge models with single orbital Xu et al. (2018) on honeycomb lattice, the AB sublattice antiferromagnetic insulating phase is suppressed in our phase diagram even at very large $U/t$, which may be understood by a perturbation theory in the large $U/t$ limit Zhou et al. (2016), but the existence of two VBS phases in such a simple model is unexpected and has not been found in a local Hubbard model before.

To study the DSM-pVBS transition, we measure the bond-bond structure factor,

[TABLE]

where bond operator $B_{i,\delta}=\sum_{l,\alpha}(c_{i,l,\alpha}^{\dagger}c_{i+\delta,l,\alpha}+h.c.)$ with $\delta$ standing for one of the three nearest-neighbor bond directions ($\hat{e_{1}}$, $\hat{e}_{2}$ and $\hat{e}_{3}$) and $\hat{e}_{1}$ is chosen in the calculation. Results show that $C_{B}(\mathbf{k})$ is peaked at momenta $\mathbf{K}$ and $\mathbf{K^{\prime}}$ $(\pm\frac{4\pi}{3\sqrt{3}a_{0}},0)$ of the BZ suggesting the VBS patterns shown in the inset of Fig. 1.

To locate the DSM-pVBS transition point, we plot the correlation ratio $R_{B}(U,L)=1-\frac{C_{B}(\mathbf{K}+\delta\mathbf{q})}{C_{B}(\mathbf{K})}$ for different $U$ and system size $L$ with $|\delta\mathbf{q}|\sim\frac{1}{L}$. This quantity approaches to one (zero) in an ordered (disordered) phase, and implies a crossing for different $L$ at an critical point Kaul (2015) as showed in Fig. 2(a), where one reads $U_{c}/t=25.1(2)$. We further collapse the bond-bond structure factor with scaling relation $C_{B}(\mathbf{K},U,L)=L^{-(1+\eta)}f(L^{1/\nu}(U-U_{c})/U_{c})$ (dynamical exponent $z$ is set to 1 due to Lorentz symmetry of massless Dirac fermions), as shown in Fig. 2 (b). and obtain the critical exponents $\eta=0.80(2)$ and $\nu=1.01(3)$. Since our Dirac fermions acquire 4 degrees of freedom per site and the pVBS phase contains an emergent $U(1)$ symmetry close to the DSM-pVBS transition as shown in Refs. Xu et al. (2018); Zhou et al. (2016), we identify this transition in the 3D $N=4$ Gross-Nevue chiral XY universality class Gross and Neveu (1974); Hands et al. (1993); Rosenstein et al. (1993); Zerf et al. (2017a); Li et al. (2017); Zhou et al. (2016); Scherer and Herbut (2016); Mihaila et al. (2017); Jian and Yao (2017); Classen et al. (2017); Torres et al. (2018); Ihrig et al. (2018). The critical exponents obtained here ($\eta=0.80(2),\nu=1.01(3)$) are comparable with those calculated theoretically or numerically in literatures, as showed in Table.I of SM sup .

Gapped phases: cVBS and pVBS —  The DSM is known to possess robust massless linear dispersion at weak interaction ($U<U_{c}$) Sorella and Tosatti (1992); Assaad and Herbut (2013); Otsuka et al. (2016); Zhou et al. (2016); Sato et al. (2017); Xu et al. (2018); Chen et al. (2019), and the Dirac fermion will be gapped out in the pVBS insulator. To monitor the opening of the single-particle gap across the DSM-pVBS transition, we measure the dynamical single-particle Green’s function and follow its decay in imaginary time $G(\mathbf{k},\tau)\propto e^{-\Delta_{\text{sp}}(\mathbf{k})\tau}$ at momentum $\mathbf{K}$ for increasing system size $L$ and imaginary time displacement $\tau$, with $G(\mathbf{k},\tau)=\frac{1}{4L^{2}}\sum_{i,j,l,\sigma}e^{i\mathbf{k}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})}\langle c_{i,l,\sigma}(\frac{\tau}{2})c^{\dagger}_{j,l,\sigma}(-\frac{\tau}{2})\rangle$. The obtained $\Delta_{\text{sp}}$ for different interaction $U$ and $L$ are shown in Fig. 3 (a). It is clear that when $U<U_{c}$, $\Delta_{\text{sp}}\to 0$ and when $U>U_{c}$, $\Delta_{\text{sp}}$ goes to a finite value, which validate the picture that the DSM-pVBS transition is accompanied by the opening of the single-particle gap.

The VBS phase is related with the formation of spin singlet, either within the hexagon plaquette (pVBS) or along the nearest-neighbor bond (cVBS), and to see that one can examine the spin excitation gap, obtained from the imaginary time decay of dynamical spin-spin correlation function $C_{S}(\mathbf{q},\tau)=\frac{1}{N}\sum_{i,j}e^{i\mathbf{q}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})}\langle\mathbf{S}_{i}(\frac{\tau}{2})\mathbf{S}_{j}(-\frac{\tau}{2})\rangle$ as $C_{S}(\mathbf{K},\tau)\propto e^{-\Delta_{\text{spin}}(\mathbf{K})\tau}$. The spin operator is defined as $\mathbf{S}_{i}=\frac{1}{2}\sum_{l,\alpha,\beta}c^{\dagger}_{i,l,\alpha}(\bm{\sigma})_{\alpha,\beta}c_{i,l,\beta}$ for $t_{2}/t\neq 0$ case with $\bm{\sigma}=(\sigma^{1},\sigma^{2},\sigma^{3})$, and $(\mathbf{S}_{i})_{\mu}^{\nu}=c_{i,\mu}^{\dagger}c_{i,\nu}-\frac{\delta_{\mu\nu}}{4}\sum_{\rho=1}^{4}c_{i,\rho}^{\dagger}c_{i,\rho}$ for $t_{2}/t=0$ case where $\mu$, $\nu$ and $\rho$ denote combination of indexes of spin and orbital as components of SU(4) generator. Similar to the case of the single-particle gap, after extrapolation of $\Delta_{\text{spin}}(\mathbf{K})$ ($\Delta_{\text{spin}}$ is the smallest and degenerate at momenta $\mathbf{K}$, $\mathbf{K^{\prime}}$ and $\Gamma$) for various $L$ and $U$, one can see from Fig. 3 (b) that the spin gap is also zero when $U<U_{c}$ in the DSM phase and becomes finite when $U>U_{c}$ as the system enters the pVBS phase.

Further increase $U/t$ from the pVBS phase, a kinetic energy jump at $U/t\approx 46$ is observed, as shown in Fig. 4 (a). At the same time, the VBS correlation $C_{B}(\mathbf{K})$ also acquires a jump at the same $U$ as shown in Fig. 4 (b). These results point out that, besides $U_{c}/t=25.1(2)$ there is another first order phase transition at $U_{\text{VBS}}/t\approx 46$ between two different VBS phases. There are three non-equivalent VBS configurations, but only two of them, the pVBS and cVBS as depicted in the inset of Fig.1, can perturbatively open the single-particle gap to VBS phases. And the jumps observed in Fig. 4 (a) and (b), might be the transition between these two VBS phases.

To verify this idea, we follow Refs. Lang et al. (2013); Zhou et al. (2016) and make use of the nearest-neighboring bonds $B_{i,\delta}$ originated from sublattice A to construct a complex order parameters $D_{\mathbf{K}}=\frac{1}{L^{2}}\sum_{i}\left(B_{i,\hat{e}_{1}}+\omega B_{i,\hat{e}_{2}}+\omega^{2}B_{i,\hat{e}_{3}}\right)e^{i\mathbf{K}\cdot\mathbf{r}_{i}}$ with $\omega=e^{i\frac{2\pi}{3}}$. The Monte Carlo histogram of $D_{\mathbf{K}}$ can reveal the difference between the two VBS phases. As shown in Fig. 4 (c). The angular distribution of pVBS will peak at ${\rm arg}(D_{\mathbf{K}})=\frac{\pi}{3},\pi,\frac{5\pi}{3}$, whereas that of cVBS will peak at ${\rm arg}(D_{\mathbf{K}})=0,\frac{2\pi}{3},\frac{4\pi}{3}$. Fig. 4 (d), (e) and (f) show the corresponding histograms at three representative interaction strengths $U=44t<U_{\text{VBS}}$, $U=46t\approx U_{\text{VBS}}$, $U=48t>U_{\text{VBS}}$. It is clear that Fig. 4 (d) and (f) are inside pVBS and cVBS, respectively, Fig. 4(e) on the other hand depicts the distribution of both characters, a typical example of the coexistence at the first order transition point. In this way, the two VBS phases and their first order transition are clearly established. Similar scenario, between two Kekulé patterned superconducting states, in the context of attractive interaction on honeycomb lattice has been discussed in Ref. Roy and Herbut (2010).

Pseudo spin Hall effect in zigzag domain of pVBS and cVBS —  Our model provides the unique opportunity that pVBS and cVBS all appear in the phase diagram due to spontaneous Dirac mass generation, and one can reveal their connection with the following analysis. In the VBS phase we consider a mean field description with bond charge order, and suppress the spin and orbital degrees of freedom for the moment. Then the static VBS order becomes a modulation in the nearest neighbor hopping. As showed in the inset of Fig. 1, there are two kinds of bonds and the hopping magnitude is defined as $(1+\delta)t$ and $(1-\delta)t$. From a tight binding Hamiltonian with such bond modulation, a $4\times 4$ $k\cdot p$ Hamiltonian at $\Gamma$ point can be derived,

[TABLE]

where momentum $\tilde{\mathbf{k}}\equiv(\frac{3}{2}k_{y},\frac{\sqrt{3}}{2}ik_{x},\sqrt{3}k_{x})$, the vector $\tilde{\mathbf{s}}=(s^{1},s^{2},s^{3})$ and mass term $m=2\delta$. Here $s^{i}$ and $\tau^{i}$ are Pauli matrices in two different spaces. From the $k\cdot p$ Hamiltonian, it is clear that the bond modulation plays the role of a mass term, and there is a sign change in it across the pVBS-cVBS transition.

The sign change in the mass term motivates us to study the domain walls between pVBS and cVBS phases. Both VBS phases coexist at the first order transition point and may go through a gap close at the domain wall. Interestingly, the calculation of the spectrum of the pVBS and cVBS zigzag domain wall on a torus shows that, there exist robust helical edge states as depicted in the Fig. 5. Such edge states are protected by the combined symmetry of sublattice (chiral) and mirror (along a bond with mirror plane denoted as blue dashed line in Fig. 5(b)) and can be interpreted as the edge states of quantum pseudo spin Hall effect Kariyado and Hu (2017); Liu et al. (2017); Sorella et al. (2018).

Experimental relevance —  Recent STM experiments Kerelsky et al. (2018); Choi et al. (2019) find significant gap opening near charge neutrality point (CNP) indicating symmetry breaking order. The insulating phase at CNP is also found in many other experiments Yankowitz et al. (2019); Cao et al. (2019); Lu et al. (2019). The spontaneous symmetry breaking VBS phases found in our simulations hence provide good candidates, and a Fourier transform of the large-area STM topograph may detect addtional features at the wavelength of VBS Gutiérrez et al. (2016). Although exact estimation of cluster charge interaction $U/t$ is not easy, our results hint the possibility that magic-angle TBG is close to a QCP where the Dirac mass is spontaneously generated and opens the charge neutrality gap. Other exciting observation shows a wide range of strange metal behavior, whose existence seems robust against experimental details Cao et al. (2019). Such universal transport behavior is the hallmark of quantum critical phenomena, possibly originated from the 3D Gross-Neveu chiral XY transition between DSM and pVBS discovered here. Moreover, the strange metal behavior is pronounced near $\pm 1/4$ filling, and not well established near charge neutrality point, this is consistent with the picture that the charge neutrality point is gapped but close to the QCP, when gated, the quantum critical fluctuation kicks in and generates strange metal behavior. A final remark is that the model we studied assumed well defined valley degrees of freedom (in terms of two orbitals). Actually, if the valleys coupling is considered in TBG, it will give a single orbital model Po et al. (2018). The interesting thing is that we still find a QCP at moderate interaction strength in the single orbital model which may be related with the intriging physics near $\pm 1/4$ filling of magic-angle TBG Xu et al. (2018).

I Symmetries and band structure for different $t_{2}/t$

The tight binding part of our Hamiltonian has the following form

[TABLE]

on a honeycomb lattice, where $t$ is the hopping amplitude between the nearest neighbor sites $\langle\cdots\rangle$, $t_{2}$ is the hopping amplitude between the fifth neighbor sites $\langle\cdots\rangle^{\prime}$ and $i^{2l-1}$ sets the correct phase factor for orbital $l$, according to Ref. Koshino et al. (2018). As there are two orbitals $l=1,2$ and two sites in one unit cell, there are four bands (do not count spin) in the Brillouin zone (BZ).

For generic $t_{2}/t$, $H_{0}$ acquires the symmetry group $G=D_{3}\times U(1)\times SU(2)\times T\times\mathcal{Z}_{2}$, where $U(1)$ stems from charge conservation of each orbital, $SU(2)$ is the spin rotational symmetry, $T$ is time reversal symmetry, and $\mathcal{Z}_{2}$ is a combination of twofold rotation in real space and chirality flip in orbital space. At the limit of $t_{2}/t=0$, $H_{0}$ acquires an even higher symmetry $D_{3}\times SU(4)\times T\times\mathcal{Z}_{2}$.

The band structure of $H_{0}$ is shown in Fig. S2. When $t_{2}/t=0$, the orbital components are degenerate, so the band structure in Fig.S2(a) is the same as that of the nearest hopping honeycomb lattice, i.e. graphene, with Dirac cones at momenta $\mathbf{K}$ and $\mathbf{K^{\prime}}$ in BZ. Increasing $t_{2}/t$, upper and lower bands within $\Gamma$-M break the orbital degeneracy and the four bands manifest in this segment of high-symmetry path, with two inner ones move towards and eventually touch each other at $t_{2}/t\approx 0.31$, as shown in Fig. S2 (b). Further increasing $t_{2}/t$ these two band crosses, as shown in Fig.S2 (c), and results in a uncommon Fermi surface (FS), which are the blue streched circles close to M points in two neighboring BZs, as shown in Fig.S2(e). The parallel segments of these FS-s contribute high density of states to the tight-binding model.

Under the cluster charge interaction $U$, as shown in Fig. S2 (d), the crossing points along $\Gamma$-M will be gapped out firstly, and consequently these FS-s will disappear and result in a jump in kinetic energy of the system, with the corresponding $U$ values hightlighted as the dash line in $t_{2}/t$–$U/t$ phase diagram in Fig.1 of the main text. After the jump, the FS-s will only contain the Dirac cones at $\mathbf{K}$ and $\mathbf{K^{\prime}}$, and only when further increasing $U/t$, the Dirac cones will be gapped out as well through the $N=4$ Gross-Neveu chiral XY transition at $U_{c}$, as discussed in the main text.

II Projection QMC method and absence of sign-problem

Since we are interested in the ground state properties of the system, the projection QMC is the method of choice Assaad and Evertz (2008). In PQMC, one can obtain a ground state wave function $|\Psi_{0}\rangle$ from projecting a trial wave function $|\Psi_{T}\rangle$ along the imaginary axis $|\Psi_{0}\rangle=\lim\limits_{\Theta\to\infty}e^{-\frac{\Theta}{2}\mathbf{H}}|\Psi_{T}\rangle$, then observable can be calculated as

[TABLE]

To evaluate overlaps in above equation, we performed Trotter decomposition and $\Theta$ is discretized into $L_{\tau}$ slices ($\Theta=L_{\tau}\Delta\tau$). Each slices $\Delta\tau$ is small and the systematic error is $\mathcal{O}(\Delta\tau^{2})$. After the Trotter decomposition, we have

[TABLE]

where the non-interacting and interacting parts of the Hamiltonian is separated. To treat the interacting part, one usually employ a Hubbard Stratonovich (HS) transformation to decouple the interacting quartic fermion term to fermion bilinears coupled to auxiliary fields. For the cluster charge interaction in the TBG model, we make use of a fourth order $SU(2)$ symmetric decoupling Assaad and Evertz (2008); Xu et al. (2018)

[TABLE]

with $\alpha=\sqrt{-\Delta\tau U}$, $\gamma(\pm 1)=1+\sqrt{6}/3$, $\gamma(\pm 2)=1-\sqrt{6}/3$, $\eta(\pm 1)=\pm\sqrt{2(3-\sqrt{6})}$, $\eta(\pm 2)=\pm\sqrt{2(3+\sqrt{6})}$ and the sum is taken over the auxiliary fields $s_{\varhexagon}$ on each hexagon which can take four values $\pm 2$ and $\pm 1$. After tracing out the free fermions degrees of freedom, we obtain the following formula with a constant factor omitted

[TABLE]

where $P$ is the coefficient matrix of trial wave function $|\Psi_{T}\rangle$. In the simulation, we choose the ground state wavefunction of the half-filled non-interacting system (described by $H_{0}$) as the trial wave function. In the above formula, the $B$ matrix is defined as

[TABLE]

and has properties $B(\tau_{3},\tau_{1})=B(\tau_{3},\tau_{2})B(\tau_{2},\tau_{1})$, where we have written the coefficient matrix of interaction part as $V[\{s_{\varhexagon,\tau}\}]$ and $K$ is the hopping matrix from the $H_{0}$. The Monte Carlo sampling of auxiliary fields are further performed based on the weight defined in the sum of Eq. (S5). The measurements are performed near $\tau=\Theta/2$. Single particle observables are measured by Green’s function directly and many body correlation functions are measured from the products of single-particle Green’s function based on their corresponding form after Wick-decomposition. The equal time Green’s function are calculated as

[TABLE]

with $R(\tau)=B(\tau,0)P$, $L(\tau)=P^{\dagger}B(\Theta,\tau)$. More technique details of PQMC method, please refer to Refs Blankenbecler et al. (1981); Hirsch (1985); Assaad and Evertz (2008); Xu et al. (2018).

At half-filling, the model is sign-problem-free, this can be seen from the following analysis. Since the model is particle-hole symmetric at half-filling, one can perform a particle-hole transformation only for the orbital $l=2$, such transformation changes the cluster charge operator $Q_{\varhexagon}\equiv Q_{\varhexagon}^{1}+Q_{\varhexagon}^{2}$ to $Q_{\varhexagon}^{\prime}\equiv Q_{\varhexagon}^{1}-Q_{\varhexagon}^{2}$ in the Hamiltonian. Here $Q_{\varhexagon}^{1}$/$Q_{\varhexagon}^{2}$ is cluster charge operator for orbital-1/orbital-2. Due to the specific form of Eq. S4, the fermion bilinears after HS transformation is invariant under the antiunitary transformation

[TABLE]

where $\sigma_{0}$ is identity in spin space, $\tau_{2}$ is second Pauli matrix in orbital space, $\mathcal{K}$ is complex conjugate. As all the fermion bilinears obeys above antiunitary symmetry, and $\mathcal{U}^{2}=-1$, these properties of our model and HS decoupling form give rise to a sufficient condition for sign-problem-free Monte Carlo simulations, as proved in Ref. Wu and Zhang (2005).

III $k\cdot p$ Hamiltonian for VBS phase

In the VBS phase, we consider a mean field description with bond charge order, while forget the spin and orbital degrees of freedom. Then we can define the static VBS order as a modulation in the nearest neighbor hopping. As showed in the inset of Fig.1 of main text, there are two kinds of bonds and the hopping magnitude is defined as $(1+\delta)t$ and $(1-\delta)t$. As such kind of bond modulation enlarges unit cell to $\sqrt{3}\times\sqrt{3}$, the tight binding Hamiltonian with enlarged unit cell writes

[TABLE]

where $m\equiv 2\delta$, $\tilde{t}=(1-\delta)t$. Fig. S3 shows band structure with different mass term $m$. We can see the bond modulation folds Dirac point K and K’ of original BZ to $\Gamma$ point in VBS BZ (small BZ) and open a gap when $m\neq 0$ and the gap is exactly $m$. This can be seen more explicitly if we expand above tight binding Hamiltonian around $\Gamma$ point and perform a further down folding from six bands to four bands with unitary transformation

[TABLE]

Then we get a $k\cdot p$ Hamiltonian $H_{\text{eff}}(\mathbf{k})=-t\left(\tilde{\mathbf{k}}\cdot\tilde{\mathbf{s}}\tau^{2}+ms^{0}\tau^{1}\right)$ where only linear terms of $k$ and $m$ are reserved. We change the basis a little bit ($\tau^{1}\rightarrow\tau^{3}$ and $\tau^{2}\rightarrow\tau^{2}$), and get the final form as shown in Eq.(4) of the main text. For convenient, we also show it here

[TABLE]

where momentum $\tilde{\mathbf{k}}\equiv(\frac{3}{2}k_{y},\frac{\sqrt{3}}{2}ik_{x},\sqrt{3}k_{x})$, the vector $\tilde{\mathbf{s}}=(s^{1},s^{2},s^{3})$. Here $s^{i}$ and $\tau^{i}$ are Pauli matrices in two different spaces. From the $k\cdot p$ Hamiltonian, it is clear that the bond modulation plays the role of a mass term, and there is a sign change in it across the pVBS-cVBS transition.

IV Critical exponents for 3D $N=4$ Gross-Nevue chiral XY universality class

We compare the critical exponents for 3D $N=4$ Gross-Nevue chiral XY universality class calculated by different methods. As showed in Table S1, the critical exponents we get are comparable with existed results. We want to remark, there are still differences between the critical exponents got from different methods. Much more efforts in the direction of developments of advanced analytical techniques for higher order calculations, as well as the developments of technically sound numeric methods to achieve much larger system sizes are requiring.