Valence Bond Orders at Charge Neutrality in a Possible Two-Orbital Extended Hubbard Model for Twisted Bilayer Graphene

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 semimetal to valence bond solid (VBS) insulator is established. The quantum critical point is discovered to belong to the 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.

Figure 1

Ground state phase diagram. As a function of $U/t$, the DSM, pVBS, and cVBS phases reveal themselves. The fifth neighbor hopping strength $t_2$ is marked as an orange line in the inset. Inside DSM, the $t_2/t$ ratio modifies the band degeneracy and Fermi surface topology, and the black dashed line signifies the corresponding band structure crossover, where excitons formed between orbitals might condense. The different patterns of the valence bonds are shown in the insets. The DSM-pVBS transition is continuous and shown to belong to the 3D $N=4$ chiral $XY$ Gross-Neveu universality. The pVBS-cVBS transition is first order, but could carry topological edge states in terms of the quantum pseudospin Hall effect in the pVBS-cVBS zigzag domain wall.

Figure 2

(a) The bond-bond correlation ratio $R_B$ and (b) data collapse analysis of the structure factor $C_B(\mathbf{K})$ at $t_2/t=0$ as a function of $U/t$ with $L=12,15,\dots ,24$. The crossing of $R_B$ in (a) gives the DSM-pVBS critical point $U_c/t=25.1(2)$. The data collapse in (b) gives the 3D $N=4$ Gross-Neveu chiral $XY$ exponents $\eta =0.80(2)$, $\nu =1.01(3)$.

Figure 3

(a) The $1/L$ extrapolation of the single-particle gap ${\mathrm{\Delta }}_{\mathrm{sp}}(\mathbf{K})$, the gap opens between $U/t=22$ and $U/t=28$, consistent with the $U_c/t$ obtained from the bond correlation ratio in Fig. 2. (b) The $1/L$ extrapolation of spin gap ${\mathrm{\Delta }}_{\mathrm{spin}}(\mathbf{K})$, the spin gap opens hand in hand with the single-particle gap as the establishment of pVBS order.

Figure 4

(a) Kinetic energy per site of the system for $U$ at large values. The sharp jump signifies a first order transition. (b) $C_B(\mathbf{K})$ for the same process, a jump in VBS order is also observed, suggesting this is a transition between different VBS phases. (c) Angular dependence of the complex order parameter $D_{\mathbf{K}}$. Black dots represent ideal pVBS order, and red dots represent ideal cVBS order. (d)–(e) Histogram of $D_{\mathbf{K}}$ at different interaction strengths $U<U_{\mathrm{VBS}}$, $U\approx U_{\mathrm{VBS}}$, and $U>U_{\mathrm{VBS}}$.

Figure 5

(a) Spectrum of the zigzag domain wall of pVBS and cVBS. The red line denotes the helical edge states. In the calculation, both the width of pVBS and cVBS strips are set to 24 times the honeycomb unit cell. (b) The zigzag domain wall of pVBS and cVBS with the bottom part the pVBS phase and top part the cVBS phase. The hopping strength is $-1.1t$ for strong bonds and $-0.9t$ for weak bonds. The pVBS and cVBS are connected by vertical bonds with hopping strength $-t$. The domain wall has periodic boundary condition along the $x$ direction, and is put on a torus (the upper boundary and lower boundary are also connected by vertical bonds with hopping strength $-t$). The shaded region denotes the super unit cell of the domain wall.