Metal-insulator transition and quantum magnetism in the SU(3) Fermi-Hubbard Model: Disentangling Nesting and the Mott Transition

,

However, in this geometry and at half-filling the properties of the SU(2) Hubbard Hamiltonian are determined simultaneously by the special features of the noninteracting dispersion -perfect nesting at k = (π, π)and by the on-site repulsion U .This leads to the anomalous feature that the critical interaction strength for the metal-insulator transition is U c = 0, and a blurring of 'Slater' insulating behavior associated with the opening of an antiferromagnetic gap at (π, π), and the Mott insulator which does not rely on long range magnetic symmetry breaking but only on the high energy cost for local double occupancy.While one can separate these effects, for example by adding a next-near-neighbor hopping t ′ , it is hard to tune such parameters in experimental realizations.
Simultaneously, the rapid development of experiments with ultracold alkaline-earth-like atoms (AEAs) that exhibit a natural SU(N ) symmetric interaction [47][48][49][50][51][52][53][54][55][56][57] provides a concrete and exciting setting to explore these intriguing phases, including the explicit separation of the non-generic features of the band structure on a square lattice from the effect of interactions.Ongoing experiments aiming to pair quantum gas microscopes [58][59][60][61] with AEAs are expected to soon enable real-space imaging of correlations, which demands a deeper understanding of the order in the SU(N ) Hubbard model.
It is a major challenge to obtain reliable answers to the question we address here, namely the ground-state properties of the SU(3) FHM on a 2D square lattice at 1/3-filling.A high accuracy treatment is required in a strongly correlated system, as well as large supercell sizes to extract the thermodynamic limit.Unlike its SU(2) counterpart, a fermion sign problem [62] is present in the SU(3) model at ⟨n⟩ = 1.In this paper we use a novel implementation of the constrained path (CP) auxiliaryfield quantum Monte Carlo (AFQMC) methodology [63] which iteratively refines the constraining Slater determinant [64].This allows us to quantify the location of the metal-to-insulator transition and determine a phase diagram.Above a critical interaction strength U c which is roughly half the bandwidth, the system develops AFM order.We also discover two intermediate AFM orders (which we label 3-2 and 3-4 AFM) along the route to the large U Heisenberg limit.

