Chiral superconductivity in the doped triangular-lattice Fermi-Hubbard model in two dimensions

The triangular-lattice Fermi-Hubbard

for the understanding of molecular materials of the κ-ET family [10,11,12,13,14]. Besides geometrical frustration, charge fluctuations play a role in the stabilization of quantum spin liquids [15,16]. The simulation of charge fluctuations can be accomplished by introducing highorder ring-exchange coupling to the effective spin Hamiltonian [17,18,19,20] or by considering the Hubbard model itself. Additionally, the observation of Hubbard-model physics in triangular lattices engineered in WeSe 2 /WeS 2 moiré superlattices [21,22] and in quantum simulators [23,24] has inspired scientific interest in such systems. In this work, we focus on the triangular-lattice FH model described by the Hamiltonian where α =↑, ↓ is the electron spin, ⟨ij⟩ indicates a sum over nearest-neighbour sites, c iα and c † iα respectively annihilates and creates an electron with spin α at the i-th lattice site, n iα = c † iα c iα is the number operator, t is the hopping strength, and U is the intensity of the onsite interaction. The non-interacting system is metallic and, at half-filling (⟨n iα ⟩ = 1/2), the strongly interacting FH Hamiltonian can be mapped into the antiferrognetic (AFM) Heisenberg one, whose ground state contains 120 • long-range spin order [25]. Away from these two cases (U/t = 0 and U/t ≫ 1), the precise nature of the system is still under debate. For weak interactions, numerical simulations assume an adiabatic connection with the non-interacting regime. Nonetheless, they might be overlooking a transition to a phase with a small but non-vanishing gap [1]. In fact, renormalization-group calculations predict a d + id superconductor at U/t ≪ 1 in 2D [26,27] and weak coupling analyses argue that, at weak interactions, the quasi-1D system is a Luther-Emery liquid with time-reversal symmetry breaking [28]. For intermediate interactions, the quest for spin liquids has been a subject of significant interest. Still, there is not even theoretical agreement on whether the spin liquid state exists in the FH model. While calculations ranging from variational cluster approximation [29,30,31], path integral renormalization group [32,33], strong coupling expansion [19], dual fermion approach [34] and exact diagonalization [35] to density matrix renormalization group (DMRG) [36,37,38,39] and variational Monte Carlo (VMC) [40,41] agreed in the existence of a spin liquid state, dynamical cluster approximation studies [42] and earlier VMC computations [43, 44] detected a direct transition from a metallic state to a magnetic ordered phase. Among the theories that support the existence of a spin liquid, its nature remains controversial. Infinite-DMRG calculations predict a gapped chiral spin liquid (CSL) [37,38], while VMC simulations on full 2D systems and finite-DMRG [36] support a gapless spin liquid that preserves timereversal symmetry. Another finite-DMRG study also supports the gapped CSL [39]. On the other hand, a multi-method approach finds that, at intermediate interactions, there is a competition between chiral and two distinct magnetic orders: collinear and 120 • order [45]. DMRG simulations of the extended AFM-Heisenberg model with four-spin interactions that arise naturally from Mott-insulator physics corroborate the existence of a CSL in lattice geometries closer to 2D than the ones used in DMRG simulations of the Hubbard model [20]. For the hole-doped system, the quest for unconventional superconductivity in the Hubbard model is a matter of current scientific interest due to its connection to High-T c superconductors [46,47,48,49]. A recent DMRG study of the doped triangular FH predicts a rich phase diagram with fractionalized excitations, spin and charge deconfinement and enhanced Cooper-pair correlations [50]. Another DMRG study estimates the spectral function of one single hole doped in the triangular-lattice CSL and observes spinon dynamics [51]. Also, DMRG simulations of the extended t − J model with three-spin chiral interactions in the triangular lattice predicted chiral superconductivity in the system, evidenced by quasi-long-range-order in the Cooper-pair correlations upon doping [52].
Shortly thereafter, emergent topological superconductivity was also reported in the simpler t − J model [53]. However, DMRG performed in quasi-1D lattices does not display true longrange order, and one has to rely on a slow decay of Cooper-pair correlations. Numerical simulations via the Linked-Cluster Expansion algorithm provide several benchmarks to the finite temperature triangular FH at intermediate to strong interactions [54], but a clear description of the weakly interacting regime, the classification of the spin liquid state, and whether or not a superconducting phase would appear upon doping is still elusive.
In this work, we numerically investigate the ground state of the doped triangular FH in 2D. Upon doping a non-magnetic chiral spin state (CSS) we observe true long-range order in the Cooper-pair correlations while the chiral order parameter remains finite, i.e. a chiral superconductor. To simulate the CSS, we first locate the MIT and the transition to the AFM phase.

