# Universality and Critical Exponents of the Fermion Sign Problem

R. Mondaini<EMAIL_ADDRESS>S. Tarat<EMAIL_ADDRESS>Beijing
Computational Science Research Center, Beijing 100193, China  R.T. Scalettar
<EMAIL_ADDRESS>Department of Physics and Astronomy, University
of California, Davis, CA 95616, USA

###### Abstract

Initial characterizations of the fermion sign problem focused on its evolution
with spatial lattice size $L$ and inverse temperature $\beta$, emphasizing the
implications of the exponential nature of the decay of the average sign
$\langle{\cal S}\rangle$ for the complexity of its solution and associated
limitations of quantum Monte Carlo studies of strongly correlated materials.
Early interest was also on the evolution of $\langle{\cal S}\rangle$ with
density $\rho$, either because commensurate filling is often associated with
special symmetries for which the sign problem is absent, or because particular
fillings are often primary targets, e.g. those densities which maximize
superconducting transition temperature (the top of the ‘dome’ of cuprate
systems). Here we describe a new analysis of the sign problem which
demonstrates that the spin-resolved sign $\langle{\cal S}_{\sigma}\rangle$
already possesses signatures of universal behavior traditionally associated
with order parameters, even in the absence of symmetry protection that makes
$\langle{\cal S}\rangle=1$. When appropriately scaled, $\langle{\cal
S}_{\sigma}\rangle$ exhibits universal crossings and data collapse. Moreover,
we show these behaviors occur in the vicinity of quantum critical points of
three well understood models, exhibiting either second order or Kosterlitz-
Thouless phase transitions. Our results pave the way for using the average
sign as a minimal correlator that can potentially describe quantum criticality
in a variety of fermionic many-body problems.

##  I Introduction

The sign problem is the fundamental obstacle that prevents accurate
computations in a variety of problems of quantum correlated matter. In evading
the ‘exponential wall’ that precludes the application of unbiased methods,
such as exact diagonalization [1, 2] and matrix-product-states-based
algorithms [3], for large systems or arbitrary dimensions, quantum Monte Carlo
techniques [4, 5] have in principle the potential to solve fundamental
questions, including understanding $d$-wave pairing mechanisms of repulsive
fermions, for example [6, 7]. Yet, the fact that the importance sampling of
quantum configurations is not constrained to render positive weights severely
limits its applicability in the most salient class of quantum problems of
interest.

Since generic solutions are not always available, a common approach relies in
restricting computations to regimes where the sign problem is still well
behaved, allowing the extraction of statistically convergent quantities.
Recent developments based on finding a local basis that mitigates it [8, 9,
10], or fine-tuned Hubbard-Stratonovich transformations to delay its onset
[11], have been very useful but do not provide an overarching circumvention
scheme. Although solutions to this conjectured NP-hard problem [12] are
unlikely to be discovered anytime soon, a precise investigation of the onset
of the sign problem, and its relation to the physics of the underlying studied
Hamiltonian, is much less explored, and potentially of great impact. Further
investigations have advanced the idea that this connection does exist, and in
different models the appearance of the sign problem seems coupled to the
manifestation of quantum critical behavior [13, 14].

When several non-hybridizing fermionic (spin) species are present, the sign
problem involves the product of contributions from each component. Previous
studies have focused on this product, both because it is required to weight
physical observables, and also because in several important cases the product
is better behaved than its constituents. Furthermore, existing work has
generally concentrated on the ‘scaling behavior’ in the sense of large space-
time systems, i.e., how the sign evolves as the inverse temperature
$\beta\rightarrow\infty$ and spatial size $L\rightarrow\infty$.

In this manuscript, we introduce two novel aspects of the study of the sign
problem, and show that they constitute a powerful new approach to using
quantum simulations to explore many-body physics. First, we analyze the ‘spin
resolved sign’, and argue that the usual approach of examining the average
product of the sign of the individual weights can actually obscure physical
content inherent when spin-resolution is used. Second, we examine the behavior
of the sign near phase transitions, both those which occur as quantum critical
points, through the variation of a parameter in the Hamiltonian, and thermal
phase transitions which occur as temperature $T$ is lowered.

Taken together, we demonstrate that the spin-resolved sign can be used to
locate phase transitions and determine critical exponents. Furthermore, it has
the potential to do so even more accurately than traditional observables
${\cal A}$ such as spin, charge, and pairing correlations. The reason is that
these latter quantities require a precise measurement of ratios $\langle{\cal
A}{\cal S}\rangle/\langle{\cal S}\rangle$ of quantities with increasing
fluctuations originating both inherently in the physics (response functions
are themselves measurements of fluctuations) and in the vanishing of the sign.
The (spin-resolved) sign, by itself, is thus a less noisy ‘observable’, if it
can be shown, as we do here, that it holds information about criticality.

##  II The SU(2) Honeycomb Hubbard model

We initially investigate the spinful Hubbard model on a honeycomb lattice with
$N=2L^{2}$ sites,

$\hat{H}=-t\sum_{\langle
ij\rangle\,\sigma}\hat{c}^{\dagger}_{i\sigma}\hat{c}^{\phantom{\dagger}}_{j\sigma}+U\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}-\mu\sum_{i,\sigma}\hat{n}_{i\sigma},$
(1)

where $\hat{c}^{\dagger}_{i\sigma}$ ($\hat{c}^{\phantom{\dagger}}_{i\sigma}$)
creates (annihilates) a fermion at site $i$ with spin $\sigma$, and
$\hat{n}_{i\sigma}$ is the number density operator. With an increasing
magnitude of the ratio of the interactions amplitude to the hopping scale,
$U/t$, the ground-state at half-filling (chemical potential $\mu=U/2$)
exhibits a continuous phase transition from a Dirac semi-metal to a Mott
insulator featuring antiferromagnetic order. This transition, described by an
effective quantum field-theory model (Gross-Neveu) [15], belongs to the chiral
Heisenberg universality class, and has been characterized in numerics in a
variety of fermionic lattice models [16, 17, 18, 19, 20, 21].

High precision computation of the critical interaction in (1) yields
$U_{c}/t=3.78-3.87$, with critical exponent in the range $\nu=0.84-1.02$ [16,
17, 18, 19] associated to the divergence of the correlation length in the
vicinity of the critical point $\xi\propto|U-U_{c}|^{-\nu}$. These values are
obtained via the scaling of physical observables: the staggered magnetization
order parameter [16, 17, 19], single particle gap [17], and quasiparticle
weight [19, 21]. Here, instead, we propose an analysis based on the average
sign. Although a non-physical observable, tied to the computational method
used, the average sign is however required to compute any physical observable
in a quantum Monte Carlo (QMC) simulation. Our results thus suggest that the
sign problem is inextricably linked to determination of the physics of the
model. In the following two sections, we explore generic aspects of the sign
problem and how they apply to this particular Hamiltonian. Subsequently, we
build on this knowledge to understand quantum and thermal phase transitions in
other fermionic models.

##  III The sign problem

We start by recalling a known scaling form of the average sign in QMC
calculations. It originates from considering the definition in terms of the
weights $W$ of the configurations $\\{x\\}$ sampled in $D$ spatial and 1
imaginary-time dimensions as [22, 12],

$\langle{\cal
S}\rangle=\frac{\sum_{\\{x\\}}W(\\{x\\})}{\sum_{\\{x\\}}|W({\\{x\\}})|}=\frac{{\cal
Z}_{W}}{{\cal Z}_{|W|}}.$ (2)

Here, $W(\\{x\\})=\det M_{\uparrow}(\\{x\\})\cdot\det M_{\downarrow}(\\{x\\})$
is a product of weights of individual fermionic flavors in the case of Eq. (1)
for determinant QMC calculations (see Appendix A for specific definitions in
the various models we study); ${\cal Z}_{W}$ is the partition function of the
original problem in its formulation in $D+1$ dimensions [23, 24], whereas
${\cal Z}_{|W|}$ instead uses the positive definite absolute value of the
weight in order to proceed with the importance sampling in the simulations.
Written in terms of the corresponding free energy densities, $f=[-1/(\beta
N)]\log{\cal Z}$, the average sign thus reduces to $\langle{\cal
S}\rangle=\exp[-\beta N(f_{W}-f_{|W|})]$. Given that
$\sum_{\\{x\\}}W(\\{x\\})\leq\sum_{\\{x\\}}|W({\\{x\\}})|$ and that the free
energy is extensive, it follows that $f_{W}\geq f_{|W|}$, and thus the average
sign exponentially decreases in terms of both real-space and imaginary-time
dimensions [22, 25], if not protected by some symmetry of the problem [26,
27].

Figure 1: Dependence of the average spin resolved sign on the space-imaginary
time dimensions. $\langle{\cal S}_{\sigma}\rangle$ vs. $L_{\tau}$ (left) and
$L$ (right) for the SU(2) honeycomb Hubbard model, before ($U/t=3.6$, top) and
after ($U/t=5.0$, bottom) the putative transition $U_{c}$. For $U<U_{c}$
($U>U_{c}$) the average spin resolved sign grows (reduces) with increasing
lattice size at low temperatures. Here, and elsewhere, error bars denote the
standard error of the mean of independent realizations. The imaginary-time
discretization is $\Delta\tau\equiv\beta/L_{\tau}=1/10$. Figure 2: Scaling
analysis of the spin resolved sign in the honeycomb SU(2) Hubbard model.
Scaling in the vicinity of the best known estimations of the critical point
$U_{c}/t$ [(a) $U/t=3.6$ and (b) $U/t=3.85$] using a re-scaled x-axis
$L_{\tau}/L^{z}$. The insets display a cost function that determines collapse
quality of $\langle{\cal S}_{\sigma}\rangle(L_{\tau},L)$ at different values
of the dynamic critical exponents $z$ (see text). A compilation is given in
the inset of (c) at a range of $U/t$ values; estimations for the critical
interactions from Refs. [18, 19] are marked by the shaded region. (c)
$\langle{\cal S}_{\sigma}\rangle$ vs. $U/t$ and different lattice sizes at
half-filling; the number of imaginary time-slices used roughly preserves the
ratio $L_{\tau}/L^{z}$ fixed, $L_{\tau}=240$, 220, 196, 170 for $L=18,$ 15,
12, 9, respectively; that is we use $z\simeq 0.5$. (d) Scaling using a
functional form $g[uL^{1/\nu}]$ whose critical exponent $\nu$, as obtained by
minimizing the error $\chi^{2}/{\rm d.o.f}$ of a high order polynomial fitting
in the space of parameters $(U_{c},\nu)$. (e) The contour plot of
$\chi^{2}/{\rm d.o.f}$, where the minimum is at $U_{c}/t=3.765$ and $\nu=0.84$
as shown by the star symbol. Recent estimations using physical quantities for
the same model [18, 19] are annotated by the cross markers. Here,
$t\,\Delta\tau=0.1$ is used.

