††thanks: These authors contributed equally to the work.††thanks: These
authors contributed equally to the work.

# Fermi surface symmetric mass generation: a quantum Monte-Carlo study

Wei-Xuan Chang  Beijing National Laboratory for Condensed Matter Physics and
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China  Sibo
Guo  Beijing National Laboratory for Condensed Matter Physics and Institute of
Physics, Chinese Academy of Sciences, Beijing 100190, China  Yi-Zhuang You
<EMAIL_ADDRESS>Department of Physics, University of California, San Diego, CA
92093, USA  Zi-Xiang Li<EMAIL_ADDRESS>Beijing National Laboratory for
Condensed Matter Physics and Institute of Physics, Chinese Academy of
Sciences, Beijing 100190, China  University of Chinese Academy of Sciences,
Beijing 100049, China

###### Abstract

The symmetric mass generation (SMG) phase is an insulator in which a single-
particle gap is intrinsically opened by the interaction, without involving
symmetry spontaneously breaking or topological order. Here, we perform
unbiased quantum Monte-Carlo simulation and systematically investigate a
bilayer fermionic model hosting Fermi surface SMG in the strongly interacting
regime. With increasing interaction strength, the model undergoes a quantum
phase transition from an exciton insulator to an SMG phase, belonging to the
(2+1)-dimensional O(4) universality class. We access the spectral properties
of the SMG phase, resembling a Mott insulating phase with relatively flat
dispersion and pronounced spectral broadening. The dispersion of Green’s
function zeros is extracted from spectral function, featuring a surface at
zero frequency precisely located at the original non-interacting Fermi
surface, which constitutes a hallmark of the Fermi surface SMG phase. The
bilayer model we study is potentially relevant to the newly discovered
high-$T_{c}$ superconductor $\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$. Our results in
SMG phase qualitatively capture the salient features of spectral function
unveiled in recent ARPES experiments, shedding new insight on the underlying
physics of $\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$.

Introduction: Symmetric mass generation (SMG) is a mechanism by which a metal
acquires excitation gaps without breaking symmetry or developing topological
order[1]. In recent years, SMG has attracted increasing attention in various
areas, including particle physics[2, 3, 4, 5, 6, 7, 8, 9] and condensed matter
physics. In condensed matter physics, SMG is investigated in the framework of
interacting topological insulator or superconductor[10, 11, 12, 13, 14, 15,
16, 17, 18], particularly the interaction reduced classification of symmetry
protected topological phases. Additionally, it is demonstrated that SMG
transition point in the Dirac fermions constitutes a fermionic version of
deconfined quantum critical point (QCP)[19, 20, 21, 22], enormously enriching
the understanding of exotic QCP beyond Landau’s paradigm[23, 24, 25, 26, 27,
28].

Figure 1:  (a) The ground-state phase diagram of Eq. (1). The system undergoes
a quantum phase transition from exciton insulator to SMG phase at $J=3.24$.
For the exciton insulator, the yellow (green) bonds denote inter-layer exciton
order with momentum $(\pi,\pi)$. For the SMG phase, the inter-layer electrons
form a spin singlet, preserving all the lattice and spin rotational
symmetries. (b) The results of correlation-length ratio for exciton order as a
function of $J$ for different system sizes.

More recently, Fermi surface SMG[29, 30, 31], denoting the gap opening without
SSB or topological order in the Fermi surface, emerges as a fascinating
subject due to its fundamental interests relevant to Fermi surface anomaly[32,
33] and Luttinger’s theorem[34]. The bilayer model featuring strong anti-
ferromagnetic (AFM) inter-layer spin interaction is proposed as a simple
microscopic model to host the Fermi surface SMG phase. Intriguingly, from
field theory analysis and perturbative numerical calculation, it is suggested
that the poles of Green’s function in the pristine gapless fermion state are
replaced by the Green’s function zeros, which encloses a Fermi volume
fulfilling the requirement of Luttinger’s theorem[30]. The surface of Green’s
function zeros at low energy is the hallmark of Fermi surface SMG.
Nevertheless, the unbiased demonstration is still lacking, and it is immensely
desirable to investigate Fermi surface SMG by unbiased theoretical approach,
unraveling the property of Green’s function zeros and other prominent features
relevant to experimental observables.

From the perspective of realistic materials, the bilayer model hosting SMG
insulator is potentially relevant to the newly discovered high-$T_{c}$
superconductor $\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$, which exhibits
superconductivity with $T_{c}$ up to 80K under pressure[35, 36, 37, 38].
$\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$ superconductor harbors a layered structure,
with each unit cell containing two layers of $\rm{Ni}\rm{O}_{2}$. The DFT
calculation[39, 40, 41, 42, 43, 44] and ARPES measurements [45] reveal that
the low-energy band structure of $\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$ under
pressure is dominated by the electrons in Ni-$3d_{z^{2}}$ and
Ni-$3d_{x^{2}-y^{2}}$ orbitals. With increasing pressure, the material
undergoes a structure transition, rendering a strong enhancement of inter-
layer AFM Heisenberg interaction between Ni-$3d_{z^{2}}$ electrons. Distinct
from Cuprates superconductors[46, 47, 48], although more experimental efforts
are required, the magnetic orders have not been detected so far in
$\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$. Hence, $\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$
is a prospective platform for experimental realizations of Fermi surface SMG.
Moreover, the theoretical understanding of the underlying mechanism of
high-$T_{c}$ SC remains under development[49, 50, 51, 52, 53, 54, 55, 56, 57,
58, 59, 60, 61, 62, 63, 64, 65, 66]. Fathoming the accurate properties of the
bilayer model with strong inter-layer spin interaction and the emergent SMG
phase by an intrinsically unbiased theoretical approach offers a promising
route to demystify the microscopic origin of high-$T_{c}$ SC in
$\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$.

To address the aforementioned questions, we perform numerically-exact quantum
Monte-Carlo (QMC) simulation[67, 68, 69, 70, 71] to systematically investigate
a bilayer fermionic model with inter-layer AFM Heisenberg interaction on
square lattice. Remarkably, the model is free from the notorious sign
problem[72, 73, 74, 75, 76, 77, 78, 79, 71, 70], so the ground-state
properties of the model with large system size are accessible. We achieve the
ground-state phase diagram of the model. In the weak coupling regime, the
ground state possesses inter-layer exciton order, degenerated with inter-layer
spin-singlet pairing order. With increasing interaction strength, the exciton
order is suppressed and the system undergoes a quantum phase transition to a
SMG phase without any SSB. The large-scale QMC simulation reveals the
transition belongs to the $(2+1)$-dimensional O(4) universality class.
Furthermore, we investigate the spectral properties of the SMG phase. The
dispersion of the spectral peak resembles of Mott insulator, featuring a
relatively flat band and pronounced broadening. The surface of Green’s
function zeros at zero frequency is extracted from the results of the spectral
function, which resides on the original Fermi surface in the non-interacting
limit, fulfilling the requirement of Luttinger’s theorem. Lastly, we discuss
the implication of our numerical results on the
$\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$ high-$T_{c}$ superconductors and potential
detection of Fermi surface SMG phase in future experiments.