Methods
We report the implementation of state-ofthe-art Auxiliary-Field Quantum Monte Carlo (AFQMC) to simulate the ground state of the full 2D triangular lattice FH model. By imaginarytime projection to the ground state we intend to reduce the bias from VMC calculations. A constrained-path (CP) approximation is required to restore polynomial convergence (otherwise plagued by the sign problem [55]) and the bias from the variational ansatz is not completely removed. However, simulations made in the past for square lattices away from half-filling have been shown to be accurate and provided several benchmarks [56,57,58,59,60]. See Appendix A for further details of the method and for a comparison of the CP-AFQMC estimates of the triangularlattice ground-state energy with exact diagonalization.
For the non-magnetic phases, we consider the generalized Hartree-Fock ansatz (GHF) for the imaginary-time projection. The mean-field Hamiltonian associated to the GHF state is obtained considering a partial particle-hole transformation on a BCS Hamiltonian, where M i = U eff ⟨c † i↑ c i↓ ⟩. The GHF ground state is obtained via a self-consistent diagonalization of Eq. (2) [60]. As done with Unrestricted Hartree Fock (UHF) wave functions [59], we consider U eff = min(U, U max eff ). We noticed that U max eff /t = 4 gives meaningful estimates of the non-magnetic states of the system. Our CP-AFQMC simulations with the GHF ansatz became unstable for strong interactions, see Appendix A for more explanation concerning the instability. With this ansatz, we were able to locate the MIT but not the transition to the magnetic phase. To see the transition to the 120 • AFM phase, we also perform simulations starting from a full Hartree-Fock (FHF) ansatz. The ansatz that provides the smaller energy after imaginary-time evolution is our best representative of the system ground state, see Appendix B for details.
In our simulations, we mainly consider lattices with N x = N y = 12 sites along theê x = (a, 0) andê y = (a/2, √ 3a/2) directions respectively (see FIG. 1 for the lattice vectors and a visual representation of the triangular FH model). Periodic boundary conditions are considered along the horizontal and vertical directions. We also run simulations with different lattice sizes to analyse the finite-size effects. See the Appendix C for the dependence of total energy and spin correlations on the lattice size. Finally, to address the effect of the dimensionality, we investigate quasi-1D lattices with N x = 36, N y = 3 and 4, PBC alongê y and open boundary condition alongê x . From now on, if not specified otherwise, we are considering the 12 × 12 lattice. We study spin-balanced systems at half filling (n = N/M = 1, where N is the number of electrons and M = N x N y ) and with hole doping (n < 1).

