Fermionic Wave Functions from Neural-Network Constrained Hidden States

We introduce a systematically improvable family of variational wave functions for the simulation of strongly correlated fermionic systems. This family consists of Slater determinants in an augmented Hilbert space involving"hidden"additional fermionic degrees of freedom. These determinants are projected onto the physical Hilbert space through a constraint which is optimized, together with the single-particle orbitals, using a neural network parametrization. This construction draws inspiration from the success of hidden particle representations but overcomes the limitations associated with the mean-field treatment of the constraint often used in this context. Our construction provides an extremely expressive family of wave functions, which is proven to be universal. We apply this construction to the ground state properties of the Hubbard model on the square lattice, achieving levels of accuracy which are competitive with state-of-the-art variational methods.

Many-body quantum systems are computationally challenging because of the exponential dependence of the size of the Hilbert space on the number of particles.Variational approaches address this problem by considering a class of wave functions depending on a set of parameters over which an optimization is performed.In this way, the computationally intractable search over the full Hilbert space is reduced to a search over a submanifold of dimension merely polynomial in the number of particles.Variational approaches have proven successful in providing qualitative and quantitative insights into the nature of the ground state and the low-energy excited states of a number of interacting quantum systems.For example, in the case of spin systems with arbitrary pairwise interactions, it has been proven 2,3 that the ratio between the energy of optimized mean-field states and the true ground state energy approaches a finite constant in the limit of large system size.The subsequent development of systematically improvable variational wave functions has led to quantitative agreement with exact energies of one-dimensional systems using matrix product states, and recently also in two dimensions using neural-network and tensor-network states 4,5 .
The remarkable success of variational states in the description of quantum spin systems unfortunately does not have a parallel in correlated systems of fermions, however.It is known, for example, that the natural mean-field analogue of direct-product states, the so-called Slater determinant (SD) states, fail to even qualitatively describe the thermodynamic limit of Fermi-Hubbard type Hamiltonians 3 and the devel-opment of systematically improvable neural-network-based trial wave functions is currently an active field of research both in second quantization [6][7][8] , and first quantization 9-15 .In the latter approach, the wave function amplitudes must be anti-symmetric functions of the particle configurations, while being able to capture correlations beyond the singleparticle Slater determinants.This is typically achieved either by considering determinants of multi-particle orbitals 10,12,14 (backflow transformations), or by Slater determinants of single-particle orbitals multiplied by a neural-network Jastrow factor that depends on the lattice occupations 9,11 .Despite being universal in the lattice, the Slater neural-network Jastrow wave functions seem to struggle to get competitive energies in the strong coupling regime.
The Hubbard model on the square lattice has been the subject of intense theoretical scrutiny, and constitutes the most iconic 'simple' model of an interacting quantum system.Despite this simplicity, a full computational solution is still to be achieved.For this model, as well as related lattice models of interacting fermions such as the t-J and Kondo lattice models, significant insight has been obtained using hidden particle approaches 1 .Although a number of different formulations are available [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] all such approaches share a basic concept which consists in augmenting the physical Hilbert space by auxiliary degrees of freedom and subsequently performing a projection back to the subspace of physical states.This projection can be regarded as a constraint that selects the representative states in the augmented space that are identified with the basis of the physical Hilbert space.In many cases, a mean-field saddle point approximation is applied both to the auxiliary particle hamiltonian and to the treatment of the constraint, which is implemented with static and uniform Lagrange multipliers.This mean-field approximation is uncontrolled in general, except when the saddle-point is associated with the limit of large number of flavors 18,20 .Even in those cases, going beyond the saddle point level is challenging and no systematic improvements beyond the meanfield variational wave functions are available, especially in view of the approximate treatment of the constraint.
In this article, we draw inspiration from hidden particle approaches to construct a systematically improvable family of variational fermionic wave functions.These states are obtained as the exact projection of Slater determinant states in a Hilbert space augmented by hidden fermionic degrees of freedom.One of the major novelties of the proposed method is that the constraint is parametrized by neural networks, giving rise to an extremely flexible family of wave function ansätze.The constraint is optimized together with the orbitals in the enlarged Hilbert space with the goal of minimizing the energy.The expressive power of this new class of wave functions is demonstrated in a variational Monte Carlo (VMC) setting, obtaining an accuracy which is competitive with the state-of-the-art for the ground-state properties of the Hubbard model in the square and rectangular lattices.
The paper is structured as follows: we begin (Sec.I) by introducing the Hamiltonian and the physical degrees of freedom of the problem.In Sec.II, we introduce the fundamentals of the hidden fermion representation, describe the Slater determinant in the augmented space together with the fully parametrized constraint function and prove the universality of this representation.This section also contains details on the VMC implementation.In Sec.III, we present groundstate energy benchmarks for the Hubbard model with increasingly large system sizes, and demonstrate that we can stabilize competing orders of charge and spin stripes for the Hubbard model on rectangular lattice geometries.

I. BACKGROUND: STATES AND HAMILTONIAN
In this paper we develop a general technique for approximating the ground state of interacting fermionic hamiltonians with discrete degrees of freedom-as defined for example by discrete orbitals or spatial coordinates.As a specific application, we focus here on the Fermi-Hubbard model, whose Hamiltonian reads (1) where the binary index σ ∈ {↑, ↓} labels two species of fermionic modes satisfying the canonical anti-commutation relations, The fermionic modes ĉiσ are the physical (electronic) degrees of freedom (DOF).The fermion dynamics is described by the lattice with V sites defined by the non-zero entries of the t ij hopping matrix, as well as by the onsite coulomb repulsion U i .In the following, we exclusively focus on the square and rectangular lattices with uniform hopping (t ij = 1) and onsite repulsion, leaving more general geometries to future studies.
< l a t e x i t s h a 1 _ b a s e 6 4 = " g k s 3 r Y g 6 a y e Q e v W + q m a P B Q E s 2 g Y = " > A A A B / H i c b V B N S 8 N A E J 3 U r 1 q / o j 1 6 W S y C p 5 K I q M e i F 4 8 V 7 A c 0 o W w 2 m 3 b p Z h N 2 N 0 I I 9 a 9 4 8 a C I V 3 + I N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q c q Z 0 o 7 z b V X W 1 j c 2 t 6 r b t Z 3 d v f 0 D + / C o q 5 J M E t o h C U 9 k P 8 C K c i Z o R z P N a T + V F M c B p 7 1 g c j v z e 4 9 U K p a I B 5 2 n 1 I / x S L C I E a y N N L T r A n k x T p V O k K c Z D 2 k h p k O 7 4 T S d O d A q c U v S g B L t o f 3 l h Q n J Y i o 0 4 V i p g e u k 2 i + w 1 I x w O q 1 5 m a I p J h M 8 o g N D B Y 6 p 8 o v 5 8 V N 0 a p Q Q R Y k 0 J T S a q 7 8 n C h w r l c e B 6 Y y x H q t l b y b + 5 w 0 y H V 3 7 B R N p p q k g i 0 V R x p F 5 d Z Y E C p m k R P P c E E w k M 7 c i M s Y S E 2 3 y q p k Q 3 O W X V 0 n 3 v O l e N i / u L x q t m z K O K h z D C Z y B C 1 f Q g j t o Q w c I 5 P A M r / B m P V k v 1 r v 1 s W i t W O V M H f 7 A + v w B z e C U 5 A = = < / l a t e x i t > n 7 !ñ ( M > 0) < l a t e x i t s h a 1 _ b a s e 6 4 = " K f t S Q O X v N c 2 x a Q R g 4 2 E s J V M 3 3 w 0 = " > A A A C N H i c b V B N S x x B E O 1 R E 3 U T k 1 W P u T R Z A u a y z I h E Q Q Q x E A L J Q S G r w s 5 m q e m p n W 2 2 u 2 f o r h G X Y X 6 U l / y Q X C S Q Q y R 4 9 T f Y u + 7 B r w f d P N 6 r o q p e U i j p K A z / B H P z C y 9 e L i 4 t N 1 6 9 X n n z t r m 6 d u z y 0 g r s i F z l 9 j Q B h 0 o a 7 J A k h a e F R d C J w p N k 9 H n i n 5 y h d T I 3 P 2 h c Y E 9 D Z u R A C i A v 9 Z v f Y s J z s r r 6 k o s R d w U I n P z G Y M q T M a 9 5 P A S q o P 4 Z p 5 B l a P u V j J 3 M N H h n d y M m q V K s v t d 8 j 4 c f + 8 1 W 2 A 6 n 4 E 9 J N C M t N s N h v / k 7 T n N R a j Q k F D j X j c K C e h V Y k k J h 3 Y h L h 3 6 j E W T Y 9 d S A R t e r p k f X / I N X U j 7 I r X + G + F S 9 3 1 G B d m 6 s E 1 + p g Y b u s T c R n / O 6 J Q 1 2 e p U 0 R U l o x N 2 g Q a k 4 5 X y S I E + l R U F q 7 A k I K / 2 u X A z B g i C f c 8 O H E D 0 + + S k 5 3 m x H n 9 p b R 1 u t / Y N Z H E v s H X v P N l j E t t k + + 8 o O W Y c J d s E u 2 T 9 2 F f w K / g b / g + u 7 0 r l g 1 r P O H i C 4 u Q V p E K s 3 < / l a t e x i t > Fock space spanned by â † i ( L = 0) 1 0 i S c P J J n 8 k r e r C f r x X q 3 P m a l B W v e c 0 j + w P r 8 A X 5 F m t k = < / l a t e x i t >

Target correlated state
< l a t e x i t s h a 1 _ b a s e 6 4 = "