An example of this protection is the case of Eq. (1) at half-filling. Via a
$\downarrow$-spin particle-hole transformation,
$c_{i\downarrow}=(-1)^{i}c_{i\downarrow}^{\dagger}$, where $(-1)^{i}=+1(-1)$
on the $A(B)$ sublattice of the bipartite honeycomb geometry (or any other
bipartite lattice), the weight simplifies to ${\rm const.}\times[\det
M_{\uparrow}(\\{x\\})]^{2}$ for whichever configuration $\\{x\\}$, when using
a spin-decomposed Hubbard-Stratonovich transformation [28, 24]. In general,
however, symmetries that preclude the onset of the sign problem are not
available for most models of interest.

##  IV The spin resolved sign

Although the ‘total’ sign problem has been investigated in detail [22, 25],
the properties of the sign of individual determinants that compose the weight
in models with larger number of local degrees of freedom were much less
explored. Moreover, past work focused on the behavior as $\beta\to\infty$ and
not near the critical point. By systematically computing the average sign of
the determinant of a single spin-species, $\langle{\cal
S}_{\sigma}\rangle\equiv\sum_{\\{x\\}}{\rm sgn}(\det
M_{\sigma}\\{x\\})|W(\\{x\\})|/\sum_{\\{x\\}}|W(\\{x\\})|$, for the problem in
Eq. (1), we have earlier demonstrated [13] (see corresponding Supplementary
Materials) a behavior reminiscent of an order parameter undergoing a typical
phase transition: it displays its maximum value $\langle{\cal
S}_{\sigma}\rangle=1$ in the quantum disordered phase while $\langle{\cal
S}_{\sigma}\rangle\to 0$ in the ordered region ($U>U_{c}$) at sufficiently low
temperatures [see Fig. 2(c)]. The latter occurs in spite of the fact that
$\langle{\cal S}\rangle$ is pinned at one since we take $\mu=U/2$, dictating
thus that in the ordered regime the most likely configurations $\\{x\\}$
display random signs of $\det M_{\sigma}(\\{x\\})$.

This behavior, including a crossing of the curves for different system sizes
at $U\simeq U_{c}$, is suggestive that a scaling function $g$, for the spin-
resolved sign exists, similar in motivation to those used for traditional,
physical observables to characterize quantum criticality:

$\langle{\cal S}_{\sigma}\rangle(u,L,L_{\tau})=g(uL^{1/\nu},L_{\tau}/L^{z}),$
(3)

where $L$ is the linear system size, $L_{\tau}=\beta/\Delta\tau$ is the number
of imaginary-time slices of the inverse temperature in writing down the path
integral ${\cal Z}$, $u=(U-U_{c})/U_{c}$ is the reduced coupling, and $z$ is
the dynamic critical exponent. The second argument comes from the fact that
the $D+1$ lattice is anisotropic in its dimensions (and eventual effective
couplings) [29, 30, 31]. The previous empirical observation determined that
for $u>0$, $g(x,y)\to 0$ when both $x,y$ diverge.

To better understand the limits of $\langle{\cal S}_{\sigma}\rangle$, we
display in Fig. 1 its dependence on both $L_{\tau}$ and $L$. We notice that
the previous expectation that the average total sign exponentially decreases
with the system size is also valid for $\langle{\cal S}_{\sigma}\rangle$, but
provided that $U\gtrsim U_{c}$ and temperatures are sufficiently low
($L_{\tau}\gg L$). If in the quantum disordered phase, however, the spin-
resolved sign increases with growing $N$. As we will see, this contrasting
behavior is fundamental for the identification of the critical interactions
using the average sign of individual determinants.

Figure 3: Scaling analysis of the spin resolved sign in the square lattice
square lattice SU(2) Ionic Hubbard. (a,b) Similar to Fig. 2 (a,b) but for the
SU(2) Ionic Hubbard model on the square lattice, for values of $U/t=2$ (a) and
2.1 (b). (c) $\langle{\cal S}_{\sigma}\rangle$ vs. $U/t$ and different lattice
sizes at half-filling; the number of imaginary time slices is chosen such that
$L_{\tau}/L^{0.5}\simeq 50$. (d) Scaling of the data in (c) using a functional
form $g[uL^{1/\nu}]$, whose critical exponent $\nu=0.97$ and critical
interaction $U_{c}/t=2.058$ are extracted from an analysis (e) as done in Fig.
2(e). Here, all results are obtained at $\Delta/t=0.5$; other values lead to
similar results, but with different $U_{c}/t$ critical values.Here, we use
$t\,\Delta\tau=0.1$.

At first sight, owing to the known Lorentz invariance that emerges at the
critical point [32, 33], the dynamic critical exponent is surmised as $z=1$.
Yet, we do not take this as a starting point, relaxing this assumption to show
that a smaller value in $\langle{\cal S}_{\sigma}\rangle$ actually gives an
optimal scaling. Motivation for $z\neq 1$ is provided in the outlook. We thus
try to scale $\langle{\cal S}_{\sigma}\rangle$ in the vicinity of the critical
point with a functional form $L_{\tau}/L^{z}$, as displayed in Figs. 2(a) and
(b), a procedure which has been argued to improve the scaling of related
quantum models [30, 31]. By defining a cost function ${\cal
C}(z)=\sum_{j}(|y_{j+1}-y_{j}|)/(\max\\{y_{j}\\}-\min\\{y_{j}\\})-1$ [34, 35],
where $y_{j}$ are the values of $\langle{\cal S}_{\sigma}\rangle(L_{\tau},L)$,
ordered according to their $L_{\tau}/L^{z}$ ratio, the dynamic critical
exponent $z^{*}$ that minimizes ${\cal C}$ can be extracted, see insets in
Figs. 2(a) and (b). A compilation of the $z^{*}(U)$ values is given in the
inset of Fig. 2(c), accompanied by a range of recent known predictions of
$U_{c}$ [18, 19]. It is clear that close to the critical point, the scaling
with the second argument of the function $g$ should be taken with $z\simeq
1/2$, for the current range of imaginary-time slices $L_{\tau}$ (or inverse
temperatures $\beta$ with $t\Delta\tau=0.1$) used.

Hence, we use this current estimation to proceed with a scaling in order to
simultaneously obtain the critical exponent $\nu$ and the critical interaction
$U_{c}$. With the $L_{\tau}/L^{z}$ ratio fixed, a clear crossing of the
average sign of individual weights when increasing the lattice size can be
seen [Fig. 2(c)], accurately determining the critical interaction. By using
the functional form of Eq. (3), we obtain the collapse of the average spin-
resolved sign [Fig. 2(d)], yielding a critical exponent $\nu\simeq 0.84$ and
$U_{c}/t\simeq 3.77$. This estimation, obtained by minimizing the error of a
high-order polynomial fit to the data in the space of parameters $(U,\nu)$, is
shown in Fig. 2(e) [See the Supplemental Materials (SM) [36] for a different
method of scaling analysis]. One can contrast these results with recent
estimations using the same model, as in Ref. 18 with $U_{c}/t=3.80(1)$ and
$\nu=0.84(4)$, while in Ref. 19, $U_{c}/t=3.85(1)$ and $\nu=1.02(1)$, both
using a zero-temperature version of the QMC method employed here [4, 37].
While larger system sizes and other finite-corrections may improve our result,
they are already in quite remarkable agreement with the best estimations to
date.

##  V The SU(2) Ionic Hubbard model

The preceding discussion provided quantitative evidence that the average sign
of a single determinant contains precise information about the quantum
criticality in a well studied model; it remains an open question whether this
is general. Here we provide compelling further validation by looking at one of
the simplest models that bypasses the symmetry that prevents the onset of the
sign problem, the ionic Hubbard model on the square lattice [38, 39, 40, 41,
42, 43, 44]. Here $\hat{H}_{\rm
Ionic}=\hat{H}+\Delta\sum_{i\sigma}(-1)^{i}\hat{n}_{i,\sigma}$, adds a
staggered onsite potential proportional to $\Delta$ to the Hamiltonian (1),
which we investigate again at half-filling.

Figure 4: KT-scaling analysis of the SU(2) square lattice attractive Hubbard
model. Color plot of the $s$-wave pair structure factor $P_{s}$ (a) and the
spin resolved sign $\langle{\cal S}_{\sigma}\rangle$ (b) in the $T$ vs. $\mu$
plane on a $L=16$ lattice; (c) and (d) display the same if using the extracted
average density $\langle\hat{n}\rangle$ instead, with black markers depicting
the outputs on the regular grid in $\mu$. (e) and (f) show a cut along the
$\mu/t=-0.75$ line for $P_{s}$ and $\langle{\cal S}_{\sigma}\rangle$ and
different lattice sizes. (g) and (h) use a KT-scaling form to collapse the
curves with (i) and (j) displaying the corresponding cost function ${\cal C}$
of the scaling in the $(T_{c},b)$ parameter space. The KT-scaling is performed
such that $L/\xi_{\rm BKT}$ for $T>T_{c}$ and $-L/\xi_{\rm BKT}$ if $T<T_{c}$.
All data are extracted at $U/t=-4$ and $t\Delta\tau=1/16$.

The qualitative physics of this model at finite $\Delta$ is generally agreed
to display a competition of band insulating ($\Delta\gg U$), Mott insulating
($U\gg\Delta$) and metallic behavior when both interactions and staggered
potential magnitudes are comparable [41, 42]. A recent investigation [13] has
indicated that this correlated metal phase can be qualitatively tracked by the
regime where the average sign of the QMC weights vanishes. We now employ our
new sign scaling method to understand the critical behavior at the transition
from the band-insulator to the metallic phase, the quantitative details of
which are still under debate in the community.