Results at half filling
We start by defining the charge structure factor with n i = n i↑ + n i↓ , around the origin k = 0 to determine whether the system is gapped or not. Since the charge gap ∆ c is proportional to lim k→0 k 2 /N (k) [61,62], a linear behavior around k = 0 indicates that the system is metallic (gapless) and a quadratic behavior indicates that the system is insulating (gapped). More precisely, for small k we have N (k) = ak 2 + bk + O(k 3 ), and whenever b ̸ = 0, ∆ c = 0. We compute the coefficients a and b by performing a quadratic fit to our data. We considered the three allowed momenta closer to the origin along k = (k x = 0, k y ), results are shown in FIG. 2 where we located the MIT around 7 < U c1 /t ≤ 8. Our critical interaction is in agreement with DMRG simulations [36,37] and a multi-method study [45]. The Linked-Cluster-Expansion calculations, after an extrapolation to zero-temperature regime, predict a critical interaction U c1 /t ∼ 7   To further corroborate our findings, we com-pute the doublon density, d = i ⟨n i↑ n i↓ ⟩/M , for which we expect different behaviors in the metallic and insulating phases [63,36]. In the metallic phase the Brinkman-Rice picture predicts that d decreases linearly with U [64], while for strong interactions the doublon density shows Heisenberg behavior [65], , where the sum runs over the nearest neighbours. In Fig. 3 we display results for d as a function of U/t. Our data for the doublon density show a deviation from the linear behavior near the MIT.  (7), see Ref. [66]. Error bars are smaller than the markers size.
Analogously, the presence of a spin gap can be accessed by the spin structure factor with S z i = n i↑ −n i↓ . We do not see the emergence of quadratic behavior in S(k) (Fig. 4), which indicates the absence of a spin gap. The excitations of the 120 • Heisenberg antiferromagnet are gapless magnons [67]. Therefore the presence of a spin gap is not a good measure to locate the AFM order in our system. On the other hand, the presence of peaks in S(k) is a signature of longrange magnetic order; for the 120 • AFM those peaks appear on the K points of the Brillouin zone. In our simulations, we see the formation of peaks on the K points which indicates the transition to the 120 • phase with critical interaction 10 < U c2 /t ≤ 11 (see Appendix D). Our estimate for the critical interaction is in agreement with recent DMRG simulations [37], which locates the transition at U/t = 10.6. We investigate the chiral order parameter where the sum is over every triangular plaquette of the lattice with vertexes i, j and k taken in the clockwise direction and S ℓ = α,β=↑,↓ c † ℓα σ αβ c ℓβ , with σ = (σ x , σ y , σ z ) a vector of Pauli matrices.
Our results for χ are show in Fig. 5 alongside the S max of S(k). We observe a competition between chiral and magnetic orders as the interaction increases as reported for quasi-1D lattices [45]. We see a sharp increase in S max in the transition to the AFM phase.
To analyse the effect of dimension in our results, we compute S(k) and χ for the quasi-1D 36 × 3 lattice. We see that the results in 2D and quasi-1D agree reasonably well, but in quasi-1D we see the AFM insulator with the peaks of the spin structure factor at the M points, which is a signature of a collinear AFM phase. We also analyse the effect of the width of the quasi-1D systems on the results by the simulation of a 36 × 4 lattice, where we see competition between 120 • and collinear orders. See Appendix D for the results of the quasi-1D simulations. The existence of a CSS was predicted in the non-magnetic insulating phase by DMRG [37,39,50]. We also detect finite chiral order before the transition to the insulating phase (U/t ≤ 7). This metallic chiral state is in disagreement with the predictions of Ref. [37], where the authors perform infinite-DMRG and extrapolate the chiral-order parameter χ on the bond dimension. A extrapolation of this sort is not needed in CP-AFQMC simulations since the summation over auxiliary-field paths restores all the correlations in the system. Moreover, other DMRG simulations of this system at weak doping also report the observation of a chiral metal [50]. Another difference, between our results and the ones of Ref. [37] after extrapolation, is the finite value of χ near the transition to the magnetic phase (11 ≤ U/t ≤ 12).Nevertheless, the computation of the chiral susceptibility in Ref.
[45] also shows evidence of chirality in the magnetic phase.

Doping the non-magnetic CSS
We search for signals of superconductivity in the doped system by looking at the long-range behavior of Cooper-pair correlations. Following the procedure described in Ref.
[50] we define the correlation function The superconductivity order parameter ∆ is given by where α =↑, ↓ andᾱ is the opposite spin of α. For pairs in the singlet state, f (↑) = −f (↓) = 1, and for triplets f (↑) = f (↓) = 1. We compute angle-averaged C(r) in the non-magnetic CSS (U/t = 9) for three distinct fillings of the hole doped system, n = 17/18, 5/6 and 11/18, covering the weak-, intermediate-, and strong-doping regimes [50]. The results are shown in Fig. 6. Quasi-longrange-order observed in previous DMRG calculations was a hint to the existence of a superconducting phase in the 2D system [50]. Our simulations provide, for the first time, clear evidence of true long-range order in the 2D triangular FH model. Furthermore, such long-range correlations increase with the hole concentration. A striking observation is that the chiral order parameter, although reduced, remains finite at the intermediate hole density (n = 5/6), χ = 0.22 (1). At the highest hole concentration, the chiral order is further suppressed. At n = 5/6, the triplet correlations at large distances are slightly stronger than the singlet ones. Upon further doping the 2D system, the singlet and triplet correlations are (almost) degenerate.
We considered the effect of reduced dimensionality by computing the Cooper-pair correlations in our 36 × 3 lattice. As expected, we do not see long-range order in quasi-1D. Still, in the quasi-1D simulations, the triplet and singlet correlations showed to be degenerate (see Appendix D). Notably, recent DMRG simulations showed a strong finite-size dependence. Indeed, singlet-pair correlations dominate for small system sizes, whereas triplet-pair correlations are enhanced for wider cylinders [50].
The source of bias in our results is twofold: finite size and the CP approximation. We analyzed the effect of finite size in the chiral order parameter χ and performed extrapolation to the thermodynamic limit in four situations where we observed chirality: i) in the metallic phase at half filling (U/t = 6); ii) in the non-magneticinsulating phase at half filling (U/t = 9); iii) in the superconducting phase at filling n = 5/6 and interaction U/t = 9; and iv) in the AFM phase at half filling (U/t = 12).   To measure the bias from the CP approximation, we computed χ for the system described by the mean-field Hamiltonian of Eq. (2). Our re-sults for U eff /t = 4 is zero, showing that chirality for U/t ≥ 4 emerges from the imaginary-time evolution regardless of the system size.

