# Exact Diagonalization of SU(N) Fermi-Hubbard Models

Thomas Botzung, Pierre Nataf

► To cite this version:

Thomas Botzung, Pierre Nataf. Exact Diagonalization of SU(N) Fermi-Hubbard Models. 2023. hal-04292727

HAL Id: hal-04292727

https://hal.science/hal-04292727v1

Preprint submitted on 17 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Exact Diagonalization of SU(N) Fermi-Hubbard Models

Thomas Botzung¹ and Pierre Nataf¹

¹Laboratoire de Physique et Modélisation des Milieux Condensés, Université Grenoble Alpes and CNRS, 25 avenue des Martyrs, 38042 Grenoble, France

(Dated: September 19, 2023)

We show how to perform exact diagonalizations of SU(N) Fermi-Hubbard models on L-sites clusters separately in each irreducible representation (irrep) of SU(N). Using the representation theory of the unitary group U(L), we demonstrate that a convenient orthonormal basis, on which matrix elements of the Hamiltonian are very simple, is given by the set of semi-standard Young tableaux (or equivalently the Gel'fand-Tsetlin patterns) corresponding to the targeted irrep. As an application, we study the robustness of some SU(N) phases predicted in the Heisenberg limit upon decreasing the on-site interaction on various lattices of size L ≤ 12 and for 2 ≤ N ≤ 6. In particular, we show that a long range color ordered phase emerges for moderate U for N=4 at filling 1/4 on the triangular lattice.

The Fermi-Hubbard model (FHM) is among the most important models in condensed matter [1–3]. In particular, the SU(2) FHM on the square lattice might describe the physics of electrons carrying a spin one half in cuprates superconductors [3–5] and has motivated numerous theoretical investigations [6, 7]. Considering N, the number of degenerate orbitals, as an integer parameter of the models, a natural extension of the SU(N = 2) FHM is the SU(N) FHM [8–10].

This higher symmetry group Hamiltonian, firstly introduced as a theoretical tool to provide an asymptotic description of spins 1/2 in the large N limit [11–14], can also describe some condensed matter systems like transition metal compounds [15, 16] or graphene with SU(4) spin valley symmetry or in twisted bilayer [17]. Alternatively, Alkaline-earth cold atoms like $^{173}$Yb or $^{87}$Sr can simulate SU(N) invariant FHMs for N up to 10 on various engineered optical lattices [18, 19]. Besides, the continuous experimental achievements in this field [20–25] brought theoreticians to investigate these systems in order to look for exotic phases that would generalize their N = 2 counterpart.

Most of the theoretical investigations focused on the large on-site repulsion limit, where the atoms, in the Mott insulating phase, are described by SU(N) Heisenberg models (HM) [26–34]. Depending on the lattices and on the number of colors N, different two-dimensional phases are predicted at T = 0 among which the SU(N) plaquette phases [35–40] are cousins of the valence band states for spins 1/2, the Neel long-range color ordered (LRO) states [30, 31] are analogous to the famous (π, π) (resp. 120°) Neel states existing on the square [41, 42] (resp. the triangular [43, 44]) lattice for N = 2, and diverse kinds of SU(N) spin liquids [45–53] generalizing the Anderson Resonating Valence Band (RVB) states [4, 54]. A part from the one-dimensional system where there is a Bethe ansatz solution [55], the theoretical investigation of these models, based on advanced numerical tools is challenging mainly because the dimension of the full Hilbert space on finite-size lattices increases exponentially, being equal to N^L where L is the number of sites of the cluster, for filling 1/N (exactly one particle per site) in the

SU(N) HM.

However, it was realized that working in the SU(N) singlets subspace, which contains the ground state (GS) in the antiferromagnetic case, is very advantageous as its dimension is much smaller than N^L. For instance for N = 6 and L = 12, such a dimension is equal to 132 while N^L ≡ 6^12 ≈ 2 × 10^9. In addition, the exact diagonalization (ED) of the HM directly in the SU(N) singlets subspace can be made easy on the basis of Standard Young Tableaux (SYT) using the orthogonal representation of the group of permutations S_L [56, 57]. It is crucial to extend this theory to the SU(N) FHM as the dimension of the full Hilbert space is even larger, i.e equal to 2^NL. In fact, for N = 6, L = 12 at filling 1/6 the dimension of the singlets subspace for the SU(6) FHM is ≈ 14 × 10^6 while the full Hilbert space has dimension 2^72 ≈ 5 × 10^21.

In this letter, we use the representation theory of the Lie group U(L) to show how to perform ED of the SU(N) FHM directly in each irrep. After the description of the method, we apply the procedure to show some ED results on square/triangular clusters of size up to L = 12 and for N up to N = 6 to see how robust are some SU(N) Mott insulating phases while decreasing the on-site repulsion.

The Hamiltonian for the SU(N) FHM reads:

$$H = \sum_{\langle i,j \rangle} \left( -t_{ij}E_{ij} + \text{h.c.} \right) + \frac{U}{2} \sum_{i=1}^{L} E_{ii}^2, \quad (1)$$

where the t_{ij} are the (possibly complex) hopping amplitude between sites i and j of a L-sites finite cluster, the on-site interaction amplitude is U. The SU(N) invariant hopping terms E_{ij} = E_{ji}^\dagger = \sum_{\sigma=1}^{N} c_{i\sigma}^\dagger c_{j\sigma} satisfy the commutation relation of the U(L) generators (∀1 ≤ i, j, k, l ≤ L):

$$[E_{ij}, E_{kl}] = \delta_{jk}E_{il} - \delta_{li}E_{kj}, \quad (2)$$

so that the Hamiltonian in Eq. (1), where the integer parameter N is hidden, can be seen as an element of the Lie algebra of the unitary group U(L) [58]. It should be considered as the counterpart of the quantum permutation Hamiltonian for the SU(N)-invariant HM, i.e
$H = \sum_{\langle i,j \rangle} J_{ij} P_{ij} + \text{h.c.}$, with $P_{ij}$ (resp. $J_{ij}$) the permutation (resp. coupling constant) between interacting sites i and j, for which the representation theory of the algebra of the group of permutations was used to perform ED directly and separately in each irrep of SU(N) [56].

We remind that an irrep of SU(N) is labelled by a Young tableau (YT), or shape α (See Fig. 1) the N rows of which representing N integers α = [α₁, α₂, ..., αN] such that α₁ ≥ α₂ ≥ ... ≥ αN ≥ 0 and ΣN i=1 αᵢ = M, where M is the number of particles. For L = 1 site, the M-fermions wave-functions are fully antisymmetric in the exchange of colors, and belong to the irrep α which is a M-boxes single column α = [1, 1, .., αM = 1], and its dimension $D_{L=1}^{M,N} = \binom{N}{M}$ simply corresponds to all the ways to choose M colors among N.

More generally, for L sites and M SU(N) fermions, the Hilbert space dimension $D_L^{M,N}$ is equal to:

$$ \sum_{\alpha} h_{L}^{\bar{\alpha}} \prod_{i=1}^{L} \binom{N}{\bar{\alpha}(i)} = \sum_{\alpha} d_{N}^{\alpha} d_{L}^{\bar{\alpha}}, \quad (3) $$

where the sums on both side run over all the YT α of maximum L columns and N rows. $\bar{\alpha} = [\bar{\alpha}(1), ..., \bar{\alpha}(L)]$ is defined as the conjugate YT of α, transforming rows into columns (i.e flipping the shape by 45°), cf Fig. 1 for some examples. In the LHS of Eq. (3), $\bar{\alpha}$ is a distribution of fermions: $\bar{\alpha}(j)$ being the number of fermions (necessarily $\le N$) on site j for $1 \le j \le L$, $\prod_{i=1}^{L} \binom{N}{\bar{\alpha}(i)}$ is the number of states for such a distribution. The factor $h_L^{\bar{\alpha}}$, defined as $h_L^{\bar{\alpha}} = L! / \prod_{k=0}^{N} (n_k^{\bar{\alpha}})!$, where $n_k^{\bar{\alpha}} = \text{Cardinal}\{j \in [1; L] / \bar{\alpha}(j) = k\}$, is the number of distributions corresponding to a given partition $\bar{\alpha}$ while permuting the $\bar{\alpha}(j)$ (or the site indices j) for $1 \le j \le L$.

In the RHS of Eq. (3), $d_N^\alpha$ (resp. $d_L^{\bar{\alpha}}$) stands for the dimension of the SU(N) irrep α (resp. the U(L) irrep $\bar{\alpha}$)[59], that we can calculate using existing formulas, like the hook length formulas [60, 61]. These dimensions are equal to the number of semi-standard Young tableaux (ssYT) of shape α (resp. $\bar{\alpha}$) filled with numbers from 1 to N (resp. L), since these latter form a basis of the SU(N) or U(L) irrep. Given a U(L) irrep represented by some YT, a ssYT is filled up with integers numbers from 1 to L in non-descending order from left to right in any row (repetitions allowed), and in strictly ascending order (repetitions not allowed) from top to bottom in any column, cf Fig. 1 and Eq. (6) for some examples.

As we detail below, for the color-invariant SU(N) FHM model in Eq. (1), the Hilbert space for the set of M SU(N) fermions on a L-sites cluster can be decomposed, or color-factorized, following the RHS of the equation for the dimension $D_L^{M,N}$, i.e Eq. (3). In particular, targeting a given collective SU(N) irrep α, we will need to diagonalize a matrix of dimension $d_L^\alpha$, and there will be $d_N^\alpha$ independent copies (some multiplicity) of the corresponding spectrum in the full energy spectrum of the model. For instance, when M is a multiple of N, one important sector is the SU(N) singlets sector as it usually

<figure><img src="image_2.png" /><figcaption>Figure 1. Examples of SU(N = 2) irrep α for a L = 4-sites cluster and M = 6 fermions in a) and M = 4 fermions in b). The Young tableaux (YT) are α = [33] (resp. α = [31]), representing singlets (resp. spin 1) states. We associate the conjugate YT, flipping the shape by 45°, α = [33] → $\bar{\alpha}$ = [222] for a) and α = [31] → $\bar{\alpha}$ = [211] for b). We fill up $\bar{\alpha}$ to get the Highest Weight State (|hws⟩) to which we associate $d_N^{\bar{\alpha}}$ physical states spanning the SU(N = 2) irrep α ($d_{N=2}^{\bar{\alpha}}$ different values of $Q_z$) and which generate $d_L^{\bar{\alpha}}$-dimensional and independent sectors under the application of the operators $E_{ij}$. See text for details.</figcaption></figure>

