The Hubbard model on triangular $N$-leg cylinders: chiral and non-chiral spin liquids

The existence of a gapped chiral spin liquid has been recently suggested in the vicinity of the metal-insulator transition of the Hubbard model on the triangular lattice, by intensive density-matrix renormalization group (DMRG) simulations [A. Szasz, J. Motruk, M.P. Zaletel, and J.E. Moore, Phys. Rev. X {\bf 10}, 021042 (2020)]. Here, we report the results obtained within the variational Monte Carlo technique based upon Jastrow-Slater wave functions, implemented with backflow correlations. As in DMRG calculations, we consider $N$-leg cylinders. For $N=4$ and in the presence of a next-nearest neighbor hopping, a chiral spin liquid emerges between the metal and the insulator with magnetic quasi-long-range order. Within our approach, the chiral state is gapped and breaks the reflection symmetry. By contrast, for both $N=5$ and $N=6$, the chiral spin liquid is not the state with the lowest variational energy: in the former case a nematic spin liquid is found in the entire insulating regime, while for the less frustrated case with $N=6$ the results are very similar to the one obtained on two-dimensional clusters [L.F. Tocchio, A. Montorsi, and F. Becca, Phys. Rev. B {\bf 102}, 115150 (2020)], with an antiferromagnetic phase close to the metal-insulator transition and a nematic spin liquid in the strong-coupling regime.


