Exact diagonalization of $\mathrm{SU}(N)$ Heisenberg and AKLT chains using the full $\mathrm{SU}(N)$ symmetry

TL;DR

This paper introduces a method for exact diagonalization of $ ext{SU}(N)$ spin chains using Young tableaux, enabling analysis of larger systems and investigation of edge states in $ ext{SU}(N)$ models.

Contribution

It develops a generalized scheme for exact diagonalization of $ ext{SU}(N)$ chains using Young tableaux, improving efficiency and scalability over previous methods.

Findings

Efficient basis construction for $ ext{SU}(N)$ subsectors

Analysis of edge states in $ ext{SU}(N)$ Heisenberg and AKLT models

Ability to study larger systems than before

Abstract

We present a method for the exact diagonalization of the $\mathrm{SU}(N)$ Heisenberg interaction Hamiltonian, using Young tableaux to work directly in each irreducible representation of the global $\mathrm{SU}(N)$ group. This generalized scheme is applicable to chains consisting of several particles per site, with any $\mathrm{SU}(N)$ symmetry at each site. Extending some of the key results of substitutional analysis, we demonstrate how basis states can be efficiently constructed for the relevant $\mathrm{SU}(N)$ subsector, which, especially with increasing values of $N$ or numbers of sites, has a much smaller dimension than the full Hilbert space. This allows us to analyze systems of larger sizes than can be handled by existing techniques. We apply this method to investigate the presence of edge states in $\mathrm{SU}(N)$ Heisenberg and AKLT Hamiltonians.

Equations

$$H=\sum _i{H}_{(i)}+\sum _{i<j}{H}_{(i,j)},$$

$$H_{(i)}=\sum _{k_i<l_i}{J}_{k_i{l}_i}{P}_{k_i{l}_i}$$

$$H_{(i,j)}=J_{ij}\sum _{k_i,l_j}{P}_{k_i,l_j},$$

$$[H,H_{(i)}]=0,\forall i=1,...,N_s.$$

$$H_{(i,j)}=J_{ij}\sum _{\mu ,\nu }{\hat{S}}_{\mu \nu }^i{\hat{S}}_{\nu \mu }^j,$$

$$\left[{\hat{S}}_{\alpha \beta }^i,{\hat{S}}_{\mu \nu }^i\right]={\delta }_{\mu \beta }{\hat{S}}_{\alpha \nu }^i-{\delta }_{\alpha \nu }{\hat{S}}_{\mu \beta }^i.$$

$${\hat{S}}_{\mu \nu }^i=\sum _{k_i=m(i-1)+1}^{mi}|{\mu }_{k_i}⟩⟨{\nu }_{k_i}|-\frac{m{\delta }_{\mu ,\nu }}{N}$$

$$\sum _{\mu ,\nu }{\hat{S}}_{\mu \nu }^i{\hat{S}}_{\nu \mu }^j$$

$$\displaystyle=\sum_{k_{i},l_{j}}\Big{\{}\sum_{\mu,\nu}|\mu_{k_{i}}\rangle\otimes|\nu_{l_{j}}\rangle\langle\nu_{k_{i}}|\otimes\langle\mu_{l_{j}}|\Big{\}}-\frac{m^{2}}{N}.$$

$$H_{(i,j)}=J_{ij}\sum _{k_i,l_j}{P}_{k_i,l_j},$$

$$H_{\text{interaction}}=\sum _{i<j}{H}_{(i,j)}.$$

$$H_{\text{interaction}}=\sum _{i=1}^{N_s-1}{H}_{(i,i+1)},$$

$$H_{\text{interaction}}=H_{(1,2)}+H_{(N_s-1,N_s)}+\sum _{i=1}^{N_s-2}{H}_{(i,i+2)}.$$

$$\mathrm{Proj}=\prod _{i=1}^{N_s}{\mathrm{Proj}}^{\beta (i)}(i).$$

$${\mathbf{Proj}}^{\beta (i)}(i)={\mathbf{v}}_j^{\beta }{\left({\mathbf{v}}_j^{\beta }\right)}^T,$$

$${\mathrm{Proj}}^{\beta (i)}(i)=\sum _{r,s}[{\mathbf{Proj}}^{\beta (i)}(i){]}_{rs}{o}_{rs}^{\beta (i)},$$

$$H_{(i)}=J_{1_i{2}_i}{P}_{1_i{2}_i}+J_{1_i{3}_i}{P}_{1_i{3}_i}+J_{2_i{3}_i}{P}_{2_i{3}_i},$$

$${\mathbf{H}}_{(i)}^{[3]}$$

$${\mathbf{H}}_{(i)}^{[2,1]}$$

$${\mathbf{H}}_{(i)}^{[1,1,1]}$$

$${\text{Proj}}_{\pm }^{[2,1]}(i)=$$

$$+(J_{12}-2J_{13}+J_{23})P_{1_i{3}_i}$$

$$+(J_{12}+J_{13}-2J_{23})P_{2_i{3}_i}]$$

$$+\frac{1}{6}(2Id-P_{1_i{2}_i}{P}_{2_i{3}_i}-P_{2_i{3}_i}{P}_{1_i{2}_i}),$$

$${\mathrm{Proj}}^{[2,1]}(1)=\frac{1}{6}(2Id+2P_{12}-P_{13}-P_{23}-P_{12}{P}_{13}-P_{13}{P}_{12})$$

$${\mathrm{Proj}}^{[2,1]}(2)=\frac{1}{6}(2Id+2P_{45}-P_{46}-P_{56}-P_{45}{P}_{46}-P_{46}{P}_{45}).$$

$$\text{ytableaushort}1234,568,79\ast [\ast (lightblue)]3\ast [\ast (lightpink)]3+1,2\ast [\ast (lightgrey)]0,2+1,2$$

$$\text{ytableaushort}1235,468,79\ast [\ast (lightblue)]3\ast [\ast (lightpink)]3+1,2\ast [\ast (lightgrey)]0,2+1,2$$

$$\text{ytableaushort}1236,459,78\ast [\ast (lightblue)]3\ast [\ast (lightpink)]3+1,2\ast [\ast (lightgrey)]0,2+1,2$$

$$\text{ydiagram}[\ast (lightblue)]3\ast [\ast (lightpink)]3+1,2\ast [\ast (lightgrey)]0,2+1,2\text{b) }$$

$$\text{ydiagram}[\ast (lightblue)]2,1\ast [\ast (lightpink)]2+2,1+1\ast [\ast (lightgrey)]0,2+1,2\text{f) }$$

$$\text{ydiagram}[\ast (lightblue)]2,1\ast [\ast (lightpink)]0,1+1,2\ast [\ast (lightgrey)]2+2,2+1\text{j) }$$

$$\mathrm{Proj}|\Psi ⟩=|\Psi ⟩.$$

$$\prod _{i=1}^{N_s}{d}_N^{\beta (i)}=\sum _{\alpha }{d}_N^{\alpha }{D}^{\alpha }(\mathcal{B}).$$

$$\underset{i=1}{\overset{N_s}{⨂}}\beta (i)$$

$${\mathrm{Proj}}^{\beta (i)}(i)|{\Psi }_i⟩=|{\Psi }_i⟩$$

$$|\Psi ⟩=\underset{i=1}{\overset{N_s}{⨂}}|{\Psi }_i⟩$$

Full text

