SU(3) Fermions on the Honeycomb Lattice at 1/3-Filling

TL;DR

This paper maps the phase diagram of SU(3) fermions on a honeycomb lattice at 1/3 filling, revealing transitions from plaquette order to dimerized states and a Mott insulator, using tensor network calculations.

Contribution

It provides the first detailed ground state phase diagram of the SU(3) Hubbard model on a honeycomb lattice at 1/3 filling with tensor network methods.

Findings

Plaquette order in the strongly interacting limit

First order transition at U/t=7.2(2) to a dimerized state

Mott transition at U/t=4.5(5) leading to a uniform ground state

Abstract

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.

Equations

$$\hat{H}=-t\sum _{⟨i,j⟩,\alpha }{\hat{c}}_{i\alpha }^{†}{\hat{c}}_{j\alpha }+H.c.+U\sum _{i,\alpha <\beta }{\hat{n}}_{i\alpha }{\hat{n}}_{i\beta },$$

$$m_i=\sqrt{\frac{3}{2}\sum _{\alpha =1}^3{\left(⟨n_{i,\alpha }⟩-\frac{1}{3}\right)}^2},$$

Full text

I Introduction

Quantum simulators based on ultracold fermionic alkaline earth atoms in an optical lattice offer a very exciting route to realize exotic phases of SU($N$) symmetric fermions which can be described by an SU($N$) Hubbard model Cazalilla et al. (2009); Gorshkov et al. (2010); Taie et al. (2012); Scazza et al. (2014); Zhang et al. (2014); Cazalilla and Rey (2014); Hofrichter et al. (2016); Ozawa et al. (2018). 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 SU($N$) symmetry. Another route to realization of SU($N$) systems (or SO(5) Wu et al. (2003)) is provided by exploiting different hyperfine states of alkali atoms Wu et al. (2003); Honerkamp and Hofstetter (2004); Rapp et al. (2008).

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 SU($N$) Heisenberg model. On two-dimensional (2D) lattices, a rich variety of exotic ground states has been found, including various states with color order Tóth et al. (2010); Bauer et al. (2012); Corboz et al. (2011a), generalized valence-bond solids Read and Sachdev (1989, 1990); Harada et al. (2003); Arovas (2008); Hermele and Gurarie (2011); Zhao et al. (2012); Corboz et al. (2012a, 2013); Nataf et al. (2016), algebraic $N$-flavor spin liquids Affleck and Marston (1988); Marston and Affleck (1989); Assaad (2005); Corboz et al. (2012b), and chiral spin liquids Hermele et al. (2009); Hermele and Gurarie (2011); Szirmai et al. (2011).

An even richer behavior can be expected in SU($N$) 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 theoryGorelik and Blümer (2009); Inaba et al. (2010); Inaba and Suga (2012); Yanatori and Koga (2016); Koga and Yanatori (2017); Del Re and Capone (2018), variational Monte Carlo Assaad (2005); Rapp and Rosch (2011); Boos et al. (2018), Quantum Monte Carlo for exceptional cases without a sign problem Lang et al. (2013); Cai et al. (2013); Wang et al. (2014), and other approximate approaches Honerkamp and Hofstetter (2004); Hasunuma et al. (2016), obtaining controlled and systematic predictions using unbiased numerical approaches remains a central challenge in general. In contrast, in the one-dimensional case SU($N$) Hubbard models have been accurately studied using matrix product states Buchta et al. (2007); Zhao et al. (2007); Szirmai et al. (2008); Capponi et al. (2008); Molina et al. (2009); Roux et al. (2009); Ulbricht et al. (2010); Nonne et al. (2010); Manmana et al. (2011); Barcza et al. (2015); Bois et al. (2015); Barbarino et al. (2016); Szirmai et al. (2017); Weichselbaum et al. (2018).