I. INTRODUCTION
The quest for spin-liquid states has fascinated the condensed matter physics community since the first proposal of the resonating-valence bond (RVB) theory by Fazekas and Anderson [1,2]. This approach has been one of the first attempts to describe a Mott insulator without any sort of symmetry breaking, even at zero temperature. In recent years, spin liquids have been reported in an increasing number of materials. Examples are given by Herbertsmithite, which is well described by the Heisenberg model on the kagome lattice [3], organic compounds, like κ(ET) 2 Cu 2 (CN) 3 and Me 3 EtSb[Pd(dmit) 2 ] 2 [4,5], or transition-metal dichalcogenides, like 1T-TaS 2 , whose low-temperature behavior could be captured by the Hubbard model on the triangular lattice [6].
An important open question concerns the nature of the insulating phase of the two-dimensional Hubbard model on the triangular lattice at half filling. Most of the investigations have been concentrated to its strong-coupling regime, where only spin S = 1/2 degrees of freedom are left. Here, spin liquids can be systematically classified, according to the projective-symmetry group (PSG) theory [7][8][9]. In particular, one can distinguish between Z 2 and U (1) spin liquids, according to the low-energy symmetry of the emerging gauge fields [10]. Starting from the Heisenberg model with nearest-neighbor (NN) super-exchange J, spin-liquid phases are expected to be stabilized when including either a next-nearest-neighbor (NNN) coupling J or a four-spin ring-exchange term K. The latter one can be justified within the fourthorder strong-coupling expansion in t/U and is usually considered for an effective description of density fluctuations close to the Mott transition [11]. As far as the J − J model is concerned, a gapless U (1) spin liquid has been proposed by both variational Monte Carlo (VMC) [12,13] and recent density-matrix renormalization group (DMRG) calculations [14], while older DMRG results suggested the presence of a gapped spin liquid [15,16]. In addition, also ring-exchange terms may stabilize a gapless spin liquid (with a spinon Fermi surface), as proposed by earlier VMC studies [17] and confirmed by later DMRG simulations [18,19]. Further VMC investigations suggested two other gapless spinliquid states, none of them possessing a spinon Fermi surface [20]. However, more recent tensor-network approaches, implemented from Gutzwiller-projected wave functions, do not support the existence of a gapless spin liquid [21].
Recently, chiral spin liquids attracted much attention, because of their similarities with quantum Hall states [22,23]. Interestingly, chiral states may exist not only when the Hamiltonian explicitly breaks timereversal symmetry (as in the quantum Hall effect) [24], but also as a result of a spontaneous symmetry breaking phenomenon [25]. On the triangular lattice, some evidence of this exotic phase has been obtained by adding a scalar chiral interaction to the Heisenberg Hamiltonian [26,27], or even in a fully symmetric Heisenberg model with super-exchange couplings up to the third neighbors [28]. A PSG classification of chiral states is possible, as worked out for the fermionic case for different lattices [29]. In particular, two simple Ansätze can be constructed [30]: The first one (dubbed CSL1) is a U (1) chiral spin liquid, with complex hoppings defined on a 2 × 1 unit cell and no pairing; the second one (dubbed CSL2) is a Gutzwiller projected d + id superconductor.
The situation in the Hubbard model, characterized by a kinetic term t and an on-site Coulomb repulsion U , is much less clear. The main difficulty comes from the presence of density fluctuations, whose energy scale is related to U that is much larger than the typical energy scale of spin fluctuations, i.e., J = 4t 2 /U . Therefore, it is not simple to detect tiny effects related to spin degrees of freedom when density fluctuations are present. In addition, numerical methods like exact diagonalization or DMRG suffer from the fact that the local Hilbert space is doubled with respect to the case of S = 1/2. Nevertheless, this effort is necessary in order to capture density fluctuations that are inevitably present in real materials. The possibility that a spin-liquid phase may exist not only in the strong-coupling regime U t, but also close to the metal-insulator transition has been discussed by different theoretical and numerical approaches in the past [31][32][33][34][35][36][37]. The term weak-Mott insulator has been used in this case, namely when a spin liquid intrudes between the weak-coupling metal and the strong-coupling antiferromagnetic insulator [11]. In particular, recent extensive DMRG calculations [38,39] highlighted the possibility for a gapped chiral spin liquid close to the Mott transition. A possible description of such a state has been proposed within a bosonic RVB description [40], as well as within a spin model with the four-spin ring-exchange term [41]. Calculations of Ref. [39] are limited to 4-leg cylinders, which highly frustrate the 120 • magnetic pattern, since the corresponding k vectors are not allowed by the quantization of momenta. Instead, in Ref. [38] also 6-leg cylinders have been considered, even if the presence of two almost degenerate momentum sectors at intermediate U/t can make the interpretation of the results not completely trivial. In addition, a recent study at finite temperature, still focusing on 4-leg cylinders, highlighted the concomitant presence of chiral correlations and nematic order at finite, but low, temperature and intermediate coupling [42]. Instead, a DMRG investigation on 3-leg cylinders suggested that a gapless and nonchiral spin liquid appears close to the Mott transition [43]. A more conventional picture, with a direct transition between a metal and an insulator with magnetic order, has been found in Refs. [44][45][46]. We would like to remark that the analysis of the insulating phase in the vicinity of the Mott transition is complicated by the significant difference in locating the Mott transition, as observed with the different methods. Finally, the effect of NNN hopping has been addressed in Ref. [47], using the VCA method with few (12) sites, leading to a large spin-liquid region for t /t > 0, and in a VMC study that suggests always a direct transition between a metal and a magnetic insulator, even if with an asymmetry between positive and negative values of t /t [48].
In this work, we present variational Monte Carlo results, based upon Jastrow-Slater wave functions and backflow correlations, for the Hubbard model on the triangular lattice on N -leg cylinders, with N = 4, 5, and 6. The role of the NNN hopping term t is also discussed. On 4-leg cylinders, we find a chiral spin liquid on a relatively extended region in the vicinity of the metalinsulator transition, when t /t < 0. This intermediate phase is gapped, in analogy with DMRG results. In addition, a gapless nonchiral state appears at larger values of the Coulomb repulsion. For t /t ≥ 0, the Mott insulator is always nonchiral, with quasi-long-range 120 • magnetic order, even if the chiral spin-liquid state is quite close in energy, at least in the vicinity of the Mott transition.
On 5-and 6-leg cylinders, the chiral spin liquid does not give the best variational energy. For N = 5, the whole insulating regime is described by a gapless nonchiral spin liquid, while for N = 6 the phase diagram is similar to the one found in the two-dimensional case [48], with antiferromagnetic order close to the metal-insulator transition and a gapless nonchiral spin-liquid phase at strong coupling. The paper is organized as follows: in section II, we describe the model and the various variational wave functions, as well as the quantities that have been used to obtain the important information; in section III, we present the numerical results; finally, in section IV, we draw our conclusions.