I I. Introduction

In recent years, considerable progress has been made in experiments with ultracold atoms 1, enabling the realization of sophisticated quantum many-body systems. In particular, degenerate gases of strontium and ytterbium loaded in optical lattices have been used to simulate the $\mathrm{SU}(N)$ Fermi-Hubbard model 2; 3; 4; 5; 6; 7; 8; 9, a generalization of the familiar $\mathrm{SU}(2)$ spin-$1/2$ Fermi-Hubbard model. When the number of particles per site is an integer and the on-site repulsion is sufficiently large, the systems are expected to be in Mott insulating phases, which are well-described by $\mathrm{SU}(N)$ Heisenberg models. This class of models is a unique playground for strongly correlated systems as it encompasses a wide variety of quantum ground states with different physical properties. In fact, even for the simplest cases with interactions limited to nearest neighbours, the zero-temperature quantum phases can be very diverse and can depend on the geometry of the lattices (one-dimensional chain, two-dimensional bipartite or frustrated lattice), the number of colours (i.e., the value of $N$), and the local $\mathrm{SU}(N)$ symmetry of the wave function.

At each site, the local $\mathrm{SU}(N)$ symmetry corresponds to a specific irreducible representation (“irrep”) of $\mathrm{SU}(N)$, and, for $m$ particles per site, can be encoded by Young diagram with $m$ boxes and no more than $N$ rows. For $m=1$, the Young diagram is a single box, representing the fundamental irrep. In this case, the $\mathrm{SU}(N)$ chain, for which a general Bethe ansatz solution exists 10, is gapless, with algebraic decaying correlations. However, if a second particle is added to each site in such a way that the resultant local wave function is fully symmetric, then, in the case of $N=2$, the system can open a Haldane gap11, while for $N>2$, the chain should be critical with universality class $\mathrm{SU}(N)_{1}$ ( although this issue has not yet been completely solved from a numerical point of view12). In two dimensions, the ground state of a square lattice with $m=1$ particle per site has been shown to be characterized by some Néel-type ordering for $\mathrm{SU}(2)$, $\mathrm{SU}(3)$13; 14, $\mathrm{SU}(4)$15 and $\mathrm{SU}(5)$16, whereas when there are $m>1$ particles per site in an antisymmetric representation, the ground state is predicted by mean-field theory to be a chiral spin liquid, provided that $m/N>5$ 17; 18.

From an experimental perspective, the study of $\mathrm{SU}(N)$ Heisenberg models with $N>2$ and non-fundamental irreps (i.e., $m>1$) on each site is indeed relevant. It has been shown, for instance, that certain many-body systems exhibit $\mathrm{SU}(N)$ symmetry with $N$ as large as 10, as in the case of strontium-879; 8, and may be implemented using optical lattices with two atoms on each site 5; 19. Furthermore, the realization of exotic phases of matter is naturally expected to involve irreps of mixed symmetry, that is, neither fully symmetric (corresponding to Young diagrams with one row) nor fully antisymmetric (corresponding to Young diagrams of one column). According to mean-field caculations17; 18, in order to obtain non-Abelian chiral spin liquids on the square lattice, the number of columns of the local irreps should be at least two. Another example is the $\mathrm{SU}(N)$ symmetry-protected topological phases in one dimension20; 21, which are generalizations of the $\mathrm{SU}(2)$ spin-1 Haldane phases22; 23; 24. These are gapped phases with non-trivial edge states, and the paradigmatic Hamiltonians of those states are the $\mathrm{SU}(N)$ version of the AKLT chain25; 26.They involve irreps with multiple rows and columns at each site and can lead to $N$ distinct topological phases, classified using group cohomology27.

The theoretical study of such systems can be extremely challenging, due in no small part to the inherent limitations of current numerical methods. The Density Matrix Renormalization Group (DMRG) technique has proven rather efficient in the investigation of $\mathrm{SU}(N)$ Hamiltonians in 1D 23; 20; 28; 29; 30 as well as Infinite Projected Entangled Pair States (iPEPS) in 2D 31; 32; 33; 15. However, as the local Hilbert space dimensions increases as a result of increasing the number of colours or the number of particles per site, the performance of DRMG significantly deteriorates. Quantum Monte Carlo methods, on the other hand, are usually able to accommodate large Hilbert spaces, but can only be used in very specific cases (to avoid the sign problem), namely, for chains with one particle per site, or for bipartite lattices on which pairs of interacting sites correspond to conjugate irreducible representations34; 35; 36; 37; 38. For configurations where the wave function on each site is completely antisymmetric, variational Monte Carlo simulations based on Gutzwiller projected wave functions have been found to produce remarkably accurate results 39; 40; 41; 42; 43, but it is not clear how this approach can be generalized to other local irreps. Finally, exact diagonalization (ED) is limited to small clusters.

Recently, we developed a procedure to exactly diagonalize the Hamiltonian for one particle per site independently in each global irrep of $\mathrm{SU}(N)$, using standard Young tableaux (SYTs) and the orthogonal representation of the symmetric group16. This method obviates the use of Clebsch-Gordan coefficients, for which the computation complexity increases dramatically with $N$ 44. The method was then extended to chains with a fully symmetric or antisymmetric irrep at each site12—a rather straightforward extension, as the number of multiplets per site is still one. However, the number of multiplets is greater than one for systems where local irreps have more than one row and more than one column, and it was not clear how SYTs could be used to solve such systems.

The purpose of this article is to proceed to such a development, and the method presented in the following two sections can be applied to the most general configuration, with one or more particles per site and any given irrep at each site (not necessarily the same from site to site). Following a brief introduction on the theory behind the method, in Section II we show how to derive a projection operator which imposes the local $\mathrm{SU}(N)$ symmetry at each site. Then, in Section III, we describe an efficient algorithm for constructing suitable basis states for a for a given global $\mathrm{SU}(N)$ subsector, using this projection operator. Further simplifications to this algorithm, which take full advantage of the inherent symmetry of the problem, are detailed in the Appendix. In Section IV, the method is used to investigate the $\mathrm{SU}(4)$ Heisenberg chain with irrep $[2,2]$ at each site: After demonstrating how to systematically express an $\mathrm{SU}(N)$ AKLT Hamiltonian in terms of permutations, we calculate the energy of the edge states directly in their corresponding irreps and determine whether they remain in the lowest part of the spectrum as we move from the AKLT point to the Heisenberg point in an interpolating Hamiltonian. Finally, conclusions are drawn and future directions are discussed in Section V.

II The Method

II.1 The Hamiltonian as a sum of permutations

As in Ref. 16, the construction of the Hamiltonian matrix will rely on the very simple representation of permutations in the basis of standard Young tableaux. It is therefore most convenient to define the model directly in terms of permutations. The equivalence of this to other formulations will be established in the following subsection.

Consider a general $\mathrm{SU}(N)$ model where, at each site, there are $m$ particles in the fundamental representation. Denoting the number of sites by $N_{s}$, we have $mN_{s}$ particles in total. We assign the number $k_{i}\equiv m(i-1)+k$ to the $k$-th particle of site $i$ ($i=1,\dots,N_{s};k=1,\dots,m$). The general Hamiltonian we will consider is a linear combination of intra- and inter-site permutations:

[TABLE]

