Dian Wu

Riccardo Rossi

Filippo Vicentini

Nikita Astrakhantsev

Federico Becca

Xiaodong Cao

Juan Carrasquilla

Francesco Ferrari

Antoine Georges

Mohamed Hibat-Allah

Masatoshi Imada

Andreas M Läuchli

Guglielmo Mazzola

Antonio Mezzacapo

Andrew Millis

Javier Robledo Moreno

Titus Neupert

Yusuke Nomura

Jannes Nys

Olivier Parcollet

Rico Pohle

Imelda Romero

Michael Schmid

J Maxwell Silvester

Sandro Sorella

Luca F Tocchio

Lei Wang

Steven R White

Alexander Wietek

Qi Yang

Yiqi Yang

Shiwei Zhang

Giuseppe Carleo

Variational Benchmarks for Quantum Many-Body Problems

The continued development of novel many-body approaches to ground-state problems in physics and chemistry calls for a consistent way to assess its overall progress. Here we introduce a metric of variational accuracy, the V-score, obtained from the variational energy and its variance. We provide the most extensive curated dataset of variational calculations of many-body quantum systems to date, identifying cases where state-of-the-art numerical approaches show limited accuracy, and novel algorithms or computational platforms, such as quantum computing, could provide improved accuracy. The V-score can be used as a metric to assess the progress of quantum variational methods towards quantum advantage for ground-state problems, especially in regimes where classical verifiability is impossible.

I. INTRODUCTION

A key aspect of the quantum many-body problem, for systems ranging from the subatomic to molecules and materials, is determining the ground state properties and energy. With the ground state, one can predict which systems are stable and whether these systems exhibit useful and exotic phases, such as superconductivity or spin liquids. However, due to the exponential complexity of the quantum wave function, finding the ground state of a many-body system can be very challenging, limiting exact numerical studies to a small number of particles. Efficiently solving the general ground-state problem is largely believed to be intractable. However, this does not apply to any particular system or class of systems, which may admit powerful approximations for ground states. Decades of research have focused on devising computational methods to find approximate solutions for specific cases of interest.

These computational methods have widely varying degrees of accuracy, and typically each method is much more successful on some systems than on others. Some of the most widely used methods include quantum Monte Carlo (QMC) [1][START_REF] Becca | Quantum Monte Carlo Approaches for Correlated Systems[END_REF][START_REF] Zhang | Auxiliary-field quantum Monte Carlo at zeroand finite-temperature[END_REF], tensor networks (TN) [START_REF] White | Density matrix formulation for quantum renormalization groups[END_REF][START_REF] Orús | Tensor networks for complex quantum systems[END_REF], and dynamical mean field theory and its extensions (DMFT) [START_REF] Georges | Hubbard model in infinite dimensions[END_REF][START_REF] Georges | Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions[END_REF]. It is known that the applicability of the numerical techniques is negatively affected by the frustration of the quantum system and particle statistics in the case of QMC methods [START_REF] Troyer | Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations[END_REF], by high entanglement for TN [START_REF] Cirac | Matrix product states and projected entangled pair states: Concepts, theorems[END_REF], and by large correlation lengths for DMFT [START_REF] Georges | Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions[END_REF]. Variational approaches based on physically motivated ansatzes [START_REF] Mcmillan | Ground state of liquid He 4[END_REF][START_REF] Ceperley | Monte Carlo simulation of a many-fermion study[END_REF] or neural networks [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF] are not explicitly affected by the aforementioned issues. However, assessing their applicability and accuracy for a given quantum many-body system is more difficult.

Quantum computers provide an alternative platform to attack quantum many-body problems [START_REF] Feynman | Simulating physics with computers[END_REF]. Notably, the dynamics of quantum many-body systems can be efficiently simulated by a digital quantum computer when the initial states are easy to prepare [START_REF] Lloyd | Universal quantum simulators[END_REF]. Besides dynamics, significant attention has been devoted to preparing ground states that are difficult to study with classical algorithms. Quantum algorithms for this task include phase estimation [START_REF] Kitaev | Quantum measurements and the abelian stabilizer problem[END_REF], variational approaches [START_REF] Peruzzo | A variational eigenvalue solver on a photonic quantum processor[END_REF][START_REF] Kandala | Hardwareefficient variational quantum eigensolver for small molecules and quantum magnets[END_REF][START_REF] Cerezo | Variational quantum algorithms[END_REF][START_REF] Bharti | Noisy intermediate-scale quantum algorithms[END_REF], adiabatic passage [START_REF] Farhi | Quantum computation by adiabatic evolution[END_REF], imaginary-time evolution [START_REF] Motta | Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution[END_REF], and subspace and Lanczos methods [START_REF] Parrish | Quantum filter diagonalization: Quantum eigendecomposition without full quantum phase estimation[END_REF][START_REF] Kirby | Exact and efficient Lanczos method on a quantum computer[END_REF].

A fundamental challenge in assessing newly established computational methods, either based on classical or quantum computing, is defining a consistent accuracy metric. Especially for ground-state problems, such a metric is necessary to clearly identify target Hamiltonians of broad interest, which cannot be solved with sufficient accuracy by existing methods. Also, this metric is crucial to quantify the improvements of computational approaches with time. In the context of assessing quantum computing-based methods, this issue pertains to the broader problem of determining in what cases quantum computers have an advantage over classical ones [START_REF] Bravyi | Quantum advantage with shallow circuits[END_REF][START_REF] Bouland | On the complexity and verification of quantum random circuit sampling[END_REF].

Determining a consistent metric for physically and chemically relevant ground-state problems is one of the goals of this work. To this end, we provide the largest to date curated collection of variational and numerically exact results on strongly correlated lattice models obtained by both state-of-the-art and baseline methods. The data we provide include multiple approaches such as exact diagonalization (ED), QMC [1] in the auxiliary field algorithm [START_REF] Blankenbecler | Monte Carlo calculations of coupled boson-fermion systems. I[END_REF][START_REF] Sorella | A novel technique for the simulation of interacting fermion systems[END_REF][START_REF] Imada | Numerical studies on the Hubbard model and the t-J model in one-and twodimensions[END_REF][START_REF] Zhang | Constrained path Monte Carlo method for fermion ground states[END_REF][START_REF] Shi | Ground-state properties of strongly interacting Fermi gases in two dimensions[END_REF], matrix product states (MPS) [START_REF] White | Density matrix formulation for quantum renormalization groups[END_REF], variational wave function formulated on a lattice [START_REF] Horsch | Exact and Monte Carlo studies of Gutzwiller's state for the localised-electron limit in one dimension[END_REF], and neural network-based methods [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF]. In addition to providing the data, we introduce an indicator of the variational accuracy of these results, named V-score, that is suitable for directly comparing classical and quantum computing-based variational approaches. The V-score, obtained as a combination of the mean energy and its variance of a given variational state, allows us to identify what Hamiltonians and regimes are hard to approximate with classical variational methods without prior knowledge of the exact solution. Furthermore, we argue that the V-score can be used as a controlled benchmark to quantify the continued progress of quantum algorithms and quantum hardware to simulate those challenging target Hamiltonians.

II. RESULTS

We focus our study on benchmarking classical and quantum variational algorithms in approximating ground states of quantum many-body systems. On the classical side, these algorithms involve explicitly maintained variational representations of wave functions, such as TN or variational Monte Carlo (VMC)-based approaches. On the quantum side, the variational methods of major interest involve parameterized quantum circuits (PQC) or other state preparation techniques based on local unitary transformations. In all cases, we assume that the methods to be benchmarked allow unbiased estimates of expectation values for Hamiltonians with few-body interactions (k-local operators, in the language of quantum information). Such expectation values can possibly be obtained with a controllable statistical error, as in the case of classical Monte Carlo-based techniques, or as a result of statistical noise due to measurements on quantum hardware.

A. Choice of problems

There is large freedom in the choice of many-body quantum problems that can be used to benchmark computational techniques. In this work, we have decided to focus on lattice Hamiltonians. These are minimal models of strong correlations and typically capture the essence of many physical systems. Lattice models first rose to prominence within classical statistical mechanics with the definition of the Ising model [START_REF] Niss | History of the Lenz-Ising model 1950-1965: From irrelevance to relevance[END_REF]. Within solid-state Figure 1. Sketch of the V-score as a metric of simulation hardness. In this work, we present an extensive dataset of computational results for quantum many-body ground-state problems (for this sketch we have selected a spin-1/2 system on a chain, a kagome, and a pyrochlore lattice). For each Hamiltonian in the dataset, we compute the mean energy and its variance with different variational techniques, including physically motivated ansatzes, neural networks, tensor networks, and parameterized quantum circuits. The energy and the variance are combined into the V-score, a metric of variational accuracy we introduce in the main text. A low V-score is associated with high accuracy. The V-score is then used to rank the Hamiltonians in terms of simulation accuracy, highlighting which quantum many-body models are hard to simulate with current methods.

physics, they find their roots in tight binding approaches to describe the electronic band structure [START_REF] Slater | Simplified LCAO method for the periodic potential problem[END_REF]. More recently, within the second quantization formalism, they are routinely used in different areas of physics to understand the low energy behavior of unconventional quantum phases and transitions among them [START_REF] Sachdev | Quantum Phase Transitions[END_REF][START_REF] Wen | Quantum Field Theory of Many-Body Systems: From the Origin of Sound to an Origin of Light and Electrons[END_REF]. In this regard, the transverse-field Ising model (TFIM) provides the simplest example of a zero-temperature phase transition purely driven by quantum fluctuations between a paramagnet and a ferromagnet as seen, e.g., in the Ising ferromagnet LiHoF 4 [START_REF] Sachdev | Quantum Phase Transitions[END_REF][START_REF] Bitko | Quantum critical behavior for a model magnet[END_REF]. Another prominent example are the various quantum impurity models, in which a localized interacting degree of freedom is embedded into a non-interacting bulk, such as the Anderson impurity model [START_REF] Anderson | Localized magnetic states in metals[END_REF]. Quantum impurity models are central to quantum embedding methods such as DMFT [START_REF] Georges | Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions[END_REF] and also have applications to nanoelectronic devices [START_REF] Kouwenhoven | Revival of the Kondo effect[END_REF]. Their lattice generalizations, such as the Kondo lattice model, describe heavy fermion systems with 4f or 5f atoms, such as Ce or U [START_REF] Hewson | The Kondo Problem to Heavy Fermions[END_REF].

Similarly, the Hubbard model [START_REF] Hubbard | Electron correlations in narrow energy bands[END_REF][START_REF] Kanamori | Electron correlation and ferromagnetism of transition metals[END_REF][START_REF] Gutzwiller | Effect of correlation on the ferromagnetism of transition metals[END_REF] has been widely used to capture the essence of strong correlation in solids and has been proven relevant to the study of hightemperature superconductivity in cuprate compounds, e.g., La 2-x Sr x CuO 4 [START_REF] Lee | Doping a Mott insulator: Physics of high-temperature superconductivity[END_REF], and the Mott metal-insulator transitions in a variety of compounds [START_REF] Imada | Metal-insulator transitions[END_REF]. A descendant of the Hubbard model, the Heisenberg model describes a wide range of magnetic phases, e.g., with ferromagnetic or antiferromagnetic orders [START_REF] Manousakis | The spin-1/2 heisenberg antiferromagnet on a square lattice and its application to the cuprous oxides[END_REF]. In addi-tion, when defined on geometrically frustrated lattices, possibly with anisotropic super-exchange couplings, the Heisenberg model gives rise to a wealth of phenomena, including spin liquid phases with topological order and exotic critical points [START_REF] Balents | Spin liquids in frustrated magnets[END_REF][START_REF] Savary | Quantum spin liquids: A review[END_REF]. In this respect, the rare earth compound YbMgGaO 4 [START_REF] Li | Rare-earth triangular lattice spin liquid: A single-crystal study of YbMgGaO4[END_REF] and the mineral herbertsmithite ZnCu 3 (OH) 6 Cl 2 [START_REF] Norman | Colloquium: Herbertsmithite and the search for the quantum spin liquid[END_REF] have offered examples for unconventional quantum phases on triangular and kagome lattices.

B. V-score

To quantify the accuracy of two or more variational methods applied on the same ground state approximation task, a key indicator is the expectation value of the energy E = Ĥ , an unbiased metric to assess the relative accuracy of variational methods. Given, for example, two independent methods preparing approximate ground states with variational energies E a and E b , the one providing the lower energy can be considered more accurate. From a practical point of view, however, it is preferable to have an absolute metric capable of predicting the accuracy of a method without comparing it with other methods. This would, for instance, allow comparing the performance of a given method on different tasks. Nonetheless, it is unlikely to find such a metric that is provably applicable in all cases, since its existence would also allow the solution of NP-hard problems [START_REF] Barahona | On the computational complexity of ising spin glass models[END_REF]. We are therefore forced to settle for an empirically applicable metric. Moreover, the metric should be easy to estimate with variational methods.