Model: We consider the bilayer model on a square lattice with inter-layer AFM
Heisenberg interaction, as described by the following Hamiltonian[80, 81, 82,
83, 84, 85, 86, 87, 30]:

$\displaystyle H=-t\sum_{\left\langle
ij\right\rangle{\sigma}\alpha}(c^{\dagger}_{i{\sigma}\alpha}c_{j{\sigma}\alpha}+\mathrm{h.c.})+J\sum_{i}\vec{S}_{i,1}\cdot\vec{S}_{i,2},~{}~{}~{}$
(1)

where $c^{\dagger}_{i{\sigma}}$ creates an electron on site $i$ with spin
polarization ${\sigma}=\uparrow,\downarrow$ and layer index $\alpha=1,2$. Here
$t$ is the hopping amplitude of electrons between nearest-neighbor sites and
$J$ denotes the strength of inter-layer AFM Heisenberg interaction. Hereafter,
we set $t=1$ as the energy unit and focus on the model at half-filling.
Remarkably, the model is free from the notorious sign problem such that we can
access the ground-state properties of the model with large system sizes
through the projective version of QMC simulation[88]. The details of the
projector QMC algorithm are included in the Supplementary Materials.

The model preserves total spin SU(2) symmetry. Additionally, the model at
half-filling respects pseudospin SU(2) symmetry for each layer[30]. In the
non-interacting limit, namely $J=0$, the model at half-filling features a
Fermi surface with nesting momentum $(\pi,\pi)$. In the weak coupling regime,
the instability with spontaneous symmetry breaking appears due to the Fermi
surface nesting. Conversely, in the strong coupling limit
$J\rightarrow\infty$, the ground state wave function is
$\left|\psi\right\rangle=\prod_{i}(c^{\dagger}_{i1\uparrow}c^{\dagger}_{i2\downarrow}-c^{\dagger}_{i1\downarrow}c^{\dagger}_{i2\uparrow})\left|0\right\rangle$.
The state is the product of the inter-layer spin singlet state on each site,
respecting all symmetries, including the lattice symmetry, the spin rotational
symmetry, and the charge conservation symmetries for each layer. In the
following, we perform a sign-problem-free QMC simulation to study the accurate
ground-state phase diagram of the model for generic coupling $J$.

Figure 2: (a) The data collapse analysis of correlation-length ratio for
excitonic order. The resulting correlation-length exponent $\nu=0.735\pm
0.024$. The rescaled correlation-length ratios for different system sizes are
collapsed into a single curve. (b) Logarithm of correlation function for
exciton order versus $\ln(r)$ at QCP $J=3.24$, giving rise to the anomalous
dimension $\eta=0.051\pm 0.036$ from the scaling function
$C(r,r)\sim\frac{1}{r^{1+\eta}}$.  Figure 3: The results of spectral function
$A(k,\omega)$ along the high-symmetry momentum line for (a) $J=2.4$, (b)
$J=2.8$, (c) $J=3.25$, (d) $J=4.0$, (e) $J=4.8$, (f) $J=6.0$.

Ground-state phase diagram: We investigate the ground-state properties of the
model Eq. (1) with varying interaction strength $J$. To unveil the possible
symmetry-breaking phase, we calculate the structure factors of distinct
orderings and the corresponding correlation-length ratios, the detailed
descriptions of which are included in the Supplementary Materials. The results
of structure factors and correlation-length ratios, as included in SM,
unambiguously reveal that the inter-layer exciton condensation (EC) is the
dominant one in the weak coupling regime. The order parameter of the EC order
is given by $\phi_{\rm EC}=\frac{1}{N}\sum_{i\sigma}c^{\dagger}_{i\sigma
1}c_{i\sigma 2}(-)^{i}$ with $(-)^{i}=\pm 1$ for $i\in\text{A,B}$ sublattices,
which carries momentum $(\pi,\pi)$ due to Fermi surface nesting. As verified
by the scaling analysis of structure factor and correlation-length ratios, the
ground state possesses EC long-range order when $J$ is small.

Further increasing interaction strength $J$, the exciton order is suppressed
and vanishes at a critical $J$. To identify the quantum critical point
$J_{c}$, we investigate the behaviors of correlation length ratios for
different system sizes. The results are depicted in Fig. 1(b), unequivocally
revealing a quantum phase transition between EC long-range ordered and
disordered phases occurring at $J_{c}=3.24$. In the parameter regime $J>3.24$,
we confirm the absence of SSB in the ground state by evaluating the
correlation-length ratio for different orderings, as shown in SM. Moreover,
the single-particle gaps are finite in the whole regime of $J$ (as shown in
SM), substantiating that the model Eq. (1) hosts Fermi surface SMG insulating
phase when $J>3.24$. Hence, we arrive at the ground-state phase diagram as
shown in Fig. 1(a).

Quantum criticality: As aforementioned, in addition to the total spin SU(2)
symmetry, the model at half filling also respects the
$\rm{SU}(2)_{1}\times\rm{SU}(2)_{2}$ pesudospin symmetry, where
$\rm{SU}(2)_{l}$ (for each layer $l=1,2$) is generated by
$\bm{K}_{l}=\sum_{i}\frac{1}{2}(-)^{i}\tilde{c}_{il}^{\dagger}\bm{\sigma}\tilde{c}_{il}$
with $\tilde{c}_{il}=(c_{il\uparrow},c_{il\downarrow}^{\dagger})^{\intercal}$
being the Nambu spinor. Owing to the enlarged pseudospin symmetry, EC order is
degenerate with inter-layer singlet SC order. The degeneracy is explicitly
explained by the fact that order parameter of EC $\phi_{\rm
EC}=\frac{1}{N}\sum_{i\sigma}c^{\dagger}_{i\sigma 1}c_{i\sigma 2}(-1)^{i}$ is
related to the inter-layer SC order $\phi_{\rm
SC}=\frac{1}{N}\sum_{i\sigma}(-)^{\sigma}c^{\dagger}_{i\sigma
1}c^{\dagger}_{i\bar{\sigma}2}$ by the particle-hole transformation
$c_{i\sigma 2}\rightarrow(-)^{i}(-)^{\sigma}c^{\dagger}_{i\bar{\sigma}2}$ (a
subgroup of $\rm{SU}(2)_{1}\times\rm{SU}(2)_{2}$), where $(-)^{\sigma}=\pm 1$
for $\sigma=\uparrow,\downarrow$ and $\bar{\sigma}$ represents the opposite
spin of $\sigma$. EC and inter-layer SC order preserve total spin rotational
symmetry, but break pesudospin symmetries
$\rm{SU}(2)_{1}\times\rm{SU}(2)_{2}$, which is equivalent to $\rm{SO}(4)$
symmetry. Consequently, the quantum phase transition from exciton insulator/SC
to SMG phase is presumably an O(4) Wilson-Fisher transition.

To decipher the critical properties of the quantum phase transition, we
perform a finite-size scaling analysis to extract the critical exponents. The
RG-invariant correlation-length ratios for different system sizes satisfy the
following scaling form for sufficiently large $L$ and $J$ close to $J_{c}$:

$\displaystyle R_{\rm{EC}}=\mathscr{F}(L^{1/\nu}(J-J_{c})),~{}~{}~{}$ (2)