We fix the staggered potential at $\Delta/t=0.5$, while increasing the
interactions $U$ to overcome the externally imposed (i.e. by the one-body
potential in $\hat{H}$) charge density wave induced by $\Delta$. As before, we
estimate the value of the dynamic critical exponent in the vicinity of the
transition, $U_{c}(\Delta=0.5)/t\simeq 2.0$ [41, 42] in Fig. 3 (a,b) by
scaling the $\langle{\cal S}_{\sigma}\rangle$ with $L_{\tau}/L^{z}$, resulting
in $z\simeq 0.5$. Using thus a roughly fixed ratio $L_{\tau}/L^{0.5}$, we
provide the scaling of the spin resolved sign in Fig. 3(c,d). Here
fluctuations are small, and the scaling renders an accurate determination of
the critical interactions driving the band-insulator-to-metal transition
$U_{c}/t\simeq 2.05$ with related critical exponent $\nu=0.97$. The one-
dimensional version of this model has been extensively studied via numerics
[45, 46, 47], and a field-theory close to the quantum critical points exists
[48]. The transition where the band-insulating phase ends, with its externally
imposed charge density wave giving way to a dimerized bond-ordered wave
insulator, belongs to the 2D Ising universality class in that case. Here in
the two-dimensional model, QMC results point out to a band-insulator to
correlated metallic transition [41, 42], whose universality class is unknown
and where field-theories describing it are currently not available, precluding
a direct comparison of the calculated exponent $\nu$, obtained from the
scaling of $\langle{\cal S}_{\sigma}\rangle$, with existing knowledge.

##  VI The attractive Hubbard model

We now generalize these two results for quantum critical behavior in the
ground-state, to finite-temperature transitions. A well studied example is the
onset of superconductivity in the two-dimensional negative-$U$ SU(2) Hubbard
model, for which, at finite chemical potentials, there is a Kosterlitz-
Thouless (KT) transition at temperature $T_{c}\neq 0$ [49, 50, 51] to a
superconducting phase. As a direct consequence of the often used (charge-
decomposed) Hubbard-Stratonovich transformation [28], the weight matrices
$M_{\sigma}$ are identical, resulting in the complete absence of sign problem,
since the remaining single-particle part of the Hamiltonian is identical for
both spin species. That the total sign is always positive does not preclude
that the average sign of individual weights converge to zero; this can be seen
in Fig. 4(c), which shows this quantity in the $T$ vs. $\mu$ parameter space
on an $L=16$ square lattice. The regime $\langle{\cal S}_{\sigma}\rangle\to 0$
is directly related to the one where the $s$-wave equal-time pair structure
factor
$P_{s}=(1/L^{2})\sum_{i,j}\langle\hat{\Delta}_{i}\hat{\Delta}_{j}^{\dagger}\rangle$
($\hat{\Delta}_{i}\equiv\hat{c}_{i\uparrow}\hat{c}_{j\downarrow}$) is also
large [Fig. 4(a)].

To make this comparison quantitative, Figs. 4(e) and (f) display a finite-size
analysis at a fixed $\mu/t=-0.75$, which gives densities close to
$\langle\hat{n}\rangle=0.5$ at low temperatures [49, 36]. The onset of the
regime at which there is significant size dependence is largely coincident for
both quantities when sweeping down the temperature. Given that the pairing
correlations below $T_{c}$ have an algebraic decay, $C(r)\sim r^{-\eta(T)}$
with $\eta(0)=0$ and $\eta(T_{c})=0.25$, one obtains that the two-dimensional
pair structure factor in a finite system of linear size $L$ scales as
$P_{s}(L)=L^{2-\eta(T)}f(L/\xi_{\rm BKT})$, with $\xi_{\rm
BKT}\propto\exp[b/\sqrt{|T-T_{c}|}]$ [50, 51], where $b$ is a non-universal
constant.

Figure 4(g) gives this scaling using a temperature-adjusted $\eta(T)$ [49,
36]. The best parameters $b$ and $T_{c}$ are extracted via the cost function
${\cal C}_{P_{s}}$, as before, shown in Fig. 4(i). Based on this, we similarly
show a scaling of the average sign of individual weights $\langle{\cal
S}_{\sigma}\rangle$ to a KT form, see Fig. 4(h), with corresponding cost
function ${\cal C}_{\langle{\cal S}_{\sigma}\rangle}$ displayed in Fig. 4(j).
Both quantities scale remarkably precisely; the minor discrepancy in $T_{c}$
[$T_{c}^{P_{s}}=0.15(2)$ and $T_{c}^{\langle{\cal
S}_{\sigma}\rangle}=0.18(3)$] can be accounted by the relatively wide
temperature-region in which ${\cal C}$ is small. For both quantities, we take
the smallest $b$ given the constraint of best collapse in a smooth curve.
Lastly, we note that allowing for the possibility that the non-universal
parameter $b$ takes different values below and above the transition when
lowering the temperature (requiring thus a multi-dimensional minimization
procedure), may improve convergence of the estimations of $T_{c}$ [35]

##  VII Discussion and Outlook

We have shown that the spin-resolved sign of auxiliary-field QMC simulations
can be used as a quantitative marker of quantum critical behavior. The total
sign also exhibits a similar role, as suggested by the Ionic Hubbard model
results (see SM [36]), but the former has the benefit of being useful when
symmetries prevent the occurrence of an overall sign problem. Our work lays
the foundation for similar investigations of other models, especially ones
that give rise to a phase problem. This can arise either from the presence of
imaginary-terms in the Hamiltonian, as in the Kane-Mele Hubbard model [52,
53], or from the particular decoupling scheme used. That is precisely the case
of SU(2) symmetric Hubbard-Stratonovich transformations [54], but
investigations in Appendix D show that the averaged spin-resolved phase
similarly tracks the onset of the ordered regime when approaching the
thermodynamic limit for the SU(2) honeycomb Hubbard model. Furthermore, other
Hamiltonians, such as the spinless fermion Hubbard model in either the
honeycomb [55] or square-lattice with a $\pi$-flux, which in the Majorana
basis evade the sign problem [56, 57], can be studied by examining the average
sign of the Pfaffian of a single weight in that basis 111The original total
weight given by a single determinant in the single-particle basis is converted
to a product of weights related by complex conjugation when using the Majorana
basis to describe the Hamiltonian, hence no sign (phase) problem., similar to
what we have done here [59]. In our results, the investigation of these three
important models emphasizes that the sign of the determinants, interpreted as
a minimal correlation function, is sufficient to assess critical properties,
circumventing what is usually employed to determine scaling properties of
physically motivated quantities.

While our investigation leads to the conclusion that the average (spin-
resolved) sign displays scaling properties associated to critical behavior, it
is less clear to understand why this happens. One possible interpretation
comes from the fact that $\langle{\cal S}\rangle$ has the physical
interpretation of the exponential of the difference of free energy densities
of the original fermionic problem, and the one used with the weights in the
Monte Carlo sampling taken to be positive-definite [22]. A similar rationale
can be derived in the case of the spin-resolved sign. For example, in a
bipartite lattice without chemical potential (the first model we investigate,
the SU(2) honeycomb Hubbard model), it is then easy to show that weights
associated to each fermionic flavor are related:
$W_{\overline{\sigma}}(\\{x\\})=C_{\\{x\\}}W_{\sigma}(\\{x\\})$, where $C$ is
a ${\\{x\\}}$-dependent positive constant
($=e^{\lambda\sum_{i\tau}x_{i\tau}}$, with ${\rm
cosh}\,\lambda=e^{|U|\Delta\tau/2}$) [24]. Therefore, the average sign of
either of the weights reads

$\displaystyle\langle{\cal S}_{\sigma}\rangle$ $\displaystyle=$
$\displaystyle\frac{\sum_{\\{x\\}}{\rm
sgn}(W_{\sigma}(\\{x\\}))C_{\\{x\\}}[W_{\sigma}(\\{x\\})]^{2}}{\sum_{\\{x\\}}C_{\\{x\\}}[W_{\sigma}(\\{x\\})]^{2}}$
(4) $\displaystyle\equiv$ $\displaystyle\frac{\sum_{\\{x\\}}{\rm
sgn}(\sqrt{\rho(x)})\rho(x)}{\sum_{\\{x\\}}\rho(x)}\equiv\frac{\sum_{\\{x\\}}\rho^{\prime}(x)}{\sum_{\\{x\\}}\rho(x)}\
.$

From here the logic follows the same as the one with a standard sign problem:
One can define this ratio as the ratio of exponentials of corresponding free
energy densities of probability distributions $\rho$ and $\rho^{\prime}$:

$\langle{\cal S}_{\sigma}\rangle=e^{-\beta V(f_{\rho^{\prime}}-f_{\rho})}\ .$
(5)

For the case of the quantum phase transitions we have investigated, we
demonstrated numerically that this quantity satisfies the scaling ansatz in
the vicinity of the quantum critical point, i.e., $\langle{\cal
S}_{\sigma}\rangle=g(uL^{1/\nu},L_{\tau}/L^{z})$. Consequently, the difference
in free energy densities reads:

$\Delta f\equiv f_{\rho^{\prime}}-f_{\rho}=-\frac{1}{L_{\tau}\Delta\tau\cdot
L^{D}}\log[g(uL^{1/\nu},L_{\tau}/L^{z})]\ .$ (6)

From Fig. 2(d), we notice that $g(uL^{1/\nu},L_{\tau}/L^{z})$ goes from 1 to 0
when $u\simeq 0$, consequently $\Delta f$ shows a departure from zero at this
same point. That is, the difference in free energy densities, initially zero
in the non-interacting regime and in the weakly correlated one, turns finite
when approaching the Mott phase, as if the free energy densities of the models
with probability distributions $\rho^{\prime}$ and $\rho$ undergo a
‘transition’ to distinct values. Finally, as we fix the ratio
$L_{\tau}/L^{z}=a$,

$\Delta f=-\frac{1}{a\Delta\tau\cdot L^{D+z}}\log[g(uL^{1/\nu},a)].$ (7)

Figures 5(a) and 5(b) summarize this reasoning for the SU(2) honeycomb Hubbard
model, showing the scaling of the difference of free energy densities, where
we emphasize that $\Delta F\equiv\Delta f\cdot L^{D+z}$ is a function of the
interaction strength 222From a historic perspective on the study of the sign
problem in DQMC it is easy to see how this was missed: in both the one-
dimensional and two-dimensional square lattice Hubbard models, the critical
interaction is $U_{c}/t=0^{+}$, hence at any finite interaction the difference
in free energy densities is positive. The honeycomb Hubbard model (and perhaps
others with $U_{c}>0$) is unique in this aspect in that it allows one to
precisely understand the regime where the free energy densities of the two
systems differentiate (start to diverge), and its relation to the onset of the
ordered phase.. Similar logic applies to the Ionic Hubbard model [Figs. 5(c)
and 5(d)], in spite of the partial weights no longer being trivially related.
That is, within the band-insulating regime, the difference in free energy
$\Delta F$ of the two distributions is zero, deviating from each other once
the the correlated metal phase at $U=U_{c}$ is approached.