where $J_{k_{i}l_{i}}$ and $J_{ij}$ are the intra- and inter-site coupling constants, respectively. Note that the inter-site operator $H_{(i,j)}$ couples all the particles of site $i$ to all the particles of site $j\neq i$ with the same coupling $J_{ij}$. Consequently, $H_{(i,j)}$ commutes with all intra-site permutations $P_{k_{i}l_{i}}$ and $P_{k_{j}l_{j}}$, and of course with all $P_{k_{n}l_{n}}$ for $n\neq i,j$. It follows that $H_{(i,j)}$ commutes with all intra-site operators $H_{(i)}$. Moreover, the intra-site operators commute with each other. We thus obtain the following fundamental property:

[TABLE]

Therefore, we can diagonalize $H$ and all of the $H_{(i)}$ in a common basis. Accordingly, the Hilbert space can be partitioned into sectors corresponding to the set of eigenvalues $e_{i},(i=1,...,N_{s})$ of the operators $H_{(i)}$, and the total Hamiltonian $H$ is block diagonal in this basis.

Since the Hamiltonians $H_{(i)}$ are linear combinations of permutations, they are $\mathrm{SU}(N)$ invariant. Hence, assuming that the spectrum of $H_{(i)}$ has no accidental degeneracy, each eigenvalue $e_{i}$ defines a local subspace of the Hilbert space at site $i$ that belongs to a local irrep, which we will denote by $\beta(i)$. The Hamiltonian $H$ restricted to the corresponding block is the $\mathrm{SU}(N)$ model we want to study, with irrep $\beta(i)$ at site $i$, and a reference energy $\sum_{i}e_{i}$.

To write the Hamiltonian in this block, we construct a projection operator $\mathrm{Proj}$ that maps onto the corresponding subspace. As we will show in subsection II.3, the explicit form of such a projection operator can be derived exactly using general properties of the permutation group, without diagonalizing the local Hamiltonians $H_{(i)}$. Hence, our model can alternatively be characterized by the Hamiltonian $H^{\prime}=\sum_{i<j}H_{(i,j)}$, restricted to the projected Hilbert space.

Note that the construction can easily be extended to the case where the number of particles $m_{i}$ at site $i$ varies from site to site. For simplicity in the examples to follow, we will focus on chains with the same number of particles at each site, i.e., $m_{i}=m$ independent of $i$.

II.2 The $\mathrm{SU}(N)$ formulation of the problem

In the most general case, an $\mathrm{SU}(N)$ Heisenberg-like interaction between two sites $i$ and $j$ can be written as:

[TABLE]

where the $\mathrm{SU}(N)$ generators satisfy at each site $i$ the commutation relation:

[TABLE]

With $m$ particles per site, the $\mathrm{SU}(N)$ generators for site $i$ are

[TABLE]

where the operator $|\mu_{k_{i}}\rangle\langle\nu_{k_{i}}|$ acts on the particle $k_{i}$ to change its colour from $|\nu_{k_{i}}\rangle$ to $|\mu_{k_{i}}\rangle$. The symbols $\mu_{k_{i}}$ and $\nu_{k_{i}}$ each stands for one of the $N$ colours, $A,B,C$, etc. The term $-\frac{m\delta_{\mu,\nu}}{N}$ renders the generators traceless. This set of local generators satisfies the commutation relation in Eq. (6). The $\mathrm{SU}(N)$ Heisenberg interaction between two sites $i$ and $j$ shown in Eq. (5) can then be rewritten as:

[TABLE]

The constant $-m^{2}/N$ may be dropped.

The term inside braces is the permutation operator $P_{k_{i}l_{i}}$, which interchanges the $k$-th particle of site $i$ and the $l$-th particle of site $j$. Thus, the interaction Hamiltonian between sites $i$ and $j$ couples each of the $m$ particles at site $i$ to each of the $m$ particles at site $j$ [as shown in Fig. 1(a)   for the case of $m=3$ particles per site]:

[TABLE]

and the general $\mathrm{SU}(N)$ Heisenberg Hamiltonian is of the form

[TABLE]

In particular, the Hamiltonian for an entire chain of $N_{s}$ sites, in which each site interacts with the site(s) adjacent to it, is simply

[TABLE]

in the case of open boundary conditions, where we have indexed adjacent sites with consecutive numbers, as demonstrated in Fig. 1(b)   for $N_{s}=7$. We can treat periodic boundary conditions by adopting the computationally convenient indexing convention shown in Fig. 1(c) , and writing

[TABLE]

II.3 Projection operators

We now describe the procedure for writing the requisite projection operator as a linear superposition of permutations. Using a projection operator is essential: For a system of $N_{s}$ sites with $m$ particles per site, the full Hilbert space has dimension $N^{mN_{s}}$ and can be very large even for a small number of sites. Solving for the energies of an interaction Hamiltonian in this space would require diagonalizing a matrix of size $N^{mN_{s}}\times N^{mN_{s}}$, and the eigenvalues would include the spectra for all possible combinations of local irreps. Furthermore, this set of eigenenergies would encompass all of the different global irreps, since the full (reducible) Hilbert space can be decomposed as $\oplus_{\alpha}V^{\alpha}$, where the $\alpha$ are $\mathrm{SU}(N)$ irreps [refer to subection VII.1 for a review]. To obtain only the spectrum associated with a given global irrep $\alpha$ and a combination of specific local irreps, i.e., a given $\beta(i)$ at each site $i$, we must apply on the sector $V^{\alpha}$ the projection operator

[TABLE]

$\mathrm{Proj}$ is formulated as a product of $N_{s}$ operators, one for each site, where $\mathrm{Proj}^{\beta(i)}(i)$ imposes the symmetry associated with $\beta(i)$ at site $i$.

Each $\mathrm{Proj}^{\beta(i)}(i)$ is a linear superposition of the $m!$ permutations among the $m$ particle numbers of site $i$. Its exact form can be determined analytically as follows.

The local eigenstates at site $i$ that belong to irrep $\beta(i)$ correspond to the eigenvectors of $\mathbf{H}_{(i)}^{\beta(i)}$, the matrix of the intra-site coupling Hamiltonian defined in Eq. (2). We can write $\mathbf{H}_{(i)}^{\beta(i)}$ using the rules for the construction of Young’s orthogonal representation of the symmetric group, provided in appendix VII.2. In this representation, the basis consists of the orthogonal units $o_{rs}^{\beta(i)}$ of $\beta(i)$. This basis of orthogonal units $o_{rs}^{\beta(i)}$ is directly related to the $f^{\beta(i)}$ standard Young tableaux (SYTs) of shape $\beta(i)$ [cf. Fig. 5 in VII.2]. It suffices for $\mathrm{Proj}^{\beta(i)}(i)$ to project onto just one of the eigenstates of $\beta(i)$. Thus, if we denote the orthonormal set of eigenvectors of $\mathbf{H}_{(i)}^{\beta(i)}$ by $\{\mathbf{v}_{1}^{\beta},\dots,\mathbf{v}_{f^{\beta(i)}}^{\beta}\}$, the matrix representation for $\mathrm{Proj}^{\beta(i)}(i)$ can be very simply calculated as

[TABLE]

using any eigenvector $\mathbf{v}_{j}$ ($1\leq j\leq f^{\beta(i)}$).

Then, to convert $\mathrm{Proj}^{\beta(i)}(i)$ from its matrix form above to a linear combination of permutation operators, we substitute the explicit form of the orthogonal units $o_{rs}^{\beta(i)}$ [given in appendix VII.3] into