contains the lowest energy eigenstates (in the large U > 0 limit for instance); it is labelled by the perfectly rectangular N-rows YT α = $\alpha_{S,M} \equiv [M/N, M/N, ..., M/N]$. In this case, $d_N^{\alpha_{S,M}} = 1$ and for L = 12 at filling 1/N ($M = L$), one has for instance $d_{L=12}^{\alpha_{S,M=L}} = 13026013$ for N = 4 and $d_{L=12}^{\alpha_{S,M=L}} = 14158144$ for N = 6. This should be compared to the dimensions of the sector usually addressed in standard ED with a fixed number of fermions of each color (conserving the U(1) symmetry) which is $d_{L,N}^{U(1)} = \binom{L}{N}^N$: one has $d_{L=12,N=4}^{U(1)} \approx 2.34 \times 10^9$ (resp. $d_{L=12,N=6}^{U(1)} \approx 8.27 \times 10^{10}$). As N increases, it is more and more advantageous to implement the full SU(N) symmetry, working in the SU(N) singlets sector, and more generally in a sector of a given irrep α.

For a given U(L) irrep $\bar{\alpha}$, the Highest Weight State (|hws⟩), uniquely (up to some similarity) and fully determines the irrep as one can generate the entire basis by applications of the generators $E_{ij}$ (for $1 \le i, j \le L$). It is represented by the shape $\bar{\alpha}$ filled with 1 for the first row (of length $\bar{\alpha}(1)$), 2 for the second row of length $\bar{\alpha}(2)$), etc...(cf Fig. 1). It is defined by the following properties: $E_{ii}|hws\rangle = \bar{\alpha}(i) \forall i \in [1; L]$ and $E_{ij}|hws\rangle = 0$ for $i < j$. Crucially, in the $D_L^{M,N}$-dimensional space of M-fermions SU(N) wave-functions, there are $d_N^\alpha$ orthonormal states $|\phi_{\alpha,k}^{\text{hws}}\rangle$ ($k = 1 \cdots d_N^\alpha$) which have these properties and can then be represented by the same ssYT associated with the |hws⟩. For example, for the SU(N) singlets irrep $\alpha_{S,M}$, there is only one state $|\phi_{\alpha,1}^{\text{hws}}\rangle$ and it is the product of SU(N) singlets for sites 1, 2, ..., M/N, with no particles on sites M/N + 1, ..., L (cf Fig. 1).