where $\nu$ is correlation-length exponent and $\mathscr{F}$ is a non-
universal scaling function. Upon choosing the optimal value of $\nu$, the data
points of relatively large $L$ should collapse on a single curve. The results
of data collapse analysis are shown in Fig. 2(a), giving rise to the
correlation-length exponent $\nu=0.735\pm 0.024$. Fixing $J$ at QCP
$J_{c}=3.24$, we evaluate the correlation function of exciton order
$C(r)=\left\langle\phi_{\rm EC}(0)\phi_{\rm EC}(r)\right\rangle$. The
anomalous dimension of exciton order is determined by the scaling form
$C(r)\sim r^{-1-\eta}$, the results of which are shown in Fig. 2(b), yielding
the anomalous dimension $\eta=0.051\pm 0.036$. The obtained critical exponents
$\nu=0.735\pm 0.024$ and $\eta=0.051\pm 0.036$ are reasonably consistent with
previous results of 3D O(4) transition[89], in which $\nu=0.7479(90)$ and
$\eta=0.0258(22)$. Hence, our large-scale QMC simulation provides convincing
evidence that the quantum phase transition between the exciton insulator and
SMG phase belongs to (2+1)D O(4) Wilson-Fisher universality class.

Figure 4:  The real part of Green’s function extracted from the spectral
function by KK relation for $J=4.0$ (a) and $J=6.0$ (b). The dashed lines
represent the zeros of Green’s function.

Spectral properties and Green’s function zeros: One prominent feature of the
SMG insulator is the zeros of single-particle Green’s function at low energy.
The celebrated Luttinger’s theorem dictates the volume enclosed by the poles
of Green’s function in momentum space yields the density of electrons. In the
SMG state, the surface of Green’s function poles is replaced by the zeros of
Green’s function, which defines Luttinger volume to match the particle
density. To delve into such a prominent feature of the SMG phase hosted in Eq.
(1), we explore the spectral properties and extract the Green’s function zeros
in the model.

Employing QMC simulation, we evaluate the imaginary-time single-particle
Green’s function and access the spectral function through stochastic
analytical continuation (SAC)[90, 91]. The results of spectral function for
different values of $J$ are presented in Fig. 3. In the exciton insulator
phase with weak interaction, for instance, $J=2.4$, the spectral function
reveals the dispersion of a symmetry-breaking band insulator, displaying the
minimum band gap located near the Fermi surface in the non-interacting limit.
The dispersion exhibits the band-folding with momentum $(\pi,\pi)$, consistent
with the ordering momentum of exciton order — a clear indication of
translational symmetry breaking. With increasing interaction strength $J$,
exciton order is suppressed and the system evolves into a SMG phase. Deep in
the SMG phase, the peaks of spectral function form two bands with relatively
flat dispersion, resembling the spectrum of Mott insulator. Notably, with
increasing $J$, the spectral peak in SMG phase exhibits remarkably pronounced
broadening, manifesting an insulator triggered by strong interaction without
well-defined quasi-particle.

Furthermore, we extract the real part of Green’s function from the single-
particle spectral function through the Kramers-Kronig (KK) relation:

$\displaystyle\text{Re}\,G(k,\omega)=\frac{1}{2\pi}\int
d\omega^{\prime}\frac{A(k,\omega^{\prime})}{\omega-\omega^{\prime}},$ (3)

and investigate the structures of Green’s function zeros, defined by
$G(k,\omega)=0$. The results for $J=4$ and $J=6$ are presented in Fig. 4(a)
and (b), respectively, in which the zeros of Green’s function real part are
dictated by the dashed line. The dispersion of Green’s function’s zeros
appears and is reversed compared with the original non-interacting band
dispersion. The “bandwidth” of the Green’s function zeros is suppressed with
increasing interaction strength $J$ [30]. Crucially, at $\omega=0$ the zeros
of Green are precisely located in the Fermi surface of the model in the non-
interacting limit. As a consequence, the volume enclosed by the surface of
zeros is equivalent to the corresponding pole Fermi surface in the non-
interacting limit, fulfilling the requirement of Luttinger theory.

Relevance to $La_{3}Ni_{2}O_{7}$ superconductor: As aforementioned, the model
in Eq. (1) is potentially relevant to the high-$T_{c}$ superconductor
$\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$, capturing the essential physics of
Ni-$3d_{z^{2}}$ electrons in $\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$. As revealed
by the recent ARPES measurements, the $3d_{z^{2}}$ electrons exhibit nearly
flat dispersion and pronounced strong mass renormalization, in sharp contrast
with the band dispersion of $3d_{x^{2}-y^{2}}$ electrons[45]. The nearly flat
dispersion and large mass renormalization resemble the spectral function in
the SMG insulating phase as shown in Fig. 3, implying that the strong inter-
layer spin interaction offers a possible origin of the electronic structure
observed in $\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$[45]. The Fermi surface SMG
phase is a plausible parent insulating state from which SC arises, similar to
the AFM Mott insulating phase in cuprates. More theoretical and experimental
investigations in this direction are needed to corroborate this scenario.

Moreover, our unbiased QMC simulation suggests that it is feasible to extract
the surface of Green’s function zeros from the spectral function, albeit in
the presence of the pronounced broadening of spectral peak and noise arising
from analytical continuation. The numerical results shed new light on the
observation of Fermi surface zeros in experiments. It is promising to
implement the procedure in our study to extract the information of Green’s
function zeros from the spectral function measured in ARPES experiments and
detect the experimental evidence of the Fermi surface SMG insulating phase.

Concluding remarks: In this paper, we perform large-scale QMC simulation to
investigate the SMG phase in an interacting bilayer model on a square lattice.
The model features exciton/SC order breaking $\rm{SO}(4)$ pseudospin symmetry
in the weak coupling regime and transitions to the SMG insulator in the
presence of strong Heisenberg interaction. We determine the ground-state phase
diagram and measure the critical properties of the quantum phase transition
between the exciton insulator and SMG phase, which confirms an $\rm{O}(4)$
Wilson-Fisher universality class. We systematically study the spectral
properties in the SMG insulator and exciton insulator. In the SMG insulator,
the spectral peaks feature a relatively flat dispersion with pronounced
broadening, resembling a Mott insulator without well-defined quasi-particles.
More essentially, we determine the zeros of Green’s function from the spectral
function. The Green’s function zeros in the zero frequency enclose the Fermi
surface of the non-interacting band structure, replacing the roles of Green’s
function poles in the non-interacting limit. The Fermi volume enclosed by the
surface associated with Green’s function zeros is equivalent to Luttinger
volume in the non-interacting limit, thus retaining Luttinger’s theorem in the
SMG insulating phase. The spectral properties revealed in the numerical
simulation capture the prominent features of ARPES measurements of the
Ni-$3d_{z^{2}}$ electrons in $\rm{La}_{3}\rm{Ni}_{2}\rm{O}_{7}$. Remarkably,
the model Eq. (1) is still sign-problem-free away from half-filling, hence it
is interesting to investigate the possible SC emerging from doping the SMG
phase by unbiased QMC simulation, which is left for our future work.

Acknowledgement: We are indebted to Kun Jiang, Yuan-Yao He and Yingfei Gu for
stimulating discussions. Z.X.L acknowledges support from the start-up grant of
IOP-CAS. Y.Z.Y. is supported by the National Science Foundation Grant No.
DMR-2238360.