II. MODEL AND METHOD
We consider the single-band Hubbard model on the triangular lattice: where c † i,σ (c i,σ ) creates (destroys) an electron with spin σ on site i and n i,σ = c † i,σ c i,σ is the electronic density per spin σ on site i. The NN and NNN hoppings are denoted as t and t , respectively; U is the on-site Coulomb interaction. We define three vectors connecting NN sites, a 1 = (1, 0), a 2 = (1/2, √ 3/2), and a 3 = (−1/2, √ 3/2); in addition, we also define three vectors for NNN sites, b 1 = a 1 + a 2 , b 2 = a 2 + a 3 , and b 3 = a 3 − a 1 . We consider clusters with periodic boundary conditions defined by T 1 = L 1 a 1 and T 2 = L 2 a 2 , in order to have L = L 1 × L 2 sites. We focus on cylinders with four (L 2 = 4), five (L 2 = 5), and six (L 2 = 6) legs, see the case with L 2 = 4 in Fig. 1. Most of the calculations have been done with L 1 = 30, which is large enough not to suffer from significant finite-size effects. The half-filled case, where the Mott transition takes place, is considered. In this case, only the sign of the ratio t /t is relevant and not the individual signs of t and t .
Our numerical results are obtained by means of the VMC method, which is based on the definition of suitable wave functions to approximate the ground-state properties beyond perturbative approaches [49]. In particular, we consider the so-called Jastrow-Slater wave functions that include long-range electron-electron correlations via the Jastrow factor [50,51], on top of an uncorrelated Slater determinant (possibly including electron pairing). In addition, the so-called backflow correlations will be applied to the Slater determinant, in order to sizably improve the quality of the variational state [52,53]. Thanks to Jastrow and backflow terms, these wave functions can reach a very high degree of accuracy in Hubbard-like models, for different regimes of parameters, including frustrated cases [54]. Therefore, they represent a valid tool to investigate strongly-correlated systems, competing with state-of-the-art numerical methods, as DMRG or tensor networks.
Our variational wave function for describing the spinliquid phase is defined as: where J d is the density-density Jastrow factor and |Φ MF is a state where the orbitals of an auxiliary Hamiltonian are redefined on the basis of the many-body electronic configuration, incorporating virtual hopping processes, via the backflow correlations [52,53]. The density-density Jastrow factor is given by where n i = σ n i,σ is the electron density on site i and v i,j are pseudopotentials that are optimized for every independent distance |R i − R j |. The density-density Jastrow factor allows us to describe a nonmagnetic Mott insulator for a sufficiently singular Jastrow factor v q ∼ 1/q 2 (v q being the Fourier transform of v i,j ) [50,51]. The auxiliary Hamiltonian is then defined as follows: where ξ k =˜ k − µ defines the free-band dispersion (including the chemical potential µ) and ∆ k is the singlet pairing amplitude. In our previous work on the twodimensional lattice [48], we found that the best spin liquid has a nematic character, the hopping terms being given by: Instead, the pairing amplitudes are: which possess a d-wave symmetry on the two bonds with hopping t. In two dimensions, we foundt d ≈ 0 and ∆ d ≈ 0, while on cylinders they may assume finite values. Remarkably, this choice (with different couplings along a 2 and a 3 ) gives the best variational energy also on cylinders, this implying an explicit breaking in pointgroup symmetries.
In addition, we focus on chiral spin-liquid states, which have been claimed to be relevant both in the Heisenberg limit [30] and in the Hubbard model close to the Mott transition [38,39]. The CSL2 state is a projected d + id superconductor characterized by uniform (real) hopping along NN and NNN bonds and a pairing where ω = e 2iπ/3 . Another chiral state (dubbed here CSL3) may be defined by the hopping amplitude of Eq. (5) and a different d + id pairing structure: (8) Finally, a chiral spin liquid with U (1) symmetry has been proposed (dubbed CSL1 in Ref. [30]), with magnetic fluxes piercing the elementary plaquettes. In presence of density fluctuations, this state breaks the translational symmetry and does not give a competitive variational energy. Therefore, in the following, this Ansatz is not reported.
Within the two-dimensional case, antiferromagnetically ordered wave functions represent an important class of states, since a large portion of the phase diagram corresponds to phases that spontaneously break the SU(2) spin symmetry. Cylinders are quasi-one dimensional systems, in which a continuous symmetry cannot be broken. Nevertheless, variational wave function can be still constructed from a magnetically ordered Slater determinant. Then, density and spin correlations may be inserted by Jastrow factors: here, J d is the density-density term of Eq. (3) and J s is the spin-spin Jastrow factor, which is written in terms of a pseudopotential u i,j that couples the z-component of the spin operators on different sites Finally, |Φ AF is obtained, after taking into account the backflow corrections, from the following auxiliary Hamiltonian: where, k is the free dispersion of Eq. (1), S i = (S x i , S y i , S z i ) is the spin operator at site i and M i is defined as M i = [cos(Q · R i ), sin(Q · R i ), 0], where Q is the pitch vector. The three-sublattice 120 • order has Q = ( 4π 3 , 0) or ( 2π 3 , 2π √ 3 ), while the stripe collinear order with a two-sublattice periodicity has Q = (0, 2π √ 3 ) or Q = (π, π √ 3 ). On 6-leg cylinders, the pitch vector corresponding to the 120 • order is allowed by the quantization of momenta; instead, on 4-and 5-leg cases it is not allowed and we take the closest possible momentum. On 5 legs, also the pitch vector of the stripe collinear order is not allowed.
In general, the effect of the spin-spin Jastrow factor J s is to reduce the value of magnetic order of the uncorrelated Slater determinant [55,56]. In purely onedimensional systems, the presence of a long-range Jastrow factor is able to completely destroy magnetic order, leading to the correct behavior of the spin-spin correlations [57]. On cylinders with a finite number of legs N , a residual magnetic order persists, thus giving rise to a spurious wave function that breaks the SU(2) symmetry. Here, we interpret the possibility to stabilize this kind of variational state as the tendency to develop magnetic order in the two-dimensional system. For simplicity, in the following, the Ansatz of Eq. (9) will be denoted by "antiferromagnetic". We remark that, in principle, it would be possible to restore the SU(2) symmetry by projecting on the S = 0 subspace [58]. However, this procedure is rather computationally expensive, whenever the computational basis has a definite value of S z = i S z i but not of S 2 = ( i S i ) 2 .
All the pseudopotentials in the Jastrow factors, the parameters in the auxiliary Hamiltonian, as well as the backflow corrections are optimized with the stochastic reconfiguration method [49].
In order to assess the metallic or insulating nature of the ground state we can compute the static densitydensity structure factor: where . . . indicates the expectation value over the variational wave function. Indeed, density excitations are gapless when N (q) ∝ |q| for |q| → 0, while a gap is present whenever N (q) ∝ |q| 2 for |q| → 0 [53,59]. Analogously, the presence of a spin gap can be checked by looking at the small-q behavior of the static spin-spin correlations [60]:

III. RESULTS
Here, we discuss the results for the variational energy of different states on the 4-leg cylinder geometry. Let us start from the case with t /t = 0, see Fig. 2 (upper panel). In this case the Mott transition occurs between U/t = 9 and 9.5, as extracted from the low-q behavior of the density-density correlations, see that the conducting phase is a standard metal, with neither magnetic nor superconducting order. Instead, in the insulating phase, the optimal wave function is the antiferromagnetic one with the pitch vector corresponding to approximately 120 • order, i.e., Q = ( 2π 3 , 7π . The overall situation is not much different from what has been obtained, within the same approach, in the two-dimensional limit [48] (except the fact that in the latter case, a true antiferromagnetic order settles down). We also remark that the energy gain of the antiferromagnetic state with respect to the spin-liquid one is smaller on four legs than in two dimensions.
Then, a large spin-liquid region appears immediately above the Mott transition, by including a finite NNN hopping t /t = −0.3, see Fig. 2 (middle panel). Here, the metal-insulator transition takes place at U/t = 11.5±0.5, see Fig 3 (middle panel). The best variational state, between U/t = 12 and 16, is given by the CSL3, even though the other spin-liquid states are very close in energy. By increasing the ratio U/t, the ground state passes through an intermediate phase where the best variational state is the antiferromagnetic one (with collinear order), before entering a further (strong-coupling) spin-liquid region that has no chiral features, in analogy with the results previously obtained in two dimensions [48]. In Fig. 2, we do not report the flux phase CSL1, since its variational energy is always significantly higher than the other states. 3 ), up to U/t ≈ 20. For larger values of the electron-electron repulsion, a nonchiral spin-liquid state emerges. Note that the energies reported for the antiferromagnetic state with collinear order below the Mott transition correspond to a local minimum with insulating features.
In order to determine the nature of the chiral spinliquid state, we analyze the spin-spin correlations by computing the spin-spin structure factor of Eq. (13). In Fig. 4, we report calculations for t /t = −0.3 and values of U/t across the Mott transition. The main result is that the chiral spin liquid, realized close to the Mott transition, has a spin gap, since S(q) ∝ |q| 2 for small values of the momentum q. This is in agreement with recent DMRG studies [38,39]. We remark that this feature is solid, since it is also shared by the other two spinliquid states with nearby energies, i.e., the CSL2 and the nonchiral one parametrized by Eqs. (5) and (6). On the contrary, the large-U state is gapless. In this regime the optimal parameterst d ≈ 0 and ∆ d ≈ 0 lead to a gapless spectrum in the auxiliary Hamiltonian (4), thus indicating that the nature of the unprojected state is not changed when including the Jastrow factor. The optimal chiral spin liquid (close to the Mott transition), as well as the nonchiral one (in the strongcoupling regime) are very anisotropic, as shown by computing the nearest-neighbor spin-spin correlations D j = S z Ri S z Ri+aj , with j = 1, 2, 3. For example, for t /t = −0.3 and U/t = 12, the CSL3 state has D 1 = D 3 = −0.029(1) and D 2 = −0.069 (1). For U/t = 20, the nonchiral spin liquid has D 1 = D 3 = −0.101(1) and D 2 = +0.041 (1). Within the error bar, these results are the same from L = 18 × 4 to L = 30 × 4. As discussed in section II, this anisotropy follows directly the parametrization of the spin-liquid state, see Eqs. (5) and (6) for the nonchiral Ansatz and Eqs. (5) and (8) for the CSL3.
Then, we show the stability of the chiral spin liquid when going from N = 4 to N = 6. Results are shown in Fig. 5, together with the ones for a truly two-dimensional cluster (with L = 18 × 18 sites), which has been already discussed in our previous work [48]. On a twodimensional cluster, the CSL3 state is a local minimum, with energy higher than the other states. Instead, on 6leg cylinders, the CSL3 state is not reported, since, upon optimization, it converges to the non-chiral state. The most important fact is that no chiral phases are present in the insulating region close to the metal-insulator transition (which, for N = 6 appears between U/t = 11 and U/t = 12). Here, the insulating phase is either an antiferromagnet with collinear order, in the vicinity of the Mott transition, or a gapless nonchiral spin liquid, in the strong-coupling regime. Note that, also in this case, the antiferromagnetic state with collinear order becomes a local minimum below the Mott transition. The reason for the stabilization of the chiral state on 4-leg cylinders comes from its remarkable energy gain when going from N = 6 (or equivalently two dimensions) to N = 4; by contrast, the variational energies of the antiferromagnetic state do not change much when varying N . Overall, the resulting phase diagram for N = 6 is qualitatively similar to the one obtained in two dimensions. Therefore, within our approach, the chiral spin liquid exists only for particular values of N , like on 4-leg cylinders.
Finally, we have also considered cylinders with an odd number of legs, i.e., with L 2 = 5. This is a particularly frustrated case, since both 120 • and stripe collinear magnetic correlations are not allowed by the quantization of transverse momenta. Results for the energies of the different variational states are reported in Fig. 6. The Mott transition is determined also in this case by looking at the static structure factor of Eq (12). In this case, the insulator is always a gapless spin liquid, with no chiral features. Indeed, the best variational state is the one defined by Eqs. (5) and (6), with optimal variational parameters ∆ d ≈ 0 andt d ≈ 0, which is the same as the large-U spin liquid reported on the 4-leg case. The two magnetic states are now both disfavored because of the 5-leg geometry and they are approximated by the pitch vectors Q = (π, 3π 5 √ 3 ) (for the stripe collinear order) and Q = ( 2π 3 , 26π 15 √ 3 ) (for the 120 • order). The two chiral states (CSL2 and CSL3) have also higher energies with respect to the nonchiral one. Our finding is in agreement with what reported by DMRG in Ref. [38], where no chiral features are observed on the 5-leg cylinder when using periodic boundary conditions.

