⇦Chinese Physics Letters, 2021, Vol. 38, No. 7, Article code 077305Express Letter Momentum Space Quantum Monte Carlo on Twisted Bilayer Graphene Xu Zhang (张栩)1†, Gaopei Pan (潘高培)2,3†, Yi Zhang (张燚)4, Jian Kang (康健)5, and Zi Yang Meng (孟子杨)1,2* Affiliations 1Department of Physics and HKU-UCAS Joint Institute of Theoretical and Computational Physics, The University of Hong Kong, Pokfulam Road, Hong Kong SAR, China 2Beijing National Laboratory for Condensed Matter Physics and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China 3School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100190, China 4Kavli Institute for Theoretical Sciences, University of Chinese Academy of Sciences, Beijing 100190, China 5School of Physical Science and Technology and Institute for Advanced Study, Soochow University, Suzhou 215006, China Received 28 May 2021; accepted 2 June 2021; published online 4 June 2021 Supported by the RGC of Hong Kong SAR of China (Grant Nos. 17303019, 17301420, and AoE/P-701/20), the National Key Research and Development Program of China (Grant No. 2016YFA0300502), the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB33000000 and XDB28000000), the National Natural Science Foundation of China (Grant Nos. 11674278, 12004383, 12074276, and 12074276), the Fundamental Research Funds for the Central Universities, and the Priority Academic Program Development (PAPD) of Jiangsu Higher Education Institutions.
†These two authors contributed equally to this work.
*Corresponding author. Email: zymeng@hku.hk Citation Text: Zhang X, Pan G P, Zhang Y, Kang J, and Meng Z Y 2021 Chin. Phys. Lett. 38 077305    Abstract We report an implementation of the momentum space quantum Monte Carlo (QMC) method on the interaction model for the twisted bilayer graphene (TBG). The long-range Coulomb repulsion is treated exactly with the flat bands, spin and valley degrees of freedom of electrons taking into account. We prove the absence of the minus sign problem for QMC simulation when either the two valleys or the two spin degrees of freedom are considered. By taking the realistic parameters of the twist angle and interlayer tunnelings into the simulation, we benchmark the QMC data with the exact band gap obtained at the chiral limit, to reveal the insulating ground states at the charge neutrality point (CNP). Then, with the exact Green's functions from QMC, we perform stochastic analytic continuation to obtain the first set of single-particle spectral function for the TBG model at CNP. Our momentum space QMC scheme therefore offers the controlled computation pathway for systematic investigation of the electronic states in realistic TBG model at various electron fillings. DOI:10.1088/0256-307X/38/7/077305 © 2021 Chinese Physics Society Article Text Twisted bilayer graphene (TBG) and other Moiré systems have attracted great theoretical[1–6] and experimental[7–20] interests in condensed matter and 2D quantum material communities. As experiments have discovered the correlated insulating phases at various integer fillings and the proximate superconductivity (SC) phases, there appears a key question how to model and understand the properties of the insulating phases, both for their own sake and that could eventually provide a clue for understanding the mechanism of the superconductivity in TBG systems. Many experimental and theoretical works have indicated the interplay between the nontrivial topology and strong interaction as the essential ingredients for the understanding of the electronic correlations in such materials, therefore pointing out a proper model for TBG system shall be significantly different from the typical Hubbard-type Hamiltonian with on-site interactions.[21–29] However, the nature of the insulating states discovered in the material is still under debate. On the one hand, the analytical and the Hartree–Fock calculations at various integer fillings have found the quantum anomalous hall (QAH) and the intervalley coherent (IVC) states as the ground states without breaking the translation symmetry, suggesting that the physics is similar to the quantum Hall ferromagnetism at the lowest Landau level (LLL).[30–40] Such similarity also led to the proposal of the skyrmion SC for the mechanism of SC discovered near the insulating phases at $\nu = \pm 2$.[41,42] On the other hand, recent numerical calculations based on the density matrix renormalization group (DMRG)[41,43–45] and exact diagonalization (ED)[46–49] have discovered a larger manifold of nearly degenerate states strongly competing with each other even in the strong coupling regime, hinting that such systems could contain much more complicated physics than that of the LLL. As a consequence, there is a crying need for applying more extensive numerical methodology that can unbiasedly solve larger system sizes to fully settle the mechanism of the insulating phases at various integer fillings. However, this is by no means an easy task. In TBG each Moiré superlattice unit cell contains more than $10^4$ carbon atoms (in the case close to the first magic angle), resulting in the same number of bands in the Moiré Brillouin zone (mBZ). It is impossible to include such a huge number of bands in any realistic calculations for strongly correlated electrons. Fortunately, all experiments have found that the correlated physics emerges only when the chemical potential lies inside the flat bands[7–9] and the large band gap that separates the flat bands and remote bands allows one to focus on the flat bands only to study the electronic correlations.[13] While the parameters of the Hamiltonian has been changed by integrating out the states on the remote bands in the presence of the Coulomb interactions, the effective Hamiltonian is still given by the Bistritzer–MacDonald (BM) model in momentum space[50] with the projected Coulomb interactions onto the flat bands, and it has been significantly simplified (yet still difficult) for realistic analytical and numerical calculations.[31,36,43,46,51–56] In light of the situation, the large-scale quantum Monte Carlo (QMC) method presents itself as the ideal choice of method to solve these TBG models at integer fillings. QMC solves the correlated electron lattice models in path-integral such that both static and dynamic properties, at finite temperature and ground state, can be obtained in unbiased manner with only statistical errors. The extrapolation to the thermodynamic limit is also possible, when the computation complexity increases polynomially with the system sizes, i.e., absence of minus sign problem. Many important features of correlated electron system such as the antiferromagnetic Mott insulator in square lattice Hubbard model,[57] non-Fermi liquid at quantum critical points,[58,59] to name a few, have been discovered from QMC simulations. In the case of TBG, by now there have been few QMC simulations in a real-space lattice model with extended interactions where interaction-driven topological state, IVC and translational symmetry-breaking insulators are found.[60–64] However, generic and systematic QMC analysis for BM-type models with flat bands, spin and valley degrees of freedom and in particular, the long-range Coulomb interaction to be fully respected in momentum space, is still missing. This is the knowledge gap we want to fill in. In this work, we develop a momentum space QMC method[59,65–67] for the aforementioned TBG models. We first prove the absence of the minus sign problem for QMC simulation at integer fillings when either the two valleys or the two spin degrees of freedom are considered. Then, by taking the realistic parameters of the twist angle and interlayer tunnelings into account, we benchmark the QMC data with the exact band gap obtained with a $6\times6$ momentum mesh in the mBZ, to reveal the insulating ground states at the charge neutrality point (CNP). Finally, by combining the QMC simulation with the stochastic analytic continuation, we obtain the first set of single-particle spectra at chiral limit and realistic parameter at CNP. Our momentum space QMC scheme therefore offers the controlled computation pathway for systematic investigation of the electronic states in the realistic TBG model at various electron fillings. Continuum Model. We start from the BM model[1–6] in plane wave basis: $$H^{\tau}_{\rm BM} \left({\boldsymbol k}\right)=\sum_{{\boldsymbol k'}}H^{\tau}_{\rm BM} {}_{{\boldsymbol k},{\boldsymbol k'}}\; e^{-i{\boldsymbol k} \cdot {{\boldsymbol r}}}e^{i{\boldsymbol k \prime} \cdot {{\boldsymbol r}}}, $$ where $$\begin{align} &H^{\tau}_{\rm BM} {}_{{\boldsymbol k},{\boldsymbol k'}}\\ =\,&\delta_{{\boldsymbol k},{\boldsymbol k'}}\begin{pmatrix} -\hbar v_F ({{{\boldsymbol k}}}-{\boldsymbol K}_1^{\tau}) \cdot \boldsymbol{\sigma}^{\tau} & U_0 \\ U_0^† & -\hbar v_F ({{{\boldsymbol k}}}-{\boldsymbol K}_2^{\tau}) \cdot \boldsymbol{\sigma}^{\tau} \end{pmatrix} \\ & +\begin{pmatrix} 0 & U_1^{\tau} \delta_{_{\scriptstyle {\boldsymbol k},{\boldsymbol k'}+\tau {\boldsymbol G}_1 }} \\ U_1^{\tau †} \delta_{_{\scriptstyle {\boldsymbol k},{\boldsymbol k'}-\tau {\boldsymbol G}_1 }} & 0 \end{pmatrix} \\ &+\begin{pmatrix} 0 & U_2^{\tau} \delta_{_{\scriptstyle {\boldsymbol k},{\boldsymbol k'}+\tau({\boldsymbol G}_1+{\boldsymbol G}_2) }} \\ U_2^{\tau †} \delta_{_{\scriptstyle {\boldsymbol k},{\boldsymbol k'}-\tau({\boldsymbol G}_1+{\boldsymbol G}_2) }} & 0 \end{pmatrix} , \end{align} $$ with $\tau=\pm 1$ denoting the valley index, $\boldsymbol{\sigma}^{\tau}=(\tau \sigma_x,\sigma_y)$ defines the A,B sublattices of the monolayer graphene. ${\boldsymbol K}^{\tau}_1$ and ${\boldsymbol K}^{\tau}_2$ are the corresponding Dirac points of the bottom and top layers of graphene that are now twisted by angles $\mp\frac{\theta}{2}$, and ${{\boldsymbol k}}\in$ mBZ and ${\boldsymbol G}_1$ and ${\boldsymbol G}_2$ are the reciprocal vectors of mBZ, as shown in Fig. 1(a). The interlayer tunneling between the Dirac states is described by the matrix $$\begin{align} U_0=\,& \begin{pmatrix} u_0 & u_1 \\ u_1 & u_0 \end{pmatrix} ,~~U_1^{\tau}= \begin{pmatrix} u_0 & u_1 e^{-\tau \frac{2\pi}{3}i} \\ u_1 e^{\tau \frac{2\pi}{3}i} & u_0 \end{pmatrix},\\ U_2^{\tau}=\,& \begin{pmatrix} u_0 & u_1 e^{\tau \frac{2\pi}{3}i} \\ u_1 e^{-\tau \frac{2\pi}{3}i} & u_0 \end{pmatrix} , \end{align} $$ where $u_0$ and $u_1$ are the intra-sublattice and inter-sublattice interlayer tunneling amplitudes. The flatness of the lowest two bands per spin per valley in the chiral limit ($u_0=0$) is determined by the dimensionless parameter $\alpha=\frac{u_1}{\hbar v_F k_{\theta}}$ with $k_{\theta}=8\pi\sin(\theta/2)/(3a_0)$ and the lattice constant of the monolayer graphene $a_0=0.246$ nm. In this Letter, we choose $\hbar v_F/a_0=2.37745$ eV, the twist angle $\theta=1.08^{\circ}$ and $u_1=0.11$ eV which leads to $\alpha=0.586$, the value corresponding to the first magic angle where the lowest two bands become completely flat in the chiral limit. We perform the QMC simulation at both the chiral limit $(u_0=0)$ and the more realistic case $(u_0=0.06$ eV), which leads to a bandwidth of 1.08 meV. The eigenstate of $H_{\rm BM}^{\tau}$ can be written in the Bloch wavefunction form $$ \psi_{m,\tau,{{\boldsymbol k}}}^X({{\boldsymbol r}})=\sum_{{\boldsymbol G}}u_{_{\scriptstyle m,\tau;{\boldsymbol G},X }}({{\boldsymbol k}}) e^{i({{\boldsymbol k}}+{\boldsymbol G})\cdot {{\boldsymbol r}}}. $$ Here $X=\{A_1, B_1, A_2, B_2\}$ denotes the layer and sublattice indices, and $u_{_{\scriptstyle m,\tau;{\boldsymbol G}, X}}({{\boldsymbol k}})$ is the Bloch wave-function with the eigen-energy $\epsilon_{m{{\boldsymbol k}}\tau}$. Here, $m$ and $\tau$ are the band and valley indices and we omit the spin index $s$ now since the Hamiltonian is spin independent. The range of $m$ can be large [consider the couplings $\{ {{\boldsymbol k}},{{\boldsymbol k}}+{\boldsymbol G}_1,{{\boldsymbol k}}+{\boldsymbol G}_1+{\boldsymbol G}_2, \cdots \}$ in Fig. 1(a) which has $m\in 1,2,\ldots,M$ elements, then $H^{\tau}_{\rm BM} {}_{{\boldsymbol k},{\boldsymbol k'}}$ is a $4M \times 4M$ matrix], we select the two flat bands and denote them as $m=1,2$ and consider the projected Coulomb interactions onto these bands in this work.
cpl-38-7-077305-fig1.png
Fig. 1. (a) Upper panel shows the Brillouin zone foldings in TBG, with the small hexagon representing the mBZ. ${\boldsymbol K}^{\pm}_{1,2}$ are the Dirac points for valley $\pm$ and layers 1 and 2, $\theta$ is the rotating angle. In lower panel, ${\boldsymbol G}_1$ and ${\boldsymbol G}_2$ are the reciprocal lattice vectors of the mBZ. The yellow dots are the allowed momentum transfer ${\boldsymbol q}$ points in QMC. Here we consider momentum transfer up to ${\boldsymbol G}$, which means that for a $6\times 6$ mesh in mBZ, the allowed number of ${\boldsymbol q}$ is 126. [(b), (c)] The comparison at the CNP of the single-particle excitation gap at the chiral limit ($u_0=0$) and realistic case ($u_0=0.06$ eV) between exact solution and the QMC results. The QMC gaps are obtained from fitting the imaginary time decay of the fermion Green's function.
Interaction Model and QMC Implementations. We develop a momentum space QMC scheme to solve the interaction Hamiltonian in the band basis: $$ H_{\rm int} = \frac{1}{2 \varOmega} \sum_{{\boldsymbol q},{\boldsymbol G},|{\boldsymbol q}+{\boldsymbol G}|\neq0} V({\boldsymbol q}+{\boldsymbol G}) \delta \rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}} \delta \rho_{_{\scriptstyle -{\boldsymbol q}-{\boldsymbol G}}},~~ \tag {1} $$ as shown in Fig. 1(a), ${\boldsymbol q}\in$ mBZ and ${\boldsymbol q}+{\boldsymbol G}$ represents a vector in extended mBZ, and we consider the momentum transfer upto the distance of ${\boldsymbol G}_1$ and ${\boldsymbol G}_2$,[52,53] as denoted by the yellow dots in the figure. The single gate Coulomb interaction is long-ranged: $$\begin{align} V({\boldsymbol q}) =\,&\frac{e^{2}}{4 \pi \varepsilon} \int d^{2} {{\boldsymbol r}}\left(\frac{1}{{{\boldsymbol r}}}-\frac{1}{\sqrt{{{\boldsymbol r}}^{2}+d^{2}}}\right) e^{i {\boldsymbol q} \cdot {{\boldsymbol r}}}\\ =\,&\frac{e^{2}}{2 \varepsilon} \frac{1}{q}\left(1-e^{\rm -q d}\right), \end{align} $$ with $\frac{d}{2}$ being the distance between graphene layer and single gate and $\epsilon$ the dielectric constant. The definition of $\delta \rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}}$ is $$\begin{align} \delta \rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}} = \,&\sum_{{{\boldsymbol k}} \in mBZ, m_{1}, m_{2},\tau, s} \lambda_{m_1,m_2,\tau}({{\boldsymbol k}},{{\boldsymbol k}}+{\boldsymbol q}+{\boldsymbol G})\\ \,&\cdot\Big(d_{{{\boldsymbol k}}, m_{1}, \tau, s}^† d_{{{\boldsymbol k}}+{\boldsymbol q}, m_{2}, \tau, s} - \frac{1}{2} \delta_{{\boldsymbol q},0} \delta_{m_1,m_2}\Big)\\ ={}&(\delta\rho_{_{\scriptstyle -{\boldsymbol q}-{\boldsymbol G}}})^†,~~ \tag {2} \end{align} $$ with the form factor $\lambda$ defined as $ \lambda_{m_1 m_2, \tau}({{\boldsymbol k}},{{\boldsymbol k}}+{\boldsymbol q}+{\boldsymbol G})=\sum_{{\boldsymbol G}',X}u^{*}_{_{\scriptstyle m_1,\tau;{\boldsymbol G}',X}}({{\boldsymbol k}})u_{_{\scriptstyle m_2,\tau;{\boldsymbol G}'+{\boldsymbol G},X}}({{\boldsymbol k}}+{\boldsymbol q})$. Physically, $\delta \rho$ is the electron density operator relative to the decoupled bilayer graphene at the CNP. The interaction in Eq. (1) differs from the normal ordered version as it can properly account for the renormalization effect of the remote bands to the flat bands.[53] The symmetry properties of $\lambda_{m_1,m_2,\tau}({{\boldsymbol k}},{{\boldsymbol k}}+{\boldsymbol q}+{\boldsymbol G})$ are discussed in the Supplemental Material (SM).[68] One can see that $\delta \rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}}$ in single particle basis is a block diagonal matrix according to the $\tau,s$ indices. Since $V({\boldsymbol q}+{\boldsymbol G})=V(-{\boldsymbol q}-{\boldsymbol G})$, terms in $H_{\rm int}$ can be written in a form ready to be QMC decoupled as $$\begin{align} &\sum_{{\boldsymbol q},{\boldsymbol G},|{\boldsymbol q}+{\boldsymbol G}|\neq0} V({\boldsymbol q}+{\boldsymbol G}) \delta \rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}} \delta \rho_{_{\scriptstyle -{\boldsymbol q}-{\boldsymbol G}}}\\ = \,&\sum_{|{\boldsymbol q}+{\boldsymbol G}|\neq0} \frac{V({\boldsymbol q}+{\boldsymbol G})}{2} \Big[\big(\delta\rho_{_{\scriptstyle -{\boldsymbol q}-{\boldsymbol G}}}+\delta\rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}}\big)^{2}\\ &-\big(\delta\rho_{_{\scriptstyle -{\boldsymbol q}-{\boldsymbol G}}}-\delta\rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}}\big)^{2}\Big],~~ \tag {3} \end{align} $$ where $\sum_{|{\boldsymbol q}+{\boldsymbol G}|\neq0}$ means summation to half of the allowed values of ${\boldsymbol q}$ and ${\boldsymbol G}$. According to the discrete Hubbard–Stratonovich transformation,[60,61,69] $e^{\alpha \hat{O}^{2}}=\frac{1}{4} \sum_{l=\pm 1,\pm 2} \gamma(l) e^{\sqrt{\alpha} \eta(l) \hat{o}}+O\left(\alpha^{4}\right)$, where $\gamma(\pm 1)=1+\frac{\sqrt{6}}{3}$, $\gamma(\pm 2)=1-\frac{\sqrt{6}}{3}$, $\eta(\pm 1)=\pm \sqrt{2(3-\sqrt{6})}$ and $\eta(\pm 2)=\pm \sqrt{2(3+\sqrt{6})}$, we can rewrite the partition function of Hamiltonian (1) in the imaginary time discretization as $$\begin{align} Z={}&{\rm Tr}\Big\{\prod_{t}e^{-\Delta \tau H_{\rm i n t}(t)}\Big\}\\ ={}&{\rm Tr}\Big\{\prod_{t} \exp\Big\{-\Delta \tau \frac{1}{4 \varOmega}\sum_{|{\boldsymbol q}+{\boldsymbol G}|\neq0} V({\boldsymbol q}+{\boldsymbol G})\\ &\cdot\left[\left(\delta\rho_{_{\scriptstyle -{\boldsymbol q}-{\boldsymbol G}}}+\delta\rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}}\right)^{2}-\left(\delta\rho_{_{\scriptstyle -{\boldsymbol q}-{\boldsymbol G}}}-\delta\rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}}\right)^{2}\right] \Big\}\Big\} \\ \approx \,& \sum_{\{l_{|{\boldsymbol q}|,t}\}} \prod_{t} \Big[ \prod_{|{\boldsymbol q}+{\boldsymbol G}|\neq0}\frac{1}{16} \gamma\left(l_{|{\boldsymbol q}|_1,t}\right) \gamma\left(l_{|{\boldsymbol q}|_2,t}\right)\Big]\\ &\cdot {\rm Tr}\Big\{\prod_{t}\Big\{\prod_{|{\boldsymbol q}+{\boldsymbol G}|\neq0}\exp \big[i \eta\left(l_{|{\boldsymbol q}|_1,t}\right) A_{{\boldsymbol q}}\left(\delta\rho_{-{\boldsymbol q}}+\delta\rho_{{\boldsymbol q}}\right)\big]\\ &\cdot \exp \big[\eta\left(l_{|{\boldsymbol q}|_2,t}\right) A_{{\boldsymbol q}}\left(\delta\rho_{-{\boldsymbol q}}-\delta\rho_{{\boldsymbol q}}\right)\big]\Big\}\Big\},~~ \tag {4} \end{align} $$ where $t$ is the imaginary time index with step $\Delta\tau$, $A_{{\boldsymbol q}+{\boldsymbol G}} =\sqrt{\frac{\Delta \tau}{4} \frac{V({\boldsymbol q}+{\boldsymbol G})}{\varOmega}}$ and $\{l_{|{\boldsymbol q}|_1,t},l_{|{\boldsymbol q}|_2,t},l_{0,t}\}$ are the four-component auxiliary field that lives in the space-time configuration of the path-integral. For each realization of the auxiliary field configuration, the fermion determinant can be evaluated exactly as the configurational weight, the QMC simulation is performed along a Markov chain of such configurations and the important sampling can be carried out with the physical observables (such as single-particle Green's function) computed through ensemble average.[70] One shall be careful about the approximation “$\approx$” in Eq. (4), since $$\begin{align} &\delta\rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}} \delta\rho_{_{\scriptstyle {\boldsymbol q}^{\prime}+{\boldsymbol G}^{\prime}}}-\delta\rho_{_{\scriptstyle {\boldsymbol q}^{\prime}+{\boldsymbol G}^{\prime}}} \delta\rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}}\\ =\,&\sum_{{{\boldsymbol k}}, m_{1}, m_{2},\tau,s}[\lambda_\tau({{\boldsymbol k}}, {{\boldsymbol k}}+{\boldsymbol q}+{\boldsymbol G}) \lambda_\tau\left({{\boldsymbol k}}+{\boldsymbol q},{{\boldsymbol k}}+{\boldsymbol q}+{\boldsymbol q}^{\prime}+{\boldsymbol G}^{\prime}\right)\\ &-\lambda_\tau\left({{\boldsymbol k}}, {{\boldsymbol k}}+{\boldsymbol q}^{\prime}+{\boldsymbol G}^{\prime}\right) \lambda_\tau\left({{\boldsymbol k}}+{\boldsymbol q}^{\prime}, {{\boldsymbol k}}+{\boldsymbol q}^{\prime}+{\boldsymbol q}+{\boldsymbol G}\right)]_{m_{1}, m_{2}}\\ &\cdot d_{{{\boldsymbol k}},m_1,\tau,s}^† d_{{{\boldsymbol k}}+{\boldsymbol q}+{\boldsymbol q}^{\prime},m_2,\tau,s}, \end{align} $$ which means $\left[\delta\rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}}, \delta\rho_{_{\scriptstyle {\boldsymbol q}^{\prime}+{\boldsymbol G}^{\prime}}}\right] \neq 0$. However, when the number of ${\boldsymbol q}$ is limited as shown in Fig. 1(a), i.e., allowing momentum transfer upto ${\boldsymbol G}$, our results show that the systematic discretization errors are acceptable. By setting $\varepsilon=7 \varepsilon_{0}$, the Moiré lattice vector $L_{M}=\frac{a_{0}}{2 \sin \left(\frac{\theta}{2}\right)}$, the area of system $\varOmega=N_{{{\boldsymbol k}}}\frac{\sqrt{3}}{2} L_{M}^{2}$ with $N_{{{\boldsymbol k}}}$ the number of ${{\boldsymbol k}}$ points in mBZ (here we have $N_{{{\boldsymbol k}}}=36$ for the $6\times6$ mesh), and the gate distance $d=40$ nm, mBZ reciprocal lattice vector $|{\boldsymbol G}|=\frac{4 \pi}{\sqrt{3} L_{M}}$, we have $ \frac{V(\bar{{\boldsymbol q}})}{\varOmega} \approx 0.01585 \frac{1}{\sqrt{N_{k}}\bar{{\boldsymbol q}}}\left(1-e^{-22.36 \frac{1}{\sqrt{N_{k}}} \bar{{\boldsymbol q}}}\right)$ eV, where $\bar{{\boldsymbol q}}$ is the distance between momenta in mBZ by setting two nearest ${{\boldsymbol k}}$ points with unit length. Discussion of the Sign-Problem. We note in Eq. (4), the exponential parts of decoupled interaction are anti-Hermitian, the following three statements about the sign structure of the QMC fermion determinants are in order: Statement 1. Considering single valley and single spin without kinetic terms at half filling, the sign of the determinant is always real. Proof. In our decoupled Hamiltonian, the configurational probability is proportional to $\exp[-\frac{1}{2}\sum_{j}\operatorname{Tr}(M_{j})]\operatorname{{\det}}\left(I+e^{M_{1}}e^{M_{2}}\dots e^{M_{n}}\right)$. Here $M_{j}$ are anti-Hermitian matrices in single particle basis and $e^{-\frac{1}{2}\sum_{j}\operatorname{Tr}(M_{j})}$ comes from constant terms in $\delta\rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G}}}$. Since $e^{M_j}$ are unitary matrices, $U=e^{M_{1}}e^{M_{2}}\dots e^{M_{n}}$ is also unitary with eigenvalue $e^{i\lambda_j}$. Set $\operatorname{{\det}}\left(U\right)=\exp[\sum_{j}\operatorname{Tr}(M_{j})] =\exp[\sum_{j}i\lambda_j]=e^{i\varGamma}$, $e^{-i\frac{\varGamma}{2}}\operatorname{{\det}}\left(I+U\right)=e^{-i\frac{\varGamma}{2}}\prod_{j}\left(1+e^{i \lambda_{j}}\right) .$ For any term $\exp[i (\Sigma_{k\in A} \lambda_{k} -\frac{\varGamma}{2})]$, $A\subseteq \left\lbrace 1,2,\dots ,n\right\rbrace $, there is always a term $\exp[i (\Sigma_{k\notin A} \lambda_{k}-\frac{\varGamma}{2})]=\exp[-i (\Sigma_{k\in A} \lambda_{k}-\frac{\varGamma}{2})]$, so adding all terms together will always be real. Statement 2. Considering single valley and double spin without kinetic terms at half filling, there is no sign problem. Proof. It is straightforward to see this result according to Statement 1 by noticing the other spin just gives a copy so that the real sign will become a non-negative sign. However, kinetic terms could change this result. See SM[68] for details. Statement 3. Considering single spin and double valley with flat band kinetic terms at half filling, there is no sign problem. Proof. One can relabel $d_{k,m,-\tau,s}$ as $-m* \tilde{d}_{k,-m,-\tau,s}^†$ in $-\tau$ subspace, then prove that the single-particle matrixes between two valleys satisfy $\delta \rho_{_{\scriptstyle {\boldsymbol q}+{\boldsymbol G},-\tau}} = -\delta \rho_{_{\scriptstyle -{\boldsymbol q}-{\boldsymbol G},\tau}}$. Thus $M_{j,-\tau}=M_{j,\tau}^*$. This transformation will keep flat band kinetic matrices intact between two valleys. Thus, the determinant of valley $-\tau$ is complex conjugated with that of valley $\tau$. We note that the similar observation has also been pointed out in Ref. [71], i.e., the TBG Hamiltonian at CNP after QMC decoupling is invariant under anti-unitary particle-hole symmetry and thus free of the sign problem. We organize these statements in Table 1. And we noted that Ref. [72] also shows some similar results.
Table 1. List of the sign structure for TBG Hamiltonian.
Degrees of freedom Kinetic terms Sign structure
Single valley single spin No Real
Single valley double spin No Non-negative
Double valley single spin Flat bands Non-negative
Double valley double spin Flat bands Non-negative
cpl-38-7-077305-fig2.png
Fig. 2. Single-particle spectra obtained from QMC+SAC at CNP with (a) chiral limit ($u_0=0$) and [(b), (c)] $u_0=0.06$ eV. The exact gaps are the same as in Figs. 1(b) and 1(c).
QMC Results and Discussions. With such understanding, we carried out QMC simulations at CNP for TBG Hamiltonian in momentum space, on the grids of $6\times 6$ in mBZ, both at chiral limit $(u_0=0)$ and with realistic parameter ($u_0=0.06$ eV). To make sure our QMC have converged to the ground state, we set the temperature $T=0.667$ meV in the simulations, which turns out to be magnitudes smaller than the obtained gap, and divide the inverse temperature $\beta=1/T$ to 150 pieces with $\Delta\tau=0.01$ such that the Trotter error negligible. Figures 1(b) and 1(c) show the comparison of the single-particle excitation gap at the chiral limit ($u_0=0$) and $u_0=0.06$ eV, obtained from fitting the imaginary decay of Green's function in QMC simulations, respectively. At the chiral limit, the two bands are degenerate, whereas in the realistic case, we diagonalize Green's function at every ${{\boldsymbol k}}$ point in the $2\times2$ band basis. The agreement between the exact gaps[55] and the QMC ones is perfect. Figure 2 shows the single-particle spectra, obtained from applying stochastically analytic continuation (SAC)[73–77] upon the imaginary time Green's function from QMC simulations. Such a QMC+SAC scheme has been shown to reliably reveal many interesting dynamical features in various strongly correlated systems.[77–86] Figure 2(a) shows the single-particle spectrum at the chiral limit ($u_0=0$) where the two bands are degenerate and Figs. 2(b) and 2(c) the spectral of the two bands at $u_0=0.06$ eV where they are not degenerate. In all cases, the spectral are particle-hole symmetric. With the establishment of such unbiased QMC computational framework, as summarized in Table 1, the controlled computation pathway for the interaction effects in the realistic TBG model is clearly opening up. The questions of the nature of the ground states at different integer fillings and the TBG phase diagram at different twist angle, chiral ratio, hNB alignment, the skyrmion SC and the dynamic and spectral properties, etc, can now be investigated as those having been investigated with QMC simulations in other exotic strongly correlated electron systems.[87–90] Acknowledgments. ZYM thanks Xi Dai for the insightful discussion and continuous encouragement for addressing the momentum-space solutions of TBG. We thank the Computational Initiative at the Faculty of Science and the Information Technology Services at the University of Hong Kong for their technical support and generous allocation of CPU time. References Localization of Dirac Electrons in Rotated Graphene BilayersNumerical studies of confined states in rotated bilayers of grapheneMoire bands in twisted double-layer grapheneElectronic properties of graphene-based bilayer systemsGraphene Bilayer with a Twist: Electronic StructureContinuum model of the twisted graphene bilayerUnconventional superconductivity in magic-angle graphene superlatticesCorrelated insulator behaviour at half-filling in magic-angle graphene superlatticesTunable correlated Chern insulator and ferromagnetism in a moiré superlatticeMaximized electron interactions at the magic angle in twisted bilayer grapheneElectronic Compressibility of Magic-Angle Graphene SuperlatticesSuperconductors, orbital magnets and correlated states in magic-angle bilayer grapheneSpectroscopic signatures of many-body correlations in magic-angle twisted bilayer grapheneCorrelated states in twisted double bilayer grapheneStrongly correlated Chern insulators in magic-angle twisted bilayer grapheneUnconventional sequence of correlated Chern insulators in magic-angle twisted bilayer grapheneObservation of superconductivity in bilayer graphene/hexagonal boron nitride superlatticesEntropic evidence for a Pomeranchuk effect in magic angle grapheneSpectroscopy of a Tunable Moiré System with a Correlated and Topological Flat BandEmergence of Chern Insulating States in Non-Magic Angle Twisted Bilayer GrapheneFragile Topology and Wannier ObstructionsFaithful tight-binding models and fragile topology of magic-angle bilayer grapheneMechanism for Anomalous Hall Ferromagnetism in Twisted Bilayer GrapheneOrigin of Mott Insulating Behavior and Superconductivity in Twisted Bilayer GrapheneOrigin of Magic Angles in Twisted Bilayer GrapheneModel for the metal-insulator transition in graphene superlattices and beyondSymmetry, Maximally Localized Wannier States, and a Low-Energy Model for Twisted Bilayer Graphene Narrow BandsMaximally Localized Wannier Orbitals and the Extended Hubbard Model for Twisted Bilayer GrapheneUnconventional superconductivity in nearly flat bands in twisted bilayer grapheneCorrelated insulating phases of twisted bilayer graphene at commensurate filling fractions: A Hartree-Fock studyGround State and Hidden Symmetry of Magic-Angle Graphene at Even Integer FillingHybrid Wannier Chern bands in magic angle twisted bilayer graphene and the quantized anomalous Hall effectNature of the Correlated Insulator States in Twisted Bilayer GraphenePseudo Landau level representation of twisted bilayer graphene: Band topology and implications on the correlated insulating phaseMachine learning enabled autonomous microstructural characterization in 3D samplesTheories for the correlated insulating states and quantum anomalous Hall effect phenomena in twisted bilayer grapheneBand structure and insulating states driven by Coulomb interaction in twisted bilayer grapheneNematic topological semimetal and insulator in magic-angle bilayer graphene at charge neutralityExact continuum model for low-energy electronic states of twisted bilayer grapheneKekulé spiral order at all nonzero integer fillings in twisted bilayer grapheneSkyrmion Superconductivity: DMRG evidence for a topological route to superconductivityCharged skyrmions and topological origin of superconductivity in magic-angle grapheneNon-Abelian Dirac node braiding and near-degeneracy of correlated phases at odd integer filling in magic-angle twisted bilayer grapheneEfficient simulation of moiré materials using the density matrix renormalization groupQuasi-flat-band physics in a two-leg ladder model and its relation to magic-angle twisted bilayer grapheneTwisted bilayer graphene. VI. An exact diagonalization study at nonzero integer fillingPossible correlated insulating states in magic-angle twisted bilayer graphene under strongly competing interactionsPhases of a phenomenological model of twisted bilayer grapheneExact Diagonalization for Magic-Angle Twisted Bilayer GrapheneRenormalization Group Study of Hidden Symmetry in Twisted Bilayer Graphene with Coulomb InteractionsStrong Coupling Phases of Partially Filled Twisted Bilayer Graphene Narrow BandsTwisted bilayer graphene. II. Stable symmetry anomalyTwisted bilayer graphene. III. Interacting Hamiltonian and exact symmetriesTwisted bilayer graphene. IV. Exact insulator ground states and phase diagramTwisted bilayer graphene. V. Exact analytic many-body excitations in Coulomb Hamiltonians: Charge gap, Goldstone modes, and absence of Cooper pairingFerromagnetism and its stability from the one-magnon spectrum in twisted bilayer grapheneTwo-dimensional Hubbard model: Numerical simulation studyNon-Fermi Liquid at ( 2 + 1 ) D Ferromagnetic Quantum Critical PointItinerant quantum critical point with fermion pockets and hotspotsValence Bond Orders at Charge Neutrality in a Possible Two-Orbital Extended Hubbard Model for Twisted Bilayer GrapheneCorrelation-Induced Insulating Topological Phases at Charge Neutrality in Twisted Bilayer GrapheneCorrelated insulating phases in the twisted bilayer graphene*Kekulé valence bond order in an extended Hubbard model on the honeycomb lattice with possible applications to twisted bilayer grapheneAntiferromagnetically ordered Mott insulator and d + i d superconductivity in twisted bilayer graphene: a quantum Monte Carlo studyElective-momentum ultrasize quantum Monte Carlo methodPhases of the ( 2 + 1 ) Dimensional SO(5) Nonlinear Sigma Model with Topological TermHalf-filled Landau levels: A continuum and sign-free regularization for three-dimensional quantum critical pointsRevealing fermionic quantum criticality from new Monte Carlo techniquesFermionic Monte Carlo study of a realistic model of twisted bilayer grapheneStochastic method for analytic continuation of quantum Monte Carlo dataIdentifying the maximum entropy method as a special limit of stochastic analytic continuationConstrained sampling method for analytic continuationUsing the average spectrum method to extract dynamics from quantum Monte Carlo simulationsNearly Deconfined Spinon Excitations in the Square-Lattice Spin- 1 / 2 Heisenberg AntiferromagnetDynamical Signature of Symmetry Fractionalization in Frustrated MagnetsDynamical signature of fractionalization at a deconfined quantum critical pointDynamics of Topological Excitations in a Model Quantum Spin IceThe bulk-corner correspondence of time-reversal symmetric insulatorsKosterlitz-Thouless melting of magnetic order in the triangular quantum Ising material TmMgGaO4Evidence of the Berezinskii-Kosterlitz-Thouless phase in a frustrated magnetAmplitude Mode in Quantum Magnets via Dimensional CrossoverVestigial anyon condensation in kagome quantum spin liquidsPseudogap and superconductivity emerging from quantum magnetic fluctuations: a Monte Carlo studyIdentification of non-Fermi liquid fermionic self-energy from quantum Monte Carlo dataPhase diagram of the spin- 1 2 Yukawa–Sachdev-Ye-Kitaev model: Non-Fermi liquid, insulator, and superconductorYukawa-SYK model and self-tuned quantum criticalityFermi arcs and pseudogap in a lattice model of a doped orthogonal metal
[1] Trambly de Laissardière G T, Mayou D, and Magaud L 2010 Nano Lett. 10 804
[2] Trambly de Laissardière G T, Mayou D, and Magaud L 2012 Phys. Rev. B 86 125413
[3] Bistritzer R and MacDonald A H 2011 Proc. Natl. Acad. Sci. USA 108 12233
[4] Rozhkov A V, Sboychakov A O, Rakhmanov A L, and Nori F 2016 Phys. Rep. 648 1
[5] Lopes dos Santos J M B, Peres N M R, and Neto A H C 2007 Phys. Rev. Lett. 99 256802
[6] Lopes dos Santos J M B, Peres N M R, and Neto A H C 2012 Phys. Rev. B 86 155449
[7] Cao Y, Fatemi V, Fang S, Watanabe K, Taniguchi T, Kaxiras E, and Jarillo-Herrero P 2018 Nature 556 43
[8] Cao Y, Fatemi V, Demir A, Fang S, Tomarken S L, Luo J Y, Sanchez-Yamagishi J D, Watanabe K, Taniguchi T, Ashoori Ray C, and Jarillo-Herrero P 2018 Nature 556 80
[9] Chen G, Sharpe A L, Fox E J, Zhang Y H, Wang S, Jiang L, Lyu B, Li H, Watanabe K, Taniguchi T, Shi Z, Senthil T, Goldhaber-Gordon D, Zhang Y, and Wang F 2020 Nature 579 56
[10] Kerelsky A, McGilly L J, Kennes D M, Xian L, Yankowitz M, Chen S, Watanabe K, Taniguchi T, Hone J, Dean C, Rubio A, and Pasupathy A N 2019 Nature 572 95
[11] Tomarken S L, Cao Y, Demir A, Watanabe K, Taniguchi T, Jarillo-Herrero P, and Ashoori R C 2019 Phys. Rev. Lett. 123 046601
[12] Lu X, Stepanov P, Yang W, Xie M, Aamir M A, Das I, Urgell C, Watanabe K, Taniguchi T, Zhang G et al. 2019 Nature 574 653
[13] Xie Y, Lian B, Jäck B, Liu X, Chiu C L, Watanabe K, Taniguchi T, Bernevig B A, and Yazdani A 2019 Nature 572 101
[14] Shen C, Chu Y, Wu Q, Li N, Wang S, Zhao Y, Tang J, Liu J, Tian J, Watanabe K, Taniguchi T, Yang R, Meng Z Y, Shi D, Yazyev O V, and Zhang G 2020 Nat. Phys. 16 520
[15] Nuckolls K P, Myungchul O, Wong D, Lian B, Watanabe K, Taniguchi T, Bernevig B A, and Yazdani A 2020 Nature 588 610
[16] Pierce A T, Xie Y, Park J M, Khalaf E, Lee S H, Cao Y, Parker D E, Forrester P R, Chen S, Watanabe K, Taniguchi T, Vishwanath A, Jarillo-Herrero P, and Yacoby A 2021 arXiv:2101.04123 [cond-mat.mes-hall]
[17] Moriyama S, Morita Y, Komatsu K, Endo K, Iwasaki T, Nakaharai S, Noguchi Y, Wakayama Y, Watanabe E, Tsuya D, Watanabe K, and Taniguchi T 2019 arXiv:1901.09356 [cond-mat.supr-con]
[18] Rozen A, Park J M, Zondiner U, Cao Y, Rodan-Legrain D, Taniguchi T, Watanabe K, Oreg Y, Stern A, Berg E, Jarillo-Herrero P, and Ilani S 2020 arXiv:2009.01836 [cond-mat.mes-hall]
[19] Liu X, Chiu C L, Lee J Y, Farahi G, Watanabe K, Taniguchi T, Vishwanath A, and Yazdani A 2020 arXiv:2008.07552 [cond-mat.mes-hall]
[20] Shen C, Ying J, Liu L, Liu J, Li N, Wang S, Tang J, Zhao Y, Chu Y, Watanabe K, Taniguchi T, Yang R, Shi D, Qu F, Lu L, Yang W, and Zhang G 2021 Chin. Phys. Lett. 38 047301
[21] Po H C, Watanabe H, and Vishwanath A 2018 Phys. Rev. Lett. 121 126402
[22] Po H C, Zou L, Senthil T, and Vishwanath A 2019 Phys. Rev. B 99 195455
[23] Bultinck N, Chatterjee S, and Zaletel M P 2020 Phys. Rev. Lett. 124 166601
[24] Po H C, Zou L, Vishwanath A, and Senthil T 2018 Phys. Rev. X 8 031089
[25] Tarnopolsky G, Kruchkov A J, and Vishwanath A 2019 Phys. Rev. Lett. 122 106405
[26] Yuan N F Q and Fu L 2018 Phys. Rev. B 98 045103
[27] Kang J and Vafek O 2018 Phys. Rev. X 8 031088
[28] Koshino M, Yuan N F Q, Koretsune T, Ochi M, Kuroki K, and Fu L 2018 Phys. Rev. X 8 031087
[29] Roy B and Juričić V 2019 Phys. Rev. B 99 121407
[30] Zhang Y, Jiang K, Wang Z, and Zhang F 2020 Phys. Rev. B 102 035136
[31] Bultinck N, Khalaf E, Liu S, Chatterjee S, Vishwanath A, and Zaletel M P 2020 Phys. Rev. X 10 031034
[32] Hejazi K, Chen X, and Balents L 2021 Phys. Rev. Res. 3 013242
[33] Xie M and MacDonald A H 2020 Phys. Rev. Lett. 124 097601
[34] Liu J, Liu J, and Dai X 2019 Phys. Rev. B 99 155415
[35] Liu J and Dai X 2020 npj Comput. Mater. 6 1
[36] Liu J and Dai X 2021 Phys. Rev. B 103 035427
[37] Cea T and Guinea F 2020 Phys. Rev. B 102 045107
[38] Liu S, Khalaf E, Lee J Y, and Vishwanath A 2021 Phys. Rev. Res. 3 013033
[39] Carr S, Fang S, Zhu Z, and Kaxiras E 2019 Phys. Rev. Res. 1 013001
[40] Kwan Y H, Wagner G, Soejima T, Zaletel M P, Simon S H, Parameswaran S A, and Bultinck N 2021 arXiv:2105.05857 [cond-mat.str-el]
[41] Chatterjee S, Ippoliti M, and Zaletel M P 2020 arXiv:2010.01144 [cond-mat.str-el]
[42] Khalaf E, Chatterjee S, Bultinck N, Zaletel M P, and Vishwanath A 2021 Sci. Adv. 7 eabf5299
[43] Kang J and Vafek O 2020 Phys. Rev. B 102 035161
[44] Soejima T, Parker D E, Bultinck N, Hauschild J, and Zaletel M P 2020 Phys. Rev. B 102 205111
[45] Huang Y, Hosur P, and Pal H K 2020 Phys. Rev. B 102 155429
[46] Xie F, Cowsik A, Song Z D, Lian B, Bernevig B A, and Regnault N 2021 Phys. Rev. B 103 205416
[47] Ochi M, Koshino M, and Kuroki K 2018 Phys. Rev. B 98 081102
[48] Dodaro J F, Kivelson S A, Schattner Y, Sun X Q, and Wang C 2018 Phys. Rev. B 98 075154
[49] Potasz P, Xie M, and MacDonald A H 2021 arXiv:2102.02256 [cond-mat.str-el]
[50] Vafek O and Kang J 2020 Phys. Rev. Lett. 125 257602
[51] Kang J and Vafek O 2019 Phys. Rev. Lett. 122 246401
[52] Song Z D, Lian B, Regnault N, and Bernevig B A 2021 Phys. Rev. B 103 205412
[53] Bernevig B A, Song Z D, Regnault N, and Lian B 2021 Phys. Rev. B 103 205413
[54] Lian B, Song Z D, Regnault N, Efetov D K, Yazdani A, and Bernevig B A 2021 Phys. Rev. B 103 205414
[55] Bernevig B A, Lian B, Cowsik A, Xie F, Regnault N, and Song Z D 2021 Phys. Rev. B 103 205415
[56] Alavirad Y and Sau J 2020 Phys. Rev. B 102 235123
[57] Hirsch J E 1985 Phys. Rev. B 31 4403
[58] Xu X Y, Sun K, Schattner Y, Berg E, and Meng Z Y 2017 Phys. Rev. X 7 031058
[59] Liu Z H, Pan G, Xu X Y, Sun K, and Meng Z Y 2019 Proc. Natl. Acad. Sci. USA 116 16760
[60] Liao Y D, Meng Z Y, and Xu X Y 2019 Phys. Rev. Lett. 123 157601
[61] Liao Y D, Kang J, Breiø C N, Xu X Y, Wu H Q, Andersen B M, Fernandes R M, and Meng Z Y 2021 Phys. Rev. X 11 011014
[62] Liao Y D, Xu X Y, Meng Z Y, and Kang J 2021 Chin. Phys. B 30 017305
[63] Xu X Y, Law K T, and Lee P A 2018 Phys. Rev. B 98 121406
[64] Huang T, Zhang L, and Ma T 2019 Sci. Bull. 64 310
[65] Liu Z H, Xu X Y, Qi Y, Sun K, and Meng Z Y 2019 Phys. Rev. B 99 085114
[66] Wang Z, Zaletel M P, Mong R S K, and Assaad F F 2021 Phys. Rev. Lett. 126 045701
[67] Ippoliti M, Mong R S K, Assaad F F, and Zaletel M P 2018 Phys. Rev. B 98 235108
[68]The symmetry properties of the form factor and the proofs of the sign structure of the fermion determinant, the QMC measurements and brief description of the stochastic analytic continuation, are presented in the Supplemental Material.
[69]Assaad F and Evertz H 2008 World-Line and Determinantal Quantum Monte Carlo Methods for Spins, Phonons and Electrons, in Computational Many-Particle Physics, ed Fehske H, Schneider R and A. Weiße (Berlin: Springer) pp 277–356
[70] Xu X Y, Liu Z H, Pan G, Qi Y, Sun K, and Meng Z Y 2019 J. Phys.: Condens. Matter 31 463001
[71]Lee J Y, Hofmann J, Khalaf E, Vishwanath A, and Berg E 2021 APS March Meeting S43 00013
[72] Hofmann J S, Khalaf E, Vishwanath A, Berg E, and Lee J Y 2021 arXiv:2105.12112 [cond-mat.str-el]
[73] Sandvik A W 1998 Phys. Rev. B 57 10287
[74] Beach K 2004 arXiv:cond-mat/0403055 [cond-mat.str-el]
[75] Sandvik A W 2016 Phys. Rev. E 94 063308
[76] Syljuåsen O F 2008 Phys. Rev. B 78 174429
[77] Shao H, Qin Y Q, Capponi S, Chesi S, Meng Z Y, and Sandvik A W 2017 Phys. Rev. X 7 041072
[78] Sun G Y, Wang Y C, Fang C, Qi Y, Cheng M, and Meng Z Y 2018 Phys. Rev. Lett. 121 077201
[79] Ma N, Sun G Y, You Y Z, Xu C, Vishwanath A, Sandvik A W, and Meng Z Y 2018 Phys. Rev. B 98 174421
[80] Huang C J, Deng Y, Wan Y, and Meng Z Y 2018 Phys. Rev. Lett. 120 167202
[81] Yan Z, Wang Y C, Ma N, Qi Y, and Meng Z Y 2021 npj Quantum Mater. 6 1
[82] Li H, Liao Y D, Chen B B, Zeng X T, Sheng X L, Qi Y, Meng Z Y, and Li W 2020 Nat. Commun. 11 1111
[83] Hu Z, Ma Z, Liao Y D, Li H, Ma C, Cui Y, Shangguan Y, Huang Z, Qi Y, Li W et al. 2020 Nat. Commun. 11 5631
[84] Zhou C, Yan Z, Wu H Q, Sun K, Starykh O A, and Meng Z Y 2020 arXiv:2007.12715 [cond-mat.str-el]
[85] Wang Y C, Yan Z, Wang C, Qi Y, and Meng Z Y 2021 Phys. Rev. B 103 014408
[86] Jiang W, Liu Y, Klein A, Wang Y, Sun K, Chubukov A V, and Meng Z Y 2021 arXiv:2105.03639 [cond-mat.str-el]
[87] Xu X Y, Klein A, Sun K, Chubukov A V, and Meng Z Y 2020 npj Quantum Mater. 5 65
[88] Wang W, Davis A, Pan G, Wang Y, and Meng Z Y 2021 Phys. Rev. B 103 195108
[89] Pan G, Wang W, Davis A, Wang Y, and Meng Z Y 2021 Phys. Rev. Res. 3 013250
[90] Chen C, Yuan T, Qi Y, and Meng Z Y 2021 Phys. Rev. B 103 165131