Intertwined orders from symmetry projected wavefunctions of repulsively interacting Fermi gases in optical lattices

Unconventional strongly correlated phases of the repulsive Fermi–Hubbard model, which could be emulated by ultracold vapors loaded in optical lattices, are investigated by means of energy minimizations with quantum number projection before variation and without any assumed order parameter. Using a tube-like geometry of optical plaquettes to realize the four-leg ladder Hubbard Hamiltonian, we highlight the intertwining of spin-, charge-, and pair-density waves embedded in a uniform d-wave superfluid background. As the lattice filling increases, this phase emerges from homogenous states exhibiting spiral magnetism and evolves towards a doped antiferromagnet. A concomitant enhancement of long-ranged d-wave pairing correlations is also found. Numerical tests of the approach for two-dimensional clusters are carried out, too.


Introduction
Low-dimensional interacting quantum matter generally exhibits several phases at low energy that challenge the ability to distinguish between competing orders and their intertwining within one single correlated state [1]. Ultracold atoms provide an ideal playground to capture the essence of this problem by their potential to properly emulate the fundamental mechanisms of quantum many-body physics [2]. In the fermionic sector, the BCS-to-BEC crossover [3,4] and the question of Stoner's itinerant ferromagnetism in repulsive gases [5,6] have been investigated. By trapping atomic vapors in optical lattices, a mimic of ideal crystalline matter can also be achieved [7]. By now, direct images of Fermi surfaces in the non-interacting limit [8] as well as s-wave superfluidity near unitary scattering [9] have been reported. Away from a Feshbach resonance, one is able to engineer almost perfectly the celebrated Hubbard model that had been first considered to describe the magnetism of metallic systems [10]. More generally, it aims to grasp the generic properties of spin-1/2 fermions moving on a lattice by hopping between neighboring sites r r , ¢   and experiencing a local two-body interaction of strength U. In second-quantized form, the Hamiltonian is given by respectively. In the attractive regime, spin-polarized systems could exhibit several exotic superfluid phases [11] while the BCS-to-BEC transition has been addressed in the spinbalanced model [12]. Otherwise, the on-site repulsion can stand for a perfectly screened Coulomb interaction, which received considerable renewed interest in two-dimensional (2D) geometry after Anderson's proposal [13] in connection to the spectacular properties of the high-T c cuprates. However, there is still no consensus about the adequacy of the positive-U Hubbard model to capture the interplay between d-wave superconductivity, magnetism and inhomogeneous phases of copper oxides. This challenging issue is even more relevant since the latest condensed-matter experiments seem to be consistent with an intriguing scenario where spin, density and long-ranged pair correlations develop cooperatively and are spatially modulated [14,15].
The exact answer to the question of whether the 2D repulsive Hubbard model supports such intertwining of multiple orders will probably be provided only through quantum emulators like ultracold atoms. Indeed, the exact low energy properties of Hamiltonian equation (1) are only accessible in one dimension [16] and for the infinitely connected Bethe lattice through the dynamical mean-field theory [17]. In other cases, computational methods to recover exact ground-states are generally marred by exponential complexity [18,19]. Nevertheless, diagrammatic quantum Monte Carlo (QMC) simulations in continuous time have recently allowed for a determination of the phase diagram at weak coupling for small to intermediate filling [20]. Even if ultracold fermions in optical lattices already enabled us to monitor the Mott transition [21] and the development of antiferromagnetic correlations at half-filling [22], the knowledge of the phase diagram at low temperature and up to the strongly repulsive limit remains a long-term goal. In the spirit of the compelling example provided by unitary Fermi gases [23], it is highly desirable to introduce theoretical approximate schemes that could guide experiments and benefit from the progressive results of this emulation. In order to embrace the full complexity of the repulsive Hubbard model, we set up in this paper a variational approach where the ground-state is progressively reconstructed from an expansion on symmetry adapted wavefunctions without any a priori input on the relevant correlations. The key features of the method are presented in section 2. Its reliability against other numerical simulations is discussed in section 3. Finally, we proceed in section 4 to a systematic application in a four-leg ladder geometry motivated by recent experimental achievement of optical lattice plaquettes [24]. The obtained quantum phase diagram in the lattice filling-interaction strength plane highlights the intertwining of magnetic, density and pairing channels.

Methodology: the symmetry projected Hartree-Fock/Bogoliubov-de Gennes scheme
For weak coupling strength U/t the determination of correlations that spontaneously emerge from the Hubbard Hamiltonian equation (1) can be achieved by identifying the channels in which instabilities develop through self-consistent perturbative or functional renormalization group methods [25,26]. In the strongly correlated regime, the problem could ideally be tackled with Gutzwiller-type wavefunctions P )ˆˆˆp artially suppresses the double occupancy entailed in a mean-field state Fñ through the real parameter g [27]. Yet, the energy minimization has to be performed in a variational Monte Carlo framework, rendering unrestricted calculations beyond reach. Hence, the reference wavefunction must be parameterized with a limited number of relevant variables to describe specific phases, such as d-wave superfluids [28], spirals [29] or stripes [30]. A step towards unbiased Gutzwiller calculations has been recently achieved [31]. However, orders exhibiting a periodicity larger than a few lattice spacings were forbidden, in contradiction to approximate QMC results [32] revealing long wavelength modes in ground states. Alternatively, correlations beyond mean field can be generated by restoring deliberately broken symmetries through quantum number projection. In fact, the Hamiltonian equation (1) is invariant under local U(1) gauge transformations, lattice translations, spin rotations and discrete symmetries of the lattice. Thus, exact eigenstates are characterized by the number of fermions N, the total pseudo-momentum K,  the total spin S and its zcomponent S z , as well as an irreducible representation of the lattice symmetry group. All these labels will be collectively denoted by Γ in the following. Their restoration on top of a single Hartree-Fock (HF) wavefunction and before energy minimization recently yielded encouraging results for 2D clusters [33]. In particular, the exact ground state of the four-site model has been analytically recovered irrespective of the interaction strength [34]. The approach, and its analog with several Slater determinants [35,36], also proved capable to evidence interplay between spin, charge and pair degrees of freedom. Potential superfluid features would nevertheless require a very large number of Hartree-Fock (HF) basis states to be accurately captured, whereas Bogoliubov-de Gennes (BdG) ansätze are well known to be more appropriate. Hence, we focus on a more entangled trial state Y G obtained through the coherent superposition of symmetry projected HF and BdG wavefunctions: is the most general quasi-particle vacuum for a lattice with N r  sites. The Peierls-Yoccoz operator P Ĝ [37] ensures the projection on quantum numbers Γ and, according to group theory, may be expressed as a specific linear combination of (unitary) symmetry transformations T ĝ : where the coefficients λ Γ,g are proportional to the characters of the irreducible representation associated with Γ.
where  stands for the energy functional obtained with Wick's theorem. However, the normal contractions c c , now correspond to matrix elements between the non-orthogonal wavefunctions a F ( ) and ( ) which can be easily expressed in terms of quasi-particle states, occupied and unoccupied HF wavefunctions [40]. Stationarity of E Γ equation (4) with respect to the amplitudes x (HF) and x (BdG) immediately leads to a generalized eigenvalue equation: On the other hand, the energy minimization with respect to the spin-orbitals r n , f s  and Bogoliubov coefficients U , is much more involved and will be detailed in a forthcoming paper [40]. It leads to a set of self-consistent equations that reads where the matrices a b ,  G ( ) are obtained with the help of the HF/BdG mean-field Hamiltonian ij as: , , The system of equations (5)-(7) allows us to determine the optimal symmetry-projected HF/BdG wavefunction through a numerical solution in which the HF and BdG states are parameterized according to the Thouless theorem [38]. No initial assumption on the ground-state is required and the method is thus able to reveal the physics embedded in the Hubbard model equation (1) at low energy.