[TABLE]

where $[\mathbf{Proj}^{\beta(i)}(i)]_{rs}$ is the matrix element in the $r$-th row and $s$-th column.

In the example of $m=3$, the local Hamiltonian of site $i$ is

[TABLE]

and by applying the rules in subsection VII.2 to each of the three irreps $[3]=$ \ydiagram3, $[2,1]=$ \ydiagram2,1, and $[1,1,1]=$ \ydiagram1,1,1  ($N\geq 3$), the matrix representations of $H_{(i)}$ are:

[TABLE]

where, for the sake of brevity, we have omitted the subscripts $i$ in the coupling constants. In particular, we can address the following question: given a set of constants $J_{k_{i}l_{i}}$, to which irrep does the ground state at site $i$ belong? This entails determining the irrep $\beta(i)$ for which the matrix $\mathbf{H}_{(i)}^{\beta(i)}$ has the lowest eigenvalue. Figure 2 depicts an algebraic example.

Applying the above procedure [Eqs. (14) - (15)] to the local irrep $\beta(i)=[2,1]$, we find:

[TABLE]

where $P_{k_{i}l_{i}}$ is the permutation operator between particles $k_{i}$ and $l_{i}$, $Id$ is the identity operator, and $\lambda^{[2,1]}=\sqrt{{J_{12}}^{2}+{J_{13}}^{2}+{J_{23}}^{2}+J_{12}J_{13}+J_{12}J_{23}+J_{13}J_{23}}$.

This formula contains two projection operators, one obtained from each of the two eigenvectors of $\mathbf{H}_{(i)}^{[2,1]}$ via Eq. (14): the positive sign corresponds to $\mathbf{v}_{+}^{[2,1]}(\mathbf{v}_{+}^{[2,1]})^{T}$, where $\mathbf{v}_{+}^{[2,1]}$ is the eigenvector with the positive eigenvalue $+\lambda^{[2,1]}$, while the negative sign corresponds to $\mathbf{v}_{-}^{[2,1]}(\mathbf{v}_{-}^{[2,1]})^{T}$. In general, for a local irrep $\beta(i)$ at site $i$, there are $f^{\beta(i)}$ local projection operators for site $i$, one for each of the $f^{\beta(i)}$ eigenstates of $\mathbf{H}_{(i)}^{\beta(i)}$. In practice, however, when we consider a Hamiltonian for interactions between different sites [Eq. (10)], only one of these $f^{\beta(i)}$ operators is required for each site, and the specific choice is inconsequential. Since the Hamiltonian for the internal coupling between particles of the same site, $H_{\textup{internal}}=\sum_{i=1}^{N_{s}}H_{(i)}$, commutes with the Hamiltonian $H=H_{\textup{internal}}+H_{\textup{interaction}}$ [cf. Eq. 4], the spectrum of $H$ is equal to the spectrum of $H_{\textup{interaction}}$ shifted by that of $H_{\textup{internal}}$.

Therefore, if we have, for instance, the irrep $\beta=[2,1]$ at every site in a chain, the eigenvalues of $H_{\textup{internal}}$ can be found directly from the matrices $\mathbf{H}_{(i)}^{[2,1]}$ [Eq. (17)], while to study $H_{\textup{interaction}}$, we can take either one of the two local operators in (19)—say, $\mathrm{Proj}_{+}^{[2,1]}(i)$—and write $\mathrm{Proj}=\prod_{i=1}^{N_{s}}\mathrm{Proj}_{+}^{[2,1]}(i)$ as our projection operator for the chain.

Furthermore, since $H_{\textup{interaction}}$ has no dependence on the local coupling constants $J_{k_{i}l_{i}}$, artificial values can be substituted for the coupling constants $J_{k_{i}l_{i}}$ into the general formula for $\mathrm{Proj}^{\beta(i)}(i)$ (such as that in Eq. (19) for $\beta(i)=[2,1]$) without affecting the final spectrum obtained for $H_{\textup{interaction}}$. In particular, if $J_{12}=1$ and $J_{13}=J_{23}=0$, the operators in Eq. (19) reduce to the same forms as the orthogonal units $o_{11}^{[2,1]}$ and $o_{22}^{[2,1]}$ [cf. appendix VII.3 for an explicit expression]. This result is general: given any local irrep $\beta(i)$ at site $i$, if we choose $J_{12}=1$ as the only non-vanishing constant, the $f^{\beta(i)}$ projection operators are equivalent to the explicit form of the orthogonal units $o_{rr}^{\beta(i)}(\forall r=1,\dots,f^{\beta(i)})$, with the indices of the permutation operators adjusted to match the particle numbers of the site (i.e., $k\to k_{i}=m(i-1)+k$). Any one of these operators can be used in the subsequent calculations for $H_{\textup{interaction}}$. In our implementation, for instance, we simply use $o_{11}^{\beta(i)}$ at every site $i$.

III The Algorithm

III.1 Equivalence classes of Young tableaux

By construction, $\mathrm{Proj}$ is a product of $N_{s}$ operators $\mathrm{Proj}(i)$ ($i=1,\dots,N_{s}$), each of which acts on one of the sites in isolation, i.e., each $\mathrm{Proj}(i)$ is composed of permutation operators that permute only the $m$ particles at site $i$. For example, if $\beta(1)=[2,1]=\beta(2)$, using $\mathrm{Proj}^{[2,1]}(i)=o_{11}^{[2,1]}$ for both sites, we obtain, for site 1,

[TABLE]

and, for site 2,

[TABLE]

The product of these two operators does not ever interchange a particle of site $1$ with a particle of site $2$.

Accordingly, we start by partitioning the $f^{\alpha}$ SYTs of shape $\alpha$ into equivalence classes: Two SYTs are said to belong to the same equivalence class if for each of the $N_{s}$ sites, the $m$ numbers labelling the particles of that site occupy the same $m$ locations on both tableaux. In other words, for each equivalence class, the locations of the blocks for each site are fixed, and SYTs of that class differ only by rearrangements of the particle numbers of the same site within those fixed locations.

For instance, for $\alpha=[4,3,2]$ and $m=3$ particles per site, all of the SYTs belonging to one of the equivalence classes are listed below. For each site $i$ ($i=1,\dots,3$), the three particle numbers $3(i-1)+1$, $3(i-1)+2$, and $3(i-1)+3$ occupy the same locations on all of the SYTs.

[TABLE]

III.2 Basis states and matrix representation

We will use $\mathcal{B}$ to denote the orthonormal basis spanning the projected Hilbert space we are looking for. Every basis state $|\Psi\rangle\in\mathcal{B}$ (sometimes called a multiplet 45) satisfies:

[TABLE]

Since the projection operator does not interchange particles belonging to different sites, it follows that each of these basis states can be identified with a superposition of SYTs that are in the same equivalence class. Indeed, if we generate all of the SYTs belonging to a certain class, we can directly extend Young’s rules [cf. appendix VII.2] to write $\mathrm{Proj}$ as a matrix in terms of this subset of SYTs. The basis state(s) belonging to this class can subsequently be found by solving for $|\Psi\rangle$ in Eq. (20). Repeating this for all of the equivalence classes with shape $\alpha$, we obtain a full basis for the particular subspace of $V^{\alpha}$ associated with the given combination of local irreps $\beta(i)$ $(i=1,\dots,N_{s})$.

