{'fraction_non_alphanumeric': 0.037608103917582604, 'fraction_numerical': 0.012679598938772697, 'mean_word_length': 4.33471188309898, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The phase diagram of the Hubbard model by Variational Auxiliary Field quantum Monte Carlo

A systematically improvable wave function is proposed for the numerical solution of strongly correlated systems. With a stochastic optimization method, based on the auxiliary field quantum Monte Carlo technique, an effective temperature Teff is defined, probing the distance of the ground state properties of the model in the thermodynamic limit from the ones of the proposed correlated mean-field ansatz. In this way their uncertainties from the unbiased zero temperature limit may be estimated by simple and stable extrapolations well before the so called sign problem gets prohibitive. At finite Teff the convergence of the energy to the thermodynamic limit is indeed shown to be possible in the Hubbard model already for relatively small square lattices with linear dimension L ~10, thanks to appropriate averages over several twisted boundary conditions. Within the estimated energy accuracy of the proposed variational ansatz, two clear phases are identified, as the energy is lowered by spontaneously breaking some symmetries satisfied by the Hubbard Hamiltonian: a) a stripe phase where both spin and translation symmetries are broken, and b) a strong coupling d-wave superconducting phase when the particle number is not conserved and global U(1) symmetry is broken. On the other hand the symmetric phase is stable in a wide region at large doping and small coupling.


I. INTRODUCTION
The accurate numerical solution of the Schroëdinger equation, namely determining the ground state of a many-body Hamiltonian H, remains the most challenging unsolved problem since the Dirac's formulation in 1931. Historically, important progress, as well as too optimistic promises, occurs periodically at least twice a decade, starting from quantum computers, the Density Functional Theory Kohn-Sham formulation [1], the Density Matrix Renormalization Group (DMRG) [2,3] and its translation within the tensor network quantum information language [4][5][6], wrong claims [7][8][9] about the solution of the sign problem in Quantum Monte Carlo (QMC), systematically improvable wave functions (WFs) based on multi reference expansions [10][11][12] or on machine learning assumptions [13,14], just to mention a few of them.
Unfortunately, determining the exact ground state of a strongly correlated Hamiltonian H in the thermodynamic limit, remains an open issue apart for particular cases as in one dimension [15] or with very particular couplings [16]. Also for this reason, the Hubbard model has been historically used to benchmark new techniques, because its simplicity represents an ideal playground for advanced computational and also experimental methods such as optical lattices [17].
As far as the 2D Hubbard model is concerned only a few results have been established by numerical techniques, namely a clear antiferromagnetic phase at one electron per site filling [18][19][20] is present and a strong evidence that this phase remains even when this condition is not fulfilled [21]: the added holes (unoccupied sites) are expelled from the antiferromagnet, and essentially fill in equally spaced vertical lines of the lattice: the so called stripes.
Instead, the question of superconductivity in the Hubbard model remains highly debated and controversial since the dis- * sorella@sissa.it covery of high-temperature superconductivity [22][23][24] to the last few months [25]. Probably the most clear evidence of superconductivity was reported in the strong coupling limit of the Hubbard model, the so called t-J model. Several years ago, this model has been studied [26] with an almost standard variational Monte Carlo method (VMC) that simply relies on the variational principle: the best wave function ansatz, defined by a set of variational parameters, is the one that minimizes the expectation value of the Hamiltonian studied. At that time, soon after the discovery of high-temperature superconductivity, it was reported that superconductivity did not need electron-phonon interaction but the driving force was rather the superexchange spin interaction J. The approach and especially the claim has been highly debated and remains controversial until now, even though it is worth mentioning that several years later, by tensor network, P. Corboz et al. have reported, without emphasizing it, almost the same value [27] of the off-diagonal d − wave superconducting long range order.
Despite this success and its simplicity, the standard VMC based on a single reference mean-field ansatz, e.g. the Hartree-Fock (HF), corrected by a simple correlation term, e.g. the Gutzwiller factor, is certainly limited as compared with the most recent and advanced variational methods, that allow a systematically improvable ansatz.
In this work, a different strategy is proposed to overcome the limitations of correcting a simple mean-field ansatz for reaching accurate ground state properties. The approach will be dubbed in the following variational auxiliary field quantum Monte Carlo (VAFQMC), that takes advantage of the enormous progress done in the last decades by two old but well established techniques: the variational quantum Monte Carlo (VMC) and the auxiliary field quantum Monte Carlo (AFQMC) [28][29][30][31][32]. From the former technique its simplicity and clarity in interpreting the results in the thermodynamic limit as well as its ability to optimize a mean-field state ψ MF in presence of a simple correlated factor are taken. Conversely from AFQMC a systematically convergent and non perturbative expression of the correlation term is used.
Indeed AFQMC allows the application of exp(−Hτ/2) to any mean-field state, where τ is the total imaginary time used in this formulation. Thus, the projection exp(−Hτ/2) filters out exactly the ground state component of ψ MF for large imaginary time τ, yielding that the ansatz |ψ τ 0 = exp(−Hτ/2)|ψ MF , basically used in this work, represents a systematically improvable-with increasing τ-correlated mean-field ansatz, retaining all advantages of the two mentioned formalisms.
The simple working hypothesis of this study is the following. Suppose that, in the thermodynamic limit, the mean-field wave function has acquired the lowest possible energy upon breaking some symmetry of the Hamiltonian in presence of an enough accurate correlation term. In this case, the implicit VMC assumption is that the exact ground state should eventually show this phase. On the other hand, with the proposed method one can in principle tune the accuracy of the electron correlation at the desired level and verify explicitly the systematic evolution of the phase diagram, derived in this way, from the simplest mean-field HF theory (τ = 0) to the converged one. This limit can be approached very closely at half-filling. Unfortunately, for the finite doping case, the exact (τ = ∞) limit cannot be met in principle because the method presented here is vexed by the so called sign problem. Nevertheless extremely accurate variational energies can be obtained even with short τ projection.
In short the high energy degrees of freedom (e.g. the average number of doubly occupied sites) of the variational ansatz are settled almost immediately with short τ projection (e.g. they are filtered out with the projection with an exponetial decay proportional to a large gap U). As long as the low energy degrees of freedom and the high energy ones are almost decoupled, an hypothesis implicitly assumed in several theoretical and computational methods in physics and chemistry, the phase diagram derived with this ansatz should converge fast with τ, because only the high energy degrees of freedom have to be settled by the projection.