SDs in physical space
FIG. 1. Depiction of the geometrical interpretation of the hidden fermion formalism.The Fock space spanned by the visible fermionic modes â † iσ is represented by the green horizontal line.The augmented Fock space is represented by the light orange plane (plane of the paper).The orange diagonal line represents the subspace in the augmented Fock space that is isomorphic to the physical Hilbert space after applying the constraint function (black arrows).The collection of Slater determinants (SDs) in the augmented space is represented by the blue shape, and the intersection with the subspace of just visible DOFs is marked in yellow.This intersection corresponds to the physical Hartree-Fock states.The constraint function changes the collection of states that represent the physical Hilbert space bringing the target correlated state close to a Slater determinant in the enlarged space.
In this work we are concerned with the subspace of definite particle numbers N ↑ and N ↓ of the individual spin species, in which case the two species are distinguishable from each other.However, it is convenient to impose full anti-symmetry between the spin species to enhance the expressivity of the family of trial wave functions, like in the so called unrestricted Hartree-Fock (HF).The projection to definite N ↑ and N ↓ subspace is imposed in the sampling of the wave function amplitudes.
+ F q 0 5 J 5 s 5 h j 9 w P n 8 A E C K N q g = = < / l a t e x i t > x 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " + F q 0 5 J 5 s 5 h j 9 w P n 8 A E C K N q g = = < / l a t e x i t > x 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " c s < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 a 3 j t c U q n y U p 6 j 2 x N < l a t e x i t s h a 1 _ b a s e 6 4 = " a k X w 6 w 6 z X f f 4 e 0 e e h t Z t P Y 6 p q + s = " d T j k I i e + n 0 h w w e a q U P g 0 j o 4 g p O 1 e 8 T K Q q l H I e + 7 g y R G s r f 3 k T 8 y + s k K n B 7 K e V x o g j H s 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " W 0 e H j g i 8 A W 6 M 0 J q I H 9 7 E / E v r 5 u q 0 P E y y p N U E Y 5 n H 4 U p g y q G k 0 B g j w q C F R t p g r C g e l e I B 0 g g r H R s B R 3 C 1 6 X w f 9 K 2 z e q x W b + q l x p n 8 z j y Y B 8 c g D K o g h < l a t e x i t s h a 1 _ b a s e 6 4 = " R t u c l k l 7 e l p e i y l w u G < l a t e x i t s h a 1 _ b a s e 6 4 = " U X < l a t e x i t s h a 1 _ b a s e 6 4 = " p P S v w 1 Y F C 5 J j z O y n 8 B 4 4 q / B 9 I G w = " > A < l a t e x i t s h a 1 _ b a s e 6 4 = " W g 7 L y J 2 A 4 U d q N t M p L N 8 3 A I e w B N 4 N j L j 0 X g x X p e t a 8 Z q 5 g T 8 g P H 2 C X X d k u s = < / l a t e x i t > N (f 1 (n)) < l a t e x i t s h a 1 _ b a s e 6 4 = " O t N g t u 1 c q + 5 S q v M 3 < l a t e x i t s h a 1 _ b a s e 6 4 = " + m o S D k S + / / W 3 S 2 q q P Y J l w b D + b g U = " > A A A B 7 X i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 B I v g q S R S 1 G P R i 8 c K 9 g P a U D a b T b t 2 s x t 2 J 0 I p / Q 9 e P C j i 1 f / j z X / j t s 1 B W x 8 M P N 6 b Y W Z e m A p u 0 8 5 h j + w P n 8 A b j t j z 0 = < / l a t e x i t > . . .< l a t e x i t s h a 1 _ b a s e 6 4 = " + m o S D k S + / / W 3 S 2 q q P Y J l w b D + b g U = " > A A A B 7 X i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 B I v g q S R S 1 G P R i 8 c K 9 g P a U D a b T b t 2 s x t 2 J 0 I p / Q 9 e P C j i 1 f / j z X / j t s 1 B W x 8 M P N 6 b Y W Z e m A p u 0 2. Hidden-fermion determinant state amplitudes with a neural-network parametrized constraint function.The top part of the determinant is constructed by slicing N rows from the top M rows of the Φ matrix, according to visible particle configuration x.Each row of the bottom sub-matrix φ h (f (x)), χ h (f (x)) (hidden sub-matrix) is parametrized by the outputs of a separate neural network (indicated by different colors), whose input is the flattened visible-lattice occupancy n.

II. HIDDEN FERMION FORMALISM AND WAVE FUNCTION ANSATZ
A. States in the augmented Hilbert space: constraint function Recall that the multi-particle physical Hilbert (Fock) space is spanned by M := 2V creation operators ĉ † iσ applied in all possible ways to the Fock vacuum |0 .The strategy of this paper is to define an augmented Fock space, constructed by M tot > M fermionic modes.
We partition the mode operators of the augmented Fock space into two species of auxiliary fermionic degrees of freedom â † µ and d † ν , referred to as visible and hidden modes respectively.We note that, although most hidden particle approaches enlarge the Hilbert space with bosonic degrees of freedom, fermionic hidden sectors have been considered in recent works [29][30][31][32][33][34] (see also 23,24,35 ).We require 1 ≤ µ ≤ M , and 1 ≤ ν ≤ M with M a free hyper-parameter.Of course M tot = M + M .The occupancy of the visible modes â † µ is identified one-to-one with the occupancy of the physical modes ĉ † iσ , establishing a direct correspondence between the index µ of the visible modes and the position-spin multiindex of the physical modes (i, σ).
Thus, the basis for the augmented Fock space is spanned by the set of states where n and ñ label the occupancy of the visible and hidden modes respectively.Note that this basis does not have a definite number of hidden fermions, even if the visible occupations are constrained to have definite particle number.Since the augmented Fock space defines a superset of the physical many-body fermionic states, a collection of "repre-sentative" states is chosen to span the basis of the physical Hilbert space within the augmented space; similarly to the constraint applied in the hidden rotor, spin or boson formalism 19,[25][26][27] .This choice produces a basis of the correct dimension, eliminating so called unphysical states.This constraint is applied by the following procedure: for each visible fermion occupancy n a particular hidden fermion occupancy ñ is chosen.The arbitrary choice of the population of the hidden modes can be summarized by a constraint function In the physical subspace the probability amplitude of the spinful fermion occupancy n is given by the overlap between the augmented basis states (Eq. 3) and a given trial state vector |Ψ of the augmented Fock space, where the hidden occupancy ñ is controlled by the visible occupancy n via the constraint: In this work we consider the search of the optimal constraint function, contrary to previous hidden-particle formulations where the constraint is a fixed physically-motivated rule.The resulting wave function ansatz is thus parametrized by both the choice of the state in the augmented space |Ψ and the constraint function F (n). Since there exist doubly exponentially many constraint functions, an extremely flexible family of correlated trial wave functions is obtained.
It should be noted that the hidden single-particle orbitals d † µ define an abstract space of hidden particle configurations.While one may be concerned by the nature of this abstract space and the form of the orthogonal set of single particle orbitals that define its basis, in practice, we work in the basis of particle configurations in the orbitals d † µ .Consequently, the relevant quantity is the combination of single-particle orbitals in the abstract space.
Fig. 1 illustrates geometrically the general concept of the hidden fermion formalism.The constraint function can be interpreted as a non-trivial rotation of the collection of states that constitute the basis of the physical Fock space embedded in the augmented space (light green horizontal line is rotated to the orange segment in Fig. 1).The goal of this transformation is to bring the target correlated state close to the parametrized family of states in the enlarged space.In Fig. 1 the chosen family of parametrized states is the family of SDs.We also show that in the particular case of M = 0 (light green subspace) the physical Fock space is directly spanned by the visible modes, and that standard SD states can be recovered in that limit.

Generalities
In order to demonstrate the versatility of the hidden fermion approach, we consider in this work the special case where |Ψ is the uncorrelated Slater determinant state |Ψ SD , which is characterized by a total number N tot ≥ N of orbital functions φ n : {1, . . .M tot } → C where 1 ≤ n ≤ N tot and Ñ = N tot − N is the number of added hidden fermions.In particular, |Ψ SD is obtained from the Fock vacuum as follows, where each φ † α is a linear combination of the original creation operators, whose coefficients are determined by the corresponding orbital.In terms of the row-vectors where Φ is the M tot × N tot matrix whose columns correspond to the orbital functions.It will be convenient to write the matrix of orbitals in the following block form, Where φ v is the M ×N matrix representing the amplitudes of the visible orbitals evaluated in the visible modes, χ v is the M × Ñ matrix representing the amplitude of the hidden orbitals evaluated in the visible modes, φ h is the M ×N matrix representing the amplitude of the visible orbitals evaluated in the hidden modes and χ h is the M × Ñ matrix representing the amplitude of the hidden orbitals evaluated in the hidden modes.
Since the SD state is an eigenstate of the total number operator, as are both the visible and hidden sectors, and anticipating that particle configurations are sampled in the VMC framework, we can represent the constraint as a mapping between the visible-particle configuration x = (x 1 , ..., x N ) and the hidden-particle configuration x = (x 1 , ..., xN ) In order to respect the Fermi statistics, it is sufficient to choose the function f to be of bosonic nature; that is, invariant under permutations of the visible configuration.The amplitudes of the wave function ansatz in the configuration basis are thus given by where φ v (x), χ v (x) and φ h (f (x)), χ h (f (x)) denote the N × (N + Ñ ) and Ñ × (N + Ñ ) sub-matrices obtained from φ v , χ v and φ h , χ h respectively by slicing the row entries corresponding to x and f (x).For convenience we denote the φ h (f (x)), χ h (f (x)) matrix as the hidden sub-matrix.

Universality and connection to other wave function ansätze
This ansatz is universal in the lattice.The proof relies on the ability of the determinant in Eq. 9 to represent a universal lookup table of amplitudes that are matched with the amplitudes of an arbitrary target state.In the particular case of φ h = 0 and χ v = 0 the flexibility of the ansatz is relied upon χ h , as φ v leads to amplitudes that correspond to an uncorrelated state in the physical space.It is possible to construct the lookup table for Ñ ≥ 1, requiring M to grow combinatorially fast with the number of physical fermionic modes.See the Supplementary Information (SI) 36 for a detailed discussion.It follows that in the general case where φ h = 0 and χ v = 0, the determinant in the enlarged Fock space does not inherit the nodes of the φ v orbitals.
Our construction bears some similarities with backflow [37][38][39][40] , in which orbitals are taken to be functions of the coordinates of all particles.In contrast to regular backflow, only the restriction of orbitals to hidden states have multiparticle position dependence (see (9)).
Jastrow-like wave function ansätze of the form where J(n) is an arbitrary function of the visible lattice occupations, also appear naturally in this formalism.This connection clearly materializes by considering Ñ = 1 and χ v = φ h = 0.In this case the amplitudes of the wave function ansatz are the product of det[φ v (x)] and a symmetric function of the visible particle configuration χ h (f (x)).Note that this class of wave functions includes the physically motivated Gutzwiller and Jastrow factors, as well as generalized neural-network Jastrow factors 9,11 applied to Slater determinants.The constraint function reproducing the Gutzwiller state can be found in the SI 36 .Configuration-interaction (CI) wave functions are also explicitly connected to the hidden fermion determinant state.Using the Laplace expansion of the determinant in Eq. 9 along its last Ñ rows yields a linear combination of Nparticle Slater determinants.If φ V and χ V are chosen to be the N lowest HF orbitals and the first Ñ virtual orbitals respectively, then a CI wave function is obtained, containing all possible (single to Ñ -tuple) excitations to the first Ñ virtual orbitals.See SI 36 for the detailed derivation.