Introducing the SU(N) total generators $Q_a = \sum_{i=1}^{L} \sum_{\sigma, \sigma'} c_{i, \sigma'}^{\dagger} \lambda_{\sigma'}^{\sigma} c_{i, \sigma}$ for $a = 1 \cdots N^2 - 1$, with
($\lambda_a^{\sigma'\sigma}$)$_{\sigma'\sigma}$ the generalized Gell-Mann matrices for SU(N), the states $|\phi_{\alpha,k}^{\text{hws}}\rangle$ are eigenvectors of the $N-1$ diagonal $Q_a$ (spanning the Cartan subalgebra) with different set of $N-1$ eigenvalues. Then, the $d_N^\alpha$ sectors generated by the (repeated) applications of the $E_{ij}$ (for $1 \le i, j \le L$) on each $|\phi_{\alpha,k}^{\text{hws}}\rangle$ (for $k \in [1; d_N^\alpha]$) will be invariant under the action of H in Eq. (1) since $[Q_a, E_{ij}] = 0, \forall 1 \le i, j \le L$ and $a \in [1; N^2 - 1]$. Moreover, each sector will independently represent the U(L) irrep $\bar{\alpha}$, i.e being in one to one correspondence with the set of ssYT of shape $\bar{\alpha}$. On such a basis, the matrix elements of the infinitesimal generators $E_{p,p}$, which are the occupation numbers on site $p$ for $p = 1 \cdots L$ and of $E_{p-1,p}$ ($p = 2 \cdots L$), which generalize the lowering operator $J_-$ for SU(2), are very simple. They were found by Gelfand and Tsetlin [62, 63] and expressed in the basis of what is now called the "Gelfand-Tsetlin" patterns (GTp) which are triangular arrangements of integers fully equivalent to the ssYT, the first row of which being the lengths $\bar{\alpha}(1)\bar{\alpha}(2)\cdots\bar{\alpha}(L)$ (See. Fig. 1 for some examples and [64] for the systematic conversion from GTp to ssYT for a given shape $\bar{\alpha}$). Firstly, for $p = 1 \cdots L$, the operators $E_{p,p}$ are diagonal in the basis of the ssYT, with coefficients equal to the number of occurrences of $p$ inside each ssYT. It also means that we know readily from its content, the occupation numbers of a ssYT on each site (cf Fig. 1).

Secondly, calling $|\nu\rangle$ a ssYT, one has for $p = 2 \cdots L$:

$$E_{p-1p}|\nu\rangle = \sum_{j=1}^{p-1} a_{p-1}^j F_{p-1}^j |\nu\rangle, \quad (4)$$

where the tableau operators $F_{p-1}^j$ transform the number $p$ in the $j^{th}$ row in $|\nu\rangle$ into $p-1$. As for the coefficients $a_{p-1}^j$, which vanishe in case such a transformation is not possible either because there is no $p$ in the $j^{th}$ row of $|\nu\rangle$, either because the resulting tableau is not a proper ssYT, they read:

$$a_{p-1}^j = \left| \frac{\prod_{i=1}^{p} (l_{i,p} - l_{j,p-1}) \prod_{i=1}^{p-2} (l_{i,p-2} - l_{j,p-1} - 1)}{\prod_{i \neq j} (l_{i,p-1} - l_{j,p-1}) \prod_{i \neq j} (l_{i,p-1} - l_{j,p-1} - 1)} \right|^{1/2}, \quad (5)$$

where $l_{k,q} = m_{k,q} - k$ with $m_{k,q}$ the length of the $k^{th}$ row of the sub-tableau that remains when we delete all the boxes containing numbers $> q$ in $|\nu\rangle$. Thus, for the SU(4) adjoint irrep at filling 1/4 for L=12 (the basis has then 57972915 elements), we have for instance:

$$E_{23} \begin{bmatrix} 1 & 1 & 2 & 3 \\ 2 & 3 & 3 & 4 \\ 4 & 5 & 6 & 5 \\ 5 & 5 & 6 & 5 \end{bmatrix} = \sqrt{\frac{5}{6}} \begin{bmatrix} 1 & 1 & 2 & 2 \\ 2 & 3 & 3 & 4 \\ 4 & 5 & 6 & 5 \\ 5 & 5 & 6 & 5 \end{bmatrix} + \sqrt{\frac{16}{6}} \begin{bmatrix} 1 & 1 & 2 & 3 \\ 2 & 2 & 3 & 4 \\ 4 & 5 & 6 & 5 \\ 5 & 5 & 6 & 5 \end{bmatrix}. \quad (6)$$

Moreover, from successive applications of the commutation relations Eq. (2), and from $E_{ij} = E_{ji}^\dagger$, one gets the matrix representing H defined in Eq. 1 in the irrep $\bar{\alpha}$, which corresponds to the SU(N) irrep $\alpha$.

Finally, let's mention the large U limit at filling 1/N. Then, the ssYT with strictly more than one occurrence of a site index are too high in energy, and we are left with the subset of SYTs, in agreement with the ED protocole developed for the HM [56]. Then, for each SYT $|\mu\rangle$ and for fixed $i = 2 \cdots L$, $a_{i-1}^j$ is either 0, either $a_{i-1}^j = \sqrt{|d-1|/d}$ for $j$ the unique row where $p$ stands in $|\mu\rangle$, with $d$ the Manhattan distance between $p$ and $p-1$. Then, the transposition $P_{ii-1}$ being equal to $-E_{ii-1}E_{i-1i} + 1$, we refind the Young rules for the orthogonal representations of the algebra of the group of permutations on the basis of SYT in a given irrep $\alpha$ [56, 65].

<figure><img src="image_3.png" /><figcaption>Figure 2. Correlations patterns of the ground states of the FHM at filling 1/N for U = 10 in a) c) and d) and for U = 12 in b). They are defined as $\langle P_{ij} \rangle - 1/N$ where $P_{ij} = -1 + E_{ij}E_{ji}$, with the reference site 1 in black, and j the site indices being blue (resp. red) for positive (resp. negative) correlation, with area proportional to its absolute value. At the top, triangular lattice with L = 12 sites, for SU(3) (in a)) compatible with the three-sublattice Neel order [31, 66, 67], and for SU(4) (in b)) compatible with the four-sublattice Neel order [28]. c) SU(5) on the cluster $\sqrt{10} \times \sqrt{10}$, pattern compatible with the (chess) knight move LRO [56]; d) SU(6) on the 3 × 4 cluster, compatible with the SU(6) plaquette state [39].</figcaption></figure>