Reliability of the HF/BdG approach
We now address the accuracy of the wavefunction equation (2) against exact diagonalization (ED) for small clusters or QMC simulations. We focus on autocorrelation functions r r ,     ( ) ( ) and r   ( ) in the magnetic, charge and d-wave pairing channels, respectively: is the spin operator at lattice node r  (with t  the usual Pauli matrices); denotes the singlet pair-field in the d x y For a 4×4 cluster in the strong coupling regime U/t=10, 12, the HF/BdG approximation reproduces very accurately the ED data, as shown in figures 1 and 2. The determination of exact ground states for larger cells is still limited by the NP-hardness of QMC simulations, except for limited parameter spaces where the stochastic sampling is protected from the notorious sign problem. This is the case at half-filling and we present in figure 3 a [37]) and a pfaffian otherwise (see [39]).
comparison between QMC and HF/BdG spin-spin correlations for a 6×6 cluster at U/t=4. No significant difference is found, especially for the largest separation distances which are essential to indicating the development of a magnetic order. In the hole-doped regime with repulsive interactions, a sign-free stochastic sampling of the ground state is certainly possible [41,42], but it remains generally plagued by systematic errors whose origin is not totally elucidated [43]. Nevertheless, it seems that these new QMC algorithms can be accurate for closed-shell fillings and moderate interaction strengths when supplemented by quantum number projection [44,45]. Superfluid correlations in the d-wave channel have been investigated in such a framework [45] and we show in figure 4 a representative result from [31] in the intermediate coupling regime U/t=4 on a 8×8 cell. The symmetry-adapted HF/BdG wavefunction essentially yields the same pairing response r ,   ( ) as shown in figure 4. Besides, the variational energies E Γ originating from equations (5)-(7) are summarized in table 1 for the clusters and on-site interactions previously considered. The agreement is excellent for 4×4 cells with a relative error smaller than 0.5%. The quality of the approximation is quite similar for their doped counterparts, even when a negative next-nearest neighbor hopping t′ is introduced to induce frustration. As the size increases, the HF/BdG energy becomes generally less accurate and the deterioration is more pronounced if the cell is doped and/or the coupling U/t is strong. Indeed, while the discrepancy for the half-filled 6×6 cluster at U/t=4 does  = Periodic-periodic boundary conditions are imposed. ED results as well as a recent VMC calculation with a symmetry restored BCS-Gutzwiller wavefunction (mVMC) are extracted from figure 8(b) of [31]. The HF/BdG ansatz, with quantum number projection before variation, correctly reproduces the shape and magnitude of the exact pairing correlations. = Sign-free QMC calculations (GQMC) and VMC results with the BCS-Gutzwiller wavefunction (mVMC), both including symmetry restoration, are extracted from figure 10 of [31]. They are compared to the pairing correlations r   ( ) originating from the symmetry-adapted HF/BdG scheme.  between next-nearest neighbors. Exact diagonalization has been performed with ALPS [46]. QMC data are borrowed from [31,43]. The VMC calculations correspond to the original (o) [28] or improved (i) [31,49] BCS-Gutzwiller wavefunction.  not exceed 0.5%, we only recover an energy E Γ very close to the one obtained with the Gutzwiller projection on top of an optimized BCS wavefunction [28] for 64-site cells and up to U/t = 12.
For the largest clusters, the above contradictory findings regarding the accuracy of HF/BdG correlation functions and energies could be reconciled provided that improving the ansatz equation (2) only has a noticeable effect on the energy. This scenario has been validated by enlarging the variational subspace through the inclusion of several HF/BdG pairs of states via a modification of the wavefunction equation (2) according to: F ( ) denote the ith HF, BdG wavefunction in the basis, respectively. The full energy minimization would then require the simultaneous variation of all HF and BdG states. This scheme is beyond our computational facilities and we therefore limit ourselves to a sequential process. In this case, HF/BdG pairs are progressively introduced and each of them is optimized while keeping unchanged the previous basis states.
The amplitudes x (a i ) are obtained through a generalized eigenvalue problem similar to equation (5) and the set of self-consistent equations determining the structure of the current HF/BdG pair now reads As expected, such superpositions of several HF/BdG wavefunctions notably improve the energy. In the case of a doped 6×6 cluster with N=24 atoms and a coupling strength U/t=8, the symmetry-adapted HF/BdG approximation with one pair of states leads to a variational energy E Γ =−34.93 t. The use of fifty HF/BdG pairs allows us to reach an energy of −36.9 t which is comparable to approximate QMC estimates depending on the constraining state chosen to avoid the sign problem [47]. As shown in figure 5, such an improvement of the trial wavefunction equation (2) induces minimal changes in the spin, charge and d-wave pair correlation functions. A similar behavior is also reported in the next section for larger cells at various fillings.
Manifestly, the ansatz equation (9) is not the only way to progressively reconstruct the exact ground state in a subspace spanned by symmetry-projected wavefunctions. As of now, such variational strategies have only been developed with HF states, which were either stochastically sampled [48] or optimized [35,36]. The inclusion of BdG wavefunctions in the basis through equation (9) yields a notable acceleration of the convergence towards the ground-state. For instance, for N=56 interacting atoms on a 16×4 cell at U/t=12, E G is decreased to −36.018 t with ten HF/BdG pairs while a subspace twice as large is required to reach a similar energy without BdG states [36]. Efficient energy lowering may also be achieved by tuning the numbers of BdG and HF wavefunctions in the basis and the order in which they are introduced as long as the optimization is reduced to a sequential process. As an example, we consider the case of N=62 atoms loaded in a 8×8 cell for U/t=8 through the ansatz While the simple HF/BdG state equation (2) gives an energy E Γ =−34.736 t not competitive with extended BCS-Gutzwiller schemes [49], the expansion equation (11) allows us to reach a similar accuracy with E Γ =−35.961 t (for N BdG =15 and N HF =35). Figure 6 reconfirms that the physical content embedded in r ,   ( ) r   ( ) and r   ( ) is unaffected against the enlargement of the HF/BdG subspace. However, noticeable changes in the values of the order parameters extracted from the long-ranged parts of the considered correlation functions are obtained.
Finally, the present calculations tend to support the HF/BdG approximation with full symmetry restoration before variation as a reliable starting point to capture the essence of correlations entailed in the repulsive Hubbard model, at least in the magnetic, density, and superfluid channels and for moderate size clusters.

