  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.