We have applied this theory to study the FHM of Eq. (1) for uniform nearest neighbors hopping $t_{ij} \equiv t = 1$ as a function of U at filling 1/N, starting from the Heisenberg limit ($U \to \infty$) and diminishing U. The ground states of the SU(3) HM on the triangular lattice (TL) and of the SU(5) HM on the square lattice (SL) are both Neel LRO states, with a three sublattice ordering pattern for SU(3) [31, 66, 67] and a (chess) knight move pattern for SU(5) [56]. We give evidence of such orders by calculating the simple correlations patterns of the exact ground states of the FHM in Fig. 2 for U = 10. Moreover, the energy spectra plotted as a function of the quadratic Casimir $C_2$ [56, 68] of the different irreps $\alpha$
<figure><img src="image_4.png" /><figcaption>Figure 3. Energy spectra of the FHM in Eq. 1 (with $t_{ij} \equiv 1$ at filling $1/N$) as a function of the quadratic casimir $C_2$ for various values of $U$ and $N$ for the $L = 12$ periodic triangular lattice (in a,b,c,d) and for the $L = 4 \times 3$ periodic square lattice in (e) where we focused on $N = 6$. Note that the constant $LU/2$ has been withdrawn. In (c), the Atos is reminiscent of the one revealing the four sublattice order in the HM with nnn couplings [28]. Inset of (e): spin (resp. singlet) gap in blue (resp. red).</figcaption></figure>

exhibit an Anderson tower of states (Atos) which reveals the continuous symmetry breaking of SU(3) (resp SU(5)) as shown in Fig. 3 (resp. in [69]). We have checked the convergence in the limit $U \to \infty$ within each irrep $\alpha$ of the eigenergies towards those of the HMs with the factor $2/U$.