The number of basis states depends on the global irrep $\alpha$. Suppose we have local irrep $\beta(i)$ at site $i$ ($i=1,\dots,N_{s})$. Let $D^{\alpha}(\mathcal{B})$ denote the number of states spanning basis $\mathcal{B}$ associated with a given $\alpha$. Then,

[TABLE]

Here, the sum $\sum_{\alpha}$ runs over all shapes $\alpha$ with $mN_{s}$ boxes and at most $N$ rows, and $d_{N}^{\alpha}$ and $d_{N}^{\beta(i)}$ are the dimensions of, respectively, the global irrep $\alpha$ and the local irrep $\beta(i)$ [cf. appendix VII.1].

A major simplification can be made by noting that, given a specific combination of local irreps, only a certain subset of equivalence classes constitute basis states of the basis $\mathcal{B}$ for that combination of irreps. In the most general configuration with irrep $\beta(i)$ at site $i$, these ”viable” classes are those that can be obtained when the Itzykson-Nauenberg46 rules are applied to form the tensor product

[TABLE]

in terms of Young diagrams. These rules can thus be implemented into an iterative scheme to obtain all of the viable classes for a given global irrep $\alpha$ and combination of local irreps $\beta(i)$ $(i=1,\dots,N_{s})$.

Furthermore, we can infer some general consequences of the rules that can be used to immediately identify invalid classes. For example, with irrep $\beta(i)$ at site $i$, a simple condition must always be satisfied: Let $r_{i}$ and $c_{i}$ denote the number of rows and number of columns, respectively, in the shape $\beta(i)$. Then, for any viable class, the $m$ blocks associated with the particles of site $i$ are situated on the tableaux such that there are no more than $c_{i}$ blocks in the same row, and no more than $r_{i}$ blocks in the same column, for all $i=1,\dots,N_{s}$. The SYTs of classes for which this condition does not hold give zero upon projection onto the local irreps $\beta(i)$, and hence do not need to be considered.

As an example, if each site is to be projected onto the irrep $[2,1]$, it follows from this condition that any class of tableaux on which the three blocks for some site are in the same row or in the same column are invalid. In the above diagram, then, tableaux a)-d), and j)-l) can all be neglected for this particular problem. By integrating this constraint with an optimized recursive algorithm, we are able to generate all 867893 of the viable equivalence classes of the global irrep $\alpha=[12,12,12]$ (i.e., the $\mathrm{SU}(3)$ singlet sector for 12 sites) with local irrep $\beta=[2,1]$ at each site in two minutes on single core of a standard CPU.

In the case of symmetric or anti-symmetric irreps at each site, there is only one basis state for each viable equivalence class, and it is possible to project a representative state of each class via a single, pre-determined formula 12. For other symmetries, there may be more than one basis state associated with a given class, and the number of basis states may vary from class to class. In the general case, therefore, we simply solve Eq. (20) for each class as a matrix equation (equivalent to finding the kernel of $(\mathrm{Proj}-\mathbf{I})|\Psi\rangle=\mathbf{0}$, where $\mathbf{I}$ and $\mathbf{0}$ are the identity matrix and the zero matrix).

Since the $\mathrm{Proj}$ does not permute particles between two different sites, the states $|\Psi\rangle$ in Eq. (20) can be calculated more easily by finding the local states $|\Psi_{i}\rangle$ that satisfy

[TABLE]

for each site $i$, then taking the tensor product

[TABLE]

(and renormalizing). In solving equation (21) for $|\Psi_{i}\rangle$, the largest matrix involved has dimension no greater than $m!$.

After finding the states of basis $\mathcal{B}$, it remains to write the interaction Hamiltonian with respect to this basis. This involves using the rules provided in appendix VII.2 to write the matrix coefficients of the permutation appearing in the interaction terms of Eq. 9, then taking matrix-vector products to make the basis transformation. In appendix VII.4, we give a provide a example and detail several useful technical simplifications that may be implemented to optimize the algorithm for this computation.

IV Results

We apply the method outlined above to investigate the presence of edge-states in $\mathrm{SU}(N)$ Heisenberg and AKLT Hamiltonians. We first describe a procedure for building an $\mathrm{SU}(N)$ AKLT-like Hamiltonian, starting a priori with the following irrep at each site:

[TABLE]

where $1\leq p\leq N/2$. $p=1$ corresponds to the adjoint irrep, while $p=N/2$ represents the perfectly rectangular self-conjugate irrep. At each site, the AKLT Hamiltonian favours the virtual decomposition of this local irrep into two fully antisymmetric irreps (i.e., two one-column irreps of lengths $N-p$ and $p$) and their recombination into a singlet over a link of the chain, as depicted in Fig. (3). Consequently, we construct the AKLT Hamiltonian in such a way that it yields the minimal energy (arbitrarily set to zero) for any state whose local wave function over a link of the lattice lives in one of the $p+1$ irreps appearing in the tensor product of the two virtual irreps (the one-column irreps of lengths $N-p$ and $p$). As an example, we consider such a tensor product for $N=6$ and $p=2$:

[TABLE]

Moreover, the AKLT Hamiltonian should give a strictly positive energy for any irrep not present in the right-hand side of Eq. (23) but present in the tensor product of two interacting sites:

[TABLE]

For instance, in the latter decomposition, the Hamiltonian should give strictly positive energy for the irrep $[3,3,2,2,2]$ as well as for all other irreps not listed, but it should give 0 for the first three irreps which, upon removing the first irrelevant column of length $N=6$, correspond to the right hand side of Eq. (23).

To fulfill this set of requirements, one possibility is to take $H_{\text{AKLT}^{N}_{p}}=\sum_{i=1}^{N_{s}-1}H^{\text{AKLT}^{N}_{p}}_{(i,i+1)}$ where the link $(i,i+1)$ Hamiltonian $H^{\text{AKLT}^{N}_{p}}_{(i,i+1)}$ is a product of $p+1$ terms, one per irrep (denoted $\beta_{j}$ below) appearing in the tensor product of the virtual irreps [cf. Eq. 23]:

[TABLE]

$H^{\text{tot}}(i,i+1)$ is the sum of all the transpositions between any pairs of particles living in the ensembles of two sites $i$ and $i+1$:

[TABLE]

In other words, it is the sum of Hamiltonian appearing in Eq. (11) and of the local atomic Hamiltonian for sites $i$ and $i+1$ of Eq. (2) with all coupling constants $J_{k_{i},l_{i}}$ and $J_{i,i+1}$ set to 1. $C^{2}$ is the quadratic Casimir of the $\mathrm{SU}(N)$ irrep.

Usually, a quadratic Casimir for $\mathrm{SU}(N)$ is written as a sum of products of generators of $\mathrm{SU}(N)$, but this approach is inconvenient here and we favour a permutation-like form of the quadratic Casimir, provided in the following equation. For a general irrep $\alpha=[\alpha_{1},\alpha_{2},...,\alpha_{k}]$ with $n$ boxes, the quadratic Casimir can be calculated from the shape $\alpha$ as

[TABLE]