II. MODEL AND WAVE FUNCTION
In this work, the square lattice Hubbard model H = K µ +V is studied, where K µ = ∑ k,σ c † k,σ c k,σ (−2t(cos k x +cos k y )− µ) is the kinetic energy, defined here by the hopping t and the chemical potential µ and V = U ∑ i n ↑,i n ↓,i is the total number of doubly occupied sites operator scaled by the coupling U > 0 of the model, where standard second quantization notations are assumed. We introduce the following ansatz [33][34][35], also inspired by similar wave functions quite popular in quantum computation [36,37], with the purpose to reach accurate ground state properties in the shortest possible projection time τ: As shown later, the above variational wave function represents a straightforward improvement of the standard auxiliary field QMC one |ψ τ 0 , mentioned in the introduction.
The mean-field wave function |ψ MF can be a generic quasifree electron state, from a simple Slater determinant to BCS pairing functions, including singlet and/or triplet correlations. Here, |ψ MF is allowed breaking some symmetries of the Hamiltonian, and is defined as the ground state of a mean-field Hamiltonian H MF (α 0 ), i. e. H MF (α 0 )|ψ MF = E 0 (α 0 )|ψ MF with a set of variational parameters indicated by the vector α 0 . For the projection operator, we take the advantage of the variational formulation, so that any extension of the variational ansatz containing a given one-i.e. |ψ τ 0 -as a particular case should necessarily improve it, namely after energy optimization, acquires a lower variational energy. The bare kinetic energy operator in exp(−Hτ/2) is therefore generalized (i.e. K µ → H MF (α)), by allowing a generic operator quadratic in the fermion ones c, c † , including for instance also a d-wave BCS pairing field. H MF (α) is therefore parameterized by a set of variational parameters indicated by another independent vector α, such that for α = 0 H MF (α) = K µ . In other words the idea is that if a symmetry is broken in the thermodynamic limit also the projection operator, i.e. the correlation factor, and not only the mean-field state |ψ MF may break the symmetry and therefore H MF (α) is conveniently parametrized in a way similar to ψ MF .
This original formulation improves the quality of the ansatz and the smooth convergence to the thermodynamic limit, as compared with the simpler ansatz |ψ τ 0 and the optimal energy is therefore generally obtained after the simultaneous optimization of both α and α 0 .
We anticipate that, in this variational formulation of the auxiliary field quantum Monte Carlo (VAFQMC), τ plays the role of an effective inverse temperature that is kept fixed during the minimization of the energy expectation value corresponding to |ψ τ .
The optimization techniques known in standard variational Monte Carlo [38] and machine learning [39] will be generalized here to the auxiliary field QMC. Before that, it is worth to emphasize simple but important properties of this ansatz: 1. It is systematically improvable. In order to realize this property, it is enough to take α = 0, when the proposed ansatz coincides with the simpler one |ψ τ 0 = exp(−Hτ/2)|ψ MF . Thus, let τ → ∞ for a mean field |ψ MF of the chosen form, that is not orthogonal to the ground state. In this limit |ψ τ 0 is obviously converging to the exact ground state. Thus, after turning on optimization, by allowing both α = 0 and, independently, α 0 different from the initial guess, a lower energy is necessarily implied for each τ, yielding that the ansatz of Eq. 1 is systematically improvable, as well as and even better than |ψ τ 0 as far as its energy is concerned.
2. It is size extensive. As discussed in Ref.28 the WF is defined directly in terms of an exponential of an extensive operator, hence the statement. In practice this means that, at a given τ, approximately the same accuracy for intensive quantities is expected, e.g. the energy per site or bulk correlation functions.
3. For finite clusters the convergence is exponential in τ due to the finite size gap between the ground state manifold (which may be also degenerate) and the first excitation with non zero energy gap. However in the thermodynamic limit this gap is probably always vanishing in this model, even for the half filled insulator, where gapless spin-wave excitations are expected due to the occurrence of antiferromagnetic order for any U > 0 [19]. This situation, as it will be shown in the following, makes the extrapolation to the 1 τ = 0 unbiased limit much simpler than the corresponding finite size case. Indeed this can be obtained by simple and stable low order polynomial extrapolations in something that can be considered an effective temperature T eff = 1/τ.
The key idea of this work is to converge first the results to the thermodynamic limit with large enough cluster size simulations and appropriate boundary conditions. Given this, simple and very stable extrapolations in T eff → 0 are employed, thus achieving, with this simple minded strategy, very accurate results of the model or at least an estimate of the accuracy of the lowest T eff variational ansatz. As presented later, this is often possible because, after the optimization, the physical properties of the ansatz given in Eq.(1) are already in the very low temperature regime, where the mentioned extrapolations are indeed stable and reliable within a given broken symmetry phase.
It is well established [40][41][42] that, when a standard type of order that breaks a continuous symmetry sets in, the corresponding gapless low energy excitations (i.e. typically bosons with a density of states ρ(ε) ∝ ε D−1 ) induce T 3 ( T 2 ) energy corrections (i.e. ∆E ∝ T 0 ερ(ε) T D+1 ) in the limit T → 0 for 2D (quasi 1D system like a finite cylinder with infinite length), and indeed it turns out that the fit is generally very appropriate for the systems studied here, because the available T eff appear always quite close to the correct asymptotic behavior. A formal derivation of the finite effective temperature corrections is given in App.A. Given the above arguments, particularly useful boundary conditions are adopted such that the convergence to the thermodynamic limit at fixed τ is as fast as possible. To this purpose twisted averaged boundary conditions TABC [43,44] are used in both the Hubbard Hamiltonian and the mean-field ones on rectangular L x × L y clusters (L x = L y = L is adopted for square lattices). This is obtained by imposing opposite twists for opposite spins: c r x ,r y +L y ,↑ = c r x ,r y ,↑ exp(i2πθ y ) c r x +L x ,r y ,↓ = c r x ,r y ,↓ exp(−i2πθ x ) c r x ,r y +L y ,↓ = c r x ,r y ,↓ exp(−i2πθ y ) with θ x = −1/2 + (i − 1/2)/N T and θ y = −1/2 + ( j − 1/2)/N T for integers 1 ≤ i, j ≤ N T , while r x and r y here and henceforth indicate integer Cartesian coordinates of the lattice 1 ≤ r x ≤ L x , 1 ≤ r y ≤ L y . All the results are then averaged on a mesh of N T × N T twists in the Brillouin zone (BZ), with N T large enough to have converged energies within statistical errors. When H MF conserves the number of particles N and has a gap in the one particle spectrum, N remains unchanged for each twist, whereas when a BCS pairing is present, grand canonical ensemble is adopted as discussed in Ref. 45 and the expectation value of H − µN is minimized by the proposed variational ansatz, whereN = ∑ r x ,r y ,σ c † r x ,r y ,σ c r x ,r y ,σ is the particle number operator. The use of opposite twists for opposite spin electrons is particularly important in this case because it allows the conservation of the translation symmetry in the BCS mean-field Hamiltonian, at least in a simple way.
Before explaining how to optimize the variational parameters it is useful to appreciate in Fig. 1 the fast and smooth convergence of the grand potential Ω = ψ τ |H−µN|ψ τ ψ τ |ψ τ in the thermodynamic limit as a function of the number N s of sites for a value of the chemical potential corresponding to doping δ = 1 − N/N s ≈ 1/8. The other important ingredient for an efficient implementation of Eq. (1) is the use of a particularly suited Trotter decomposition for the corresponding propagator: where in principle h i and t i can be independently optimized to minimize the Trotter error, as proposed in a recent work [36].
Henceforth it is assumed that the operators in Eq. (4) are ordered from left to right according to increasing values of the integer i. In this work τ is kept fixed, and therefore the following constraint is imposed: Moreover, in order to minimize the number of variational parameters and the number n, the following parameterization is adopted: with h 0 = 0 and γ > 1 chosen in a way to satisfy the constraint in Eq. (5). This expression is based on the conventional small ∆τ symmetric Trotter decomposition exp , where a nonuniform time step according to Eq. (6) has been chosen because, for high accuracy, it is important only that the first time step is small and γ 1. In this way the efficiency and the number of operators involved is significantly optimized, without using too many variational parameters as in Ref. 36. In all forthcoming calculations only a single parameter ∆τ is optimized with the convenient choice: leading to a Trotter time step error on the energy that, for all U values reported, can be hardly distinguished from the error bars. Here and henceforth the square brackets indicate the integer part of a real number.