Parametrized constraint function, practical implementation
In contrast to physically motivated constraint functions, a more general approach involves considering f (x) to belong to a parametrized family, whose parameters are optimized, together with the orbitals, in the energy minimization.The variational Monte Carlo seeks optimal parameters θ for a variational family of wave functions ψ θ (x), which are assumed to be differentiable with respect to θ.Although the requirement of differentiability appears to be in tension with the combinatorial nature of f (x), this obstacle is easily overcome by parametrizing instead the composition of functions φ(f (x)) and χ(f (x)), which appears in the hidden sub-matrix of the enlarged determinant (see ( 9)).This parametrization is connected to the notion of a continuous set of orthogonal hidden modes d † µ , which accounts to M → ∞.Remarkably, this automatically satisfies the condition that M must grow combinatorially with M for the determinant in Eq. 9 to be a universal lookup table of amplitudes.The hidden fermion configurations f (x) are thus represented by some internal state of the parametric function.However, in practice we are never interested in such internal state.
Since the hidden sub-matrix is a matrix-valued function which by construction is a permutation-invariant function of the visible configuration x, we choose to represent it by neural networks, taking as an input the visible occupation numbers n, without loss of generality.Neural networks are the perfect candidate to reduce the intractable complexity of choosing the optimal constraint, as they define an extremely flexible family of functions.Furthermore, sufficiently large neural networks can represent arbitrary contrainst functions, since they satisfy a universal approximation theorem 41 .The set of variational parameters θ of our ansatz consists on the the matrices φ v and χ v together with the weights and biases parametrizing the corresponding neural network.Fig. 2 details the precise parametrization used in this work.In practice, each row of the hidden sub-matrix is parametrized by its own neural network, as shown in Fig. 2 by different colored neural-network blocks.We consider multilayer perceptrons with hyperbolic-tangent activations.The hyperparameters of the ansatz include the neural network architecture as well as the number of added hidden fermions Ñ .
The cost of evaluating the enlarged determinant and its derivatives with respect to the variational parameters scales with the number of visible and hidden fermions as O (N + Ñ ) 3 , coming from the used LU factorization.Typically we choose Ñ ∼ O N , and therefore the asymptotic cost of the evaluation of the hidden fermion determinant state is O N 3 .The computation of the wave function amplitudes and gradients is the only step in the VMC algorithm where the required resources are larger, by a constant factor, than the computation of a N -fermion determinant.

C. Methods
Both the amplitudes of the matrices φ v and χ v , together with the weights and biases of the neural networks parametrizing the rows of the hidden sub-matrix, are jointly optimized using the stochastic reconfiguration method 43 , an extension of the classical natural gradient optimization method 44 to variational quantum states.Given that we are interested in the approximation of the ground-state wave function, we rely on the variational principle and use the expectation value of the Hamiltonian with respect to the variational state as the objective function to be optimized.For every Hamiltonian parameter choice a new trial state is optimized from scratch.
General expectation values and gradients of the objective function are computed using Markov Chain Monte Carlo sampling according to the probability distribution defined by the square of the wave function amplitudes |ψ θ (x)| 2 , working in the basis of particle configurations.We use the Python library NetKet 45 for the implementation (see the SI 36 for details), where gradients of the wave function amplitudes with respect to the variational parameters are computed by the so called automatic differentiation implemented in the Python library Jax 46 .

III. NUMERICAL EXPERIMENTS
In this section we benchmark the hidden fermion determinant wave function ansatz with a fully parametrized constraint function.We first study the square lattice at average site occupation n = 1/2 and n = 5/8.We use the 4 × 4 square lattice as a test-bed to study the accuracy (compared to exact diagonalization (ED)) of the proposed ansatz.The accuracy is quantified by the difference between the ED ground-state energy and the variational energy, relative to the ED ground-state energy.We analyze the effect of the neural network complexity and compare against relevant results in the literature.Lastly we focus on rectangular geometries of size 4 × L, where we consider 1/8 hole-doping (n = 7/8).Periodic boundary conditions are set in the short side of the rectangle in all cases.We study the case of both open and periodic boundary conditions on the long side.In the former, we compare our energies with Density Matrix Renormalization Group (DMRG) results and study the competing stripe orders of the system.In the latter and in the smallest system size (4 × 4) we analyse the relative error in the ground-state energy, obtained from ED.In the larger sizes (L = 8 and L = 16) we compare the ground state energy with the results obtained using a Slater-Jastrow ansatz and the neural network backflow wave function from Ref. 10 .In all cases we focus on the zero magnetization and fixed visible and hidden particle subspaces.

A. Benchmarks in the square lattice
We begin by considering the particular case of Ñ = N , which provides a good trade-off between computational complexity and accuracy; and a single-hidden-layer neural network parametrizing each row of the hidden sub-matrix.This architecture is a good starting point to study the effect of the neural network expressive power in the accuracy of the r w X q 2 X j 5 b U 9 Z o Z p n 8 g P X 6 A Y u 9 l c 8 = < / l a t e x i t >

Slater-RBM:
< l a t e x i t s h a 1 _ b a s e 6 4 = " 9 l h u k 8 N 3 8 0 A W w v 3 v X Z W c L / r 5 6 g g

Slater-RBM:
< l a t e x i t s h a 1 _ b a s e 6 4 = " 9 4 w t T 8 F S 8 e F P H q d 3 j z b 5 x J c t D E B w W P 9 6 q o q u d G U m i 0 r G 9 j Y X F p e W U 1 t 5 Z f 3 9 j c 2 j Z 3 d u s 6 j B W H G g 9 l q J o u 0 y