Discussion and Conclusion.
We implemented for the first time AFQMC simulations of the zero-temperature triangular-lattice FH model. The method of choice is a stateof-the-art technique that provides unbiased estimates subjected to a small systematic error arising from constraints that eliminate the fermionic sign problem. We located the transition to the insulating phase with a critical interaction between 7 < U c1 /t ≤ 8. The transition to the magnetic phase was located at 10 < U c2 /t ≤ 11. We observe a finite value of the chiral order parameter in all phases observed, and in the insulating phase there is competition between chiral and magnetic orders. Finally, we considered hole doping the non-magnetic CSS near the phase transition at U/t = 9. Our results show long-range order in the Cooper-pair correlations.
Topological states that break time-reversal symmetry and have quasi-long-range order in the Cooper-pair correlations were observed in DMRG simulations of the t − J model extended to include three-spin chiral interactions [52]. Here we reported the simulation of a chiral superconductor in the 2D triangular Hubbard model. For filling n = 5/6, we see a finite superconducting order parameter with chiral order and enhanced triplet-pair correlations. A recent DMRG study also finds emergent chiral superconductivity in the t − J model [53], which further supports our results.
A natural followup of this work would include: (i) The characterization of our CSS via spin-flux insertion to compute the Chern number [37,52]; (ii) The evaluation of nonlocal correlations to uncover string patterns [68]; (iii) The investigation of quasi-periodic lattices [69,70] where the interplay between aperiodicity and long-range order leads to exotic phases, analogously to supersolidity [71] and glass physics [72,73] in Bose systems; (iv) The study of the extended FH model. Longrange interactions can be engineered via Rydberg dressing [74,75,76], as recently realised in square lattices with unidirectional hopping [77].

Acknowledgments
We acknowledge Peter Schauss