Apart from the mean energy, for most variational methods, we also have a controllable estimate of the energy variance Var E = Ĥ2 -Ĥ 2 . It has the important property that it exactly vanishes if computed on the exact ground state. Therefore, Var E can be used to infer some information about the distance of the variational energy E from the exact, and a priori unknown, groundstate energy E 0 . After early empirical observations [START_REF] Kwon | Effects of three-body and backflow correlations in the twodimensional electron gas[END_REF], it has been shown that Var E scales linearly with the deviation E -E 0 [START_REF] Imada | Path-integral renormalization group method for numerical study of strongly correlated electron systems[END_REF][START_REF] Kashima | Path-integral renormalization group method for numerical study on ground states of strongly correlated electron systems[END_REF][START_REF] Sorella | Generalized Lanczos algorithm for variatonal quantum Monte Carlo[END_REF], so it can be used as a measure of the accuracy of the variational state.

We can use E and Var E to create a dimensionless, intensive combination:

V-score := N Var E (E -E ∞ ) 2 , ( 1
)
where N is the number of degrees of freedom, which is the number of spins for spin models, and the number of particles for fermionic models. The constant E ∞ serves as a zero point of the energy, compensating for any global shift of the energy in the definition of the Hamiltonian. The V-score is dimensionless in energy units and system size for the variational states we consider and is also invariant under energy shifts by construction. A detailed discussion of the definition of the V-score is presented in the Supplementary Materials.

To further justify the definition of the V-score, in Fig. 2, we present a comparison of this quantity against the energy relative error (E -E 0 )/(E ∞ -E 0 ) for a wide range of Hamiltonians and variational methods, where the ground-state energy E 0 is obtained by ED or numerically exact QMC. Despite the great diversity of Hamil-tonians and variational methods considered, the V-score is a remarkably consistent and reliable estimator for the order of magnitude of the energy relative error, as shown by the linear fit in Fig. 2. In the inset of Fig. 2, we show that the same linear fit also well describes classically simulated PQCs, optimized with the variational quantum eigensolver (VQE) algorithm. These results validate the V-score as an absolute performance metric for both classical and quantum variational algorithms, at least for the Hamiltonians and the techniques we consider in the paper.

C. Identifying hard problems

We can now discuss which Hamiltonians are hard for the state-of-the-art variational methods presented in our collection. Given the intrinsic exponential complexity of the problem, it is no longer possible to obtain ED results on larger system sizes. Thus, using the V-score as a guide in this task is crucial. In Fig. 3 we show the Vscore of all methods and models in our dataset. We first select the best method for each Hamiltonian by choosing the one with the lowest variational energy. We then use this method's V-score as an absolute hardness metric in the ground state approximation task, which we refer to as the V-score of this Hamiltonian. The V-score of the best-performing method is marked in bold in Fig. 3. For additional clarity, in Fig. 4, we classify all those results by Hamiltonian types and lattice geometries.

It is well-known that 1D (chain) geometries are easy to solve with density matrix renormalization group (DMRG). The small values of their V-scores in Fig. 3 and Fig. 4 clearly label them as accessible, particularly for spin models. Unfrustrated spin models typically also have small V-scores, ranging from 10 -6 to 10 -4 for TFIM and the Heisenberg model on square lattices with open boundary conditions (OBC). Moreover, these models can be efficiently simulated with unbiased stochastic techniques like QMC, and thus easy to study on classical computers. On the other hand, the V-scores show that frustrated geometries like triangular, kagome, pyrochlore, and the J 1 -J 2 square lattice, as well as fermionic models like the Hubbard model, are the most demanding for variational algorithms. With the help of Fig. 2, we can infer the order of magnitude of the energy relative error from the V-score. When this correspondence is applied to those hard problems, it predicts that we cannot expect an accuracy on the energy better than one or two digits.

D. A perspective on quantum advantage

Many recent theoretical and experimental efforts have been dedicated to showing the computational advantage of quantum computers over classical computers. Informed by theoretical computer science arguments, random circuit sampling has been proposed as a specific task < 10 16 10 8 10 4 10 2 10 1 10 0.5

Harder Easier

V-score 3. V-score rankings of the complete dataset of variational benchmark results. Each radial slice shows the results for a given Hamiltonian. The bold marker in the slice shows the method with the lowest variational energy (rather than the lowest V-score) for that Hamiltonian, and its V-score determines the clockwise ranking of that Hamiltonian. The radial axis uses the doubly log scale to show the results across a wide range of V-score magnitude. The results with V-scores less than 10 -16 are indistinguishable from the exact solutions due to limited numerical precision. The label outside each slice is the name of the Hamiltonian in the dataset, which describes the lattice geometry, the lattice size Ns, the boundary conditions (BC), and other parameters such as the interaction strength and the number of fermions. The '+' marker near the label indicates that an ED solution is available for that Hamiltonian, and the ' * ' marker indicates a numerically exact QMC solution.

-4x4x4_256_ P -3x3x 3_108_P _10 0_P_1 -8 x8_1 92_P + -4x4 _48 _P _ 196 _O
_ 8 _ t 1 2 _ 2 5 6 _ P A _ 1 0 7 _ 8 _ t 1 2 _ U V 1 V 2 _ 1 9 6 _ P _ 0 . 5 _ 1 0 0 _ O _ 2
to show such advantage [START_REF] Boixo | Characterizing quantum supremacy in near-term devices[END_REF][START_REF] Arute | Quantum supremacy using a programmable superconducting processor[END_REF]. Still, it is unclear if current noisy experimental quantum computing platforms can provide a solid and scalable advantage over classical ones [START_REF] Pan | Solving the sampling problem of the Sycamore quantum circuits[END_REF][START_REF] Gao | Limitations of linear crossentropy as a measure for quantum advantage[END_REF][START_REF] Aharonov | A polynomial-time classical algorithm for noisy random circuit sampling[END_REF]. In addition to tasks of purely theoretical interest, there is also growing interest in finding practical quantum advantage [START_REF] Daley | Practical quantum advantage in quantum simulation[END_REF], where quantum devices show a speedup for problems of scientific or technological relevance. The benchmarks introduced in this work belong to the family of approaches that can be useful to assess a quantum advantage that is practically relevant to physics.

In the context of variational ground state algorithms, the V-score can readily be used as a good quality metric to assess quantum advantage. Furthermore, the Vscore can be used as an absolute indicator of hardness for Hamiltonians. In this respect, Hamiltonians with large classical V-scores are identified as hard problems that are not yet satisfactorily solved by classical computers and can be targeted by quantum computations. Last, in the absence of classical verifiability of the quantum solutions, the V-score can benchmark the progress of variational quantum computing-based approaches in solving ground-state problems that are relevant to physics.

In Fig. 3, we show the V-scores of the classical variational methods we have analyzed. From these results, we can infer that there is little room for quantum advantage in one-dimensional geometries, where DMRG is very effective. In higher dimensions, unfrustrated spin models, such as the TFIM or the Heisenberg model on the square lattice, are similarly well approximated by exist- ing classical methods. On the contrary, specific regimes of higher dimensional frustrated spin models constitute a clear challenge for existing classical methods. For example, the pyrochlore or kagome Heisenberg models typically present V-scores significantly larger than their unfrustrated counterparts. A similar scenario emerges for the Hubbard model in two dimensions. In the specific regimes of interaction strengths, geometries, and frustration identified by large V-scores, these models represent natural targets for variational quantum algorithms. Impurity models with multiple bands also represent an ideal terrain for practical quantum advantage as the ability of classical algorithms to simulate them rapidly degrades upon increasing the number of bands, as Fig. 3 shows already for three-band models, and because of their importance for material science. We are also in position to provide an early assessment of the V-scores obtained by the type of variational states that can be efficiently prepared on quantum computers. In this respect, it is encouraging to remark that PQC perform well compared to classical variational methods, as shown in the inset of Fig. 2, at least for the small system sizes we consider, where PQC can be classically simulated and ideally optimized. Applications on quantum hardware are much more challenging because of stochastic fluctuations and noise. However, the baseline of ideally optimized PQC is nonetheless remarkable and promising for applications.

III. DISCUSSION

In this work, we have introduced the V-score, an empirical metric to quantify the absolute accuracy of varia-tional solutions to strongly interacting quantum models. Supplemented with state-of-the-art results obtained by a large variety of numerical methods, this metric allows us to clearly identify models, geometries, and regimes in which existing approaches are currently less accurate. With the introduction of novel computational techniques and improved computing architectures, the outcomes of this analysis will naturally evolve in time, revealing in a certifiable manner the continuous improvements happening in the field. In this respect, the dataset presented here can be a standardized way of taking snapshots of the evolution of quantum many-body techniques with time.

Besides the importance of these benchmarks for future developments in computational techniques based on classical computers, it will be especially interesting to use the V-score to measure the impact of quantum computingbased approaches directly. The hardest classical problem instances identified by large V-scores can be good candidates for studies based on quantum algorithms. In that context, the V-score can be used as a metric to assess progress in quantum variational state preparation, in the absence of classical verifiability.

Supplementary Materials for "Variational Benchmarks for

Quantum Many-Body Problems"

S1. OVERVIEW OF THE MANY-BODY HAMILTONIANS

Here we provide the definitions of the many-body quantum Hamiltonians used for benchmarking purposes.

A. Spin models

The transverse-field Ising model (TFIM) is

Ĥ = J i,j σz i σz j + Γ i σx i , (S1)
where i, j runs over nearest neighbors, Γ is the transverse field strength, and σx , σy , σz are Pauli matrices. In the dataset we use J = 1. The Heisenberg model is

Ĥ = J i,j ,a σa i σa j , (S2)
where a runs in {x, y, z}.

The J 1 -J 2 model is Ĥ = J 1 i,j ,a σa i σa j + J 2 i,j ,a σa i σa j , (S3)
where i, j runs over next-nearest neighbors (diagonals if the lattice is 2D square), and J 2 is the next-nearest neighbor interaction. In the dataset we use J 1 = 1.

B. Fermions

The t-V Hamiltonian is

Ĥ = -t i,j ĉ † i ĉj + ĉ † j ĉi + V i,j ni nj , ( S4
)
where V is the Coulomb repulsive interaction strength, and the number of fermions is fixed to N f . In the dataset we consider t = 1. The Hubbard model is

Ĥ = -t i,j ,σ ĉ † iσ ĉjσ + ĉ † jσ ĉiσ + U i ni↑ ni↓ , (S5)
where U is the on-site interaction strength, and the numbers of fermions are fixed to N ↑ and N ↓ . We only consider the case of N ↑ = N ↓ .

C. Impurity models

A typical Anderson impurity Hamiltonian ĤA contains two parts ĤA := Ĥloc + Ĥbath , (S6)