## References

  * Wang and You [2022] J. Wang and Y.-Z. You, Symmetry 14, 10.3390/sym14071475 (2022).
  * BenTov [2015] Y. BenTov, Journal of High Energy Physics 2015, 34 (2015).
  * Wen [2013] X.-G. Wen, Chinese Physics Letters 30, 111101 (2013).
  * You and Xu [2015] Y.-Z. You and C. Xu, Phys. Rev. B 91, 125147 (2015).
  * Zeng  _et al._ [2022] M. Zeng, Z. Zhu, J. Wang, and Y.-Z. You, Phys. Rev. Lett. 128, 185301 (2022).
  * Ayyar and Chandrasekharan [2016] V. Ayyar and S. Chandrasekharan, Phys. Rev. D 93, 081701 (2016).
  * Ayyar and Chandrasekharan [2017] V. Ayyar and S. Chandrasekharan, Phys. Rev. D 96, 114506 (2017).
  * Wang and Wen [2020] J. Wang and X.-G. Wen, Phys. Rev. Res. 2, 023356 (2020).
  * Tong [2022] D. Tong, Journal of High Energy Physics 2022, 1 (2022).
  * Fidkowski and Kitaev [2010] L. Fidkowski and A. Kitaev, Phys. Rev. B 81, 134509 (2010).
  * Ryu and Zhang [2012] S. Ryu and S.-C. Zhang, Phys. Rev. B 85, 245132 (2012).
  * Qi [2013] X.-L. Qi, New Journal of Physics 15, 065002 (2013).
  * Yao and Ryu [2013] H. Yao and S. Ryu, Phys. Rev. B 88, 064507 (2013).
  * Wang  _et al._ [2014] C. Wang, A. C. Potter, and T. Senthil, Science 343, 629 (2014), https://www.science.org/doi/pdf/10.1126/science.1243326 .
  * Wang and Senthil [2014] C. Wang and T. Senthil, Phys. Rev. B 89, 195124 (2014).
  * Slagle  _et al._ [2015] K. Slagle, Y.-Z. You, and C. Xu, Phys. Rev. B 91, 115121 (2015).
  * Cheng  _et al._ [2018] M. Cheng, Z. Bi, Y.-Z. You, and Z.-C. Gu, Phys. Rev. B 97, 205109 (2018).
  * Xu and Xu [2021] Y. Xu and C. Xu, Green’s function zero and symmetric mass generation (2021), arXiv:2103.15865 [cond-mat.str-el] .
  * You  _et al._ [2018a] Y.-Z. You, Y.-C. He, C. Xu, and A. Vishwanath, Phys. Rev. X 8, 011026 (2018a).
  * You  _et al._ [2018b] Y.-Z. You, Y.-C. He, A. Vishwanath, and C. Xu, Phys. Rev. B 97, 125112 (2018b).
  * Hou and You [2023] W. Hou and Y.-Z. You, Phys. Rev. B 108, 125130 (2023).
  * He  _et al._ [2016] Y.-Y. He, H.-Q. Wu, Y.-Z. You, C. Xu, Z. Y. Meng, and Z.-Y. Lu, Phys. Rev. B 94, 241111 (2016).
  * Senthil  _et al._ [2004a] T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P. A. Fisher, Science 303, 1490 (2004a), https://www.science.org/doi/pdf/10.1126/science.1091806 .
  * Senthil  _et al._ [2004b] T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and M. P. A. Fisher, Phys. Rev. B 70, 144407 (2004b).
  * Levin and Senthil [2004] M. Levin and T. Senthil, Phys. Rev. B 70, 220403 (2004).
  * Sandvik [2007] A. W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007).
  * Shao  _et al._ [2016] H. Shao, W. Guo, and A. W. Sandvik, Science 352, 213 (2016).
  * Li  _et al._ [2017] Z.-X. Li, Y.-F. Jiang, S.-K. Jian, and H. Yao, Nature Communications 8, 314 (2017).
  * Lu  _et al._ [2023a] D.-C. Lu, M. Zeng, J. Wang, and Y.-Z. You, Phys. Rev. B 107, 195133 (2023a).
  * Lu  _et al._ [2023b] D.-C. Lu, M. Zeng, and Y.-Z. You, Green’s function zeros in fermi surface symmetric mass generation (2023b), arXiv:2307.12223 [cond-mat.str-el] .
  * Lu  _et al._ [2023c] D.-C. Lu, J. Wang, and Y.-Z. You, Definition and classification of fermi surface anomalies (2023c), arXiv:2302.12731 [cond-mat.str-el] .
  * Else  _et al._ [2021] D. V. Else, R. Thorngren, and T. Senthil, Phys. Rev. X 11, 021005 (2021).
  * Else and Senthil [2021] D. V. Else and T. Senthil, Phys. Rev. Lett. 127, 086601 (2021).
  * Luttinger [1960] J. M. Luttinger, Phys. Rev. 119, 1153 (1960).
  * Sun  _et al._ [2023] H. Sun, M. Huo, X. Hu, J. Li, Z. Liu, Y. Han, L. Tang, Z. Mao, P. Yang, B. Wang, J. Cheng, D.-X. Yao, G.-M. Zhang, and M. Wang, Nature 10.1038/s41586-023-06408-7 (2023).
  * Hou  _et al._ [2023] J. Hou, P. T. Yang, Z. Y. Liu, J. Y. Li, P. F. Shan, L. Ma, G. Wang, N. N. Wang, H. Z. Guo, J. P. Sun, Y. Uwatoko, M. Wang, G. M. Zhang, B. S. Wang, and J. G. Cheng, Emergence of high-temperature superconducting phase in the pressurized $\mathrm{La_{3}Ni_{2}O_{7}}$ crystals (2023), arXiv:2307.09865 [cond-mat.supr-con] .
  * Wang  _et al._ [2023] G. Wang, N. Wang, J. Hou, L. Ma, L. Shi, Z. Ren, Y. Gu, X. Shen, H. Ma, P. Yang, Z. Liu, H. Guo, J. Sun, G. Zhang, J. Yan, B. Wang, Y. Uwatoko, and J. Cheng, Pressure-induced superconductivity in polycrystalline $\mathrm{La_{3}Ni_{2}O_{7}}$ (2023), arXiv:2309.17378 [cond-mat.supr-con] .
  * Zhang  _et al._ [2023a] Y. Zhang, D. Su, Y. Huang, H. Sun, M. Huo, Z. Shan, K. Ye, Z. Yang, R. Li, M. Smidman, M. Wang, L. Jiao, and H. Yuan, High-temperature superconductivity with zero-resistance and strange metal behavior in la3ni2o7 (2023a), arXiv:2307.14819 [cond-mat.supr-con] .
  * Luo  _et al._ [2023a] Z. Luo, X. Hu, M. Wang, W. Wú, and D.-X. Yao, Phys. Rev. Lett. 131, 126001 (2023a).
  * Geisler  _et al._ [2023] B. Geisler, J. J. Hamlin, G. R. Stewart, R. G. Hennig, and P. J. Hirschfeld, Structural transitions, octahedral rotations, and electronic properties of $\mathrm{A_{3}Ni_{2}O_{7}}$ rare-earth nickelates under high pressure (2023), arXiv:2309.15078 [cond-mat.supr-con] .
  * Zhang  _et al._ [2023b] Y. Zhang, L.-F. Lin, A. Moreo, and E. Dagotto, Electronic structure, orbital-selective behavior, and magnetic tendencies in the bilayer nickelate superconductor $\mathrm{La_{3}Ni_{2}O_{7}}$ under pressure (2023b), arXiv:2306.03231 [cond-mat.supr-con] .
  * Sakakibara  _et al._ [2023] H. Sakakibara, N. Kitamine, M. Ochi, and K. Kuroki, Possible high $\mathrm{T_{c}}$ superconductivity in $\mathrm{La_{3}Ni_{2}O_{7}}$ under high pressure through manifestation of a nearly-half-filled bilayer hubbard model (2023), arXiv:2306.06039 [cond-mat.supr-con] .
  * Gu  _et al._ [2023] Y. Gu, C. Le, Z. Yang, X. Wu, and J. Hu, Effective model and pairing tendency in bilayer Ni-based superconductor $\mathrm{La_{3}Ni_{2}O_{7}}$ (2023), arXiv:2306.07275 [cond-mat.supr-con] .
  * Zhang  _et al._ [2023c] Y. Zhang, L.-F. Lin, A. Moreo, T. A. Maier, and E. Dagotto, Electronic structure, magnetic correlations, and superconducting pairing in the reduced ruddlesden-popper bilayer $\mathrm{La_{3}Ni_{2}O_{6}}$ under pressure: different role of $d_{3z^{2}-r^{2}}$ orbital compared with la3ni2o7 (2023c), arXiv:2310.17871 [cond-mat.supr-con] .
  * Yang  _et al._ [2023] J. Yang, H. Sun, X. Hu, Y. Xie, T. Miao, H. Luo, H. Chen, B. Liang, W. Zhu, G. Qu, C.-Q. Chen, M. Huo, Y. Huang, S. Zhang, F. Zhang, F. Yang, Z. Wang, Q. Peng, H. Mao, G. Liu, Z. Xu, T. Qian, D.-X. Yao, M. Wang, L. Zhao, and X. J. Zhou, Orbital-dependent electron correlation in double-layer Nickelate $\mathrm{La_{3}Ni_{2}O_{7}}$ (2023), arXiv:2309.01148 [cond-mat.supr-con] .
  * Keimer  _et al._ [2015] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006).
  * Sachdev [2003] S. Sachdev, Rev. Mod. Phys. 75, 913 (2003).
  * Jiang  _et al._ [2023] K. Jiang, Z. Wang, and F.-C. Zhang, High temperature superconductivity in $\mathrm{La_{3}Ni_{2}O_{7}}$ (2023), arXiv:2308.06771 [cond-mat.supr-con] .
  * Luo  _et al._ [2023b] Z. Luo, B. Lv, M. Wang, W. Wú, and D.-X. Yao, High-tc superconductivity in $\mathrm{La_{3}Ni_{2}O_{7}}$ based on the bilayer two-orbital t-j model (2023b), arXiv:2308.16564 [cond-mat.supr-con] .
  * Lu  _et al._ [2023d] C. Lu, Z. Pan, F. Yang, and C. Wu, Interlayer coupling driven high-temperature superconductivity in $\mathrm{La_{3}Ni_{2}O_{7}}$ under pressure (2023d), arXiv:2307.14965 [cond-mat.supr-con] .
  * Lu  _et al._ [2023e] D.-C. Lu, M. Li, Z.-Y. Zeng, W. Hou, J. Wang, F. Yang, and Y.-Z. You, Superconductivity from doping symmetric mass generation insulators: Application to $\mathrm{La_{3}Ni_{2}O_{7}}$ under pressure (2023e), arXiv:2308.11195 [cond-mat.str-el] .
  * Qin and Yang [2023] Q. Qin and Y.-F. Yang, Phys. Rev. B 108, L140504 (2023).
  * feng Yang  _et al._ [2023] Y. feng Yang, G.-M. Zhang, and F.-C. Zhang, Minimal effective model and possible high-$\mathrm{T_{c}}$ mechanism for superconductivity of $\mathrm{La_{3}Ni_{2}O_{7}}$ under high pressure (2023), arXiv:2308.01176 [cond-mat.supr-con] .
  * Shen  _et al._ [2023] Y. Shen, M. Qin, and G.-M. Zhang, Effective bi-layer model hamiltonian and density-matrix renormalization group study for the high-tc superconductivity in la3ni2o7 under high pressure (2023), arXiv:2306.07837 [cond-mat.str-el] .
  * Christiansson  _et al._ [2023] V. Christiansson, F. Petocchi, and P. Werner, Correlated electronic structure of $\mathrm{La_{3}Ni_{2}O_{7}}$ under pressure (2023), arXiv:2306.07931 [cond-mat.str-el] .
  * Shilenko and Leonov [2023] D. A. Shilenko and I. V. Leonov, Phys. Rev. B 108, 125105 (2023).
  * Oh and Zhang [2023] H. Oh and Y.-H. Zhang, Type II t-J model and shared antiferromagnetic spin coupling from hund’s rule in superconducting la3ni2o7 (2023), arXiv:2307.15706 [cond-mat.str-el] .
  * Qu  _et al._ [2023] X.-Z. Qu, D.-W. Qu, J. Chen, C. Wu, F. Yang, W. Li, and G. Su, Bilayer $t$-$j$-$j_{\perp}$ model and magnetically mediated pairing in the pressurized nickelate la3ni2o7 (2023), arXiv:2307.16873 [cond-mat.str-el] .
  * Liao  _et al._ [2023] Z. Liao, L. Chen, G. Duan, Y. Wang, C. Liu, R. Yu, and Q. Si, Electron correlations and superconductivity in $\mathrm{La_{3}Ni_{2}O_{7}}$ under pressure tuning (2023), arXiv:2307.16697 [cond-mat.supr-con] .
  * Liu  _et al._ [2023] Y.-B. Liu, J.-W. Mei, F. Ye, W.-Q. Chen, and F. Yang, The s±-wave pairing and the destructive role of apical-oxygen deficiencies in $\mathrm{La_{3}Ni_{2}O_{7}}$ under pressure (2023), arXiv:2307.10144 [cond-mat.supr-con] .
  * Chen  _et al._ [2023a] X. Chen, P. Jiang, J. Li, Z. Zhong, and Y. Lu, Critical charge and spin instabilities in superconducting $\mathrm{La_{3}Ni_{2}O_{7}}$ (2023a), arXiv:2307.07154 [cond-mat.supr-con] .
  * Zhang  _et al._ [2023d] J.-X. Zhang, H.-K. Zhang, Y.-Z. You, and Z.-Y. Weng, Strong pairing originated from an emergent $\mathbb{Z}_{2}$ berry phase in $\mathbf{La}_{3}\mathbf{Ni}_{2}\mathbf{O}_{7}$ (2023d), arXiv:2309.05726 [cond-mat.str-el] .
  * Schlömer  _et al._ [2023] H. Schlömer, U. Schollwöck, F. Grusdt, and A. Bohrdt, Superconductivity in the pressurized nickelate $\mathrm{La_{3}Ni_{2}O_{7}}$ in the vicinity of a BEC-BCS crossover (2023), arXiv:2311.03349 [cond-mat.str-el] .
  * Lu  _et al._ [2023f] C. Lu, Z. Pan, F. Yang, and C. Wu, Interplay of two $\mathrm{E_{g}}$ orbitals in superconducting $\mathrm{La_{3}Ni_{2}O_{7}}$ under pressure (2023f), arXiv:2310.02915 [cond-mat.supr-con] .
  * Chen  _et al._ [2023b] J. Chen, F. Yang, and W. Li, Orbital-selective superconductivity in the pressurized bilayer nickelate $\mathrm{La_{3}Ni_{2}O_{7}}$: An infinite projected entangled-pair state study (2023b), arXiv:2311.05491 [cond-mat.str-el] .
  * Blankenbecler  _et al._ [1981] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
  * Hirsch  _et al._ [1981] J. E. Hirsch, D. J. Scalapino, R. L. Sugar, and R. Blankenbecler, Phys. Rev. Lett. 47, 1628 (1981).
  * Hirsch [1985] J. E. Hirsch, Phys. Rev. B 31, 4403 (1985).
  * Assaad and Evertz [2008] F. Assaad and H. Evertz, World-line and determinantal quantum monte carlo methods for spins, phonons and electrons, in  _Computational Many-Particle Physics_, edited by H. Fehske, R. Schneider, and A. Weiße (Springer Berlin Heidelberg, Berlin, Heidelberg, 2008) pp. 277–356.
  * Li and Yao [2019] Z.-X. Li and H. Yao, Annual Review of Condensed Matter Physics 10, 337 (2019).
  * Loh  _et al._ [1990] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino, and R. L. Sugar, Phys. Rev. B 41, 9301 (1990).
  * Troyer and Wiese [2005] M. Troyer and U.-J. Wiese, Phys. Rev. Lett. 94, 170201 (2005).
  * Wu and Zhang [2005] C. Wu and S.-C. Zhang, Phys. Rev. B 71, 155115 (2005).
  * Berg  _et al._ [2012] E. Berg, M. A. Metlitski, and S. Sachdev, Science 338, 1606 (2012).
  * Li  _et al._ [2015] Z.-X. Li, Y.-F. Jiang, and H. Yao, Phys. Rev. B 91, 241117 (2015).
  * Li  _et al._ [2016] Z.-X. Li, Y.-F. Jiang, and H. Yao, Phys. Rev. Lett. 117, 267002 (2016).
  * Wang  _et al._ [2015] L. Wang, Y.-H. Liu, M. Iazzi, M. Troyer, and G. Harcos, Phys. Rev. Lett. 115, 250601 (2015).
  * Wei  _et al._ [2016] Z. C. Wei, C. Wu, Y. Li, S. Zhang, and T. Xiang, Phys. Rev. Lett. 116, 250601 (2016).
  * Zhai  _et al._ [2009] H. Zhai, F. Wang, and D.-H. Lee, Phys. Rev. B 80, 064517 (2009), arXiv:0905.1711 [cond-mat.supr-con] .
  * Ma  _et al._ [2022] T. Ma, D. Wang, and C. Wu, Phys. Rev. B 106, 054510 (2022).
  * Zhang and Sachdev [2020a] Y.-H. Zhang and S. Sachdev, Physical Review Research 2, 023172 (2020a), arXiv:2001.09159 [cond-mat.str-el] .
  * Zhang and Sachdev [2020b] Y.-H. Zhang and S. Sachdev, Phys. Rev. B 102, 155124 (2020b), arXiv:2006.01140 [cond-mat.str-el] .
  * Nikolaenko  _et al._ [2021] A. Nikolaenko, M. Tikhanovskaya, S. Sachdev, and Y.-H. Zhang, Phys. Rev. B 103, 235138 (2021), arXiv:2103.05009 [cond-mat.str-el] .
  * Bohrdt  _et al._ [2021] A. Bohrdt, L. Homeier, C. Reinmoser, E. Demler, and F. Grusdt, Annals of Physics 435, 168651 (2021), arXiv:2107.08043 [cond-mat.quant-gas] .
  * Bohrdt  _et al._ [2022] A. Bohrdt, L. Homeier, I. Bloch, E. Demler, and F. Grusdt, Nature Physics 18, 651 (2022), arXiv:2108.04118 [cond-mat.str-el] .
  * Hirthe  _et al._ [2023] S. Hirthe, T. Chalopin, D. Bourgund, P. Bojović, A. Bohrdt, E. Demler, F. Grusdt, I. Bloch, and T. A. Hilker, Nature (London) 613, 463 (2023), arXiv:2203.10027 [cond-mat.quant-gas] .
  * Feldbacher and Assaad [2001] M. Feldbacher and F. F. Assaad, Phys. Rev. B 63, 073105 (2001).
  * Kanaya and Kaya [1995] K. Kanaya and S. Kaya, Phys. Rev. D 51, 2404 (1995).
  * Sandvik [1998] A. W. Sandvik, Phys. Rev. B 57, 10287 (1998).
  * Shao and Sandvik [2023] H. Shao and A. W. Sandvik, Physics Reports 1003, 1 (2023), progress on stochastic analytic continuation of quantum Monte Carlo data.

