Flat-band Ferromagnetism of the SU$(N)$ Hubbard Model on the Tasaki Lattice

TL;DR

This paper studies the phase transition to ferromagnetism in the SU(N) Hubbard model on Tasaki lattices, using rigorous proofs and numerical simulations to explore how the transition depends on dimensionality and spin degrees of freedom.

Contribution

It provides the first rigorous proof of the para-ferro transition in one dimension and extensive numerical analysis for higher dimensions and larger N, revealing the impact of spin on critical behavior.

Findings

Rigorous proof of para-ferro transition in 1D lattice.

Critical density exceeds standard percolation and increases with N.

Effective repulsive interaction strengthens with larger N.

Abstract

We investigate the para-ferro transition of the repulsive SU($N$) Hubbard model on one- and two-dimensional Tasaki lattices. Under certain restrictions for constructing localized many-particle ground states of flat-band ferromagnetism, the quantum strongly correlated electrons model is mapped to a classical statistical geometric site-percolation problem, where the nontrivial weights of different configurations must be considered. We prove rigorously the para-ferro transition for the SU($N$) Hubbard model on one-dimensional Tasaki lattice by transfer-matrix method, and numerically verify it for the model with symmetries up to SU($10$). In two dimensions, we numerically investigate the phase transition of SU($3$), SU($4$) and SU($10$) Hubbard models by Metropolis Monte Carlo simulation. And we find the critical density exceeds that of standard percolation, and it increases with spin degrees of freedom, implying that the effective repulsive interaction becomes stronger for larger $N$. We further rigorously prove the flat-band ferromagnetism of the SU$(N)$ Hubbard model, when the number of particles $N_e$ equals to the degeneracy $N_d$ of the lowest band in the single-particle energy spectrum.

Figures

Equations

$$H=\sum _{i,j\in \Lambda }{t}_{i,j}\left(c_{i,\sigma }^{†}{c}_{j,\sigma }+\text{mbox}H.c.\right)+U\sum _i\sum _{\sigma \ne {\sigma }^{\prime }}{n}_{i,\sigma }{n}_{i,{\sigma }^{\prime }},$$

$$d_{\text{SU}(N)}(|C|)=(N+|C|-1)!/(|C|!(N-1)!).$$

$$W(q)=\prod _{i=1}^{M_q}{e}^{\mu |C_i|}{d}_{\text{SU}(N)}(|C_i|),$$

$$⟨O⟩=\frac{\sum _q{O}_{q}W(q)}{\sum _{q}W(q)}.$$

$$S^2(|C|)=\frac{(N-1)(N+|C|)|C|}{2N}.$$

$$S_q^2$$

$$=\sum _{l=1}^n{N}_q(l)S^2(l)=\mathcal{N}\sum _{l=1}^n{n}_q(l)S^2(l),$$

$$n(l,\text{SU}(N))=f_{\text{SU}(N)}(z)\frac{(l+N-1)!}{l!}{\left(\frac{z}{{\lambda }_{\text{max}}(z)}\right)}^l,$$

$$\frac{S^2}{{\mathcal{N}}^2}=\sum _{l=1}^n\frac{S_l^2}{\mathcal{N}}{f}_{\text{SU}(N)}(z)\frac{(l+N-1)!}{l!}{\left(\frac{z}{{\lambda }_{\text{max}}}\right)}^l,$$

$$H=H_0+H_{\text{int}},$$

$$a_{u,\sigma }^{†}=\sum _{i\in \Lambda }{\phi }_i^u{c}_{i,\sigma }^{†},$$