where the $\alpha_{i}$ are the lengths of the rows and the $\alpha^{T}_{j}$ are the lengths of the columns (which are also the rows of the transposed shape $\alpha^{T}$). To apply the last formula in Eq. (25), it is important not to forget to add one column of $N$ boxes to the shapes $\beta_{j}$ appearing in the tensor product in Eq. (23) in order for them to have exactly $m^{2}=N^{2}$ boxes in total. So, for the example with $N=6$ and $p=2$, we obtain:

[TABLE]

where we have set $H^{J\equiv 1}_{(i)}=H^{J\equiv 1}_{(i+1)}=C^{2}([2211])=-5$ since we have the irrep $[2,2,1,1]$ at each site. As a final step, one should check that the quadratic Casimir of all the other irreps (the ones appearing in Eq. (24) but not in Eq. (23)) are strictly larger than the ones appearing in Eq. (23). For the example treated above, it is true. In cases where it is not, to ensure the strict positivity of the last Hamiltonian on the states which would live locally (i.e., along a link) in one of the ”bad” irreps, one could simply square the Hamiltonian (25) (or just part of it).

We have applied this logic to a simpler case that has already attracted some attention in literature, for which $N=4$ and $p=2$. Then, since:

[TABLE]

and

[TABLE]

one can write:

[TABLE]

where we have added the normalization constant $1/56$ so that the linear term has amplitude one, in the development of the last product. A relevant question is whether the $\mathrm{SU}(4)$ Heisenberg Hamiltonian for this irrep (which we can denote $H^{\text{Heis}^{4}_{2}}$) holds the same quantum phase as $H^{\text{AKLT}^{4}_{2}}$. In particular, are the edge states, which are expected to be the states of lowest energy for $H^{\text{AKLT}^{4}_{2}}$, also the lowest energy states for $H^{\text{Heis}^{4}_{2}}$? To answer this question, we have exactly diagonalized the interpolating Hamiltonian

[TABLE]

in each relevant irrep for different chain lengths up to a chain of 10 sites, for which the full Hilbert has a dimension of the order of $10^{13}$. The lowest energies for 10 sites are shown in Fig. (4):

we see that the three lowest energy states are contained in the three irreps appearing in the RHS of Eq. (29), and that they are well separated from the rest of the spectrum. This situation remains true as we tune $\lambda$ from the AKLT point ($\lambda=0$) to the Heisenberg point ($\lambda=1$), demonstrating continuity between the two points.

Interestingly, by calculating the energy spectrum in each global irrep of the Hilbert space, our approach enables us to have a direct characterization of the edge states. In particular, the virtual decomposition of each local irrep $[2,2]$ into two virtual $[1,1]$ irreps leads to the presence of edge states living in those latter irreps, the tensor product of which equal to the sum of the irreps contained in the ground state manifold, according to Eq. (29). These edge irreps $[1,1]$ form an effective and fractionalized representation of the symmetry group: mathematically, it corresponds to some projective representation of the symmetry group. For $\mathrm{SU}(N)$, they can be classified according to the number of boxes (modulo $N$) of the virtual irreps (equal to $2$ in the example above), thus possibly giving rise to $N-1$ non-trivial topological phases27. Those phases (and the corresponding edge states), which have already received some attention in literature, mainly in the two aforementioned cases ($p=1$47; 48; 49; 50 and $p=N/2$51 in the Hamiltonian (25)) are usually characterized through other means due to the difficulty of calculating the $\mathrm{SU}(N)$ quantum numbers when other numerical methods are used. For instance, one can discriminate those topological states by calculating some non-local string order parameter28; 52; 49; 51, or the $Z_{n}$ Berry phase53. Another method makes use of the entanglement spectrum54 whose structure and degeneracies have been shown to reproduce that of the physical edge states 55. In particular, in the case $N=4$ and $p=2$, the degeneracy of the edge irrep $[1,1]$ equal to $6$ has been observed in the entanglement spectrum of an Hamiltonian interpolating between an AKLT point and an Heisenberg point51, consistent with our results. Our method reveals complementary and additional information.

V Conclusion and Perspectives

In this paper, we have shown that the description of the basis states of $\mathrm{SU}(N)$ lattice models in terms of standard Young tableaux, which has previously been implemented for the fundamental representation16 and for totally symmetric or antisymmetric representations12, can be extended to arbitrary irreducible representations. The main difficulty as compared to the case of fully symmetric or antisymmetric irreps lies in the fact that the basis for each equivalence class of Young tableaux is no longer of dimension 1, and special emphasis has been placed on the construction of this basis using projection operators defined in terms of elementary permutations. As in previous cases, this approach allows one to work directly in irreps of the global $\mathrm{SU}(N)$ symmetry. Since for antiferromagnetic interactions the low lying states are typically found in irreps whose dimensions are much, much smaller than that of the full Hilbert space, this enables us to reach sizes that are inaccessible to other formulations.

We have applied the method to the study of $\mathrm{SU}(4)$ Heisenberg chain with the irrep $[2,2]$ at each site, demonstrating that the ground state is of VBS type, and explicitly characterizing the irreps of the edge states, a very interesting by-product of this approach.

In the future, we aim to further develop our technique in three directions. First, we hope to improve its efficiency in treating longer range interactions as well as three-site ring exchange interactions, in order to investigate 2D systems with exotic properties. With one particle per site, $\mathrm{SU}(N)$ Heisenberg models on 2D systems have already been shown to possibly host chiral spin liquids with chiral edge states described by the $\mathrm{SU}(N)_{1}$ Wess-Zumino-Novikov-Witten conformal field theory 56; 57; 43. To create non-Abelian Chiral Spin Liquids described by the $\mathrm{SU}(N)_{k}$ WZNW CFT (with $k>1$) requires more sophisticated irreps at each site of a 2D lattice 17; 58; 59. Secondly, we would like to generalize the use of SYTs to DMRG in order to take advantage of the $\mathrm{SU}(N)$ symmetry, an issue of considerable interest45; 60; 61; 62 at present. Finally, it would be interesting to determine whether it is possible to integrate both the complete $\mathrm{SU}(N)$ group and spatial symmetries in the same algorithm, as has been done for the special case of $N=2$ 63; 64.

VI Acknowledgements

We thank Akira Furusaki, Philippe Lecheminant, and Thomas Quella for useful discussions. This work has been supported by the Swiss National Science Foundation, and by the Research Internship program at EPFL.

VII Appendix

VII.1 Irreps of $\mathrm{SU}(N)$

In general, for a system of $n$ particles, each irrep of $\mathrm{SU}(N)$ can be associated with a Young diagram composed of $n$ boxes arranged in at most $N$ rows. This represents a particular set of $\mathrm{SU}(N)$-symmetric $n$-particle wave functions. The shape $\alpha$ of the Young diagram is specified by a partition $\alpha=[\alpha_{1},\alpha_{2},\dots,\alpha_{k}]$ (with $1\leq k\leq N$ and $\sum_{j=1}^{k}\alpha_{j}=n$), where the row lengths $\alpha_{j}$ satisfy $\alpha_{1}\geq\alpha_{2}\geq\dots\geq\alpha_{k}\geq 1$. The diagram can be filled with numbers 1 to $n$, and the resultant tableau is said to be standard if the entries are increasing from left to right in every row and from top to bottom in every column. Standard Young tableaux (SYTs) play a central role in representation theory.