Results: quantum phase diagram of ultracold fermions loaded in optical four-leg tubes
Though the symmetry-projected HF/BdG wavefunction equation (2) displays a polynomial complexity with the number N r  lattice sites, the numerical optimization remains challenging by requiring the simultaneous determination of around N 3 r 2  parameters. Unfortunately, very large square cells are needed to support both the emergence of an off-diagonal long-ranged order linked with superfluidity and the development of the longwavelength collective modes expected in the density and magnetic channels [32]. Therefore, we now restrict ourselves to four-leg ladders which are natural steps in the dimensional crossover from the exactly solvable chain to the unknown 2D limit. This geometry can, in fact, be indirectly emulated with ultracold vapors by loading the atoms in optical tubes of plaquettes created from four wells arranged in a square pattern, as depicted in figure 7.
A three-dimensional array of such independent clusters has already been realized using an optical superlattice configuration along two orthogonal directions [24]. By tuning the laser potentials' parameters to allow for the tunneling between adjacent planes, a collection of uncoupled identical tubes could be obtained. When unfolded, each of them realizes the Hubbard Hamiltonian on a rectangular cell with four legs and periodic boundary conditions along the y-direction.
The variational state equation (2), free of symmetry breaking, is now systematically determined to unravel the relevant orders and their potential intertwining in the low-lying energy states. We focus on tubes of length L16 loaded with slightly less than one atom per site, so that they are characterized by their hole-doping δ=1−n with n the lattice filling factor. Coupling strengths ranging from the moderate U/t=4 to the strongly U/t=12 interacting regime are considered. We stress that all energy minimizations have been independently carried out, thereby allowing for crosschecking the results. Moreover, full periodic boundary conditions on finite-size clusters could bias pairing correlations by disadvantaging the d x y 2 2 channel: the corresponding wavefunction in momentum space would indeed be zero for a non-negligible fraction of wavevectors in the first Brillouin zone. Therefore, antiperiodic boundary conditions along the legs (x-direction) are chosen. Their influence is discussed in the appendix, where it is more generally shown that a tube length L=16 is large enough to ensure a weak sensitivity of the physical content to boundary conditions.
We first address a system of L=16 four-sites plaquettes with an on-site interaction U=12t and investigate the hole-doping dependence of relevant correlation functions. Each optimization involves the determination of around 10 4 complex parameters that enter the variational wavefunction equation (2)    = A rectangular 16×4 cell is considered. Spin and density autocorrelation functions are calculated from the numerical solution of the symmetry projected HF/BdG scheme. All symmetries are restored through projections on the number of atoms N , a zero total pseudo-momentum K ,  the spin-singlet subspace and the irreducible representation A 1 of the C v 2 lattice symmetry group. The latter is physically associated to a many-body wavefunction invariant under horizontal and vertical mirrors. Note that these quantum number projection are also included during the energy minimization, except for the total spin where only its z-component and parity are imposed. In both parts (a) and (b), 3D histograms in the front are obtained with one HF/BdG pair of states, while those in the back result from an enlarged subspace spanned by several sequentially optimized HF/BdG wavefunctions (five couples for N 50, 54, 58, 62 = and ten couples for N 56.  ( )in the spin-spin correlations corresponds to an antiferromagnet with a staggered magnetization oscillating in amplitude with a period of λ m =8 lattice spacings in the x-direction. Similarly, the density-density correlation function reveals inhomogeneities distributed with a period λ c =4 along the x-axis in the variational ground state. Note that these orders and their symmetry-related counterparts are necessarily superimposed to respect all invariances of the Hamiltonian. Furthermore, with λ c even, the relation λ m =2λ c characterizes stripes at the boundaries of antiferromagnetic domains separated by a π phase shift. Their intertwining with d-wave superfluidity is eventually proved by highlighting in figure 9(b) a non-zero average of the pairing correlation function r   ( ) at large separation distance r. The non-decaying tail observed for r>4 is consistent with the off-diagonal longranged order that signs superfluidity. Besides, the 4-period small oscillations of r   ( ) around its averaged value indicate the existence of pairs at a finite momentum equal to the charge-order wavevector. Such stripes with a dwave superfluidity spatially modulated in phase with the density profile have also been proposed in recent simulations [50,51] of the t-J Hamiltonian that approximates the Hubbard model in the limit U/t→∞. Superfluid domain wall states in four-leg ladders also find support from density-matrix renormalization group calculations [52] of the t-J model, despite a possible contamination by Friedel oscillations stemming from open boundary conditions [53].
Stripe-like states are robust against a decrease of the hole number as shown in figure 8 for the 16×4 cluster considered here. However, the shift of the peaks in the structure factors S m and S c reflects a doubling of the period when crossing δ=1/8. In addition, pairing correlations at large distance are totally suppressed in the 8and 4-hole systems corresponding to the perfectly filled and half-filled 8-period vertical stripes, respectively (see figure 9(c)). The spin and charge pattern associated with domain walls separated by eight lattice spacings is also realized for N=58. As for N=54, these stripes are neither filled nor half-filled and again the behavior of r   ( ) at large distance is consistent with the development of a pair-density wave of period λ c . When moving towards the half-filling limit, antiferromagnetism no longer exhibits amplitude modulation and a uniform density profile is recovered. Finally, a pure d-wave off-diagonal long-ranged order is unambiguously supported as long as such a background is doped with few holes (see figure 9(b) for N=62).
Another scenario emerges when considering an increase of the hole doping from δ≈16%. While peaks related to charge-density waves disappear, incommensurate spin-spin correlations persist. At the same time, its associated wavevector leaves the side of the Brillouin zone to its diagonal. The nature of the underlying incommensurate magnetic ordering is not unambiguously revealed by such peaks, as they are compatible with both collinear spins or spirals [29]. One way to test whether spins rotate on the lattice is to detect a non-decaying four-body correlation function between spin chirality vectors V S S S r r r u r u for different densities is shown in figure 10. The long-ranged (r>4) part systematically displays an oscillating behavior reflecting significant quantum fluctuations. Two regimes are however clearly distinguished: spiral correlations averaged over large distances vanish in the striped and antiferromagnetic states (N54), while they are non-zero and positive at larger dopings. This signal remains of small amplitude and thus rather characterizes a spiral ordering component embedded in a spindensity wave (SDW). Note that despite this, pure spiral ground states are not expected in the large-U Hubbard model considered here [54]. As shown in figure 9(a), the d-wave pairing correlation function in the SDW/spiral state displays a complex behavior at large distance, yet free of a rapid decay to zero as was found at half-filling or in the stripes at commensurate dopings. It can be viewed as the precursor of the dwave superfluidity that is better established for larger lattice fillings.
The energy minimization with the symmetry projected HF/BdG wavefunction essentially exhibits all the above features from the intermediate coupling U/t=6 to the strongly correlated regime U/t=12. The results are summarized in the quantum phase diagram shown in figure 11 for hole doping δ smaller than 1/4. The stripe-like states are stabilized in the intermediate doping range and once U/t exceeds a critical value. The latter is suppressed with decreasing δ. A similar feature has also been obtained with inhomogeneous dynamical mean-field [55] and constrained-path QMC approaches [32]. In addition, the change of the charge period from λ c =4 to λ c =8 takes place for U10t when crossing δ=1/8. Close to the half-filling limit, only antiferromagnetic correlations persist, while stripes melt for larger doping. Instead, incommensurate antiferromagnetism in the form of coexisting spiral and spin-density waves is found. It develops along the xdirection for intermediate interaction strengths (6U/t8) and tends towards the diagonal direction at large U/t. Furthermore, the spiral component appears for couplings that increase with the doping. Finally, longranged d-wave pairing correlations are systematically evidenced, except when all the holes are perfectly trapped into filled or half-filled vertical stripes. These trends are altered at smaller U/t. In particular, for U/t=4, charge inhomogeneities are missing and a clear tendency towards magnetic ordering is obtained for doping δ<16% only, in agreement with the latest diagrammatic QMC calculations [20]. Eventually, the superfluid signal is rather erratic, though this non-monotonicity proved stable against changes of boundary conditions to investigate the influence of shell effects, commonly invoked at small coupling in the attractive regime [56]. Further details are presented in the appendix.
The symmetry-projected HF/BdG phase diagram ( figure 11) in four-leg ladder geometry confirms the emergence of correlations proposed separately for the hole-doped 2D Hubbard model in the spin, charge and pairing channels. While the scenario of a competition between the resulting orders is usually retained, our findings rather point towards a subtle entanglement of the associated degrees of freedom. It induces the wide variety of strongly correlated states observed in figure 11 as a function of the hole doping. Their robustness requires persistence of the observed correlations when refining the grid of available densities by increasing the tube length L. Some representative examples are shown in figures 12 and 13 for different U/t regimes to explore additional parts of the phase diagram originally obtained at L=16. Either close to the half-filling point or on both sides of the 1/8 hole doping, no qualitatively new features appear in the spin and density autocorrelation functions. Not even the stripe period changes, when relevant. Besides, long-ranged pair correlations are still evidenced whether or not they are intertwined with antiferromagnetism or stripes. In the latter case, the increase of L allows us to grasp the oscillations of r   ( ) which clearly match the charge period λ c . It is remarkable that Figure 11. Phase diagram arising from the symmetry-projected HF/BdG approach for repulsively interacting cold fermions loaded in optical four-leg tubes of length L 16.
= Colors refer to different magnetic (charge) orders revealed by a peak at the wavevector q m  (q c  ) in the Fourier transform of the spin (density) autocorrelation function. For each hole doping d and interaction strength U t, the d-wave superfluidity symbol is made more visible when the pair correlation function exhibits off-diagonal long-ranged order.
only the quantum number projection on top of the mean-field like wavefunctions and before variation remains efficient to generate such unconventional superfluid signals for clusters with a hundred lattice sites. Indeed, we recall that the standard BdG approach fails to stabilize superfluid states for the repulsive Hubbard model. However, the value of the long-ranged tail in r   ( ) tends to decrease as compared to the length L=16 that was previously considered. While this feature could indicate the establishment of a quasi-long-ranged order in the d-wave pairing channel, the benchmark comparisons presented in section 3 point towards a deterioration of the symmetry-projected HF/BdG approximation to the ground-state when enlarging the cluster. So, one cannot exclude the need to significantly increase the dimension of the HF/BdG subspace to recover the accuracy on

Conclusion
Summarizing, we have highlighted insights into generic features of repulsively interacting ultracold fermions loaded in optical four-leg ladders through their description by the Hubbard model. First, we have shown that such systems are ideal candidates to realize a whole sequence of magnetic phases that may be tailored by varying the filling of the lattice or the ratio t/U. Above all, the long since proposed scenario of d-wave superfluidity emerging from a doped Mott insulator has been put forward thanks to energy minimizations with no physical assumption on the relevant orders. Nevertheless, such intertwining of magnetic and pair degrees of freedom manifests itself under various facets depending on whether antiferromagnetic correlations grow from homogenous collinear spins, spatially modulated spin-density waves or spirals. It also involves the charge degree of freedom as stripes that either destroy or support superfluidity, depending on their filling. These features have been extracted from symmetry-adapted states originating from quantum number projection that also induce correlations beyond mean field. Furthermore, magnetic, charge and superfluid correlations remain robust against improvements of this wavefunction. The quantum phase diagram in the four-leg tube geometry therefore provides an additional reference for the cross-validation between theory and quantum emulation from experiments that is necessary to face the exponential complexity of low-dimensional quantum matter.