SU(3) Fermions on the Honeycomb Lattice at 1/3-Filling
Sangwoo S. Chung Philippe Corboz









 Institute for Theoretical Physics Amsterdam and Delta Institute for Theoretical Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

  Institute for Theoretical Physics Amsterdam and Delta Institute for Theoretical Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands


   SU(N) symmetric fermions on a lattice, which can be realized in ultracold-atom-based quantum simulators, have very promising prospects for realizing exotic states of matter. Here we present the ground state phase diagram of the repulsive SU(3) Hubbard model on a honeycomb lattice at 1/3 filling obtained from  infinite projected entangled pair states tensor network calculations. In the strongly interacting limit the ground state has plaquette order. Upon decreasing the interaction strength U/t we find a first order transition at U/t=7.2(2) into a dimerized, color-ordered state, which extends down to U/t=4.5(5) at which the Mott transition occurs and the ground state becomes uniform. Our results may serve as a prediction and benchmark for future quantum simulators of SU(3) fermions.


SU(3) fermions on the honeycomb lattice at 1/3-filling
    Philippe Corboz
    December 30, 2023
======================================================





§ INTRODUCTION

Quantum simulators based on ultracold fermionic alkaline earth atoms in an optical lattice offer a very exciting route to realize exotic phases of symmetric fermions which can be described by an Hubbard model <cit.>.
By exploiting different nuclear spin states of the alkaline earth atoms, systems with up to N=10 different flavors (or colors) of fermions can be obtained. Thanks to an almost-perfect decoupling of the nuclear spin from the electronic angular momentum, the interactions are essentially independent of the nuclear spin, giving rise to a large symmetry. 

Another route to realization of systems (or SO(5) <cit.>) is provided by exploiting different hyperfine states of alkali atoms <cit.>.

On the theory side, substantial progress has been made in recent years in predicting the ground states of these systems in the Mott insulating state in the strongly interacting limit, where the system is effectively described by an Heisenberg model. On two-dimensional (2D) lattices, a rich variety of exotic ground states has been found, including various states with color order <cit.>,  generalized valence-bond solids <cit.>, algebraic N-flavor spin liquids <cit.>, and chiral spin liquids <cit.>. 

An even richer behavior can be expected in  Hubbard models, where enhanced charge fluctuations in the intermediate interaction range may give rise to new interesting phases. Accurate calculations would also be desirable for a more direct comparison and prediction for future quantum simulators. However,  the Hubbard model is also substantially more challenging than the Heisenberg model, especially in the strongly correlated regime in two dimensions. While there has been progress in recent years based on dynamical mean-field theory<cit.>,  variational Monte Carlo <cit.>, Quantum Monte Carlo for exceptional cases without a sign problem <cit.>, and other approximate approaches <cit.>,  obtaining controlled and systematic predictions using unbiased numerical approaches remains a central challenge in general. In contrast, in the one-dimensional case  Hubbard models have been accurately studied using matrix product states <cit.>.

In this paper we demonstrate that the 2D Hubbard model has become within reach of state-of-the-art 2D tensor networks, which are a generalization of matrix product states to higher dimensions. Two-dimensional tensor networks have already been successfully applied in the Heisenberg case <cit.> and the SU(2) Hubbard model <cit.> (and many other strongly correlated systems; see, e.g., Refs. Corboz2014b,Niesen2017a,Liao2017,Haghshenas2018,chen18,Lee2018,Jahromi2018 and references therein). 

Here we consider the SU(3) Hubbard model on a honeycomb lattice and determine the nature of the Mott insulating phases as a function of the interaction strength, summarized in the phase diagram shown in Fig. <ref>. Our main findings include  a plaquette ground state in the strongly interacting limit, consistent with the result for the Heisenberg case <cit.>, whereas in the intermediate interaction range another state which has a dimer order with a coexisting color order gets stabilized by enhanced charge fluctuations.














§ MODEL
 
The Hubbard model is given by 


    Ĥ=-t ∑_⟨ i,j⟩ ,αĉ_iα^†ĉ_jα+H.c.+U∑_i,α<βn̂_iαn̂_iβ,


where α, β = 1...N label the N different colors (or flavors) of fermions, i and j denote the lattice sites, and ⟨ i,j⟩ runs over pairs of nearest neighbor sites. The operator ĉ_iα^† (ĉ_iα) creates (annihilates) a fermion of type α at site i, and n̂_iα≡ĉ_iα^†ĉ_iα is the particle number operator of type α.  Here we consider the model for N=3 at 1/3 filling (one particle per site) on a honeycomb lattice and study its phase diagram as a function of U/t.








§ METHOD

An infinite projected entangled pair state (iPEPS) is a variational tensor network ansatz to efficiently represent 2D ground states of local Hamiltonians in the thermodynamic limit <cit.>.  It consists of a unit cell of tensors, with one tensor per lattice site, which is periodically repeated on the lattice. Each tensor has one physical index of dimension d representing the local Hilbert space of a lattice site (here d=8), and Z auxiliary indices with a bond dimension D where Z is the coordination number of the lattice (Z=3 for the honeycomb lattice), and D controls the accuracy of the ansatz. 
 