A Constrained path AFQMC
Let a given initial state |Ψ(0)⟩ be nonorthogonal to the ground-state of the Hamiltonian H, |Ψ 0 ⟩. The imaginary-time evolution |Ψ(τ )⟩ = exp(−τ H)|Ψ(0)⟩ will asymptotically converge to |Ψ 0 ⟩ as τ (a real number) increases. In the AFQMC method the anti-symmetric wave function is written as a linear combination of Slater determinants In our simulations, ξ(Φ k ) is not considered explicitly. As the imaginary-time evolution proceeds, Slater determinants are replicated or killed. The number of a given |Φ k ⟩ in the sum resembles ξ(Φ k ) [56, 57]. In practice one starts at τ = 0 with all Slater Determinants equal to a given trial state |Φ T ⟩, an approximation to the ground state usually obtained from mean-field theories, and update them by the application of the imaginary-time evolution operator in a stochastic process. For large systems the diagonalization of the FH Hamiltonian becomes impractical. It is then usual to consider the Trotter formula [79] to factor the evolution operator into the product of three therms, where, for the Fermi Hubbard (FH) Hamiltonian, K (V ) is the hopping (interaction) therm. If we chose δτ small enough, the error introduced by neglecting the O(δτ 2 ) therms in Eq. (8) can be made smaller than the statistical uncertainty inherent to Monte Carlo calculations; the method remains numerically exact. The desired limit τ = nδτ ≫ t −1 , with t being the hopping strength, is obtained after n successive applications of this small-δτ approximation to |Ψ(0)⟩. A given iteration on the Slater Determinants is |Φ n+1 where the superscript n indicates the imaginary time τ = nδτ . The application of one-body operators in |Φ n k ⟩ results in another Slater Determinant. Therefore, both exp(−δτ K/2) only propagates |Φ n k ⟩. On the other hand, since V is a sum of two-body operators, the remaining exponential imposes a challenge. To handle it, one can use the Hubbard-Stratonovich decomposition to transform the twobody interaction into one-body ones between each electron and auxiliary fields x. We choose the spin discrete decomposition [3], with p(x) = 1/2 and γ being given by the relation cosh(γ) = exp(δτ U/2). Considering the FH Hamiltonian H, one writes where ⃗ x = (x 1 , x 2 , ..., x M ) is a configuration of auxiliary fields, with M the number of lattice sites, to be sampled within Monte Carlo calculations, p(⃗ x) = (1/2) M is a probability distribution function (pdf) and B V (⃗ x) is a product of one-body exponentials. Explicitly, The stochastic simulation of Eq. (9) is very inefficient since p(⃗ x) is constant, and an importance sampling technique is required [56,57]. The importance sampling is also useful to define an estimator to the properties of the system and to determine the constraints that eliminate the sign problem. The importance function we implement is O T (Φ n k ) = ⟨Φ T | Φ n k ⟩, which leads to the modified imaginary-time evolution with the modified pdfp(⃗ ). The pdf now depends on the overlap with the approximated ground state after and before the diffusion procedure. Since the pdfp(⃗ x) is usually not normalized, we define the normalization factor of each Slater determinant N (Φ n k ), and the iterative projection equation becomes To handle the normalization we introduce weights to each Slater Determinant |Ψ n ⟩ = k ω n k |Φ n k ⟩, and in each iteration the weights are updated as ω n k = N (Φ n k )ω n−1 k , with ω 0 k = 1 . In practice, the sampling of the pdfp(⃗ x) is done considering individually each auxiliary field in the configuration ⃗ x. A detailed description of how to implement the aforementioned sampling and how to updated the weights is given in Ref. [58].
The equivalence between |Φ n k ⟩ = − |Φ n k ⟩ causes a sign problem that prevents numerical convergence. To eliminate the sign problem, auxiliary-field paths are constrained to a region of the configuration space were O T (Φ n k ) > 0 (akin to the fixed-node approximation [80]). Due to the equivalence between the positive and negative regions of the ground state, the method is numerically exact if the nodal structure of the trial wave function is equivalent to the ground-state one. Unfortunately the latter is unknown in the majority of cases and |Φ T ⟩ is used as an approximation to it. For that reason, the constrained-path approximation has a systematic error shown to be small [56,57,58,60].
Ground-state estimates of the system total energy can be obtained using the mixed estimator, with E n k = ⟨Φ T |H|Φ n k ⟩/O T (Φ n k ) and sufficiently large n. As can be noticed, the mixed estimator is only exact if the operator in the numerator of Eq. (15) commutes with the Hamiltonian H.
Estimates of other physical observables require the back-propagation technique [56, 81,82]. The back-propagation estimator is constructed from the formula which asymptotically reaches the average value of the observable O in the ground state for τ − τ bp and τ bp ≫ t −1 . The numerical evaluation of Eq. (16) is implemented efficiently by storing the auxiliary fields sampled in the forward propagation exp(−τ bp H)|Ψ(0)⟩ and using them to back propagate ⟨Φ T |. For a detailed description of this estimator see Ref. [82]. Additional simulation details. Our time step is δτ = 10 −2 /t, with 10 3 Slater determinants. Also, we see convergence of our estimates with τ bp = 1.6/t. Instability in the strong interaction regime. We noticed that for U/t ≥ 10, the weights of the Slater Determinants suffer from strong oscillations, which cause higher variances in our final results. These oscillations are controlled by the importance function, i.e., the mean-field ansatz being projected in imaginary time.
Comparison with exact diagonalization. To have an estimate of the effect of the CP approximation in our simulations, we considered a 3 × 3 triangular lattice filled with 6 electrons under open boundary conditions. We computed the relative difference between the ground-state energy calculated with exact diagonalization and CP-AFQMC. The results are shown in FIG. 8 where we also show the relative difference associated with the mean-field GHF ansatz with U eff = min(U, 4t).