IV. CONCLUSIONS
In summary, we have studied the Hubbard model on cylinders with a triangular lattice geometry by means of the VMC approach. Both a NN hopping t and a NNN hopping t are considered in the model. First, we fo-cused on the 4-leg case, with different values of the ratio t /t. For t /t < 0, a spin liquid is stabilized in the vicinity of the Mott transition. This state is a gapped chiral spin liquid that also breaks the point-group symmetry. At larger values of U/t, a further gapless spin liquid appears. For t = 0, the insulating region is always antiferromagnetic (with approximately 120 • order), while, for t /t > 0, we observe a gapless spin liquid in the strong-coupling regime. However, the chiral spin liquid disappears on cylinders with 5 and 6 legs, as well as in the truly two-dimensional case. In these cases, a gapless spin liquid survives in the large-U region. These results are summarized in Fig. 7.
Our calculations convey two main messages: On one side, the spin liquid that we obtain on the 4-leg cylinder, close to the Mott transition, is chiral and spin gapped, in agreement with recent DMRG calculations [38,39]. In addition, the best chiral state breaks the reflection symmetry, as also suggested in the finite-temperature tensornetwork method of Ref. [42]. Nevertheless, within variational Monte Carlo, an additional NNN hopping is necessary to stabilize the chiral state. On the other side, our results suggest that a chiral spin liquid exists only in particular geometries (e.g., the 4-leg cylinder). Instead, on cylinders with 5 and 6 legs (as well as in two dimensions), the chiral spin liquid either is not stable upon optimization or has a variational energy that is quite higher than the optimal state. Finally, we observe that chiral flux phases (defined on the 2 × 1 unit cell) have a variational energy that is not competitive with other wave functions. As every variational calculations, our results suffer from an intrinsic bias, given by the choice of the variational Ansatz; still, the Jastrow-Slater state possess a large flexibility, being able to describe a wide variety of different phases, including quantum spin liquids, with or without chiral order. The fact that we do not observe a chiral spin liquid in 5-and 6-leg cylinders and in two dimensional clusters, suggests that either this state is not present or it cannot be represented by the Ansätze that have been considered here.