Using $[1]=\Box$ to denote the fundamental irrep, the set of all $n$-particle wave functions live in the full Hilbert space $\Box^{\otimes n}$. The multiplicity $f^{\alpha}$ of irrep $\alpha$ in this space is equal to number of SYTs with shape $\alpha$. This number can be calculated from the hook length formula,

[TABLE]

where the hook length $l_{i}$ of the $i$-th box is defined as the number of boxes to the right of it in the same row, plus the number of boxes below it in the same column, plus one (for the box itself). The dimension $d_{N}^{\alpha}$ of the irrep can also be calculated from the shape as

[TABLE]

where $\gamma_{i}$ is the algebraic distance from the $i$-th box to the main diagonal, counted positively (resp. negatively) for each box above (resp. below) the diagonal. The full Hilbert space can be decomposed as

[TABLE]

where $V^{\alpha}$ is the sector corresponding to irrep $\alpha$ (and if $d_{N}^{\alpha}>1$, $V^{\alpha}$ can itself be decomposed into $d_{N}^{\alpha}$ equivalent subsectors, $V^{\alpha}=\oplus_{i=1}^{d^{\alpha}_{N}}V^{\alpha}_{i}$). The equation for the dimension of the full Hilbert space thus reads

[TABLE]

where the sum runs over all Young diagrams with $n$ boxes and no more than $N$ rows.

As an example, for $n=2$, $\Box^{2}$ can be decomposed as a sum of the subspace spanned by the symmetric two-particle wave functions and that spanned by the antisymmetric two-particle wave functions:

[TABLE]

There is only one SYT for each of the diagrams [2] and [1,1], so $f^{[2]}=f^{[1,1]}=1$, while Eq. (34) gives $d_{N}^{[2]}=\frac{N(N+1)}{2}$ and $d_{N}^{[1,1]}=\frac{N(N-1)}{2}$. It is easy to check that Eq. (36) is satisfied.

VII.2 Young’s orthogonal representation of the symmetric group

This subsection summarizes some useful results concerning the orthogonal representation of the symmetric group.

For a given tableau shape $\alpha$, a convenient representation of the symmetric group $\mathcal{S}_{n}$ can be formulated using Young’s orthogonal units $\{o_{rs}^{\alpha}\}_{r,s=1\dots f^{\alpha}}$. These are specific linear combinations of permutations, whose explicit forms are provided in section VII.3 (for $n=3$ and $n=4$). They satisfy orthonormality:

[TABLE]

as well as completeness:

[TABLE]

and form a basis in which any linear superposition $\eta$ of permutations belonging to $\mathcal{S}_{n}$ can be uniquely decomposed as

[TABLE]

where $\mu_{rs}^{\alpha}(\eta)$ are real coefficients.

An important result we will make frequent use of is that successive transpositions $P_{k,k+1}$, i.e., permutations between the consecutive numbers $k$ and $k+1$ ($1\leq k\leq n-1$), takes an extremely simple form in the basis of orthogonal units. If we write $P_{k,k+1}=\sum_{\alpha,t,q}\mu_{tq}^{\alpha}(P_{k,k+1})o_{tq}^{\alpha}$, then, for a given shape $\beta$, the matrices $\bar{\bar{\mu}}^{\beta}(P_{k,k+1})$ defined by

[TABLE]

are symmetric and orthogonal, and very sparse, with at most two nonzero entries in each row and in each column. These entries can be explicitly calculated as follows. We assign some fixed order (named last letter order sequence) to the $f^{\alpha}$ SYTs and label them $S_{1},\dots,S_{f^{\beta}}$. If $k$ and $k+1$ are in the same row (resp. column) on the tableau $S_{t}$, then $\mu_{tt}^{\beta}(P_{k,k+1})=+1$ (resp. $-1$), and all other matrix elements involving $t$ vanish. If $k$ and $k+1$ are not in the same column nor the same row on $S_{t}$, and if $S_{q}$ is the tableau obtained from $S_{t}$ by interchanging $k$ and $k+1$, then the only non-vanishing matrix elements involving $t$ or $u$ are given by

[TABLE]

Here, $\rho$ is the inverse of the axial distance from $k$ to $k+1$ on $S_{t}$, which is computed by counting $+1$ (resp. $-1$) for each step made downward or to the left (resp. upward or to the right) to reach $k+1$ from $k$.

This simple yet incredibly useful formula is clarified in Figure 5. Moreover, since every permutation can be factorized into successive transpositions, we can use these rules to write the exact matrix representation of any permutation or linear superposition thereof via a few elementary calculations.

VII.3 Explicit form of the orthogonal units

We provide below the explicit expression of some orthogonal units for the irreps $[2,1]$ and $[2,2]$, as examples.

To construct the orthogonal units $\left\{o^{\alpha}_{rs}\right\},\forall r,s=1\dots f^{\alpha},$ for a given shape $\alpha$ with $n$ boxes, one can use the procedure of Thrall65, which involves recursively constructing a product of antisymmetrizers and symmetrizers, adding one box at a time, for a total of $n$ boxes.

We provide an alternative method that is useful when $n$ is not too large (typically, for $n\lesssim 10$). For a given shape $\alpha$, we first generate all of the SYTs (using either the algorithm NEXYTB from Chap. 14 of Ref. (66) or the algorithm outlined in Appendix 5 of Ref. (12)). We then compute the $n!$ matrices of size $f^{\alpha}\times f^{\alpha}$ that represent the permutations $\eta$ of the symmetric group $\mathcal{S}_{n}$. The rules given in VII.2 can be used to write the $n-1$ matrices for transpositions between consecutive numbers, $P_{k,k+1}$ ($1\leq k\leq n-1$), and the remaining permutations of $\mathcal{S}_{n}$ can be obtained by decomposing them into transpositions and calculating the matrix products. Finally, we use the relation in Ref. (67), which is the inverse of that shown in Eq. (38):

[TABLE]

where $\mu^{\alpha}_{sr}(\eta^{-1})$ is the element in the $s^{th}$ row and $r^{th}$ column of the inverse of the matrix corresponding to permutation $\eta$.

The four orthogonal units associated for the shape [2,1] are thus linear combinations of the $3!=6$ permutations of $\mathcal{S}_{3}$:

[TABLE]

while for the shape $[2,2]$, for instance, $o_{11}^{[2,1]}$ (which may directly be used as the projection operator onto the irrep associated with $[2,2]$) is

[TABLE]

VII.4 Matrix representation of the interaction Hamiltonian

We discuss below several technical simplifications that may be implemented to optimize the construction of the matrix representing the permutation Hamiltonian. We will focus on the specific example of the SU($3$) open chain of length $N_{s}=5$ sites with local irrep $\beta=[2,1]$ at every site. We will concentrate on the singlet sector (i.e of global irrep $[5,5,5]$). There are 16 viable equivalence classes, listed below. It will be found that each class may be spanned by one, two, or eight states, for a total of 32 basis states in $\mathcal{B}$.

[TABLE]

For instance, to write the the matrix of the Hamiltonian for the interaction between site 2 (particles 4, 5, and 6) and 3 (particles 7, 8, and 9),

[TABLE]

[cf. Eq. (9)], we consider, for each class, only the pink and grey blocks.

For class a), the set of possible standard subtableaux for site 2  $\left(\ydiagram{1+2,1}\,\right)$ is

[TABLE]