Slater-RBM
< l a t e x i t s h a 1 _ b a s e 6 4 = " 1 B 0 e s N g o F a x N s 9 p I K 4 T T e a 6 f S p p g M s E j 2 t M 0 w i G V 3 m x 5 7 h x d a G W I g l j o i h R a q t 8 n Z j i U c h r 6 u j P E a i x / e w v x L 6 + X q q D i z < l a t e x i t s h a 1 _ b a s e 6 4 = " S S m 9 m g t h m n + B G z r / T 8 a 9 W P T j / b f 7 n Y P X m z i 2 2 V P 2 j H W Z Y A N 2 w N 6 w I z Z i k n 1 m X 9 k 3 d h l 9 i S 6 i q + j 7 d e t W t J l 5 w v 5 A 9 O M n N f a l 1 w = = < / l a t e x i t > (Variance extrapolated Rel.error: 5.52(6) • 10 5 ) < l a t e x i t s h a 1 _ b a s e 6 4 = " L z 2 y p 6 t G q 3 y e 3 r h r / T C g Q p x 5 Y f g D a n U V u 2 T R n p + b 0 a e e U m u p K M o u g l W H q 0 + f v J 0 7 V n t + Y u X r 9 b r r 9 8 M n C m s w L 4 w y t j T B B w q q b F P k h S e 5 h Y h S x S e J B e f K  FIG. 3. Exact diagonalization benchmarks of the ground-state energy in the 4 × 4 lattice with periodic boundary conditions.(a) Relative error in the ground-state energy as a function of the inverse of the width density α of the single-hidden-layer neural networks parametrizing the rows of the hidden sub-matrix.Average physical site occupation is n = 1/2 and Ñ = 8.Different values of U are considered, as indicated by each color.The error for a Slater-RBM ansatz (see main text) with hidden neuron density α = 32, at the same values of U is included for comparison.Indicated is also the relative error from the variance-extrapolated energy for each value of U (see SI 36 for details).(b) Relative error in the ground-state energy as a function of the coupling constant U , at n = 5/8 average site occupancy (first closed shell) and Ñ = 10.The rows of the hidden sub-matrix are given by single-hidden-layer neural networks with α = 64.The error from Slater-Jastrow and Slater-RBM ansätze are included for comparison.The green diamond is the relative error found with the stateof-the-art, tensor-network-based ansatz from Ref. 42 .Shown is also the relative error according to the projection of the converged hidden fermion determinant state to the subspace of invariant wave functions under the action of π/2 rotations (C4) and the group of all possible translations T with K = 0 momentum, separately and together.
ansatz.In this case, the expressive power is only determined by the number of hidden units.More hidden units improve the flexibility of the neural network.Furthermore, this single hidden layer architecture is the minimal architecture that satisfies the universal approximation theorem 41 .
Panel (a) of Fig. 3 shows the relative error in the ground state energy as a function of the ratio between the number of hidden units and input features (α), at n = 1/2 average site occupation.Different values of U are shown, including challenging cases (U = 7.75 and U = 10) where the ground state is strongly correlated.There is a systematic trend to decrease the error in the energy as α is increased, providing a clear and controllable pathway to obtaining more expressive wave function ansätze.Moreover, for the largest values of α, and contrary to what is observed on typical wave function ansätze, the error does not significantly increase with U as the correlations in the ground state increase.Remarkably, the relative error at U = 10, typically the most challenging case, is orders of magnitude lower than the error of the Slater-RBM wave function ansatz.The Slater-RBM ansatz is a particular case of the wave function in Eq. 10, where J(n) is a restricted Boltzmann Machine of complex weights.
The direct extrapolation of the relative error to the α → ∞ limit is challenging as the asymptotic scaling of the accuracy with the neural network complexity is not a well understood matter in the field.However, from the different energy and variance estimates obtained for each α we perform an energy-variance extrapolation procedure 47 to obtain better estimates to the ground-state energy.See SI 36 for details.The relative error corresponding to the variance extrapolated energies is shown in Fig. 3 (a), where the relative error is in this case defined as the difference between the variance extrapolated and the ground-state energies, relative to the ground-state energy.
The improvement of the accuracy with the increase of α is accompanied by a gentle increase in the computational complexity of the determinant amplitudes.The scaling with α is linear, as the evaluation of the elements hidden submatrix requires O Ñ (M • αM + αM • (N + Ñ )) operations, coming from the two afine transformations of the fully connected neural networks with a single hidden layer.For reference, the scaling of the evaluation of the neuralnetwork backflow from 10 is O(N 3 ), from the evaluation of the determinant of multi-particle orbitals, while the evaluation of the matrix elements that enter the determinant requires to store M distinct fully connected neural networks and O N (M • αM + αM • (N ) operations.This makes the asymptotic scaling of the hidden fermion determinant state with N , M and α identical to the scaling of the neuralnetwork backflow In principle, deeper architectures provide a greater expressive power than their shallower counterparts 48 , at the expense of a higher computational cost.We observe that, while deeper architectures provide marginal gains in the energy error, increasing the number of hidden fermions yields a greater impact on the accuracy of the ansatz (see the SI 36 for a detailed study of the effect of increasing Ñ and the depth of the neural networks in the accuracy of the ansatz).
Benchmarks on physically motivated constraint functions were also performed (see the SI 36 for details).Our experiments reveal that parametrizing f is advantageous compared to the physically inspired rigid rules, which show a marginal improvement in accuracy as compared to the Slater-Jastrow state.
At n = 5/8 average site occupation we can compare the relative error in the energy against the state-of-the-art ansatz from Ref. 42 .n = 5/8, which corresponds to N = 10, is the first closed shell for the model under consideration in the non-interacting limit.Fig. 3 (b) shows the relative error in the ground-state energy as a function of U .The error does not increase monotonically with the value of U , as standard wave function ansätze do.This is shown by the reference errors displayed by the Slater-Jastrow and Slater-RBM ansätze.Remarkably, the single Slater determinant ansatz with the parametrized constraint function outperforms (by a factor of two in the relative error) the relative error reported in Ref. 42 that uses the state-of-the-art ansatz consisting on a pairing reference state multiplied by Jastrow, Gutzwiller and doublon-holon correlation factors as well as fat tree tensor network of bond dimension 16, all projected into the zero momentum singlet subspace, with enforced C 4 rotational symmetry.Remarkably, while the result from Ref. 42 relies on the projection of the trial state onto given symmetry sectors, our ansatz achieves better accuracy with no symmetry projections.Symmetry projections are an independent avenue to improve the accuracy.So far, we have only considered the increase of the neural-network complexity to obtain better trial states.Not surprisingly, our ansatz is further improved when projected to relevant symmetry subspaces after its convergence, as shown in panel (b) of Fig. 3.
See the SI 36 for more benchmarks in square geometries at half filling, where we compare the variational energies from the hidden fermion determinant state at increasingly larger system sizes with Auxiliary Field Quantum Monte Carlo (AFQMC) calculations 51 .Our energies are in better agreement to AFQMC than those obtained with the neuralnetwork Jastrow wave function from Ref. 9 .

B. Increasing system size and stripe order at 1/8 hole-doped rectangular geometries
To conclude this work, we investigate the validity of the proposed wave function ansatz on increasingly larger system sizes in rectangular geometries.In particular, we choose rectangular lattices of dimensions L×4, with L = {4, 8, 16}.We focus on the 1/8 hole doped and zero total magnetization subspace, where, in the strong coupling regime, the ground state is expected to show hole stripes every eight lattice sites across the long side of the rectangle (λ = 8).The high hole density regions coincide with domain walls in the antiferromagnetic order 49,52 .In this section we study the particular case of U = 8.For this particular choice of coupling constant and filling, previous works have found different competing orders close in energy to the λ = 8 stripe order 49 .
To guide the wave function ansatz towards the λ = 8 stripe order, we add a soft mean-field constraint to the φ v matrix of orbitals.A M × N matrix is added to the variational φ v .This matrix has zeros everywhere except for a single entry in every column.These entries are filled up with a constant factor S that multiplies max(|φ v |), following the charge and spin order described above.This forces each of the N visible orbitals to peak in a certain position of the physical lattice.The value of S is a hyper-parameter.The wave function is optimized with this constraint until its convergence, then the guiding matrix is merged into φ v as part of the variational parameters and the energy optimization is continued.
A good trade-off between accuracy and computational resource use is achieved by the addition of Ñ = 16 hidden fermions for all system sizes.We also consider a two layer fully connected neural network of hidden-unit density α = {60, 14, 6} for the L = {4, 8, 16} sizes respectively, to parametrize the rows of the hidden sub-matrix.
We investigate the case with periodic boundary conditions on the short side of the rectangle and open boundary conditions on the long side (PBC-OBC).The left panel of Fig. 4 (a) shows the energy per site as a function of L and the comparison with DMRG variational energies used in Ref. 49 .The DMRG algorithm finds two meta-stable solutions, one with half-filled stripes and one with filled stripes.In the L = 8 case, we can stabilize both meta-stable arrangements by tuning the value of S. S = 0 or small values of S lead to a halffilled stripe configuration of higher energy.In this system size the charge distribution shows high hole density every four sites, coinciding with domain walls in the antiferromagnetic order.Larger values of S yield a filled stripe configuration, showing only one stripe of high hole density in the system, that coincides with a domain wall in the antiferromagnetic order.The charge and spin configurations for the two competing orders are shown on the right panel of Fig. 4 (a).They are in good agreement with the DMRG hole distributions.The variational energy from the hidden fermion determinant state is in good agreement with the DMRG energies.At L = 16, the best variational energy found by tuning S lies between the DMRG energies that correspond to the half and filled stripes.Our method, not being specifically tailored to quasi one-dimensional problems, does not outperform DMRG in this particular lattice geometry.
A more interesting case is the addition of periodic boundary conditions along the long side of the rectangle, a situation that is not amenable to DMRG calculations due to its computational cost.The left panel of Fig. 4 (b) shows the energy per site as a function of L and the comparison with the ED energy in the 4 × 4 system.The hidden fermion determinant with a parametrized constraint function achieves a significantly lower energy than both the standard Slater-Jastrow ansatz and the state-of-the-art neural network backflow ansatz from Ref. 10 .In the 4 × 4 lattice, the relative error in the energy is reduced by over one order of magnitude compared to the neural network backflow wave function.In the 4 × 8 and 4 × 16 lattices the energy per site is noticeably lower than it from the neural network backflow and Slater-Jastrow ansätze.These results demonstrate the escalability of the proposed formalism, which outperforms existing state- g h G y y y E + E h Z E 1 S s P q c g k M x V A T y i T X t 1 q s T y V l O g t V 0 C E 4 s y / P k / p J y T k t l W / K x c r F N I 4 8 2 S c H 5 I g 4 5 I x U y B W p k h p h J C P P 5 J W 8 G U / G i / F u f E x a c 8 Z 0 Z p f 8 g f H 5 A 5 8 3 l x Q = < / l a t e x i t > Slater-Jastrow < l a t e x i t s h a 1 _ b a s e 6 4 = " P J O c 6 U x T h Z 5 V 7 7 m 0 L j 3 Z I y w r a q c = " > A A A B / X i c b V D L S g M x F M 3 U V 6 2 v 8 b F z E y y C q z I j R V 0 W 3 b g q F e w D 2 q F k 0 k w b m p k M y R 2 1 D s V f c e N C E b f + h z v / x r S d h b Y e u H A 4 5 9 7 k 3 u P H g m t w n G 8 r t 7 S 8 s r q W X y 9 s b G 5 t 7 9 i 7 e w 0 t E 0 V Z n U o h V c s n m g k e s T p w E K w V K 0 Z C X 7 C m P 7 y a + M 0 7 p j S X 0 S 2 M Y u a F p B / x g F M C R u r a B x 1 g D 6 D C t F r

DMRG (filled stripe)
< l a t e x i t s h a 1 _ b a s e 6 4 = " 3 o 5 u L h T m a 4 0 1 l f 9 R 3 k H e n Y 9 5 a 8 H J Z w 7 h D 5 z P H / L 0 j P 0 = < / l a t e x i t > 15 < l a t e x i t s h a 1 _ b a s e 6 4 = " O I 1 6 a i C F A D + g J P V t 3 1 q P 1 Y r 0 u W 9 e s 1 c w J + g H r 7 R P s 8 4 0 L < / l a t e x i t > 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 a d C 7 L 9 2 v j 4 a a i C F A D + g J P V t 3 1 q P 1 Y r 0 u W 9 e s 1 c w J + g H r 7 R P s 8 4 0 L < / l a t e x i t > 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 a d C 7 L 9 2 v j 4 a

2
< l a t e x i t s h a 1 _ b a s e 6 4 = " S v E w 9 j + G 6 / n N S Z e b h l v Q w R v l S i I = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m k q M e i F 4 8 t 2 F p o Q 9 l s J + 3 a z S b s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 a a T 2 J A t s Z U T P S q 9 5 c / M / r p i a 8 8 a d c J q l B y Z a L w l Q Q E 5 P 5 1 2 T A F T I j J p Z Q p r i 9 l b A R V Z Q Z m 0 3 B h u C t v r x O W l d l 7 7 p c a V R K t d s s j j y c w T l c g g d V q M E 9 1 K E J D B C e 4 R X e n E f n x X l 3 P p a t O S e b O Y U / c D 5 / A I a j j M Q = < / l a t e x i t > 7 < l a t e x i t s h a 1 _ b a s e 6 4 = " S v E w 9 j + G 6 / n N S Z e b h l v Q w R v l S i I = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m k q M e i F 4 8 t 2 F p o Q 9 l s J +  of-the-art wave function ansätze used in the field.For these cases we set S = 3.
In addition, we analyze the hole density and staggered spin density distributions in the largest system size (4 × 16) in the right panel of Fig. 4 (b).The hole density distribution shows repeating maxima separated by 8 lattice sites.Coinciding with the maxima in the hole density, the antiferromagnetic order displays a domain wall.The amplitude of the staggered magnetization is modulated along the long side of the rectangles.These features are consistent with observations from previous studies 49,52 coming from different many-body numerical methods, further validating the accuracy of the hidden fermion determinant state to find good approximations to the highly correlated ground states of complex Hamiltonians.

IV. CONCLUSIONS
In this paper we have shown that the variational treatment of interacting electrons in an augmented Fock space can be highly beneficial to improve the generality of the wave function ansatz, especially in the strong correlation limit.We found that key elements for the success of this approach are to optimize the constraint function relating the enlarged and physical Hilbert spaces, as well as treating this constraint exactly.A simple Slater determinant state in the augmented Fock space is found to provide an extremely expressive wave function ansatz, which we proved to be universal.The optimization of the constraint function and of the hidden sector orbital amplitudes was performed using a neural network representation.We presented numerical experiments which show that the proposed wave function is competitive with the state-of-the-art variational accuracy in the ground state of the Hubbard model in the square lattice.Furthermore, in contrast to standard variational approaches 9,42 , the accuracy of this ansatz does not rely on imposing symmetries, potentially allowing for a great level of accuracy on systems with a small number of symmetries.This also opens the possibility of applying our approach to systems without an underlying lattice, with potential applications to quantum chemistry, nuclear physics and materials science.
In particular, we envision the accurate calculation of the ground-state properties of molecular Hamiltonians, which in the Molecular Orbital basis lacks both an underlying lattice and exploitable symmetries.The connection between a compact representation of a CI wave function with all possible single-, double-, up to Ñ -tuple excitations provides an accurate post Hatree-Fock starting point for the trial state.The formalism introduced in this manuscript is also well-suited for models with quenched disorder, such as models of interacting fermions on fully connected lattices with random exchange interactions 53 .While exact solutions of such models are available in the large-M limit when the spin symmetry is extended to SU (M ), the physical SU (2) case requires computational approaches.Similarly to Molecular Hamiltonians, this class of Hamiltonians also lacks any translational symmetries that can be exploited to improve the accuracy of traditional wave-function ansätze.
consequently, det χ h (f (x)) can be interpreted as the determinant of multiple-visible-particle orbitals.Therefore, in order to have ψ(x) = ψ tar (x), the determinant det χ h (f (x)) must satisfy: which can only be accomplished iff det[φ v (x)] = 0 for all x.This requirement is satisfied by choosing the rows of the (M × N ) φ v matrix to be non-colinear N -dimenstional vectors.Note that both ψ tar (x) and det φ v (x) are anti-symmetric functions, and in consequence, det χ h (f (x)) must be a symmetric function of x, which is achieved by taking f (x) to be symmetric.Now we provide an explicit construction of f (x) and χ h such that requirement A2 is satisfied for an arbitrary state ψ tar (x).This construction is based on identifying det χ h (f (x)) with a lookup table with α max = M N entries, that provides the corresponding amplitudes for the different α max occupation configurations n of the visible particles on the visible modes.Therefore, in the worst case scenario we require M = α max + ( Ñ − 1) hidden modes.A possible choice for the M × Ñ matrix of amplitudes χ h that achieves the desired lookup table is given by with I ( Ñ −1) the ( Ñ − 1) × ( Ñ − 1) identity matrix.We identify χ 1 , ..., χ αmax with the α max different amplitudes of the ratio ψ tar (x)/ det φ v (x) .The position of hidden fermion i is determined by i th component of the constraint function, i.e. f i (x), which must satisfy 0 < f i (x) ≤ M and f i (x) ∈ Z.We choose where 0 < α(x) ≤ α max and α ∈ Z provides a distinct label to each of the M N different occupation configurations of the visible particles on the visible modes.
With the construction we obtain det χ h (f (x)) = χ m = ψ tar (x)/ det φ v (x) with 1 ≤ m ≤ α max , thus concluding the proof that the hidden fermion determinant state can match the amplitudes of an arbitrary state in the fixed particle subspace of Fock space.
It is noteworthy that the explicit construction of a combinatorially-large number of hidden modes can be circumvented by directly parametrizing the function composition χ h (f (x)) by a universal function approximator.
Following the Quantum Chemistry notation |Φ µ α with α ≤ N and µ > N , is the excitation of the fermion in occupied orbital α to the virtual orbital µ: Analogously |Φ µ1...µn α1...αn labels the promotion of the n electrons in occupied orbitals α 1 through α n to virtual orbitals µ 1 through µ n .
We choose to identify the visible orbitals φ v with the N occupied HF orbitals φ HF α (1 ≤ α ≤ N ) and the Ñ hidden orbitals χ h with the Ñ first virtual orbitalsφ HF µ (N < µ ≤ N + Ñ ).
1. Ñ = 1 and single excitations According to the above identification of orbitals, the wave function amplitudes become: The cofactor of the above determinant C ij (x) is defined as where A ij (x) is the minor of the matrix of orbitals, obtained by eliminating row i and column j from the matrix that enters the original determinant.Using the Laplace expansion of the determinant along its last row: We realize that C Ntot,α (x) can be identified with the amplitudes of the N -particle Slater determinant obtained from promoting particle in occupied orbital α to the first virtual orbital.In the particular case where α = N tot , C Ntot,Ntot (x) are the amplitudes of the HF ground state.The prefactors to the cofactors can be set to constants . For this choice of constraint function and orbital amplitudes, the projected augmented determinant to the physical space corresponds to the CI wavefunction with single-particle excitations: (B5)

Ñ = 2 and double excitations
In this case, the amplitudes of the hidden fermion determinant state are given by the determinant: We now proceed to expand the above determinant on its cofactors by its last row, obtaining the linear combinations of (N +1)particle determinants: C i,j is the cofactor of the (N + 2)-particle determinant C ij = det A ij where A ij is the minor obtained by removing row i and column j from the matrix entering the determinant in Eq.B6.Now C Ntot,α are (N + 1)-particle determinants.The cofactors C Ntot,α can also be expanded on their cofactors by their last row, obtaining: where A (kl)(ij) is the minor obtained by removing row k and column l from the previously defined minor A ij .Now C (Ntot−1,β)(Ntot,α) are N -particle Slater determinants obtained by promoting the particles in the occupied orbitals α and β to the two lowest virtual orbitals φ HF N +1 and φ HF N +2 .It must be noted that the limiting cases C (Ntot−1,Ntot−1)(Ntot,Ntot) , C (Ntot−1,Ntot−1)(Ntot,α) and C (Ntot−1,β)(Ntot,Ntot) are the HF ground-state determinant, the excitation of the particle in occupied orbital α to the second virtual orbital and the excitation of the particle in occupied orbital β to the first virtual orbital receptively.As before, setting the prefactors to the cofactors to be the constants: the augmented determinant, projected to the physical space corresponds to the the CI wavefunction with single and doubleparticle excitations: )   3. General Ñ and Ñ -tuple excitations The generalization of the two cases above to arbitrary Ñ implies that the hidden fermion determinant state of Ñ hidden fermions can be chosen to be a compact representation of a CI wavefunction with all possible single through Ñ -tuple particle excitations to the Ñ lowest virtual orbitals: (B11) Note that in the above expression the coefficients c µ1...µn α1...αn are given, up to a sign, by the product 1≤i,j≤n φ αi (f j (x)).