$$\varphi_{i}^{u}=\left\{\begin{array}[]{rcl}1,&&{i=u+\delta;}\\ -\sqrt{c},&&{i=u;}\\ 0,&&\textrm{otherwise.}\end{array}\right.$$

$$\left[H_0,a_{u,\sigma }^{†}\right]={\epsilon }_0{a}_{u,\sigma }^{†}$$

$$|{\Phi }_{A,N}⟩=\prod _{u\in A}{a}_{u,N}^{†}|0⟩$$

$$|{\Phi }_{A,\{m_i\}}⟩=\left\{\prod _{i=1}^M{\left[\prod _{k=1}^{m_i}{\mathbf{S}}_{k,i}^{-}\right]}^{\prime }\right\}|{\Phi }_{A,N}⟩,$$

$$|\Phi {⟩}_{\text{GS}}=\sum _{A_1,\cdots A_N}f(A_1,\cdots A_N)\prod _{u_1\in A_1}{a}_{u_1,1}^{†}\cdots \prod _{u_N\in A_N}{a}_{u_N,N}^{†}|0⟩,$$

$$(G{)}_{u,v}=\sum _{i\in \Lambda }{\phi }_i^u{\phi }_i^v,$$

$${\kappa }_i^u=\sum _{v\in V}(G^{-1}{)}_{u,v}{\phi }_i^v,$$

$$\sum _{i\in \Lambda }{\kappa }_i^u{\phi }_i^v$$

$$\sum _{u\in V}{\kappa }_i^u{\phi }_j^u$$

$${\psi }_{ij}={\psi }_{ji}={\delta }_{i,j}-\sum _{u,v\in V}{\phi }_i^u(G^{-1}{)}_{u,v}{\phi }_j^v.$$

$$b_{u,\sigma }=\sum _{i\in \Lambda }{\kappa }_i^u{c}_{i,\sigma }$$

$$d_{i,\sigma }=\sum _{j\in \Lambda }{\psi }_{ij}{c}_{j,\sigma },$$

$$\left\{b_{u,\sigma },a_{u^{\prime },{\sigma }^{\prime }}^{†}\right\}$$

$$\left\{d_{i,\sigma },a_{u,{\sigma }^{\prime }}^{†}\right\}$$

$$c_{i,\sigma }^{†}$$

$$c_{i,\sigma }$$

$$f(A_1,\cdots ,A_N)=0$$

$$\left|A_{\sigma }\cap C_k\right|=\left|A_{\sigma }^{\prime }\cap C_k\right|$$

$$f(A_1,\cdots ,A_N)=f(A_1^{\prime },\cdots ,A_N^{\prime }).$$

$$\mathbf{T}$$

$$\displaystyle=\left(\begin{array}[]{ccccc}1&1&1&\cdots&1\\ z&z&z&\cdots&z\\ z&0&z&\cdots&z\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ z&0&0&\cdots&z\\ \end{array}\right)$$

$$\Xi (z,\mathcal{N})=\text{Tr}{\mathbf{T}}^{\mathcal{N}}.$$

Full text

Keywords : SU($N$), flat-band ferromagnetism, percolation, para-ferro transition

**1.Introduction

**

In condensed matter physics, the interplay of the Coulomb interaction with Pauli principle is a fundamental problem, from where many intriguing phases can appear. One important problem is the emergence of ferromagnetism in two or higher dimensions, which can date back to the early age of quantum theory and reveal the importance of many-body effect. After more than one century of investigation, a few special cases have been solved by either exact results Nagaoka ; Lieb-theorem ; Mielke1 ; Mielke2 ; TasakiPRL ; Mielke-Tasaki ; Tasaki-review ; YiLi ; Li-Nagaoka or sign-problem-free quantum Monte Carlo simulation Li-QMC-FM . However, a general solution to this problem is still lack, partly due to the notorious sign problem of fermions book-Lieb ; book-Flliott . As an interesting and important limiting case, a system with a flat band in the single-electron energy spectrum may favor the stabilization of ferromagnetism as it reduces the kinetic energy cost for spin polarization to zero due to the band flatness. For instance, Mielke and Tasaki prove flat-band ferromagnetism in a class of lattices, based on a systematic route called “cell construction” Mielke1 ; Mielke2 ; TasakiPRL ; Mielke-Tasaki . After that, flat-band ferromagnetism has attracted a lot of attention and has been connected with the discussion of topological and magnetic phases ShengNatComm ; SunPRL ; Tang2011 ; Katsura-Tasaki ; Scaffidi ; Essafi ; Tsai ; Goldman ; Shizhong-percolation .

The previous investigation on flat-band ferromagnetism mainly focus on fermionic systems possessing the SU($2$) spin-rotational symmetry. The SU($N>2$) generalization of flat-band ferromagnetism is not only of theoretical interests Affleck1 ; Affleck2 ; Sachdev1 ; Sachdev2 ; Wu-SO5 ; Wu-NatPhy ; Wu-view ; Cazalilla ; NiePRA , but also of experimental relevance in the context of cold atomic gases  Takahashi_PRL ; Takahashi_Nat ; Taie ; Bloch , where the high degrees of freedom of nuclear spin provides a fascinating playground to simulate the SU($N$) fermionic Hubbard model Cappellini ; ZhangSci ; Pagano2014 ; Scazza . By loading ultracold alkaline-earth-metal(-like) atoms in different optical lattices, the quantum simulation of various models has been theoretically proposed and further experimentally realized Sr-Rice ; Fallani ; Takahashi_PRL ; Takahashi_Nat , where the spin symmetries can be as high as SU($6$) for 173Yb TaiePRL20010 ; Takahashi_Nat ; Taie ; Hofrichter and SU($10$) for 87Sr Sr-Rice ; Sr-Desalvo .

In this paper, we investigate the para-ferro magnetic transition and flat-band ferromagnetism of the SU($N$) Hubbard model on a special class of decorated cubic lattices in one and two dimensions, which are illustrated in Fig. 1. The existence of ferromagnetism with SU($2$) spin symmetry on this lattice, later referred as Tasaki lattice for short, is first considered by Mielke and Tasaki Mielke1 ; Mielke2 ; TasakiPRL ; Mielke-Tasaki , and further by Maksymenko et al. MaksymenkoPRL , who introduced a mapping from the quantum Hubbard model to a classical statistical geometric site-percolation problem. Here, we first prove a generalized version of this mapping from the SU($N$) Hubbard model to an $N$-state Pauli-correlated percolation (PCP) problem, where an SU($N$) correlated weight for each configuration is considered. Then we obtain the exact results in one-dimensional ($1$D) Tasaki lattice for the SU($N$) Hubbard model, and further verify the para-ferro magnetic transition numerically. For the two-dimensional ($2$D) case, we implement importance samplings by Metropolis Monte Carlo simulation and determine the para-ferro magnetic transition numerically. Moreover, we prove rigorously the existence of flat-band ferromagnetism of the SU$(N)$ Hubbard model, when the number of particles $n$ equals to the degeneracy $\mathcal{N}$ of the lowest band in the single-particle energy spectrum.

**2.Mapping from the SU($N$) Hubbard model to an $N$-state Pauli-correlated percolation problem

**

We start with the SU$(N)$ Hubbard model with repulsive interaction $U>0$ on a Tasaki lattice $\Lambda$,

[TABLE]

where ${\sigma}=\{1,\cdots,N\}$ labels the spin (or equivalently flavor or color) of a fermion. The operator $c_{i,\sigma}$ and its adjoint $c_{i,\sigma}^{{\dagger}}$ satisfy the fermionic anticommutation relations $\{c_{i,\sigma}^{{\dagger}},c_{j,\sigma^{\prime}}\}=\delta_{ij}\delta_{\sigma\sigma^{\prime}}$ and $\{c_{i,\sigma},c_{j,\sigma}\}=\{c_{i,\sigma}^{\dagger},c_{j,\sigma}^{\dagger}\}=0$, and the number operator $n_{i,\sigma}\equiv c_{i,\sigma}^{\dagger}c_{i,\sigma}$. The Tasaki lattice $\Lambda$ is a $d$-dimensional hypercubic lattice with decorated sites located at the middle of the nearest links of the hypercubic lattice Tasaki-review . The hopping pattern, therefore, is composed by the nearest-neighbor ones ($t_{1}$) on the original hypercubic lattice and the nearest-neighbor ones ($t_{2}$) after the decoration, as shown in Fig. 1.

The lowest band of single-particle energy spectrum of the Tasaki lattice is completely dispersionless as long as $t_{2}=\sqrt{c}\,t_{1}>0$, in which $c$ is the coordination number of the corresponding undecorated lattice ($c=2$ and $c=4$ for $1$D and $2$D lattices, respectively). The existence of flat band is due to destructive interferences of particle hoppings, implying that we can construct the localized single-particle eigenstate in each trapping cell Schulenburg-PRL2002 ; Zhitomirsky-Tsunetsugu ; Derzhko2007summary ; DerzhkoPRB2007 , whose wave function only overlaps with that of particles in adjacent cells, as shown by dark grey areas in Fig. 1. When non-interacting particles are filled into the flat band, the ground states are highly degenerate respective to the spin of individual particles. The ferromagnetic state with all particles possess the same spin is obviously one of these degenerate ground states.

Since the repulsive interaction is positive semidefinite, the ground state energy may not decrease as one turns on $U>0$. In the limit of low particle density, each trapping cell is occupied by at most one particle to avoid energy penalty. When the density increases, the wave function of a trapping cell can not avoid overlapping with that of adjacent one. But the energy penalty can be avoided as well if the particles in neighboring trapping cells are in the same symmetric spin state MaksymenkoPRL . By doing so, particles on these linked trapping cells are aligned in spin, and form a polarized “cluster” as shown in Fig. 2a. Namely, the minimization of energy requires all spins of particles belong to the same cluster being aligned, and this is the origin of ferromagnetism.

We then establish a one-to-one correspondence between the states composed of ferromagnetic clusters and the geometric configurations of $n$ $(\leq\mathcal{N})$ particles distributed over $\mathcal{N}$ traps. In comparison to the standard percolation Stauffer ; Isichenko , the ferromagnetic state of a cluster $C$ of size $|C|$ (or the number of particles in cluster $C$) acquires an additional spin degeneracy, which is given by the dimension of the irreducible representation of the SU$(N)$ group Georgi-book ; Katsura ,

[TABLE]

For example, the degeneracy of the ferromagnetic state in the SU($2$) Hubbard model is $2S+1$, where $S$ is the total spin of this cluster. Therefore, for a ferromagnetic cluster of length three and with total spin $S=3/2$ as shown on the left side Fig. 2a, we need to assign a factor of $4$ to the geometric configuration of the cluster composed by three linked trapping cells depicted on the left of Fig. 2b. Note that the irreducible representation is the “permutation symmetric” representation of the SU$(N)$ group. As the ferromagnetic cluster carries the fully symmetric representation of the SU($N$) group, the electron wave function is fully antisymmetric in the spatial space, thereby avoiding the interaction energy penalty. Thus, for a given configuration $q$ composed by a set of linked clusters, a nontrivial weight $W(q)$ must be assigned by taking into account the spin degeneracy of each cluster

[TABLE]

where $e^{\mu}$ is the fugacity to tune the number of particles in a grand-canonical ensemble. Further analysis shows that for a given number of particles, the relative probability for them to decompose into smaller clusters is higher than forming larger clusters, which implies that this nontrivial weight gives rise to an effective repulsive interaction, tending to break up large clusters MaksymenkoPRL . In the following model calculation, it will be shown that the repulsive effect becomes stronger with $N$, such that the formation of clusters with size proportional to the system size, i.e., the percolation transition, becomes more difficult in systems with higher SU$(N)$ spin symmetry.

The key theorem then needs to be proved is that the ground states of the SU$(N)$ Hubbard model when the flat band is partially filled can be represented by the states of ferromagnetic clusters, which can further be mapped into geometric configurations (see Supplementary Material accompanying this paper). Then the quantum many-body system of SU$(N)$ flat-band ferromagnetism is reduced to a classical site-percolation problem with a nontrivial weight. The expectation value of an operator $O$ is the average over all degenerate ground states,

[TABLE]

For the sake of magnetism, we consider the total spin operator $\mathbf{S}$ and focus on the expectation value $\langle\mathbf{S}^{2}\rangle=\sum_{q}\mathbf{S}_{q}^{2}W(q)/\sum_{q}W(q)$, where $\mathbf{S}^{2}$ is the quadratic Casimir operator of the SU$(N)$ group. Because the particles belonging to the same cluster are aligned in the ferromagnetic state, the eigenvalue of $\mathbf{S}^{2}$ for a cluster with size $|C|$ is Ma-book

[TABLE]

The total spin $\mathbf{S}_{q}^{2}$ for a given geometric configuration $q$ can be calculated in two equivalent ways

[TABLE]

where $N_{q}(l)$ is the number of clusters of size $l$, and $n_{q}(l)\equiv N_{q}(l)/\mathcal{N}$ is the normalized fraction with respect to the number of traps $\mathcal{N}$. The average value of the fraction $n(l)=\sum_{q}n_{q}(l)W(q)/\sum_{q}W(q)$ is a characteristic quantity in percolation as it signals the information of cluster size distribution Stauffer ; Isichenko .

3.**Exact results in one dimension

**

As we have discussed in the previous section, the many-body ground states of the SU$(N)$ Hubbard model on flat bands can be mapped to a geometric site-percolation problem with a nontrivial weight. In this section, we consider the 1D Tasaki lattice as an example, and present some exact results obtained by the transfer-matrix method MaksymenkoPRL . The most important step is to construct a one-to-one mapping between the enumeration of all ground states of the Hubbard model and the counting of configurations when populating $n$ particles with spin component $1,2,\cdots,N$ into $\mathcal{N}$ trapping cells. Here, we establish such a mapping by introducing the following rules to populate particles with spin component $\sigma=1,2,\cdots,N$ to the trapping cells DerzhkoPRB2010 :

(1) Each trap $j$ can be either empty or occupied by only one particle with spin component $\sigma$. In other words, for each trap $j=0,1,\cdots,\mathcal{N}-1$, there are $N+1$ trap states $\sigma_{j}=0,1,\cdots,N$, where [math] indicates the trap is empty. Then a ground state of the chain is presented as a sequence of the trap states.

(2) We take the convention $\sigma_{i}\leq\sigma_{i+1}$ to avoid over counting, because only one sequence of two trap states corresponding to the neighboring cells being occupied by different spin-component particles is allowed when moving along adjacent traps.

After employing the two rules above, we can construct a transfer matrix $\mathbf{T}$ to represent the percolation problem, where the grand partition function is given by $\Xi(z,\mathcal{N})=\text{Tr}\,\mathbf{T}^{\mathcal{N}}$ (see Supplementary Material accompanying this paper). Figure 3 shows the cluster distribution $n(l)$ and the macroscopic magnetic moment $S^{2}/S^{2}_{\text{max}}$ of the Hubbard model with SU($3$), SU($4$) and SU($10$) symmetries for instances, at the particle density $p\equiv n/\mathcal{N}=0.99$. For the cluster distribution $n(l)$ depicted in Fig. 3a, a non-monotonic behavior with a peak structure at some finite $l$ and vanishing value for large enough cluster is observed for all cases. With increasing $N$, the peak becomes higher and narrower, with the central position shifted to a larger value of $l$. This tendency indicates that the effective repulsion induced by the nontrivial weight is more prominent for cases of higher spin symmetry. From Fig. 3b, we find that the macroscopic magnetic moment scales to zero for all cases, implying that the ground state is still paramagnetic when the flat band is nearly $1/N$ filled. Besides, the macroscopic magnetic moment for SU$(10)$ is smaller than that of SU$(3)$ at the same concentration. This observation suggests that the formation of large enough clusters is more difficult in SU$(10)$ system than the SU$(3)$ case, as can also be inferred from the tails on the right side of Fig. 3a. We will find these results consistent with following analytic discussions.

In the thermodynamic limit ($\mathcal{N}\to\infty$ with $n/\mathcal{N}$ finite), we obtain the cluster distribution

[TABLE]

where $z=e^{\mu}$ is the fugacity, $f_{\text{SU}(N)}(z)$ is a function depends only on $z$ and $N$, and $\lambda_{\text{max}}(z)$ is the eigenvalue with maximal modulus of the $\mathbf{T}$-matrix. The exact expression of $f_{\text{SU}(N)}(z)$ is in general complicated, but fortunately we can treat it as a coefficient function such that some important properties of the system, including the shape of the cluster distribution $n(l)$, can be obtained qualitatively. From Eq. (7), one can easily conclude that as long as ${z/\lambda_{\text{max}}(z)}<1$, clusters of large enough size must vanish. The percolation transition thus occurs when the condition ${z_{c}/\lambda_{\text{max}}(z_{c})}=1$ is satisfied, which results in $z_{c}\to+\infty$ (see Supplementary Material accompanying this paper).

For the macroscopic magnetic moment, we have

[TABLE]

which also proves the ground state is paramagnetic for finite $z$. Further calculation shows that in the thermodynamic limit, the particle density $p$ indeed approaches unity as $z\to+\infty$. Hence the ground state is paramagnetic at any particle density $p<1$ for arbitrary $N$, with an onset of ferromagnetism at $p=1$. These results are consistent with our rigorous proof of flat band ferromagnetism for SU$(N)$ Hubbard model on the Tasaki lattice at filling $p=1$ (see Supplementary Material accompanying this paper). One should remind that the conclusions stated above does not contradict with the Lieb-Mattis theorem Lieb-Mattis , which proves the absence of ferromagnetism of 1D Hubbard model on lattices without loops. However, here we consider a special type of lattices with loops composed by nearest $t_{1}$ and next-nearest neighbor hoppings $t_{2}$ as depicted in Fig. 1. In fact, as the existence of loop is crucial to validate the construction of the many-body ground states Tasaki-review ; Nie2013 ; Nie2018 , the Lieb-Mattis theorem does not apply in this circumstance.

4.**Para-ferro magnetic transition in two dimensions

**

For 2D Tasaki lattice, the transfer matrix formalism employed above can not be directly applied 2DTMM . For instead, we make use of the Metropolis Monte Carlo method with importance samplings for the nontrivial weights to simulate the percolation problem on $2$D $L\times L$ Tasaki lattices with periodic boundary conditions. A new configuration is accepted with the Metropolis probability min$[1,W(q^{\prime})/W(q)]$, in which $q$ denotes the old configuration and $q^{\prime}$ denotes a new configuration. Specifically, to obtain the density of particles versus fugacity $z=e^{\mu}$, we perform grand-canonical simulations with exchange Monte Carlo. We randomly choose a site and change its occupying status, i.e., if it is (not) occupied, we (populate) remove one particle.

The percolation transition can be identified as a first order jump between concentrations $p_{-}$ and $p_{+}$ shown in Fig. 4. For fillings below $p_{-}$, the possibility to find a cluster of the same color is vanishingly small in the thermodynamic limit, suggesting a paramagnetic state in the original Hubbard model. When $p>p_{+}$, the clusters of the same color can be linked together to reach a macroscopic size. In the SU($N$) Hubbard model, this regime corresponds to a system with a macroscopic ferromagnetic domain.

In order to determine the para-ferro magnetic transition point, we perform canonical simulations with fixed density of particles for various lattice sizes $L$ to obtain the macroscopic magnetic moments, and extrapolate the results to the thermodynamic limit. Numerically, we generate a new configuration $q^{\prime}$ from the old one $q$ by randomly permuting the occupied status of two sites to guarantee the constraint that the number of particles is fixed. As shown in Figs. 5a, 5d and 5g, the macroscopic magnetic moments scale to zero at $p=0.64$ for the SU$(3)$ Hubbard model, at $p=0.65$ for the SU$(4)$ model, and at $p=0.69$ for the SU$(10)$ case. This implies that the corresponding ground states are paramagnetic below such fillings. In Figs. 5b, 5e and 5h, the particle densities are increased to $p=0.74$, $0.77$, and $0.82$ for the SU$(3)$, SU$(4)$, and SU$(10)$ models, respectively. The macroscopic magnetic moments in all cases scale to a finite values in the thermodynamic limit, indicating the emergence of ferromagnetism. In this regime, the snapshot (see Supplementary Material accompanying this paper) shows the concomitant phase-separated states, where macroscopic ferromagnetic domains can be observed. In Figs. 5c, 5f and 5i, the particle density is further enhanced to $p=0.83$, $0.88$, and $0.95$ for the SU$(3)$, SU$(4)$, and SU$(10)$ cases, respectively. The macroscopic magnetic moments scale to finite values with $S^{2}/S_{\text{max}}^{2}\approx 1$ in the thermodynamic limit, and are independent of the system size. In this regime, the snapshot shows the ferromagnetic phase spans across the entire system.

Based on the numerical results in both grand canonical and canonical ensembles, we obtain the phase diagrams for SU$(N)$ Hubbard models on 2D Tasaki lattice. By increasing particle density $p$, the para-ferro transition as a first order jump is found to lay between concentrations $p_{-}$ and $p_{+}$, whose exact values are both enhanced with $N$. When $p<p_{-}$, the ground state is paramagnetic with zero macroscopic magnetic moment. When $p_{-}<p<p_{+}$, the ground state falls into the phase-separation regime of ferromagnetic domains. As $p_{+}<p<1$, the ground state is of unsaturated ferromagnetism, where the ferromagnetic phase spans across the entire system. The dependence of the exact values of $p_{-}$ and $p_{+}$ on the value of $N$ is rooted from the stronger effective repulsion in systems with larger $N$. Finally, as the flat band is $1/N$ filled with $p=1$, the ground state for arbitrary $N$ is rigorously proved to be the ferromagnetic state (see Supplementary Material accompanying this paper).

It is interesting to consider the limit $N\to\infty$, where a numerical calculation becomes difficult by directly using the method discussed above due to the largely extended space of spin degrees of freedom. This difficulty may be overcome by employing the very large symmetry of the system and mapping to other more tractable models. From our Monte Carlo simulation results and in the view that the effective repulsive interaction that leads to the breakup of the large clusters increases with larger $N$, the para-ferro transition is expected to occur only at $p=1$ when $N$ goes to infinity. The conjecture coincides with the rigorous proof that the ground state at $p=1$ is of saturated flat-band ferromagnetism (see Supplementary Material accompanying this paper). In this limit $N\to\infty$, the ground state is paramagnetic as long as $p<1$ for both $1$D and $2$D cases.

6.**Conclusions

**

We study the Hubbard models with SU$(N)$ spin symmetry and repulsive interaction on a type of decorated cubic lattices in one and two dimensions, where flat bands are present at the bottom of the band structures. When the filling density $p\equiv n/\mathcal{N}\leq 1$ with $n$ the number of particles and $\mathcal{N}$ the single particle flat band degeneracy, we prove that the ground states of the SU$(N)$ Hubbard model can be mapped to the geometrical configuration of an $N$-state Pauli-correlated percolation problem with a nontrivial weight. Based on this exact mapping, we solve for the ground states for 1D Tasaki lattice by a transfer-matrix method, and conclude that the system is paramagnetic at density $p<1$ for arbitrary $N$. For 2D cases, we consider the SU$(N=3,4,10)$ symmetries as typical examples, and use the Metropolis Monte Carlo method to obtain the ground states numerically. The para-ferro magnetic transition point $p_{c}$ is estimated to lay between $p_{-}<p_{c}<p_{+}$ with $p_{-}(N=3)=0.64(1)$, $p_{+}(N=3)=0.83(0)$, $p_{-}(N=4)=0.65(1)$, $p_{+}(N=4)=0.87(7)$, $p_{-}(N=10)=0.69(4)$, $p_{+}(N=10)=0.94(5)$. We also prove rigorously the existence of flat-band ferromagnetism for all values of $N\geq 2$ in arbitrary dimension, provided that the number of particles equals to the degeneracy of the flat band in the single-particle energy spectrum. The present problem reduces to the standard site percolation problem by taking $N=1$.

**Acknowledgements

**

We would like to thank Zheng-Xin Liu, Deping Zhang and Zi-Xiang Li for helpful discussions. R.L. and W.Z. are supported by the National Key R&D Program of China (Grant No. 2018YFA0306501), the National Natural Science Foundation of China (Grants No. 11434011, No. 11522436, and No. 11774425), the Beijing Natural Science Foundation (Grant No. Z180013), and the Research Funds of Renmin University of China (Grants No. 10XNL016 and No. 16XNLQ03). W.N. is supported in part by the National Key R&D Program of China (Grant No. 2017YFB0405700), the National Natural Science Foundation of China (Grant No. 11704267) and start-up funding from Sichuan University (Grant No. 2018SCU12063).

**Author contributions

**

W. Nie and W. Zhang conceived the project. All authors contributed equally to the analytical calculation. R. Liu and W. Nie performed the Metropolis Monte Carlo simulation for the 2D case. W. Nie and W. Zhang wrote the manuscript with inputs from all coauthors.

Appendix A Characterization of the ground states

In this section, we prove that the many-body ground states of the SU$(N)$ repulsive Hubbard model on lattices where the single-particle dispersion acquires flat-band ground states can be characterized by linear combinations of ferromagnetic clusters, provided the filling density $p\equiv n/\mathcal{N}\leq 1$ with $n$ the number of particles and $\mathcal{N}$ the number of flat-band degeneracy. As the single-particle states in the non-dispersive flat band correspond to localized states, the value of $\mathcal{N}$ equals to the number of those localized states, which are also referred as trapping cells. We emphasize that although the Tasaki lattice are considered as a typical example, the proof is general for all flat-band lattices in arbitrary dimensions.

We denote the set of sites of the Tasaki lattices by $\Lambda$ and the undecorated hypercubic lattice by $V$, then $|V|=\mathcal{N}$. We rewrite the Hamiltonian in the main text as

[TABLE]

with $H_{0}=\sum_{i,j\in\Lambda}t_{i,j}\left(c_{i,\sigma}^{{\dagger}}c_{j,\sigma}+\mbox{H.c.}\right)$, $H_{\textrm{int}}=U\sum\limits_{i}\sum\limits_{\sigma\neq\sigma^{\prime}}n_{i,\sigma}n_{i,\sigma^{\prime}}$. When there exists a flat band, the eigenstates of single particle are localized states and can be represented as $\{a_{u,\sigma}^{\dagger}|0\rangle\}$, where $u=\{1,2,\cdots,\mathcal{N}\}$ denotes the trapping cells, $\sigma=\{1,2,\cdots,N\}$ labels the spin component, and $|0\rangle$ represents the vacuum state. Note that the trapping cells and the sites of the undecorated lattice have a one-to-one correspondence, such that $u$ also denotes the undecorated lattice sites. The expression of the creation operator $a_{u,\sigma}^{\dagger}$ reads

[TABLE]

where $\varphi_{i}^{u}$ is defined as

[TABLE]

Here, $\delta$ is the unit vector between the undecorated and decorated sites and $c$ is the coordination number of the corresponding undecorated lattice. The commutation relation between $H_{0}$ and $a_{u,\sigma}^{\dagger}$ can be easily obtained

[TABLE]

with $\varepsilon_{0}$ the single-particle eigenenergy of the flat band.

When the filling density $p<1$, the number of particles is less than that of the trapping cells, such that each trapping cell is occupied by at most one particle to avoid energy penalty of repulsion. Obviously, the ferromagnetic state

[TABLE]

with all $n$ particles aligned at the spin $N$ component is a ground state of $H$ given by Eq. (9). Here, $A$ (with $|A|=n\leq\mathcal{N}$) is an arbitrary subset of the undecorated lattice $V$. It can be verified that $H_{0}|\Phi_{A,N}\rangle=n\varepsilon_{0}|\Phi_{A,N}\rangle$ and $H_{\textrm{int}}|\Phi_{A,N}\rangle=0$.

From the fully-polarized ferromagnetic state $|\Phi_{A,N}\rangle$, we can construct states composed of ferromagnetic clusters

[TABLE]

where $i=\{1,2,\cdots M\}$ labels the linked clusters of occupied sites, $\mathbf{S}^{-}_{k,i}$ is one of the lowering operators of the complete lowering/raising operator set $\{\mathbf{S}^{\pm}_{1},\cdots,\mathbf{S}^{\pm}_{C_{N}^{2}}\}$ for the SU($N$) group acting on the i-th cluster, $m_{i}$ is a non-negative integer among $\{0,1,\cdots,d_{\textrm{SU}(N)}(|C_{i}|)-1\}$ labeling the spin degeneracy, and $C_{N}^{2}$ is the binomial coefficient. The action of operators $[\mathbf{O}]^{\prime}$ is arranged in such a way that $\mathbf{O}|\Phi_{A,N}\rangle\neq 0$. It can be verified that by selecting all possible combinations of $m_{i}$ and $\{\mathbf{S}^{\pm}_{k,i}\}$, we can obtain a number of $d_{\textrm{SU}(N)}(|C_{i}|)$ linearly independent degenerate states of the i-th cluster. Notice that these $d_{\textrm{SU}(N)}(|C_{i}|)$ degenerate ferromagnetic states can be represented by the Young tableau representation $:\bm{\sigma}_{i}:=D\sum_{P}\prod_{u_{i}\in C_{i}}a_{u_{i},\sigma_{u_{i}}}^{\dagger}|0\rangle$, where $D$ is the normalized coefficient, $\bm{\sigma}_{i}=(\sigma_{u_{1}},\cdots,\sigma_{u_{|C_{i}|}})$ is the sector of spin configuration of the $i$-th cluster. The arrangement of spin components in $:\quad:$ is non-decreasing from left to right, and the summation over permutation $P$ takes account of all possible non-repeated spin configurations $(\sigma_{u_{1}},\cdots,\sigma_{u_{|C_{i}|}})$ of the i-th cluster. Because of the commutation relation $[\mathbf{S}^{-}_{k,i},H]=0$, one can easily conclude that the states $|\Phi_{A,\{m_{i}\}}\rangle$ expressed in Eq. (14) are also ground states of $H$. In the following, we will demonstrate that an arbitrary ground state of the SU$(N)$ Hubbard mode can be expanded by states defined by Eq. (14) with various $A$, $\{m_{i}\}$ and $\{\mathbf{S}_{k,i}^{-}\}$.

An arbitrary ground state $|\Phi\rangle_{\textrm{GS}}$ of the SU$(N)$ Hubbard model can be represented with the aid of trapping cell states

[TABLE]

where $A_{1},\cdots,A_{N}$ are subsets of $V$ satisfying the number constraints $|A_{1}|+|A_{2}|+\cdots+|A_{N}|=n$. For any given subsets $\{A_{1},\cdots,A_{N}\}$, we can decompose them into sets of connected clusters $C_{k}$, such that $\bigcup\limits_{\sigma=1}^{N}A_{\sigma}=\bigcup\limits_{k=1}^{M}C_{k}$. By using the commutation relation $[H_{0},a_{u,\sigma}^{\dagger}]=\varepsilon_{0}a_{u,\sigma}^{\dagger}$, we obtain $H_{0}|\Phi\rangle_{\textrm{GS}}=n\varepsilon_{0}|\Phi\rangle_{\textrm{GS}}$.

The ground state $|\Phi\rangle_{\textrm{GS}}$ must satisfy another condition $H_{\textrm{int}}|\Phi\rangle_{\textrm{GS}}=0$, which is equivalent to $c_{i,\sigma}c_{i,\sigma^{\prime}}|\Phi\rangle_{\textrm{GS}}=0$. The fermionic site operators $c_{i,\sigma}^{\dagger}$ and $c_{i,\sigma}$ can be expressed by the trapping cell operators $a_{u,\sigma}^{\dagger}$ and its corresponding adjoint. By introducing the Gramm matrix $G$

[TABLE]

and defining

[TABLE]

we have the completeness relation

[TABLE]

where

[TABLE]

Thus, we can define the operators $b_{u,\sigma}$ adjoint to $a_{u,\sigma}^{\dagger}$

[TABLE]

and the auxiliary operators

[TABLE]

which satisfy the anti-commutation relations

[TABLE]

With the aid of the completeness relations (18) and (19), we obtain the expansion formulae

[TABLE]

for arbitrary $i\in\Lambda$ and $\sigma=1,2,\cdots,N$. By substituting Eq. (25) into the ground state condition $c_{i,\sigma}c_{i,\sigma^{\prime}}|\Phi\rangle_{\textrm{GS}}=0$, and employ the quasi-locality and local connectivity conditions 1 ; 2 ; 3 ; 4 satisfied by the Tasaki lattices, we can finally obtain two constraints for the coefficients $f(A_{1},\cdots A_{N})$:

(i) One must have

[TABLE]

if $A_{\sigma}\cap A_{\sigma^{\prime}}\neq\emptyset$ for any $\sigma\neq\sigma^{\prime}$. This indicates that in each particle distribution $\{A_{1},\cdots,A_{N}\}$ in the ground states, the subsets $A_{\sigma}$ and $A_{\sigma^{\prime}}$ do not overlap. In other words, a lattice site is occupied by at most one particle.

(ii) For any two subsets $\{A_{1},\cdots,A_{N}\}$ and $\{A_{1}^{\prime},\cdots,A_{N}^{\prime}\}$ which can be decomposed into the same set of clusters $\bigcup\limits_{\sigma=1}^{N}A_{\sigma}=\bigcup\limits_{\sigma^{\prime}=1}^{N}A_{\sigma^{\prime}}^{\prime}=\bigcup\limits_{k=1}^{M}C_{k}$ and satisfy $A_{\sigma}\cap A_{\sigma^{\prime}}=A_{\sigma}^{\prime}\cap A_{\sigma^{\prime}}^{\prime}=\emptyset$ for any $\sigma\neq\sigma^{\prime}$, if the number of particles of the $\sigma$-th spin component are the same within each cluster, i.e.,

[TABLE]

for all $\sigma$ and $k$, the coefficients associated with the two subsets are equivalent

[TABLE]

Thus, the coefficients $f(A_{1},\cdots,A_{N})$ are symmetric under permutation of particles within an individual cluster. The summation of all possible permutation states then resembles the ferromagnetic state obtained from the fully-polarized cluster after a sequence of spin lowering operator.

The constraints (i) and (ii) prove that any ground state $|\Phi\rangle_{\textrm{GS}}$ is a linear combination of the states given by Eq. (14). This characterization of ground states is important when mapping to a percolation problem as well as using the transfer-matrix method.

Before concluding this section, we emphasize that for the special case with $n=\mathcal{N}$, the flat-band is $1/N$ filled and there can only exist one cluster with size equals to the number of trapping cells. Thus, the states by linear combination of Eq. (14) are fully ferromagnetic states with a spin degeneracy associated with the SU$(N)$ symmetry, which proves the flat-band ferromagnetism of the ground state (When we are submitting this paper, we are informed that flat band ferromagnetism of SU($N$) Hubbard model at filling $p=1$ was mentioned in the talk of H. Katsura: https://pcs.ibs.re.kr/PCS_Workshops/PCS_Flatband_Networks_in_Condensed_Mater_and_Photonics.html). Notice that this statement is valid for arbitrary $N$ in any dimensions.

Appendix B Transfer-matrix method for 1D Tasaki lattice

As discussed in the main text, the populating rules of $N$-component particles over $\mathcal{N}$ traps can be expressed by a transfer matrix $\mathbf{T}$, which reads

[TABLE]

with $z=e^{\mu}$. The matrix element $T(\sigma_{i},\sigma_{i+1})$ takes the value $1$ ($z$) if the site $i$ is empty (occupied) for spin configurations satisfying the two rules defined in the main text. The contribution of all possible configurations to the grand-canonical partition function can be calculated as

[TABLE]

To calculate some physical quantities such as the average occupation number of a given site $\langle n_{i}\rangle$, we introduce the $\mathbf{N}$-matrix,

[TABLE]

and we have

[TABLE]

Besides, we introduce the $\mathbf{S},\mathbf{C},\mathbf{F}$-matrix,

[TABLE]

which means a sequence starting with an empty site, both neighboring traps in a sequence being occupied, and the last site of a sequence being followed by an empty one respectively. Note that the $\mathbf{C}$-matrix satisfies

[TABLE]

For convenience, we diagonalize the $\mathbf{T}$-matrix with $\mathbf{\Lambda}=\mathbf{X}^{-1}\mathbf{T}\mathbf{X}$, where the diagonal matrix reads

[TABLE]

with $\lambda_{1},\lambda_{2},\cdots,\lambda_{N}$ the non-zero and non-repeated roots of the following equation

[TABLE]

and the transformation matrix

[TABLE]

Before proceeding the calculation, it is necessary to analyze the properties of these $N$ roots. It can be verified that if $\lambda_{i}$ is a root, then $\lambda_{i}^{*}$ is also a root. In other words, there exits a pair of complex conjugate roots. Besides, we point out that the root with the maximal modulus is always real as demonstrated in Fig. 6.

The cluster distribution $n(l)$ can be obtained as

[TABLE]

We can find that the elements of $\mathbf{S}\mathbf{C}^{l-1}\mathbf{F}$ are all zeros except

[TABLE]

Thus, in the thermodynamic limit $\mathcal{N}\to\infty$ we obtain

[TABLE]

where $f_{\text{SU}(N)}(z)=\frac{(\mathbf{X})_{1,\text{max}}(\mathbf{X}^{-1})_{\text{max},1}}{\lambda_{\text{max}}(z)}$, $\lambda_{\text{max}}(z)$ is the root with the maximal modulus of the $N$ roots, and $\text{max}(\neq 1)$ is the index which satisfies $(\mathbf{\Lambda})_{\text{max},\text{max}}=\lambda_{\text{max}}$. We emphasize again that $\lambda_{\text{max}}$ is real and if not, we will obtain a cluster distribution related to $\cos[\mathcal{N}\text{arg}(\lambda_{\cdots})]$, which will oscillate in the limit of $\mathcal{N}\to\infty$ rather than converging to a fixed value.

The percolation transition occurs when the condition ${z_{c}/\lambda_{\text{max}}(z_{c})}=1$ is satisfied. To determine the transition point $z_{c}$, we rewrite the roots of Eq. (78) as $\lambda=\alpha(z)z$ with complex $\alpha(z)$, such that the percolation transition takes place when $\alpha=1$. By substituting the representation of $\lambda$ back into Eq. (78), we obtain

[TABLE]

with $i=1,2,\cdots,N$. From the expression above, one can easily find that $z\to+\infty$ when $\alpha\to 1$. Indeed, it can be verified that there exists at least one real $\alpha(z)$ satisfying Eqs. (78) and (87) for arbitrary $z>0$, and the absolute value of $\alpha$ is no less than 1 (only if $z\to+\infty$). Hence the percolation transition occurs at $z_{c}\to+\infty$. We further point out that at the critical point $z_{c}\to+\infty$, the real component of $\alpha_{i}$ approaches unity while the imaginary part goes to zero.

In order to determine the results in terms of the average occupation number $p$, in the thermodynamic limit we have,

[TABLE]

We further let $z\to+\infty$ and obtain with the aid of Eq. (87) that $\lim\limits_{z\to\infty}(\lambda_{i}-z)/{z}=\lim\limits_{z\to\infty}(\alpha_{i}(z)-1)\to 0$ for $i=1,2,\cdots,N$. Thus we have $(\mathbf{X})_{1,i}=-1\delta_{i,1}$ and $(\mathbf{X}^{-1})_{i,1}=-1\delta_{i,1}$ for $i=1,2,\cdots,N+1$, which leads to $(\mathbf{X}\mathbf{\Lambda}^{\mathcal{N}}\mathbf{X}^{-1})_{11}/(\textrm{Tr}\,\mathbf{X}\mathbf{\Lambda}^{\mathcal{N}}\mathbf{X}^{-1})\to 0$ and $p\to 1$.

Appendix C Supplementary material for 2D Monte-Carlo simulation

To show the features of the ground state at different particle concentrations in canonical simulations, we plot as a typical example the snapshots of Pauli-correlated percolation configurations for SU($3$) Hubbard model in Fig. 7. When particle density $p_{1}=0.64<p_{c}$ as shown in Fig. 7a, the ground state is paramagnetic with no macroscopic clusters. When particle density exceeds $p_{c}$ as illustrated in Fig. 7b, the ground state becomes ferromagnetic with macroscopic domains which are phase-separated from each other. If the filling factor is further increased, ferromagnetic phase spans the entire system as depicted in Fig. 7c.

Appendix D Connection to the standard percolation problem in the limit of $N=1$

From Eq.(86) of the Supplementary Material, we obtain the cluster distribution for 1D Tasaki lattice in the limit of $N=1$,

[TABLE]

With the aid of Eq. (88) of the Supplementary Material, we can further get the average occupation number $p(z)$ in the thermodynamic limit,

[TABLE]

The result above determines the fugracity $z=\frac{p}{1-p}$. By substituting this expression into Eq. (89) in the Supplementary Material, we finally obtain

[TABLE]

which is consistent with the standard percolation 5 .

For 2D case, according to Eq. (2) in the main text, the spin degeneracy $d_{\text{SU}(N)}=1$ in the limit of $N=1$. Thus, we recover the standard percolation problem.