We write $\mathrm{Proj}^{[2,1]}(2)=\frac{1}{6}(2Id+2P_{45}-P_{46}-P_{56}-P_{45}P_{46}-P_{46}P_{45})$ with respect to these three subtableaux using Young’s orthogonal representation, forming a $3\times 3$ matrix $\mathbf{Proj}^{[2,1]}(2)$.

Computing the kernel of $\mathbf{Proj}^{[2,1]}(2)-\mathbf{I}$ yields a single nonzero 3-dimensional vector:

[TABLE]

Performing a similar routine for site 3, with standard subtableaux $S_{3}=\left\{\ytableaushort{\none\none\none\none 7,\none\none 8,9}\,,\,\ytableaushort{\none\none\none\none 7,\none\none 9,8}\,,\,\ytableaushort{\none\none\none\none 8,\none\none 7,9}\,,\,\right.$

$\left.\ytableaushort{\none\none\none\none 8,\none\none 9,7}\,,\,\ytableaushort{\none\none\none\none 9,\none\none 7,8}\,,\,\ytableaushort{\none\none\none\none 9,\none\none 8,7}\right\}$ yields two basis states represented as 6-dimensional vectors:

[TABLE]

We then take the tensor products: Let

[TABLE]

and

[TABLE]

These are the vectors we need to work with to write $H_{(2,3)}$ with respect to our basis. They each have 18 entries, corresponding to the set of 18 subtableaux formed by combining those of sites 2 and 3:

[TABLE]

In a similar fashion, the vectors representing the states $|\Psi_{2\otimes 3}\rangle$ can be found for the other classes. Each will be at most $6!^{2}=36$-dimensional.

Then, the matrix elements of $\mathbf{H}_{(2,3)}$ in our 32-dimensional basis can be found the following way: We first write $H_{(2,3)}$ in the basis of subtableaux (please refer to the second bullet below) and then take the products

[TABLE]

where $x$ and $y$ run over all of the equivalence classes.

The computational time and memory required can be reduced (by several orders of magnitude, for large systems) by taking advantage of the following key observations:

Certain matrix elements are zero and need not be calculated. (In fact, the matrix representations for the $H_{(i,j)}$ in basis $\mathcal{B}$ are extremely sparse.)

Condition 1: In order for Eq. (43) to be nonzero, the two classes $x$ and $y$ must be such that 6 blocks corresponding to sites 2 and 3 are in the same 6 locations on both classes. This follows from the fact that $H_{(2,3)}$ consists only of permutations between sites 2 and 3. Therefore, products between states of between class a) and those of class d) would automatically be zero, since the overall locations of blocks for sites 2 and 3 on class a), \ydiagram2+3,1+2,1, are different from those on class d), $\ydiagram{2+2,1+2,2}$. It can be seen that only classes g), j), n), and p) have the same total shape for sites 2 and 3 as class d). Additionally, between two classes that give the same subtableaux shape, there must be at most one interchange between the three blocks of site 2 and those of site 3. For example, the subtableaux of class d), \ydiagram[(lightgrey)]0, 2+1, 2[(lightpink)]2+2, 1+1, and of class g), \ydiagram[(lightgrey)]0, 1+2,1+1*[*(lightpink)]2+2,0,1, differ only by an interchange of the leftmost pink block with the leftmost grey block. However, the subtableaux of class p), $\ydiagram[*(lightgrey)]{2+2,2+1}*[*(lightpink)]{0,1+1,2}$, differs from that of classes d) as well as g), j), and n) by two interchanges. Hence, Eq. (43) does not need to be calculated between the states of class p) and those of any other class.

Condition 2: Moreover, classes $x$ and $y$ must be such that the locations of the blocks for each and every one of the other sites are identical. To illustrate this more clearly, we consider a few classes of a slightly larger shape:

[TABLE]

Here, the blocks of sites 2 and 3 on all of above classes satisfy Condition 1. However, the blocks for site 5 and for site 6 are situated differently on class ii) than on classes i) and iii). As a result, products analagous to that in Equation (43) between class ii) and class i) or between class ii) and iii) will all be zero.

These two conditions follow from the nature of the interaction Hamiltonian and the orthogonality of the basis states in $\mathcal{B}$. 2. 2.

Two vectors $|\Psi^{x}_{2\otimes 3}\rangle$, $|\Psi^{y}_{2\otimes 3}\rangle$ from different classes, i.e., $x\neq y$, are with respect to completely different bases. For class a), we found the set $S_{2}\otimes S_{3}$ [Eq. (42)], which includes 18 subtableaux, but the analogue for class e) consists of 9 completely different subtableaux:

[TABLE]

Then, to calculate $\langle\Psi^{e}_{2\otimes 3,1}|H_{(2,3)}|\Psi^{a}_{2\otimes 3,1}\rangle$ and $\langle\Psi^{e}_{2\otimes 3,1}|H_{(2,3)}|\Psi^{a}_{2\otimes 3,2}\rangle$ (there are two vectors for class a), as found above, and one for class e)), we first observe that $H_{(2,3)}$ is a sum of permutations between particles of sites 2 and 3, each of which can be decomposed into a product of transpositions with $P_{67}$ in the middle. For instance, $P_{69}=P_{79}P_{67}P_{79}$ and $P_{59}=P_{56}P_{79}P_{67}P_{79}P_{56}$.

Since matrix-vector products are linear, we take a single permutation in the sum at a time. Consider the product $\langle\Psi^{e}_{2\otimes 3,1}|P_{59}|\Psi^{a}_{2\otimes 3,2}\rangle=\langle\Psi^{e}_{2\otimes 3,1}|P_{56}P_{79}P_{67}P_{79}P_{56}|\Psi^{a}_{2\otimes 3,2}\rangle$. We find matrices for $P_{56}$ and $P_{79}$ with respect to each of the sets of subtableaux $S^{e}_{2}\otimes S^{e}_{3}$ and $S^{a}_{2}\otimes S^{a}_{3}$, by applying Young’s rules. Then, instead of writing a full matrix representation for $P_{67}$, we place $S^{a}_{2}\otimes S^{a}_{3}$ on the columns and $S^{e}_{2}\otimes S^{e}_{3}$ on the rows (see below) to find a $9\times 18$ matrix that can be used in the above product.

[TABLE] 3. 3.

For the interaction Hamiltonian between two given sites, we only need to work with the subtableaux composed of the $2m$ blocks of those two sites. Some of these sets of subtableaux will be the same for multiple classes. For example, among the 20 classes in our example, the blocks for site 1 and site 2 are the same on classes a) through d): \ydiagram[(lightblue)]2,1[*(lightpink)]2+2,1+1. Hence, when calculating $H_{(1,2)}$, it suffices to find the vectors and resultant matrix elements for such a configuration just once, instead of repeating the computation. 4. 4.

Finding the kernel of $\mathbf{Proj}(i)-\mathbf{I}$ for each site $i$ individually and forming the tensor product(s) between the resultant vectors of two sites may not account for all of the basis states, since, as we have seen, some sites give rise to more than one vector. However, when writing the matrix for $H_{i,i+1}$, it suffices to find the vectors for the sites $i$ and $i+1$ and repeating their matrix elements in certain positions in the matrix (so that the same ordered basis is used for all the $H_{i,i+1}$). These positions depend on the number of vectors that correspond to the sites $j$ for $j<i$ and for $j>i+1$.

It should be noted that the list of simplifications provided above is vital, but by no means exhaustive.