For our simulations, we use the fermionic iPEPS ansatz from Ref. corboz2010 (see also Refs. kraus2010,Barthel2009,Corboz10_fmera,pineda2010,Corboz09_fmera,shi2009) with unit cells up to 18 sites. The corner transfer matrix method <cit.> is used to contract the infinite 2D network of tensors, where we map the honeycomb lattice onto a brick-wall square lattice <cit.>. The optimal variational parameters to approximate the ground state are obtained based on an imaginary time evolution <cit.>. The truncation of a bond index, required at each imaginary time step, is done using both the so-called simple-update <cit.> and the more accurate (but computationally more expensive)  full-update scheme <cit.>.  To increase the efficiency we use  tensors with  U(1)× U(1)× U(1) symmetry, which significantly reduces the number of variational parameters, enabling us to reach large bond dimensions up to D=28. 










§ RESULTS

We first focus on the strongly interacting limit (U/t=100) where we find results which are consistent with those from the previously studied SU(3) Heisenberg model, <cit.> namely, a ground state with plaquette order in a six-site unit cell depicted in Fig. <ref>(a). The thickness of the bonds is proportional to the square magnitude of the corresponding local energies (called bond energies E_b), and the pie charts show the local densities of each color, which sum to n=1. In this case, all color densities are equal (i.e. n_α=1/3) showing the absence of color order [i.e. SU(3) symmetry is not broken], and strong energy bonds are found around plaquettes consistent with  SU(3) singlets which are predominantly formed around the hexagons and thereby break translational symmetry. 

As in the Heisenberg case <cit.> we find a competing low-energy state in a 18-site unit cell which is dimerized and has color order, shown in Fig. <ref>(b). Each dimer (in orange)  has two colors that are dominant at each end, and at the four sites surrounding a dimer the remaining third color is dominant. 


 
 

In Fig. <ref>(c) we compare the energies of the two states as a function of 1/D using both simple- and full-update optimizations. For large bond dimensions we clearly find that the plaquette state is energetically  favorable. The same conclusion is reached by plotting the energy as a function of the truncation error w, shown in Fig. <ref>(d), which often provides a better parameter for an extrapolation of the energy to the exact infinite-D limit <cit.>. The extrapolated energy of the dimer state, E_s=-0.043936(9), is clearly higher than that of the plaquette state, E_s=-0.044184(9).

In Fig. <ref>(e) we present  results for the local ordered moment m, given by 

    m_i=√(3/2∑_α=1^3(⟨ n_i,α⟩ -1/3)^2),

averaged over all lattice sites i in the unit cell, which is finite in the dimer state and vanishing in the plaquette state. Figure <ref>(f) shows the difference in energy between the highest and the lowest bond energy in the unit cell, Δ E = max(E_b) - min(E_b), which is finite for both states. Since the large-D results between simple and full updates are found to be similar for both states, we continue using the former, since it allows us to reach larger D values with substantially smaller computational effort.


We next move away from the Heisenberg limit by reducing U/t which enhances charge fluctuations. Interestingly, we observe that the energy difference between the two competing states gets smaller with decreasing U/t [see Figs. <ref>(a)-<ref>(d), and around a critical value of U_c/t = 7.2(2)[This value, which is based on an extrapolation in the truncation error w, is compatible with the one obtained using a 1/D extrapolation, U/t=7.4(7).]  the energies of the two states intersect [Fig. <ref>(e)] so that the dimer state actually becomes the ground state for U<U_c [We also tested other unit cell sizes but could not find other competing states, except deformed versions of the color-ordered, dimerized state which has a higher energy]. Consequently, at the critical point, both m and Δ E  exhibit a discontinuity [Fig. <ref>(f)], characteristic for the first-order nature of the transition.

 Thus, enhanced charge fluctuations present at a lower U/t are found to favor the dimerized state over the plaquette state. Intuitively this is because in the latter the kinetic energy is mostly concentrated locally around the plaquettes, whereas in the dimer state charges are delocalized more uniformly in the system [which can be seen from the smaller value of Δ E in the dimerized state than in the plaquette state; see Fig. <ref>(f)], which is favored as we decrease the interaction (with the fermions being completely delocalized in the non-interacting limit).







To get more insights into the energetics of the two states we present a comparison of the different energy contributions in Fig. <ref>  for U/t=100, U/t=9, and U/t=6. The kinetic  term is split into two parts. The first part is the one relevant for the superexchange processes  in the Heisenberg (large-U/t) limit, i.e., matrix elements between two singly occupied sites |1_α,1_β⟩ and an empty and doubly occupied site  |0,1_α 1_β⟩ (or |1_α 1_β ,0⟩), where α≠β denote the colors of the two fermions. The second term E_kin2 includes all remaining kinetic contributions. In the Heisenberg limit the latter can be neglected (see top panel in Fig. <ref>), and the plaquette state wins over the dimerized state thanks to a slightly lower exchange energy. When U/t is lowered the energy cost to form doubly occupied sites decreases, and enhanced charge fluctuations lead to a gain in kinetic energy, including also contributions which are not included in the exchange term. These contributions are clearly stronger in the dimer state, at the expense of a larger positive contribution from the on-site repulsion (see middle panel in Fig. <ref>). Eventually, for even smaller U/t these kinetic contributions become dominant compared to the on-site repulsion, and the dimerized state wins over the plaquette state (bottom panel in Fig. <ref>). 






When U/t is lowered even further one can expect that the system will eventually undergo a transition from the Mott insulating dimer state into a conducting (or superconducting) state at weak U/t. In the non-interacting case, U/t=0, the model reduces to a tight-binding model of three independent species of fermions. In contrast to the half-filled case where the Fermi level crosses the Dirac nodes (as, e.g., in graphene), the Fermi surface at 1/3 filling is one-dimensional manifold in  momentum space, which is known to lead to a multiplicative logarithmic correction to the area law of the entanglement entropy <cit.>. Thus, in the low-U/t limit the states are strongly entangled  and therefore very challenging to  accurately  represent by a tensor network ansatz. This is also reflected in the increasing truncation error with decreasing U/t (for a fixed value of D; see Fig. <ref>).




To obtain an estimate of the Mott transition we study the stability of the dimer phase upon lowering U/t. The  energy difference between the highest and the lowest bond energies shown in Fig. <ref>(a) gets strongly suppressed with increasing D for U/t ≤ 5. Taking a linear extrapolation in 1/D as a rough estimate for the value in the infinite-D limit, we find that Δ E — where Δ E > 0 indicates breaking of translational symmetry — vanishes for U/t=4. Also the square of the local ordered moment — where m^2>0 indicates breaking of SU(3) color symmetry — shown in Fig. <ref>(b) extrapolates to a vanishingly small value in the infinite-D limit for U/t=4.  From these results we conclude that the dimer state gets unstable between U/t=4 and U/t=5, i.e. that the ground state becomes uniform for U/t<U_c/t=4.5(5). Furthermore, we also find that the charge gap of the uniform state, obtained using a two-site unit cell iPEPS ansatz, clearly vanishes for U/t=4 [see Fig. <ref>(c)]. For U/t=4.5 and D=28 the charge gap is still visible by the plateau in n(μ), but its size gets suppressed with increasing D. This suggests that the Mott transition also takes place at U_c (i.e.,  there is no clear indication of an additional intermediate quantum spin liquid insulating phase). [We note that determining the nature of the transition point between the dimer and uniform state is challenging. Due to the increased entanglement with decreasing U/t it is difficult to distinguish a continuous transition from a weak first-order transition here]. 




 











§ CONCLUSIONS

Using iPEPS simulations pushed to large bond dimensions, up to D=28, we have investigated the ground-state phase diagram of the SU(3) Hubbard model on a honeycomb lattice at 1/3 filling, focusing on the Mott insulating phases. 
In the large-U/t (Heisenberg) limit we found a plaquette ground state which breaks lattice symmetry but preserves symmetry. For an intermediate interaction strength, 4.5(5)≤ U/t ≤ 7.2(2), the ground state is a dimerized, color-ordered state which breaks both and lattice symmetry. The order parameters of the dimerized state vanish for U ≤ 4.5(5), indicating that the ground state becomes uniform, compatible with a conducting (or superconducting) state in the weakly interacting limit. 

This work demonstrates that the 2D Hubbard model in the strongly correlated regime has become within reach of state-of-the-art tensor network methods, offering a systematic and controlled way to predict the phases of these challenging models. By making use of global symmetries it is also  possible to study systems with a larger N, despite the  large local dimension of d=2^N.
 
Our results may serve as a prediction for future quantum simulators based on SU(3) ultra-cold atoms in an optical honeycomb lattice. Since the plaquette phase breaks a discrete lattice symmetry we can expect that it extends also to finite temperatures, so that the phase can be realized in experiments at sufficiently low temperatures,  offering also an interesting possibility to study a finite-temperature phase transition in a 2D quantum simulator. While the finite color order of the dimer phase can only exist at zero temperature in two dimensions, the dimerization may in principle set in already at finite temperature (without coexisting color order), or occur simultaneously with the color order at zero temperature. Recently developed tensor network approaches for finite-temperature simulations <cit.> may provide further insights into the critical temperatures in the future.
  
 
 






S.S.C. acknowledges helpful discussions with I. Niesen. This project received funding from the European Research Council (ERC) under the European Union's Horizon 2020  research and innovation program (Grant Agreement No. 677061). This work is part of the  consortium, a program of the Netherlands Organization for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW).







apsrev4-1