##  I Supplemental Material for squarelattice-SMG

###  I.1 Section I: Projector determinant Quantum Monte-Carlo algorithm

In this work, we implement projector Quantum Monte-Carlo (PQMC) simulation,
which is an unbiased numerical algorithm for studying the ground state of
interacting fermionic systems. The key scheme of PQMC is that one can access
the accurate ground-state properties through imaginary-time evolution on a
trial wave function, as long as the trial wave function is not orthogonal to
the real ground-state wave function which is generically satisfied for a
quantum many-body Hamiltonian. The ground-state wave function is accessed as:
$|\psi_{G}\rangle=\lim_{\Theta\rightarrow\infty}e^{-\Theta\hat{H}}|\psi_{T}\rangle$,
where $\Theta$ is the projective parameter, $|\psi_{G}\rangle$ is the ground
state of Hamiltonian under consideration and $|\psi_{T}\rangle$ is trail wave
function. For simplicity, the ground-state wave function of given non-
interacting Hamiltonian is generally chosen as the trial wave function, which
is a Slater-determinate. Under this scheme, the corresponding ground-state
expectation value of observable is achieved as:

$\displaystyle\langle\hat{O}\rangle=\frac{\langle\psi_{G}|\hat{O}|\psi_{G}\rangle}{\langle\psi_{G}|\psi_{G}\rangle}=\frac{\langle\psi_{T}|e^{-\Theta\hat{H}}\hat{O}e^{-\Theta\hat{H}}|\psi_{T}\rangle}{\langle\psi_{T}|e^{-2\Theta\hat{H}}|\psi_{T}\rangle}$
(S1)

