$${H}_{0}=\t \\sum\\limits_{s,\\eta,{{{\\bf{k}}}},{{{\\bf{G}}}},X,{{{{\\bf{G}}}}}^{{\\prime} },{X}^{{\\prime} }}{d}_{s,\\eta,{{{\\bf{k}}}},{{{\\bf{G}}}},X}^{{{\\dagger}} }{H}_{{{{\\bf{G}}}},X;{{{{\\bf{G}}}}}^{{\\prime} },{X}^{{\\prime} }}^{s,\\eta,{{{\\bf{k}}}}}{d}_{s,\\eta,{{{\\bf{k}}}},{{{{\\bf{G}}}}}^{{\\prime} },{X}^{{\\prime} }}\\\\=\t \\sum\\limits_{s,\\eta,{{{\\bf{k}}}}}{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}^{{{\\dagger}} }{H}^{s,\\eta,{{{\\bf{k}}}}}{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}.$$<\/p>\n
\n (6)\n <\/p>\n
Here \\({d}_{s,\\eta,{{{\\bf{k}}}},{{{\\bf{G}}}},X}^{{{\\dagger}} }\\) is the creation operator for spin s, valley \u03b7, momentum k in mBZ, momentum folding G to account for the contributions from extended Brillouin zones, layer (1, 2) and sublattice (A, B), and \\({\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}^{{{\\dagger}} }\\) is the vector of \\({d}_{s,\\eta,{{{\\bf{k}}}},{{{\\bf{G}}}},X}^{{{\\dagger}} }\\). We use X\u2009\u2208\u2009{1A, 1B, 2A, 2B} to denote the layer and sublattice degrees of freedom. So, Hs,\u03b7,k is a matrix with the dimension of NGNX, which is actually the BM model7<\/a>, and with the sub-blocks as<\/p>\n $$\\begin{array}{l}{H}_{{{{\\bf{G}}}},{{{{\\bf{G}}}}}^{{\\prime} }}^{s,\\eta,{{{\\bf{k}}}}} \\hfill \\\\={\\delta }_{{{{\\bf{G}}}},{{{{\\bf{G}}}}}^{{\\prime} }}\\hslash {v}_{{{{\\rm{F}}}}}\\left(\\begin{array}{ll}-({{{\\bf{k}}}}+{{{\\bf{G}}}}-{{{{\\bf{K}}}}}_{1}^{\\eta })\\cdot {\\sigma }^{\\eta }&0\\\\ 0&-({{{\\bf{k}}}}+{{{\\bf{G}}}}-{{{{\\bf{K}}}}}_{2}^{\\eta })\\cdot {\\sigma }^{\\eta }\\end{array}\\right)\\\\+\\left(\\begin{array}{ll}0&{T}_{1}^{\\eta }\\\\ {{T}_{2}^{\\eta }}^{{{\\dagger}} }&0\\end{array}\\right) \\hfill \\end{array}$$<\/p>\n \n (7)\n <\/p>\n with<\/p>\n $$\\begin{array}{l}{T}_{l}^{\\eta }=\\left(\\begin{array}{ll}{u}_{0}&{u}_{1}\\\\ {u}_{1}&{u}_{0}\\end{array}\\right){\\delta }_{{{{\\bf{G}}}},{{{{\\bf{G}}}}}^{{\\prime} }}+\\left(\\begin{array}{ll}{u}_{0}&{u}_{1}{{{{\\rm{e}}}}}^{-i\\frac{2\\pi }{3}\\eta }\\\\ {u}_{1}{{{{\\rm{e}}}}}^{i\\frac{2\\pi }{3}\\eta }&{u}_{0}\\end{array}\\right){\\delta }_{{{{\\bf{G}}}},{{{{\\bf{G}}}}}^{{\\prime} }+{(-1)}^{l}\\eta {{{{\\bf{G}}}}}_{1}}\\\\+\\left(\\begin{array}{ll}{u}_{0}&{u}_{1}{{{{\\rm{e}}}}}^{i\\frac{2\\pi }{3}\\eta }\\\\ {u}_{1}{{{{\\rm{e}}}}}^{-i\\frac{2\\pi }{3}\\eta }&{u}_{0}\\end{array}\\right){\\delta }_{{{{\\bf{G}}}},{{{{\\bf{G}}}}}^{{\\prime} }+{(-1)}^{l}\\eta \\left({{{{\\bf{G}}}}}_{1}+{{{{\\bf{G}}}}}_{2}\\right)},\\hfill \\end{array}$$<\/p>\n \n (8)\n <\/p>\n where \\({{{\\bf{G}}}},{{{{\\bf{G}}}}}^{{\\prime} }\\in \\{{n}_{1}{{{{\\bf{G}}}}}_{1}+{n}_{2}{{{{\\bf{G}}}}}_{2}\\}\\) defining the lattices of mBZs, with n1 and n2 being integers, and a cutoff, \u2223G\u2223, \\(| {{{{\\bf{G}}}}}^{{\\prime} }| \\le 6| {{{{\\bf{G}}}}}_{1}|\\), is applied. vF is the Fermi velocity and \\(\\hslash {v}_{{{{\\rm{F}}}}}\/(\\sqrt{3}a)=2.37745\\) eV. \\({\\sigma }^{\\eta }=\\left(\\eta {\\sigma }_{x},{\\sigma }_{y}\\right)\\), where \u03c3x and \u03c3y are Pauli matrices in sublattice space. The first part of \\({H}_{{{{\\bf{G}}}},{{{{\\bf{G}}}}}^{{\\prime} }}^{s,\\eta,{{{\\bf{k}}}}}\\) stands for intra-layer hopping of the top (1) layer and the bottom (2) layer, respectively, and the second part stands for inter-layer hopping u0 and u1, which are actually the moir\u00e9 potentials.<\/p>\n Fig. 5: Momentum-space configuration of BZs and mBZs. Twisting of graphene BZs as blue and orange hexagons leads to the mBZs as green hexagon or black rhombus with lattices of G1 and G2. \u0398 is the twist angle and \\({{{{\\bf{K}}}}}_{l}^{\\eta }\\) is a Dirac point of valley \u03b7 and layer l. \u0393 and M are two high-symmetric points.<\/p>\n The dispersion of H0 for several twist angles are shown in Supplementary Fig.\u20092<\/a> with u0\u2009=\u200960\u2009meV, where enlarging of the bandwidth of two low-energy bands is seen with increasing \u0398 from MA, and the two low-energy bands are always isolated from the remote bands by \u00a0~ 60\u2009meV. For the cases considered in this study, the minimum gap between low-energy bands and remote bands is not less than 20 meV while the maximum Coulomb potential is less than 6\u2009meV, as shown in Supplementary Fig.\u20093<\/a>. Thus, the low-energy physics can be well captured by these two low-energy bands.<\/p>\n The single-gated Coulomb interaction after Fouriers transformation is<\/p>\n $${H}_{{{{\\rm{I}}}}}=\\frac{1}{2\\Omega }\\sum \\limits_{{{{\\bf{q}}}}\\in {{{\\rm{mBZ}}}},{{{\\bf{G}}}}}V({{{\\bf{q}}}}+{{{\\bf{G}}}})\\delta {\\rho }_{{{{\\bf{q}}}}+{{{\\bf{G}}}}}\\delta {\\rho }_{-{{{\\bf{q}}}}-{{{\\bf{G}}}}},$$<\/p>\n \n (9)\n <\/p>\n where<\/p>\n $$\\begin{array}{l}\\delta {\\rho }_{{{{\\bf{q}}}}+{{{\\bf{G}}}}} \\hfill \\\\=\\sum\\limits_{s,\\eta,{{{\\bf{k}}}},{{{{\\bf{G}}}}}^{{\\prime} },{X}^{{\\prime} }}\\left({d}_{s,\\eta,{{{\\bf{k}}}},{{{{\\bf{G}}}}}^{{\\prime} },{X}^{{\\prime} }}^{{{\\dagger}} }{d}_{s,\\eta,{{{\\bf{k}}}}+{{{\\bf{q}}}}+{{{\\bf{G}}}},{{{{\\bf{G}}}}}^{{\\prime} },{X}^{{\\prime} }}-\\frac{\\nu+4}{8}{\\delta }_{{{{\\bf{q}}}},0}{\\delta }_{{{{\\bf{G}}}},0}\\right)\\\\=\\sum\\limits_{s,\\eta,{{{\\bf{k}}}}}\\left({\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}^{{{\\dagger}} }{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}+{{{\\bf{q}}}}+{{{\\bf{G}}}}}-\\frac{\\nu+4}{8}{\\delta }_{{{{\\bf{q}}}},0}{\\delta }_{{{{\\bf{G}}}},0}\\right).\\hfill\\end{array}$$<\/p>\n \n (10)\n <\/p>\n Here \u03b4\u03c1q+G is the electron density operator with respect to (\u03bd\u2009+\u20094)\/8 filling, and \u03bd is the filling parameter with \u03bd\u2009=\u20090 corresponding to the charge neutrality. \\(\\Omega=N| {{{{\\bf{a}}}}}_{{{{\\rm{M1}}}}}| | {{{{\\bf{a}}}}}_{{{{\\rm{M2}}}}}| \\sqrt{3}\/2\\) is the total area of mori\u00e9 superlattice in real space, while aM1 and aM2 are the lattice vectors of the real-space unit cell of the superlattice with \\(| {{{{\\bf{a}}}}}_{{{{\\rm{M1,M2}}}}}|=\\sqrt{3}a\/[2\\sin (\\Theta \/2)]\\). For the notational simplicity, from here on, q\u2009+\u2009G is denoted as \\({{{\\mathcal{Q}}}}\\).<\/p>\n Denoting the eigenvalues and eigenvectors of the BM-model matrix Hs,\u03b7,k in Eq.\u2009(7<\/a>) as \\({\\epsilon }_{{{{\\bf{k}}}},m}^{s,\\eta }\\) and \\(\\left\\vert {u}_{{{{\\bf{k}}}},m}^{s,\\eta }\\right\\rangle\\), \\({H}^{s,\\eta,{{{\\bf{k}}}}}\\left\\vert {u}_{{{{\\bf{k}}}},m}^{s,\\eta }\\right\\rangle={\\epsilon }_{{{{\\bf{k}}}},m}^{s,\\eta }\\left\\vert {u}_{{{{\\bf{k}}}},m}^{s,\\eta }\\right\\rangle\\), or<\/p>\n $${H}^{s,\\eta,{{{\\bf{k}}}}}{U}^{s,\\eta,{{{\\bf{k}}}}}={U}^{s,\\eta,{{{\\bf{k}}}}}{E}^{s,\\eta,{{{\\bf{k}}}}},$$<\/p>\n \n (11)\n <\/p>\n with each column of Us,\u03b7,k consisting of \\(\\left\\vert {u}_{{{{\\bf{k}}}},m}^{s,\\eta }\\right\\rangle\\) and Es,\u03b7,k is a diagonal matrix consisting of \\({\\epsilon }_{{{{\\bf{k}}}},m}^{s,\\eta }\\), then H0 can be projected to the band basis as<\/p>\n $${H}_{0}=\t \\sum\\limits_{s,\\eta,{{{\\bf{k}}}}}{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}^{{{\\dagger}} }{H}^{s,\\eta,{{{\\bf{k}}}}}{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}\\\\=\t \\sum\\limits_{s,\\eta,{{{\\bf{k}}}}}{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}^{{{\\dagger}} }{U}^{s,\\eta,{{{\\bf{k}}}}}{{U}^{s,\\eta,{{{\\bf{k}}}}}}^{{{\\dagger}} }{H}^{s,\\eta,{{{\\bf{k}}}}}{U}^{s,\\eta,{{{\\bf{k}}}}}{{U}^{s,\\eta,{{{\\bf{k}}}}}}^{{{\\dagger}} }{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}\\\\=\t \\sum\\limits_{s,\\eta,{{{\\bf{k}}}}}{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}^{{{\\dagger}} }{U}^{s,\\eta,{{{\\bf{k}}}}}{E}^{s,\\eta,{{{\\bf{k}}}}}{{U}^{s,\\eta,{{{\\bf{k}}}}}}^{{{\\dagger}} }{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}.$$<\/p>\n \n (12)\n <\/p>\n Defining the projection<\/p>\n $${\\vec{c}}_{s,\\eta,{{{\\bf{k}}}}}^{{{\\dagger}} }={\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}^{{{\\dagger}} }{U}^{s,\\eta,{{{\\bf{k}}}}},$$<\/p>\n \n (13)\n <\/p>\n where \\({\\vec{c}}_{s,\\eta,{{{\\bf{k}}}}}\\) is an operator vector with indices of G and X, with the element,<\/p>\n $${c}_{s,\\eta,{{{\\bf{k}}}},m}^{{{\\dagger}} }=\\sum\\limits_{{{{\\bf{G}}}},X}{d}_{s,\\eta,{{{\\bf{k}}}},{{{\\bf{G}}}},X}^{{{\\dagger}} }{\\left| {u}_{{{{\\bf{k}}}},m}^{s,\\eta }\\right\\rangle }_{{{{\\bf{G}}}},X},$$<\/p>\n \n (14)\n <\/p>\n then<\/p>\n $${H}_{0}=\t \\sum\\limits_{s,\\eta,{{{\\bf{k}}}}}{\\vec{c}}_{s,\\eta,{{{\\bf{k}}}}}^{{{\\dagger}} }{E}^{s,\\eta,{{{\\bf{k}}}}}{\\vec{c}}_{s,\\eta,{{{\\bf{k}}}}}\\\\=\t \\sum\\limits_{s,\\eta,{{{\\bf{k}}}},m}{\\epsilon }_{{{{\\bf{k}}}},m}^{s,\\eta }{c}_{s,\\eta,{{{\\bf{k}}}},m}^{{{\\dagger}} }{c}_{s,\\eta,{{{\\bf{k}}}},m}.$$<\/p>\n \n (15)\n <\/p>\n This is the projection from the plane-wave basis, d\u2020, to the band basis, c\u2020. To make the projection to two low-energy bands, m relating to these two bands are kept.<\/p>\n Similarly, for the interaction,<\/p>\n $${\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}^{{{\\dagger}} }{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}+{{{\\mathcal{Q}}}}}=\t {\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}}^{{{\\dagger}} }{U}_{{{{\\bf{k}}}}}^{s,\\eta }{{U}_{{{{\\bf{k}}}}}^{s,\\eta }}^{{{\\dagger}} }{U}_{{{{\\bf{k}}}}+{{{\\mathcal{Q}}}}}^{s,\\eta }{{U}_{{{{\\bf{k}}}}+{{{\\mathcal{Q}}}}}^{s,\\eta }}^{{{\\dagger}} }{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}+{{{\\mathcal{Q}}}}}\\\\=\t {\\vec{c}}_{s,\\eta,{{{\\bf{k}}}}}^{{{\\dagger}} }{{U}_{{{{\\bf{k}}}}}^{s,\\eta }}^{{{\\dagger}} }{U}_{{{{\\bf{k}}}}+{{{\\mathcal{Q}}}}}^{s,\\eta }{\\vec{c}}_{s,\\eta,{{{\\bf{k}}}}+{{{\\mathcal{Q}}}}}\\\\=\t \\sum\\limits_{m,n}{c}_{s,\\eta,{{{\\bf{k}}}},m}^{{{\\dagger}} }{\\left({{U}_{{{{\\bf{k}}}}}^{s,\\eta }}^{{{\\dagger}} }{U}_{{{{\\bf{k}}}}+{{{\\mathcal{Q}}}}}^{s,\\eta }\\right)}_{m,n}{c}_{s,\\eta,{{{\\bf{k}}}}+{{{\\mathcal{Q}}}},n}\\\\=\t \\sum\\limits_{m,n}{c}_{s,\\eta,{{{\\bf{k}}}},m}^{{{\\dagger}} }{\\lambda }_{m,n}^{s,\\eta }({{{\\bf{k}}}},{{{\\mathcal{Q}}}}){c}_{s,\\eta,{{{\\bf{k}}}}+{{{\\mathcal{Q}}}},n},$$<\/p>\n \n (16)\n <\/p>\n where the form factor is defined as<\/p>\n $$\\begin{array}{rcl}{\\lambda }_{m,n}^{s,\\eta }({{{\\bf{k}}}},{{{\\mathcal{Q}}}})&=&{\\left({{U}_{{{{\\bf{k}}}}}^{s,\\eta }}^{{{\\dagger}} }{U}_{{{{\\bf{k}}}}+{{{\\mathcal{Q}}}}}^{s,\\eta }\\right)}_{m,n}\\\\ &=&\\left\\langle {u}_{{{{\\bf{k}}}},m}^{s,\\eta }| {u}_{{{{\\bf{k}}}}+{{{\\mathcal{Q}}}},n}^{s,\\eta }\\right\\rangle .\\end{array}$$<\/p>\n \n (17)\n <\/p>\n Note that we have made \\({U}_{{{{\\bf{k}}}}}^{s,\\eta }={U}^{s,\\eta,{{{\\bf{k}}}}}\\) to make the expressions compact. So,<\/p>\n $$\\delta {\\rho }_{{{{\\mathcal{Q}}}}}=\\sum\\limits_{s,\\eta,{{{\\bf{k}}}}}\\left(\\sum\\limits_{m,n}{\\lambda }_{m,n}^{s,\\eta }({{{\\bf{k}}}},{{{\\mathcal{Q}}}}){c}_{s,\\eta,{{{\\bf{k}}}},m}^{{{\\dagger}} }{c}_{s,\\eta,{{{\\bf{k}}}}+{{{\\mathcal{Q}}}},n}-\\frac{\\nu+4}{8}{\\delta }_{{{{\\bf{q}}}},0}{\\delta }_{{{{\\bf{G}}}},0}\\right)$$<\/p>\n \n (18)\n <\/p>\n Moreover, since<\/p>\n $$\\begin{array}{rcl}{\\delta }_{{{{\\bf{G}}}},0}&=&\\sum\\limits_{m}\\left\\langle {u}_{{{{\\bf{k}}}},m}^{s,\\eta }| {u}_{{{{\\bf{k}}}}+{{{\\bf{G}}}},m}^{s,\\eta }\\right\\rangle \\\\ &=&\\sum\\limits_{m,n}{\\lambda }_{m,n}^{s,\\eta }({{{\\bf{k}}}},{{{\\bf{G}}}}){\\delta }_{m,n},\\end{array}$$<\/p>\n \n (19)\n <\/p>\n and further,<\/p>\n $${\\delta }_{{{{\\bf{q}}}},0}{\\delta }_{{{{\\bf{G}}}},0}=\\sum\\limits_{m,n}{\\lambda }_{m,n}^{s,\\eta }({{{\\bf{k}}}},{{{\\mathcal{Q}}}}){\\delta }_{{{{\\bf{q}}}},0}{\\delta }_{m,n},$$<\/p>\n \n (20)\n <\/p>\n so that<\/p>\n $$\\delta {\\rho }_{{{{\\mathcal{Q}}}}}=\\sum\\limits_{s,\\eta,{{{\\bf{k}}}},m,n}{\\lambda }_{m,n}^{s,\\eta }({{{\\bf{k}}}},{{{\\mathcal{Q}}}})\\left({c}_{s,\\eta,{{{\\bf{k}}}},m}^{{{\\dagger}} }{c}_{s,\\eta,{{{\\bf{k}}}}+{{{\\mathcal{Q}}}},n}-\\frac{\\nu+4}{8}{\\delta }_{{{{\\bf{q}}}},0}{\\delta }_{m,n}\\right).$$<\/p>\n \n (21)\n <\/p>\n Since the set of \\(\\{| {{{\\mathcal{Q}}}}| \\ne 0\\}\\) is inversely symmetric, then<\/p>\n $${H}_{{{{\\rm{I}}}}}=\t \\sum\\limits_{| {{{\\mathcal{Q}}}}| \\ne 0}\\frac{1}{2\\Omega }V({{{\\mathcal{Q}}}})\\delta {\\rho }_{{{{\\mathcal{Q}}}}}\\delta {\\rho }_{-{{{\\mathcal{Q}}}}}\\\\=\t \\sum\\limits_{{{{\\bf{Q}}}}}\\frac{1}{2\\Omega }V({{{\\bf{Q}}}})(\\delta {\\rho }_{{{{\\bf{Q}}}}}\\delta {\\rho }_{-{{{\\bf{Q}}}}}+\\delta {\\rho }_{-{{{\\bf{Q}}}}}\\delta {\\rho }_{{{{\\bf{Q}}}}})\\\\=\t \\sum\\limits_{{{{\\bf{Q}}}}}\\frac{1}{4\\Omega }V({{{\\bf{Q}}}})\\left({\\left(\\delta {\\rho }_{-{{{\\bf{Q}}}}}+\\delta {\\rho }_{{{{\\bf{Q}}}}}\\right)}^{2}-{\\left(\\delta {\\rho }_{-{{{\\bf{Q}}}}}-\\delta {\\rho }_{{{{\\bf{Q}}}}}\\right)}^{2}\\right)\\\\=\t \\sum\\limits_{{{{\\bf{Q}}}}}\\frac{1}{4\\Omega }V({{{\\bf{Q}}}})\\left({A}_{{{{\\bf{Q}}}}}^{2}-{B}_{{{{\\bf{Q}}}}}^{2}\\right),$$<\/p>\n \n (22)\n <\/p>\n where the set of {Q} is any half of \\(\\{| {{{\\mathcal{Q}}}}| \\ne 0\\}\\). H0 in Eq.\u2009(15<\/a>) and HI in Eq.\u2009(22<\/a>) are exactly the terms used in the Hamiltonian (1<\/a>) in the main text.<\/p>\n The KIVC order parameter is defined in the plane-wave basis as<\/p>\n $${O}_{{{{\\rm{KIVC}}}}}({{{\\bf{q}}}})\\equiv \t \\sum\\limits_{s,{{{\\bf{k}}}}}{\\vec{d}}_{s,{{{\\bf{k}}}}+{{{\\bf{q}}}}}^{{{\\dagger}} }{\\eta }_{x}{\\mu }_{y}{\\sigma }_{x}{\\vec{d}}_{s,{{{\\bf{k}}}}}\\\\=\t \\sum\\limits_{s,\\eta,{{{\\bf{k}}}}}{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}+{{{\\bf{q}}}}}^{{{\\dagger}} }{\\mu }_{y}{\\sigma }_{x}{\\vec{d}}_{s,-\\eta,{{{\\bf{k}}}}},$$<\/p>\n \n (23)\n <\/p>\n and likewise projected to the band basis as<\/p>\n $${O}_{{{{\\rm{KIVC}}}}}({{{\\bf{q}}}})\\\\=\t \\sum\\limits_{s,\\eta,{{{\\bf{k}}}}}{\\vec{d}}_{s,\\eta,{{{\\bf{k}}}}+{{{\\bf{q}}}}}^{{{\\dagger}} }{U}_{{{{\\bf{k}}}}+{{{\\bf{q}}}}}^{s,\\eta }{{U}_{{{{\\bf{k}}}}+{{{\\bf{q}}}}}^{s,\\eta }}^{{{\\dagger}} }{\\mu }_{y}{\\sigma }_{x}{U}_{{{{\\bf{k}}}}}^{s,-\\eta }{{U}_{{{{\\bf{k}}}}}^{s,-\\eta }}^{{{\\dagger}} }{\\vec{d}}_{s,-\\eta,{{{\\bf{k}}}}}\\\\=\t \\sum\\limits_{s,\\eta,{{{\\bf{k}}}}}{\\vec{c}}_{s,\\eta,{{{\\bf{k}}}}+{{{\\bf{q}}}}}^{{{\\dagger}} }{{U}_{{{{\\bf{k}}}}+{{{\\bf{q}}}}}^{s,\\eta }}^{{{\\dagger}} }{\\mu }_{y}{\\sigma }_{x}{U}_{{{{\\bf{k}}}}}^{s,-\\eta }{\\vec{c}}_{s,-\\eta,{{{\\bf{k}}}}}\\\\=\t \\sum\\limits_{s,\\eta,{{{\\bf{k}}}},m,n}{c}_{s,\\eta,{{{\\bf{k}}}}+{{{\\bf{q}}}},m}^{{{\\dagger}} }{\\left({{U}_{{{{\\bf{k}}}}+{{{\\bf{q}}}}}^{s,\\eta }}^{{{\\dagger}} }{\\mu }_{y}{\\sigma }_{x}{U}_{{{{\\bf{k}}}}}^{s,-\\eta }\\right)}_{m,n}{c}_{s,-\\eta,{{{\\bf{k}}}},n}\\\\=\t \\sum\\limits_{s,\\eta,{{{\\bf{k}}}},m,n}{c}_{s,\\eta,{{{\\bf{k}}}}+{{{\\bf{q}}}},m}^{{{\\dagger}} }\\langle {u}_{{{{\\bf{k}}}}+{{{\\bf{q}}}},m}^{s,\\eta }| {\\mu }_{y}{\\sigma }_{x}| {u}_{{{{\\bf{k}}}},n}^{s,-\\eta }\\rangle {c}_{s,-\\eta,{{{\\bf{k}}}},n}.$$<\/p>\n \n (24)\n <\/p>\n Similarly, the VP order parameter from the plane-wave basis to the band basis is<\/p>\n $${O}_{{{{\\rm{VP}}}}}({{{\\bf{q}}}}) \\equiv \t \\sum\\limits_{s,{{{\\bf{k}}}}}{\\vec{d}}_{s,{{{\\bf{k}}}}+{{{\\bf{q}}}}}^{{{\\dagger}} }{\\eta }_{z}{\\mu }_{0}{\\sigma }_{0}{\\vec{d}}_{s,{{{\\bf{k}}}}}\\\\=\t -\\sum\\limits_{s,\\eta,{{{\\bf{k}}}},m,n}\\eta {c}_{s,\\eta,{{{\\bf{k}}}}+{{{\\bf{q}}}},m}^{{{\\dagger}} }\\langle {u}_{{{{\\bf{k}}}}+{{{\\bf{q}}}},m}^{s,\\eta }| {\\mu }_{0}{\\sigma }_{0}| {u}_{{{{\\bf{k}}}},n}^{s,\\eta }\\rangle {c}_{s,\\eta,{{{\\bf{k}}}},n}\\\\=\t -\\sum\\limits_{s,\\eta,{{{\\bf{k}}}},m,n}\\eta {c}_{s,\\eta,{{{\\bf{k}}}}+{{{\\bf{q}}}},m}^{{{\\dagger}} }\\langle {u}_{{{{\\bf{k}}}}+{{{\\bf{q}}}},m}^{s,\\eta }| {u}_{{{{\\bf{k}}}},n}^{s,\\eta }\\rangle {c}_{s,\\eta,{{{\\bf{k}}}},n}.$$<\/p>\n \n (25)\n <\/p>\n Gauge fixing<\/p>\n As can be seen from Eqs. (12<\/a>), (16<\/a>), the projections from the plane-wave basis to the band basis are gauge-independent within each spin-valley flavor, so that it is unnecessary to fix gauges within each spin-valley flavor. For different valleys, we fix \\({{{\\mathcal{PT}}}}\\), which is \u03b7x\u03bcy\u03c3x in plane-wave basis or \u03b7yny in the band basis with ny acting on bands. \\({{{\\mathcal{PT}}}}\\) leads to a complex conjugate between valleys.<\/p>\n In addition, to use the periodic condition, cs,\u03b7,k+G,m\u2009=\u2009cs,\u03b7,k,m, we impose<\/p>\n $${\\left\\vert {u}_{{{{\\bf{k}}}}+{{{\\bf{G}}}},m}^{s,\\eta }\\right\\rangle }_{{{{{\\bf{G}}}}}^{{\\prime} },X}={\\left\\vert {u}_{{{{\\bf{k}}}},m}^{s,\\eta }\\right\\rangle }_{{{{{\\bf{G}}}}}^{{\\prime} }-{{{\\bf{G}}}},X}.$$<\/p>\n \n (26)\n <\/p>\n The CFMC algorithm<\/p>\n The Hamiltonian in Eq.\u2009(1<\/a>) has been solved with momentum-space QMC method with discretized auxiliary fields38<\/a>,39<\/a>. Similar to other determinantal QMC methods, it applies the Hubbard-Stratonovich transformation to decouple the Coulomb interaction, following the Trotter decomposition127<\/a>,128<\/a>,129<\/a>,130<\/a>,131<\/a>. Due to \\({{{\\mathcal{PT}}}}\\) symmetry, the system at charge neutrality is free of sign problem38<\/a>,39<\/a>. The simulation temperature is set to be 1\/T\u00a0~\u00a0L with T\u00a0= 0.5 meV for L\u00a0= 6.<\/p>\n Original version of the momentum-space QMC method38<\/a>,42<\/a> uses the discrete version of the Hubbard-Stratonovich transformation following a Trotter decomposition:<\/p>\n $$\\exp \\left(-{\\alpha }_{1}({{{\\bf{Q}}}}){A}_{{{{\\bf{Q}}}}}^{2}\\right)=\t \\frac{1}{4}\\sum\\limits_{{l}_{{{{\\bf{Q}}}},1}}\\gamma \\left({l}_{{{{\\bf{Q}}}},1}\\right)\\exp \\left({{{\\rm{i}}}}\\xi \\left({l}_{{{{\\bf{Q}}}},1}\\right)\\sqrt{{\\alpha }_{1}({{{\\bf{Q}}}})}{A}_{{{{\\bf{Q}}}}}\\right)\\\\ \t+{{{\\mathcal{O}}}}\\left({{\\alpha }_{1}({{{\\bf{Q}}}})}^{4}\\right),$$<\/p>\n \n (27)\n <\/p>\n and<\/p>\n $$\\exp \\left({\\alpha }_{1}({{{\\bf{Q}}}}){B}_{{{{\\bf{Q}}}}}^{2}\\right)=\t \\frac{1}{4}\\sum\\limits_{{l}_{{{{\\bf{Q}}}},2}}\\gamma \\left({l}_{{{{\\bf{Q}}}},2}\\right)\\exp \\left(\\xi \\left({l}_{{{{\\bf{Q}}}},2}\\right)\\sqrt{{\\alpha }_{1}({{{\\bf{Q}}}})}{B}_{{{{\\bf{Q}}}}}\\right)\\\\ \t+{{{\\mathcal{O}}}}\\left({{\\alpha }_{1}({{{\\bf{Q}}}})}^{4}\\right).$$<\/p>\n \n (28)\n <\/p>\n Here lQ,1(2)\u2009=\u2009\u00b1\u00a01,\u00a0\u00b1\u00a02 is a discrete auxiliary field, and the coeficients, \\(\\gamma (\\pm 1)=1+\\sqrt{6}\/3\\), \\(\\gamma (\\pm 2)=1-\\sqrt{6}\/3\\), \\(\\xi (\\pm 1)=\\pm \\sqrt{6-2\\sqrt{6}}\\), and \\(\\xi (\\pm 2)=\\pm \\sqrt{6+2\\sqrt{6}}\\) are taken from the eighth-order approximation as in 129<\/a>. \u03b11(Q)\u2009=\u2009\u0394\u03c4V(Q)\/4\u03a9\u2009<\u20098\u2009\u00d7\u200910\u22123 with the imaginary time \u03b2\u2009=\u20091\/(kBT) evenly sliced into N\u03c4 pieces, as \u0394\u03c4\u2009=\u2009\u03b2\/N\u03c4. The partition function, \\(Z={{{\\rm{Tr}}}}\\{\\exp \\left(-\\beta H\\right)\\}\\), yields<\/p>\n $$Z={\\sum}_{\\{{l}_{\\tau,{{{\\bf{Q}}}},1},{l}_{\\tau,{{{\\bf{Q}}}},2}\\}}{{{\\rm{Tr}}}}\\,\\{{U}_{C}\\}$$<\/p>\n \n (29)\n <\/p>\n with<\/p>\n $${U}_{C}=\t {\\prod }_{\\tau=\\bigtriangleup \\tau }^{\\beta }\\exp \\left(-\\bigtriangleup \\tau {H}_{0}\\right)\\prod\\limits_{{{{\\bf{Q}}}}}\\frac{1}{16}\\gamma \\left({l}_{\\tau,{{{\\bf{Q}}}},1}\\right)\\gamma \\left({l}_{\\tau,Q,2}\\right)\\\\ \t \\times \\exp \\left({{{\\rm{i}}}}\\xi \\left({l}_{\\tau,{{{\\bf{Q}}}},1}\\right)\\sqrt{{\\alpha }_{1}({{{\\bf{Q}}}})}{A}_{{{{\\bf{Q}}}}}\\right)\\\\ \t \\times \\exp \\left(\\xi \\left({l}_{\\tau,{{{\\bf{Q}}}},2}\\right)\\sqrt{{\\alpha }_{1}({{{\\bf{Q}}}})}{B}_{{{{\\bf{Q}}}}}\\right)\\\\ \t+{{{\\mathcal{O}}}}(8\\times 1{0}^{-3}),$$<\/p>\n \n (30)\n <\/p>\n where C stands for a configuration of the auxiliary fields, and \\({{{\\mathcal{O}}}}(8\\times 1{0}^{-3})\\) stands for an error less than 8\u2009\u00d7\u200910\u22123, which is usually 5\u2009\u00d7\u200910\u22123, of the exact value.<\/p>\n With the partition function at hand, measurements can be done by sampling the discrete auxiliary fields with the Metropolis-Hastings scheme. Local update of one auxiliary field lQ,1(2) yields an acceptance rate around 0.4, and cluster update of multiple fields yields a vanishing acceptance rate. Each local update of a field l\u03c4,Q,1(2) leads to the computation of matrix determinant with the complexity \\({{{\\mathcal{O}}}}({N}^{3})\\), due to that AQ or BQ is non-diagonal. Combining with the absence of the efficient cluster update scheme for discrete fields, it leads to an overall complexity of at least \\({{{\\mathcal{O}}}}(\\beta {N}^{4})\\) in the original momentum-space QMC for a sweep to update 2N\u03c4NQ fields. Note that the number of the combinations (Q,\u00a01(2)) is 2NQ\u2009\u2248\u20093N (the cutoff at \u2223G1\u2223 makes the number of momentum exchanges roughly 3 times of the number of k points in mBZ). Therefore, the previous studies on TBG by QMC were heavily constrained by the accessible system sizes, with the largest one being L\u2009=\u2009939<\/a>,42<\/a>.<\/p>\n To reduce the computational complexity, we develop a continuous-field Monte Carlo (CFMC) scheme to carry out cluster updates of multiple auxiliary fields at once. First of all, we switch to continuous auxiliary fields using the Gaussian integrals<\/p>\n $${{{{\\rm{e}}}}}^{\\frac{1}{2}{\\alpha }_{2}{A}^{2}}=\\frac{1}{\\sqrt{2\\pi }}\\int_{-\\infty }^{\\infty }d\\phi {{{{\\rm{e}}}}}^{-\\frac{1}{2}{\\phi }^{2}}{{{{\\rm{e}}}}}^{-\\phi \\sqrt{{\\alpha }_{2}}A},$$<\/p>\n \n (31)\n <\/p>\n to rewrite the partition function of Eq.\u2009(29<\/a>) in the form<\/p>\n $$Z=\t \\int\\left(\\prod\\limits_{\\tau,{{{\\bf{Q}}}}}d{\\phi }_{\\tau,{{{\\bf{Q}}}},1}d{\\phi }_{\\tau,{{{\\bf{Q}}}},2}\\right){{{{\\rm{e}}}}}^{-\\frac{1}{2}\\sum\\limits_{\\tau,Q}\\left({\\phi }_{\\tau,{{{\\bf{Q}}}},1}^{2}+{\\phi }_{\\tau,{{{\\bf{Q}}}},2}^{2}\\right)}\\\\ \t \\times {{{\\rm{Tr}}}}\\left(\\prod\\limits_{\\tau }{{{{\\rm{e}}}}}^{{{{\\rm{i}}}}\\sum\\limits_{{{{\\bf{Q}}}}}\\left(-{\\phi }_{\\tau,{{{\\bf{Q}}}},1}\\sqrt{{\\alpha }_{2}({{{\\bf{Q}}}})}{A}_{{{{\\bf{Q}}}}}+{{{\\rm{i}}}}{\\phi }_{\\tau,{{{\\bf{Q}}}},2}\\sqrt{{\\alpha }_{2}({{{\\bf{Q}}}})}{B}_{{{{\\bf{Q}}}}}\\right)}{{{{\\rm{e}}}}}^{-\\bigtriangleup \\tau {H}_{0}}\\right)\\\\ \t+{{{\\mathcal{O}}}}(8\\times 1{0}^{-3}),$$<\/p>\n \n (32)\n <\/p>\n where \u03b12(Q)\u00a0=\u00a0\u0394\u03c4V(Q)\/2\u03a9. Note that the Gaussian decoupling is exact, so that the error of Z here comes from the Trotter decomposition, same as in Eq.\u2009(29<\/a>).<\/p>\n We employ the same Metropolis-Hastings scheme for the updates of the continuous fields \u03d5\u03c4,Q,1(2). Small updates of \u03d5\u03c4,Q,1(2) would make an acceptance rate around one, while large updates would lead to a vanishing acceptance rate. The changes of \u03d5\u03c4,Q,1(2), \u0394\u03d5\u03c4,Q,1(2), can be proposed from a Gaussian distribution or from the Hamiltonian\u2019s dynamics.<\/p>\n Gaussian update scheme<\/p>\n For the Gaussian proposal, we tune the magnitudes of \u0394\u03d5\u03c4,Q,1(2) to make the acceptance rate around 0.4. The difference between the new and the old value of the auxiliary field is written as \u0394\u03d5\u03c4,Q,1(2)\u00a0=\u00a0rx with r being the control parameter to tune the magnitudes of \u0394\u03d5\u03c4,Q,1(2) and x generated from a standard normal distribution, \\({{{\\mathcal{N}}}}(0,1)\\). The r is updated by r\u2009\u2192\u2009r\u2009+\u2009(Ra\u2009\u2212\u20090.4)\u0394r where \u0394r is a positive constant and Ra is the measured acceptance rate. Ra smaller than 0.4 means \u0394\u03d5\u03c4,Q,1(2) are larger than the desirable ones, and thus r is decreased, and vice versa. Ra\u2009~\u20090.4 should be achieved before the Markov process reaches equilibrium distribution.<\/p>\n All \u03d5\u03c4,Q,1(2) of a specific imaginary time \u03c4, with the number of \u00a0\u2248\u20093N in total, are updated each time. Thus, only one calculation of matrix determinant with complexity \\({{{\\mathcal{O}}}}({N}^{3})\\) is needed for each imaginary time \u03c4.<\/p>\n The full protocol of CFMC with Gaussian proposals can be described as follows<\/p>\n 1.<\/p>\n Generating all \u03d5\u03c4,Q,1(2) fields using normal distribution \\({{{\\mathcal{N}}}}(0,1)\\) and the current value of r parameter (initially we set r and \u0394r to 0.1).<\/p>\n 2.<\/p>\n Updating all \u03d5\u03c4,Q,1(2) of \u03c4\u2009=\u2009\u03b2 by \u0394\u03d5\u03c4,Q,1(2)\u2009=\u2009rx and accept-rejecting the update according to the Metropolis-Hastings scheme.<\/p>\n 3.<\/p>\n Repeating the updates for all time slices from \u03b2 to \u0394\u03c4, then reversely to complete a full sweep, in the process of which Ra is calculated. After each sweep altering r using the rule r\u2009\u2192\u2009r\u2009+\u2009(Ra\u2009\u2212\u20090.4)\u0394r.<\/p>\n 4.<\/p>\n Repeating the sweep until sufficient statistics is gained for all observables.<\/p>\n Hamiltonian update scheme<\/p>\n Since the auxiliary fields are continuous, it is easy to adopt more sophisticated updating schemes like Hamiltonian dynamics to improve the proposed \u0394\u03d5\u03c4,Q,1(2). This approach can significantly boost the acceptance rate, even for large global updates\u2014a key feature of the so-called Hybrid Monte Carlo (HMC) algorithm, which has been successfully applied in condensed matter systems, especially those with long-range Coulomb interactions 132<\/a>.<\/p>\n The partition function in Eq.\u2009(32<\/a>) can be further expressed as<\/p>\n $$Z=\t \\int\\left(\\prod\\limits_{i}d{\\phi }_{i}d{p}_{i}\\right){{{{\\rm{e}}}}}^{-\\frac{1}{2}{\\sum }_{i}\\left({\\phi }_{i}^{2}+{p}_{i}^{2}\\right)}\\det \\left({{{\\rm{I}}}}+{B}_{C}(\\beta,0)\\right)\\\\ \t+{{{\\mathcal{O}}}}(8\\times 1{0}^{-3}),$$<\/p>\n \n (33)\n <\/p>\n where the index i\u2009\u2261\u2009(\u03c4,\u00a0Q,\u00a0m) with m\u2009=\u20091,\u00a02 is employed for compactness and an artificial momentum pi is introduced for \u03d5i. Moreover,<\/p>\n $$\t \\det \\left({{{\\rm{I}}}}+{B}_{C}(\\beta,0)\\right)=\\\\ \t \\quad {{{\\rm{Tr}}}}\\left(\\prod\\limits_{\\tau }{{{{\\rm{e}}}}}^{{{{\\rm{i}}}}{\\sum }_{{{{\\bf{Q}}}}}\\left(-{\\phi }_{\\tau,{{{\\bf{Q}}}},1}\\sqrt{{\\alpha }_{2}({{{\\bf{Q}}}})}{A}_{{{{\\bf{Q}}}}}+{{{\\rm{i}}}}{\\phi }_{\\tau,{{{\\bf{Q}}}},2}\\sqrt{{\\alpha }_{2}({{{\\bf{Q}}}})}{B}_{{{{\\bf{Q}}}}}\\right)}{{{{\\rm{e}}}}}^{-\\bigtriangleup \\tau {H}_{0}}\\right),$$<\/p>\n \n (34)\n <\/p>\n which is nonnegative for charge-neutral TBG. Furthermore, the partition function can be expressed as<\/p>\n $$Z=\\int\\left(\\prod\\limits_{i}d{\\phi }_{i}d{p}_{i}\\right){{{{\\rm{e}}}}}^{-{{{\\mathcal{H}}}}}+{{{\\mathcal{O}}}}(8\\times 1{0}^{-3}),$$<\/p>\n \n (35)\n <\/p>\n with the effective Hamiltonian<\/p>\n $${{{\\mathcal{H}}}}\\equiv \\frac{1}{2}{\\sum}_{i}\\left({\\phi }_{i}^{2}+{p}_{i}^{2}\\right)-\\ln \\left(\\det \\left({{{\\rm{I}}}}+{B}_{C}(\\beta,0)\\right)\\right).$$<\/p>\n \n (36)\n <\/p>\n The Hamiltonian dynamics is applied as<\/p>\n $$\\frac{d{p}_{i}}{dt}=\t -\\frac{\\partial {{{\\mathcal{H}}}}}{\\partial {\\phi }_{i}}\\\\=\t -{\\phi }_{i}+4{{{\\rm{Tr}}}}\\left({J}_{{{{\\bf{Q}}}}}({{{\\rm{I}}}}-G(\\tau ))\\right)+{{{\\mathcal{O}}}}(1{0}^{-2})\\\\ \\frac{d{\\phi }_{i}}{dt}=\t \\frac{\\partial {{{\\mathcal{H}}}}}{\\partial {p}_{i}}\\\\=\t {p}_{i},$$<\/p>\n \n (37)\n <\/p>\n where G(\u03c4) is the equal-time Green\u2019s function for s\u2009=\u2009\u2191 and \u03b7\u2009=\u2009+\u00a0, and JQ is the matrix form of AQ or BQ for m\u2009=\u20091 or 2. The Hamiltonian dynamics facilitates the increased acceptance rates through the conservation of the Hamiltonian along equal-energy trajectories.<\/p>\n The outperformance of the CFMC algorithm<\/p>\n The comparison of performances of CFMC and discretized-field QMC (DQMC) algorithms is shown in Fig.\u20096<\/a> using the same physical settings. It can be clearly seen that CFMC indeed reduces the computational complexity from O(\u03b2N4) to O(\u03b2N3), as the system size increases.<\/p>\n Fig. 6: Computational complexities of DQMC and CFMC with the same physical settings. We measure the CPU time consumed by the two methods for 100 sweeps in the units of CPU core hours, and show the time in a linear-linear scale (a) and a log-log scale (b). Note that inverse temperature also sclaes with the systems size: \u03b2\u2009~\u2009L. The data can be very well described by the fits of \u03b2N4 and \u03b2N3 with N\u2009=\u2009L\u2009\u00d7\u2009L for DQMC and CFMC, respectively, especially in the limit of large system sizes.<\/p>\n The CPU time per update is not the only metric which defines the performance of the Monte Carlo scheme. It can be that different updates lead to substantially different distributions of the observables (despite the average values being fixed)133<\/a>. Thus, the faster algorithm can potentially suffer from larger statistical fluctuations, which can even cancel the gain in the speed of update. In order to check this, we plot the Monte Carlo time series of an IVC structure factor from DQMC and CFMC simlations. Results for \u0398\u2009=\u20091.08\u00b0, 1.\u00a02\u00b0, and 1.\u00a03\u00b0 are shown in Supplementary Fig.\u20096<\/a>. The histograms for the time series are drawn on top of them using the red and blue lines, respectively. The overlap of histograms indicates that fluctuations of DQMC and CFMC are very similar, thus the overall increase in performance from O(\u03b2N4) scaling to O(\u03b2N3) sclaing still holds.<\/p>\n The autocorrelations for the CFMC and DQMC methods are also compared, as illustrated in Supplementary Fig.\u20097<\/a>. Since the measurements are done after the updates of the fields in each time slice, we compare autocorrelation for measurements separated by a Euclidean time interval \u0394t. The autocorrelation time is indeed larger for CFMC, which offsets the factor of N in speedup in simulations at MA and at the smallest lattice size. However, since the autocorrelation time does not depend on the lattice size, moving away from L\u2009=\u20096 and \u0398\u2009=\u20091.08\u00b0 reveals a greater advantage: the difference in scaling (O(\u03b2N4) for DQMC versus O(\u03b2N3) for CFMC) quickly outweighs the impact of the increased autocorrelation time (see the concrete examples in the caption of the Supplementary Fig.\u20097<\/a>. We therefore conclude that CFMC significantly enhances the performances of the Monte Carlo simulations.<\/p>\n Stochastic analytic continuation<\/p>\n The stochastic analytic continuation (SAC)41<\/a>,42<\/a>,43<\/a>,44<\/a>,85<\/a>,86<\/a>,134<\/a> method is also employed to extract the real frequency spectral functions. The Green\u2019s function relates to the real frequency spectral function as 42<\/a><\/p>\n $$G({{{\\bf{k}}}},\\tau )=\\int_{-\\infty }^{\\infty }d\\omega \\left(\\frac{{{{{\\rm{e}}}}}^{-\\omega \\tau }}{1+{{{{\\rm{e}}}}}^{-\\beta \\omega }}\\right)A({{{\\bf{k}}}},\\omega ).$$<\/p>\n \n (38)\n <\/p>\n A variational ansatz of A(k,\u00a0\u03c9) simulating an annealing process is applied. For each momentum k, \\(A({{{\\bf{k}}}},\\omega )={\\sum }_{i=1}^{{N}_{\\omega }}{A}_{i}\\delta (\\omega -{\\omega }_{i})\\) where frequencies \u03c9i initially form regular grid with length of 0.5kBT of N\u03c9\u2009=\u20094000 points symmetrically distributed around 0. Then A(k,\u00a0\u03c9) is estimated by sampling over the grid of \u03c9i, using Eq.\u2009(38<\/a>) and the effective action<\/p>\n $${\\chi }^{2}=\t {\\sum}_{{\\tau }_{1},{\\tau }_{2}}\\left(\\overline{G}({\\tau }_{1})-\\int_{-\\infty }^{\\infty }d\\omega \\left(\\frac{{{{{\\rm{e}}}}}^{-\\omega \\tau }}{1+{{{{\\rm{e}}}}}^{-\\beta \\omega }}\\right)A(\\omega )\\right){({C}^{-1})}_{{\\tau }_{1},{\\tau }_{2}}\\\\ \t \\times \\left(\\overline{G}({\\tau }_{2})-\\int_{-\\infty }^{\\infty }d\\omega \\left(\\frac{{{{{\\rm{e}}}}}^{-\\omega \\tau }}{1+{{{{\\rm{e}}}}}^{-\\beta \\omega }}\\right)A(\\omega )\\right),$$<\/p>\n \n (39)\n <\/p>\n where<\/p>\n $$\\begin{array}{l}{C}_{{\\tau }_{1},{\\tau }_{2}}\\\\=\\frac{1}{{N}_{b}({N}_{b}-1)}{\\sum }_{b=1}^{{N}_{b}}\\left({G}^{b}({\\tau }_{1})-\\overline{G}({\\tau }_{1})\\right)\\left({G}^{b}({\\tau }_{2})-\\overline{G}({\\tau }_{2})\\right),\\end{array}$$<\/p>\n \n (40)\n <\/p>\n with \u03c41 and \u03c42 referring to selected imaginary times where the errors of Green\u2019s functions are less than 0.1. Gb(\u03c41,2) is the Green\u2019s function of b-th bin at \u03c41,2, and \\(\\overline{G}({\\tau }_{1,2})\\) is the average of Gb(\u03c41,2) over all bins. Importance sampling of Metropolis-Hastings type is employed, and the sampling weight is \\(W=\\exp \\left(-{\\chi }^{2}\/(2{T}_{1})\\right)\\) where T1 is the artificial temperature, which decreases from 5 by a factor of 0.6 in each step. The final value of T1 is defined by the average effective action satisfying the relation\\(\\langle {\\chi }^{2}\\rangle={\\chi }_{\\min }^{2}+2\\sqrt{{\\chi }_{\\min }^{2}}\\)42<\/a>, and this final value of T1 is used to obtain the spectra.<\/p>\n Computation of the free energy<\/p>\n The free energy is difficult to obtain numerically, since exponentially large or small numbers are usually involved. Let us denote the partition function as Z\u2009=\u2009\u2211CWCPC, where WC corresponds to the bosonic part of the action and PC corresponds to the fermionic determinant. Since the kinetic energy is included as \\({{{{\\rm{e}}}}}^{-\\bigtriangleup \\tau {H}_{0}}\\) in Eq.\u2009(32<\/a>), PC becomes exponentially large as the kinetic energy grows with twist angle (see Supplementary Fig.\u20092<\/a>). One might think that positive and negative terms in the matrix \\({{{{\\rm{e}}}}}^{-\\bigtriangleup \\tau {H}_{0}}\\) form pairs, so that they cancel each other and the magnitude of PC does not change a lot with the twist angle. However, the identity matrix in the relation, \\({{{\\rm{Tr}}}}\\left({{{{\\rm{e}}}}}^{{\\sum }_{i,j}{c}_{i}^{{{\\dagger}} }{A}_{i,j}{c}_{j}}{{{{\\rm{e}}}}}^{{\\sum }_{i,j}{c}_{i}^{{{\\dagger}} }{B}_{i,j}{c}_{j}}\\right)=\\det \\left({{{\\rm{I}}}}+{{{{\\rm{e}}}}}^{A}{{{{\\rm{e}}}}}^{B}\\right)\\), forbids this cancellation.<\/p>\n The free energy F can be theoretically measured by<\/p>\n $$F \\equiv \t -\\frac{1}{\\beta }\\ln (Z)\\\\=\t -\\frac{1}{\\beta }\\ln \\left(\\sum\\limits_{C}{W}_{C}{P}_{C}\\right)\\\\=\t -\\frac{1}{\\beta }\\ln \\left(\\left(\\sum\\limits_{C}{W}_{C}\\right)\\frac{{\\sum }_{C}{W}_{C}{P}_{C}}{{\\sum }_{C}{W}_{C}}\\right)\\\\=\t -\\frac{1}{\\beta }\\ln \\left(\\sum\\limits_{C}{W}_{C}\\right)-\\frac{1}{\\beta }\\ln \\left(\\frac{{\\sum }_{C}{W}_{C}{P}_{C}}{{\\sum }_{C}{W}_{C}}\\right),$$<\/p>\n \n (41)\n <\/p>\n which simply means the importance sampling of PC over the weight defined by WC. However, PC is exponentially large in our case, so that direct sampling of a moderate-size system is inaccessible.<\/p>\n Recently, an integral algorithm to efficiently sample exponential observables such as free energy and entanglement entropy was developed 87<\/a>. Dumping the constant number \\(-\\ln \\left({\\sum }_{C}{W}_{C}\\right)\/\\beta\\) for a specific T, we can rewrite F as<\/p>\n $$F=-\\frac{1}{\\beta }\\int_{0}^{1}dt\\frac{\\partial \\ln (f(t))}{\\partial t},$$<\/p>\n \n (42)\n <\/p>\n if \\(f(1)=\\left({\\sum }_{C}{W}_{C}{P}_{C}\\right)\/\\left({\\sum }_{C}{W}_{C}\\right)\\) and f(0)\u00a0=\u00a01. Subsequently, we define f(t) as f(t)\u2009=\u2009\u3008etX\u3009, so that \\(\\langle {{{{\\rm{e}}}}}^{X}\\rangle=\\left({\\sum }_{C}{W}_{C}{P}_{C}\\right)\/\\left({\\sum }_{C}{W}_{C}\\right)\\) and \\(X=\\ln (P)\\). Substituting it in (42<\/a>), we get<\/p>\n $$F=\t -\\frac{1}{\\beta }\\int_{0}^{1}dt\\frac{\\langle X{{{{\\rm{e}}}}}^{tX}\\rangle }{\\langle {{{{\\rm{e}}}}}^{tX}\\rangle }\\\\=\t -\\frac{1}{\\beta }\\int_{0}^{1}dt\\frac{{\\sum }_{C}{W}_{C}{P}_{C}^{t}\\ln ({P}_{C})}{{\\sum }_{C}{W}_{C}{P}_{C}^{t}}.$$<\/p>\n \n (43)\n <\/p>\n Thus, the exponentially large values of observables are avoided and we need only the importance samplings of \\(\\ln ({P}_{C})\\) over the weight of \\({W}_{C}{P}_{C}^{t}\\).<\/p>\n Pad\u00e9 approximants and scaling exponents<\/p>\n The critical exponents characterize a particular universality class and depend on the dimensionality, symmetry and the relevant degrees of freedom. The transition to the KIVC order in TBG falls in the universality classe of the chiral XY GN model. This was considered previously via perturbative renormalization up to four-loop order in an \u03f5 expansion54<\/a>, where \u03f5 is the deviation from the the upper critical spacetime dimension D\u2009=\u20094\u2009\u2212\u2009\u03f5 of the model In ref. 54<\/a>, a scale dependence is introduced in the fields \u03d5 and \u03c8, as well as for the mass term m, quartic coupling \u03bb, and Yukawa coupling g from Eq.\u2009(3<\/a>), and their renormalization is calculated. Based on this, the inverse correlation-length exponent \u03bd\u22121, the correction-to-scaling exponent \u03c9, the fermion and the boson anomlous dimension \u03b7\u03c8,\u00a0\u03b7\u03d5 are determined and their expression for general fermion flavor number Nf is provided up to fourth order in \u03f5. We use these results and apply them to our case of interest, the KIVC-DSM transition in TBG. To obtain an estimate for the critical exponents for the case of (2+1) dimensions, we calculate the different Pad\u00e9 approximant<\/p>\n $${P}_{[m\/n]}(\\epsilon )=\\frac{{a}_{0}+{a}_{1}\\epsilon+\\cdots+{a}_{m}{\\epsilon }^{m}}{1+{b}_{1}\\epsilon+\\cdots+{b}_{n}{\\epsilon }^{n}}$$<\/p>\n \n (44)\n <\/p>\n for a general fermion flavor number Nf for the exponents \u03bd\u22121,\u00a0\u03b7\u03d5. Selecting only those that do not exhibit any singularities or poles, which in our case are P[2\/2](\u03f5) and P[3\/1](\u03f5), we obtain the numerical values for Nf\u2009=\u200916 and \u03f5\u2009=\u20091 and use the hyperscaling relation to determine the order parameter exponent \u03b2, \u03b2\u2009=\u2009\u03bd(D\u2009\u2212\u20092\u2009+\u2009\u03b7\u03d5)\/2. We subsequently average over these values and based on those provide the numerical uncertainties.<\/p>\n","protected":false},"excerpt":{"rendered":"Continuum model of TBG The schematic twist of graphene Brillouin zone and the consequent moir\u00e9 Brillouin zone are…\n","protected":false},"author":2,"featured_media":61878,"comment_status":"","ping_status":"","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[24],"tags":[49,48,3881,1099,1100,3883,314,66],"class_list":{"0":"post-61877","1":"post","2":"type-post","3":"status-publish","4":"format-standard","5":"has-post-thumbnail","7":"category-physics","8":"tag-ca","9":"tag-canada","10":"tag-electronic-properties-and-materials","11":"tag-humanities-and-social-sciences","12":"tag-multidisciplinary","13":"tag-phase-transitions-and-critical-phenomena","14":"tag-physics","15":"tag-science"},"_links":{"self":[{"href":"https:\/\/www.newsbeep.com\/ca\/wp-json\/wp\/v2\/posts\/61877","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.newsbeep.com\/ca\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.newsbeep.com\/ca\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.newsbeep.com\/ca\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/www.newsbeep.com\/ca\/wp-json\/wp\/v2\/comments?post=61877"}],"version-history":[{"count":0,"href":"https:\/\/www.newsbeep.com\/ca\/wp-json\/wp\/v2\/posts\/61877\/revisions"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.newsbeep.com\/ca\/wp-json\/wp\/v2\/media\/61878"}],"wp:attachment":[{"href":"https:\/\/www.newsbeep.com\/ca\/wp-json\/wp\/v2\/media?parent=61877"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.newsbeep.com\/ca\/wp-json\/wp\/v2\/categories?post=61877"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.newsbeep.com\/ca\/wp-json\/wp\/v2\/tags?post=61877"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}![\"figure](\"https:\/\/www.newsbeep.com\/ca\/wp-content\/uploads\/2025\/08\/41467_2025_62461_Fig5_HTML.png\")
<\/a><\/p>\n![\"figure](\"https:\/\/www.newsbeep.com\/ca\/wp-content\/uploads\/2025\/08\/41467_2025_62461_Fig6_HTML.png\")
<\/a><\/p>\n