Concluding remarks
It must be noted that for arbitrary Ñ not all of the coefficients of the above expansion can be independent from each other.This claim is supported by a simple counting argument.There are Ñ • ( Ñ + N ) of the φ HF α (f (x)) coefficients from the bottom hidden sub-matrix.However, in general , the number of free parameters in the above CI expansion is: which except for Ñ = 1 is larger than Ñ • ( Ñ + N ).While this imposes a constraint on the type of CI wavefunctions that can be obtained with this construction, it must be noted that other successful trial states like unitary coupled clusters also suffer from similar constraints.
While the above construction assumed the staring orbitals to be the occupied and virtual orbitals from the HF solution of the given Hamiltonian, a more general case concerns the choice of arbitrary linear combinations of those orbitals.In practice we allow such arbitrary linear combinations, and therefore, the trial state we use in this work contains, as a particular case, the CI wavefunction.
Finally it is important to remark that in this derivation we have set φ HF α (f (x)) to be constants.However, in the general case they will be bosonic functions of the visible particle positions, which leads to Jastrow-like factors in the CI expansion, thus leading to a more general wave function.

Appendix C: Physically motivated constraint functions
In this appendix we first show an explicit construction of the Gutzwiller wave function ansatz from the hidden fermion determinant state.This construction serves as a motivation to showcase the fact that a simple constraint function leads to a correlated wave function.Inspired by this, we present other physically motivated (non-optimizable) constraint functions, and benchmark them in the 4 × 4 Hubbard model at quarter occupation.We also study the possibility of adding neural-networkbased correlation factors to the hidden fermion Slater determinant state.

Explicit construction of the Gutzwiller wave function in the Slater determinant hidden fermion state
In order to gain intuition about the nature of the hidden fermion determinant state, consider as a trivial warmup, the molecular orbital wave function, which has the following Fock-space representation corresponding to N = 2 fermions in M = 4 modes, by introducing one hidden fermion ( Ñ = 1) with two hidden states ( M = 2) and using the following disjointly supported orbitals, In particular, with the following choice of constraint function we find the probability amplitudes which is of the required form.By similar reasoning, one can construct the Gutzwiller projection of an arbitrary Slater determinant state with 2M modes described by orbitals φ v , simply by choosing disjointly supported hidden and visible orbitals φ h = 0, χ v = 0, χ h = 1, e −g , . . ., e −M g T together with the constraint function k t p N N T x 7 j Y 6 f 0 c K y N K w V 4 q v 6 e y I i 0 d i Q j 1 y k J D O y 8 N x H / 8 8 I U 4 u t O x l W S A l N 0 t i h O B Q a N J / / j H j e M g h g 5 Q q j h 7 l Z M B 8 Q Q C i 6 l o g s h m H 9 5 k T T O K s F l 5 f z + o l y 9 y e M o o E N 0 h E 5 Q g K 5 Q F d 2 h G q o j i j R 6 R q / o z Q P v x X v 3 P m a t S 1 4 + c 4 D + w P v 8 A Z D V k X I = < / l a t e x i t > a) < l a t e x i t s h a 1 _ b a s e 6 4 = " / l 0 D + U s < l a t e x i t s h a 1 _ b a s e 6 4 = " p 6 P

Benchmarks with physically motivated constraint functions
Motivated by the ability to reproduce correlated wave function ansätze as described above, we propose and test a collection of constraint functions built from the physical intuition about the ground state physics of the Hamiltonian in Eq. (1) in the main text.
a. Gutzwiller-inspired constraint function Take Ñ = 1 and M = N/2, with the constraint function f (x) = M i=1 n i↑ n i↓ .The position of the hidden fermion in the hidden modes is determined by the number of doubly occupied sites in the x configuration.The wave function ansatz is the Slater determinant in the enlarged Hilbert space of Eq. (5) in the main text.The variational parameters are the coefficients (orbitals) of the non-orthogonal change of single-particle basis in Eq. (6) in the main text.For a particular choice of orbitals this ansatz can reproduce the Gutzwiller wave function as described above.