Hamiltonian and Methodology. The SU(3) FHM is defined by the Hamiltonian
where c † iσ (c iσ ) is the creation (annihilation) operator for a fermion with spin flavor σ = A, B, C on site i on a 2D square lattice.N s = L x × L y denotes the number of lattice sites, n iσ = c † iσ c iσ is the number operator for flavor σ on site i, ⟨i, j⟩ denotes nearest-neighbor pairs, t is the hopping amplitude which will serve as the energy unit, and U is the interaction strength.We study the ground state of the Hamiltonian in Eq. ( 1) as a function of U in the spin-balanced case at 1/3-filling (where there is an average of one fermion per lattice site), that is, in the sector of Hilbert space with To characterize the ground-state properties, we measure energies, correlation functions, and structure factors.For example the spin-resolved two-body correlation functions, . This is complemented by the corresponding momentumspace measures, for example the equal time, equal species, density structure factor S(k , where r i and r j are the coordinates of sites i and j. We perform unrestricted Hartree-Fock (UHF) calculations as preliminary explorations to suggest candidate phases of matter, and to generate trial wavefunctions for AFQMC (see below).In the UHF treatment we approximate the Hamiltonian as H HF = σ H σ HF , where H σ HF = −t ⟨i,j⟩ c † iσ c jσ + h.c.+ U eff i,τ ̸ =σ ⟨n iτ ⟩n iσ − 1 2 ⟨n iτ ⟩⟨n iσ ⟩ , and determine the mean-fields ⟨n iσ ⟩ self-consistently from values initialized randomly or constructed from several possible ordered patterns proposed for the 1/3 filling SU(3) FHM [15,16].Our method reproduces the SU(2) UHF phase diagram [6,65].In H σ HF above we have allowed an effective interaction strength U eff , which does not have to equal the "physical" U .This turns out to be important in allowing a self-consistent procedure between UHF and AFQMC to determine a best U eff value for producing the trial wave function for the CP constraint [64].As discussed in more detail later, the four phases depicted in Fig. 1 emerge as mean-field ground states for different values of the interaction strength; however, estimates of the critical values of U or even the sequence of phases could differ significantly from the exact result.
In order to compute the ground-state properties beyond the mean field level, we use the state-of-the-art CP-AFQMC method [63,66].AFQMC is a projection quantum Monte Carlo approach, based on the fact that the ground state can be obtained by acting an imaginary time evolution operator on a trial wave function, |ϕ 0 ⟩ ∝ lim τ →∞ e −τ H |ϕ T ⟩, as long as the trial wave function is non-orthogonal to the ground state ⟨ϕ T |ϕ 0 ⟩ ̸ = 0. Within AFQMC, a combination of the Trotter decomposition and the Hubbard-Stratonovich transformation maps the imaginary time evolution onto a random walk in the manifold of Slater determinants, and ground state properties are averages over this random walk.Observables are measured by back-propagation [63,67] in the open-ended random walk approach of CP-AFQMC.
AFQMC yields exact ground state correlations in special situations, like the SU(2N ) Hubbard model at halffilling [36,41,43], and electron phonon Hamiltonians including the Holstein [68][69][70][71] and the Su-Schrieffer-Heeger models [72][73][74][75].In the SU(3) model, the infamous fermion sign problem [62,[76][77][78] necessitates a constrained path approximation which imposes a sign or gauge condition on the sampling in auxiliary-field space [63,79] relying on a trial wave function.In this work we leverage the recent methodological advances [64], to implement self-consistency loops which minimize the biases related to the trial wave function.The very promising success of such techniques in the SU(2) model, in both the repulsive [80][81][82][83] and attractive regimes [84], together with several cross-checks we did in this work (Supplemental Material, Fig. 5, and analysis below), demonstrates the high accuracy of the method.
Results.Our main results are summarized in Fig. 1, as a ground-state phase diagram of the 1/3-filling (⟨n⟩ = 1) SU(3) FHM.Panels (a)-(d) illustrate the patterns in the magnetic correlations in each of the phases.Black and white represent high and low density correlations for a single flavor, respectively.The many-body phase diagram is given in panel (f).(The UHF result is shown in panel (e) for reference and comparison.)At small U , the system is in a paramagnetic (PM) (metal) phase, while at large U , multiple occupancies are suppressed and superexchange causes adjacent sites to favor different spin species, and therefore generates an anti-ferromagnetic (AFM) pattern.In contrast to the SU(2) case where the AFM correlation determines a unique (up to translations) spatial pattern, many patterns are possible for SU (3).We find multiple ordering vectors occur, representing distinct ground states for different U values.∼ 3.25 than in UHF, reflecting the fact that the latter overestimates the ordering tendency due to its absence of fluctuations.As in UHF, AFQMC also gives the 3-3 AFM as the ground state at large U , consistent with the ground state found in prior work in the Heisenberg limit [25].However, the energy of this state, obtained from AFQMC, is nearly degenerate with the 3-4 pattern for the larger U values investigated here.
To further characterize the phase in different U regions, the ground state structure factor at the three characteris- tic momenta k = (2π/3, π/2), (2π/3, π), (2π/3, 4π/3) are shown in Fig. 3 as functions of U .The vertical lines separate the regimes where the different AFM orders are stable.There is no peak in the PM metal; in the 3-2 phase, a peak is present at k = (2π/3, π); in the 3-4 phase, peaks are seen at both k = (2π/3, π) and (2π/3, π/2); finally the 3-3 AFM phase shows a peak at k = (2π/3, 4π/3).In a UHF calculation, the phase boundary can be deter- mined by either searching for the global ground state at each U , or by comparing the energies of different selfconsistent ordered solutions -ensuring in either cases that the thermodynamic limit is reached.This is challenging to execute in many-body computations where statistical error bars and (much larger) finite-size effects can exceed energy differences.We use a combination of strategies to resolve the different boundaries as illustrated in Fig. 4.
The transition between the PM metal and the 3-2 AFM phases can be observed with robust and unambiguous self-consistent AFQMC calculations.In Fig. 4(a), we show the computed structure factor S(2π/3, π) as a function of U for a number of lattice sizes.At each U , we start the calculation from a trial wave function generated from UHF with an essentially arbitrary U eff , and perform self-consistency via the natural orbitals [64].The selfconsistent iteration consistently converges to the same final value (see Fig. 5 in SM) which is shown in the plot.The existence of the quantum critical point U c1 is evident from the different size dependencies of S(2π/3, π) as U is varied across it.Extrapolating the finite-temperature compressibilities leads to a U c1 [85] consistent with our results.A sharp contrast is seen between this behavior and that at 1/2-filling, shown in the inset of panel (a), emphasizing the result of decoupling of nesting from the Mott transition: At ⟨n⟩ = 3/2, the spins form an alternating pattern "A(B&C)A(B&C)..." which singles out one species (e.g.A) with the other two (B and C) behaving interchangeably [13].This spin arrangement has period 2 in both the x and y directions, hence giving rise to a structure factor peak at (π, π).In this case, the dependence of the structure factor on lattice size starts as soon as the interaction U is turned on, implying that the long-range order develops at U c1 = 0, similar to the half filling case in the SU(2) Hubbard model.This is expected since the SU(3) case also has the same perfectly nested Fermi surface at = (π, π), as argued in [13].
To determine the nature of the ordered AFM phases, we apply complementary measurements, with additional non-self-consistent calculations using multi-determinant trial wave functions.As discussed above, the UHF solutions using different U eff can produce distinct patterns, which can be used as trial wave function and the initial Slater determinant to start our AFQMC calculations at any U .However, this turns the whole procedure from a linear process (projection) into a nonlinear process, which can get stuck in local minima.In this situation, we determine the true ground state by comparing the energies obtained from different constraints, as illustrated in panel (b) of Fig. 4. Three regions are seen, with the lowest energies given by the 3-2, 3-4, and 3-3 states, respectively.Their boundaries, where two energy curves cross, give two more transition points, U c2 and U c3 .AFQMC with single Slater determinant constraints typically achieves energy accuracy of well over sub-percent level [79][80][81]83].When the energy difference becomes much smaller, we augment our process -in the vicinity of U c2 we use a superposition of two Slater determinants corresponding to 3-2 & 3-4 AFM as the trial wave function.This leads to a finite structure factor S(2π/3, π) (3-2 AFM order) on one side, in contrast with finite S(2π/3, π) and S(2π/3, π/2) (3-4) on the other, as illustrated in Fig. 4(b2).Near U ∼ U c3 , we conduct the same procedures, as shown in Fig. 4(b3).Beyond U/t ∼ 8.25, however, even this procedure leads to ambiguous results, where both 3-4 and 3-3 structure factors can become finite, depending on the random number seed.This combined with the tiny energy differences indicate that, to within the resolution of the current calculations, the 3-4 and 3-3 orders are essentially degenerate for the larger U values investigated here.
Discussion.Using high-accuracy many-body numerical computations, we have mapped out the ground-state phase diagram of the two-dimensional SU(3) Hubbard model at 1/3 filling (one particle per site).We find a metal-insulator phase transition at the quantum critical point U c1 ∼ 5.5, above which long range magnetic orders develop.Compared to the more ubiquitous SU(2) counterpart, the SU(3) case allows the separation of strong interaction effects at integer filling from non-generic features of band structure -nesting of the Fermi surface and a van Hove singularity in the density of states.An immediate consequence is the non-zero value of U c .
Above U c1 , our results suggest two novel intermediate magnetic orders (3-2 and 3-4 AFM phases) before the final, strong coupling 3-3 phase.We determined the critical interaction strengths for the ordering wave vectors by comparing their energies and also by using a superposition of two states as the initial and trial wave functions in our AFQMC calculations.The results given by these two approaches are consistent.Beyond U c3 ∼ 8.25, the 3-3 AFM energy appears to be slightly lower than 3-4 AFM, but they are almost degenerate at the interaction strengths considered in the study.
The metal-insulator transitions and novel magnetic ordered phases predicted in our work provide important steps to understand SU(N ) physics.This is especially timely with the intense experimental efforts to emulate such systems using ultracold AEAs.How these magnetic orders evolve at finite temperature and other potentially exotic phases in the doped systems emerge are interesting questions for future work, especially in the experimental context.

Supplemental Materials
In these supplemental materials, we present additional details concerning: (1) Self-consistency procedures we conduct in the constrained path Auxiliary Field Quantum Monte Carlo (CP-AFQMC) method; (2) Comparisons of energies in different phases and finite lattice size effect analysis in the Hartree-Fock.FIG. 5. Illustration of convergence of real space densitydensity correlation.We have chosen Uc1 < U = 7.5 < Uc2, where the system is in the 3-2 AF phase.Spin correlations converge from lower and higher values when the starting point is in the 3-2 phase using an initial Slater determinant with U eff = 3.5 and 10 respectively.The legend indicates the number of self-consistent iterations.
Self consistency -Real space spin correlations.In this paper we have implemented a self-consistent CP-AFQMC approach, described in [64].An example of the convergence of the same-spin density-density correlations ⟨n A i n A j ⟩ is given in Fig. 5 to illustrate this process.The calculation is initiated from a Slater determinant given by the U eff = 3.5 Hartree-Fock 3-2 AFM solution and the one-body density matrix is obtained from QMC.We diagonalize the one-body density matrix and find the eigenvectors which correspond to the largest N σ eigenvalues.Those eigenvectors are used to construct a new Slater determinant, which is then used as the initial Slater determinant for the next iteration of simulation.The number in the legend of Fig. 5 indicates the self-consistent iteration value.Blue right triangles ("5" in the legend) and blue stars ("6") are almost on top of each other, implying the density-density correlations converge after ∼ 6 self-consistent iterations.
As a check that the results are independent of initialization, we also start our simulation with another initial Slater determinant, obtained from the U eff = 10, 3-2 AFM Hartree-Fock solution.The results (red stars) agree with the one starting from the U eff = 3.5 Hartree-Fock (blue star).All results presented in the main text are given by the last self-consistent (fully converged) iteration.
Finite lattice size effects -Energy.In Fig. 6, we explore finite lattice size effects in the Hartree-Fock calculations.First, similar to Fig. 4(b), we show the total energy difference per site ∆E/N s given by Hartree-Fock as a function of U for a 12 × 12 square lattice in panel (a).The 3-2, 3-4, and 3-3 AFM phase has the lowest energy for  2)).We can extrapolate the energy to the thermodynamic limit 1/N s → 0. There is no crossing between different AFM phases, implying the energy ordering between the three magnetic phases on the 12 × 12 lattice in Fig. 6(a) still holds in the thermodynamic limit.We have also verified that this conclusion is robust to the application twisted boundary conditions to the lattice and taking the average of energies for different twist angles.This method effectively samples a finer mesh of momentum points and hence larger spatial sizes [86][87][88].