Figure 5: Scaling analysis of the spin resolved sign and the corresponding
difference of free energies for the SU(2) honeycomb Hubbard model [(a) and
(b)] and the square lattice SU(2) Ionic Hubbard model [(c) and (d)].
Parameters used are indicated, with imaginary-time discretization
$t\Delta\tau=0.1$.

One of the aspects of the scaling analysis of $\langle{\cal
S}_{\sigma}\rangle$ that defies current expectations relates to the value of
the dynamical critical exponent $z$ we have used. In particular for the SU(2)
honeycomb Hubbard model, field-theory predictions assert $z=1$ [32, 33], and
numerical simulations using projective quantum Monte Carlo methods that
directly tackle the $T=0$ limit often use this as a starting point [16, 17,
19]. Our simulations, on the other hand, employ the corresponding finite-
temperature version of this algorithm [23, 24], such that the ground-state
physics is only obtained asymptotically when $\beta\to\infty$ or when the
typical correlation lengths $\xi$ are sufficiently large such that they are
comparable to the linear system size $L$ [24]. A verification of the latter is
possible by examining the $\beta$-dependence of the antiferromagnetic
structure factor

$S_{\rm
AF}=\frac{1}{2L^{2}}\sum_{i,j}{(-1)}^{\delta}\langle(\hat{n}_{i,\uparrow}-\hat{n}_{i,\downarrow})(\hat{n}_{j,\uparrow}-\hat{n}_{j,\downarrow})\rangle$
(8)

with $\delta=0$ ($\delta=1$) if sites $i$ and $j$ are in the same (different)
sublattice.

Figure 6: Dependence of the antiferromagnetic structure factor on the inverse
temperature $\beta$ (number of imaginary-time slices $L_{\tau}$, fixing
$t\Delta\tau=0.1$) for the SU(2) honeycomb Hubbard model in the vicinity of
the quantum critical point: (a) $U/t=3.60$, (b) $U/t=3.77$ and (c) $U/t=3.85$.
The vertical arrows mark the value of $L_{\tau}$ used in the scaling analysis
of $\langle{\cal S}_{\sigma}\rangle$, color-matching with the corresponding
linear system size $L$.

A saturation of $S_{\rm AF}$ with increasing $\beta$ indicates that $\xi\simeq
L$ [61], and is readily obtained deep in the ordered phase. Close to $U_{c}$,
however, the observation of such a saturation demands (numerically)
prohibitive values of $L_{\tau}$, as indicated in Fig. 6. As a result, the
currently employed values of the imaginary-time slices in the scaling analysis
of $\langle{\cal S}_{\sigma}\rangle$ [marked by the arrows in Fig. 6]
inevitably lead to the conclusion that finite-temperature effects are at play
here, and the observed scaling relates to a low-but-finite temperature
crossover that emanates from the quantum phase transition. Consequently, the
dynamical critical exponent need not be pinned at $z=1$, and the value we use,
obtained after a scaling of $L_{\tau}/L^{z}$ for the current range of
imaginary-time slices employed, endows the ability to study the quantum
criticality. In other words, the fact that we adjust the dynamical critical
exponent for the current range of temperatures is what allows one to obtain
numerically accurate values of the pair $(U_{c},\nu)$.

While this may come as a surprise, it becomes more clear after performing a
scaling of a physical quantity, in particular the one which dictates the onset
of magnetic ordering at the quantum critical point, the antiferromagnetic
structure factor, as shown in Appendix C. There one finds that a $z<1$ (in
practice $z\simeq 0.5$) gives the best data collapse, and that the same
combination of $(U_{c},\nu)$ which scales the spin-resolved sign is seen to
similarly scale $S_{\rm AF}$.

Turning to the thermal transitions, a final, more speculative, line of inquiry
is to investigate the potential existence of a Kosterlitz-Thouless transition
in the repulsive Hubbard model (and its variants) away from half-filling, via
the analysis of the average sign when entering the $\langle{\cal S}\rangle\to
0$ phase in this Hamiltonian [13]. A preliminary study is given in the SM
[36], affirmatively indicating this connection.

In summary, our determination of the quantum criticality via the scaling of
the sign in many models is within the range of existing investigations.
Because the literature shows considerable variation in the quantitative
location of the different QCP’s [18, 19, 41, 42, 50, 51], the results here
offer an alternative, and potentially more accurate, route of resolving a
challenging problem in correlated electron models. Similar investigation can
be carried out with methods that tackle directly the ground-state limit, and
we defer this analysis for a future study.

###### Acknowledgements.

We acknowledge discussions with I. Herbut, S. Sachdev, Z.Y. Meng, and R.
Singh. R.M. acknowledges support from the National Natural Science Foundation
of China (NSFC) Grants No. U1930402, 12050410263, 12111530010 and No.
11974039. R.T.S. was supported by the grant DE‐SC0014671 funded by the U.S.
Department of Energy, Office of Science. Computations were performed on the
Tianhe-2JK at the Beijing Computational Science Research Center.

##  Appendix A Methods

In all calculations we make use of the finite-temperature determinant quantum
Monte Carlo method [23, 24]. Via a sequence of Trotterization, Hubbard-
Stratonovich decoupling of quartic terms by means of the introduction of an
auxiliary field [28], and final fermionic integration, the partition function
is written in terms of the determinants of $N\times N$ matrices $M_{\sigma}$
[where $N$ is the number of sites in Eq. (1)] for each spin-component
$\sigma$. These are the weights $W(\\{x\\})$ referred to in the main text.
Instead of summing over all configurations of the field
$\\{x\\}\to\\{x_{i\tau}\\}$, importance sampling is performed while observing
the statistical convergence of both physical observables (when possible) and
the average sign of the weights. The only approximation used is the imaginary-
time discretization $\Delta\tau$ which we take as $1/10$ for the quantum
transitions or $1/16$ for the thermal ones. Statistical sampling varies among
the different models, but in all cases an average of the results are taken for
each individual set of parameters over dozens of independent samplings
(typically from 20 to 48), with thousands of Monte Carlo sweeps for each run.

In order to decouple the interactions in all SU(2) models we investigate, we
apply either the spin-decomposed Hubbard-Stratonovich transformation [28, 4],

$\displaystyle e^{-\Delta\tau
U(\hat{n}_{i\uparrow}-\frac{1}{2})(\hat{n}_{i\downarrow}-\frac{1}{2})}=\frac{1}{2}e^{-U\Delta\tau/4}\sum_{x_{i}=\pm
1}e^{\lambda x_{i}(\hat{n}_{i\uparrow}-\hat{n}_{i\downarrow})}\,,$

for repulsive interactions ($U>0$), or its counterpart (charge decomposition)

$\displaystyle e^{-\Delta\tau
U(\hat{n}_{i\uparrow}-\frac{1}{2})(\hat{n}_{i\downarrow}-\frac{1}{2})}=\frac{1}{2}e^{-|U|\Delta\tau/4}\sum_{x_{i}=\pm
1}e^{\lambda x_{i}(\hat{n}_{i\uparrow}+\hat{n}_{i\downarrow}-1)}\,,$

in the case that $U<0$. In both transformations, ${\rm
cosh}\,\lambda=e^{|U|\Delta\tau/2}$. Finally, the matrices $M_{\sigma}$
entering in the weights read

$\displaystyle
M_{\sigma}=\mathds{1}+B_{\sigma,L_{\tau}}B_{\sigma,L_{\tau}-1}\ldots
B_{\sigma,1},$ (9)

with $B_{\sigma,\tau}=e^{K}e^{V_{\sigma,\tau}}$. Here, $K$ is an imaginary-
time independent $N\times N$ matrix containing all one-body terms in the
Hamiltonian (including hopping and chemical potential), whose entries are
multiplied by $-\Delta\tau$. In turn, $V_{\sigma,\tau}$ is diagonal with
entries that depend on the Hubbard-Stratonovich transformation used. For the
repulsive case, $V^{ii}_{\sigma,\tau}=\lambda\sigma x_{i\tau}$ ($\sigma=\pm 1$
for $\uparrow$ and $\downarrow$), while
$V^{ii}_{\uparrow,\tau}=V^{ii}_{\downarrow,\tau}=\lambda x_{i\tau}$ for
attractive interactions.

##  Appendix B Eigenvalues of the $M_{\sigma}$ matrices

A possibility to infer numerically that the signs of the determinants can
track phase transitions is via the analysis of the spectrum of the
$M_{\sigma}$ matrices, as defined in Eq. 9, whose determinant gives the
partial weight of a certain configuration $\\{x\\}$. Similarly, one can define
this matrix in its space-time formulation [23],

$M_{\sigma}(\\{x\\})=\begin{bmatrix}\mathds{1}&&&&B_{\sigma,{L_{\tau}}}\\\
-B_{\sigma,1}&\mathds{1}&&&\\\ &-B_{\sigma,2}&\mathds{1}&&\\\
&&\ddots&\ddots&\\\ &&&-B_{\sigma,{L_{\tau}-1}}&\mathds{1}\end{bmatrix},$ (10)

by structuring the $N\times N$ matrices $B_{\sigma,\tau}$, the single-particle
propagators, as defined in the Appendix A. A drawback is that
$M_{\sigma}(\\{x\\})$ has now dimensions $(NL_{\tau})\times(NL_{\tau})$, but
one of the benefits of this representation is that the range of eigenvalues is
now shrunken, while preserving the value of the determinant. Besides that, it
is numerically stable since no matrix multiplications among $B_{\sigma}$’s are
necessary to build it.

Figure 7: The eigenspectrum of the $M_{\sigma}(\\{x\\})$ matrix represented in
the complex plane for a single Hubbard-Stratonovich configuration
$\\{x_{i,\tau}\\}$ extracted in the Monte Carlo sampling, for $U/t=2$ (a) and
$U/t=6$ (b). (c) and (d) exhibit a two-dimensional histogram of the
eigenspectrum when combining results of 96 configurations – brighter colors
display a larger counting. In (a–d), the linear lattice size is $L=9$ and
$L_{\tau}=200$, such that each configuration leads to $32\ 400$ eigenvalues.
(e) shows both the (normalized) number of eigenvalues which are in the
negative real axis (these dictate whether there is a spin-resolved sign
problem or not) and the fraction of the configurations that possess an odd
number of eigenvalues in $\mathbb{R}^{-}$. Empty (full) symbols refer to $L=6$
($L=9$); vertical shaded region gives the confidence interval of the quantum
critical point location. As in the main text, here we use $t\,\Delta\tau=0.1$.

