The two-dimensional Hubbard model at finite doping hosts competing or intertwined orders, resulting in conflicting conclusions from different computational approaches regarding its ground state. We show that a key source of such discrepancies is the bias encoded in the variational ansatz. We consider three different Transformer backflow fermionic wave functions based on a Slater determinant, its particle-hole counterpart, and a Pfaffian, initialized without any mean-field pretraining. We show that, despite achieving nearly degenerate, state-of-the-art variational energies, each ansatz converges to a state with qualitatively different spin, charge, and pairing correlations. Upon improving accuracy via symmetry restoration and variance reduction, however, all three converge to the same physical picture: coexisting superconducting and stripe orders. These results demonstrate that variational energy alone is insufficient to identify the ground state in the presence of competing phases, and highlight the importance of tracking how correlation functions evolve as the wave function is systematically improved before drawing physical conclusions.
The two-dimensional Hubbard model was originally proposed as the minimal model of uniform
d-wave superconductivity in the cuprate high-temperature superconductors, but recent calculations with relevant parameters in a realistic doping regime have often found other ground states. Further, uncertainty remains over such results due to the intrinsic bias of different computational methods. Here we address this uncertainty through state-of-the-art neural tensor network and neural backflow state simulations, enabled by large-scale GPU computation, that together reach an unprecedented variational accuracy on lattices with up to 1024 sites. At
1/8 doping in the pure Hubbard model, as well as with next-nearest-neighbor hopping relevant to the hole-doped cuprates, we reveal a near degeneracy among a low-energy uniform
d-wave state, superconducting stripes, and non-superconducting stripe states on the milli-
t energy scale. Our work rebalances the narrative that the Hubbard model supports a single dominant order in relevant parameter regimes towards a picture of extreme competition between stripes and superconductivity.
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 already be possible in the Hubbard model for relatively small square lattices with linear dimension $L\ensuremath{\simeq}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.
We apply the Tensor-Backflow method to investigate the Fermi-Hubbard model on two-dimensional lattices up to 256 sites, exploring various interaction strengths
U, electron fillings
n, next-nearest-neighbor hopping
t′, and boundary conditions. By considering backflow terms from nearest- or next-nearest-neighbor sites, we achieve competitive results without enforcing geometric symmetries on the variational wave-function. The optimizations were stable from a prior unrestrictied Hartree-Fock state, followed by adding backflow corrections. Meanwhile, changing interaction strengths in the prior unrestrictied Hartree-Fock state is helpful to bypass the local minima. When
t′=0, by considering nearest-neighbor backflow terms, linear stripe order emerges successfully for the case of
n=0.875 and
U=8 on a
16×16 lattice with periodic boundary conditions. In a similar case with open boundary conditions, the energy obtained is only
4.5×10−4 higher than the state-of-the-art method fPEPS with bond dimension
D=20. Compared to state-of-the-art neural network methods, the energies obtained using the Tensor-Backflow approach are competitive, with relative errors below
5×10−3. For
n=0.8 and
n=0.9375, direct optimizations yield results consistent with the phase diagram from AFQMC. When
t′=-0.2, considering next-nearest-neighbor backflow terms leads to energies that are either competitive with or even lower than those from state-of-the-art neural network approaches. For instance, for
n=0.875 and
U=8 on a
12×12 lattice with periodic boundary conditions, the energy obtained is
8.1×10−4 lower than that from the neural network result. Thus, the Tensor-Backflow method demonstrates strong representational capabilities for solving the Fermi-Hubbard model.
Building on the recently developed tensor-backflow method, which accurately describes energies and stripe states in the Fermi-Hubbard model, we investigate
d-wave pairing in the two-dimensional
t-
t′-
U Hubbard model. For the
8×8 lattice, we improve the wave function using symmetry projections followed by one Lanczos step, whereas for the
12×12 and
16×16 lattices, we apply one Lanczos step without symmetry projection. We focus primarily on the parameter regime with filling
n=0.875, on-site repulsion
U=8 and next-nearest-neighbor hopping
t′=−0.2t, a regime in which enhanced
d-wave pairing has been reported. The resulting energies are competitive with those obtained using state-of-the-art neural quantum states. For example, on a
16×16 lattice with periodic boundary conditions, the tensor-backflow wave function after one Lanczos step without explicit symmetry enforcement achieves a relative energy difference of
4.4×10−3 from the symmetry-preserving neural quantum state result. From the real space pair correlations, for
t′=−0.2t, we observe stronger
d-wave pair correlations together with the characteristic relative sign between horizontal and vertical bond-pair components. Compared with
t′=0, finite negative
t′ increases the spectral weight of the corresponding pair-density matrix, which indicates the enhanced overall
d-wave pair fluctuations. The leading eigenvalue is enhanced on the
8×8 lattice for both with and without symmetry projections. For larger lattices such as
12×12 and
16×16, the leading pair-density matrix spectral weight is distributed among several nearly degenerate eigenmodes rather than being concentrated in a single dominant mode. This indicates enhanced and more broadly distributed
d-wave pair fluctuations, but no single eigenvalue exhibits the extensive scaling required for ODLRO.
Simulating strongly correlated fermions remains a long-standing challenge due to the exponential complexity of the Hilbert space and the intricate sign structure of many-body wavefunctions. We introduce a variational framework centered on a neural network backflow transformation that combines deep learning with variational Monte Carlo. The proposed ansatz employs a multi-determinant expansion with low-rank shifts to capture non-local correlations and complex sign structures. Applied to the two-dimensional Hubbard model at both half-filling and
1/8 doping, the method achieves energies within
0.45% of auxiliary-field quantum Monte Carlo at half-filling and captures intertwined charge- and spin-density stripe patterns at
1/8 doping. These results demonstrate the potential of this framework as a scalable and interpretable approach to variational simulations of strongly correlated fermionic systems.