In PQMC, one usually employs Trotter decomposition to discretize imaginary
time and Hubbard-Stratonovich (H-S) transformation to decouple the interacting
terms in the Hamiltonian. The Trotter decomposition is expressed as:

$\displaystyle e^{-\Theta\hat{H}}=\lim_{\Delta\tau\rightarrow
0}\left[e^{-\Delta\tau\hat{H}_{t}}e^{-\Delta\tau\hat{H}_{I}}\right]^{N\tau}$
(S2)

where $\Delta\tau=\frac{\Theta}{N_{\tau}}$ is the imaginary-time slice,
$\hat{H}_{t}$ and $\hat{H}_{I}$ is the kinetic and interacting part of
Hamiltonian, respectively. In our model, the kinetic and interacting parts of
Hamiltonian are written as:

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle\hat{H}_{t}+\hat{H}_{I}$
(S3) $\displaystyle\hat{H}_{t}$ $\displaystyle=$
$\displaystyle-t\sum_{\left\langle
ij\right\rangle\sigma\alpha}(c^{\dagger}_{i\sigma\alpha}c_{j\sigma\alpha}+\mathrm{h.c.})$
(S4) $\displaystyle\hat{H}_{I}$ $\displaystyle=$ $\displaystyle
J\sum_{i}\vec{S}_{i,1}\cdot\vec{S}_{i,2}$ (S5)

