Intertwining of spin, charge and pairing correlations in the repulsive two-dimensional Hubbard model is shown through unrestricted variational calculations, with projected wavefunctions free of symmetry breaking. A crossover from incommensurate antiferromagnetism to stripe order naturally emerges in the hole-doped region when increasing the on-site coupling. Although effective pairing interactions are identified, they are strongly fragmented in several modes including d-wave pairing and more exotic channels related to an underlying stripe. We demonstrate that the entanglement of a mean-field wavefunction by symmetry restoration can largely account for interaction effects.Transition metal oxides are prototypical strongly correlated systems that exhibit a rich variety of quantum phenomena such as antiferromagnetism (AF), incommensurate charge and spin ordering or high-T c superconductivity. The understanding of their subtle competition at low-temperature remains one of the most challenging topics in condensed matter physics. According to Anderson's proposal [1], the two-dimensional (2D) single-band Hubbard model is expected to provide a minimal framework for addressing these issues in rare earth cuprates. It describes d-electrons hopping between the Wannier states of neighboring lattice sites in a copper-oxygen plane experiencing a purely local Coulomb repulsion. At half-filling, the interaction strength drives a Mott transition to an insulator [2] with long-range AF order [3]. When the lattice is doped away from half-filling, the exact form of the phase diagram is still controversial. Of central interest, in this regime, is whether the ground state supports unconventional fermion-pair condensates or charge inhomogeneities, and if so, how their order parameters are intertwined with magnetic properties. Only a partial answer can be currently obtained through approximate many-body techniques, such as cluster extensions [4,9] of the dynamical mean-field theory [5,9], the two-particle self-consistent approximation [6,9], Gutzwiller variational schemes [7,9] or slave-boson approaches [8,9]. Standard quantum Monte-Carlo simulations (QMC) are also restricted [10] owing to the notorious sign problem that is particularly severe for doped Hubbard models. Although new sign-free stochastic reformulations have been recently introduced [11,12], they are not, for now, immune to systematic errors [13]. To overcome these theoretical difficulties, a direct quantum simulation has been suggested by loading ultracold mixtures of two interacting fermionic species in optical lattices. Indeed, such atomic systems allow for an almost perfect implementation of the Hubbard model with tunable parameters [14]. The crossover from a metallic into a Mott insulating regime has already been observed [15], while magnetic ordering and potential exotic superfluidity remains to be achieved using, e.g., new cooling techniques [16].