While diminishing $U$, the structure of the energy spectra stays the same up to $U \sim 2.5$ (resp. $U \sim 1$) for SU(3) on the $L = 12$ TL (resp. SU(5) on the $L = 10$ SL). Then, some energy plateau as a function of $C_2$ appears for smaller $U$, which is also true for SU(2) on the TL, as shown in Fig. 3. Such a system should be in the metallic phase for $U \lesssim 8.5$, as expected from DMRG simulations on large cylinders [70], so that the plateau could be a signature of the metallic phase in the weak coupling limit. To further characterize the metallic phase and to locate its boundary, we show in Fig. 4 and in [69] the charge gap defined by: $\Delta_c = E_0(M = L + 1) + E_0(M = L - 1) - 2E_0(M = L)$, where $E_0(M)$ is the minimal energy for the lattice with M fermions, which implies the diagonalization over all the relevant $M = L, L \pm 1$ boxes irreps $\tilde{\alpha}$. It suggests that the metallic phase develops for $U \le U_c = 9.8(\pm 0.4)$ for SU(3) on a TL and for $U \le U_c = 8.75(\pm 0.15)$ for SU(5) on a SL, with apparently no intermediate phase between the latter and the LRO in the large $U$ limit.

The scenario of successive LRO phases, with different antiferromagnetic orders as what has been pointed out in the SU(3) FHM on the SL [71], does not seem to occur here, since the correlation patterns are monotonic [69].

<figure><img src="image_5.png" /><figcaption>Figure 4. (a) Charge gaps ($\Delta_c$) for the FHM on the $L = 12$ sites triangular lattice for $N=2,3$ and 4. The dashed lines (cf [69] for the fitting procedure) cross the x axis at $U \simeq 8.6$ for $N=2$, $U \simeq 9.8$ for $N=3$ and $U \simeq 11.2$ for $N=4$ separating the small U metallic phase from the Mott insulators. (b) Spin gaps ($\Delta_s$) defined as the difference between the SU($N$) adjoint and the SU($N$) singlets minimal energies.</figcaption></figure>

The presence of an in-between phase might occur when the HM limit is not a LRO, like for SU(3) on the Honeycomb lattice [72]. The SU(4) FHM on the TL might enter into this category as the HM limit is a gapless QSL [28, 40, 52]. In fact, while there are four low lying energy SU(4) singlet states for $U = 50$ on the $L = 12$ sites TL, when decreasing the interaction, the singlets gap starts increasing and an emerging Atos appears for $U \lesssim 15$ (cf Fig. 3 with both $U = 12$ and $U = 50$), which is similar to the one occurring in a pseudo HM with next-nearest neighbour (nnn) couplings on the same $L = 12$ TL cluster [28]. The nnn couplings, which are present at order 4 in $t/U$ in the large $U$ limit [51], were shown to stabilize a four sublattice order [28], also apparent in the correlation patterns of our Fig. 2. From the charge gap shown in Fig. 4, the boundary between this phase and the metallic phase occurs around $U = U_c = 11.2 \pm 0.2$. Like for the other systems, such a value changes for larger $L$ due to finite size effects, as illustrated for SU(2), $L = 12$ and $L = 16$ in [69].

Finally, we have investigated the SU(6) FHM on the SL, as the HM limit is also not a LRO but a plaquette state [39]. Through the correlations pattern of the GS on a $4 \times 3$ periodic SL for $U = 10$ in Fig. 2, we found some evidence of the 6 sites plaquette state in the Mott phase, a feature confirmed by the presence of two low lying energy SU(6) singlet states, compatible with the periodic boundary conditions, for large $U$ (i.e $U = 20$ cf Fig. 3). When $U$ decreases, the spin gap becomes smaller than the singlet gap, possibly revealing the Mott transition. However, with the current version of our code and with the limitation of our computational resources, the necessary calculation of the charge gap was too demanding,
leaving open both the question of the size of the metallic
phase and the presence of some intermediate phase.

To conclude, we found an efficient protocole to per-
form ED of the FHM on L-sites clusters directly in each
SU(N) irrep, which uses the set of ssYT (or GT patterns)
as a convenient basis with matrix elements of the U(L)
group generators. This approach, which generalizes the
use of SYT for SU(N) in HMs [56], dramatically reduces
the dimension of the matrices to diagonalize. We applied
our method to study the survival of the SU(N) Mott
phases for N = 3 to N = 6 on TL and SL when the on-

site interaction *U* decreases. In particular, we found an
emerging intermediate LRO phase for SU(4) on the TL,
reminiscent of the four sublattice order in the HM with
nnn couplings [28]. In the future, one possible perspec-
tive would be the implementation of the ssYT basis in
tensor networks/DMRG algorithm, in a fashion similar
to what was done with the SYT for the HM [73, 74].

We acknowledge F. Mila for useful discussions. This
work has been supported by an Emergence grant from
the Institute of Physics of CNRS.

[1] J. Hubbard, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences **276**, 238 (1963).

[2] M. C. Gutzwiller, Phys. Rev. Lett. **10**, 159 (1963).

[3] D. J. Scalapino, Rev. Mod. Phys. **84**, 1383 (2012).

[4] P. W. Anderson, Science **235**, 1196 (1987).