b. Unrestricted hidden spin ansatz
Consider a pair of hidden fermionic modes d † iβ with β = 0, 1, associated to each site of the visible Hubbard lattice i, and a single hidden fermion populating each pair of hidden modes.Thus Ñ = M/2 and M = M .The position of the hidden fermion depends on the visible site occupancy via the local constraint function f i (x) = n i↑ n i↓ = β.The wave function ansatz is the Slater determinant in the enlarged Hilbert space of Eq. (5) in the main text.The variational parameters are the coefficients (orbitals) of the non-orthogonal change of single-particle basis in Eq. (6) in the main text.This ansatz is closely related to the hidden spin formalism 27 .The difference is to consider the hidden fermions to be indistinguishable particles amongst each other and the visible ones, as opposed to the hidden spin formalism where the added spins are distinguishable particles.This unrestricted hidden spin ansatz provides more flexibility, the same way the unrestricted Hartree-Fock does compared to factorised HF.It is easy to see that the unrestricted hidden spin ansatz can reproduce the Gutzwiller ansatz for a particular choice of orbitals.Furthermore, it was proven that the hidden spin formalism can capture the Mott transition at the mean field level.Therefore, this ansatz can capture correlations beyond the Gutzwiller single-site correlations.Table III provides the variational (for the hidden fermion determinant state with a fully parametrized hidden sub-matrix) and exact diagonalization (ED) energies in the square lattice Hubbard model of size 4×4 and average site occupation n = 5/8 (first closed shell of the model) and periodic boundary conditions, for various values of U .The hidden sub-matrix is parametrized by single-hidden-layer neural networks of width density α.Projections to different symmetry subspaces are applied to the converged trial state: rotation of π/2 (C 4 ), translations of momentum K = 0 and the intersection of both subspaces.See the main text for more details.The data corresponds to the relative errors displayed in Fig. 3    Here we present the converged ground-state energies of the hidden fermion determinant state with a fully parametrized hidden sub-matrix in the Hubbard model in rectangular geometries of dimension 4 × L at 1/8 hole doping (n = 7/8) and U = 8.Periodic boundary conditions are considered along the short side of the rectangle, while both open and periodic boundary conditions are taken along the long side of the rectangle.In Table IV both configurations are labelled as PBC-PBC and PBC-OBC respectively.Table IV shows the variational energy per site obtained with the hidden fermion determinant state with a fully parametrized hidden sub-matrix.The value of α depends on the system size (see Section 3. B in the main text for details).The energies are those displayed in Fig. 4  Table V shows the energy per site on the square lattice Hubbard model of size L × L with periodic boundary conditions along one of the sides of the square and anti-periodic boundary conditions along the other side.The trial state is the hidden fermion determinant state with a fully parametrized hidden sub-matrix.Details on the number of hidden fermions and hidden unit densities can be found on Section 6 of the Supplementary Information.
< l a t e x i t s h a 1 _ b a s e 6 4 = " d z M 8 q s M Q a Z l N H U 2 0B I u L G L G V h d o = " > A A A C A X i c b V C 7 S g N B F J 2 N r x h f q z a C z W A Q r M K u B L W M 2 l h G M A 9 I Q p i d 3 C R D Z h / M 3 B X D E h t / x c Z C E V v / w s 6 / c T b Z Q h M P D B z O u Y + 5 x 4 u k 0 O g 4 3 1 Z u a X l l d S 2 / X t j Y 3 N r e s X f 3 6 j q M F Y c a D 2 W o m h 7 T I E U A N R Q o o R k p Y L 4 n o e G N r l O / c Q 9 K i z C 4 w 3 E E H Z 8 N A t E X n K G R u v Z B G + E B l Z 9 c x g M f A o Q e 1 R H j M O n a R a f k T E E X i Z u R I s l Q 7 d p f 7 V 7 I 4 3 Q I l 0 z r l u t E 2 E m Y Q s E l T A r t W I O Z P G I D a B k a M B 9 0 J 5 l e M K H H R u n R f q j M C 5 B O 1 d 8 d C f O 1 H v u e q f Q Z D v W 8 l 4 r / e a 0 Y + x e d R A R R j B D w 2 a J + L C m G N I 2 D 9 o Q C j n J s C O N K m L 9 S P m S K c T S h F U w I 7 v z J i 6 R + W n L P S u X b c r F y l c W R J 4 f k i J w Q l 5 y T C r k h V V I j n Dy S Z / J K 3 q w n 6 8 V 6 t z 5 m p T k r 6 9 k n f 2 B 9 / g A 1 R J d i < / l a t e x i t >Augmented space< l a t e x i t s h a 1 _ b a s e 6 4 = "1 8 M D o v n s U V U D J O N O z 6 L l A 9 J P L z k = " > A A A B + H i c b V B N S 8 N A E N 3 4 W e t H o x 6 9 L B a h X k o i R T 1 J 0 Y s X o Y L 9 g D a U z W b T L t 1 s w u 5 E q K G / x Is H R b z 6 U 7 z 5 b 9 y 2 O W j r g 4 H H e z P M z P M T w T U 4 z r e 1 s r q 2 v r F Z 2 C p u 7 + z u l e z 9 g 5 a O U 0V Z k 8 Y i V h 2 f a C a 4 Z E 3 g I F g n U Y x E v m B t f 3 Q z 9 d u P T G k e y w c Y J 8 y L y E D y k F M C R u r b p U o P u A h Y d j f B V 9 g 5 7 d t l p + r M g J e J m 5 M y y t H o 2 1 + 9 I K Z p x C R Q Q b T u u k 4 C X k Y U c C r Y p N h L N U s I H Z E B 6 x o q S c S 0 l 8 0 O n + A T o w Q 4 j J U p C X i m / p 7 I S K T 1 O P J N Z 0 R g q B e 9 q f i f 1 0 0 h v P Q y L p M U m K T z R W E q M M R 4 m g I O u G I U x N g Q Q h U 3 t 2 I 6 J I p Q M F k V T Q j u 4 s v L pH V W d c + r t f t a u X 6 d x 1 F A R + g Y V Z C L L l A d 3 a I G a i K K U v S M X t G b 9 W S 9 W O / W x 7 x 1 x c p n D t E f W J 8 / 8 3 K R / A = = < / l a t e x i t >