III. VARIATIONAL AUXILIARY FIELD METHOD
Since H MF (α) is a one body operator and |ψ MF is a meanfield state (e.g. a Slater determinant) there is no difficulty to apply exp [−t i H MF (α)] to |ψ MF . Moreover, taking into account that exp [−t i H MF (α)| ψ MF remains a mean-field state, this operation can be also performed iteratively.
The application of exp(−h i V ) is instead more complicated and can nevertheless be implemented by using the auxiliary field technique. To this purpose the following discrete Hubbard-Stratonovich transformation [18] is used: where m j = n j,↑ − n j,↓ ,N is the total number of particles operator, cosh λ i = exp(Uh i ) and the Ising variables σ ji = ±1 for each site j of the lattice. Notice that the factor exp(−Uh iN /2) represents only an irrelevant change of the chemical potential in H MF (α). Therefore this term is omitted for simplicity here and henceforth.
We get therefore that the expectation value of Hamiltonian H (or any relevant correlation function [46]) for the wave func-tion in Eq.(1) is given by: where σ indicates the N s × n dimensional vector with components σ ji , and: E n can be therefore computed by Monte Carlo, by sampling the Ising fields σ and σ according to the weight |W n (σ , σ)|, where: and, quite generally in the complex case, W n has a non trivial phase that is determined by S n (σ , σ) = W n (σ ,σ) |W n (σ ,σ)| , a complex number with unit modulus, that plays the role of the QMC sign.
Finally E n can be computed by : namely by evaluating the ratio of the means corresponding to two real random variables ℜ[e n (σ , σ)S n (σ , σ)] and ℜ[S n (σ , σ)] [47] over the configurations generated by Monte Carlo sampling according to the probability p n (σ , σ) = |W n (σ ,σ)| ∑ σσ |W n (σ ,σ)| , using the standard technique described in Ref. 28, and: is a sort of local energy, namely an estimate of E n for a given configuration of the Ising fields σ and σ . Indeed both e n (σ , σ) and S n (σ , σ), as well as W n (σ , σ) can be computed in ∝ N 3 s n operations because they involve essentially imaginary time propagations of mean-field states under time dependent one body propagators U n (σ) and U n (σ ).
The basic ingredient introduced in this work is the possibility to compute energy derivatives of E n with respect to all the parameters defining the WF, i.e. α and α , henceforth assumed to be defined by the 2p variational parameters α 1 , α 2 , · · · α 2p and the minimum time step used ∆τ = α 2p+1 at fixed τ, according to Eq.(5). A simple algebra, very similar to the one known for VMC, implies that any energy derivative ∂ E n ∂ α j with respect to an arbitrary variational parameter α j , for j = 1, 2, · · · 2p + 1, can be computed by means of corresponding derivatives of two complex functions ∂ e n ∂ α j , O j = ∂ ln(W n ) ∂ α j = ∂W n ∂ α j /W n of the local energy and the logarithm of the weight, respectively: where here and henceforth the symbol * indicates the average of the generic random variable * over the probability distribution p n defined before, and for shorthand notations the dependence on σ and σ of all the quantities involved is not explicitly shown. The differentiation of the complex quantities lnW n and e n , required for the ∂ E n ∂ α j evaluation, at given values of σ and σ , may appear very cumbersome and involved especially considering that, it is often necessary, as in VMC, to optimize several parameters. This task can be easily achieved in a computational time equal to the one required to compute the complex quantities ln(W n ) and e n , remarkably only up to a small prefactor regardless how large is the number of variational parameters involved. This is possible by using Adjoint Algorithmic Differentiation (AAD) [48], a technique that is becoming popular in the field of Machine Learning with another name, i.e. "back propagation", but was certainly known before in applied mathematics [49], and only recently has been appreciated in physics [50,51]. Once all energy derivatives ∂ j E are known, the usual scheme adopted in VMC can be applied also here. Variational parameters are changed according to the equation: where rate learning is a suitable small constant, determining the speed of convergence to the minimum and F is the so called Fisher-information matrix, given by: where the symbol AB = AB − A B indicates here the covariance of two random variables over the probability p n . We adopt here rate learning 6 and the same regularization with ε = 0.01 described in Ref. [28] to avoid instabilities in the calculation of the inverse matrix. Typically convergence is reached with a few hundreds iterations and, variational parameters are averaged after convergence for about 50 steps.

IV. TESTING THE METHOD
In this section the proposed method is tested against known benchmark results on infinite systems. Unfortunately there are only a few results available in the thermodynamic limit, and those ones are mainly limited to ground state energies. Nevertheless they are extremely relevant for a variational method, because its predictions can be supported by a good estimate of the energy.
First VAFQMC is tested at U/t = 8 on cylinders with finite width L y and periodic boundary conditions (PBC) in the short y− direction, in order to compare with the very accurate results determined by DMRG with open boundary conditions (OBC) in the x−direction, and extrapolated to the one dimensional L x = ∞ thermodynamic limit. As it is seen in Fig. 3 VAFQMC is in very good agreement with the known results for L y ≤ 6, especially when the T 2 eff → 0 extrapolation is employed in a relatively small L y × 16 cluster with 32 twists in the long direction [52]. In this picture we can also appreciate the importance to use appropriate boundary conditions to reach the thermodynamic limit, as OBC require very large clusters for this purpose, a clear limitation that prevents a very accurate extrapolation. Notice that VAFQMC extrapolated results, though affected by a visible non linear term, are obtained with variational energies that are better (below) than the available [53] DMRG ones also for the smallest width L y = 4. This is remarkable because DMRG has obviously best performances in almost one dimensional systems. The small discrepancy between VAFQMC and DMRG extrapolated values for L y = 6 does not necessarily imply that VAFQMC is less accurate in this case, but may be due to the limitation of DMRG in obtaining reliable extrapolations as soon as we approach the 2D case, as also shown by the large error bar reported. In any case it is clear that the proposed extrapolations, can be used to estimate the quality of the VAFQMC best (lowest T eff ) variational energy estimates.
In all these VAFQMC calculations, at doping δ = 1/8 it has been found that the W = 8 stripe is the most favorable meanfield, as shown in Fig. 2 for the L y = 4 case. Here this state is parametrized by the most general mean-field Hamiltonian with local and nearest neighbor couplings independent of r y (see also the Supplementary Information): where n ↑,r x ,r y = c † ↑,r x ,r y c ↑,r x ,r y , n ↓,r x ,r y = c † ↓,r x ,r y c ↓,r x ,r y , n r x ,r y = n ↑,r x ,r y + n ↓,r x ,r y . Thus the local magnetic antiferromagnetic field satisfies ∆ AF (r x ) = −∆ AF (r x +W ) and the corresponding local chemical potential µ(r x ) = µ(r x +W ). Moreover nearest neighbor hoppings t x (r x ) = t x (r x +W ) and t y (r x ) = t y (r x +W ) are also included, amounting to a total of 4W variational parameters defining H MF (α 0 ). All their values are assumed independent of the twists [45]. The few parameter choice adopted in Ref. 54 is used only to initialize the present most general 8W +1 (by including also the optimization of H MF (α) and ∆τ) parameter ansatz |ψ τ .
For the uniform solution (see Sec. VI) antiferromagnetic order is allowed only at half filling while at finite doping a four (three) parameter ansatz is adopted in H MF (α) (H MF (α 0 ), here t = 1 can be left unchanged, as it sets the scale of the mean-field energy, irrelevant for |ψ MF ) including nearest and next nearest neighbor hoppings as well has a uniform chemical potential µ 0 and d x 2 −y 2 pairing. In all the forthcoming calculations, when the particle number is not defined neither in the mean-field wave function nor in the projection, the energy per site e(δ ) at fixed doping δ is accurately estimated by an appropriate choice of the chemical potential µ. This is obtained by simple and stable interpolations with a few calculations in the grand canonical ensemble or, in other words, by inverting the Legendre transform from chemical potential dependence to the conjugate density one 1 − δ .
With this method it is therefore possible to compute the energy in the thermodynamic limit of true 2D clusters without particular effort as the average (complex-) sign S n is always larger than 0.3 for all the simulations reported in this work. It is important to remark that, when the 2Dthermodynamic limit is approached, a small but systematic increase of the energy can be appreciated, confirming that, in order to determine the correct two-dimensional thermodynamic limit, true two-dimensional clusters have to be used.
For L y small the stripe solution is clearly favored as compared with the uniform solution, as found by DMRG. In Fig. 3 -0. the uniform solution is optimized with a non zero d x 2 −y 2 BCS pairing but remains clearly above the stripe solution even after extrapolation to T eff = 0. However the situation is quite different when the twodimensional thermodynamic limit is required. In this case, in order to compare the stripe with periodicity W = 8 and W = 7 and the uniform d − wave solution, calculations on a 16 × 16 and 14 × 16 clusters are carried out with 16 × 16 different twists in the BZ. For the d−wave ansatz several calculations at different chemical potentials µ 1.8t are attempted in order to fulfill a doping δ = 1/8 for each effective temperature T eff .
It is seen from the Fig. 4 that in 2D the uniform solution turns out to have an energy competitive with the stripe ones, though the W = 7 stripe appears to be the most stable solution also because it provides always the lowest variational energy for each T eff . Nevertheless the stripe-uniform phase energy difference, reached at the lowest T eff shown, is very small 0.001t, though the two variational states remain with qualitatively different types of correlation functions, i.e. a sizable magnetic moment in the former and non negligible d − wave pairing correlations in the latter (see inset in Fig. 4). In Tab. I the various extrapolated and best variational energies are compared with the available benchmark results [20,21,54]  ods, because DMRG becomes rapidly inaccurate with increasing L y , and the best variational energy obtained by iPEPS is significantly higher than the VAFQMC one. This may explain why the iPEPS extrapolated energy has a large uncertainty and is much below the VAFQMC and DMRG 6×∞ ones. It is also possible that the extrapolation in the bond dimension is quite inaccurate, and can be substantially improved [55]. Finally the energy and magnetic order parameter are compared at half-filling. In this case there is no sign problem and in principle large imaginary times can be employed without particular difficulties. However, τt ≤ 1.3 has been chosen also in this case, in order to show the strength of this approach even when short time projections are employed. For the order parameter finite effective temperature scaling anal-ysis implies convergence linear in T eff much slower than the energy (see App.A). Nevertheless remarkably accurate results can be obtained in both cases, as shown in Fig. 5. At finite effective temperature T eff very weak finite size effects are seen, showing once more the great advantage to use TABC with this approach. The error in the T eff → 0 extrapolated energy is then compatible with the reference one within the error bars (this work extrapolation E h = −0.52443 (14), reference −0.5247 (2)) whereas the extrapolated order parameter (0.2894 (16)) is in very good agreement with the benchmark one (0.2991 (2)).

V. THE HARTREE-FOCK STRIPE PHASE
Also within the Hartree-Fock, namely by considering a simple Slater determinant or also uncorrelated BCS mean-field wave functions, very few established results are known, if we allow also nonuniform solutions [58,59], i.e. non translation invariant solutions with large unit cell. At half filling it is clear that the simple antiferromagnetic solution, with two site per unit cell, is stable for any U > 0 as it has been shown analytically by Hirsch [18]. On the other hand it is well-known that no superconducting solution is possible for U > 0. Restricting the solution to a small unit cell containing only two sites leads to a phase separation instability [60], namely the energy per hole: acquires a minimum at a finite doping δ c and for any doping δ ≤ δ c it is more convenient to expel the holes from a pure antiferromagnet with energy per site e(0) in a region with appropriate size containing an hole rich phase at doping δ c . After this construction (see e.g. App.B) e h (δ ) will be constant in the thermodynamic limit and equal to e h (δ c ) for all dopings δ ≤ δ c . This represents a fingerprint of phase separation and the study of the energy per hole can be also considered an implicit way to build a nonuniform solution (when a minimum is found) within not only HF but any variational approach based on a translationally invariant ansatz. However a better way to expel holes from an antiferromagnet is found by means of the stripe solution already reported in Fig. 2. In particular when the doping δ = 1/W the HF bands show a clear insulating behavior because the unit cell 2×2W contains 4W doubly-degenerate bands [61] and 2W −2 are fully occupied, as is the case for δ = 1/8 shown in Fig. 6. The insulating nature of this stripe solution was also pointed out in weak coupling HF theory by H. Schulz [58] who discovered first the stability of incommensurate magnetic states with finite wavevector Q = (π ± πδ , π) close to the antiferromagnetic one Q = (π, π), and the opening of a full gap away from half filling within HF. The occurrence of a finite gap in this case is easily understood not only because the mean-field HF Hamiltonian turns out to have a gap, but also for the following simple general argument holding also in the correlated case. This incommensurate state is adiabatically connected to the insulator having equally spaced empty (i.e. with no-electrons) vertical lines of sites separating half-filled antiferromagnetic insulating regions.
In this way the HF solution can avoid phase separation but the corresponding energy per hole is almost flat at small doping (see Fig. 7a) that is almost equivalent to phase separation. Indeed at small doping the stripes are very far apart W = 1 δ and do not interact at all, as can be appreciated in Fig. 7b where the interaction between two vertical stripes at distance W is given by I(W ) = e h (δ = 1/W ) − e h (δ = 1/∞). I(W ) here defined represents the energy cost per hole for two stripes being at finite distance W rather than at infinite distance, where they do not interact.
At non commensurate doping, for instance 1/8 < δ < 1/7 it is possible to verify that such kind of insulating solutions at δ = 1/7 and δ = 1/8 can be joined together by forming a smooth doping dependent insulating stripe phase at any intermediate doping (see for instance Fig.7a for U/t = 8). Many stripes at positions r x i of the lattice can be thought to interact by means of a pairwise repulsive interaction ∝ ∑ i< j I(r x i − r x j ), with interaction I(W ) as the one computed in Fig. 7b. Thus the rule for obtaining the minimum HF energy is as follows: in order to get the appropriate hole density δ one alternates energetically more expensive defect stripes of length W smaller than [1/δ + 1] (W = 7 in this case) by placing them as far as possible.

VI. RESULTS FOR THE CORRELATED PHASE DIAGRAM
The rule determined in the previous section to identify the minimum energy insulating stripe solution has been verified also in the correlated case up to the maximum τ possible with essentially no sign problem, namely τt = 1.3, 1.73, 2.6 for U/t = 8, 6, 4, respectively. Therefore first the stripe wave function (breaking translation and spin symmetry) is computed at commensurate doping δ = 1/W , where W , being the distance between equally spaced stripes, is an integer. Then, in order to account for intermediate dopings, the corresponding energies are interpolated. Moreover, in order to estimate the error to the exact T eff = 0 limit, the results are extrapolated to this limit. Hence the lowest energy of the uniform solution is estimated, that is parametrized by the following mean-field Hamiltonian: whereK = ∑ i, j,σ t R,R c † R,σ c R ,σ is the translation invariant kinetic energy defined here by the nearest and next nearest neighbor hopping t and t , namely t R,R = −t (t R,R = −t ) if R and R are (next) nearest neighbor sites, ∆ x 2 −y 2 determines the gap function for a d x 2 −y 2 BCS superconductor, ∆ AF is the meanfield gap due to a commensurate antiferromagnetic state, and µ 0 is the mean-field chemical potential value.
By the proposed energy optimization method, as the doping is decreased clear evidence of an instability towards d x 2 −y 2 BCS pairing is found, because as soon as ∆ x 2 −y 2 turns out to be non zero in the thermodynamic limit, the variational wave function |ψ τ , breaks the U(1) symmetry c R,σ → exp(iθ )c R,σ of the Hamiltonian, acquiring the best possible energy for a uniform superconducting phase. In this approach ∆ AF turns out to be non zero only at half filling but it is important to emphasize that, despite the uncorrelated HF case, a non zero d x 2 −y 2 BCS pairing is possible at finite doping for T eff lower than a critical effective temperature, and this represents one of the most important effect of electron correlation, as it is not present in the infinite T eff HF theory.
In this way we can determine the transition between the translationally invariant phase and the stripe phase as reported in Fig. 8 for the three representative values of U/t. Notice also that the location of the transition points does not depend much on the extrapolation, clearly supporting that the phase diagram converges quite fast by lowering T e f f .
As in HF theory, away from the commensurate dopings 1/W , an insulating phase covering a continuum of different dopings can be constructed by appropriately joining such commensurate solutions, i.e. a doping 2/13 with insulating properties can be easily obtained by alternating a W = 7 stripe with a W = 6 one. This is actually the case in a 26 × 12 cluster, within VAFQMC, namely the 6 + 7 + 6 + 7 alternating stripe length solution has an energy slightly lower than the corresponding one 6 + 6 + 7 + 7, thus satisfying the HF rule stated at the end of the previous section. This phase, competing with the uniform d x 2 −y 2 BCS one, has a full charge gap in H MF (α 0 ), a property that cannot be changed by a small τ−projection. Moreover no coexistence of stripe and BCS order was found within VAFQMC.
In summary, according to Fig. 8, it is possible to conclude that an insulating stripe phase acquires an energy lower than the corresponding uniform phase already at quite large doping around δ 20%. Nevertheless, it is important to remark that the uniform phase has a very good energy, competing with the lowest possible ones. Indeed in Fig. 4 the uniform d − wave phase looses only a tiny energy (less than 0.001t per site at the lowest T eff ) vs the nonuniform stripe phases. Interestingly just around the doping δ = 1/8 equally spaced stripes with length W < 1/δ appear lowest in energy, as can be seen in Fig. 9, where the energy of such metallic (at least within the mean-field Hamiltonian H MF (α 0 )) stripes are compared with the commensurate insulating ones (black lines). For instance the W = 9 stripe at doping around 1/13 (i.e. a metallic stripe with W < 1/δ ) has an hole energy more than 0.01t below the corresponding insulating stripe with W = 13, a huge gain that can be hardly explained by artifacts of the approximations (finite T eff and/or finite clusters). In this plot we also see that phase separation (a minimum of the hole energy) cannot be excluded for doping δ 5% but the results strongly depend on the extrapolation and requires therefore more accuracy to be settled. By judging from the extrapolated values it is plausible to expect a similar flat behavior of the energy per hole at low doping as in the corresponding HF plot (Fig.7a). However, within HF, a metallic solution with W < 1/δ was never found to have the lowest energy at low doping, again in agreement with Ref. 58. In the following it is worth to discuss rather extensively the main result of this work, namely that, at large doping the stripe melts and a small but sizable superconducting d-wave order clearly remains within the present approach, in the sense that the energy is lowered in the thermodynamic limit by breaking the U(1) global symmetry related to number of particle conservation. This is because both the mean-field Hamiltonians H MF (α 0 ) and H MF (α) clearly support a broken symmetry solution of this type when ∆ x 2 −y 2 > 0 in Eq. 21. Indeed, in Fig. 10 the superconducting gap ∆ x 2 −y 2 parameter, corresponding to H MF (α 0 ), is displayed for U/t = 8 as a function of the doping. When the ∆ x 2 −y 2 is large there are negligible size effects. However when it approaches zero, a long slowly decaying tail shows evident size effects. Indeed with the largest affordable sizes it is possible that the calculated critical doping δ c above which a symmetric phase is stable indicates only a sharp crossover region separating a phase with sizable strong coupling superconductivity from another phase with very small order parameter, compatible with an exponentially small pairing of the Kohn-Luttinger type [24,62]. Indeed, within weak coupling theory, the d − wave superconducting phase should remain stable up to 40% doping, a doping much larger than the one detected in the present work.
As shown in Fig. 11, within this approach, the critical doping characterizing a strong coupling superconductivity (i.e. a non negligible pairing) can be identified because for µ ≤ µ c the energy clearly improves by increasing the value of ∆ x 2 −y 2 even when starting from a negligible value, representing the symmetric Fermi liquid ground state. Nevertheless there may be non negligible size effects as shown for U/t = 8, especially considering the difficulty to locate the end of a long tiny tail of the order parameter, as discussed before. Moreover, for smaller U these size effects are expected even larger, but much larger cluster simulations are prohibitive at present. For this reason for U/t = 2 (U/t = 4) a much larger mesh for the TABC was used with 64 × 64 (32 × 32) different twists in the BZ, instead of the smaller 16 × 16 grid used for U/t = 8 and U/t = 6. In a mean-field calculation a 64 × 64 TABC grid in a 16 × 16 square lattice corresponds to a 1024 × 1024 cluster with PBC, large enough to probe even a very tiny but non negligible gap parameter ∆ x 2 −y 2 . Remarkably, at U/t = 2 no evidence of such a tiny value of the order parameter is found in all relevant doping region, as can be clearly seen in Fig. 11.
This small but relevant strong coupling region where dwave superconductivity appears to be stable in the Hubbard model is in agreement with the old findings by cluster dynamical mean-field theory [23], though it is located at a slight  higher doping, namely at 10% the stripe should be stable also for U/t = 4. Moreover it is worth mentioning a recent VMC calculation [63] employed at U/t = 10 and density matrix embedding theory(DMET) [64] supporting also the stability of a d−wave superconductivity in a small doping region. On the other hand the proposed phase diagram is also in agreement with the claim [25] of absence of d−wave superconductivity in the 2D Hubbard model, because in this work the authors refer mainly to doping 1/8, where no BCS pairing was found also in this work. The final phase diagram is therefore reported in Fig. 12.

VII. CONCLUSIONS
In this work the VAFQMC method, able to exploit the power of the auxiliary field technique combined with the sim-  For all cases the initial wave function was prepared with ∆ x 2 −y 2 = 0 defining the Fermi liquid state. Then this parameter was set to a small value and the iterative stochastic optimization was started. In these calculations, in order to reduce at most the stochastic bias, at each iteration a quite large number (> 3 × 10 5 ) of samples for evaluating energy derivatives was used, in a way that they are sensitive to very small energy differences < 10 −5 t. Optimization of the grand potential Ω at fixed chemical potential µ was employed in all cases shown. The corresponding doping was estimated by ∂ Ω ∂ µ = −(1 − δ )N s using the data for Ω reported in the supplementary material SI extrapolated to the T eff → 0 limit. These estimates are consistent with the direct evaluation of the doping δ = 1 − ψ τ |N|ψ τ N s ψ τ |ψ τ at the lowest T eff within half the confidence doping interval reported for each U/t = 2. plicity and generality of the standard variational quantum Monte Carlo method, has been introduced. By combining these two successful approaches to strongly correlated systems, it is possible to estimate the evolution of the phase diagram of a lattice model, by improving systematically the accuracy of the correlation term, starting from the HF approximation. This is achieved by lowering an effective temperature T eff up to a value that turns out to be low enough to provide very accurate variational energies of the Hubbard model in the thermodynamic limit, while remaining far from any sign problem instability, so far affecting most fermion quantum Monte Carlo techniques. Within the assumption that the tiny energy gain between different broken symmetry phases of the model can be sorted out by a careful optimization of a mean-field ansatz in presence of an accurate enough correlation term, the phase diagrams obtained in this way should be considered reliable, at least for the phases studied.
It is clear that the above assumption is very important as, for instance in the Hubbard model, we have estimated an energy per site error of about 2 × 10 −3 t (4 × 10 −4 t) at U/t = 8 (U/t = 4), while energy differences of the various phases at the lowest T eff considered can be even more than an order of magnitude smaller. This is not a difficulty of the present method, providing state of the art variational energies (see supplementary information for detailed benchmark results), but clarifies once more the enormous challenge of the elec- . The blue line is also determined by the condition that for U/t = 2 no BCS pairing has been detected, while the green line by the one that the incommensurate spin order cannot show up at non zero doping for U → 0 [58]. In the stripe phase there may be a transition between an insulating stripe and a small doping metallic stripe, as suggested by Fig. 9 for U/t = 8. In this transition the mean-field VAFQMC Hamiltonian H MF (α 0 ) has a finite gap (bands are completely filled) or a vanishing one, respectively.
tron correlation problem in numerical calculations. For instance on such clusters it is basically impossible to distinguish an unbroken symmetric metallic phase from a d−wave superconducting one with small ∆ x 2 −y 2 gap parameter, by judging merely on correlation functions (e.g. pairing correlations), and the above assumption provides a much more sensitive criterium for occurrence of broken symmetry phenomena. Relying on the above assumption, dramatically helps numerical techniques, such as VAFQMC, able to work with boundary conditions fulfilling the symmetries of H, that may be possibly broken in the thermodynamic limit.
In any case, the present work provides rigorous and accurate upper bounds of the ground state energies of the Hubbard model for several U/t and doping values well converged in the thermodynamic limit (see supplementary information), that can be useful in future works on this important subject and/or for benchmarking new computational, theoretical or experimental techniques. The fact that extrapolated energies do not depend much on the ansatz adopted (all extrapolated energies agree within 0.001t per site) implies that energetic properties are settled within a reasonable accuracy. However, in order to get the right order without assuming it in the initial mean-field ansatz, much larger projection times are required (see for instance Fig.4 where two completely different phases differ in energy by less than 0.001t per site), that are at present not possible.
In this study many different mean-field solutions have been attempted with several different initializations of parameters, including also broken time reversal solutions [65] and different modulations of stripes and pairing. In all cases the smooth convergence of the energy as a function of T eff has been verified. This is not only important for achieving reasonable extrapolations but especially for excluding being trapped in spurious local minima, that is a limitation of any non linear optimization technique. The phase diagram presented here in Fig.12 shows the phases so far determined with the lowest possible variational energies among the chosen ansatzs.
An important progress in this study is the control of finite size effects within VAFQMC. The conventional approach is to attempt exact calculations (usually very hard if not impossible) at finite number of sites N s and try to extrapolate to 1 N s → 0. Within VAFQMC it is much simpler and indeed possible to extrapolate only in the effective temperature T eff → 0 limit the finite T eff t/2 results. To this purpose it is enough to perform N s 100 calculations because they are already very close to the 1 N s → 0 limit, thanks to the very effective twisted average boundary condition method.
The results shown in this work are weakly affected by the sign problem. In general the average sign, appearing for instance in the denominator of Eq. (16), should decay exponentially to zero both with 1/T eff and N s , yielding prohibitive statistical errors. However for T eff large enough (i.e. the ones used here) the exponential decay in N s of the average sign is very weak, due to a stability property of the adopted auxiliary field technique, that was discovered several years ago [66]. Hence much larger clusters could be safely simulated, that is very promising considering also that the ones reported here are not far from the state of the art simulations not vexed by the "sign problem".
In this work a large stripe phase in the Hubbard model exists as an effective way to expel the holes from a clean antiferromagnet. In this way genuine phase separation does not occur in the Hubbard model as a result of a very weak repulsive coupling between stripes at large distance. In practice, at small doping, the energy per hole is almost constant, namely the inverse charge compressibility is vanishingly small, in agreement with the phase separation scenario [60] and some previous numerical calculations [67][68][69].
An important outcome of the present work is that a sizable superconducting phase is present in the Hubbard model with a non negligible order parameter only in the strong coupling regime. For instance for U/t = 2 no evidence of d x 2 −y 2 BCS order was probed and the maximum ∆ x 2 −y 2 0.02t was obtained with the largest U/t = 8 considered at about 20% doping. The critical doping δ c separating the metal from the d−wave superconductor turns out to be much smaller than the one obtained with weak coupling analysis [24,62] and the question is how to reconcile the weak and the strong coupling limits. In principle, as suggested by the results in Fig. 10, a small tail with very weak BCS pairing could distinguish a strong coupling phase with sizable pairing and a weak coupling phase with a sharp crossover (and not a transition) clearly separating the two regimes.
The main features of cuprate superconductivity can be ex-plained with the Hubbard model, but quantitatively there remain several unsolved issues. First, the stripes are completely filled by holes and do not match the half filled ones detected experimentally [70]. Second, the d−wave superconductivity is sizable but still far to explain T c 100K because the gap function, that is growing when approaching the Mott insulating phase at zero doping (see Fig.10), is limited by the phase transition to the stripe phase. VAFQMC has been developed to study properties of strongly correlated systems in the thermodynamic limit. Obviously this method applies also to finite systems for obtaining the best ground state estimates with the largest possible projection time τ. However for finite systems, the symmetry is never broken, and a more accurate wave function should be used by restoring the symmetry with appropriate projection operators [71,72], an improvement that can be done and should be worth doing for this particular purpose. Moreover, several new applications and extensions are possible within VAFQMC. The WF |ψ τ can be further generalized to more realistic Hamiltonians, including for instance long range Coulomb interaction or electron-phonon coupling. The finite effective temperature T eff here defined for a variational wave function is the key to obtain very accurate results well converged in the thermodynamic limit, a feature that could be possibly generalized to several other computational techniques, from tensor networks to machine learning variational wave functions.

ACKNOWLEDGMENTS
Acknowledgments: I am particularly in debt with A. Parola, N. Costa, and E. Tosatti for several suggestions for improving the manuscript. I am also tgrateful for very useful discussions with G. Carleo, N. Costa, A. Tirelli, F. Becca, L. Capriotti and M. Fabrizio. I acknowledge PRACE for awarding access to Marconi at CINECA Italy, and Riken collaboration for access to HOKUSAI supercomputer in Saitama Japan, as well as financial support from PRIN 2017BZPKSZ and the Simon's foundation.

Energy and order parameter
The expectation value of the energy E(τ) over the state |ψ τ = exp(−Hτ/2)|ψ MF can be written as: Since E(τ) in the expression (A6), depends on τ only by means of ω n , it follows that: On the other hand the local magnetization:

Evaluation
As in any field theory the above expressions diverge without introducing infinite counter terms. In the expression of the energy one can subtract and add one in the integrand, whereas in the correction of the magnetic moment it is enough to deal with particular care the q = 0 mode (this mode may be also avoided by using twisted boundary conditions): The infinite RHS term in the energy corresponds to an infinite shift of the ground state energy, a typical feature of quantum field theories, as this shift becomes finite by introducing a cutoff, implicitly present in any physical model. In order to evaluate the above expressions it is useful to consider the following sum: that is valid for any complex number z. In particular for z = whereas for z = iy we get: Defining the velocity of the Goldstone's modes c = ϒ χ , the sum over n is extended to −∞ ≤ n ≤ ∞, using the even dependence on n and taking into account the extra n = 0 contribution in the definition of ε 0 . In this way, by using Eq. A15 for evaluating the infinite sums the universal finite size corrections are given by: where ε 0 is the "infinite" ground state energy corrected for the zero point harmonic energy a divergent expression depending on the cutoff but not on the boundary conditions (i.e. it is the same divergent term appearing in Ref. 40, by using the trace). In the thermodynamic limit and fixed τ the convergent contribution can be worked out, by replacing sums over momenta with integrals (i.e. 1 V ∑ q → ( dq 2π ) D ), and changing the momentum scale q → cτ q, so that: where f e is a dimensional dependent constant defined by: that is clearly a convergent integral for any dimension D ≥ 1 because coth(|q|) converges to one exponentially for large |q|. The values of f e are therefore π 12 , ζ (3) 4π (ζ (3) = 1.20205690316) and π 2 240 in D = 1, D = 2 and D = 3, respectively. In D = 1 the coefficient f e has not to be taken for granted, since no true long range order is possible (see later), but the qualitative behavior of the correction T 2 eff should hold as it is consistent with the corresponding finite temperature correction derived in D = 1 for the 1D-Hubbard model by means of the Bethe Ansatz [74]. Thus, in the thermodynamic limit, the finite τ corrections to the energy depend only on the universal constant f e , the number of components M of the order parameter, and the velocity c of the gapless Goldstone modes. In D=2 the finite temperature corrections to the internal energy at infinite volume have been determined in Ref. [41]. For this purpose we have repeated the same type of calculation and get that, taking into account the different boundary condition for the field (PBC in time for the trace τ = β = 1/T ): that is consistent with the previous calculation reported in D=2, i.e. [41]: On the other hand, for the magnetic order parameter, by noticing that sin( ω n τ 2 ) 2 = 1 for n = 2ν + 1 and zero otherwise, the sum has to be evaluated for odd integer frequencies ω ν = π(2ν+1) τ . For this purpose the sum ∞ ∑ ν=0 1 χω 2 ν +ϒ| q| 2 is extended to 1 2 ∞ ∑ ν=−∞ 1 χω 2 ν +ϒ| q| 2 because ω ν = −ω −(ν+1) for ν < 0. Therefore by applying Eq. A14 for evaluating the following infinite sums, one obtains: 1 |q| (the first term in ∆m 0 is the q = 0 contribution in Eq. A12 that is convergent for finite τ) and the dimensional dependent constant f m is defined by: that converges in D > 1. The values of f m are therefore − ln(2) π and − 1 12 for D = 2 and D = 3, respectively. Thus, in the thermodynamic limit the finite T eff corrections to the antiferromagnetic order parameter depend on the universal constant f m , the number of components M of the order parameter, the velocity c and the stiffness ϒ of the gapless Goldstone modes.
Notice that the expression above for f m converges in D = 2. Instead if periodic boundary conditions in τ are assumed for the field Π(x,t), corresponding to the standard finite temperature calculation at β = τ, by repeating the same steps, it follows that: ∆m(τ) = − M−1 coth(c|q|τ/2) |q| which blows up for D ≤ 2, thus recovering the Mermin-Wagner theorem [75]: no finite m is possible at finite temperature for D ≤ 2. The advantage of using the projection technique is therefore evident especially in D = 2 for the study of broken symmetry phases, that are possible, for continuous symmetries, only at zero temperature.
When using finite cylinders with even L y all the above results obtained for D = 1 apply because the y momentum q y = 2πn y L y is quantized with n y = 0, 1, · · · , L y − 1 and only the q y = 0 momentum value provides power law corrections in 1/τ to energy and correlation functions (if converging). This is because all the other contributions acquire a finite ∼ 1 L y gap and converge much faster.
Appendix B: Convexity of the energy as a function density One can divide a large system A + B in two regions A and B containing a macroscopic number of sites L A and L B and electrons N A and N B , respectively. Then, by defining E A+B the energy of the total system, E A and E B the ones of the corresponding parts, ρ = N A +N B L A +L B = pρ 1 + (1 − p)ρ 2 , ρ 1 = N A L A , ρ 2 = N B L B , p = L A L A +L B , e(ρ 1 ) = E A /L A and e(ρ 2 ) = E B /L B converged in the thermodynamic limit because the two parts are macroscopic, then e(ρ) = E A+B L A +L B ≤ E A +E B L A +L B = pe(ρ 1 ) + (1 − p)e(ρ 2 ) follows because one can neglect in the thermodynamic limit the surface term contribution to the energy separating the region A from the region B. This is because, by assumption, the model is short range as the hopping term connects only nearest neighbor sites and the interaction U is on site. The final inequality: e(ρ) ≤ pe(ρ 1 ) + (1 − p)e(ρ 2 ) (B1) therefore holds for arbitrary densities ρ 1 ≤ ρ ≤ ρ 2 and implies, the convexity property of the function e(ρ) in the thermodynamic limit. Whenever phase separation occurs between two densities ρ 1 and ρ 2 the inequality B1 turns to a strict equality, namely e(ρ) is a linear function of the density because p = ρ−ρ 2 ρ 1 −ρ 2 , so that in particular the energy per hole in Eq. 20 is constant for δ ≤ δ c .