Focusing on the SU(2) honeycomb Hubbard model, we start by analyzing in Fig. 7
the spectrum $\\{\varepsilon_{i}\\}$ of $M_{\sigma}(\\{x\\})$ for values of
the interactions far below ($U/t=2$) and far above ($U/t=6$) the transition
point $U_{c}$. Figures 7 (a) and 7 (b) show that a structural transition
occurs in the eigenvalue spectrum, here computed for a single typical
configuration of the auxiliary field $\\{x_{i,\tau}\\}$. This observed
structural transition is generic, as shown by the corresponding two-
dimensional histograms in Figs. 7 (c) and (d), obtained by combining
eigenvalues of four field configurations ‘visited’ over the course of the
Monte Carlo sampling for 24 independently seeded Markov chains, resulting in
96 configurations in total.

While the quantum phase transition is hinted in the eigenvalues of
$M_{\sigma}(\\{x\\})$, so far we have not drawn a connection to the spin-
resolved sign problem. Being a real matrix (for this model with the Hubbard
Stratonovich transformation highlighted in the Appendix A), its eigenvalues
come as either complex conjugate pairs, or real numbers. Since the determinant
(product of eigenvalues) does not change sign when multiplying the conjugate
pairs, a sign problem is only a function of the number of eigenvalues in the
negative real axis, $n_{\mathbb{R}^{-}}$. That is, if $n_{\mathbb{R}^{-}}$ is
odd (even), $\det M_{\sigma}(\\{x\\})<0$ ($>0$). Typically,
$n_{\mathbb{R}^{-}}$ is very small in comparison to the total number of
eigenvalues $2L^{2}L_{\tau}$ in this Hamiltonian. Yet, it is a clear function
of the interaction strength, growing at large $U/t$, as shown in Fig. 7(e) by
$\overline{n_{\mathbb{R}^{-}}}$, after averaging over different configurations
$\\{x\\}$. Finally, the percentage of those configurations that possess an odd
number of eigenvalues in the negative real axis, $P_{\rm
odd}(n_{\mathbb{R}^{-}})$, also grows reaching around 50$\%$ within the
ordered phase. As a result, the average spin-resolved sign $\langle{\cal
S}_{\sigma}\rangle$ converges to zero.

It is currently unclear to us if a physical meaning can be attributed to the
number of negative eigenvalues of $M_{\sigma}(\\{x\\})$, in terms of the
fields $\\{x_{i,\tau}\\}$, and the winding of world lines they are associated
with.

##  Appendix C The scaling of $S_{\rm AF}$

Figure 8: The scaling analysis of the antiferromagnetic structure factor. In
(a) and (b) the scaled $S_{\rm AF}$ dependence on $L_{\tau}/L^{z}$, with
$z=0.5$ and 1, respectively. Insets give the corresponding cost function
dependence on the ratio of exponents $\beta/\nu$. (c) Color-mesh plot of the
cost function ${\cal C}(z,\beta/\nu)$, where the white marker at
$(z^{*},\beta/\nu^{*})=(0.54,0.885)$ depicts its minimum value. (d) The
scaling collapse of $S_{\rm AF}$ with $uL^{1/\nu}$, with parameters extracted
from the analysis of the spin-resolved sign in the main text [Fig. 2]. The
imaginary-time discretization used is $t\Delta\tau=0.1$.

In the main text, we argue that owing to finite-temperature effects one needs
to adjust the dynamical critical exponent $z$ from its expected $z=1$ value in
order to perform a scaling of the spin-resolved sign $\langle{\cal
S}_{\sigma}\rangle$. While the scaling analysis we perform for $\langle{\cal
S}_{\sigma}\rangle$ is quantitatively precise despite its novelty, similar
constraints should apply to the case of the scaling of physical quantities.
Following this logic we perform the scaling of the antiferromagnetic order
parameter, $m_{\rm AF}=\lim_{L\to\infty}\sqrt{S_{\rm AF}/N}$, with $N=2L^{2}$
the number of sites of the SU(2) honeycomb Hubbard model, Eq. (1). This
quantity follows a scaling ansatz of the form $m_{\rm
AF}=L^{-\beta/\nu}g(uL^{1/\nu},L_{\tau}/L^{z})$ [17, 19], which in turn
implies the antiferromagnetic structure factor scales as

$\frac{S_{\rm AF}}{N}=L^{-2\beta/\nu}g(uL^{1/\nu},L_{\tau}/L^{z})$ (11)

Unlike previous studies that tackle this same transition using $T=0$ quantum
Monte Carlo methods [16, 17, 18, 19], we make the $L_{\tau}/L^{z}$ dependence
explicit in order to account for a finite-$T$ influence on the scaling.

We start by fixing $U/t=3.77$, such that the dependence on the first argument
of the scaling function in Eq. (11) is negligible, i.e., $u\simeq 0$ 333Notice
that this value is chosen based on the critical value $U_{c}=3.765$ that comes
from the scaling of the spin-resolved sign.. Figure 8(a) and 8(b) display the
scaled structure factor vs. $L_{\tau}/L^{z}$ by fixing $z=0.5$ and $1.0$,
respectively, while adjusting the ratio of exponents $\beta/\nu$ that gives
the best data collapse. The latter is obtained by the minimization of the cost
function ${\cal C}(\beta/\nu)$, whose definition is the same as given in the
main text, and is displayed as insets in Figs. 8(a) and 8(b). Notably, the
data collapse is significantly better at $z\simeq 0.5$ in comparison to the
one for $z=1$, a first indication that a dynamic critical exponent smaller
than one results in an improved scaling. Compiling the cost function in a
range of $(z,\beta/\nu)$-values, shown in Fig. 8(c) as a color-mesh plot, we
obtain the minimum cost function at $z^{*}=0.54$ and $(\beta/\nu)^{*}=0.885$.
The latter is compatible with the value $\beta/\nu\simeq 0.9$ obtained after
the first-order $\epsilon$-expansion of the Gross-Neveu model [33, 17].
Lastly, by fixing this ratio of exponents $(\beta/\nu)$ while choosing the set
of parameters used to perform the scaling of $\langle{\cal S}_{\sigma}\rangle$
in the main text, $(U_{c}/t,\nu,z)=(3.765,0.84,0.5)$, we report in Fig. 8(e)
the dependence of the scaled structure factor with respect to the first
argument of the scaling function: The exhibited collapse is remarkably good,
thus confirming that the average (spin-resolved) sign can indeed be used to
infer criticality. As far as we are aware, this is the first time that a
finite-temperature quantum Monte Carlo method was used to obtain critical
exponents of a transition pertaining to the Gross-Neveu universality class,
and the key step for its success is the tuning of the dynamical critical
exponent $z$.

##  Appendix D Other Hubbard-Stratonovich transformations

Our main results point out that the (spin-resolved) average sign carries
fundamental information about phase transitions and their universality
classes. However, given that the sign problem is basis-dependent, that is, by
choosing another Hubbard-Stratonovich transformation, the average sign of the
quantum Monte Carlo weights can change [63, 64], an immediate question that
arises is: can one still infer quantum critical points using ${\rm
sgn}(|W_{\sigma}(\\{x\\})|)$? To answer it, we report further numerical tests.

Figure 9: Dependence of the (one minus the) average spin resolved phase on the
interaction magnitude $U/t$ for the SU(2) honeycomb Hubbard model. We fix the
ratio $L_{\tau}/L^{z}$ approximately constant, with $z=0.5$. The inset is a
zoom-in of the critical region, and the vertical shaded region are the recent
estimations in the literature of the QCP location [18, 19]. The imaginary-time
discretization used is $t\Delta\tau=0.1$.

An often used Hubbard-Stratonovich transformation is one that explicitly
conserves the SU(2) symmetry for each configuration $\\{x_{i\tau}\\}$ of the
field [54, 53, 52]:

$\displaystyle e^{-\Delta\tau
U(\hat{n}_{i\uparrow}+\hat{n}_{i\downarrow}-1)^{2}/2}=$
$\displaystyle\sum_{x_{i\tau}=\pm 1,\pm
2}\gamma(x_{i\tau})\prod_{\sigma}e^{{\rm i}\sqrt{\Delta\tau
U/2}\eta(x_{i\tau})(\hat{n}_{i\sigma}-1/2)}+{\cal O}(\Delta\tau^{4})\ .$ (12)

It comes at the expense of having a four-valued discrete field $x_{i\tau}=\pm
1,\pm 2$, accompanied by a few (real) constants,

$\displaystyle\gamma(\pm 1)=1+\sqrt{6}/3\ ;\ \eta(\pm
1)=\pm\sqrt{2(3-\sqrt{6})}$ $\displaystyle\gamma(\pm 2)=1-\sqrt{6}/3\ ;\
\eta(\pm 2)=\pm\sqrt{2(3+\sqrt{6})}\ .$ (13)

In such case, the transformation is not strictly exact but brings an error
proportional to ${\cal O}(\Delta\tau^{4})$, negligible in practice in
comparison to the one introduced by the Trotter decomposition of the one and
two-body terms in the Hamiltonian, ${\cal O}(\Delta\tau^{2})$.

Given its form, the Monte Carlo weights can become complex, in principle, but
in the SU(2) honeycomb Hubbard model, owing to its bipartite nature, one can
show that the weights associated to the two-spin components are complex
conjugate pairs at half-filling [54, 53, 52]. As a result, no phase problem
emerges. Nonetheless, this does not guarantee that the average phase of each
fermionic flavor $\langle e^{{\rm i}\theta_{\sigma}}\rangle$ needs to be real
and raises the immediate question whether the sign still captures information
about the onset of an ordered phase. Numerical simulations we performed,
however, point out affirmatively to both: ${\rm Im}\langle e^{{\rm
i}\theta_{\sigma}}\rangle\rightarrow 0$ over the course of the sampling, and
that $\langle e^{{\rm i}\theta_{\sigma}}\rangle(U=U_{c})\rightarrow 0$ when
approaching the thermodynamic limit. The latter is reported in Fig. 9, using
$L_{\tau}/L^{1/2}=\frac{240}{18^{1/2}}$ approximately fixed for different
system sizes.

Unlike in the case of a spin-decomposed Hubbard-Stratonovich transformation,
where a crossing of $\langle{\cal S}_{\sigma}\rangle$ for different system
sizes leads to an immediate identification of $U_{c}$, here the nature of the
dependence of the average spin-resolved phase with $U$ makes a scaling process
more challenging. However, the trend of $\langle e^{{\rm
i}\theta_{\sigma}}\rangle$ with different system sizes clearly suggests that
the average phase converges towards zero when approaching the quantum critical
point. Future investigations with both larger sizes and improved statistics
may settle the possible determination of critical exponents in this regime,
but it is overall clear that the (partial) weights remain efficient trackers
of a quantum phase transition despite the use of a fundamentally different
Hubbard-Stratonovich transformation.