< l a t
e x i t s h a 1 _ b a s e 6 4 = " c W d a F u a S C N 7 U c 6 o A w h 0 X F c r V B D Q = " > A A A B + H i c b V B N S 8 N A E N 3 4 W e t H o x 6 9 L B a h X k o i R b 0 I R S 9 e h A r 2 A 9 p Q N p t N u 3 S z C b s T o Y b + E i 8 e F P H q T / H m v 3 H b 5 q C t D w Y e 7 8 0 w M 8 9 P B N f g O N / W y u r a + s Z m Y a u 4 v b O 7V 7 L 3 D 1 o 6 T h V l T R q L W H V 8 o p n g k j W B g 2 C d R D E S + Y K 1 / d H N 1 G 8 / M q V 5 L B 9 g n D A v I g P J Q 0 4 J G K l v l y o 9 4 C J g 2 d 0 E X 2 H n t G + X n a o z A 1 4 m b k 7 K K E e j b 3 / 1 g p i m E Z N A B d G 6 6 z o J e B l R w K l g k 2 I v 1 S w h d E Q G r G u o J B H T X j Y 7 f I J P j B L g M F a m J O C Z + n s i I 5 H W 4 8 g 3 n R G B o V 7 0 p u J / X j e F 8 N L L u E x S Y J L O F 4 W p w B D j a Q o 4 4 I p R E G N D C F X c 3 I r p k C h C w W R V N C G 4 i y 8 v k 9 Z Z 1 T 2 v 1 u 5 r 5 f p 1 H k c B H a F j V E E u u k B 1 d I s a q I k o S t E z e k V v 1 p P 1 Y r 1 b H / P W F S u f O U R / Y H 3 + A P H r k f s = < / l a t e x i t > ( M = 0) < l a t e x i t s h a 1 _ b a s e 6 4 = " v 8 4 x 0 T h y v w 8 f z b m 1 l E A q H 0 E e v D A = " > A A A C C X i c b V C 7 S g N B F J 2 N r x h f U U u b w S B Y h V 0 R t Q z a W E b I C 5 I l z M 7 e T Y b M P p i 5 K 4 Y l r Y 2 / Y m O h i K 1 / Y O f f O J u k 0 M Q D w x z O u Zd 7 7 / E S K T T a 9 r d V W F l d W 9 8 o b p a 2 t n d 2 9 8 r 7 B y 0 d p 4 p D k 8 c y V h 2 P a Z A i g i Y K l N B J F L D Q k 9 D 2 R j e 5 3 7 4 H p U U c N X C c g B u y Q S Q C w R k a q V + m P Y Q H V G H W Y G o A S H m s F E i G 4 F O N 5 p v 0 y x W 7 a k 9 B l 4 k z J x U y R 7 1 f / u r 5 M U 9 D i J B L p n X X s R N 0 M 6 Z Q c A m T U i / V k D A + Y g P o G h q x E L S b T S + Z 0 B O j + D S I l X r J 9 j 5 N 0 6 S L T T x w M D h n P u Y e / x I C o O u + + 0 s L C 4 t r 6 x m 1 r L r G 5 t b 2 7 m d 3 a o J Y 8 2 h w k M Z 6 r r P D E i h o I I C J d Q j D S z w J d T 8 w e X Y r 9 2 DN i J U d z i M o B W w n h J d w R l a q Z 0 7 a C I 8 o A 6 S 2 y t D h a I s 7 g W g E D r U R I z D q J 3 L u w V 3 A j p P v J T k S Y p y O / f V 7 I Q 8 H g / h k h n T 8 N w I W w n T K L i E U b Y Z G 7 C T B 6 w H D U s V C 8 C 0 k s k h I 3 p k l Q 7 t h t o + h X S i / u 5 I W G D M M P B t Z c Cw b 2 a 9 s f i f 1 4 i x e 9 5 K h I p i B M W n i 7 q x p B j S c S q 0 I z R w l E N L G N f C / p X y P t O M o 8 0 u a 0 P w Z k + e J 9 W T g n d a K N 4 U 8 6 W L N I 4 M 2 S e H 5 J h 4 5 I y U y D U p k w r h 5 J E 8 k 1 f y 5 j w 5 L 8 6 7 8 z E t X X D S n j 3 y B 8 7 n D + e E m e k = < / l a t e x i t > SDs in augmented space < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 h l d 4 c 8 b O i / P u f C x a C 0 4 + c w x / 4 H z + A C / g j S U = < / l a t e x i t > n < l a t e x i t s h a 1 _ b a s e 6 4 = " V g i L g 5 e B u Q x S 8 u u C n 5 H s 8 P l t R y k = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e x K U I 9 B L 5 4 k o n l A s o T Z y S Q Z M j u 7 z P S K Y c k n e P G g i F e / y J t / 4 y T Z g y r a B D + P o U / k 9 a j m k f m 7 W r W q l + N o 8 j D / b B A S g D G 5 y A O r g A D d A E G N y B B / A E n o 1 7 4 9 F 4 M V 5 n r T l j P r M L f s B 4 + w S 9 r J e b < / l a t e x i t > 1 (f Ñ (n)) t e x i t s h a 1 _ b a s e 6 4 = " 0 V f e F w U F E + W P U 0 m v e 8 3 5 K H r 5 b N 1 w h r P r J E f s F 4 / A I X 0 l c s = < / l a t e x i t >Slater-RBM:< l a t e x i t s h a 1 _ b a s e 6 4 = " c C P C 6 H y W 4 B g sL P 6 l X K 2 e O U l O 8 m 0 = " > A A A C A H i c d V D L S g N B E J y N 7 / i K e v D g Z T A I X l x 2 4 2 r U U 4 g X L 0 J 8 J A p J C L O T T j I 4 + 2 C m V w x L L v 6 K F w + K e P U z v P k 3 7 i Y R V L S g o a j q p r v L D a X Q a F k f R m Z i c m p 6 Z n Y u O 7 + w u L S cW 1 m t 6 S B S H K o 8 k I G 6 d p k G K X y o o k A J 1 6 E C 5 r k S r t y b 4 9 S / u g W l R e B f Y j + E p s e 6 v u g I z j C R W r n 1 B s I d a D d p L R C 0 N 6 k C o d 6 o U q r Q D p S P 0 9 k T B f 6 4 H v p p 0 + w 5 6 e 9 j L x P 6 8 Z o 3 f m J C K I Y o S A j x d 5 s a Q Y 0 i w P 2 h E K O M p B S h h X I r 2 V 8 h 5 T j K d p 6 C w E e / r l W V I 7 L t o n x d J 1 q V A + n 8 S R I 3 t k n x w S m 5 y S M r k k F V I l n D y S Z / J K 3 o w n 4 8 V 4 N z 7 G r X P G Z G a H / I H x + Q M 8 Y p d S < / l a t e x i t > Slater-Jastrow < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 o 5 u L h T m a 4 0 1 l f 9 R 3 kC 9 v + e o 4 V o = " > A A A C C X i c b V A 9 S w N B E N 2 L 3 / E r a m m z G A S r c C e i l q I W K S O a K C Q h 7 O 3 N 6 Z L d v W N 3 T g x H W h v / i o 2 F I r b + A z v / j Z v k C r 8 e D D z e m 2 F m X p h K Y d H 3 P 7 3 S 1 P T M 7 N z 8 Q n l x a X l l t b K 2 3 r J J Z j g 0 e S I T c x U y C 1 J o a K J A C V e p A a Z C C Z d h / 2 T k X 9 6 C s S L R F z h I o a v Y t R a x 4 A y d 1 K v Q D s I d G p X X R R S B p j E Y 5 Q x 6 f k o t M o R h r 1 L 1 a / 4 Y 9 C 8 J C l I l B R q 9 y k c n S n i m Q C O X z N p 2 4 K f Y z Z l Bw S U M y 5 3 M Q s p 4 n 1 1 D 2 1 H N F N h u P v 5 k S L e d E t E 4 M a 4 0 0 r H 6 f S J n y t q B C l 2 n Y n h j f 3 s j 8 T + v n W F 8 2 M 2 F T j M E z S e L 4 k x S T O g o F h o J A x z l w B H G j X C 3 U n 7 D D O P o w i u 7 E I L f L / 8 l r d 1 a s F / b O 9 u r H h 0 X c c y T T b J F d k h A D s g R q Z M G a R J O 7 s k j e S Y v 3 o P 3 5 L 1 6 b 5 P W k l f M b J A f 8 N 6 / A I / o m k A = < / l a t e x i t > Hidden fermion SD state < l a t e x i t s h a 1 _ b a s e 6 4 = " T o H d A a K 7 e a s u 3 d Y U S l V F 4 p U 3 p 9 M = " > A A A B 6 n i c d V D L S s N A F L 3 x W e u r 6 t L N Y B F c h a Q k a d 0 V u 3 F Z 0 T 6 g D W U y n b R D J w 9 m J k I J / Q Q 3 L h R x 6 x e 5 8 2 + c P g Q V P X D h c M 6 9 3 H t P k H I m l W V 9 G G v r G 5 t b 2 4 W d 4 u 7 e / s F h 6 e i 4 L Z N M E N o i C U 9 E N 8 C S c h b T l m K K 0 2 4 q K I 4 C T j v B p D H 3 O / d U S J b E d 2 q a U j / C o 5 i F j G C l p d v G w B m U y p Z p u 5 d e p Y Y s 0 3 M 9 x 7 U 1 q b D F e p w D U 1 o A Y E R P M A T P B v c e D R e j N d l 6 5 q x m j m B H z D e P g F A b 4 3 O < / l a t e x i t > C 4 < l a t e x i t s h a 1 _ b a s e 6 4 = " n T 0 n y a I 3 J 3 V e n H f n Y 9 F a c P K Z Y / Q H z u c P + D i R p g = = < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " c M h o z M 3 a D k k 9 r / 0 e N Q 4 5 V n 3 R V L o = " > A A A C K n i c d V B N b x M x E P W W r x I + G u D I x S J C S g + s 1 m k 3 i T g V e u F Y E E k r J S G a 9 U 5 a q 1 5 7 Z c 8 i o l V + D x f + C p c e q C q u / B C 8 a Z A A w Z N G e n 5 v R p 5 5 W a m V p y S 5 i r Z u 3 L x 1 + 8 7 2 3 d a 9 + w 8 e 7 r Q f P R 5 7 W z m J I 2 m 1 d S c Z e N T K 4 I g U a T w p H U K R a T z O z g 8 b / / g j O q + s e U / L E m c F n B q 1 U B I o S P P 2 q + 6 U 8 B O 5 o h 6 D U 2 A k 8 u Y J p d V A m P N 3 q 0 9 A W O l 3 5 e I c L A j y 6 d Z 8 C L 8 u 5 f 8 n g 1 Y Y d 8 L 2 U b u x t 7 + M Y 4 2 9 Z Z u s y W L W Z X v s g B 2 y P h P s O / v B r t n P 4 D K 4 C m 6 C 2 4 f W l W A 5 s 8 H + Q H B 3 D x Y R p c A = < / l a t e x i t > (Variance extrapolated Rel.error: 9.86(3) • 10 5 ) < l a t e x i t s h a 1 _ b a s e 6 4 = " e d I 2 7 1 B f g p V 6 b 2 6 2 4 m L / + + 2 F I o g = " > A A A C K n i c d V B N b x M x E P W W r x K + A h y 5 W E R I 6 Y G V N 9 m k t K d C L x w L I m m l J E S z 3 k l r 1 W u v 7 N m q 0 S q / p x f + C p c e q C q u / B C c N E i A 4 E k j P b 8 3 I 8 + 8 r N T K k x D X 0 c a t 2 3 f u 3 t u 8 3 3 j w 8 N H j J 8 2 n z 4 b e r 0 k 3 p 1 6 o K / C k r O Z 9 a F M s R X 6 u 8 T N R T e z 4 s s d B Z A J / 5 v b y n + y x t V N H s z q Z U p K 0 I j b z 6 a V Z q T 5 c v c e K 4 c S t L z Q E A 6 F X b l 8 g Q c S A r p N k I I v y 7 l / y f D T p z 0 4 / R D 2 t p 7 t 4 5 j k 7 1 g L 1 m b J W y b 7 b H 3 7 I A N m G Q X 7 C v 7 x q 6 i L 9 F l d B 1 9 v 2 n d i N Y z z 9 k f i H 7 8 B B S N p c A = < / l a t e x i t > (Variance extrapolated Rel.error: 1.69(2) • 10 4 ) < l a t e x i t s h a 1 _ b a s e 6 4 = " j O 3 b a a p Q 8 z Q D 6 W p N 8 P T Z j p Q H p 4 U = " > A A A B 9 H i c b V D L S g N B E J y N r x h f U Y 9 e B o M Q L 2 F X g n o M P s B j B P O A Z A m z k 0 k y Z G Z 2 n e k N h i X f 4 c W D I l 7 9 G G / + j Z N k D 5 p Y 0 F B U d d P d F U S C G 3 D d b y e z s r q 2 R 0 9 w S d W 6 e J e q G 0 p w D P 1 9 0 R C p D F j G d h O S W B g F r 2 p + J / X i q F 3 6 S d c R T E w R e e L e r H A E O J p A r j L N a M g x p Y Q q r m 9 F d M B 0 Y S C z S l n Q / A W X 1 4 m 9 b O S d 1 4 q 3 5 c L l a s 0 j i w 6 Q s e o i D x 0 g S r o D l V R D V H 0 i J 7 R K 3 p z R s 6 L 8 + 5 8 z F s z T j p z i P 7 A + f w B X J O R 1 w = = < / l a t e x i t > (ED) < l a t e x i t s h a 1 _ b a s e 6 4 = " t o Z 9 i Z I Q 0 r C P e o g u Q J p p T u g o L X S p 7 R P T 3 a 8 3 R p v K q S e e 5 v s 5 / u N TN b N u U P Q v 8 Q b k y w Z o 9 j I f N a a i s c B h M g l M 6 b q u R H W E 6 Z R c A n 9 d C 0 2 E D H e Z W 2 o W h q y A E w 9 G T 7 V p 9 t W a d K W 0 r Z C p E P 1 5 0 T C A m N 6 g W 8 7 A 4 Y d M + k N x P + 8 a o y t 4 3 o i w i h G C P l o U S u W F B U d J E S b Q g N H 2 b O E c S 3 s r Z R 3 m G Yc b Y 5 p G 4 I 3 + f J f c n O Q 8 w 5 z + c t 8 t n A 2 j m O e b J I t s k M 8 c k Q K 5 I I U S Y l w 8 k C e y A t 5 d R 6 d Z + f N e R + 1 p p z x z A b 5 B e f j G 7 0 o n V M = < / l a t e x i t > Relative error: 9(1) • 10 4 < l a t e x i t s h a 1 _ b a s e 6 4 = " h w b + 5 R b y H n s W t j E S 6 1 g d i 7 W o 5 W Q = " > A A A C G 3 i c b Z A 9 S w N B E I b 3 4 l e M X 1 F L m 8 U g x M J w F + I H V k E b S x V j h C S G v c 3 E L N m 7 P X b n x H D k f 9 j 4 V 2 w s F L E S L P w 3 7 m k K T X x h 4 e W Z G W b n 9 S M p D L r u p 5 O Z m p 6 Z n c v O 5 x Y W l 5 Z X 8 q t r l 0 b 0 p 5 j 2 n G 0 c a Z s y F 4 4 y d P m s t y y d s r V c 4 q h e r R K I 4 s 2 S C b p E g 8 s k + q 5 I S c k h r h 5 J 4 8 k m f y 4 j w 4 T 8 6 r 8 / b T m n F G M + v k j 5 y P L x 7 T n u g = < / l a t e x i t > Relative error: 1.58(0.001)• 10 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " e b d U G Z Q q 4 D n x x J 1 Z R 8 H i 9 K y L R k w = " > A A A C G 3 i c b Z A 9 S w N B E I b 3 / D Z + R S 1 t F o M Q C 8 N d N C p W o o 1 l F B M D S Q x 7 m 0 m y u H d 7 7 M 4 F w 5 H / Y e N f s b F Q x E q w 8 N + 4 F 1 P 4 9 c L C y z M z z M 7 r R 1 I Y d N 0 P Z 2 J y a n p m d m 4 + s 7 C 4 r G t b B / p b z H N O N o 4 8 z Y E L z f J / 8 1 1 W L B 2 y / s n e / l j k / G c c y R D b J J 8 s Q j B + S Y n J E y q R B O 7 s g D e S L P z r 3 z 6 L w 4 r 1 + t E 8 5 4 Z p 3 8 k P P + C R 1 F n u c = < / l a t e x i t > Relative error: 5.35(0.001)• 10 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " J v I H r N 2 S s + w f 1 2 h m + N d R G v i N P T c = " > A A A C A H i c b V D J S g N B E O 2 J W 4 z b q A c P X g a D 4 M U w I 0 E 9 B r 2 I p 4 h m g W Q I P Z 2 a p E n P Q n e N G o a 5 + C t e P C j i 1 c / w 5 t / Y W Q 6 R M 3 p F b 9 a T 9 W K 9 W x + z 1 p y V z e y j P 7 A + f w D R i 5 V 5 < / l a t e x i t >NN backflow< l a t e x i t s h a 1 _ b a s e 6 4 = " 3 o 5 u L h T m a 4 01 l f 9 R 3 k C 9 v + e o 4 V o = " > A A A C C X i c b V A 9 S w N B E N 2 L 3 / E r a m m z G A S r c C e i l q I W K S O a K C Q h 7 O 3 N 6 Z L d v W N 3 T g x H W h v / i o 2 F I r b + A z v / j Z v k C r 8 e D D z e m 2 F m X p h K Y d H 3 P 7 3 S 1 P T M 7 N z 8 Q n l x a X l l t b K 2 3 r J J Z j g 0 e S I T c x U y C 1 J o a K J A C V e p A a Z C C Z d h / 2 T k X 9 6 C s S L R F z h I o a v Y t R a x 4 A y d 1 K v Q D s I d G p X X R R S B p j E Y 5 Q x 6 f k o t M o R h r 1 L 1 a / 4 Y 9 C 8 J C l I l B R q 9 y k c n S n i m Q C O X z N p 2 4 K f Y z Z l Bw S U M y 5 3 M Q s p 4 n 1 1 D 2 1 H N F N h u P v 5 k S L e d E t E 4 M a 4 0 0 r H 6 f S J n y t q B C l 2 n Y n h j f 3 s j 8 T + v n W F 8 2 M 2 F T j M E z S e L 4 k x S T O g o F h o J A x z l w B H G j X C 3 U n 7 D D O P o w i u 7 E I L f L / 8 l r d 1 a s F / b O 9 u r H h 0 X c c y T T b J F d k h A D s g R q Z M G a R J O 7 s k j e S Y v 3 o P 3 5 L 1 6 b 5 P W k l f M b J A f 8 N 6 / A I / o m k A = < / l a t e x i t > Hidden fermion SD state < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 I 9 c M H q W b w I 2 W B 4 k 6 K o Z D 2 U r C b Y = " > A A A B + X i c b V B N S w M x E M 3 6 W e v X q k c v w S J 4 s e x K U Y + l v X i z g v 2 A d i n Z N G 1 D k + y S z B b L 0 n / i x Y M i X v 0 n 3 v w 3 p u 0 e t P X B w O O 9 G W b m h b H g B j z v 2 1 l b 3 9 j c 2 s 7 t 5 H f 3 9 g 8 O 3 a P j h o k S T V m d R i L S r Z A Y J r h i d e A g W C v W j M h Q s G Y 4 q s 7 8 5 p h p w y P 1 C J O Y B Z I M F O 9 z S s B K X d f t A H s C L d N a p X p 5 X 6 l O u 2 7 B K 3 p z 4 F X i Z 6 S A M t S 6 7 l e n F 9 F u c P 0 T O T I g = = < / l a t e x i t > PBC-OBC < l a t e x i t s h a 1 _ b a s e 6 4 = " h S 3 P w W D / Q w v e W 6 p p K / qI z O n o b z o = " > A A A B + X i c b V B N S w M x E M 3 W r 1 q / V j 1 6 C R b B i 2 V X i n o s 7 c V j B f s B 7 V K y a d q G J t k l m S 2 W p f / E i w d F v P p P v P l v T N s 9 a O u D g c d 7 M 8 z M C 2 P B D X j e t 5 P b 2 N z a 3 s n v F v b 2 D w 6 P 3 O O T p o k S T V m D R i L S 7 Z A Y J r h i D e A g W D v W j M h Q s F Y 4 r s 3 9 1 o R p w y P 1 C N O Y B Z I M F R 9 w S s B K P d f t A n s C L d N 6 t X Z l a 9 Z z i 1 7 J W w C v E z 8 j R Z S h 3 n O / u v 2 I J p I p o I I Y 0 / G 9 G I K U a O B U s F m h m x g W E z o m Q 9 a x V B H J T J A u L p / h C 6 v 0 8 S D S t h T g h f p 7 I i X S m K k M b a c k M D K r 3 l z 8 z + s k M L g L U q 7 i B J i i y 0 W D R G C I 8 D w G 3 O e a U R B T S w j V 3 N 6 K 6 Y h o Q s G G V b A h + K s vr 5 P m d c m / K Z U f y s V K N Y s j j 8 7 Q O b p E P r p F F X S P 6 q i B K J q g Z / S K 3 p z U e X H e n Y 9 l a 8 7 J Z k 7 R H z i f P 9 K 6 k y M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " T Q z d y V U 3 d 7 5 0 K / e g 9 I i D O 6 w H 0 H T Z 5 1 A e I I z N F I r s 9 d A e E D l J x f X t 5 c 0 5 w k p o U 0 1 K h H B 4 a C V y d p 5 e w Q 6 S 5 w J y Z I J S q 3 7 s k j e S Y v 3 o P 3 5 L 1 6 b 5 P W k l f M b J A f 8 N 6 / A I / o m k A = < / l a t e x i t > Hidden fermion SD state t e x i t s h a 1 _ b a s e 6 4 = " Y 3 < / l a t e x i t > 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " t Y l D W 7 p I Z F 2 I + R H D s A N N T e o D v J k = " > A A A B 6 H i c b V D L T g J B E O z F F + I L 9 e h l I j H x R H Y N E Y 9 E L x 4 h k U c C G z I 7 9 M L I 7 O x m Z t a E E L 7 A i w e N 8 e o n e f N v H G A P C l b S S a W q O 9 1 d Q S K 4 N q 7 7 7 e Q 2 N r e 2 d / K 7 h b 3 9 g 8 O j 4 v F J S 8 e p Y t h k s 3 a z S b s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 B + f C o r e N U M W y x W M S q E 1 C N g k t s G W 4 E d h K F N A o E P g T j 2 5 n / 8 I R K 8 1 j e m 0 m C f k S H k o e c U W O l p t s v V 9 y q O w d Z J V 5 O K p C j 0 S 9 / 9 Q Y x S y O U h g m q d FIG. 4. Energy per site and competing charge and spin orders in the 4 × L rectangular lattice at 1/8 hole doping (n = 0.875) and U = 8.(a) Periodic boundary conditions on the short side of the cylinder and open on the long side (PBC-OBC).The left panel compares the hidden fermion determinant state energies with DMRG energies.The width of the DMRG symbols shows the range of converged variational energies for different bond dimensions used in Ref. 49 .For L = 8, blue points labelled as [1] and [2] correspond to filled and half filled stripes.The right panel shows the hole and staggered spin distribution for both meta-stable configurations.The diameter of the grey circles is proportional to the hole density.(b) Periodic boundary conditions along both sides of the rectangles (PBC-PBC).The left panel compares the hidden fermion determinant state energies with the Slater-Jastrow and neural network backflow ansätze (from Ref. 10 ).The dashed horizontal line marks the ED (4 × 4 with PBCs from Ref. 50 ) energy.In the 4 × 4 lattice the relative error in the ground state energy is displayed for each ansatz.The right panel shows the hole and staggered spin distributions in the 4 × 16 lattice.
t e x i t s h a 1 _ b a s e 6 4 = " p q b 4 1 H J b C e b n j z G x 0 7 p 4 V g b V w r w V P 0 9 k R F p 7 U h G r l M S G N h 5 b y L + 5 4 U p x N e d j K s k B a b o b F G c C g w a T / 7 H P W 4 Y B T F y h F D D 3 a 2 Y D o g h F F x K R R d C M P / y I m m c V Y L L y v n 9 R b l 6 k 8 d R Q I f o C J 2 g A F 2 h K r p D N V R H F G n 0 j F 7 R m w f e i / f u f c x a l 7 x 8 5 g D 9 g f f 5 A 5 J b k X M = < / l a t e x i t > b) n m P a W Z j 0 C Z r Q / B m X 5 4 n l Z O C d 1 o o 3 h T z p Y t p H B l y Q A 7 J M f H I G S m R a 1 I m P u H k k T y T V / L m P D k v z r v z M W l d c K Y z e + Q P n M 8 f G s y X T w = = < / l a t e x i t > Unrestricted SD < l a t e x i t s h a 1 _ b a s e 6 4 = " s 6 l f y z / L F 2 c g 7 s z y y l U l e y n O y g y y Y y p e W 6 E j Y R p F F x C P 1 W P D U S M 3 7 A O 1 C x V L A D T S I b v 9 G n W K i 3 a D r U t h X S o / p 5 I W G B M L / B t Z 8 C w a y a 9 gf i f V 4 u x f d R I h I p i B M V H i 9 q x p B j S Q T a 0 J T R w l D 1 L G N f C 3 k p 5 l 2 l m M 9 E m Z U P w J l + e J u X 9 v H e Q L 1 w W M s X j c R x L Z I f s k h z x y C E p k n N y Q U q E k 0 f y T F 7 J m / P k v D j vz s e o d c Y Z z 2 y T P 3 A + f w C n L J x 7 < / l a t e x i t > Unrestricted SD (Gutzwiller) < l a t e x i t s h a 1 _ b a s e 6 4 = " K q S r 9 u b v m 8 L w o 8 8 z r 5 J 7 n V z / W f Q = " > A A A C C 3 i c b V A 9 T w J B E N 3 z E / H r 1 N J m I z H B h t w Z o p Z E L Y y V R h E S I G R v G W D D 3 t 5 l d 0 4 l F 3 o b / 4 q N h c b Y + g f s / D c u S K H g S y Z 5 + 9 5 M d u Y F s R Q G P e / L m Z m d m 1 9 Y z C x l l 1 d W 1 9 b d j c 0 b E y x L e y u l H e Z Z j Y Q b b I 2 B H / y 5 G l y s 1 / w D w r F y 2 K u d D y O I 0 O 2 y Q 7 J E 5 8 c k h I 5 I x e k T D h 5 I E / k h b w 6 j 8 6 z 8 + a 8 / 7 T O O O O Z L f I H z s c 3 E f K b D g = = < / l a t e x i t > Unrestricted SD (Jastrow) < l a t e x i t s h a 1 _ b a s e 6 4 = " e F w j + n 9 s Q Y T d t Z T Y 5 X c 5 0 3 w 8 b 6 I FIG. 5. Benchmarks of physically motivated constraint functions with ED energies in the 4×4 Hubbard model at n = 1/2 average physical site occupation.Results from standard wave function ansätze are shown as dashed lines for comparison purposes.(a) Relative error in the ground-state energy as a function of the coupling constant U .The different constraint function ansätze are a single Slater determinant in the augmented Fock space with no projections.(b) Same as panel (a) including a complex RBM projection factor both in the control unrestricted HF ansatz and a E-RBM factor in the the hidden fermion ansätze.

3. 4
× L Hubbard model with n = 7/8 and U = 8 in the main text.L Energy per site HFDS (PBC-PBC) Energy per site HFDS (PBC-OBC) Variational energy per site in the Hubbard model in rectangular geometries of dimensions 4 × L at 1/8 hole doping (n = 7/8) and U = 8, obtained with the hidden fermion determinant state with a fully parametrized hidden sub-matrix.The value of α depends on the system size (see Section 3. B in the main text for details).The energies are those displayed in Fig.4in the main text.4. L × L Hubbard model with n = 1 25693(1) -0.9120(1) -0.68135(1) -0.5401(1) 6 -1.2079(1) -0.8717(2) -0.6609(4) -0.5289(3) 8 -1.1900(2) -0.8621(4) -0.6574(2) -0.5244(6) TABLE V. Variational energy per site in the L × L Hubbard model at half filling with periodic boundary conditions along one of the sides of the square and anti-periodic boundary conditions along the other side.The trial state is the hidden fermion determinant state with a fully parametrized hidden sub-matrix (see Section 6 in the Supplementary Information for more details).The energies correspond to the relative errors shown in Fig4in the Supplementary Information.