FIG. 1 .
FIG. 1. Phase diagram of the 1/3-filling SU(3) Fermi-Hubbard Model.(a-d): schematic density correlations for one of the three spin species in different phases on a N = 12 × 12 square lattice.(a) At small U < Uc1, the system resides in a paramagnetic (PM) phase, as indicated by the uniform gray color map; (b) For Uc1 < U < Uc2 a non-uniform pattern emerges which has period 3 in one direction and 2 in the other, which we denote by "3-2 AFM" phase; (c) Uc2 < U < Uc3 is similarly denoted as "3-4 AFM"; and (d) finally, for U > Uc3, the order is a "3-3 AFM".(e) Hartree-Fock phase diagram: The quantum phase transition points are U MF c1 ≈ 3.25, U MF c2 ≈ 4.75 and U MF c3 ≈ 5.65.(f) AFQMC phase diagram: Uc1 ≈ 5.5, Uc2 ≈ 7.7 and Uc3 ≈ 8.3.The question mark indicates that, although 3-3 has a slightly lower energy, it cannot be clearly distinguished from 3-4 in our calculations.

Fig. 1 (
b) has a period of 3 in one direction and 2 in the other direction, thus it is denoted as a 3-2 AFM.Similarly, patterns (c) and (d) are denoted 3-4 and 3-3 AFM respectively.Within the Hartree-Fock treatment, we have U MF c1 ∼ 3.25 below which the system is in the paramagnetic phase.When 3.25 ≲ U ≲ 4.75, the 3-2 AFM has a lower energy with respect to the other two states, while the 3-4 AFM is the ground state for the intermediate region 4.75 ≲ U ≲ 5.65.However, the energy gap to other orders in this intermediate regime is relatively small.The 3-3 AFM is the mean-field ground state for the large U region U ≳ 5.65.The corresponding phase diagram from AFQMC is shown in Fig. 1(f).As expected, the metal-insulator phase transition happens at a larger U c1 ∼ 5.5 > U MF c1