## References

  * Poilblanc [1991] D. Poilblanc, Twisted boundary conditions in cluster calculations of the optical conductivity in two-dimensional lattice models, Phys. Rev. B 44, 9562 (1991).
  * Lin  _et al._ [1993] H. Lin, J. Gubernatis, H. Gould, and J. Tobochnik, Exact diagonalization methods for quantum systems, Computers in Physics 7, 400 (1993).
  * Schollwöck [2005] U. Schollwöck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).
  * Loh and Gubernatis [1992] E. Loh and J. Gubernatis, Chapter 4 - Stable numerical simulations of models of interacting electrons in Condensed-Matter Physics, in  _Electronic Phase Transitions_, Modern Problems in Condensed Matter Sciences, Vol. 32, edited by W. Hanke and Y. Kopaev (Elsevier, 1992) pp. 177–235.
  * Gubernatis  _et al._ [2016] J. Gubernatis, N. Kawashima, and P. Werner,  _Quantum Monte Carlo Methods: Algorithms for Lattice Models_ (Cambridge University Press, 2016).
  * Dagotto [1994] E. Dagotto, Correlated electrons in high-temperature superconductors, Rev. Mod. Phys. 66, 763 (1994).
  * Qin  _et al._ [2021] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, The Hubbard model: A computational perspective (2021), arXiv:2104.00064 [cond-mat.str-el] .
  * Marvian  _et al._ [2019] M. Marvian, D. A. Lidar, and I. Hen, On the computational complexity of curing non-stoquastic Hamiltonians, Nat. Comm. 10, 1571 (2019).
  * Hangleiter  _et al._ [2020] D. Hangleiter, I. Roth, D. Nagaj, and J. Eisert, Easing the Monte Carlo sign problem, Science Advances 6, 10.1126/sciadv.abb8341 (2020).
  * Levy and Clark [2021] R. Levy and B. K. Clark, Mitigating the sign problem through basis rotations, Phys. Rev. Lett. 126, 216401 (2021).
  * Wan  _et al._ [2020] Z.-Q. Wan, S.-X. Zhang, and H. Yao, Mitigating sign problem by automatic differentiation (2020), arXiv:2010.01141 [cond-mat.str-el] .
  * Troyer and Wiese [2005] M. Troyer and U.-J. Wiese, Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations, Phys. Rev. Lett. 94, 170201 (2005).
  * Mondaini  _et al._ [2022] R. Mondaini, S. Tarat, and R. T. Scalettar, Quantum critical points and the sign problem, Science 375, 418 (2022).
  * Wessel  _et al._ [2017] S. Wessel, B. Normand, F. Mila, and A. Honecker, Efficient Quantum Monte Carlo simulations of highly frustrated magnets: the frustrated spin-1/2 ladder, SciPost Phys. 3, 005 (2017).
  * Gross and Neveu [1974] D. J. Gross and A. Neveu, Dynamical symmetry breaking in asymptotically free field theories, Phys. Rev. D 10, 3235 (1974).
  * Sorella  _et al._ [2012] S. Sorella, Y. Otsuka, and S. Yunoki, Absence of a spin liquid phase in the Hubbard model on the honeycomb lattice, Sci. Rep. 2, 992 (2012).
  * Assaad and Herbut [2013] F. F. Assaad and I. F. Herbut, Pinning the order: The nature of quantum criticality in the Hubbard model on honeycomb lattice, Phys. Rev. X 3, 031010 (2013).
  * Parisen Toldin  _et al._ [2015] F. Parisen Toldin, M. Hohenadler, F. F. Assaad, and I. F. Herbut, Fermionic quantum criticality in honeycomb and $\pi$-flux Hubbard models: Finite-size scaling of renormalization-group-invariant observables from quantum Monte Carlo, Phys. Rev. B 91, 165108 (2015).
  * Otsuka  _et al._ [2016] Y. Otsuka, S. Yunoki, and S. Sorella, Universal quantum criticality in the metal-insulator transition of two-dimensional interacting Dirac electrons, Phys. Rev. X 6, 011029 (2016).
  * Zhang  _et al._ [2019] Y.-X. Zhang, W.-T. Chiu, N. C. Costa, G. G. Batrouni, and R. T. Scalettar, Charge order in the Holstein model on a honeycomb lattice, Phys. Rev. Lett. 122, 077602 (2019).
  * Otsuka  _et al._ [2020] Y. Otsuka, K. Seki, S. Sorella, and S. Yunoki, Dirac electrons in the square-lattice Hubbard model with a $d$-wave pairing field: The chiral Heisenberg universality class revisited, Phys. Rev. B 102, 235105 (2020).
  * Loh  _et al._ [1990] E. Loh, J. Gubernatis, R. Scalettar, S. White, D. Scalapino, and R. Sugar, Sign problem in the numerical simulation of many-electron systems, Phys. Rev. B 41, 9301 (1990).
  * Blankenbecler  _et al._ [1981] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Monte Carlo calculations of coupled boson-fermion systems. I, Phys. Rev. D 24, 2278 (1981).
  * Hirsch [1985] J. E. Hirsch, Two-dimensional Hubbard model: Numerical simulation study, Phys. Rev. B 31, 4403 (1985).
  * Iglovikov  _et al._ [2015] V. Iglovikov, E. Khatami, and R. Scalettar, Geometry dependence of the sign problem in quantum Monte Carlo simulations, Phys. Rev. B 92, 045110 (2015).
  * Wu and Zhang [2005] C. Wu and S.-C. Zhang, Sufficient condition for absence of the sign problem in the fermionic quantum Monte Carlo algorithm, Phys. Rev. B 71, 155115 (2005).
  * Wei  _et al._ [2016] Z. C. Wei, C. Wu, Y. Li, S. Zhang, and T. Xiang, Majorana positivity and the fermion sign problem of quantum Monte Carlo simulations, Phys. Rev. Lett. 116, 250601 (2016).
  * Hirsch [1983] J. E. Hirsch, Discrete Hubbard-Stratonovich transformation for fermion lattice models, Phys. Rev. B 28, 4059 (1983).
  * Binder and Wang [1989] K. Binder and J.-S. Wang, Finite-size effects at critical points with anisotropic correlations: Phenomenological scaling theory and Monte Carlo simulations, Journal of Statistical Physics 55, 87 (1989).
  * Rieger and Young [1994] H. Rieger and A. P. Young, Zero-temperature quantum phase transition of a two-dimensional Ising spin glass, Phys. Rev. Lett. 72, 4141 (1994).
  * Guo  _et al._ [1994] M. Guo, R. N. Bhatt, and D. A. Huse, Quantum critical behavior of a three-dimensional Ising spin glass in a transverse magnetic field, Phys. Rev. Lett. 72, 4137 (1994).
  * Herbut  _et al._ [2009a] I. F. Herbut, V. Juričić, and B. Roy, Theory of interacting electrons on the honeycomb lattice, Phys. Rev. B 79, 085116 (2009a).
  * Herbut  _et al._ [2009b] I. F. Herbut, V. Juričić, and O. Vafek, Relativistic Mott criticality in graphene, Phys. Rev. B 80, 075432 (2009b).
  * Šuntajs  _et al._ [2020] J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, Ergodicity breaking transition in finite disordered spin chains, Phys. Rev. B 102, 064207 (2020).
  * Aramthottil  _et al._ [2021] A. S. Aramthottil, T. Chanda, P. Sierant, and J. Zakrzewski, Finite-size scaling analysis of the many-body localization transition in quasiperiodic spin chains, Phys. Rev. B 104, 214201 (2021).
  * [36] See Supplemental Material, and references therein, for more details on the scaling of the total sign, different scaling procedures, and further investigation on the thermal transitions either in attractive or repulsive Hubbard model.
  * Assaad [2002] F. F. Assaad, Quantum Monte Carlo methods on lattices: The determinantal approach, in  _Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms, Lecture Notes NIC Series_ , Vol. 10 (2002) pp. 99–156.
  * Fabrizio  _et al._ [1999a] M. Fabrizio, A. O. Gogolin, and A. A. Nersesyan, From band insulator to Mott insulator in one dimension, Phys. Rev. Lett. 83, 2014 (1999a).
  * Kampf  _et al._ [2003] A. Kampf, M. Sekania, G. Japaridze, and P. Brune, Nature of the insulating phases in the half-filled ionic Hubbard model, Journal of Physics: Condensed Matter 15, 5895 (2003).
  * Garg  _et al._ [2006] A. Garg, H. Krishnamurthy, and M. Randeria, Can correlations drive a band insulator metallic?, Phys. Rev. Lett. 97, 046403 (2006).
  * Paris  _et al._ [2007] N. Paris, K. Bouadim, F. Hebert, G. Batrouni, and R. Scalettar, Quantum Monte Carlo study of an interaction-driven band-insulator–to–metal transition, Phys. Rev. Lett. 98, 046403 (2007).
  * Bouadim  _et al._ [2007] K. Bouadim, N. Paris, F. Hébert, G. Batrouni, and R. Scalettar, Metallic phase in the two-dimensional ionic Hubbard model, Phys. Rev. B 76, 085112 (2007).
  * Craco  _et al._ [2008] L. Craco, P. Lombardo, R. Hayn, G. Japaridze, and E. Müller-Hartmann, Electronic phase transitions in the half-filled ionic Hubbard model, Phys. Rev. B 78, 075121 (2008).
  * Garg  _et al._ [2014] A. Garg, H. Krishnamurthy, and M. Randeria, Doping a correlated band insulator: A new route to half-metallic behavior, Phys. Rev. Lett. 112, 106406 (2014).
  * Lou  _et al._ [2003] J. Lou, S. Qin, T. Xiang, C. Chen, G.-S. Tian, and Z. Su, Transition from band insulator to Mott insulator in one dimension: Critical behavior and phase diagram, Phys. Rev. B 68, 045110 (2003).
  * Manmana  _et al._ [2004] S. R. Manmana, V. Meden, R. M. Noack, and K. Schönhammer, Quantum critical behavior of the one-dimensional ionic Hubbard model, Phys. Rev. B 70, 155115 (2004).
  * Tincani  _et al._ [2009] L. Tincani, R. M. Noack, and D. Baeriswyl, Critical properties of the band-insulator-to-Mott-insulator transition in the strong-coupling limit of the ionic Hubbard model, Phys. Rev. B 79, 165109 (2009).
  * Fabrizio  _et al._ [1999b] M. Fabrizio, A. O. Gogolin, and A. A. Nersesyan, From band insulator to Mott insulator in one dimension, Phys. Rev. Lett. 83, 2014 (1999b).
  * Scalettar  _et al._ [1989] R. T. Scalettar, E. Y. Loh, J. E. Gubernatis, A. Moreo, S. R. White, D. J. Scalapino, R. L. Sugar, and E. Dagotto, Phase diagram of the two-dimensional negative-U Hubbard model, Phys. Rev. Lett. 62, 1407 (1989).
  * Moreo and Scalapino [1991] A. Moreo and D. J. Scalapino, Two-dimensional negative-U Hubbard model, Phys. Rev. Lett. 66, 946 (1991).
  * Paiva  _et al._ [2004] T. Paiva, R. R. dos Santos, R. T. Scalettar, and P. J. H. Denteneer, Critical temperature for the two-dimensional attractive Hubbard model, Phys. Rev. B 69, 184501 (2004).
  * Hohenadler  _et al._ [2012] M. Hohenadler, Z. Y. Meng, T. C. Lang, S. Wessel, A. Muramatsu, and F. F. Assaad, Quantum phase transitions in the Kane-Mele-Hubbard model, Phys. Rev. B 85, 115132 (2012).
  * Zheng  _et al._ [2011] D. Zheng, G.-M. Zhang, and C. Wu, Particle-hole symmetry and interaction effects in the Kane-Mele-Hubbard model, Phys. Rev. B 84, 205121 (2011).
  * Meng  _et al._ [2010] Z. Meng, S. Wessel, A. Muramatsu, T. Lang, and F. Assaad, Quantum spin liquid emerging in two-dimensional correlated Dirac fermions, Nature 464, 847 (2010).
  * Wang  _et al._ [2014] L. Wang, P. Corboz, and M. Troyer, Fermionic quantum critical point of spinless fermions on a honeycomb lattice, New J. Phys. 16, 103008 (2014).
  * Li  _et al._ [2015] Z.-X. Li, Y.-F. Jiang, and H. Yao, Solving the fermion sign problem in quantum Monte Carlo simulations by Majorana representation, Phys. Rev. B 91, 241117 (2015).
  * Li  _et al._ [2016] Z.-X. Li, Y.-F. Jiang, and H. Yao, Majorana-time-reversal symmetries: A fundamental principle for sign-problem-free quantum Monte Carlo simulations, Phys. Rev. Lett. 117, 267002 (2016).
  * Note [1] The original total weight given by a single determinant in the single-particle basis is converted to a product of weights related by complex conjugation when using the Majorana basis to describe the Hamiltonian, hence no sign (phase) problem.
  * Götz  _et al._ [2022] A. Götz, S. Beyl, M. Hohenadler, and F. F. Assaad, Valence-bond solid to antiferromagnet transition in the two-dimensional Su-Schrieffer-Heeger model by Langevin dynamics, Phys. Rev. B 105, 085151 (2022).
  * Note [2] From a historic perspective on the study of the sign problem in DQMC it is easy to see how this was missed: in both the one-dimensional and two-dimensional square lattice Hubbard models, the critical interaction is $U_{c}/t=0^{+}$, hence at any finite interaction the difference in free energy densities is positive. The honeycomb Hubbard model (and perhaps others with $U_{c}>0$) is unique in this aspect in that it allows one to precisely understand the regime where the free energy densities of the two systems differentiate (start to diverge), and its relation to the onset of the ordered phase.
  * Paiva  _et al._ [2005] T. Paiva, R. Scalettar, W. Zheng, R. Singh, and J. Oitmaa, Ground-state and finite-temperature signatures of quantum phase transitions in the half-filled Hubbard model on a honeycomb lattice, Phys. Rev. B 72, 085123 (2005).
  * Note [3] Notice that this value is chosen based on the critical value $U_{c}=3.765$ that comes from the scaling of the spin-resolved sign.
  * Hirsch [1986] J. E. Hirsch, Connection between world-line and determinantal functional-integral formulations of the Hubbard model, Phys. Rev. B 34, 3216 (1986).
  * Batrouni and Scalettar [1990] G. G. Batrouni and R. T. Scalettar, Anomalous decouplings and the fermion sign problem, Phys. Rev. B 42, 2282 (1990).
  * Dalla Piazza  _et al._ [2012] B. Dalla Piazza, M. Mourigal, M. Guarise, H. Berger, T. Schmitt, K. J. Zhou, M. Grioni, and H. M. Rønnow, Unified one-band Hubbard model for magnetic and electronic spectra of the parent compounds of cuprate superconductors, Phys. Rev. B 85, 100508 (2012).
  * Hirayama  _et al._ [2018] M. Hirayama, Y. Yamaji, T. Misawa, and M. Imada, Ab initio effective Hamiltonians for cuprate superconductors, Phys. Rev. B 98, 134501 (2018).
  * Hirayama  _et al._ [2019] M. Hirayama, T. Misawa, T. Ohgoe, Y. Yamaji, and M. Imada, Effective Hamiltonian for cuprate superconductors derived from multiscale ab initio scheme with level renormalization, Phys. Rev. B 99, 245155 (2019).
  * Arovas  _et al._ [2021] D. P. Arovas, E. Berg, S. Kivelson, and S. Raghu, The Hubbard model (2021), arXiv:2103.12097 [cond-mat.str-el] .