where $c^{\dagger}_{i\sigma}$ creates an electron on site $i$ with spin
polarization $\sigma=\uparrow,\downarrow$ and layer index $\alpha=1,2$. Here
$t$ is the hopping amplitude of electrons between NN sites and $J$ denotes the
strength of inter-layer AFM Heisenberg interaction. To facilitate the Hubbard-
Strotonovich transformation on the interaction, the interacting term is
rewritten as:

$\displaystyle\hat{H}_{I}=\frac{J}{4}\sum_{i}\left[\left(S_{i1}^{x}+S_{i2}^{x}\right)^{2}+\left(S_{i1}^{y}+S_{i2}^{y}\right)^{2}+\left(S_{i1}^{z}+S_{i2}^{z}\right)^{2}-\left(S_{i1}^{x}-S_{i2}^{x}\right)^{2}-\left(S_{i1}^{y}-S_{i2}^{y}\right)^{2}-\left(S_{i1}^{z}-S_{i2}^{z}\right)^{2}\right]$
(S6)

In our simulation, the discrete H-S transformation is implemented[70]:

$\displaystyle e^{\Delta\tau\lambda A^{2}}=\sum_{l=\pm 1,\pm
2}\gamma\left(l\right)e^{\sqrt{\Delta\tau\lambda}\eta\left(l\right)A}+\mathcal{O}\left(\Delta\tau^{4}\right)$
(S7)

where the auxiliary field $\eta$ and $\gamma$ are taken the following values:

$\displaystyle\gamma\left(\pm 1\right)=1+\sqrt{6}/3,$ (S8)
$\displaystyle\gamma\left(\pm 2\right)=1-\sqrt{6}/3,$ (S9)
$\displaystyle\eta\left(\pm 1\right)=\pm\sqrt{2\left(3-\sqrt{6}\right)},$
(S10) $\displaystyle\eta\left(\pm
2\right)=\pm\sqrt{2\left(3+\sqrt{6}\right)}.$ (S11)

Under the transformation, the exponential of interacting term is expressed as:

$\displaystyle
e^{-\Delta\tau\hat{H}_{I}}=\prod_{{\sigma_{1}}=x,y,z}\left[\sum_{a_{1}^{\sigma_{1}}=\pm
1,\pm
2}\gamma\left(a_{1}^{\sigma_{1}}\right)e^{\hat{h}_{1}^{\sigma_{1}}}\right]\times\prod_{{\sigma_{2}}=x,y,z}\left[\sum_{a_{2}^{\sigma_{2}}=\pm
1,\pm
2}\gamma\left(a_{2}^{\sigma_{2}}\right)e^{\hat{h}_{2}^{\sigma_{2}}}\right]$
(S12)

where $\hat{h}_{1}^{\sigma_{1}}=i\sqrt{\Delta\tau
J/4}\eta(a_{1}^{\sigma_{1}})(S_{i1}^{\sigma_{1}}+S_{i2}^{\sigma_{1}})$,
$\hat{h}_{2}^{\sigma_{2}}=\sqrt{\Delta\tau
J/4}\eta(a_{2}^{\sigma_{2}})(S_{i1}^{\sigma_{2}}-S_{i2}^{\sigma_{2}})$ and
$a_{1}^{\sigma_{1}},a_{2}^{\sigma_{2}}=\pm 1,\pm 2$. Then, it is
straightforward to employ the standard procedure of PQMC to evaluate the
ground-state expectation value of observables as shown in Eq. (S1)[70, 71].

Figure S1: (a)-(d) The results of structure factors for SDW, CDW, EC and SC.
(e)-(h) The results of correlation-length ratios for SDW, CDW, EC and SC. The
system sizes are $L=8,10,12,14,16,18$.

###  I.2 Section II: The absence of sign problem

In this section, we demonstrate the absence of sign problem in the procedure
of PQMC as illustrated in the last section. After H-S transformation, the
decoupled interacting terms are expressed as:

$\displaystyle\hat{h}_{1}^{\sigma}$ $\displaystyle=$ $\displaystyle
i\sqrt{\Delta\tau J/4}\eta(a)(S_{i1}^{\sigma}+S_{i2}^{\sigma})$ (S13)
$\displaystyle\hat{h}_{2}^{\sigma}$ $\displaystyle=$
$\displaystyle\sqrt{\Delta\tau J/4}\eta(a)(S_{i1}^{\sigma}-S_{i2}^{\sigma})$
(S14)

where $\sigma=x,y,z$, $a=\pm 1,\pm 2$ and
$\vec{S}_{il}=(S^{x}_{il},S^{y}_{il},S^{z}_{il})=\frac{1}{2}c^{\dagger}_{i\alpha
l}\vec{\sigma}^{\alpha\beta}c_{i\beta l}$, with
$\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ being Pauli matrix. The
definition of $\eta(a)$ is shown in Eq. (S11). Remarkably, there exists a non-
unitary transformation $\hat{T}=i\tau_{x}\sigma_{y}K$, where $\tau_{x}$ is the
first Pauli matrix operating in the space of layer, $\sigma_{y}$ is the second
Pauli matrix operating in the spin space, and $K$ is the operation of complex
conjugation. It is straightforward to show that the non-interacting
Hamiltonian and the interacting terms after H-S decoupling are invariant under
such non-unitary transformation:

$\displaystyle\hat{T}\hat{H}_{t}\hat{T}^{-1}$ $\displaystyle=$ $\displaystyle
H_{t}$ (S15) $\displaystyle\hat{T}\hat{h}_{1}^{\sigma}\hat{T}^{-1}$
$\displaystyle=$ $\displaystyle\hat{h}_{1}^{\sigma}$ (S16)
$\displaystyle\hat{T}\hat{h}_{2}^{\sigma}\hat{T}^{-1}$ $\displaystyle=$
$\displaystyle\hat{h}_{2}^{\sigma}$ (S17)

where $\sigma=x,y,z$ and $\hat{H}_{t}$ is defined in Eq. (S4). Because
$\hat{T}$ obeys $\hat{T}^{2}=-1$, the model is sign-problem-free in QMC
simulation owing to the principle introduced in [74].