[5] F. C. Zhang and T. M. Rice, Phys. Rev. B **37**, 3759 (1988).

[6] D. P. Arovas, E. Berg, S. A. Kivelson, and S. Raghu, Annual Review of Condensed Matter Physics **13**, 239 (2022).

[7] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, Annual Review of Condensed Matter Physics **13**, 275 (2022).

[8] R. Assaraf, P. Azaria, M. Caffarel, and P. Lecheminant, Phys. Rev. B **60**, 2299 (1999).

[9] C. Honerkamp and W. Hofstetter, Phys. Rev. Lett. **92**, 170403 (2004).

[10] S. Capponi, P. Lecheminant, and K. Totsuka, Annals of Physics **367**, 50 (2016).

[11] I. Affleck, Nuclear Physics B **265**, 409 (1986).

[12] I. Affleck and J. B. Marston, Phys. Rev. B **37**, 3774 (1988).

[13] D. S. Rokhsar, Phys. Rev. B **42**, 2526 (1990).

[14] M. Marder, N. Papanicolaou, and G. C. Psaltakis, Phys. Rev. B **41**, 6920 (1990).

[15] S. K. Pati, R. R. P. Singh, and D. I. Khomskii, Phys. Rev. Lett. **81**, 5406 (1998).

[16] M. G. Yamada, M. Oshikawa, and G. Jackeli, Phys. Rev. Lett. **121**, 097201 (2018).

[17] Y.-H. Zhang, D. N. Sheng, and A. Vishwanath, Phys. Rev. Lett. **127**, 247701 (2021).

[18] A. V. Gorshkov, M. Hermele, V. Gurarie, C. Xu, P. S. Julienne, J. Ye, P. Zoller, E. Demler, M. D. Lukin, and A. Rey, Nature physics **6**, 289 (2010).

[19] M. A. Cazalilla and A. M. Rey, Reports on Progress in Physics **77**, 124401 (2014).

[20] S. Taie, R. Yamazaki, S. Sugawa, and Y. Takahashi, Nature Physics **8**, 825 (2012).

[21] C. Hofrichter, L. Riegger, F. Sazza, M. Höfer, D. R. Fernandes, I. Bloch, and S. Fölling, Phys. Rev. X **6**, 021030 (2016).

[22] B. Abeln, K. Sponselee, M. Diem, N. Pintul, K. Sengstock, and C. Becker, Phys. Rev. A **103**, 033315 (2021).

[23] S. Taie, E. Ibarra-García-Padilla, N. Nishizawa, Y. Takasu, Y. Kuno, H.-T. Wei, R. T. Scalettar, K. R. A. Hazzard, and Y. Takahashi, Nature Physics **18**, 1356 (2022).

[24] D. Tusi, L. Franchi, L. F. Livi, K. Baumann, D. Benedicto Orenes, L. Del Re, R. E. Barfknecht, T. W. Zhou, M. Inguscio, G. Cappellini, M. Capone, J. Catani, and L. Fallani, Nature Physics **18**, 1201 (2022).

[25] G. Pasqualetti, O. Bettermann, N. D. Oppong, E. Ibarra-García-Padilla, S. Dasgupta, R. T. Scalettar, K. R. A. Hazzard, I. Bloch, and S. Fölling, (2023), arXiv:2305.18967 [cond-mat.quant-gas].

[26] N. Read and S. Sachdev, Nuclear Physics B **316**, 609 (1989).

[27] Y. Q. Li, M. Ma, D. N. Shi, and F. C. Zhang, Phys. Rev. Lett. **81**, 3527 (1998).

[28] K. Penc, M. Mambrini, P. Fazekas, and F. Mila, Phys. Rev. B **68**, 012408 (2003).

[29] M. Greiter and S. Rachel, Phys. Rev. B **75**, 184441 (2007).

[30] P. Corboz, A. M. Läuchli, K. Penc, M. Troyer, and F. Mila, Phys. Rev. Lett. **107**, 215301 (2011).

[31] B. Bauer, P. Corboz, A. M. Läuchli, L. Messio, K. Penc, M. Troyer, and F. Mila, Phys. Rev. B **85**, 125116 (2012).

[32] D. Wang, Y. Li, Z. Cai, Z. Zhou, Y. Wang, and C. Wu, Phys. Rev. Lett. **112**, 156403 (2014).

[33] R. Bondesan and T. Quella, Nuclear Physics B **886**, 483 (2014).

[34] A. Weichselbaum, S. Capponi, P. Lecheminant, A. M. Tsvelik, and A. M. Läuchli, Phys. Rev. B **98**, 085104 (2018).

[35] A. Paramekanti and J. B. Marston, Journal of Physics: Condensed Matter **19**, 125215 (2007).

[36] P. Corboz, K. Penc, F. Mila, and A. M. Läuchli, Phys. Rev. B **86**, 041106 (2012).

[37] T. C. Lang, Z. Y. Meng, A. Muramatsu, S. Wessel, and F. F. Assaad, Phys. Rev. Lett. **111**, 066401 (2013).

[38] P. Nataf, M. Lajkó, P. Corboz, A. M. Läuchli, K. Penc, and F. Mila, Phys. Rev. B **93**, 201113 (2016).

[39] O. Gauthé, S. Capponi, M. Mambrini, and D. Poilblanc, Phys. Rev. B **101**, 205144 (2020).