Supplementary Materials:
Universality and Critical Exponents of the Fermion Sign Problem

In these Supplementary Materials we provide additional data to support the
scalings used, and try different scaling procedures as well. In addition, we
show the scaling of the total sign $\langle{\cal S}\rangle$ for the Ionic
Hubbard model to argue that this quantity follows $\langle{\cal
S}_{\sigma}\rangle$ as long as it is not protected by some symmetry of the
problem. Lastly, further details of the attractive Hubbard model are
presented, discussing the problematic finite-size effects that occur in
relatively small lattices and the temperature-dependence of the correlation
exponent $\eta$. We end by providing preliminary data suggestive of a KT-
scaling form on both $\langle{\cal S}_{\sigma}\rangle$ and $\langle{\cal
S}\rangle$ for the repulsive Hubbard model at densities below half-filling.

###### Contents

  1. I Introduction
  2. II The SU(2) Honeycomb Hubbard model
  3. III The sign problem
  4. IV The spin resolved sign
  5. V The SU(2) Ionic Hubbard model
  6. VI The attractive Hubbard model
  7. VII Discussion and Outlook
  8. A Methods
  9. B Eigenvalues of the $M_{\sigma}$ matrices
  10. C The scaling of $S_{\rm AF}$
  11. D Other Hubbard-Stratonovich transformations
  12. E Scaling of the total sign in the Ionic Hubbard model
  13. F Using a different scaling procedure
  14. G More details on the attractive Hubbard model
  15. H The repulsive Hubbard model away from half-filling
  16. I Further remarks

##  Appendix E Scaling of the total sign in the Ionic Hubbard model

Figure S1: Scaling analysis of the total sign $\langle{\cal S}\rangle$ of the
QMC weights in the SU(2) Ionic Hubbard model. We use the same parameters
obtained for the scaling of $\langle{\cal S}_{\sigma}\rangle$ in the main
text, $U_{c}/t=2.05$ and $\nu=0.97$. As before, the $L_{\tau}/L^{0.5}$ ratio
is fixed at 50, and $\Delta/t=0.5$.

Since we focus attention on the scaling of the sign of an individual
determinant in the main text, here we show that the same scaling properties
are valid in the case of the total sign of the QMC weights. This is possible
in the case of the Ionic Hubbard model precisely because there is no symmetry-
protection that guarantees the positive-definiteness of the total sign. We
report in Fig. S1 (a) the dependence of the average sign of the product of
determinants (total weight of quantum configurations in QMC calculations) as a
function of the interaction’s strength. Similar to the case of spin-resolved
sign, $\langle{\cal S}\rangle\to 0$ in the ordered phase, displaying a
crossing at the critical interaction $U_{c}/t$ when considering different
system sizes. That the critical interactions and correlation exponent $\nu$
extracted from the scaling of $\langle{\cal S}_{\sigma}\rangle$ similarly
scales the curves for the total sign demonstrates the generality of the
interpretation that the sign of the weights already hosts information about
the critical properties of quantum models.

##  Appendix F Using a different scaling procedure

Figure S2: Different scaling procedure. Cost function for the scaling of the
spin resolved sign in the vicinity of the quantum phase transition for the
Honeycomb Hubbard model (left) and the Ionic Hubbard model on the square
lattice (right). The critical set of parameters $(\nu,U_{c})$ is given by the
white marker and corresponds to the point where ${\cal C}_{\langle{\cal
S}_{\sigma}\rangle}$ is minimized. Cross markers in the left panel, as in Fig.
2(c), display the results of Refs. 18, 19. The results are
$(\nu,U_{c}/t)=(0.71,3.73)$ for the honeycomb case, and
$(\nu,U_{c}/t)=(0.83,2.06)$ in the ionic model.

In computing the critical exponent $\nu$ and the critical interactions $U_{c}$
that characterize and drive the Dirac semi-metal to Mott insulator transition,
we used an analysis that relies on the fitting quality with a high-order
polynomial in the vicinity of the transition. Here, we demonstrate that such
results are hardly changed if we use different approaches. In particular, we
try the cost-function method used in the extraction of the dynamical critical
exponents $z$, and of the critical temperature in the KT-transition for the
attractive Hubbard model. The benefits of this procedure arise from the fact
that ‘flat’ regions in the data set do not contribute to the goodness of the
fit, unlike in a polynomial approach, which may be challenging when dealing
with both flat and decaying regimes of the data, which often incurs in
problems of over-fitting it.

Figure S2 displays the cost function ${\cal C}_{\langle{\cal
S}_{\sigma}\rangle}=\sum_{j}(|y_{j+1}-y_{j}|)/(\max\\{y_{j}\\}-\min\\{y_{j}\\})-1$,
where $y_{j}$ are the values of $\langle{\cal S}_{\sigma}\rangle(L_{\tau},L)$,
ordered according to their $[(U-U_{c})/U_{c}]L^{1/\nu}$ value [34]. This is
equivalent to Fig. 2(e) and 3(e), but using a different collapse method. In
both cases the critical exponents are fairly close to their previous
predictions (see caption of Fig. S2). Small discrepancies highlight the
sensitivity of scaling analyses to the data treatment method. Yet, for the
honeycomb Hubbard model where accurate predictions using physical observables
are known [18, 19], a different scaling method produces similarly good
results.