FIG. 4 .
FIG. 4. Determining the transition points between phases.(a) The 3-2 AFM structure factor vs. U is shown for three lattice sizes.S(2π/3, π) is small and nearly independent of the system size in the metallic phase but grows proportionally to system size in the 3-2 AFM phase.The boundary between the two regimes gives the metal-insulator phase transition point Uc1.The inset shows S(π, π) for 1/2 filling, which exhibits a critical point at U = 0, in contrast with the finite Uc1 at 1/3 filling.(b) Energy differences with respect to the 3-4 AFM phase.Red, green curves represent 3-2, 3-4, and 3-3 AFM states respectively.Insets (b2) and (b3) show complementary results on structure factors to help validate the location of Uc2 and Uc3.In each, the calculations use a trial wave function constructed from the superposition of the two states involved, and the two corresponding structure factors are shown.Different symbol shapes represent results from different random seeds.
c3and U ≳ U M F c3 respectively.In order to detect whether the small energy difference at strong interactions is due to the finite lattice size effect, we choose U = 4 (3-2 phase), U = 5 (3-4 phase) and U = 7 (3-3 phase) and present the energy per site E/N s as a function of 1/N s in Fig.6(b).Different symbols in the legend represent the local minimum 3-2, 3-4, 3-3 AFM states.Lattice sizes are labeled next to the data points in Fig.6(b)-(