[40] H.-K. Jin, R.-Y. Sun, H.-H. Tu, and Y. Zhou, Science Bulletin **67**, 918 (2022).

[41] P. W. Anderson, Phys. Rev. **86**, 694 (1952).

[42] E. Manousakis, Rev. Mod. Phys. **63**, 1 (1991).

[43] B. Bernu, P. Lecheminant, C. Lhuillier, and L. Pierre, Phys. Rev. B **50**, 10048 (1994).

[44] L. Capriotti, A. E. Trumper, and S. Sorella, Phys. Rev. Lett. **82**, 3899 (1999).

[45] F. Wang and A. Vishwanath, Phys. Rev. B **80**, 064413 (2009).
[46] M. Hermele and V. Gurarie, Phys. Rev. B **84**, 174441 (2011).

[47] P. Corboz, M. Lajkó, A. M. Läuchli, K. Penc, and F. Mila, Phys. Rev. X **2**, 041013 (2012).

[48] H.-H. Lai, Phys. Rev. B **87**, 205131 (2013).

[49] G. Chen, K. R. A. Hazzard, A. M. Rey, and M. Hermele, Phys. Rev. A **93**, 061601 (2016).

[50] P. Nataf, M. Lajkó, A. Wietek, K. Penc, F. Mila, and A. M. Läuchli, Phys. Rev. Lett. **117**, 167202 (2016).

[51] C. Boos, C. J. Ganahl, M. Lajkó, P. Nataf, A. M. Läuchli, K. Penc, K. P. Schmidt, and F. Mila, Phys. Rev. Res. **2**, 023098 (2020).

[52] A. Keselman, B. Bauer, C. Xu, and C.-M. Jian, Phys. Rev. Lett. **125**, 117202 (2020).

[53] J.-Y. Chen, J.-W. Li, P. Nataf, S. Capponi, M. Mambrini, K. Totsuka, H.-H. Tu, A. Weichselbaum, J. von Delft, and D. Poilblanc, Phys. Rev. B **104**, 235104 (2021).

[54] P. W. Anderson, Materials Research Bulletin **8**, 153 (1973).

[55] B. Sutherland, Phys. Rev. B **12**, 3795 (1975).

[56] P. Nataf and F. Mila, Phys. Rev. Lett. **113**, 127204 (2014).

[57] K. Wan, P. Nataf, and F. Mila, Phys. Rev. B **96**, 115159 (2017).

[58] Or more precisely of its universal enveloping algebra since it contains some square of generators.

[59] The irreps of $U(N)$ or $SU(N)$ are basically the same; distinction is relevant when we include/exclude representation of generators with non vanishing trace like some combinations of $E_{ii}$.

[60] H. Weyl, Math. Z. **23**, 271 (1925).

[61] G. de B. Robinson, *Representation theory of the symmetric group* (Edinburgh: EUP, 1961).

[62] I. M. Gelfand and M. L. Tsetlin, Dokl. Akad. Nauk SSSR **71**, 825 (1950).

[63] N. J. Vilenkin and A. U. Klimyk, *Representation of Lie groups and special functions, Vol. 3* (Kluwer Academic Publishers, 1992).

[64] A. Alex, M. Kalus, A. Huckleberry, and J. von Delft, Journal of Mathematical Physics **52**, 023507 (2011).

[65] A. Young, *Proceedings of the London Mathematical Society s1-33*, **97** (1900).

[66] H. Tsunetsugu and M. Arikawa, *Journal of the Physical Society of Japan* **75**, 083701 (2006).

[67] A. Läuchli, F. Mila, and K. Penc, Phys. Rev. Lett. **97**, 087205 (2006).

[68] K. Pilch and A. N. Schellekens, *Journal of Mathematical Physics* **25**, 3455 (1984).

[69] See Supplemental Material at XXX for some complementary numerical results.

[70] A. Szasz, J. Motruk, M. P. Zaletel, and J. E. Moore, Phys. Rev. X **10**, 021042 (2020).

[71] C. Feng, E. Ibarra-García-Padilla, K. R. A. Hazzard, R. Scalettar, S. Zhang, and E. Vitali, (2023), arXiv:2306.16464 [cond-mat.str-el].

[72] S. S. Chung and P. Corboz, Phys. Rev. B **100**, 035134 (2019).

[73] P. Nataf and F. Mila, Phys. Rev. B **97**, 134420 (2018).

[74] S. Gozel, P. Nataf, and F. Mila, Phys. Rev. Lett. **125**, 057202 (2020).
# Exact Diagonalization of SU(N) Fermi-Hubbard Models

## – Supplemental Material –

Thomas Botzung¹, and Pierre Nataf¹

¹ Laboratoire de Physique et Modélisation des Milieux Condensés, Université Grenoble Alpes and CNRS, 25 avenue des Martyrs, 38042 Grenoble, France (Dated: September 19, 2023)

In this supplemental material, we show complementary numerical results obtained via ED. In Sec. I, we focus on SU(5) on the square lattice with L=10 sites. For this specific case, we show some energy spectra and we detail the fitting procedure of the charge gap that we have used for all the systems under consideration in our manuscript. Then, in Sec. II, we consider the case SU(2) on a triangular lattice of L = 12 and L = 16 sites which illustrate the finite size effects on the spin/charge gaps. Finally, in Sec. III, we show the evolution of real space correlations while increasing the interaction $U$ for SU(3) on the L=12 triangular lattice and for SU(5) on the square lattice.