##  Appendix G More details on the attractive Hubbard model

Figure S3: Severe finite-size effects on the SU(2) attractive Hubbard model.
Spin resolved sign and pair structure factor in the Hubbard model for small
lattice sizes $L=4$ (a, b) and 8 (c, d) in the space of parameters $T/t$ vs.
$\mu/t$. The dome structure in which $\langle{\cal S}_{\sigma}\rangle\to 0$
and $P_{s}$ is large tracks the regimes where a $-\vec{k},\vec{k}$ ‘nesting’
of the Fermi surfaces occurs. Here interaction strength is $U/t=-4$.

### Severe finite-size effects at small lattice sizes.—

The attractive Hubbard model (analyzed in detail in the main text concerning
the relationship between the spin-resolved average sign and superconducting
behavior) has the peculiarity of being especially prone to the manifestation
of finite-size effects in comparison to the other models investigated. In Fig.
S3, we show the dependence of $\langle{\cal S}_{\sigma}\rangle$ and $P_{s}$ in
the $T/t$ vs. $\mu/t$ plane for smaller lattice sizes than the ones considered
for the KT scaling analysis. It is clear now that only a few special values of
$\mu/t$ lead to an enhancement of the $s$-wave pairing at sufficiently low
temperatures. Nevertheless, the structures are also reproduced by the onset of
regions where $\langle{\cal S}_{\sigma}\rangle\to 0$. One can understand this
non-monotonic behavior by recalling that in small lattices, the coarse grained
discretization of the Brillouin zone implies that a few special $\mu/t$ lead
to a perfect nesting condition, where the Fermi surface crosses one of the
allowed momentum points. Remarkably this occurs even at the intermediate
values of $U/t=4$ as seen in Fig. S3. Such oscillations have been noted before
[50], and using small lattice sizes to extract the critical temperature can be
responsible for large deviations from the expected behavior in the
thermodynamic limit as shown in Ref. 51.

### Estimating $\eta(T)$.—

In the main text, we perform a scaling analysis of the $s$-wave pair structure
factor $P_{s}$ using the fact that in a finite system of linear size $L$ it
should scale as $P_{s}(L,T)=L^{2-\eta(T)}f(L/\xi_{\rm KT})$ [50, 51]. For a
KT-type scaling, the exponent $\eta$ is temperature-dependent, and related to
the decay of the real-space pair correlations $C(r=|{\bf i}-{\bf
j}|)\equiv\langle\Delta_{\bf i}^{\dagger}\Delta_{\bf
j}^{\phantom{\dagger}}+\Delta_{\bf j}^{\dagger}\Delta_{\bf
i}^{\phantom{\dagger}}\rangle\sim r^{-\eta(T)}$ [51]. At $T=0$, true long-
range order is manifest and $\eta(0)=0$. As one approaches $T_{c}$ from below,
the exponent smoothly evolves until $\eta(T_{c})=0.25$. In analogy to what was
performed in Ref. 49 for the classical two-dimensional Heisenberg model, we
here extract $\eta(T)$ by dividing the scaling form above for different system
sizes. If comparing lattices with linear sizes $L$ and $L^{\prime}$, the
exponent is thus

$\eta(T)=2-\frac{\ln[P_{s}(L,T)/P_{s}(L^{\prime},T)]}{\ln(L/L^{\prime})}.$
(S1)

Thence, selecting a reference (large) lattice $L$, one can extract $\eta(T)$
by checking the ratio in Eq. S1 at each temperature. This procedure is
illustrated in Fig. S4(b), where the reference lattice has $L=22$. While even
at the relatively low temperatures we investigate it is unclear the
extrapolation to $\eta(T\to 0)=0$ is attainable, the second known limit,
$\eta(T\to T_{c})=0.25$, is readily observed. At high temperatures, $\eta$
saturates at 2, which is consistent with the fact that the pair structure
factor no longer has size-dependence in this regime [see Fig. 4(g) in the main
text].

Figure S4: Extra results of the SU(2) attractive Hubbard model. The resulting
density as the temperature is decreased (a) and the estimation of the
temperature-dependent exponent $\eta$ (b) which governs the real-space decay
of the pairing correlations (see text). As in previous analyses, interaction
strength is $U/t=-4$ and chemical potential is set at $\mu/t=-0.75$.

### Fixed chemical potential.—

The properties of the Hubbard model are often studied in terms of a fixed
density of particles. When performing a grand-canonical, finite-temperature
simulation, this is accomplished by setting the chemical potential in the
regime of parameters that returns the desired electronic density. In
particular, such an approach was applied before in studying the attractive
Hubbard model away from half-filling [49, 50, 51]. Here, we take a different
route and fix the chemical potential (in the main text, $\mu/t=-0.75$). As
displayed in Fig. S4(a), this procedure results in a slight deviation ($\sim
1\%$) of the densities around quarter-filling ($\langle\hat{n}\rangle=0.5$) in
the salient region where the KT transition takes place, with minimal size
effects. As the critical temperature $T_{c}$ possesses minor variations around
this density range [51], such a procedure is justifiable.

##  Appendix H The repulsive Hubbard model away from half-filling

The implication that an estimation of the superconducting critical temperature
$T_{c}$ can be extracted by the regime where the average sign of single
weights in the Monte Carlo sampling drops to zero (furthermore obeying a KT-
scaling form), immediately opens room for the exploration of the same
procedure to tackle the repulsive case, which, unlike the attractive Hubbard
model, displays a sign problem for the total weights (thus preventing using
physical observables to track an eventual finite-temperature transition [13]).
We here use refinements of the standard repulsive Hubbard model [Eq. (1)], by
including a next-nearest neighbor hopping $t_{2}$. This has been shown to
provide more accurate comparisons to experiments on cuprates [65, 66, 67].

Figure S5: The scaling of the spin-resolved sign in the repulsive Hubbard
model. (a) The map of the spin resolved sign in the (a) $T$ vs. $\mu$ and in
the (b) $T$ vs. $\langle\hat{n}\rangle$ planes. In the latter, black markers
depict the densities extracted from the regular $\mu$ grid, and where the
interpolation of the data is performed. (c) A cut of $\langle{\cal
S}_{\sigma}\rangle$ at chemical potential $\mu/t=-2$, corresponding to density
$\sim 0.8$. (d) A KT-scaling form ($\xi_{\rm
BKT}=\exp\\{b/\sqrt{|T-T_{c}|}\\}$) with parameters as displayed, using
lattice sizes ranging from $L=16$ to 26. (e) The corresponding cost function
of the spin-resolved sign in the regime of parameters $(T_{c},b)$. The inset
in (c) displays the $T$-dependence of the $d$-wave pair structure factor in
the regime of parameters where $\langle{\cal S}\rangle>0.25$: system-size
dependence is negligible at this regime, and shows that the extensive regime
is not attainable up until $\langle{\cal S}_{\sigma}\rangle\to 0$. The
parameters are: interaction strength $U/t=6$; next-nearest neighbor hopping
$t_{2}/t=-0.2$, and imaginary-time discretization is $\Delta\tau=1/16$. Figure
S6: The scaling of the total average sign in the repulsive Hubbard model. The
same as in Fig. S5 but for the total sign $\langle{\cal S}\rangle$ instead of
the spin-resolved one $\langle{\cal S}_{\sigma}\rangle$. The same set of
parameters $b$ and $T_{c}$ scale both quantities equally well.

Figure S5 plots the ‘phase diagram’ of the spin resolved sign [as in Fig. 4(a)
and (b) in the main text for the attractive case] with varying temperatures
and both chemical potentials [Fig. S5(a)] and resulting electronic densities
[Fig. S5(b)]. The appearance of a $\langle{\cal S}_{\sigma}\rangle$-dome
suggests the study of criticality is possibly analogous to the attractive
case. By selecting a fixed chemical potential $\mu/t=-2$ (rendering densities
$\langle\hat{n}\rangle\sim 0.8$ across a wide span of temperatures), we
perform a KT-scaling analysis in Fig. S5(c) and S5(d); the corresponding cost
function is given in Fig. S5(e). The trend is remarkably similar to the $U<0$
model, where a finite critical temperature is known to occur. Indetermination
in the location of $T_{c}$ arises from the fact that there is a wide regime in
which $C_{\langle{\cal S}\rangle_{\sigma}}(T_{c},b)$ is small. Future studies
on larger lattices, and considering the possibility that the non-universal
parameter $b$ assumes different values above and below the putative
transition, may render a more accurate estimation.

The repulsive Hubbard model at densities below unity allows us to further
reinforce the argument that the total sign of the weights in the QMC
calculations already provides information about critical behavior. As there is
no symmetry protection that guarantees it to be pinned at 1, we can see in
Fig. S6 that the same set of parameters that scale $\langle{\cal
S}_{\sigma}\rangle$, equally well scale the total sign. Accompanied by similar
results of the Ionic Hubbard model [see Fig. S1], we are led to the conclusion
that in generic models without symmetry protection, the total sign reflects
the same critical properties of the average sign of partial weights, either in
quantum or thermal phase transitions.

Although the critical temperature obtained via the KT-scaling is remarkably
higher than one would expect for describing a cuprate system, similar to the
case of the attractive Hubbard model in the main text, the best estimation of
$T_{c}$ that matches the physical observables (here unattainable) likely comes
from the extrapolation that the parameter $b$ is small, driving it to values
close to $T_{c}\simeq 0.3t$. Yet, assuming a typical bandwidth $W\sim 2$eV for
the copper-oxide compounds [68] would result in critical temperatures much
above room temperature if using a non-interacting bandwidth estimation. There
are two explanations, the first is that the bandwidth for single-particle
excitations is much larger than that in the presence of finite interactions,
and is certainly altered by the large $U/t$ we use. Second, the starting
point, the single-layer Hubbard model with these set of parameters (here using
$U=6t$ and $t^{\prime}=-0.2t$), does not quantitatively characterize the
cuprates. Lastly, interplane hybridization is known to reduce the pairing
fluctuations, and may ultimately drive the thermal transition close to the
ones currently observed in copper-based superconductors.

##  Appendix I Further remarks

We have here demonstrated that the sign problem can be used as a tracker to
understand criticality of fermionic lattice models in four important examples.
That the sign problem depends on the basis in which one describes the system’s
Hamiltonian is well known. What is less understood is whether the average sign
in other local bases can also be used to understand critical properties.
Similarly, an in-depth study using other Hubbard-Stratonovich transformations
[64] would lend insight on the wider generality of our results, on whether
scaling behavior can be similarly inferred.