Ĥloc := {α} 0 α1α2 d † α1 dα2 + {α} U α1α2α3α4 d † α1 d † α2 dα3 dα4 , (S7) Ĥbath := {α} N b l=1 l α1α2 ĉ † lα1 ĉlα2 + {α} N b l=1 ν l α1α2 d † α1 ĉlα2 + h.c. , (
|ν l | 2 = I l dω - Im ∆(ω) π . ( S10
)
We consider two types of interactions that are frequently encountered in DMFT calculations: the singleband Hubbard interaction

Ĥloc = U n↑ n↓ + 0 (n ↑ + n↓ ) , (S11)
where nσ = d † σ dσ is the particle number operator, with σ ∈ {↑, ↓}; the three-band rotationally invariant Kanamori interaction [START_REF] Georges | Strong correlations from Hund's coupling[END_REF] Ĥloc = ĤDD + ĤSF + ĤPH , (S12)

ĤDD := U m nm↑ nm↓ + (U -2J) m >m,σ nmσ nm σ + (U -3J) m >m,σ nmσ nm σ , (S13) ĤSF := J m m d † m↑ dm↓ dm ↑ d † m ↓ + h.c. , (S14)
ĤPH := -J m >m d † m↑ d † m↓ dm ↑ dm ↓ + h.c. , (

S2. V-SCORE

A. Definition and justification of the V-score for lattice models

Given the observed energy expectation E = ψ| Ĥ|ψ and variance Var E = ψ| Ĥ2 |ψ -ψ| Ĥ|ψ 2 of a variational quantum state |ψ , we want to introduce a function of these two quantities, the V-score, as a metric quantifying how close E is to the ground-state energy. As the exact ground-state energy is unknown in general, the mean energy alone is not enough to characterize the quality of a variational optimization. The energy variance is zero for an eigenstate of the Hamiltonian, thus, also for the ground state. Therefore, assuming that the variational optimization does not converge towards an excited state, we can infer that the V-score should be a monotonic function of Var E.

To start with, we investigate how these two observables scale asymptotically with the number of degrees of freedom N . For any well-defined variational state, E scales linearly with N , as the energy is an extensive thermodynamic quantity. To analyze the scaling of Var E, we evoke the cluster property of the variational state. The Hamiltonian Ĥ is written as a sum of

N H = O(N ) local terms, Ĥ = N H i=1 ĥi . If the correlations of these local terms satisfy the cluster property ĥi ĥj -ĥi ĥj ≤ A d(i, j) D+ , ( S16
)
where A, > 0, d is a distance function on the lattice, and D is the space dimension, then Var E scales linearly with N . We can therefore construct a dimensionless number,

N Var E E 2 , ( S17
)
that does not scale with the energy unit or with N asymptotically. Moreover, for variational optimizations converging towards the ground state, Var E can be expected to scale linearly with the energy difference from the ground state [52-54, 62, 63]. Therefore, the dimensionless number in Eq. (S17) linearly quantifies the energy difference as well. However, the expression in Eq. (S17) is still prone to a shift of energy Ĥ → Ĥ + C, where C is an arbitrary constant. In particular, if we fine-tune C such that the ground-state energy is zero, Eq. (S17) can be expected to scale inversely with the energy difference from the ground state. Similarly, if C is such that E = 0, Eq. (S17) is illdefined. To solve this issue, we need to fix a zero point of energy E ∞ in the definition of the V-score:

V-score := N Var E (E -E ∞ ) 2 , (S18)
such that the average energy E of a variational groundstate wave function |ψ is always different from E ∞ . In this work we choose E ∞ to be the energy expectation of a random state (sampled uniformly on the unit sphere surface) in the Hilbert (sub)space H, because an optimized variational state typically has lower energy than a random state. E ∞ is also the energy expectation of a thermal state at infinite temperature restricted to H, which can be computed from the trace of Ĥ:

E ∞ := Tr Ĥ dim H , ( S19
)
where dim H is the dimension of the Hilbert (sub)space.

For the models we consider in this work, the dimension of the Hilbert space is finite for a finite lattice size, and E ∞ is a finite number. We now discuss our choices for the number of degrees of freedom N . For unconstrained spin-1/2 Hilbert spaces, we define N to be equal to the number of lattice sites N s . For the t-V model with fixed particle number, N is defined to be equal to the particle number N f , and we have dim H = B(N s , N f ), where B is the binomial coefficient. For the Hubbard model with fixed particle numbers, N is defined to be equal to the sum of the numbers of spin up and down fermions, N = N ↑ + N ↓ , while dim H = B(N s , N ↑ ) B(N s , N ↓ ). We remark that the Vscore can also be applied to estimations of the lowestenergy excited states in symmetry sectors different from the ground state symmetry sector.

This choice of energy shift E ∞ supposes that the model contains only relevant low-energy degrees of freedom. Actually, if we add many high-energy states, they will have the effect of artificially raising E ∞ without contributing much to the ground state. Therefore, for some models we do not consider in this work (e.g. bosonic models or quantum chemical models), a cutoff on the relevant energy scale must be set into place in order to use this definition of E ∞ . An alternative strategy to define E ∞ for these models would be to use the mean-field energy, which is not affected by the problem of having highenergy states. However, it is generally not straightforward to compute the mean-field energy as it is in itself an NP-hard problem [START_REF] Schuch | Computational complexity of interacting electrons and fundamental limitations of density functional theory[END_REF], there exist many variants of mean-field theory, and it could make weak-coupling calculations beyond mean-field theory artificially hard. Impurity models, which we have introduced in Sec. S1 C, require an adapted definition of E ∞ , see Sec. S2 E.

B. Calculation of E∞ for lattice models

Analytical formulae for specific models

For quantum spin models, we have E ∞ = 0 when the Hamiltonian is written as a sum of Pauli strings with no term proportional to the identity operator, as all Pauli matrices are traceless. For the spinless t-V model with fixed particle number, only the diagonal term V i,j ni nj contributes, and we have

E ∞ = V |E|N f (N f -1) N s (N s -1) , ( S20
)
where |E| is the number of nearest neighbor bonds. For the Hubbard model with fixed particle numbers, only the diagonal term U i ni↑ ni↓ contributes, and we have

E ∞ = U N ↑ N ↓ N s . (S21)
Apart from fixing the number of fermions, in this work we do not consider symmetries of the Hamiltonians when calculating E ∞ .

General case

It is generally efficient to get a numerical estimate of E ∞ with stochastic methods. For simplicity, we limit our discussion to spinless fermions and a shortrange translation-invariant Hamiltonian Ĥ. Considering a Hilbert space H with fixed particle number N f and lattice size N s , we estimate E ∞ by sampling uniformly a bit string x = (x 1 , . . . , x Ns ) with x i ∈ {0, 1} and the constraint i x i = N f , then taking the average

E ∞ = x| Ĥ|x x∼U Ns ,N f (x) , ( S22
)
where |x is an element of the Fock basis and U Ns,N f (x) is the aforementioned uniform distribution. The variance of the estimator of E ∞ can be written as

Var E ∞ = x| Ĥ|x 2 x∼U Ns ,N f (x) -E 2 ∞ . (S23)
Comparing to the physical variance of the energy at infinite temperature

Var H Ĥ := Tr Ĥ2 dim H -E 2 ∞ , (S24)
which can be computed by the same sampling method:

Var H Ĥ = x| Ĥ2 |x x∼U Ns ,N f (x) -E 2 ∞ , (S25)
we have Var E ∞ ≤ Var H Ĥ. As Ĥ is short-range and translation invariant, and as |x satisfies the cluster property being a product state, we have the scalings

E ∞ = O(N ), Var H Ĥ = O(N ), (S26)
which implies that the one-sample stochastic relative error on E ∞ vanishes with increasing system size:

√ Var E ∞ |E ∞ | ≤ Var H Ĥ |E ∞ | = O 1 √ N . ( S27
)
This shows that, in the thermodynamic limit N → ∞, even just one sample is enough to estimate E ∞ . Moreover, the calculation of x| Ĥ|x can be done in a computational time increasing linearly with the number of lattice sites. Therefore, we conclude that the statistical procedure we discussed is efficient.

C. Bounds on the V-score

We now consider bounds on the ratio of the V-score and the energy relative error. The lower bound is obtained when the variational state exactly coincides with an excited state, which is therefore zero. In order to prove an upper bound, we maximize Var E given a fixed mean energy E. When the spectrum is bounded from above, e.g., for finite systems, Var E is maximized, at fixed average energy E, when the variational state |ψ is a linear combination of the ground state and the maximal energy state

|ψ = √ 1 -I |E 0 + √ I |E M , (S28)
where E 0 , |E 0 and E M , |E M are the minimal and the maximal eigenvalue-eigenvector pairs respectively, I ∈ [0, 1] is an interpolation parameter such that E = (1 -I)E 0 + IE M , and we have assumed that both Ĥ and -Ĥ are non-degenerate for simplicity of notation. The energy variance of |ψ is equal to (E M -E)(E -E 0 ). When we express the variance in terms of the V-score, we reach the following bound for the ratio of the V-score and the relative energy error:

0 ≤ V-score (E -E 0 )/(E ∞ -E 0 ) ≤ N (E ∞ -E 0 )(E M -E) (E ∞ -E) 2 .
(S29) Eq. (S29) shows that, while the V-score can arbitrarily underestimate the difference from the ground-state energy, it cannot arbitrarily overestimate it. We remark that the upper bound is linear in system size, as the aforementioned variational state does not respect the cluster property for energy correlations.

D. Scaling of the V-score in the limit of vanishing ground state infidelity

In this section, we adapt the argument from Ref. [START_REF] Taddei | Iterative backflow renormalization procedure for many-body ground-state wave functions of strongly interacting normal Fermi liquids[END_REF] to show the linear scaling of the V-score with the energy relative error in the limit where the ground state infidelity goes to zero. We remark that this result is not directly applicable to large system sizes as the ground state infidelity is expected to grow exponentially with system size for a fixed energy relative error.

We consider a variational state parameterized by a control parameter c ≥ 0:

|ψ(c) := i ψ i (c) |E i , (S30)
where ψ i (c) ∈ C and |E i are energy eigenstates with eigenvalue E i . We suppose E 0 ≤ E i , and we define the set of ground state indices as G := {i | E i = E 0 }. We introduce the ground state infidelity

I(c) := i / ∈G |ψ i (c)| 2 i |ψ i (c)| 2 . (S31)
We suppose that when the control parameter c goes to infinity, the ground state infidelity goes to zero:

lim c→∞ I(c) = 0. (S32)
We also suppose that I(c) = 0 for every c.

In the following, we determine the scaling of the Vscore when c → ∞. For an operator Ô, we define

Ô c := ψ(c)| Ô|ψ(c) ψ(c)|ψ(c) , (S33) Ô ex,c := i / ∈G |ψ i (c)| 2 E i | Ô|E i i / ∈G |ψ i (c)| 2 . ( S34
)
Then we have

Ĥ c -E 0 = ( Ĥ ex,c -E 0 )I(c). (S35)
We also define

σ 2 c := Ĥ2 c -Ĥ 2 c (S36) = ( Ĥ -E 0 ) 2 c -( Ĥ c -E 0 ) 2 (S37) = ( Ĥ -E 0 ) 2 ex,c I(c) -( Ĥ ex,c -E 0 )I(c) 2 (S38) = ( Ĥ -E 0 ) 2 ex,c I(c) 1 + O I(c) , (S39)
where we used 0

< ( Ĥ ex,c -E 0 ) 2 ≤ ( Ĥ -E 0 ) 2 ex,c . We then write lim c→∞ V-score ( Ĥ c -E 0 )/(E ∞ -E 0 ) = lim c→∞ N σ 2 c (E ∞ -E 0 ) ( Ĥ c -E 0 )(E ∞ -Ĥ c ) 2 (S40) = lim c→∞ N ( Ĥ -E 0 ) 2 ex,c (E ∞ -E 0 ) ( Ĥ ex,c -E 0 )(E ∞ -Ĥ c ) 2 (S41) = lim c→∞ N ( Ĥ -E 0 ) 2 ex,c ( Ĥ ex,c -E 0 )(E ∞ -E 0 ) . ( S42
)
If we suppose that the following limits exist:

σ 2 ex := lim c→∞ Ĥ2 ex,c -Ĥ 2 ex,c , (S43)
E ex := lim c→∞ Ĥ ex,c , (S44)
then we have

lim c→∞ V-score ( Ĥ c -E 0 )/(E ∞ -E 0 ) = N σ 2 ex + (E ex -E 0 ) 2 (E ex -E 0 )(E ∞ -E 0 )
.

(S45) This limit is non-singular because

E 0 + ∆ ≤ E ex ≤ E M , (S46) 0 ≤ σ 2 ex ≤ 1 4 (E M -E 0 -∆) 2 , ( S47
)
where ∆ is the gap between the ground state subspace and the rest of the spectrum, and E M is the maximal energy. If we suppose further that E ex = E 0 + ∆ and σ 2 ex = 0, which is valid, e.g., for two-level systems or for a thermal state at temperature T = 1/c, then we have

lim c→∞ V-score ( Ĥ c -E 0 )/(E ∞ -E 0 ) = N ∆ E ∞ -E 0 . ( S48
)

E. E∞ for impurity models

As mentioned in Sec. S2 A, when applied to impurity models, the definition of E ∞ should be modified such that |x in Eq. (S22) is restricted to the low-energy subspace. For impurity models, there is a natural way to filter out high-energy states that do not contribute to the ground state. Noticing that the fast convergence of density of bath sites with negative (positive) on-site energy to be fully occupied (empty), |x can be restricted to the lowenergy subspace as |x = |o ⊗ |e ⊗ |i , where |o (|e ) denotes the product state of completely occupied (empty) bath sites, |i ∈ H imp ⊗ H a is a product state belonging to the Hilbert space composed by the impurity and the active bath site. The active bath site is determined by the one with the smallest absolute on-site energy. Following Eq. (S19), E ∞ is modified for impurity models as

E ∞ = Tr ρb ⊗ ρi ĤA dim H imp , (

S3. OVERVIEW OF THE NUMERICAL METHODS

A. Exact diagonalization

The quantum many-body problem can be solved numerically to arbitrary precision on small lattice sizes using exact diagonalization without any approximation. Typically, an iterative algorithm like the Lanczos algorithm [START_REF] Lanczos | An iteration method for the solution of the eigenvalue problem of linear differential and integral operators[END_REF] is used to solve for the eigenvalues and eigenvectors of the static Schrödinger equation

Ĥ |ψ = E |ψ . (S50)
The Lanczos algorithm has proven to be a reliable tool for computing ground-state energies up to machine precision. Hence, data retrieved from exact diagonalization can be considered exact. To achieve the currently largest system sizes, both memory and CPU time limitations must be dealt with. To avoid memory bottlenecks, the matrix-vector multiplication operations in the Lanczos algorithm are performed "on-the-fly", i.e., without storing the Hamiltonian matrix, neither in full nor in sparse format, but by implementing a matrix-vector multiplication function. Moreover, using Hamiltonian symmetries allows block-diagonalization to reduce the memory footprint further. Using a symmetry-adapted basis requires efficient algorithms to evaluate matrix elements in this basis, e.g. sublattice coding algorithms. Finally, also large-scale parallelization for distributed memory computers is required, which poses challenges in managing load balancing. The necessary techniques required to achieve the currently largest system sizes are explained in detail in Ref. [START_REF] Wietek | Sublattice coding algorithm and distributed memory parallelization for large-scale exact diagonalizations of quantum many-body systems[END_REF].

B. Tensor networks

Matrix product states and density matrix renormalization group

The density matrix renormalization group (DMRG) is a variational technique, first introduced by White in 1992 to accurately describe the ground-state properties of onedimensional (1D) quantum lattices [START_REF] White | Density matrix formulation for quantum renormalization groups[END_REF]. While exact diagonalization methods operate in an exponentially large basis, DMRG works with the degrees of freedom tied to a few sites at a time. At the heart of the DMRG algorithm is the matrix product state (MPS) ansatz [START_REF] Schollwöck | The density-matrix renormalization group in the age of matrix product states[END_REF], which represents the many-indexed wave function as a chain of tensors, one for each site, with links connecting the sites in a 1D layout. DMRG is an extremely efficient procedure for optimizing the coefficients of the MPS. The required bond dimension, i.e. dimension of the indices linking the tensors, is determined by the degree of entanglement in the state being described. In the calculations presented here, which do not utilize parallelization beyond a single node, bond dimensions are limited to about 10, 000. The required bond dimension for the ground state of a model system is thus tied to the area law of entanglement, which states that the entanglement entropy of a bipartition of a system varies as the size of the boundary rather than the volume of either subsystem [START_REF] Eisert | Area laws for the entanglement entropy -a review[END_REF]. DMRG is ideal for gapped 1D systems, where the area law implies that the bond dimension (for a fixed error) is independent of the length of the chain. For a gapless chain, there is a slow logarithmic growth of the bond dimension with length; nevertheless, spin-1/2 chains with lengths in the thousands are still rather easy on a laptop.

For two-dimensional (2D) clusters, the area law implies that the bond dimension is independent of the length but grows exponentially with the width. Despite the exponential, the general efficiency and robustness of DMRG makes it one of the most powerful and versatile methods for studying many 2D lattice models [START_REF] Stoudenmire | Studying twodimensional systems with the density matrix renormalization group[END_REF]. Success has required developing a variety of specific techniques and "tricks"; for example, the standard two-site DMRG algorithm gives a rough measure of the error associated with using a finite bond dimension, called the truncation error. Even though the truncation error is a crude approximation of the true error, protocols for extrapolating the truncation error to zero can give greatly improved energy estimates and approximate errors on those estimates. Another trick is to rely as much as possible on measurements of local quantities rather than longdistance correlation functions, which are determined by DMRG much less accurately. Similar information to correlation functions can be obtained from a local perturbation of the system followed by the measurement of local quantities away from the perturbation by following the linear response theory. Note that perturbations, e.g. a global antiferromagnetic field breaking the SU(2) symmetry of an antiferromagnet, can sometimes reduce entanglement, making the calculation easier.

The energy variance and the related V-score presented in this paper deviates from the usual DMRG protocols, but they provide a natural way to compare different algorithms. The calculation of the variance of the energy is straightforward in DMRG, but it is much more costly to compute than the truncation error. Aside from the cost, extrapolations to zero variance are a potential improvement over truncation error extrapolations. To mitigate the cost, a two-site variance has been introduced, which is potentially more robust than the truncation error, and can be used in the one-site DMRG algorithm [START_REF] Hubig | Error estimates for extrapolations with matrix-product states[END_REF]. The two-site variance has a cost similar to the rest of a twosite DMRG sweep. Here, however, we only report fullvariance calculations and forgo extrapolations to allow comparisons with other methods.

A key limitation in calculating the full variance is that one can run out of memory. We use a matrix product operator (MPO) form for the Hamiltonian, and calculating the variance can be done using the square of this MPO -except that this tends to require large amounts of memory. A useful trick is to break the Hamilonian MPO into k pieces, each with smaller bond dimension. Then one needs to sum k 2 terms which may be calculated in parallel.

Fork tensor product states for impurity models

The many-body wave function of the impurity models is parameterized by the Fork tensor product states (FTPS) [START_REF] Bauernfeind | Fork tensor-product states: Efficient multiorbital real-time DMFT solver[END_REF]. Compared to MPS, which has a chain geometry, FTPS avoids the artificial long-range interaction that is detrimental to MPS by explicitly separating bath degrees of freedom belonging to different bands. Hence, the FTPS is expected to efficiently capture the entanglement structure of multiorbital problems. Furthermore, the bipartite nature of the FTPS makes it possible to extend the efficient DMRG algorithm developed for MPS to FTPS. The ground state is found in our calculations by the single-site DMRG algorithm supplied with a subspace expansion [START_REF] Hubig | Strictly single-site DMRG algorithm with subspace expansion[END_REF]. Except for the three-band model with spin-orbital coupling, which has only U(1) symmetry for the charge sector, all calculations are performed under the global U(1) symmetries for the charge and spin sectors. We first perform 30 DMRG sweeps without symmetries with a relatively low bond dimension to find the correct charge and spin sector. Then, the charge and the spin quantum numbers are fixed to be Nα and Ŝz α for each spin-orbital α, respectively. Finally, the ground state is found by another 60 DMRG sweeps in the fixed quantum number sector with a maximum bond dimension m = 100 for single-band models and m = 350 for three-band models.

C. Variational Monte Carlo

Variational Monte Carlo (VMC) methods are a family of computational methods that do not suffer from the sign problem and whose computational cost is tractable. In particular, VMC combines a variational encoding of the wave function, to reduce the memory complexity associated with storing the wave function, with Monte Carlo techniques which lower the computational complexity. This approach was originally introduced to treat models in the continuum such as the helium atom or the electron gas [START_REF] Mcmillan | Ground state of liquid He 4[END_REF][START_REF] Ceperley | Monte Carlo simulation of a many-fermion study[END_REF], and later adapted to find the ground states of lattice systems as those discussed in the main text of this manuscript [START_REF] Horsch | Exact and Monte Carlo studies of Gutzwiller's state for the localised-electron limit in one dimension[END_REF][START_REF] Shiba | Properties of strongly correlated fermi liquid in valence fluctuation system -a variational Monte-Carlo study[END_REF]. Since then, more sophisticated methods have been proposed on various approximation levels [START_REF] Yokoyama | Variational Monte-Carlo studies of Hubbard model. I[END_REF][START_REF] Gros | Antiferromagnetic cor-relations in almost-localized fermi liquids[END_REF][START_REF] Capriotti | Resonating valence bond wave functions for strongly frustrated spin systems[END_REF][START_REF] Tahara | Variational Monte Carlo method combined with quantum-number projection and multivariable optimization[END_REF][START_REF] Misawa | mVMC -open-source software for many-variable variational Monte Carlo method[END_REF][START_REF] Nomura | Restricted Boltzmann machine learning for solving strongly correlated quantum systems[END_REF]. Nowadays, several open-source software implementations of those algorithms are available, such as mVMC [START_REF] Misawa | mVMC -open-source software for many-variable variational Monte Carlo method[END_REF] and NetKet [START_REF] Vicentini | NetKet 3: Machine learning toolbox for many-body quantum systems[END_REF].

In variational approaches such as VMC, a quantum state |ψ θ is encoded into a parameterized function ψ θ (x) usually referred to as the variational ansatz. This function takes vectors x from a certain basis of the Hilbert space as input and outputs the complex wave function amplitudes such that

|ψ θ := x ψ θ (x) |x . (S51)
If the function ψ θ is fixed, then a quantum state is uniquely identified by the vector of variational parameters θ. Variational ansatzes are usually chosen such that the number of variational parameters θ increases at most only algebraically with the number of degrees of freedom in the system (in most cases linearly or quadratically). This provides an exponential reduction in the memory complexity with respect to storing the full wave function, which is a vector in an exponentially large Hilbert space.

Even when the wave function is encoded into a small vector of a few variational parameters θ, the expectation value of physical quantities involves two sums over the entire basis of the Hilbert space, leading to exponential computational complexity. To work around this issue, expectation values are computed through statistical averages of stochastic estimators. Instead of summing over the whole Hilbert space, only a few basis elements are selected through rigorous Monte Carlo sampling techniques such as Markov chain Monte Carlo (MCMC). In practice, the expectation value of the energy over the state, E θ = ψ θ | Ĥ|ψ θ , is approximated by the following expectation value [START_REF] Becca | Quantum Monte Carlo Approaches for Correlated Systems[END_REF]:

E θ ≈ E x∼P θ (x) E loc θ (x) , (S52)
where the samples {x} are distributed according to the Born probability

P θ (x) := |ψ θ (x)| 2 ψ θ |ψ θ , (S53)
and

E loc θ (x) := x x| Ĥ|x ψ θ (x ) ψ θ (x) (S54)
is called the local energy. Assuming that the Hamiltonian Ĥ has only a polynomial number of non-zero entries x| Ĥ|x for every row x|, the computational complexity of this local estimator is algebraic.

In general, when using VMC to determine the ground state of a Hamiltonian Ĥ one updates the variational parameters θ according to the conjugate gradient of the energy [START_REF] Amari | Natural gradient works efficiently in learning[END_REF][START_REF] Sorella | Green function Monte Carlo with stochastic reconfiguration[END_REF]. If ψ θ (x) has real-valued parameters or if it is holomorphic [START_REF]For non-holomorphic ansatzes with complex-valued parameters the formula gains a second term ∝ ∂ log ψ * θ (x) ∂θ[END_REF], the gradient is estimated using the formula

∂E θ ∂θ * ≈ E x∼P θ (x) E loc θ (x) -E θ ∂ log ψ * θ (x) ∂θ * , (S55)
where

∂ log ψ θ (x)
∂θ is computed by automatic differentiation (AD) [START_REF] Baydin | Automatic differentiation in machine learning: A survey[END_REF] in modern softwares. This estimator also has the useful property that ∂E θ ∂θ * = 0 when the variational state hits an eigenstate of the Hamiltonian such as the ground state, meaning that the optimization will stop if convergence is reached.

The simplest first-order optimization scheme is the stochastic gradient descent (SGD):

θ ← θ -η ∂E θ ∂θ * , ( S56
)
where η is a sufficiently small positive number called the learning rate. More elaborate optimization schemes that involve accumulating momentums of the gradient, such as Adam [START_REF] Kingma | 3rd International Conference on Learning Representations[END_REF], can also be employed.

A more advanced optimization method that leverages information about the local geometry of the variational manifold [START_REF] Provost | Riemannian structure on manifolds of quantum states[END_REF][START_REF] Stokes | Quantum natural gradient[END_REF] is often used. This method, known as natural gradient [START_REF] Amari | Natural gradient works efficiently in learning[END_REF] in the machine learning literature and as stochastic reconfiguration (SR) [START_REF] Sorella | Green function Monte Carlo with stochastic reconfiguration[END_REF] in the VMC literature, approximates the imaginary-time evolution by exp(-τ Ĥ) for the Hamiltonian Ĥ with sufficiently large τ to reach the ground state.

After the stochastic optimization of the energy, the variational wave function can be improved further by applying the Lanczos operator 1 + α Ĥ to |ψ θ once or twice with the variational parameter α and is employed in some cases of the present benchmark.

Physically motivated ansatzes

In fermionic systems, which are one of the main focuses of this paper, the variational ansatz takes the form |ψ = Ĉ |φ , where |φ is an uncorrelated fermionic state and Ĉ denotes a generic correlator. Centering the discussion on systems of spinful electrons, the uncorrelated state |φ is given either by a Slater determinant [START_REF] Yokoyama | Variational Monte-Carlo studies of Hubbard model. I[END_REF][START_REF] Gros | Antiferromagnetic cor-relations in almost-localized fermi liquids[END_REF][START_REF] Ferrari | Charge-density waves in kagome-lattice extended Hubbard models at the van Hove filling[END_REF] or by a pair-product (PP) wave function (also known as geminal wave function) [START_REF] Capriotti | Resonating valence bond wave functions for strongly frustrated spin systems[END_REF][START_REF] Tahara | Variational Monte Carlo method combined with quantum-number projection and multivariable optimization[END_REF][START_REF] Gros | Superconductivity in correlated wave functions[END_REF][START_REF] Astrakhantsev | Brokensymmetry ground states of the Heisenberg model on the pyrochlore lattice[END_REF]. The variational PP wave function is defined by

φ(x) := x|φ = x|   ijσσ f σσ ij ĉ † iσ ĉ † jσ   Ne 2 |φ 0 , (S57)
where |φ 0 is the vacuum state, N e is the number of electrons and the pair amplitudes f σσ ij form a matrix of 4N 2 variational parameters on N sites, which depend on the spatial coordinates i and j, and the spins σ and σ . The spin-dependence of f σσ ij can be chosen such that either singlet or triplet components appear in the wave function. The PP state, in general, contains Slater determinants as a subset and offers higher flexibility and accuracy. It can accommodate the Hartree-Fock-Bogoliubov type wave function with magnetic, charge, and superconducting orders [START_REF] Tahara | Variational Monte Carlo method combined with quantum-number projection and multivariable optimization[END_REF], and paramagnetic metals as well, in a unified and flexible fashion.

The correlator Ĉ operates to save the exponentially large number of basis functions, and frequently used examples in the present benchmark are given by introducing and applying various physically adequate operators in Ĉ, such as Ĉ := LS Lq PG PJc PJs .

(S58)

The correlation factors

PG := exp i g G i ni↑ ni↓ , (S59)
PJc := exp   i<j g Jc ij ni nj   , (S60) PJs := exp   i<j g Js ij Ŝz i Ŝz j   (S61)
are the Gutzwiller factor [START_REF] Gutzwiller | Effect of correlation on the ferromagnetism of transition metals[END_REF], the long-range Jastrow correlation factor [START_REF] Jastrow | Many-body problem with strong forces[END_REF][START_REF] Capello | Variatinal description of Mott insulators[END_REF], and the long-range spin Jastrow correlation factor [START_REF] Huse | Simple variational wave functions for two-dimensional heisenberg spin-1/2 antiferromagnets[END_REF], with spatially dependent variational parameters g G i , g Jc ij , and g Js ij , respectively; Ŝz i = 1 2 (n i↑ -ni↓ ) is the z component of the local spin operator. In practice, the translational symmetry is often imposed on the Gutzwiller and the Jastrow factors in order to reduce the number of independent variational parameters.

Furthermore, since on finite sizes the exact ground state possesses all the symmetries of the Hamiltonian, while the variational state |φ may break them, the quantum number projections L are also considered to restore the symmetries. In Eq. (S58), LS and Lq are examples of projection operators which enforce the fixed total spin S and momentum q, respectively [START_REF] Tahara | Variational Monte Carlo method combined with quantum-number projection and multivariable optimization[END_REF], and Lq is used in the case of periodic boundary conditions.

In the above example, the variational parameters are contained both in the correlator Ĉ (i.e. g G i , g Jc ij , and g Js ij )

and in the uncorrelated state |φ (i.e. f σσ ij ). They are optimized to better approximate the ground state by lowering the numerically evaluated energy. By fully optimizing f σσ ij including the long-range part, it is known that not only the symmetry-broken Hartree-Fock-Bogoliubov states and simple single-particle noninteracting Slater determinants, but also correlated metallic states such as Tomonaga-Luttinger liquid can be represented [START_REF] Liang | Some new variational resonating-valence-bond-type wave functions for the spin-1/2 antiferromagnetic heisenberg model on a square lattice[END_REF][START_REF] Kaneko | Improved multivariable variational Monte Carlo method examined by high-precision calculations of one-dimensional Hubbard model[END_REF]. In order to lower the computational cost by reducing the number of independent variational parameters, one can assume that f σσ ij has a sublattice structure such that it depends only on the relative position vector r i -r j and a sublattice index of site j which is denoted as η(j).

Then one can rewrite it as f σσ η(j) (r i -r j ). In many cases, one can employ variational states which do not break any translational symmetry and assume a fully translational invariance (1 × 1 sublattice structure, in the case of 2D models), restricting to singlet pairings only (i.e.

f σσ ij = -f σ ,σ ij (1-δ σσ ))
. Antiferromagnetic states can be included by extending the sublattice structure to a 2 × 2 (or larger) unit cell. In most of the present benchmark studies, we do not restrict the sublattice size; namely, the sublattice size is the same as the full lattice. In some cases, PJs and LS are omitted. One can also improve the wave function by implementing dependence on the local density of |x (see Ref. [START_REF] Ido | Unconventional dual 1D-2D quantum spin liquid revealed by ab initio studies on organic solids family[END_REF]). The advantage of using this scheme together with Eqs. (S57) and (S58) is that the quantum entanglement beyond the area law can be taken into account.

In an alternative approach, the parameters f σσ ij for |φ are obtained starting from an auxiliary non-interacting Hamiltonian, containing hopping and pairing amplitudes [START_REF] Capriotti | Resonating valence bond wave functions for strongly frustrated spin systems[END_REF][START_REF] Hu | Direct evidence for a gapless Z2 spin liquid by frustrating néel antiferromagnetism[END_REF][START_REF] Tocchio | Metallic and insulating stripes and their relation with superconductivity in the doped Hubbard model[END_REF]. The most general form is given by

Ĥaux := ijσσ χ σσ ij ĉ † iσ ĉjσ + ∆ σσ ij ĉ † iσ ĉ † jσ + h.c. , (S62)
where χ σσ ij and ∆ σσ ij are hopping and pairing terms, respectively. Then the ground state of Ĥaux is a PP wave function in Eq. (S57), with the pair amplitudes f σσ ij determined by the parameters χ σσ ij and ∆ σσ ij . The advantage of this approach is two-fold. First, long-range pair amplitudes can be obtained within a short-range parametrization of χ σσ ij and ∆ σσ ij , thus avoiding delicate optimizations of the long-range tails of f σσ ij . Second, symmetries may be imposed directly on the auxiliary Hamiltonian, thus avoiding including further projectors, e.g., LS and Lq . For example, by restricting to "singlet" hoppings (χ σσ ij = δ σσ χ ij with χ ij = χ * ji ) and pairings (∆ σσ ij = (1 -δ σσ )∆ ij with ∆ ij = ∆ ji ), the PP wave function is already a singlet, while by imposing translational symmetry in the dependence of these parameters on the lattice sites, the PP state is translationally symmetric automatically. Within this framework, the magnetic order can be described by including on-site hoppings which mimic the presence of a (site-dependent)

Zeeman field h i = (h x i , h y i , h z i ), with h z i = 1 2 (χ ↑↑ ii -χ ↓↓ ii ), h x i = χ ↑↓ ii + χ ↓↑
ii , and h y i = i(χ ↑↓ ii -χ ↓↑ ii ) (see for example Refs. [START_REF] Iqbal | Spin liquid nature in the Heisenberg J1-J2 triangular antiferromagnet[END_REF][START_REF] Ferrari | Spectral signatures of fractionalization in the frustrated Heisenberg model on the square lattice[END_REF]).

The general parametrization in Eq. (S62) allows us to include site-dependent terms to describe charge and/or spin inhomogeneities (i.e. stripes) [START_REF] Tocchio | Metallic and insulating stripes and their relation with superconductivity in the doped Hubbard model[END_REF][START_REF] Marino | Stripes in the extended t-t Hubbard model: A variational Monte Carlo analysis[END_REF]. Jastrow factors, when included, are chosen to have translational and rotational symmetries. Furthermore, backflow correlations can be included, redefining the orbitals of Eq. (S62) on the basis of the many-body electronic configuration [START_REF] Tocchio | Role of backflow correlations for the nonmagnetic phase of the t-t Hubbard model[END_REF][START_REF] Tocchio | Backflow correlations in the Hubbard model: An efficient tool for the study of the metal-insulator transition and the large-U limit[END_REF]. Analogous variational wave functions can be employed to study localized spins systems, e.g., frustrated Heisenberg models on different lattice geometries [START_REF] Capriotti | Resonating valence bond wave functions for strongly frustrated spin systems[END_REF][START_REF] Hu | Direct evidence for a gapless Z2 spin liquid by frustrating néel antiferromagnetism[END_REF][START_REF] Iqbal | Spin liquid nature in the Heisenberg J1-J2 triangular antiferromagnet[END_REF][START_REF] Iqbal | Gapless spin-liquid phase in the kagome spin-1/2 heisenberg antiferromagnet[END_REF][START_REF] Ferrari | Dynamical structure factor of the J1-J2 heisenberg model on the triangular lattice: Magnons, spinons, and gauge fields[END_REF].

Notice that in the context of spin models, the Gutzwiller factor is replaced by the Gutzwiller projector PG ∞ := i ni (2 -ni ), which enforces single occupancy on each lattice site, on top of a fermionic uncorrelated PP state |φ . In addition, analogously to the electronic case previously discussed, a spin Jastrow factor PJs can be included, and quantum number projectors LS , Lq can be applied to enforce lattice symmetries. Also, in this case, the f σσ ij pair amplitudes of the PP state can be assumed to be direct variational parameters (with certain symmetry constraints) or defined through an auxiliary Hamiltonian like the one of Eq. (S62).

In addition, we mention that a possible bias inevitable in the original wave function can be progressively removed by adding a correlator M implemented as a restricted Boltzmann machine (RBM) [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF][START_REF] Nomura | Restricted Boltzmann machine learning for solving strongly correlated quantum systems[END_REF] on top of the variational state, i.e., |ψ = Ĉ M |φ . In this case, Jastrow factors can be omitted to save computational cost. We refer the reader to the following subsection for more details about the RBM. Recently, the VMC approach with the RBM implementation has been successfully applied to reveal the nature of the quantum spin liquids [START_REF] Ido | Unconventional dual 1D-2D quantum spin liquid revealed by ab initio studies on organic solids family[END_REF][START_REF] Ferrari | Neural Gutzwillerprojected variational wave functions[END_REF][START_REF] Nomura | Dirac-type nodal spin liquid revealed by refined quantum many-body solver using neural-network wave function, correlation ratio, and level spectroscopy[END_REF].

Neural quantum states

Neural networks are structured nonlinear functions that can universally approximate any well-behaved function. Motivated by their recent success of efficiently representing probability distributions in machine learning tasks, they have also been used as variational ansatzes to represent quantum states. They typically have more parameters and higher expressive power than traditional ansatzes. The first proposed neural quantum state (NQS) is a restricted Boltzmann machine (RBM) [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF] applied first to quantum spin Hamiltonians and then extended to fermionic systems [START_REF] Nomura | Restricted Boltzmann machine learning for solving strongly correlated quantum systems[END_REF] as is introduced in the last subsection. The RBM is also equivalent to a multilayer perceptron (MLP) with two layers. The finite temperature path integral formalism was shown to exactly map to a deep Boltzmann machine (DBM), where the Trotter-Suzuki layers correspond to the multiple hidden layers in the DBM [START_REF] Carleo | Constructing exact representations of quantum many-body systems with deep neural networks[END_REF]. The correlator M represented by NQS can also benefit from techniques developed for traditional ansatzes, such as backflow correlation [START_REF] Luo | Backflow transformations via neural networks for quantum many-body wave functions[END_REF] and Gutzwiller projection [START_REF] Ferrari | Neural Gutzwillerprojected variational wave functions[END_REF].

Inspired by the development in computer vision, NQS based on convolutional neural networks (CNN) [START_REF] Choo | Two-dimensional frustrated J1-J2 model studied with neural network quantum states[END_REF] takes advantage of the locality and the translational symmetry of physical systems on regular lattices. We can also restore symmetries by quantum number projections discussed in the previous subsection. Some of the data in this paper are obtained by RBM combined with quantum number projections [START_REF] Nomura | Helping restricted Boltzmann machines with quantum-state representation by restoring symmetry[END_REF]. Group convolutions [START_REF] Luo | Gauge equivariant neural networks for quantum lattice gauge theories[END_REF][START_REF] Roth | Group convolutional neural networks improve quantum state accuracy[END_REF] further generalize them to richer symmetries including rotations, reflections, and gauge transformations.

Another direction for developing NQS is perfect sampling from complicated many-body probability distributions, which avoids the issue of autocorrelation time in MCMC. It can be achieved by autoregressive neural networks (ARNN) [START_REF] Sharir | Deep autoregressive models for the efficient variational simulation of many-body quantum systems[END_REF], which are exact likelihood models. They rely on decomposing a joint probability into conditional probabilities that can be sequentially sampled, and they can comprise either dense or convolutional layers. Recurrent neural networks (RNN) [START_REF] Hibat-Allah | Recurrent neural network wave functions[END_REF][START_REF] Roth | Iterative retraining of quantum spin models using recurrent neural networks[END_REF] maintain a memory of past information during the sampling procedure, which makes them suitable to capture long-range correlations, as suggested by their success in natural language processing.

One interesting feature about the RNN is its capability to encode translation-invariant properties of the bulk of a quantum system [START_REF] Roth | Iterative retraining of quantum spin models using recurrent neural networks[END_REF]. Moreover, RNNs can be extended to multiple spatial dimensions. In particular, one can construct 2D RNNs that were shown to be competitive with DMRG, and also cheaper in terms of computational complexity compared to projected-entangled pair states (PEPS) [START_REF] Hibat-Allah | Recurrent neural network wave functions[END_REF]. RNNs can use much fewer parameters than other architectures to encode information in a large spatial area, thanks to weight sharing between RNN cells at all sites. In a similar fashion to tensor networks, tensorized versions of RNNs have been used to provide higher expressive power [START_REF] Hibat-Allah | Variational neural annealing[END_REF][START_REF] Hibat-Allah | Supplementing recurrent neural network wave functions with symmetry and annealing to improve accuracy[END_REF][START_REF] Wu | From tensor network quantum states to tensorial recurrent neural networks[END_REF] and more accurate estimations of ground-state energies that outperform DMRG in certain regimes [START_REF] Hibat-Allah | Supplementing recurrent neural network wave functions with symmetry and annealing to improve accuracy[END_REF][START_REF] Wu | From tensor network quantum states to tensorial recurrent neural networks[END_REF]. They were also shown to be able to encode the area law of entanglement [START_REF] Wu | From tensor network quantum states to tensorial recurrent neural networks[END_REF]. Symmetries, such as U(1) symmetry and spatial symmetries, can be imposed in RNNs to improve the quality of variational calculations [START_REF] Hibat-Allah | Recurrent neural network wave functions[END_REF][START_REF] Hibat-Allah | Supplementing recurrent neural network wave functions with symmetry and annealing to improve accuracy[END_REF]. Being the recent state of the art for many machine learning tasks, transformers have also been proposed for NQS as they have a more flexible autoregressive architecture [START_REF] Luo | Gauge equivariant neural networks for quantum lattice gauge theories[END_REF][START_REF] Zhang | Transformer quantum state: A multi-purpose model for quantum many-body problems[END_REF]. Additionally, it is worth noting that thermal-like fluctuations can be added to the training of autoregressive models in the hope of escaping local minima that can be encountered when studying disordered or frustrated systems [START_REF] Roth | Iterative retraining of quantum spin models using recurrent neural networks[END_REF][START_REF] Hibat-Allah | Variational neural annealing[END_REF][START_REF] Hibat-Allah | Supplementing recurrent neural network wave functions with symmetry and annealing to improve accuracy[END_REF].

Neural networks can be used as well to simulate fermionic systems within the second quantization formalism [START_REF] Choo | Fermionic neuralnetwork states for ab-initio electronic structure[END_REF][START_REF] Yoshioka | Solving quasiparticle band spectra of real solids using neuralnetwork quantum states[END_REF][START_REF] Bennewitz | Neural error mitigation of near-term quantum simulations[END_REF]. The fermionic modes are mapped into an interacting quantum spin model in this formalism. This can be achieved via the Jordan-Wigner [START_REF] Jordan | Über das paulische äquivalenzverbot[END_REF], the parity, or the Bravyi-Kitaev [START_REF] Bravyi | Fermionic quantum computation[END_REF] transformations. The reduction of the fermionic problem into a spin Hamiltonian makes it possible to exploit the success of NQS on spin systems. However, this approach suffers from the disadvantage of producing a spin Hamiltonian with non-local interactions. First quantization is therefore an attractive alternative as it preserves the locality of the physical interactions. In this case, the amplitudes of the variational state ψ θ (x) must be antisymmetric with respect to permutations of the particle indices. NQSbased parametrizations of fermionic wave functions borrow inspiration from traditional ansatzes like the Slater-Jastrow state, backflow correlations, and hidden particle representations.

The amplitudes of a Slater-Jastrow-inspired NQS consist of the product of a parameterized antisymmetric reference factor |φ (Slater determinant or PP) and a symmetric neural network factor M. The neural network is in charge of incorporating correlations on top of the reference wave function. First introduced in Ref. [START_REF] Nomura | Restricted Boltzmann machine learning for solving strongly correlated quantum systems[END_REF], a positive RBM was used as the correlation factor. Later works also implemented correlation factors that can alter the nodal structure of the reference state and respect the translational symmetry via the use of CNN with skip connections [START_REF] Stokes | Phases of two-dimensional spinless lattice fermions with first-quantized deep neural-network quantum states[END_REF]. It it noted that if |φ is a general Slater determinant of non-orthogonal orbitals, then the NQS is a universal parametrization in the lattice [START_REF] Moreno | Fermionic wave functions from neural-network constrained hidden states[END_REF]. Alternatively, neural networks have been used to parametrize the N -particle orbitals of an N -particle determinant in order to encode backflow correlations. This variational family has also been shown to be universal in the lattice [START_REF] Luo | Backflow transformations via neural networks for quantum many-body wave functions[END_REF].

Lastly, antisymmetric NQS has also been constructed using "hidden" additional fermionic degrees of freedom [START_REF] Moreno | Fermionic wave functions from neural-network constrained hidden states[END_REF]. In this case, the variational state is represented by a Slater determinant in the Hilbert space that is augmented by adding the "hidden" fermions. The Slater determinant in the augmented space is then projected into the physical Hilbert space, and this projection is parameterized by a neural network. The neural network parameters are optimized together with the singleparticle orbitals of the determinant. This ansatz explicitly contains the above Slater-Jastrow-inspired factorization, as well as a compact representation of configurationinteraction wave functions [START_REF] Moreno | Fermionic wave functions from neural-network constrained hidden states[END_REF].

Variational auxiliary-field quantum Monte Carlo

The variational auxiliary-field quantum Monte Carlo (VAFQMC) [START_REF] Sorella | The phase diagram of the Hubbard model by variational auxiliary field quantum Monte Carlo[END_REF] approach creates a variational ansatz for the ground state wave function of the Hubbard Hamiltonian Ĥ = K + V , by projections via an optimizable pseudo-Hamiltonian, using the formalism of AFQMC (see Section S3 E 1). In VAFQMC, a single Slater determinant |ψ MF is first constructed from an effective mean-field calculation with a set of variational parameters α 0 , such that ĤMF (α 0 ) |ψ MF = E 0 (α 0 ) |ψ MF . A variational ansatz is then constructed, by operating onto |ψ MF a projection operator:

|ψ τ = exp - τ 2 ĤMF (α) + V |ψ MF , (S63)
where α denotes a set of variational parameters. In Eq. (S63), ĤMF (α) replaces the kinetic part of the Hubbard Hamiltonian with a general quadratic operator of fermionic creation and annihilation operators, which can include, for instance, a d-wave BCS pairing field. In this work, ĤMF (α) was designed to give a |ψ MF describing AFM stripes with the expected wavelength [START_REF] Xu | Stripes and spin-density waves in the doped two-dimensional Hubbard model: Ground state phase diagram[END_REF]. The potential part V is kept as the original Hubbard on-site interaction. The parameter τ , the total imaginary time of the projection, is kept fixed and plays the role of an effective inverse temperature.

The projection in Eq. ( S63) is further broken up as

exp - τ 2 ĤMF (α) + V = n i=1 exp -t i ĤMF (α) exp -h i V × exp -t n+1 ĤMF (α) , ( S64
)
where n is the number of time steps given by n = max ([(U τ /0.4 -1)/5], 1) depending on the interaction strength U . The variable steps h i and t i are treated as additional variational parameters to optimize in order to minimize Trotter errors [START_REF] Beach | Making Trotters sprint: A variational imaginary time ansatz for quantum many-body systems[END_REF]. However, here we introduce a simple functional form for the non-uniform time steps which depends on a single variational parameter ∆τ (see Ref. [START_REF] Sorella | The phase diagram of the Hubbard model by variational auxiliary field quantum Monte Carlo[END_REF] for details). The variational ansatz |ψ τ in Eq. (S63) can thus be equivalently denoted as |ψ n , yielding a variational energy E n := ψ n | Ĥ|ψ n / ψ n |ψ n . Now following similar procedures to AFQMC, we recast the variational energy as

E n = σσ ψ MF | Û † n (σ ) Ĥ Ûn (σ)|ψ MF σσ ψ MF | Û † n (σ ) Ûn (σ)|ψ MF , (S65)
where

Ûn (σ) := n i=1 exp -t i ĤMF (α) exp λ i j σ j,i mj × exp -t n+1 ĤMF (α) , (S66)
with the vector σ := {σ ji } ∈ {±1} Ns×n denoting the collection of auxiliary fields arising from the discrete Hubbard-Stratonovich transformation [START_REF] Hirsch | Two-dimensional hubbard model: Numerical simulation study[END_REF]. The spin operator for each site j is defined as mj = nj,↑ -nj,↓ , and the parameter λ i is given by cosh λ i = exp(U h i ).

Eq. ( S65) is now rewritten in a form suitable for Monte Carlo:

E n = σσ |W n (σ , σ)|e n (σ , σ)S n (σ , σ) σσ |W n (σ , σ)|S n (σ , σ) , ( S67
)
where

W n (σ , σ) := ψ MF | Û † n (σ ) Ûn (σ)|ψ MF , (S68) S n (σ , σ) := W n (σ , σ) |W n (σ , σ)| , (S69) e n (σ , σ) := ψ MF | Û † n (σ ) Ĥ Ûn (σ)|ψ MF W n (σ , σ) (S70)
are the weight factor, the phase factor, and the local energy respectively. For each set of variational parameters, we estimate the energy by E n = S n e n Wn / S n Wn , where • Wn indicates the average of a random variable with respect to the Monte Carlo samples from the probability distribution |W n |.

The optimization of the variational parameters is then carried out by generalizing techniques from standard variational Monte Carlo [START_REF] Becca | Quantum Monte Carlo Approaches for Correlated Systems[END_REF][START_REF] Sorella | Green function Monte Carlo with stochastic reconfiguration[END_REF] and from machine learning [START_REF] Amari | Natural gradient works efficiently in learning[END_REF]. Minimizing E n (α) with respect to α requires the energy derivatives. Assuming there are 2p variational parameters α 1 , α 2 , . . . , α 2p plus an extra parameter of the minimum time step α 2p+1 = ∆τ , we can compute the derivatives by

∂E n ∂α j = S n ∂en ∂αj + (e n -E n )O j Wn S n Wn , (S71)
where O j := ∂ ln Wn ∂αj . The complex derivatives ∂en ∂αj and O j are obtained by automatic differentiation [START_REF] Griewank | Evaluating Derivatives (Society for Industrial and Applied Mathematics[END_REF]. The final optimized E n defined within the variational ansatz is dependent on the parameter τ , and provides an upper bound to the true ground-state energy of the given Hamiltonian.

D. Parameterized quantum circuits and variational quantum eigensolver

Variational states can be represented by parameterized quantum circuits (PQC) [START_REF] Cerezo | Variational quantum algorithms[END_REF]. A PQC has parameterized single and two-qubit gates, e.g. rotation gates, as well as multi-qubit entangling gates [START_REF] Nielsen | Quantum Computation and Quantum Information: 10th Anniversary Edition[END_REF]. Different hardware architectures come with different sets of native gates. The set of variational parameters θ is optimized using the variational quantum eigensolver (VQE) algorithm [START_REF] Peruzzo | A variational eigenvalue solver on a photonic quantum processor[END_REF]. When using quantum hardware, the expectation value of the Hamiltonian and the energy gradient components must be measured through repeated wave function collapses, thus introducing a sampling noise similar to the one in VMC simulations. We only consider noiseless emulations of the quantum circuits, i.e., neglecting this measurement sampling noise, as well as hardware noise due to decoherence. We do this because our goal is to test whether the relation between the V-score and the relative energy accuracy also holds in the quantum computing setting.

We employ several types of PQC ansatzes and optimization methods. We first tackle the TFIM at the criticality Γ = 1, defined on a finite chain with L sites.

The circuit is made of a series of d blocks built from single-qubit rotations ÛR (θ k ), interlayered with entangler blocks Ûent , that spans the required length of the qubit register, with k = 1, . . . , d + 1. This is made of a ladder of CNOT (also known as CX) gates with linear connectivity, such that qubit q i is target of qubit q i-1 and controls qubit q i+1 for i = 1, . . . , L -2. Since the singlequbit rotations are all local operations, ÛR (θ k ) can be written as a tensor product of rotations on single qubits:

ÛR (θ k ) = L-1 i=0 Ry (θ k qi ), ( S72
)
where Ry (θ k qi ) is a rotation around the y-axis on the Bloch sphere of qubit q i . The full unitary circuit operation is described by ÛR-CX (θ) := ÛR (θ d+1 )

1 i=d Ûent ÛR (θ i ), ( S73
)
and the final parameterized state is

|ψ(θ) = ÛR-CX (θ) |0 ⊗L . (S74)
The accuracy of the calculation is controlled by the circuit depth d, and the total number of variational parameters is L(d + 1). An alternative is to use physically motivated ansatzes, such as the Hamiltonian variational (HV) ansatz [START_REF] Wecker | Progress towards practical quantum variational algorithms[END_REF]. The unitary operator defining the HV ansatz is made of d blocks, and each block is a product of operators Ûα = exp(iθ k α Ĥα ), with α = 1, . . . , indexing the noncommuting terms of the Hamiltonian. For TFIM we only need = 2. In this case, the full unitary operator is

ÛHV (θ) := 1 i=d Û2 (θ i 2 ) Û1 (θ i 1 ), (S75)
which can be efficiently decomposed using one-and twoqubit quantum gates, and the final parameterized state is

|ψ(θ) = ÛHV (θ) |0 + |1 √ 2 ⊗L , (S76)
where the initial non-entangled state can be obtained from |0 ⊗L by placing one Hadamard gate on each qubit. The total number of parameters is d.

We use both the heuristic R-CX ansatz and the physically motivated HV ansatz on the TFIM. We have obtained a family of optimized trial states that depend on the circuit depth d, as shown in Fig. S1, confirming the linear scaling of V-score versus energy relative error proposed in the main text.

For the Heisenberg and the J 1 -J 2 models, we employ the symmetry-enhanced architecture introduced in Ref. [START_REF] Seki | Symmetry-adapted variational quantum eigensolver[END_REF]. This PQC is more sophisticated compared to the previous ones, which are included mainly for demonstration purposes, and allows one to enforce translational, point-group, and SU(2) symmetries of the variational wave function in the device noise-free case. This circuit has shown outstanding performance in representing the ground state with only a few variational parameters when applied to 2D frustrated magnets [START_REF] Astrakhantsev | Algorithmic phases in variational quantum ground-state preparation[END_REF].

The initial state to be prepared is the dimerized state, splitting the L site indices into L/2 arbitrary pairs:

|ψ D := 0≤i<L/2 1 √ 2 ↑ 2i ↓ 2i+1 -↓ 2i ↑ 2i+1 , (S77)
which is a product of L/2 dimers on the selected pairs. This state is manifestly the SU(2)-singlet. To access other total spin quantum numbers, one should replace one or more singlets in the product Eq. (S77) with triplets

1 √ 2 ↑ 2i ↓ 2i+1 + ↓ 2i ↑ 2i+1 . We then introduce the SWAP operator Pij := 1 2 Ŝi • Ŝj + 1 , (S78)
which exchanges spin states between sites i and j. We note that the SWAP operator commutes with the total spin operator Ŝ2 , therefore the variational ansatz

|ψ(θ) = α e iθα Piαjα |ψ D (S79)
has a well-defined total spin quantum number. The exponential SWAP (eSWAP) operator e iθα Piαjα can be efficiently implemented on a quantum computer, as the two-qubit SWAP gate can be compiled into single-qubit gates and CNOTs, or is a native gate in architectures alternative to superconductive qubits.

Let us now suppose that, in addition to the SU(2) symmetry, the system also has translational or point-group symmetries, which can be all represented as qubit permutations. In Ref. [START_REF] Seki | Symmetry-adapted variational quantum eigensolver[END_REF], the authors provide a way that allows one to effectively project the wave function onto any irreducible representation subspace of the spatial symmetry. To this end, let us consider the spatial symmetry projector

P := 1 |G| g∈G χ g ĝ, ( S80
)
where G is the full spatial symmetry group consisting of the elementary unitary permutations ĝ, and χ g are the characters depending on the desired projection quantum number (irreducible representation). The projected wave function |ψ

P (θ) = P √ N (θ)
|ψ(θ) is normalized with N (θ) := ψ(θ)| P |ψ(θ) , since P 2 = P . This wave function is optimized using the natural gradient approach discussed in Section S3 C. The energy gradient is preconditioned using the metric tensor, which mimics imaginary-time evolution in the allowed variational subspace. For this symmetry-enhanced ansatz, the energy gradient reads

∂ i E(θ) = 2 Re ψ(θ)| Ĥ P |∂ i ψ(θ) N (θ) -A i (θ) E(θ) , (
G(θ) ij := ∂ i ψ(θ)| P |∂ j ψ(θ) N (θ) -A * i (θ)A j (θ). ( S82
)
The parameter update reads

θ ← θ -η j Re G(θ) -1 ij ∂ j E(θ) . (S83)
The matrix elements required to construct these objects can be measured using the Hadamard test rule [START_REF] Seki | Symmetry-adapted variational quantum eigensolver[END_REF][START_REF] Astrakhantsev | Algorithmic phases in variational quantum ground-state preparation[END_REF]. During the simulation of the VQE optimization on a classical computer, we either measure these matrix elements using N s circuit shots or compute them exactly. In the former case, the metric tensor obtained with sampling should be regularized to make the matrix inversion well defined. Here, we employ the regularization G reg = √ GG + β1 as suggested in Ref. [START_REF] Gacon | Simultaneous perturbation stochastic approximation of the quantum fisher information[END_REF]. Finally, it is also important to check that the V-score can be efficiently computed with VQE. The current way to estimate the energy with VQE is by calculating the weighted sum of the expectation values of all Pauli operators that compose the Hamiltonian. For instance, the TFIM is made of L local terms, yet most of them can be measured simultaneously [START_REF] Kandala | Hardwareefficient variational quantum eigensolver for small molecules and quantum magnets[END_REF]. It turns out that only two bases are needed: the computational basis, and the rotated {σ x } ⊗L basis. A similar argument applies for Ĥ2 : While an upper bound for the number of terms is L 2 , it is possible to check that the number of groups of Pauli operators that can be measured simultaneously grows sub-linearly with the system size, in both one and two dimensions.

The 2D Hubbard model introduces the additional complication of the fermion-to-qubit mapping. In more than one dimension, the Jordan-Wigner mapping generally transforms a two-local fermionic operator, i.e. the hopping term, into a non-local qubit operator [START_REF] Peruzzo | A variational eigenvalue solver on a photonic quantum processor[END_REF]. However, numerical tests show that the number of bases, thus the number of measurements, grows only sub-quadratically with the system size.

E. Quantum Monte Carlo

In the benchmarks, we include some quantum Monte Carlo (QMC) methods. Although they are not strictly variational and the V-score is not applicable to them, they produce numerically exact results in some cases, while in other cases highly accurate energies and the bias has been certified to be small. These cases are discussed below, with the methods described in detail. The QMC results are used to assess the V-scores of variational methods when ED is not practical.

Auxiliary-field quantum Monte Carlo

There exist two different auxiliary-field Monte Carlo algorithms, based on one hand the finite temperature grand canonical ensemble [START_REF] Blankenbecler | Monte Carlo calculations of coupled boson-fermion systems. I[END_REF], and on the other hand the ground state canonical ensemble [START_REF] Sorella | A novel technique for the simulation of interacting fermion systems[END_REF][START_REF] Imada | Numerical studies on the Hubbard model and the t-J model in one-and twodimensions[END_REF]. The auxiliaryfield quantum Monte Carlo (AFQMC) method (for an overview, see e.g. Ref. [START_REF] Zhang | Auxiliary-field quantum Monte Carlo at zeroand finite-temperature[END_REF]) used in this work is based on the latter. It filters out the ground state from an initial state by an imaginary-time propagation: |ψ G = lim τ →∞ e -τ Ĥ |ψ I , where |ψ I and |ψ G are the initial and ground state respectively, τ the imaginary time, and Ĥ the many-body Hamiltonian. The initial state must satisfy ψ I |ψ G = 0 but can be otherwise arbitrary. We have typically taken it from a mean-field calculation [START_REF] Qin | Benchmark study of the twodimensional hubbard model with auxiliary-field quantum Monte Carlo method[END_REF][START_REF] Purwanto | Eliminating spin contamination in auxiliary-field quantum Monte Carlo: Realistic potential energy curve of F2[END_REF].

The projection is achieved by first discretizing the imaginary time into m small time steps ∆τ : e -τ Ĥ = (e -∆τ Ĥ ) m , and then applying Trotter-Suzuki breakup [START_REF] Trotter | On the product of semi-groups of operators[END_REF][START_REF] Suzuki | Generalized Trotter's formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems[END_REF] to each time step: e -∆τ Ĥ = e -∆τ K/2 e -∆τ V e -∆τ K/2 + O ∆τ 3 . (S84)

Here K is the kinetic part consisting of one-body operators, and V the interacting part containing two-body operators. For all systems computed in this work, we have either extrapolated ∆τ to zero or set it to a fixed value (typically ∆τ = 0.01) and verified that the Trotter error is well within our statistical error.

Hubbard-Stratonovich (HS) transformation is applied to rewrite the interacting part into a one-body form coupled with auxiliary fields (AFs) {x i }. In the Hubbard model, a spin decomposition is applied to the U > 0 case: e -∆τ U ni↑ ni↓ = e -∆τ U (n i↑ +n i↓ )/2 xi=±1 1 2 e γxi(n i↑ -n i↓ ) , (S85) where cosh γ = exp(∆τ U/2). For the U < 0 cases, we apply a charge decomposition form of the HS transformation. A systematic study of the effect of the different HS transformations can be found in Ref. [START_REF] Shi | Symmetry in auxiliary-field quantum Monte Carlo calculations[END_REF].

The projection is carried out by evaluating the propagator as an integral using the Monte Carlo method: e -∆τ Ĥ = p(x) B(x)dx, (S86) where x := {x i } = {x 1 , x 2 , . . . , x Ns } for a lattice with N s sites, and p(x) is a probability distribution in AF space, which is a uniform function in the discrete HS transformation above. The propagator B(x) now only contains one-body operators: B(x) := e -∆τ K/2b (x)e -∆τ K/2 , where b(x) := Ns i=1 bi (x i ) is the product of the one-body operators transformed from the interacting part, i.e., the right-hand side of Eq. (S85).

In this work, two different ground state AFQMC methods are used. For sign-problem-free systems, e.g., the repulsive Hubbard model (U > 0) at half-filling or the spin-balanced attractive Hubbard model (U < 0) [START_REF] Hirsch | Two-dimensional hubbard model: Numerical simulation study[END_REF][START_REF] Loh | Sign problem in the numerical simulation of many-electron systems[END_REF], numerically exact results are given by AFQMC using a generalized Metropolis algorithm with force bias [START_REF] Shi | Ground-state properties of strongly interacting Fermi gases in two dimensions[END_REF][START_REF] Shi | Infinite variance in fermion quantum Monte Carlo calculations[END_REF]. On the other hand, for a doped Hubbard model with U > 0, AFQMC is applied with a constraint path (CP) approximation [START_REF] Zhang | Constrained path Monte Carlo method for fermion ground states[END_REF][START_REF] Nguyen | CPMC-Lab: A matlab package for constrained path Monte Carlo calculations[END_REF] to control the sign problem. The two ground state AFQMC methods are described in the following two subsections, respectively.

Sign-problem-free Hubbard model: exact AFQMC with generalized Metropolis algorithm

For the half-filled repulsive Hubbard model on a square lattice, as well as the spin-balanced attractive Hubbard model, the sign problem is absent in AFQMC because of symmetry. In these cases, exact calculations are performed. (Note that in almost all sign-problem-free calculations, the standard approach has an infinite variance problem which must be properly taken care of [START_REF] Shi | Infinite variance in fermion quantum Monte Carlo calculations[END_REF].) We study these systems using the AFQMC with a generalized Metropolis algorithm [START_REF] Shi | Ground-state properties of strongly interacting Fermi gases in two dimensions[END_REF]. The ground-state energy is given by Ĥ = ψ L | Ĥ|ψ R / ψ L |ψ R , where ψ L | := ψ I | exp -τ L Ĥ and |ψ R := exp -τ R Ĥ |ψ I . Since the Hamiltonian commutes with propagators, energy measurements can be inserted between any two time steps, i.e., with any combination of τ L + τ R = τ , pro-vided τ is sufficiently larger than the equilibration time, τ > τ eq , to reach the ground state from |ψ I .

To illustrate the sampling process, we rewrite the energy expectation as a path integral in AF space: ) |ψ I . These are single Slater determinants if the initial state is chosen as a single Slater determinant. In the case of a multi-determinant |ψ I , the different terms in the linear combination can be sampled. The number M = τ /∆τ is the total number of time slices, and M L and M R correspond to τ L and τ R respectively. The variables X := {x (1) , x (2) , . . . , x (M ) } form a N s × M -dimensional vector in AF space. The probability distribution P (X) := M m=1 p(x (m) ). We adopt a force bias method [START_REF] Shi | Ground-state properties of strongly interacting Fermi gases in two dimensions[END_REF] instead of the usual single-site heat bath update, which helps reduce the autocorrelation time in the Monte Carlo sampling. In this method, we update a cluster of AFs at each time slice n simultaneously. The size of the cluster N c can be as large as N s and can be tuned according to the acceptance ratio. Below we sketch the algorithm for the spin decomposition; generalization to other HS transformations is straightforward [START_REF] Zhang | Auxiliary-field quantum Monte Carlo at zeroand finite-temperature[END_REF][START_REF] Shi | Ground-state properties of strongly interacting Fermi gases in two dimensions[END_REF]. Each AF in the cluster is proposed a new value according to i ). The optimal choice of force bias is given by niσ := φ L |n iσ |φ R / φ L |φ R , where the σ and i are the spin and lattice site indices, and leads to a P that approximates the target probability to O( √ ∆τ ).

Ĥ = φ L | Ĥ|φ R φ L |φ R P (X) φ L |φ R dX P (X) φ L |φ R dX , ( S87

Doped repulsive Hubbard model: constrained path AFQMC

All the AFQMC results on the repulsive Hubbard model away from half-filling are obtained with CP-AFQMC [START_REF] Zhang | Auxiliary-field quantum Monte Carlo at zeroand finite-temperature[END_REF][START_REF] Zhang | Constrained path Monte Carlo method for fermion ground states[END_REF][START_REF] Zhang | Constrained path quantum Monte Carlo method for fermion ground states[END_REF]. CP-AFQMC is built on an openended random walk, with an indefinite value of τ . The sign problem is removed by introducing a trial state |ψ T to guide and constrain the random walk. In this work, |ψ T is obtained from noninteracting calculations for closed-shell systems and Hartree-Fock calculations for open-shell systems. In the latter cases, the CP error is further reduced by applying a self-consistent constraint [START_REF] Qin | Coupling quantum Monte Carlo and independent-particle calculations: Selfconsistent constraint for the sign problem based on the density or the density matrix[END_REF] which couples to a generalized Hartree Fock (GHF) calculation with an effective U eff determined via the self-consistency [START_REF] Qin | Benchmark study of the twodimensional hubbard model with auxiliary-field quantum Monte Carlo method[END_REF]. Several results are also provided from constraint release [START_REF] Shi | Symmetry in auxiliary-field quantum Monte Carlo calculations[END_REF] which are essentially exact, as indicated in the main text.

Different from the Metropolis approach discussed above, CP-AFQMC is performed via an open-ended branching random walk, in which a population of walkers { |φ (n) k , w (n) k } are propagated following the time evolution |ψ (n+1) = e -τ Ĥ |ψ (n) . These walkers sample the many-body wave function in the sense that |ψ

(n) ∝ k w (n) |φ (n) k / ψ T |φ (n) k
. The ground-state energy is then given by Ĥ = ψ T | Ĥ|ψ G / ψ T |ψ G . After the random walk has equilibrated (i.e., after a sufficient number of steps n eq = τ eq /∆τ ), the walkers will sample |ψ G , and all subsequent steps n > n eq can be used to compute the ground-state energy:

Ĥ = k,n ψ T Ĥ φ (n) k ψ T φ (n) k w (n) k k,n w (n) k . ( S89
)
In the random walks, importance sampling is introduced which amounts to sampling new AF to advance the walker by a modified probability density:

|φ (n+1) ← p(x) B(x)dx |φ (n) , (S90)
where p(x) ∝ ψ T |φ (n+1) / ψ T |φ (n) builds in the knowledge from |ψ T to improve the sampling efficiency [START_REF] Zhang | Auxiliary-field quantum Monte Carlo at zeroand finite-temperature[END_REF]. The actual form of p(x) contains force bias terms similar to how we formulated the generalized Metropolis algorithm above.

The importance sampling transformation also automatically imposes a constraint [START_REF] Zhang | Auxiliary-field quantum Monte Carlo at zeroand finite-temperature[END_REF]. As ∆τ approaches zero, the random walkers will not cross the surface defined by ψ T |φ (n) k = 0, which separates two degenerate regions in AF space, each of which is overcomplete and can fully represent |ψ G . Constraining the random walks in one region of the Slater determinant space (or equivalently, the AF space) is an exact condition if |ψ T = |ψ G . The CP approximation uses an approximate |ψ T to impose this condition, which leads to a systematic error. The ground-state energy computed by the mixed estimate of Eq. (S89) is therefore not variational [START_REF] Carlson | Issues and observations on applications of the constrainedpath Monte Carlo method to many-fermion systems[END_REF]. The CP-AFQMC method has been extensively benchmarked, and the CP error was shown to be small for Hubbard-like systems (see, e.g., Refs. [START_REF] Qin | Benchmark study of the twodimensional hubbard model with auxiliary-field quantum Monte Carlo method[END_REF][START_REF] Shi | Symmetry in auxiliary-field quantum Monte Carlo calculations[END_REF][START_REF] Leblanc | Solutions of the two-dimensional Hubbard model: Benchmarks and results from a wide range of numerical algorithms[END_REF]), and the results provided here for the doped Hubbard model are expected to be very accurate, with the relative error below a few tenths of a percent.

As mentioned, the energies for the doped Hubbard model were computed using either (i) constraint release (taken from Ref. [START_REF] Shi | Symmetry in auxiliary-field quantum Monte Carlo calculations[END_REF], essentially exact), (ii) with freeelectron |ψ T (for closed-shell systems), or (iii) with a self-consistent constraint [START_REF] Qin | Coupling quantum Monte Carlo and independent-particle calculations: Selfconsistent constraint for the sign problem based on the density or the density matrix[END_REF]. In the data included in this paper, we have indicated how each energy was obtained and, in the case of (iii), including the final U eff in the GHF after convergence of the self-consistent iteration between GHF and AFQMC.

Continuous-time quantum Monte Carlo

The continuous-time quantum Monte Carlo [START_REF] Iazzi | Efficient continuous-time quantum Monte Carlo algorithm for fermionic lattice models[END_REF] For the t-V model, the initial state |ψ I is chosen as a single Slater determinant, and Ĥ0 and Ĥ1 are the noninteraction and the interaction terms respectively. Here, we write Ĥ1 = V 4 i,j e iπ(ni+nj ) to ensure that each term in the interaction expansion is evaluated as a determinant. Simulation of the spinless t-V model on the bipartite lattice at half-filling with repulsive interaction is free from the fermion sign problem [START_REF] Wang | Efficient continuous-time quantum Monte Carlo method for the ground state of correlated fermions[END_REF]. Therefore, the results are free from systematic Trotter or constrained path errors. The computational time complexity of this method scales as O(τ V N 3 s ).

Figure 2 .

 2 Figure2. Validation of V-score against exact results. We compare V-scores versus energy relative errors on various strongly correlated models for which exact results (ED or QMC) are available. The black dashed line is a least square fit of log(energy rel. err.) = log(V-score) + C, where C = -1.80 ± 0.08. The inset focuses on PQC results run on classical hardware (no shot noise included).

  S8) where a locally interacting impurity Ĥloc is coupled to a non-interacting bath Ĥbath . The indices {α} are a collection of quantum numbers denoting the impurity (or the l-th bath site) degrees of freedom of the fermionic creation d † (or ĉ † l ) and annihilation d (or ĉl ) operators, and N b is the number of bath sites per spin-orbital. The bath parameters { l , ν l } are connected to the hybridization function as ∆(ω) = N b l=1 ν l ν l † ω-l , and can be obtained by discretizing the hybridization function ∆(ω) on the real frequency axis into N b equidistant intervals {I l } of size ∆ω as l αα = min I l

2

 2 S15)with m ∈ {1, 2, 3} being the orbital index, and ĤDD , ĤSF , and ĤPH denoting the density-density, the spin-flip, and the pair-hopping interactions respectively.The following models representing a collection of typical solutions in practical DMFT calculations are considered: (SB-Imp) single-band Anderson impurity model with a semielliptic spectral function, i.e., -1 π Im ∆(ω) = with D being the half-bandwidth and U = D; (SB-DMFT-MT-HF) DMFT metal solution of the single band Hubbard model on the Bethe lattice with U = 2D at half-filling n = 1 and (SB-DMFT-MT-AHF) doped case n = 0.8; (SB-DMFT-MI-HF) DMFT Mott-insulator solution of the single band Hubbard model on the Bethe lattice with U = 4D at halffilling n = 1; three-band models with Kanamori interaction U = 2.3 eV and J = 0.4 eV that are based on the material-realistic DMFT solutions of the archetypal Hund's metal Sr 2 RuO 4 in the t 2g subspace (TB-DMFT-SOC) with and (TB-DMFT) without spin-orbit coupling.

  S49) with ρb = |e |o o| e| and ρi = |i i|. For an impurity model with N e electrons and N o occupied bath sites, we evaluate E ∞ using the stochastic sampling method presented in Sec. S2 B 2 with the constraint that each sampled state |i has a fixed particle number of N b -N o .

Figure S1 .

 S1 Figure S1. Comparison of V-scores of VQE ansatzes versus energy relative errors on a 10 sites TFIM. The black dashed line is the least square fit of all variational methods shown in Fig. 2 in the main text. Each data point corresponds to a full VQE optimization using different circuit depths d. For the HV ansatz, we report d = 8, 12, 16, 20, 24, 26, while for R-CX we report d = 4, 6, 8, 10, 12. Longer circuits provide systematically better V-scores and energy relative errors.

  S81) where A i (θ) := 1 N (θ) ψ(θ)| P |∂ i ψ(θ) is the connection. The metric tensor is defined as

  )whereφ L | := ψ I | M L m=1 B(x (M -m+1) ) and |φ R := M R m=1 B(x (m)

P

  (x i ) := e γxi(n i↑ -n i↓ ) xi=±1 e γxi(n i↑ -n i↓ ) , (S88) thus giving a probability density for a new candidate cluster P(x(n)

1 i=k- 1 e

 11 algorithm carries out a diagrammatic expansion of the imaginary-time projection operator e -τ Ĥ and samples interaction expansion termsψ I |e -τ Ĥ |ψ I = ψ I |e -(τ -τ k ) Ĥ0 Ĥ1 -(τi+1-τi) Ĥ0 Ĥ1 e -τ1 Ĥ0 |ψ I , (S91)where Ĥ = Ĥ0 + Ĥ1 . Ground state physical observables are measured asÔ = lim τ →∞ ψ I |e -τ Ĥ/2 Ô e -τ Ĥ/2 |ψ I ψ I |e -τ Ĥ |ψ I . (S92)

ACKNOWLEDGMENTS

We acknowledge discussions with A. Sandvik and M. Stoudenmire. This research was supported by the NCCR MARVEL, a National Centre of Competence in Research, funded by the Swiss National Science Foundation (grant number 205602). The Flatiron Institute is a division of the Simons Foundation. A.W. acknowledges support by the DFG through the Emmy Noether programme (WI 5899/1-1). N.A. and T.N. acknowledge support from the European Unions Horizon 2020 research and innovation program (ERC-StG-Neupert-757867-PARATOP). M.I., Y.N., R.P., and M.S. thank the support by MEXT as "Program for Promoting Researches on the Supercomputer Fugaku" (Basic Science for Emergence and Functionality in Quantum Matter -Innovative Strongly Correlated Electron Science by Integration of Fugaku and Frontier Experiments, JP MXP1020200104) together with computational resources of supercomputer Fugaku provided by the RIKEN Center for Computational Science (Project ID: hp210163 and hp220166). M.I. acknowledges the support of Grant-in-Aid for Transformative Research Areas, No.22A202 from MEXT. J.C. acknowledges support from the Natural Sciences and Engineering Research Council (NSERC), the Shared Hierarchical Academic Research Computing Network (SHARCNET), Compute Canada, and the Canadian Institute for Advanced Research (CIFAR) AI chair program. The dataset of benchmark results is available at https://github.com/varbench/varbench