B Mean-field ansätze
The FHF Hamiltonian [83] that we test in the regime of strong interactions is where we omitted the terms that do not conserve the number of particles. We tested the FHF ansatz and the GHF one, from Eq. (2), considering U eff = min(U, U max eff ) with U max eff = 2t, 4t and 6t. The results for the ground-state energy are in Table 1.  Table 1: Ground-state energy. Simulations of the 12 × 12 system with U/t = 9 ans 12 and different initial ansätze.
Considering the variational principle, the state that provides the smaller energy after imaginary-time evolution is the best representation of the ground state of the system for a given value of U/t. For example, in FIG. 9 we show the spin structure factor of the U/t = 9 case after imaginary-time evolution starting from four different ansätze. We see sharp peaks in the spin structure factor after evolution from the FHF ansatz with U max eff = 6t. Those peaks are evidence of long-range magnetic order, but since the other ansätze provided a smaller energy, the ground state of the system for U/t = 9 is non-magnetic.
We computed the spin structure factor of the FHF ansatz with U max eff = 6t itself and concluded that long-range magnetic order after imaginary-time evolution is reminiscent of this initial choice. Therefore, for U/t ≥ 9 we perform simulations with two different ansätze: a) the GHF one with U max eff = 4t and b) the FHF one with U max eff = 6t. After imaginary-time evolution, the state with smaller energy is our best representation of the ground state of the system for a particular value of U/t. In FIG. 10, we show the difference between those energies, which indicates a phase transition to the magnetic phase with critical interaction close to U/t = 11. We emphasize that, due to the CP approximation, our imaginary-time evolution is not variational. Therefore, care must be taken when comparing the ground-state energies projected from GHF and FHF ansätze. Even though, we notice  Figure 9: Spin structure factor for U/t = 9 after imaginary-time evolution from four different ansätze. All results are for allowed momenta in the considered geometry.
that the difference between those energies is always greater than our measure of the error due to the CP approximation (FIG. 8), excluding the case of interaction U/t = 11. This corroborates the transition to the magnetic phase at strong interactions. Figure 10: Difference between energies obtained evolving the GHF (U max eff = 4t) and the FHF (U max eff = 6t) ansätze in imaginary time.

C Finite system size effects
We compute the total energy of the system at half filling for two lattice sizes, with 108 sites (N x ×N y = 9 × 12) and 144 sites (N x × N y = 12 × 12).

D Comparison between 2D and quasi-1D
In FIG. 13 we compare results for S max and χ in 2D (12 × 12) and quasi-1D (36 × 3). We see that, for the quasi-1D system, S max suddenly increases for interactions between 11 < U/t ≤ 12 indicating the transition to the AFM phase. For the 2D lattice, the increment of S max happens earlier and is much more pronounced. We emphasize that all simulations in quasi-1D were made considering the GHF ansatz. We also compute the Cooper pair correlations for the quasi-1D system at U/t = 9, and as expected by the Mermin-Wagner theorem, they do not display long-range order. We considered the fillings n = 17/18 and n = 5/6, see FIG. 14. Furthermore, In FIG. 15 we show the spin structure factor S(k) for U/t = 9, 10, 11 and 12. The maxima of S(k) is around K for the 2D lattice, but it moves to the M in quasi-1D. In both cases one can see the formation of peaks indicating spin order.
Finally, we analyze the influence of the width of the quasi-1D systems in our results. We compare the chirality and the spin correlations in 36 × 3 and 36 × 4 lattices; results are shown in FIG. 16. In both systems, we see the transition to the magnetic phase marked by a sudden increase of S max , but for the quasi-1D cylinder of width four, the peaks of the spin structure factor S(k) seem to be more intense in the K rather than in the M points as was noticed in the width-three system (refer to  FIG. 15b for the width-three result). Nonetheless, we also see high values of S(k) at the M points for the width-four system. This indicates competition between two different magnetic orders as reported in Ref. [45]. Chirality is also seen in both systems.