Figure S2:  Finite-size scaling analysis of SDW (a), CDW (b), EC (c) and SC
(d). The results are fitted by a second-order polynomial function of $1/L$.
The system sizes are $L=8,10,12,14,16,18$.

###  I.3 III. The definition and results of structure factors and correlation-
length ratios

To characterize the symmetry spontaneously breaking phases and the phase
transitions, we compute structure factors of various symmetry breaking orders,
which are defined as following:

$\displaystyle
S(\vec{Q},L)=\frac{1}{L^{4}}\sum_{ij}\left\langle\hat{O}(i)\hat{O}(j)\right\rangle
e^{i\vec{Q}\cdot(\vec{r}_{i}-\vec{r}_{j})}$ (S18)

where $\hat{O}(i)$ is the symmetry breaking operator field located at site
$i$, $L$ is the linear system size and $\vec{Q}$ is the momentum of symmetry
breaking order under consideration. In our cases, we consider the spin-density
wave (SDW) and charge-density wave (CDW) orders. The corresponding operator
fields are defined as:

$\displaystyle
O_{\rm{SDW}}(i)=\sum_{\alpha}c^{\dagger}_{i\uparrow\alpha}c_{i\uparrow\alpha}-c^{\dagger}_{i\downarrow\alpha}c_{i\downarrow\alpha}\quad
O_{\rm{CDW}}(i)=\sum_{\alpha}c^{\dagger}_{i\uparrow\alpha}c_{i\uparrow\alpha}+c^{\dagger}_{i\downarrow\alpha}c_{i\downarrow\alpha}$
(S19)

where $\alpha$ denotes the index of layer, $\uparrow$/$\downarrow$ is spin
polarization. Here we only consider spin order in z-direction without the loss
of generality because the model in Hamiltonian Eq. (1) preserves spin $SU(2)$
symmetry. For the bilayer model, the interlayer exciton condensation (EC)
order is also considered in our study. The operator field is defined as:

$\displaystyle O_{\rm{EC}}(i)=\sum_{\sigma}c^{\dagger}_{i\sigma 1}c_{i\sigma
2}$ (S20)

where $\sigma=\uparrow/\downarrow$ is the spin polarization and $1/2$ is the
layer index. Additionally, we consider the order of superconductivity (SC). In
our case, the dominant SC channel is the interlayer on-site spin singlet
pairing, with the corresponding operator field defined as:

$\displaystyle O_{\rm{SC}}(i)=c^{\dagger}_{i\uparrow
1}c^{\dagger}_{i\downarrow 2}-c^{\dagger}_{i\downarrow
1}c^{\dagger}_{i\uparrow 2}$ (S21)

For SC ordering, the dominant ordering momentum is $(0,0)$. For the order
parameter field involving particle-hole instabilities, including SDW, CDW and
EC, the dominant ordering carries momentum $(\pi,\pi)$ in our cases owing to
the Fermi surface nesting of non-interacting band structure. We also confirm
them via scrutinizing the momentum distribution of structure factors for
various orders.

To reveal the possible SSB in model Eq. (1), we evaluate the structure factors
at peaked momentum for different system sizes, as presented in Fig. S1(a)-(d).
In the weakly interacting regime, the structure factors for EC and SC orders
are obviously larger than the ones for SDW and CDW. The finite-size scaling of
structure factors for SDW, CDW, EC and SC are shown in Fig. S2. The results
clearly show that SDW and CDW long-ranged order is absent in the entire
interaction regime of the model Eq. (1). The EC and inter-layer singlet SC
order exist in the weakly interacting regime and vanish when interaction is
strong. The degeneracy of EC and SC order is also explicitly witnessed by the
equivalence of structure factors.

Figure S3:  (a)-(e) Finite-size scaling of single-particle gap $\Delta_{\rm
sp}(L)$ for various interaction strength $J$. The $\Delta_{\rm sp}$ is
extracted at the momentum $(\frac{\pi}{2},\frac{\pi}{2})$, at which the gap
size is minimum. (f) The results of single-particle gap at thermodynamic limit
obtained by finite-size scaling as a function of $J$.

To further verify the conclusions, we evaluate the correlation-length ratio
for different system sizes, which is believed to bear weaker finite-size
effects compared with the direct extrapolation of structure factors. The
correlation-length ratio is defined as:

$\displaystyle R(L)=\frac{S(\vec{Q}_{m},L)}{S(\vec{Q}_{m}+\delta\vec{q},L)}-1$
(S22)

where $S(\vec{Q}_{m},L)$ is the structure factor of the order under
consideration, $\vec{Q}_{m}$ is the ordering momentum, and
$\delta\vec{q}=(\frac{2\pi}{L},\frac{2\pi}{L})$ denotes the minimum momentum
shift from ordering momentum $\vec{Q}_{m}$ on a lattice with linear system
size $L$. For the disordered phase of the order parameter under consideration,
$R(L)$ decreases with $L$ and $R(L)\rightarrow 0$ in the thermodynamic limit,
whereas for the long-ranged ordered phase, $R(L)$ increases with $L$ and
$R(L)\rightarrow\infty$ in the thermodynamic limit. At putative quantum
critical point, RG-invariant quantity $R(L)$ is independent on system size as
$L$ is sufficiently large. The results of correlation-length ratios for SDW,
CDW, EC and SC are presented in Fig. S1(e)-(h). For SDW and CDW order, $R(L)$
decreases with system size in the whole interacting regime. For EC and SC
order, a quantum phase transition from long-ranged ordered phase to the
disordered phase occurs with increasing interaction strength. The degeneracy
between EC and SC order is further verified by the equivalence of correlation-
length ratios for two orders.

###  I.4 IV. The finite-size scaling of single-particle gap

To reveal the nature of SMG insulating phase, a prerequisite is the single-
particle excitation is gaped. In this section, we present the finite-size
scaling results of single-particle gap. In QMC simulation, single-particle gap
is accessible by evaluating the time-dependent single-particle Green’s
function $G(\vec{k},\tau)=\left\langle
c_{\vec{k}}(0)c^{\dagger}_{\vec{k}}(\tau)\right\rangle$, which exhibits the
scaling behavior $G(\vec{k},\tau)\sim e^{-\Delta_{sp}(\vec{k})\tau}$ as
imaginary-time $\tau$ is long enough. Hence, the single-particle gap at
momentum $\vec{k}$ is extracted from the scaling of $G(\vec{k},\tau)$ versus
$\tau$. According to the numerical results of single-particle gap in the
Brillouin zoom, it is shown that the single-particle gap is minimum at
momentum $(\pm\frac{\pi}{2},\pm\frac{\pi}{2})$. Fixing
$\vec{k}=(\frac{\pi}{2},\frac{\pi}{2})$, we compute time-dependent Green’s
function for different system sizes and various $J$, and extract the single-
particle gap from the scaling behavior of $G(\vec{k},\tau)$. The results of
single-particle gap are included in Fig. S3. In the entire interaction
strength regime we consider, the single-particle gap at thermodynamic limit is
finite, indicating that the ground state is an insulator. Particularly, in the
strongly interacting regime $J>3.24$ where the ground state does not involve
any spontaneous symmetry breaking, the single-particle exciton features
pronounced gap, demonstrating the nature of a SMG insulating phase.