### I. SU(5) ON A SQUARE LATTICE

In this section, we show the data for SU(5) on a square lattice with L=10 sites at filling 1/5, that we have discussed in the main text. The long-range color ordered state with broken SU(5) symmetry is revealed by an Anderson tower of states (Atos), clearly visible on Fig. S1 (b) for $U = 10$. The situation is quite different at small $U$, as seen from Fig. S1 (a), where the Atos has disappeared and an energy plateau appears. In Fig. S1 (c), we plot the charge gap $\Delta_c$ to locate the boundary $U = U_c$ between the metallic and the LRO phases. We see that $\Delta_c$ opens rapidly around $U \approx 8$.

To estimate the critical point roughly, we use a fit of the form $\Delta_c = a * U - b$. The latter corresponds to the expected behavior of the charge gap at large $U$, since the energy cost to add one particle should be equal to $U$: using the definition of $\Delta_c$ one can readily show that for $U \to \infty$, $\Delta_c = U$. Then, such a linear function shown in dashed lines in the figures cross the x axis at $U = U_c$. However, due to the small curvature of $\Delta_c$, it depends on the values of $U$ used for the fit; we took different sets belonging to $U \in [40; 100]$. With this method, we obtain $8.6 < U_c < 8.9$ with $0.992 < a < 0.997$. In addition, in Fig. S1 (d), we show the spin gap, defined as the difference between the SU(5) adjoint and the SU(5) singlet minimal energies. The spin gap increases rapidly in the metallic phase, reaching its maximal value around $U \approx 8$, close to $U = U_c$, a feature also observed in the other cases studied here.

<figure><img src="image_6.png" /></figure>

Figure S1. Energy spectra of the FHM in Eq. (1) (main text) (with $t_{ij} = 1$ at filling 1/N) as a function of the quadratic casimir $C_2$ for $N = 5$, for the L=10 periodic squared lattice for various values of $U$, (a) $U = 0.5$, (b) $U = 10$. In (b), the black line is a guide for the eye to show the presence of the Anderson tower of states. (c) Charge gap for the FHM on the L = 10 squared periodic lattice as a function of the one-site repulsion $U$. The dashed line (fit of the form $\Delta_c = a * U - b$) crosses the x axis at $U_c \approx 8.75$, separating the small U metallic phase from the Mott insulators. (d) Spin gaps as a function of $U$. It exhibits a peak at $U \approx 8$.
<figure><img src="image_7.png" /></figure>

Figure S2. (a) Charge gaps for the FHM on the $L = 12$ and $L = 16$ sites triangular lattice for $N=2$. The dashed lines correspond to a fit of the form $\Delta_c = a * U - b$. They cross the x axis at $U = U_c = 8.6 \pm 0.3$ for $L = 12$ and at $U = U_c = 7.9 \pm 0.3$ for $L = 16$, separating the small U metallic phase from the Mott insulators. (b) Spin gaps for FHM on the $L = 12$ and $L = 16$ sites triangular lattice for $N=2$. A peak is present at $U \approx 8$ in both cases.

## II. SU(2) ON THE TRIANGULAR LATTICES WITH $L=12$ AND $L=16$ SITES

In this section, we analyze the case SU(2) on a triangular lattice with $L = 12$ and $L = 16$ sites at filling $1/2$. In Fig. I (a), we show the charge gap as a function of $U$ for $L = 12$ (blue diamond) and $L = 16$ (red circle). Using the fitting procedure detailed in the Sec. S1 above, we estimate the phase boundary between a metallic phase and an insulating phase: we have $U_c = 8.6 \pm 0.3$ for $L = 12$ and $U_c = 7.9 \pm 0.3$ for $L = 16$. In Fig. I (b), we show the spin gap as a function of $U$.

## III. EVOLUTION OF CORRELATIONS WITH U

In this section, we study how the real space correlations, defined as $\langle P_{1j} \rangle - 1/N$ where $P_{1j} \equiv -1 + E_{1j}E_{j1}$, evolve while increasing the on-site interaction $U$. In Fig. S3, we show the correlations for two lattice geometries and different $U$. From left to right, we have $U=1, U=5, U=10, $ and $U=20$. In (a), we consider the case SU(3) on a triangular lattice. We see that for $U = 1$, the correlations do not display the three-sublattice Neel order. The order becomes clearly visible at $U = 5$. In (b), the case SU(5) on a square lattice is shown. Here, the pattern compatible with LRO is visible at any $U$, although with smaller correlation strengths for small $U$. Note that while increasing the on-site interaction $U$, we do not observe the appearance of a new ( or differently configured) long-range orders. For completeness, we have checked that patterns do not evolve for $U$ up to 50 (not shown).
<figure><img src="image_8.png" /><figcaption>Figure S3. Evolution of the correlations patterns of the ground states of the FHM at filling 1/N for various from left to right $U = 1, 5, 10, 20$. (a) On a triangular lattice SU(3). (b) On a square lattice SU(5). They are defined as $\langle P_{1j} \rangle - 1/N$ where $P_{1j} = -1 + E_{1j}E_{j1}$, with the reference site 1 in black, and j the site indices being blue (resp. red) for positive (resp. negative) correlation, with area proportional to its absolute value.</figcaption></figure>