In this paper we demonstrate that the 2D SU($N$) 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 SU($N$) Heisenberg case Corboz et al. (2011a); Bauer et al. (2012); Zhao et al. (2012); Corboz et al. (2012a, b, 2013); Nataf et al. (2016) and the SU(2) Hubbard model Zheng et al. (2017) (and many other strongly correlated systems; see, e.g., Refs. Corboz and Mila, 2014; Niesen and Corboz, 2017; Liao et al., 2017; Haghshenas et al., 2018; Chen et al., 2018; Lee and Kawashima, 2018; Jahromi and Orús, 2018 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. 1. Our main findings include a plaquette ground state in the strongly interacting limit, consistent with the result for the Heisenberg case Zhao et al. (2012); Corboz et al. (2013), whereas in the intermediate interaction range another state which has a dimer order with a coexisting color order gets stabilized by enhanced charge fluctuations.

II Model

The SU($N$) Hubbard model is given by

[TABLE]

where $\alpha,\beta=1...N$ label the $N$ different colors (or flavors) of fermions, $i$ and $j$ denote the lattice sites, and $\left\langle i,j\right\rangle$ runs over pairs of nearest neighbor sites. The operator $\hat{c}_{i\alpha}^{\dagger}$ ($\hat{c}_{i\alpha}$) creates (annihilates) a fermion of type $\alpha$ at site $i$, and $\hat{n}_{i\alpha}\equiv\hat{c}_{i\alpha}^{\dagger}\hat{c}_{i\alpha}$ is the particle number operator of type $\alpha$. 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$.

III 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 Verstraete and Cirac (2004); Jordan et al. (2008); Nishio et al. (2004). 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. Corboz et al., 2010a (see also Refs. Kraus et al., 2010; Barthel et al., 2009; Corboz et al., 2010b; Pineda et al., 2010; Corboz and Vidal, 2009; Shi et al., 2009) with unit cells up to $18$ sites. The corner transfer matrix method Baxter (1976); Nishino and Okunishi (1996); Orús et al. (2009); Corboz et al. (2011b, 2014) is used to contract the infinite 2D network of tensors, where we map the honeycomb lattice onto a brick-wall square lattice Corboz et al. (2012b). The optimal variational parameters to approximate the ground state are obtained based on an imaginary time evolution Corboz et al. (2010a); Phien et al. (2015). The truncation of a bond index, required at each imaginary time step, is done using both the so-called simple-update Jiang et al. (2008) and the more accurate (but computationally more expensive) full-update scheme Corboz et al. (2010a); Phien et al. (2015). To increase the efficiency we use tensors with $U(1)\times U(1)\times U(1)$ symmetry, which significantly reduces the number of variational parameters, enabling us to reach large bond dimensions up to $D=28$.

IV 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, Zhao et al. (2012); Corboz et al. (2013) namely, a ground state with plaquette order in a six-site unit cell depicted in Fig. 2(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_{\alpha}=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 Lee and Yang (2012); Corboz et al. (2013) we find a competing low-energy state in a 18-site unit cell which is dimerized and has color order, shown in Fig. 2(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. 2(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. 2(d), which often provides a better parameter for an extrapolation of the energy to the exact infinite-$D$ limit Corboz (2016). 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. 2(e) we present results for the local ordered moment $m$, given by

[TABLE]

averaged over all lattice sites $i$ in the unit cell, which is finite in the dimer state and vanishing in the plaquette state. Figure 2(f) shows the difference in energy between the highest and the lowest bond energy in the unit cell, $\Delta 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. 3(a)-3(d), and around a critical value of $U_{c}/t=7.2(2)~{}$111This 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. 3(e)] so that the dimer state actually becomes the ground state for $U<U_{c}$ 222We 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 $\Delta E$ exhibit a discontinuity [Fig. 3(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 $\Delta E$ in the dimerized state than in the plaquette state; see Fig. 2(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. 4 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_{\alpha},1_{\beta}\rangle$ and an empty and doubly occupied site $|0,1_{\alpha}1_{\beta}\rangle$ (or $|1_{\alpha}1_{\beta},0\rangle$), where $\alpha\neq\beta$ 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. 4), 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. LABEL:Fig:Econtr). 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. LABEL:Fig:Econtr).

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 Wolf (2006); Gioev and Klich (2006). 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. 3).

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. 5(a) gets strongly suppressed with increasing $D$ for $U/t\leq 5$. Taking a linear extrapolation in $1/D$ as a rough estimate for the value in the infinite-$D$ limit, we find that $\Delta E$ — where $\Delta 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. 5(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. 5(c)]. For $U/t=4.5$ and $D=28$ the charge gap is still visible by the plateau in $n(\mu)$, 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). 333We 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.

V 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 SU($N$) symmetry. For an intermediate interaction strength, $4.5(5)\leq U/t\leq 7.2(2)$, the ground state is a dimerized, color-ordered state which breaks both SU($N$) and lattice symmetry. The order parameters of the dimerized state vanish for $U\leq 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 SU($N$) 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 Czarnik et al. (2012); Czarnik and Dziarmaga (2014, 2015); Czarnik et al. (2016); Kshetrimayum et al. (2019); Czarnik et al. (2019) may provide further insights into the critical temperatures in the future.

Acknowledgements.

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 Delta-ITP consortium, a program of the Netherlands Organization for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW).