# Hidden Bethe states in a partially integrable model

Zhao Zhang<EMAIL_ADDRESS>SISSA and INFN, Sezione di Trieste, via Bonomea
265, I-34136, Trieste, Italy  Tsung-Dao Lee Institute, Shanghai Jiao Tong
University, Shanghai 200240, China  Giuseppe Mussardo<EMAIL_ADDRESS>SISSA
and INFN, Sezione di Trieste, via Bonomea 265, I-34136, Trieste, Italy

###### Abstract

We present a one-dimensional multi-component model, known to be partially
integrable when restricted to the subspaces made of only two components. By
constructing fully anti-symmetrized bases, we find integrable excited
eigenstates corresponding to the totally anti-symmetric irreducible
representation of the permutation operator in the otherwise non-integrable
subspaces. We establish rigorously the breakdown of integrability in those
subspaces by showing explicitly the violation of the Yang-Baxter’s equation.
We further solve the constraints from Yang-Baxter’s equation to find
exceptional momenta that allows Bethe Ansatz solutions of solitonic bound
states. These integrable eigenstates have distinct dynamical consequence from
the embedded integrable subspaces previously known, as they do not span their
separate Krylov subspaces, and a generic initial state can partly overlap with
them and therefore have slow thermalization. However, this novel form of weak
ergodicity breaking contrasts that of quantum many-body scars in that the
integrable eigenstates involved do not have necessarily low entanglement. Our
approach provides a complementary route to arrive at quantum many-body scars
since, instead of solving towers of single mode excited states based on a
solvable ground state in a non-integrable model, we identify the integrable
eigenstates that survive in a deformation of the Hamiltonian away from its
integrable point.

##  I Introduction

Recent progress in cold atom experiments have made possible the simulation of
a larger variety of strongly correlated systems [1]. In particular, alkine-
earth atoms such as 137Yb and 87Sr in optical atoms can be used to realize the
Mott insulating phase of Fermi-Hubbard model, which can be described by the
SU(N) Heisenberg Hamiltonian [2, 3, 4, 5, 6]. This multi-component
generalization of the spin-$\frac{1}{2}$ Heisenberg model was solved based on
Yang’s generalization [7, 8] of the Bethe Ansatz method developed by
Sutherland [9]. Using the same nested Bethe Ansatz method, Babelon, de Vega,
and Viallet constructed an integrable $Z_{n+1}\times Z_{n+1}$ symmetric
generalization of the XXZ model, from the solution of Yang-Baxter’s equations
[10]. However, the anisotropy in this model was artificially designed so that
the model is exactly solvable but not realistic to be realized in experiment.

In this paper, we propose a minimal generalization of the spin-$\frac{1}{2}$
XXZ Hamiltonian, written in terms of permutation operators, whose SU(N)
symmetry is explicitly broken to the symmetric group $\mathfrak{S}_{N}$ by a
diagonal term composed of the Cartan operators. As shown in Sutherland’s
original formulation of the isotropic model [9], such a multi-component spin
model can be equivalently formulated in the language of one-dimensional hard-
core bosonic particles with density-density interactions. The model proposed
in this paper is in part motivated by a recent quantum system bearing a number
theory analogy [11], whose low energy effective theory can be seen as the
current model in certain limits of its parameters. Some of the questions
raised in the paper [11], such as the the scaling of spectral gap from ground
state, are partially answered in this paper.

It is worth to underline that, while the model proposed here is in general not
integrable, the Krylov subspaces spanned by configurations involving only two
components are nonetheless integrable, as was originally noted in a series of
papers in the past [12, 13, 14]. Hence, such a model falls in the category of
quasi-exactly solvable models, which has been studied since long in quantum
mechanics [15]. In many-body systems, there is a class of models called
frustration-free, whose zero energy ground state can be explicitly written as
the common lowest energy eigenstate of each local operator in the Hamiltonian.
Examples include the Majumdar-Ghosh model [16], AKLT model [17], Motzkin and
Fredkin spin chains [18, 19, 20, 21, 22, 23] and, in two dimensions, Kitaev’s
toric code model [24], to name a few. Many of these models boast solvable
eigenstates beyond ground state as well. For instance, stabilizer code
Hamiltonians consist of mutually commuting operators, so the entire spectrum
is solvable. Despite the fact that neither Majumdar-Ghosh [25], nor spin-1
AKLT model [26] satisfy this condition, they are known to have exact excited
states. Recently these models have attracted a lot of attention for the
experimental observation of anomalous dynamics in close related Rydberg atom
experiments [27] that interpolates between ergodicity breaking and
thermalization behaviours. Such a dynamics has been attributed to a large
overlap of the initial state with a tower of numerically solved eigenstates
with equally spaced energies that violate the eigenstate thermalization
hypothesis (ETH) in the so-called PXP Hamiltonian [28, 29], used to model the
Rydberg blockade experiment [30]. The discovery of this new form of weak
ergodicity breaking has drawn more attention to analytical efforts for solving
exact excited states in non-integrable models [31, 32, 33, 34, 35, 36, 37, 38,
39, 40, 41, 42, 43]. The tremendous success achieved in searching these
quantum many-body scar (QMBS) states is largely due to the approach which
examines the sparsity of entanglement spectrum for numerically exactly
diagonalized states [32]. In this way, these states can be analytically
represented in terms of matrix product states [36] making use of a reverse
engineering. Furthermore, a large portion of these exact eigenstates share the
common features of having energy expressed in terms of either integer or
rational numbers, and they can be identified using single mode approximation
(SMA) with $\pi$-momenta [36]. Some of these features observed in individual
examples have found a more general explanation in the unified framework of
spectrum generating algebra [44]. Concerning the partially integrability of
the model proposed in this paper, it is useful to say that it does not lead to
the slow dynamics of weak ergodicity breaking, as the integrable subspaces
studied before [12, 13, 14]. Indeed, depending on the subspace to which an
initial state belongs to, the system either thermalizes according to ETH, or
exhibits strong ETH-violation due to partial integrability and behaves as a
generalized Gibbs ensemble.

The plethora of frustration-free models (and their success in providing
insights in understanding various aspects of many-body systems and quantum
computing through rigorous mathematical theorems [45, 46, 47, 48, 49, 50]),
gives naturally rises to the question: to what extent can a frustration-free
Hamiltonian interpolate between the situation in which only the ground state
is known, and a stabilizer code that provides the entire solvable spectrum?
The examples of exact QMBS in AKLT models partially answered this question, in
that the additional eigenstates are of momenta $\pi$, allowing the diffractive
scattering to be canceled out in a momentum eigenstate [32]. In this work, we
will extend this idea to non-diffractive scattering with other momenta using
Bethe Ansatz, in a non-integrable multi-component model with a “weaker”
frustration. Such a strategy of easing frustration by enlarging local degrees
of freedom has been employed by models with SU(N)-singlet simplex solid ground
states [51], in the form of geometric frustration, which has also one-
dimensional cousin models with longer range interaction, such as SU(3) spin
chains with trimer and valance bond solid (VBS) ground states [52], and SU(N)
VBS [53]. Here, instead of designing a Hamiltonian as sum of projectors
pinning down a desired ground state, we work with a minimal multi-component
generalization to the spin-$\frac{1}{2}$ XXZ spin chain with 2-local
interaction which violates the Yang-Baxter’s equation (YBE) and therefore is
not integrable. We show that the reduced frustration with larger local Hilbert
space allows us to have additional Bethe Ansatz integrable eigenstates in any
generic non-integrable Krylov subspace. Unlike the partial integrability
studied before, these integrable eigenstates now bear a new form of weak
ergodicity breaking, as a generic initial state in these subspaces will have
overlap with both integrable and non-integrable eigenstates.

Entanglement entropy and its scaling has been a major tool in probing many-
body systems with strong correlations, see [54] for a review. For one-
dimensional systems, Hastings have given a rigorous proof that the ground
state entanglement entropy scaling with system size is bounded by a constant
for gapped systems [55]. Meanwhile, people have been searching for area-law
violating grounds states in gapless systems with logarithmic [18, 21], power
law [19], and even linear scaling [20, 22, 23, 56, 57]. In the study of
quantum many-body scars, the sparsity of entanglement spectrum and the sub-
volume law scaling of entanglement entropy has been used to show that the
exact excited states violate strong eigenstate thermalization hypothesis [32].
With the enlarged local Hilbert space, the entanglement entropy of the ground
state of our model scales not only with the system size, but with the
dimensionality of local Hilbert space as well. By rigorously establishing a
volume law scaling of the ground state in the extreme frustration-free case of
number of components equal to the system size (which also accounts for the
additional contribution to the entanglement of integrable excited states from
the antisymmetrized basis in addition to that from the corresponding Bethe
Ansatz states in the spin-$\frac{1}{2}$ XXZ model[58, 59]), we show rigorously
that the entanglement entropy scaling of our integrable eigenstates
interpolates between area-law (which violates entanglement entropy of
corresponding excited states in XXZ chain) and volume law scaling of the
totally antisymmetric bases.

The paper is organized as follows. In Sec. II, we introduce the model and
discuss its symmetries, which we will use later for achieving the
diagonalization of its Hamiltonian. Sec. III is rather articulated: in III.1
we first construct a proper basis to diagonalize the integrable subspaces of
the Hilbert space; in III.2 we illustrate the previous procedure in terms of
an elementary approach for diagonalizing the two-body case; in III.3, using
the general coordinate Bethe Ansatz, we establish the embedding of the
spectrum of spin-$\frac{1}{2}$ Heisenberg in our model. In Sec. IV, we
rigorously prove the that a generic sector with multiple identical particles
is non-integrable unless the anisotropy is turned off, and we use this result
to shed more light on the reason of the partial integrability of the sectors
discussed in the previous section. We also provide an explanation to the
$\pi$-momenta QMBSs emerging from the condition to satisfy YBE with a generic
anisotropy. In Sect. V, we address the case of a large number of components
(this makes the ground state frustration free) and we show analytically that
the entanglement entropy of the ground state scales linearly in the
thermodynamic limit. We also derive a rigorous decomposition of entanglement
entropy contribution from the Bethe Ansatz wave-function and the
antisymmetrized bases. In Sec. VI, we draw our conclusions, also pointing out
few interesting directions worth exploring in the future. The paper has also
three appendices which address particular technical points discussed in the
main text.

##  II  The Hamiltonian and its symmetry

The local Hilbert space of our model $\mathbb{C}^{N}$ is spanned by the $N$
components of the SU(N) group, or $N$ species of hard-core bosons. The model
is defined on a one-dimensional lattice of length $L$, with periodic boundary
condition $\mathbb{C}_{L+1}\equiv\mathbb{C}_{1}$. Its Hamiltonian takes the
form

$H=\sum_{i=1}^{L}\left[P_{i,i+1}+2\Delta C_{i,i+1}\right]$ (1)

where transposition operator $P_{i,i+1}$ and the diagonal anisotropy term
$C_{i,i+1}$ are defined as

$\displaystyle P_{i,i+1}=$ $\displaystyle\sum_{a,b=1}^{N}e_{i}^{(ab)}\otimes
e_{i+1}^{(ba)},$ (2) $\displaystyle C_{i,i+1}=$
$\displaystyle\sum_{a=1}^{N}e_{i}^{(aa)}\otimes e_{i+1}^{(aa)},$ (3)

where $(e^{(ab)})_{cd}=\delta_{c}^{a}\delta_{d}^{b}$ is the standard basis of
$N\times N$ matrices, such that $P_{i,i+1}(v_{i}\otimes
v_{i+1})=v_{i+1}\otimes v_{i}$ for any $v_{i},v_{i+1}$. The transposition term
is SU(N) invariant, as it can be written as traces over product of SU(N)
generators. The diagonal term however, can be written only in terms of the
generators of the Cartan subalgebra of SU(N), and therefore breaks the
symmetry down to the symmetric group $\mathfrak{S}_{N}$.

Apart from the multi-component chain interpretation, one can also look at this
model as a chain of hard-core bosons with a local color degree of freedom
subjected to a repulsive density-density interaction. Henceforward, we shall
interchange freely between these two languages for the convenience of
presentation. The formulation in terms of hard-core bosons permits to put in
correspondence the present model with the number-theory quantum model analyzed
in [11] (the so-called coprime spin ladder model), in the sense that each
boson with different color can be put in correspondence with a prime number
while the diagonal term in the Hamiltonian can be associated to the ”coprime
interaction”. There are though some differences between the two models: unlike
the one discussed in [11], we do not have here a quasi-2D lattice and,
moreover, in the present case each boson does not have a composite structure
consisting of more elementary degrees of freedom. Moreover, there is no onsite
dynamical terms in the Hamiltonian that allows bosons of one species to
transmute into others. Nonetheless, the coprime spin ladder model does
reproduce the Hamiltonian (1) when restricted to a certain sector of particle
content $\\{c_{1}^{n_{1}},c_{2}^{n_{2}},\cdots,c_{s}^{n_{s}}\\}$, where
$c_{A}$ denotes the particle of color $A$, and $n_{A}$ denotes the number of
times it appear along the chain. Therefore, our Hilbert space is fragmented
into Krylov sub-spaces specified by the particle content. In each sector, the
restricted Hamiltonian will be equivalent to [9]

$H(\\{c_{1}^{n_{1}},c_{2}^{n_{2}},\cdots,c_{s}^{n_{s}}\\})=\sum_{i=1}^{L}\left[(1+2\Delta)\sum_{A}N_{i,i+1}^{(AA)}+\sum_{A<B}P_{i,i+1}^{(AB)}\right],$
(4)

where $N^{(AA)}$ counts the number of neighboring pairs of species $A$, and
$P^{(AB)}$ transposes only neighboring pairs of species $A$ and $B$.

This restricted Hamiltonian manifests permutation symmetries in addition to
the translational invariance. Firstly, the total number of each species
$n_{A}$ is conserved. Secondly, being each particle indistinguishable, the
Hamiltonian is invariant under
$\mathfrak{S}_{n_{1}}\otimes\mathfrak{S}_{n_{2}}\otimes\cdots\otimes\mathfrak{S}_{n_{s}}$,
where $\sum_{A=1}^{s}n_{A}=L$, and $\otimes$ denote the outer product between
symmetric groups [60, 61]. Notice that this is a subgroup of the original
$\mathfrak{S}_{L}$ symmetry. Therefore, eigenvectors of this Hamiltonian can
be identified within each of its irreducible representations (irrep’s). This
step would already reduce the complexity of the problem significantly enough
so that one could diagonalize the Hamiltonian written in terms of standard
irreducible representations matrices by brute force, as was done in [62, 63].
With some more thoughts, Gaudin [64, 65] turned the Young tableaux formalism
(equivalently described in terms of Hund’s method, usually more appealing to
physicist), into a linear system of constraints from Fock’s conditions, which
by one remarkable algebraic identity after another, arrived at the same Bethe-
Yang Hypothesis states. However, the anisotropy term in our Hamiltonian mixes
irreducible representations, making this approach not applicable. So, in our
case it is more promising to follow the strategy of Yang and Sutherland [7, 8,
9] and first take advantage of the translational invariance to go to momentum
space and consider irrep’s of permutation operators.

##  III  Diagonalization of the integrable subspaces

The Hamiltonian (1) without the anisotropy term was diagonalized by Sutherland
using nested Bethe Ansatz [9] 111Although Sutherland was actually solving the
ferromagnetic model for fermionic particles, he showed the spectrum is that of
the antiferromagnetic bosonic one in reversed ordering.. As shown in the next
section, the anisotropy term breaks the integrability of the original
Hamiltonin by violating the Yang-Baxter Equation (YBE) of scattering matrices,
except for certain irrep’s. The large number of local degrees of freedom, or
species of particles, can be both a curse and a blessing. In this paper, we
take advantage of this property by constructing a particular basis, in which
the solution of relatively low-lying eigenstates using the Bethe Ansatz is
drastically simplified. In this section, we first show how this approach helps
solving eigenstates in the sector where only one specie of particle appear
multiple times, while the rest of the species each appears only once. In the
next section, we show that the YBE is still violated in this basis for a
generic sector where multiple species appear more than once, as well as in the
orthogonal subspace to the basis within the same symmetry sector.

###  III.1  Maximally antisymmetrized basis

The key observation which leads to use this method consists of realizing that
each transposition operator is minimized by the eigenvalue $-1$. So the ground
state corresponds to that state reached by antisymmetrizing as many neighbors
as possible. Of course this can only reconcile among different species, as
antisymmetrizing a symmetric pair gives $0$. Therefore, we can restrict our
searching for the ground state in the subspace of the Hilbert space spanned by
the basis

$|{i_{1},i_{2},...,i_{n}}\rangle=\sum_{\sigma\in\mathfrak{S}_{L-n}}\text{sgn}(\sigma)(-)^{\sum_{a=1}^{n}i_{a}}|{i_{1},i_{2},...,i_{n}};\sigma\left(c_{2}~{}c_{3}~{}\cdots~{}c_{L-n}\right)\rangle,$
(5)

where $i_{1}<i_{2}<\cdots<i_{n}$ labels the location of the $n$ identical
species (using the correspondence with the coprime spin ladder model, these
species are chosen to be the one corresponding to the smallest prime numbers).
The symbol sgn($\sigma$) denotes the signature of the permutation $\sigma$
acting on the reference state in a certain order of species with $c_{j}$
denoting particle of the $j$-th color, before permuting the identical species
to their final locations. An illustration of the $n=2,L=5$ case is given in
Fig. 1.

Figure 1: An example of the basis vectors in a sector with particles of $4$
different colors, the red onappearing twice.

The reason for this convenient choice of basis is such that

$P_{j,j+1}|{i_{1},...,i_{n}}\rangle=\begin{cases}-|{i_{1},...,i_{n}}\rangle,&j\neq
i_{a},i_{a}-1,\forall a\\\
-|{i_{1},...,i_{a-1},i_{a}+1,i_{a+1},...,i_{n}}\rangle,&j=i_{a}\neq
i_{a+1}-1\\\ |{i_{1},...,i_{n}}\rangle,&j=i_{a}=i_{a+1}-1\end{cases}$ (6)

for any $i_{a}<L$. Since the transposition operators act the kinetic term of
the Hamiltonian, antisymmetrizing the species is equivalent to picking a
subspace where all the particles with different species except the identical
one that appears multiple times to have momenta $\pi$. Periodic boundary
condition on the physical chain leads to (anti-)periodic boundary condition on
our artificial basis, depending on the parity of the length of the chain.

$|i_{1},...,i_{n-1},L+1\rangle\equiv(-1)^{L}|1,i_{1},...,i_{n-1}\rangle,$ (7)

where $i_{1}>1$. For simplicity, we restrict our discussion to the case of
even length from now on. The odd length case can be treated likewise. In this
basis, we can diagonalize the Hamiltonian by expanding the eigenvector as

$|v\rangle=\sum_{1\leq i_{1}<\cdots<i_{n}\leq
L}w_{i_{1},...,i_{n}}|{i_{1},i_{2},...,i_{n}}\rangle.$ (8)

Keeping in mind that our basis vectors are already eigenvectors of
transposition operators unless it acts on the identical species, the
eigenvalue equation

$H|v\rangle=E_{v}|v\rangle.$ (9)

results in a difference equation version of the Helmholtz equation of weights
$w_{i_{1},...,i_{n}}$

$\nabla^{2}w_{i_{1},...,i_{n}}=-(E_{v}+L)w_{i_{1},...,i_{n}},$ (10)

where we used $\nabla^{2}$ to denote the discrete Laplacian operator so as to
avoid confusion with the anisotropy parameter. If none of the position $i$’s
is consecutive,

$\nabla^{2}w_{i_{1},...,i_{n}}\equiv\sum_{a=1}^{n}\left(w_{i_{1},...,i_{a}-1,...,i_{n}}+w_{i_{1},...,i_{a}+1,...,i_{n}}-2w_{i_{1},...,i_{n}}\right),$
(11)

where $w_{i_{1},...,i_{n-1},L+1}\equiv(-1)^{n-1}w_{1,i_{1},...,i_{n-1}}$ due
to (7). When certain position $i$’s form consecutive strings
$l_{a},l_{a}+1,...,r_{a}$, $\nabla^{2}$ becomes instead

$\begin{split}\nabla^{2}w_{l_{1},...,r_{1},...,l_{s},...,r_{s}}=\sum_{\alpha=1}^{s}\bigl{(}w_{l_{1},...,r_{1},...,l_{\alpha}-1,l_{\alpha}+1,...,r_{\alpha},...,l_{s},...,r_{s}}\\\
+w_{l_{1},...,r_{1},...,l_{\alpha},...,r_{\alpha}-1,r_{\alpha}+1,...,l_{s},...,r_{s}}\\\
-[\sum_{a=l_{\alpha}}^{r_{\alpha}-1}(2+2\Delta)+2]w_{l_{1},...,r_{1},...,l_{s},...,r_{s}}\bigr{)},\end{split}$
(12)

Let’s now make few remarks to underline the distinction which exists between
our model and the ferromagnetic spin-1/2 Heisenberg model: in the Bethe Ansatz
solution of this latter model, the coefficient of the third term in eq. (12)
would remain $-2$, regardless of how long the string of consecutive spin downs
would be. This is the main difference between the two models which will
eventually lead to our different solution of the Bethe Ansatz equations, as
explicitly illustrated below. One might think that our derivation of Bethe
Ansatz equations, obtained from subtracting these two equations and requiring
them both to hold simultaneously, would put a restriction to our solution,
namely that the number of components must exceed half of the size of the
system for (11) to be valid. However, a closer examination reveals that the
subtraction is merely a shortcut in deriving the Bethe Ansatz equations, since
alternatively one can compare the coefficients of independent monomials on
both sides of (12) and get the same equations. That is to say, our solution
exists for arbitrary number of components $N$.

###  III.2 Elementary two-body example: the
$\\{c_{1}^{2},c_{2},\cdots,c_{L-2}\\}$ sector

####  III.2.1 The isotropic case

The following two subsections serve as a very elementary introduction to
readers who may be not familiar with the Bethe Ansatz formalism. More
experienced readers should jump directly to our main result in III.3. The case
of $n=2$ as we will see shortly corresponds to a 2-quasi-particle excitation
from the background of momentum-$\pi$ species. Considering this case is
particularly illuminating as it allows an exact diagonalization by solving a
second-order linear homogeneous recurrence relation. To see this, we perform
the transformation

$w_{i_{1},i_{2}}\rightarrow a_{r}^{(s)}=w_{r,r+s},$ (13)

where $r$ labels the site of the first occurrence of the identical species,
while $s$ denotes the distance between the two of them. The eigenvector is
then expressed as

$|v\rangle=\sum_{r=1}^{L}\sum_{s=1}^{L-r}a_{r}^{(s)}|{r,r+s}\rangle.$ (14)

And the anti-periodic boundary condition is given by

$a_{r}^{(L-r+1)}=-a_{1}^{(r-1)}.$ (15)

We first exploit the translational invariance of the Hamiltonian by
diagonalizing in the basis formed by eigenvectors of the translation operator

$T|v\rangle=t^{-1}|v\rangle.$ (16)

According to the anti-periodic boundary condition (7), we have

$T^{L}|1,2\rangle=T^{2}|L-1,L\rangle=-T|1,L\rangle=|1,2\rangle.$ (17)

So $T^{L}=1$ still applies, giving

$t=e^{2i\theta},\quad\theta=\frac{i\pi k}{L},\quad k=1,...,L.$ (18)

Therefore, (16) implies $a_{r+1}^{(s)}=ta_{r}^{(s)}$, and up to a
normalization constant, we can choose $a_{r}^{(1)}=t^{r}$ as the initial
condition for recurrence equations of $a_{r}^{(s)}$ with respect to index $s$.

The recurrence relation (10) now becomes

$\displaystyle a_{r}^{(2)}+a_{r-1}^{(2)}=$ $\displaystyle ea_{r}^{(1)}$ (19)
$\displaystyle a_{r+1}^{(s-1)}+a_{r}^{(s-1)}+a_{r}^{(s+1)}+a_{r-1}^{(s+1)}=$
$\displaystyle ea_{r}^{(s)},\quad s>1,$ (20)

where $e=-(E_{v}+N-4)$. Using translation invariance, the first equation gives
the second initial condition necessary to solve the second-order recurrence
relation $a_{r}^{(2)}=\frac{et^{r+1}}{1+t}$, while the second equation becomes

$a_{r}^{(s+1)}-\frac{et}{1+t}a_{r}^{(s)}+ta_{r}^{(s-1)}=0.$ (21)

The solution of its characteristic equation depends on its discriminant

$D=\frac{e^{2}t^{2}}{(1+t)^{2}}-4t\,\,\,.$ (22)

#####  $D=0$: Identical root:

In this scenario, the general solution is of the form
$a_{r}^{s}=(c_{0}+c_{1}s)t^{s}$, which cannot satisfy the anti-periodic
boundary condition (15) and the initial conditions simultaneously. We conclude
that the discriminate must be nonzero.

#####  $D\neq 0$: Different root:

The characteristic equation has two roots in this case

$\displaystyle\lambda_{\pm}=$ $\displaystyle\left(\frac{e}{4\cos\theta}\pm
i\sqrt{1-\left(\frac{e}{4\cos\theta}\right)^{2}}\right)e^{i\theta}$ (23)
$\displaystyle=$ $\displaystyle e^{i(\theta\pm\alpha)},$ (24)

where $\alpha=\arccos(\frac{e}{4\cos\theta})$. Here, we have assumed $|e|\leq
4|\cos\theta|$, which is necessary for the anti-periodic boundary to hold.
Together with the initial conditions given by $a_{r}^{(1)}$ and $a_{r}^{(2)}$,
(21) gives

$a_{r}^{(s)}=t^{r+(s-1)/2}\frac{\sin s\alpha}{\sin\alpha}.$ (25)

Plugging this into the anti-periodic boundary condition (15), we get

$\alpha=\theta+\theta^{\prime},\quad\theta^{\prime}=\frac{2\pi
k^{\prime}}{L},\quad k^{\prime}=1,2,...,L.$ (26)

We can parameterize the energy eigenvalues with
$\theta_{1}\equiv\theta-\alpha=-\theta^{\prime}$ and
$\theta_{2}\equiv\theta+\alpha=2\theta+\theta^{\prime}$, and the eigenvectors
with $\mu_{1,2}\equiv e^{i\theta_{1,2}}$ as

$\displaystyle E(\theta_{1},\theta_{2})=$
$\displaystyle-L+4-2\cos\theta_{1}-2\cos\theta_{2}$ (27)
$\displaystyle|\theta_{1},\theta_{2}\rangle=$
$\displaystyle\frac{1}{N}\sum_{1\leq i<j\leq
L}(\mu_{1}^{i}\mu_{2}^{j}-\mu_{1}^{j}\mu_{2}^{i})|i,j\rangle,$ (28)

where $\theta_{1,2}=i2\pi k_{1,2}/N$, with $k_{1,2}=1,2,...,N$. Since its
coefficients $a_{i,j}=-a_{j,i}$, we can alternatively express it as

$|\theta_{1},\theta_{2}\rangle=\frac{1}{N}\sum_{1\leq i,j\leq
L}\mu_{1}^{i}\mu_{2}^{j}|i,j\rangle,$ (29)

with $|i,i\rangle\equiv 0$. Notice the $\theta_{1}=\theta_{2}$ solutions give
vanishing eigenvectors, as is expected when one attempts to antisymmetrize a
symmetrized pair of identical species. In other words, when the anisotropy or
density-density interaction is absent, these identical particles of same
species behave as free fermions. The fact that our original bosonic degrees of
freedom has ended up behaving as fermions is exactly due to the unitary
transformation pointed out by Sutherland [9] of multiplying all wave functions
by completely anti-symmetric one in all objects, which flips the sign of the
Hamiltonian and the exchange statistics of the particles at the same time. The
ground state energy in this sector is therefore $-L+2-2\cos\frac{2\pi}{L}$. In
the thermodynamic limit, it gives rise to a vanishing gap that scales with
$1/L^{2}$ counting from the universal ground state of $L$ different species.

The partial spectrum and corresponding eigenstates given above agree with the
results from nested Bethe Ansatz by Sutherland [9], which starts from a
reference state of identical species and consider the different ones as quasi-
particles moving around. To map the eigenstates solved from our approach of
treating the different species as background and identical ones as quasi-
particles, to the framework of Bethe Ansatz, the different species are
understood to have $\pi$ momenta. Our solution corresponds to two particles
having momenta $\theta_{1}$ and $\theta_{2}$, with a scattering phase of
$\pi$, which says the two identical species albeit impossible to be
antisymmetrized in real space, is antisymmetrized in the momentum space.

####  III.2.2 Two-body problem with anisotropy/interaction

Since the interaction term acts on neighboring sites, it would only change the
$e$ in (19) to $e+2\Delta$ in determining the initial condition $a_{r}^{2}$,
leaving the recurrence equation unchanged. Their solution is given by

$\displaystyle a_{r}^{(s)}(\Delta)=$
$\displaystyle\frac{t^{r+(s-1)/2}}{2i\sin\alpha}\left[\left(\frac{e+2\Delta}{2\cos\theta}-e^{-i\alpha}\right)e^{i(s-1)\alpha}+\left(e^{i\alpha}-\frac{e+2\Delta}{2\cos\theta}\right)e^{-i(s-1)\alpha}\right]$
$\displaystyle=$ $\displaystyle
t^{r+(s-1)/2}\left(\frac{\sin(s\alpha)}{\sin\alpha}+\frac{\Delta}{\cos\theta}\frac{\sin(s-1)\alpha}{\sin\alpha}\right),$

which reproduces (25) when $\Delta$ is taken to be $0$. Quantization of the
energy is once again given by the solution of the anti-periodic boundary
condition, now written as

$e^{iL(\alpha-\theta)}=\frac{\cos\theta+\Delta e^{i\alpha}}{\cos\theta+\Delta
e^{-i\alpha}}$ (30)

Taking the logarithm on both sides, one realize that for each value of
$\theta=\frac{k\pi}{N}$, there are $N$ roots $\alpha_{J}$, labeling the $N$
degenerate eigenstates at energy level
$E(k,J)=-N+4-4\cos\theta\cos\alpha_{J}$,

$L(\alpha_{J}-\theta)+2\pi
J=2\arctan\frac{\Delta\sin\alpha_{J}}{\cos\theta+\Delta\cos\alpha_{J}},\quad
J=1,2,...,L.$ (31)

Parametrizing as before, $\theta_{1}\equiv\theta-\alpha$,
$\theta_{2}\equiv\theta+\alpha$, and $\mu_{1,2}\equiv e^{i\theta_{1,2}}$, we
still have

$\displaystyle E(\theta_{1},\theta_{2})=$
$\displaystyle-L+4-2\cos\theta_{1}-2\cos\theta_{2}$ (32)
$\displaystyle|\theta_{1},\theta_{2}\rangle=$
$\displaystyle\frac{1}{\mathcal{N}}\sum_{1\leq i<j\leq
L}(\mu_{1}^{i}\mu_{2}^{j}-S_{12}\mu_{1}^{j}\mu_{2}^{i})|i,j\rangle,$ (33)

where $S_{12}=\mu_{1}^{-L}\equiv\mu_{2}^{L}$, since
$(\mu_{1}\mu_{2})^{L}\equiv 1$, and normalization constant
$\mathcal{N}^{2}=L^{2}-L-Le^{iL\theta}\sin(L-1)\alpha/\sin\alpha$.

###  III.3  The multi-mode scars: $\\{c_{1}^{n},c_{2},\cdots,c_{L-n+1}\\}$
sector

When there are more than two identical particles in one species, the recursive
equations as in the previous section will involve multiple indices. Therefore,
it is necessary to construct the Bethe Ansatz wave function

$w_{i_{1},...,i_{n}}=\sum_{\sigma\in\mathfrak{S}_{n}}A_{\sigma}\prod_{a=1}^{n}\mu_{\sigma
a}^{i_{a}},$ (34)

where $\mathfrak{S}_{n}$ denotes symmetric group of order $n$, and
$\mu_{a}=e^{i\theta_{a}}$. When $i_{a+1}>i_{a}+1$ for all $a$, (10) takes the
form

$\sum_{\sigma\in\mathfrak{S}_{n}}A_{\sigma}\sum_{a=1}^{n}(\mu_{\sigma
a}^{-1}+\mu_{\sigma a}-2)\prod_{a=1}^{n}\mu_{\sigma
a}^{i_{a}}=-(E_{v}+L)\sum_{\rho\in S_{n}}A_{\rho}\prod_{a=1}^{n}\mu_{\sigma
a}^{i_{a}}.$ (35)

The common factor on the l.h.s. can be taken out of the sum over $\sigma$,
which gives $E_{v}=-L+\sum_{a=1}^{n}(2-2\cos\theta_{a})$. When
$i_{a_{0}+1}=i_{a_{0}}+1$ for some $a_{0}$, but $i_{a+1}>i_{a}+1$ for $a\neq
a_{0}$,(10) becomes

$\begin{split}\sum_{\sigma\in\mathfrak{S}_{n}}A_{\sigma}\bigl{(}\mu_{\sigma
a_{0}}^{-1}+\mu_{\sigma(a_{0}+1)}-4-2\Delta+\sum_{a\neq
a_{0}}(2\cos\theta_{\sigma a}-2)\bigr{)}\prod_{a=1}^{n}\mu_{\sigma
a}^{i_{a}}\\\ =-(E_{v}+L)\sum_{\rho\in S_{n}}A_{\rho}\prod_{a=1}^{n}\mu_{\rho
a}^{i_{a}}.\end{split}$ (36)

Noticing that (35) still holds in this case, we can subtract this equation
from it to get

$\sum_{\sigma\in\mathfrak{S}_{n}}A_{\sigma}(\mu_{\sigma
a_{0}}+\mu_{\sigma(a_{0}+1)}^{-1}+2\Delta)\prod_{a=1}^{n}\mu_{\sigma
a}^{i_{a}}=0,$ (37)

or using explicitly $i_{a_{0}+1}=i_{a_{0}}+1$,

$\sum_{\sigma\in\mathfrak{S}_{n}}A_{\sigma}\bigl{(}\mu_{\sigma
a_{0}}\mu_{\sigma(a_{0}+1)}+1+2\Delta\mu_{\sigma(a_{0}+1)}\bigr{)}\mu_{\sigma(a_{0}+1)}^{-1}\prod_{a=1}^{n}\mu_{\sigma
a}^{i_{a}}=0.$ (38)

Combining terms from permutations differing by a transposition $(a_{0}\
a_{0}+1)$, one can see that for this equation to be satisfied for any
coordinates $i_{a}$’s and momenta $k_{a}$’s, it requires

$\frac{A_{\sigma^{\prime}}(\Delta)}{A_{\sigma}(\Delta)}=-\frac{\mu_{\sigma
a_{0}}\mu_{\sigma(a_{0}+1)}+1+2\Delta\mu_{\sigma(a_{0}+1)}}{\mu_{\sigma
a_{0}}\mu_{\sigma(a_{0}+1)}+1+2\Delta\mu_{\sigma a_{0}}},$ (39)

where $\sigma^{\prime}=\tau_{a_{0}}\sigma$, with $\tau_{a_{0}}$ transposing
$a_{0}$ and $a_{0}+1$. In appendix A, we introduce different parameterization
of the momenta for different magnitudes of the anisotropy, in which the
scattering phase depends only on the difference of the rapidities of the two
scattering particles, so as to facilitate a unified treatment of nested and
algebraic Bethe Ansatz in the next section. Quantization condition comes from
the periodic boundary condition (7). Decomposing the cycle-$n$ permutation
$\tau$ into a product of transpositions between neighbors, and using (39), we
have

$\prod_{a=1}^{n}\frac{\Delta\cos\frac{1}{2}(\theta_{\sigma n}-\theta_{\sigma
a})+\cos\frac{1}{2}(\theta_{\sigma n}+\theta_{\sigma
a})+i\Delta\sin\frac{1}{2}(\theta_{\sigma n}-\theta_{\sigma
a})}{\Delta\cos\frac{1}{2}(\theta_{\sigma n}-\theta_{\sigma
a})+\cos\frac{1}{2}(\theta_{\sigma n}+\theta_{\sigma
a})-i\Delta\sin\frac{1}{2}(\theta_{\sigma n}-\theta_{\sigma
a})}=e^{i\theta_{\sigma n}L},$ (40)

where we have absorbed the minus signs into the scattering phases. Since (43)
holds for any $\sigma\in S_{n}$, we have the Bethe Ansatz equations

$2\sum_{b=1}^{n}\arctan\frac{\Delta\sin\frac{1}{2}(\theta_{a}-\theta_{b})}{\Delta\cos\frac{1}{2}(\theta_{a}-\theta_{b})+\cos\frac{1}{2}(\theta_{a}+\theta_{b})}=L\theta_{a}+2\pi
J_{a}$ (41)

where $J_{a}=1,...,L$, for any $a=1,...,n$. Taking the sum on both sides over
$a$, we get $\sum_{a=1}^{n}\theta_{a}=2\pi k/L$, with $k=1,...,L$.

Once the momenta are solved from the Bethe Ansatz equations, (39) can be used
to deduce the relative weights of the permutations among momenta, which gives
the eigenstates in this sector. If the interaction term is turned off, (39)
becomes

$A_{\sigma}=-A_{\sigma^{\prime}}.$ (42)

Therefore $A_{\sigma}=(-1)^{\sigma}$. Boundary condition now becomes

$\sum_{\sigma}A_{\sigma}\mu_{\sigma n}^{L+1}\prod_{a=1}^{n-1}\mu_{\sigma
a}^{i_{a}}=(-1)^{n-1}\sum_{\rho}A_{\rho}\mu_{\rho 1}\prod_{a=2}^{n}\mu_{\rho
a}^{i_{a}},$ (43)

for any $i_{1},...,i_{n}$. Comparing the terms on both sides differing by a
global translation $\tau=(1\ 2\ \cdots\ n)$, we have

$(-1)^{n-1}=\frac{A_{\sigma\tau^{-1}}}{A_{\sigma}}=(-1)^{n-1}\mu_{\sigma
n}^{L},\quad\forall\sigma\in S_{n}.$ (44)

Hence, $\theta_{a}=2\pi k_{a}/L$, and the ground state energy in the sector
with $n$ particles belonging to the identical species is

$E\,=\,-L+2\sum_{k=-m}^{m}\left(1-\cos 2\pi k/L\right)$ (45)

for $n=2m+1$, or

$E=-L+2\sum_{k=-m+1}^{m}\left(1-\cos 2\pi k/L\right)$ (46)

for $n=2m$, since if any two of the momenta coincide, the anti-symmetrized
coefficient (34) would vanish. We can use the Lagrange’s trigonometric
identity to estimate their gap from the ground state in the thermodynamic
limit $n\to\infty$,

$E-E_{0}\to 2L\left(\frac{n}{L}-\frac{2}{\pi}\sin\frac{\pi n}{L}\right)$ (47)

With some extra steps explained in appendix B, one can even show that the
normalization constant is fixed to be $\mathcal{N}=L^{n/2}$. The eigenstates
are then expressed as

$|\theta_{1},...,\theta_{n}\rangle=L^{-n/2}\sum_{1\leq i_{1}<\cdots<i_{n}\leq
L}\sum_{\sigma\in\mathfrak{S}_{n}}(-1)^{\sigma}e^{i\sum_{a=1}^{n}\theta_{\sigma
a}i_{a}}|i_{1},...,i_{n}\rangle.$ (48)

The careful reader might now recognize that, choosing $n$ momenta among their
$L$ possible values allowed by quantization condition, the dimension of this
eigenvector subspace $\binom{L}{n}$, is precisely the dimension of the
irreducible representation corresponding to the one-column Young tableaux with
$n$ rows, whereas the absolute ground state of the full Hilbert space
expressed by Slater determinant correspond to the $1$ dimensional irreducible
representation of the $L$ row, $1$ column tableaux.

We emphasize that the integrable eigenstates found above have finite energy
density (evaluated wrt the ground state energy) in the thermodynamic limit. To
see that, notice that the energy of eigenstates in this sector consists of two
parts, a contribution from the antisymmetrized basis, at the lowest energy
$\epsilon_{\text{GS}}=-L$, and a contribution from the solution of
spin-$\frac{1}{2}$ XXZ problem, of the form
$L\epsilon_{\text{XXZ}}(S^{z}_{\text{tot}}=\frac{L}{2}-n)$, where lower and
upper bound on $\epsilon_{\text{XXZ}}(S^{z}_{\text{tot}}=\frac{L}{2}-n)$ is
established to be finite [67]. So the total energy density in the
thermodynamic limit is

$\epsilon_{v}-\epsilon_{\text{GS}}=\epsilon_{\text{XXZ}}.$ (49)

##  IV  Violation of YBE: the
$\\{c_{1}^{n_{1}},\cdots,c_{m}^{n_{m}},c_{m+1},\cdots,c_{L-n+m}\\}$ sector

The spectrum of a generic sector, where multiple colors appear multiple times
in the particle content, is solved by recourse to the Bethe-Yang hypothesis
[7], used in Yang’s solution to the one-dimensional Fermi problem with
repulsive $\delta$ interaction. As well known, the same model was originally
and independently solved by Gaudin, using a less physical and more
algebraically involved approach [64],) based on the pioneering work of McGuire
[68, 69], Lieb and Flicker [70], and subsequently further generalized by
Sutherland [8, 9].

Yang’s approach differs from Bethe Ansatz used in the previous section since
it does not assume any constraint on the permutation symmetry of the wave
function, which is indeed the case when multiple species are involved.
Therefore, if we still label the sites of all identical species with $i$’s,
eq. (6) no longer holds. Consequently, the basis $|i_{1},\cdots,i_{n}\rangle$
no longer spans the whole Hilbert space when restricted to
$i_{1}<\cdots<i_{n}$. Instead, we should write the eigenvector of the
Hamiltonian as

$|v\rangle=\sum_{Q\in\mathfrak{S}_{n}}\sum_{1\leq x_{Q1}<\cdots,x_{Qn}\leq
L}\psi_{Q}(\mathbf{x})|\mathbf{x}\rangle,$ (50)

where $\mathbf{x}$ is an $n$-component array of the combinations of coordinate
and spin/color $\\{q_{i};c_{i}\\}$, and $x_{i}<x_{j}$ if $q_{i}<q_{j}$, or
$c_{i}<c_{j}$ if $q_{i}=q_{j}$. The simultaneous swapping of a pair both of
them gives a minus sign due to the fermionic nature of the basis,

$|\mathbf{x}\tau\rangle=\text{sgn}(\tau)|\mathbf{x}\rangle,$ (51)

where
$\mathbf{x}\tau=(x_{\tau(1)},...,x_{\tau(n)})\equiv(\\{q_{\tau(1)};c_{\tau(1)}\\},...,\\{q_{\tau(n)};c_{\tau(n)}\\})$.
Therefore, only the anti-symmetric part of the wave function multiplying an
anti-symmetric basis is meaningful

$\psi_{\tau^{-1}Q}(\mathbf{x}\tau)=\text{sgn}(\tau)\psi_{Q}(\mathbf{x}).$ (52)

While we are discussing permutation symmetry, it’s worth mentioning theorems
on the ordering of energy levels of the isotropic Hamiltonian. According to
Lieb and Mattis’ theorem on the ordering of energy levels of antiferromagnets,
the lowest energy eigenstate is determined by the so-called “pouring
principle” [71, 72]. More recently, this result has been generalized to a
lattice model with higher spin with both open and periodic boundary conditions
[73]. In short, it states that the lowest energy eigenstate corresponds to the
Young Tableaux that gives the highest weight state of the $sl(n)$ algebra. The
proof is given much in the same spirit of the one due to Lieb and Mattis,
namely, showing that the representation corresponding to the Young Tableaux of
highest weight state follows from the Perron-Frobenius theorem and the
requirement to have non-negative components: hence, it overlaps with the
ground state, but, given at the same time that the Hamiltonian does not mix
representations, one can conclude that it has to be the only representation
corresponding to the ground state. Irreducible representations of outer
product can be constructed according to a simple rule from the irreducible
representations of each group in the outer product, which in our case, has to
be the single row Young tableaux $[n_{1}],\cdots,[n_{m}]$, as one cannot
antisymmetrize among an identical species. Then the Young Tableaux in the
outcome of the outer product that every other one can “pour into”, alias
corresponding to the highest weight state, is $[n_{1},\cdots,n_{m}]$, where
each row is filled with the same species, assuming $n_{1}>n_{2}>\cdots>n_{m}$.
The above analysis relies purely on the permutation symmetry of the
Hamiltonian. As we will see shortly, Yang’s approach utilizes both permutation
symmetry and the translational invariance, and is more convenient in practice.
It starts with a first layer of Bethe Ansatz treating all particles as if they
have the same color, while keeping in mind that their colors could be
different later by keeping track of the swapping of color indices

$\psi_{Q}(\mathbf{x})=\sum_{P\in\mathfrak{S}_{n}}A_{Q,P}\prod_{i=1}^{n}\mu_{Pi}^{q_{Qi}},$
(53)

for each sector $1\leq x_{Q1}<x_{Q2}<\cdots<x_{Qn}\leq L$ in the coordinate
space labeled by a permutation $Q\in\mathfrak{S}_{n}$, which are a priori
independent. Each $Q$ corresponds to a particular sector of the whole
coordinate space, with the $Q=id$ one called the fundamental sector. $A_{Q,P}$
is an $n!\times n!$ matrix, whose columns are denoted by $\xi_{P}$. To relate
the wave functions defined in these separate sectors, we use the boundary
conditions between two sectors $Q^{\prime}=Q\tau_{a}$ that differ by a
transposition between indices $a$ and $a+1$. In appendix C, we give the
detailed derivation of Yang’s scattering matrix $Y$ from comparing the
eigenvalue equations of nonadjacent and adjacent cases, which requires

$\xi_{P\sigma_{ij}}=Y_{ij}^{a}\xi_{P},$ (54)

for $P(a)=i,P(a+1)=j$, where the $Y$ matrix as defined in Appendix C. These
$n!(n-1)$ equations are consistent for the isotropic Hamiltonian only. For
completeness, we carry out its diagonalization following Sutherland [9].

The above Yang’s $Y$ matrices are given in the reflection(-diagonal)
representation of the asymptotic wave function, according to Sutherland [74],
which is convenient for relating the a priori independent wave functions
defined in separate sectors, as the coordinate of each particle remains
invariant after scattering in this representation. However, in order to
utilize the periodic boundary condition as a quantization condition for the
momenta, we switch to the transmission(-diagonal) representation, where the
diagonal terms of scattering $S$ matrices are the transmission amplitudes. In
the transmission representation, momentum is fixed to a particle after
scattering. Picking a particular permutation of momenta $P_{0}=Q$, such that
$i=a,j=a+1$, we have

$S_{ij}=\pi(\tau_{ij})Y_{ij}^{ij}=T_{ij}\mathbb{1}+R_{ij}\pi(\tau_{ij}),$ (55)

where reflection $R_{jj}$ and the transmission $T_{jj}$ coefficients are
defined as in eqs. (57),(56) herefater. In Appendix C, we show that in the
presence of anisotropy in the Hamiltonian, the Yang-Baxter’s equation, which
is necessary for the following procedures to work, is only satisfied for sub-
sectors of the Hilbert space corresponding to two types of Young tableaux,
namely those of the shape illustrated in Fig. 2 (c) and (d). The former
corresponds to the solution of spin-$\frac{1}{2}$ XXZ model, as well as the
partially integrable subspaces studied in [12, 13, 14]. The latter corresponds
to our scarred sectors. In fact, the Young tabeaux in Fig. 2 (c) correspond to
an even broader class of integrable eigenstate: while the antisymmetric bases
requires the constituents to be all distinguishable, a symmetric basis does
not require its constituents to be identical. So we can also form symmetric
bases from those different components and obtain the upper spectrum
counterparts of the eigenstates solved in the previous section. The important
point is that we can either form totally symmetric or antisymmetric irrep’s,
but not a mixture of them in order for the YBE to be satisfied.

Figure 2: (a) Young diagram corresponding to a generic sector in the Hilbert
space of our model, integrable only in the isotropic case. (b) The Young
diagram corresponding to the sector where our ground state of a model with
only a finite number $s$ of colors divisible by the size of the chain, not
integrable as anisotropy mixes irreducible representations corresponding to
different Young tableaux. (c) Young diagram corresponding to
spin-$\frac{1}{2}$ XXZ model as well as the partially integrable subspaces
studied in [12, 13, 14]. (d) Young diagram corresponding to the scarred
sectors of our model.

Hence, for scattering between different colors, $c_{a}\neq c_{a+1}$, we have

$\displaystyle R_{ij}^{ab}(\Delta)$
$\displaystyle=\frac{-i}{\lambda_{i}-\lambda_{j}+i}\equiv R_{ij}$ (56)
$\displaystyle T_{ij}^{ab}(\Delta)$
$\displaystyle=\frac{\lambda_{i}-\lambda_{j}}{\lambda_{i}-\lambda_{j}+i}\equiv
T_{ij},$ (57)

where

$\lambda_{i}=\frac{i}{2}\frac{\mu_{i}+1}{\mu_{i}-1}.$ (58)

Yet for scattering between the same colors, $c_{a}=c_{a+1}$, we have
$\pi(\tau_{a(a+1)})=-\mathbb{1}$ (due to (52)), and

$Y_{ij}^{a(a+1)}(\Delta)=R_{ij}^{a(a+1)}(\Delta)-T_{ij}^{a(a+1)}(\Delta)=-\frac{\mu_{i}\mu_{j}+1+2\Delta\mu_{j}}{\mu_{i}\mu_{j}+1+2\Delta\mu_{i}}=\frac{\varphi^{\Delta}(\lambda_{j}^{\Delta}-\lambda_{i}^{\Delta}+i\eta)}{\varphi^{\Delta}(\lambda_{j}^{\Delta}-\lambda_{i}^{\Delta}-i\eta)}\equiv\Theta_{ij},$
(59)

with the parameterization (81). It’s easily verified that while the unitarity
relation

$Y_{ij}^{ab}(\Delta)Y_{ji}^{ab}(\Delta)=1,$ (60)

is satisfied for both cases, the Yang-Baxter’s equations

$Y_{jk}^{ab}(\Delta)Y_{ik}^{bc}(\Delta)Y_{ij}^{ab}(\Delta)=Y_{ij}^{bc}(\Delta)Y_{ik}^{ab}(\Delta)Y_{jk}^{bc}(\Delta),$
(61)

requires

$\displaystyle\Theta_{jk}R_{ik}\Theta_{ij}$
$\displaystyle=R_{ij}\Theta_{ik}R_{jk}+T_{ij}R_{ik}T_{jk},$ (62)
$\displaystyle R_{jk}T_{ik}\Theta_{ij}$
$\displaystyle=T_{ij}\Theta_{ik}R_{jk}+R_{ij}R_{ik}T_{jk},$ (63)

which only hold when $\Delta=0$ or $\Delta\to\infty$, see Fig.3.

Figure 3: Yang-Baxter Equation in terms of reflection and transmission
coefficient.

If $\Delta\neq 0$, YBE hold respectively for scattering between identical
colors, and different colors, but not when they are mixed. This means the
scattering among three particles and more, is only factorizable into
consecutive scatterings between pairs of two particles, irrespective of their
ordering, if their colors are either all the same, which leads to the solution
in III.3, or all different, corresponding to the absolute ground state. Hence,
for the generic case, where scattering between particles of both same and
different colors are present, the scattering is not non-diffractive, and
Yang’s nested Bethe Ansatz no longer applies. In fact, Sutherland has derived
a necessary consistent condition for the scattering between different species
to be non-diffractive [74], an example of such an exactly solvable anisotropic
generalization of the multi-component Heisenberg model was solved by Babelon,
et. al. [10], albeit being less physically natural than our model.

Another way to look at the YBE is to find the momenta that solves Eq. (63) for
generic $\Delta$. Multiplying the denominators on both sides and comparing the
coefficients of different orders of $\Delta$ give the following solutions

$\begin{cases}\mu_{k}=1\\\ \mu_{j}=\mu_{i};\\\ \mu_{k}=\mu_{j};\\\
\mu_{k}=\mu_{i};\\\ \mu_{j}=2-\frac{1}{\mu_{i}};\\\
\mu_{k}=2-\frac{1}{\mu_{j}};\\\ \mu_{k}=2-\frac{1}{\mu_{i}}.\end{cases}$ (64)

The first solution is trivial, meaning the particle with a different color has
momentum $\pi$, or is antisymmetrized. This can only happen is appears only
once along the chain for the wave-function not to vanish. So this is included
in the case discussed before. The next three solutions corresponds to two of
the three particles having the same momentum, in which case they move
collectively keeping their distance, so a three-body scattering never happens.
The last three solutions are non-trivial. They require two of the momenta to
be complex-valued, with the corresponding particles forming bond states.

##  V Entanglement entropy of integrable eigenstates

By now it should be clear that our integrable eigenstates are in 1-to-1
mapping with the Yang’s solution to the $\delta$-interaction problem [7] of
$n$ particles, with the antisymmetrized bases of distinguishable particles
playing the role of the vacuum there, or the pseudovaccum of all spin up state
in the XXZ chain. In this section, we try to get a sense of the additional
entanglement that these antisymmetrization contribute, by first show that in
the extreme case of $N=L$, where the antisymmetrized basis occupy the entire
system, the ground state has exactly linear scaling of entanglement entropy in
the thermodynamic limit. Then we argue that the entanglement entropy of
integrable excited states are lower-bounded by the corresponding Bethe Ansatz
states in spin-$\frac{1}{2}$ XXZ chain [58, 59].

###  V.1 Ground state entanglement entropy for $N=L$

In the case of $N=L$, the ground state can be expressed as

$|GS\rangle=\frac{1}{\sqrt{L!}}\sum_{\mathcal{P}\in\mathfrak{S}_{L}}\text{sgn}(\mathcal{P})|c_{\mathcal{P}1},c_{\mathcal{P}2},...,c_{\mathcal{P}L}\rangle,$
(65)

where $|...,c_{i},...\rangle$ denotes the configuration with the $i$-th site
having particle of color $c_{i}$.

This ground state is highly entangled as the color at each site depends on
those at all the others. This can be seen by calculation the entanglement
entropy between two halves of the system. Let $L=2l$, then the Schmidt
decomposition of the ground state is

$|GS\rangle=\sum_{\mathcal{P}_{0}\in\mathfrak{S}_{L}/(\mathfrak{S}_{l})^{2}}\text{sgn}(\mathcal{P}_{0})\frac{l!}{\sqrt{(2l)!}}|\mathcal{P}_{0}1,\mathcal{P}_{0}2,...,\mathcal{P}_{0}l\rangle\otimes|\mathcal{P}_{0}(l+1),\mathcal{P}_{0}(l+2),...,\mathcal{P}_{0}L\rangle,$
(66)

where
$|i_{1},i_{2},...,i_{l}\rangle=\frac{1}{\sqrt{l!}}\sum_{\mathcal{P}\in\mathfrak{S}_{l}}\text{sgn}(\mathcal{P})|c_{\mathcal{P}i_{1}},s_{\mathcal{P}i_{2}},...,s_{\mathcal{P}i_{l}}\rangle$.The
entanglement entropy between two subsystems is then

$S=-{2l\choose l}\frac{(l!)^{2}}{(2l)!}\log\frac{l!^{2}}{(2l)!}.$ (67)

In the thermodynamic limit $l\to\infty$, this is approximated by the sterling
formula as

$\displaystyle S$ $\displaystyle\simeq 2l(\log(2l)-1)-2l(\log l-1)$
$\displaystyle\simeq L\log 2.$

###  V.2 Entanglement entropy of integrable excited states

An integrable eigenstate labeled by the $n$ momenta
$\\{\theta_{1},\cdots,\theta_{n}\\}$ of the identical color $c_{1}$, and the
colors $\\{c_{2},\cdots,c_{L-n+1}\\}$ of the rest $L-n$ sites can be expressed
as the superposition of bipartition across the middle over different
possibility of coloring in the left half system specified by the number of
site in the identical color $n_{A}$, and the combination of the rest of the
colors $C$:

$\begin{split}|\\{\theta_{1},\cdots,\theta_{n}\\},\\{c_{1};c_{2},\cdots,c_{L-n+1}\\}\rangle=\sum_{n_{A}=n-l}^{l}\sum_{C=1}^{\binom{L-n}{l-n_{A}}}&\text{sgn}(C)\sqrt{\frac{(l-n_{A})!(l-n+n_{A})!}{(L-n)!}}\sum_{\\{i_{1},\cdots,i_{n_{A}}\\}}\sum_{\\{i_{n_{A}+1},\cdots,i_{n}\\}}\frac{w_{i_{1},\cdots,i_{n}}}{\sqrt{\mathcal{N}}}\\\
&\frac{1}{\sqrt{(l-n_{A})!}}|i_{1},\cdots,i_{n_{A}}\rangle\otimes\frac{1}{\sqrt{(l-n+n_{A})!}}|i_{n_{A}+1},\cdots,i_{n}\rangle,\end{split}$
(68)

where $\mathrm{sgn}(C)$ is the signature of the permutation of the coloring of
the whole system divided by the those of the two subsystems, and $\mathcal{N}$
takes care of the normalization for the usual Bethe wave function without
antisymmetric basis. Plugging in the definition of the bases vectors (5), and
taking partial trace of the density matrix over the right subsystem, we get
the reduced density matrix of the left subsystem

$\begin{split}\rho_{A}=&\mathrm{tr}_{B}|\\{\theta_{1},\cdots,\theta_{n}\\},\\{c_{1};c_{2},\cdots,c_{L-n+1}\\}\rangle\langle\\{\theta_{1},\cdots,\theta_{n}\\},\\{c_{1};c_{2},\cdots,c_{L-n+1}\\}|\\\
=&\sum_{n_{A}=n-l}^{l}\sum_{C=1}^{\binom{L-n}{l-n_{A}}}\frac{(l-n_{A})!(l-n+n_{A})!}{(L-n)!}\sum_{\\{i_{1},\cdots,i_{n_{A}}\\}}\frac{\sum_{\\{i_{n_{A}+1},\cdots,i_{n}\\}}|w_{i_{1},\cdots,i_{n}}|^{2}}{\mathcal{N}}|i_{1},\cdots,i_{n_{A}}\rangle\langle
i_{1},\cdots,i_{n_{A}}|\\\
=&\bigoplus_{n_{A}=n-l}^{l}\bigoplus_{C=1}^{\binom{L-n}{l-n_{A}}}\frac{(l-n_{A})!(l-n+n_{A})!}{(L-n)!}\rho_{A}^{\mathrm{XXZ}}(n_{A}),\end{split}$
(69)

where $\rho_{A}^{\mathrm{XXZ}}(n_{A})$ denotes the block of the reduced
density matrix of the XXZ chain with $n_{A}$ spin down in the left subsystem.
The entanglement entropy of the integrable excited states are computed from
the Schmidt coefficients

$p_{A}(n_{A};C;\\{i_{n_{A}+1},\cdots,i_{n}\\})=\frac{(l-n_{A})!(l-n+n_{A})!}{(L-n)!}p_{A}^{\mathrm{XXZ}}(n_{A};\\{i_{n_{A}+1},\cdots,i_{n}\\}),$
(70)

which are given in terms of the Schmidt coefficients of the corresponding XXZ
excited states $p_{A}^{\mathrm{XXZ}}(n_{A};\\{i_{n_{A}+1},\cdots,i_{n}\\})$.

$\displaystyle S_{A}=$
$\displaystyle-\sum_{n_{A}=n-l}^{l}\sum_{\\{i_{1},\cdots,i_{n_{A}}\\}}p_{A}^{\mathrm{XXZ}}(n_{A};\\{i_{n_{A}+1},\cdots,i_{n}\\})\log\big{(}\frac{(l-n_{A})!(l-n+n_{A})!}{(L-n)!}p_{A}^{\mathrm{XXZ}}(n_{A};\\{i_{n_{A}+1},\cdots,i_{n}\\})\big{)}$
(71) $\displaystyle=$ $\displaystyle
S_{A}^{\mathrm{XXZ}}+\log(L-n)!-\sum_{n_{A}=n-l}^{l}p_{A}^{\mathrm{XXZ}}(n_{A})\big{(}\log(l-n_{A})!+\log(l-n+n_{A})!\big{)},$
(72)

where
$p_{A}^{\mathrm{XXZ}}(n_{A})=\mathrm{tr}_{A}\rho_{A}^{\mathrm{XXZ}}(n_{A})$ is
the sum of the Schmidt coefficient of XXZ eigenstates with a fixed number of
down spins $n_{A}$ in subsystem A. So the entanglement entropy of our
integrable excited states factorizes nicely into the entanglement contribution
from the Bethe Ansatz state and the antisymmetrized bases. In the
thermodynamic limit of $L=2l\to\infty$,

$S_{A}-S_{A}^{\mathrm{XXZ}}\sim(2l-n)\log(2l-n)-\sum_{n_{A}=n-l}^{l}p_{A}^{\mathrm{XXZ}}(n_{A})\big{(}(l-n_{A})\log(l-n_{A})+(l-n+n_{A})\log(l-n+n_{A})\big{)}.$
(73)

The additional entanglement entropy from the bases varies with the number of
colors involved, and the Schmidt decomposition of the particular XXZ
eigenstate in question. It ranges between $0$ and $(L-n)\log 2$. So the
entanglement entropy of the integrable excited states are lower-bounded by
those of the corresponding Bethe Ansatz states in spin-$\frac{1}{2}$ model,
which are area law breaking [58, 59], and upper-bounded by volume law, as a
function of the number of components available in the Hilbert space.

##  VI  Conclusions and Outlooks

In this paper we have addressed the integrable excited states in a partially
integrable $\mathfrak{S}_{N}$-invariant antiferromagnetic multi-component
Heisenberg chain. The model is not integrable in symmetry sectors involving
more than two components, because three-body scatterings between quasi-
particles depend generally on the ordering of the factorized two-body
scatterings, as interaction between identical and different components are
different. However, in addition to the more apparently integrable Krylov
subspaces studied before [12, 13, 14], we found a large number of Bethe Ansatz
integrable excited states in the symmetry sectors involving more than two
components, by constructing an antisymmetized basis among components that
appear once in the chain. The excited states we found differ from the
previously studied set in the integrable symmetry sectors in that an initial
state in the non-integrable sector overlaps both with integrable eigenstates
and ETH-violating ones. So there will be consequences for what concerns the
slow dynamics of local observables toward thermalization as a weak form of ETH
violation, i.e. a process that interpolates between integrable and chaotic
systems.

While our discussion has referred to a local Hilbert space of arbitrary $N$
dimension (we have even used thermodynamically large $N$ as an extreme case to
illustrate the scaling of the contribution to entanglement entropy from the
antisymmetrized bases), it is worth stressing that these integrable excited
states exist even for as small a number of components as $N=3$: this implies
that, employing a mapping between $SU(3)$ and $SU(2)$ spin-1, these states can
be realized in cold atom experiment. Although as $N$ gets smaller, integrable
excited states are more outnumbered by the non-integrable ones, we speculate
that the slowdown in thermalization can still be observed with a smart choice
of initial state in an array of smaller number of cold atoms. However a
proposal for the realization of the anisotropy interaction goes beyond the
scope of the current manuscript and is deferred to future work.

Our approach to arrive at weak ergodicity breaking should be compared with the
paradigm of QMBS. In those models, one usually starts with a non-integrable
model of which little is known with the possible exception of a solvable
ground state. So the non-integrability is established by numerical evidence of
level-spacing statistics. Then by applying SMA to the solvable ground state,
one can construct a tower of exact excited states, usually with $\pi$ momenta.
In our case instead, we start with a model that is Bethe Ansatz integrable,
and adds deformation to break integrability, so we can directly check the YBE
to see the generic breakdown of integrability, as well as identify the
surviving integrable eigenstates. Although the integrable excited states have
finite energy density, they do not fall into the category of exact quantum
many-body scar states, in the following sense. First, given these states are
solved with Bethe Ansatz, their corresponding energy are not integer or
rational valued, nor are they equal distant in the spectrum. So the periodic
revival on top of a slow dynamics characteristic to QMBS will not be
observable here. Second, they exhibit severe violation of area law and are not
expressible by matrix product states with finite bond dimension [36, 75].

We emphasize that violation of YBE only implies the absence of non-diffractive
scattering. When Bethe Ansatz fails, Sommerfeld diffraction ansatz can still
apply if internal consistency depending on solution of certain Riemann-Hilbert
problem is met. Indeed, McGuire and Hurst has developed algebraic formulation
for solving the three-body problem with different scattering between particles
of different species [76, 77]. Yet so far, this approach has only been applied
to three- or four-body problems [78], except for scattering with bound states
[79]. At the end of Sec. IV, we have identified exceptional momenta that allow
such exact analytical solutions by solving the YBE as equations of momenta for
arbitrary anisotropy. We believe that this method can be more widely
applicable than model specific result. For instance, one can apply Bethe
Ansatz instead of SMA to the frustration-free ground state of AKLT chain, and
solve the YBE of scattering matrices to either find momenta different from
$\pi$ that allows multi-mode exact excited state, or have a definitive
explanation of the relation between QMBS and $\pi$-momenta. Another direction
to pursue in terms of diagonalization without going into the most general
Sommerfeld diffraction ansatz is to calculate the diffractive scattering
amplitude for the $\Delta\to 0$ case, when the Yang-Baxter equation is weakly
violated [80].

We have been able to express the entanglement entropy of our integrable
excited state as a sum of the entanglement entropy of the corresponding Bethe
Ansatz eigenstates in the XXZ chain, and an additional entanglement
contribution from the antisymmetrized bases. This gives another way to see how
the integrable excited states we found are more non-trivial than those in the
integrable symmetry sector previously studies. We believe such a sum of
entanglement entropy contributions may exist elsewhere, such as the in the
isotropic multi-component model integrable by nested Bethe Ansatz. We also
showed that when the number of components are large enough, the ground state
have volume law scaling of entanglement entropy, and the spectrum is gapless.
It would be interesting to see how these features change as the number of
local degrees of freedom becomes finite, either from numerics or a combination
of analytical and numerical methods.

###### Acknowledgements.

We thank Hosho Katsura for discussions, and especially for bringing to our
attention related results in the literature, and thank Lorenzo Piroli for
discussions. ZZ thanks Israel Klich, Pasquale Calabrese, Feng He, and Haiyuan
Zou for fruitful discussions. ZZ acknowledges the kind hospitality of the
workshop “Randomness, Integrability and Universality” at the Galileo Galilei
Institute for Theoretical Physics during the final stage of the work. The work
of ZZ was partly supported by the National Natural Science Foundation of China
Grant No. 12005129, and GM’s work was supported in part by the grant Prin
2017-FISI.

##  Appendix A Parameterization of Scattering matrices

The reparameterization of momenta into rapidities for different anisotropy
parameters is due to Orbach [81].

#####  $\Delta>1$:

Using parameterization $\Delta\equiv\cosh\eta$, and

$\mu\equiv-\frac{\sin(\lambda+i\eta/2)}{\sin(\lambda-i\eta/2)},$ (74)

(39) becomes

$\frac{A_{\sigma^{\prime}}(\eta)}{A_{\sigma}(\eta)}=\frac{\sin(\lambda_{\sigma(a_{0}+1)}-\lambda_{\sigma
a_{0}}+i\eta)}{\sin(\lambda_{\sigma(a_{0}+1)}-\lambda_{\sigma a_{0}}-i\eta)},$
(75)

#####  $\Delta=1$:

In this case, we parameterize

$\mu\equiv-\frac{\lambda+i/2}{\lambda-i/2},$ (76)

and (39) takes the form

$\frac{A_{\sigma^{\prime}}(\eta)}{A_{\sigma}(\eta)}=\frac{\lambda_{\sigma(a_{0}+1)}-\lambda_{\sigma_{a_{0}}}+i}{\lambda_{\sigma(a_{0}+1)}-\lambda_{\sigma_{a_{0}}}-i}$
(77)

#####  $0\leq\Delta<1$:

Using parameterization $\Delta\equiv\cos\eta$, and

$\mu\equiv-\frac{\sinh(\lambda+i\eta/2)}{\sinh(\lambda-i\eta/2)},$ (78)

(39) then becomes

$\frac{A_{\sigma^{\prime}}(\eta)}{A_{\sigma}(\eta)}=\frac{\sinh(\lambda_{\sigma(a_{0}+1)}-\lambda_{\sigma_{a_{0}}}+i\eta)}{\sinh(\lambda_{\sigma(a_{0}+1)}-\lambda_{\sigma_{a_{0}}}-i\eta)},$
(79)

Summarizing, we have

$\varphi^{\Delta}(\lambda^{\Delta}\pm
i\eta)=\begin{cases}\sin(\lambda^{\Delta}\pm
i\operatorname{arccosh}\Delta),&\Delta>1,\\\ \lambda^{\Delta}\pm
i,&\Delta=1\\\ \sinh(\lambda^{\Delta}\pm
i\arccos\Delta),\quad&0\leq\Delta<1,\end{cases}$ (80)

and

$\mu=\begin{cases}-\frac{\sin(\lambda^{\Delta}+i\eta/2)}{\sin(\lambda^{\Delta}-i\eta/2)},&\Delta>1,\\\
-\frac{\lambda^{\Delta}+i/2}{\lambda^{\Delta}-i/2},&\Delta=1\\\
-\frac{\sinh(\lambda^{\Delta}+i\eta/2)}{\sinh(\lambda^{\Delta}-i\eta/2)},\quad&0\leq\Delta<1.\end{cases}$
(81)

##  Appendix B Normalization of eigenstates

In the absence of the interaction term, the eigenvalues of translation
operator are roots of unity. Therefore it admits a cute calculation of the
normalization constant, as we demonstrate below.

$\displaystyle\mathcal{N}^{2}=$ $\displaystyle\sum_{1\leq i_{1}<...<i_{n}\leq
N}|w_{i_{1},...,i_{n}}|^{2}$ $\displaystyle=$ $\displaystyle\sum_{1\leq
i_{1}<...<i_{n}\leq L}\sum_{\sigma,\rho\in
S_{n}}(-1)^{\sigma+\rho}\prod_{a=1}^{n}(\mu_{\sigma a}\mu_{\rho
a}^{*})^{i_{a}}$ $\displaystyle=$ $\displaystyle\sum_{1\leq
i_{1}<...<i_{n}\leq L}\sum_{\sigma,\tau\in
S_{n}}(-1)^{\tau}\prod_{a=1}^{n}(\mu_{\sigma a}\mu_{\tau\sigma
a}^{*})^{i_{a}}\qquad$ $\displaystyle(\tau=\rho\sigma^{-1})$ $\displaystyle=$
$\displaystyle\sum_{1\leq i_{1}<...<i_{n}\leq L}\sum_{\sigma,\tau\in
S_{n}}(-1)^{\tau}\prod_{a=1}^{n}(\mu_{a}\mu_{\tau
a}^{*})^{i_{\sigma^{-1}a}}\qquad$ $\displaystyle(\text{dummy index}\ a)$
$\displaystyle=$ $\displaystyle\sum_{1\leq i_{1}<...<i_{n}\leq
L}\sum_{\sigma,\tau\in S_{n}}(-1)^{\tau}\prod_{a=1}^{n}(\mu_{a}\mu_{\tau
a}^{*})^{i_{\sigma a}}\qquad$ $\displaystyle(\text{dummy index }\sigma)$
$\displaystyle=$ $\displaystyle\sum_{1\leq i_{1}\neq...\neq i_{n}\leq
L}\sum_{\tau\in S_{n}}(-1)^{\tau}\prod_{a=1}^{n}(\mu_{a}\mu_{\tau
a}^{*})^{i_{a}}$ $\displaystyle=$ $\displaystyle\sum_{\tau\in
S_{n}}(-1)^{\tau}\prod_{a=1}^{n}\sum_{i_{a}=1}^{L}(\mu_{a}\mu_{\tau
a}^{*})^{i_{a}}-\sum_{\begin{subarray}{c}\exists b,c,\\\
i_{b}=i_{c}\end{subarray}}\sum_{\tau\in S_{n}}(-1)^{\tau}(\mu_{b}\mu_{\tau
b}^{*}\mu_{c}\mu_{\tau c}^{*})^{i_{b}}\prod_{a\neq c,d}(\mu_{a}\mu_{\tau
a}^{*})^{i_{a}}$ $\displaystyle=$ $\displaystyle\sum_{\tau\in
S_{n}}(-1)^{\tau}\prod_{a=1}^{n}(\delta_{a,\tau a}L)\qquad$
$\displaystyle(\text{as }\mu_{a}\text{'s are distinct})$
$\displaystyle-\sum_{\begin{subarray}{c}\exists b,c,\\\
i_{b}=i_{c}\end{subarray}}\sum_{\tau\in S_{n}/S_{2}}(-1)^{\tau}\prod_{a\neq
c,d}(\mu_{a}\mu_{\tau a}^{*})^{i_{a}}\sum_{\pi\in\\{id,(b\
c)\\}}(-1)^{\pi}(\mu_{b}\mu_{\tau b}^{*}\mu_{c}\mu_{\tau c}^{*})^{i_{b}}$
$\displaystyle=$ $\displaystyle\sum_{\tau\in
S_{n}}(-1)^{\tau}\delta_{\tau,id}L^{n}$ $\displaystyle=$ $\displaystyle
L^{n}.$

As a special case, when $n=L$, this norm can be easily calculated by
representing the wave function as Slater determinant, and evaluating the
product of determinants with the determinant of the product of the two
matrices.

##  Appendix C Derivation of scattering matrix

In this appendix, we give the detailed derivation of Yang’s scattering $Y$
matrix in the reflection(-diagonal) representation. Using (50), (51) and (52),
the eigenvalue equation $H|v\rangle=E|v\rangle$, when considering adjacent
identical particles case $Qa=Q(a+1)-1$ separately, gives

$\begin{split}&\sum_{Q}\biggl{(}\sum_{x_{Qi}+1<x_{Q(i+1)}}\psi_{Q}\bigl{[}\sum_{i=1}^{n}-\bigl{(}|...,x_{i}-1,...\rangle+|...,x_{i}+1,...\rangle\bigr{)}-(L-2n)|...\rangle\bigr{]}\\\
&+\sum_{x_{Qa}+1=x_{Q(a+1)}}\psi_{Q}\bigl{[}\sum_{i\neq
k,k+1}-\bigl{(}|...,x_{i}-1,...\rangle+|...,x_{i}+1,...\rangle\bigr{)}\\\
&-|...,x_{a}-1,x_{a+1},...\rangle-|...,x_{a},x_{a+1}+1,...\rangle-(L-2n+1)|...,x_{i},...\rangle\\\
&-|...,x_{a+1},x_{a},...\rangle+2\Delta\delta^{a,a+1}|...,x_{a},x_{a+1},...\rangle\bigr{]}+\cdots\biggr{)}\\\
&=\sum_{Q}\sum_{1\leq x_{Q1}<...<x_{Qn}\leq
L}E\psi_{Q}|...\rangle,\end{split}$

where $\delta^{a,a+1}=1$ if $c_{Qa}=c_{Q(a+1)}$, and $\delta^{a,a+1}=0$ if
$c_{Qa}\neq c_{Q(a+1)}$, and we have omitted the unchanged variables in the
ket vectors, and the ellipsis in the sum on the l.h.s. denote terms with
multiple adjacent identical particles. Collecting coefficients of the same
basis vectors on both sides, we have

$-\sum_{i=1}^{n}[\psi_{Q}(\mathbf{x}-\mathbf{e_{i}})+\psi_{Q}(\mathbf{x}+\mathbf{e_{i}})]-(L-2n)\psi_{Q}(\mathbf{x})=E\psi_{Q}(\mathbf{x}),$
(82)

for the non-adjacent case, where $\mathbf{e_{i}}$ denotes the unit vector in
the $i$’th coordinate component, and

$\begin{split}&-\sum_{i=1}^{n}[\psi_{Q}(\mathbf{x}-\mathbf{e_{i}})+\psi_{Q}(\mathbf{x}+\mathbf{e_{i}})]+\psi_{Q}(\mathbf{x}+\mathbf{e_{a}})+\psi_{Q}(\mathbf{x}-\mathbf{e_{a+1}})-(L-2n+1)\psi_{Q}(\mathbf{x})\\\
&-\psi_{\tau_{a(a+1)}Q}(\mathbf{x}\tau_{a(a+1)})+2\Delta\delta^{a,a+1}\psi_{Q}(\mathbf{x})=E\psi_{Q}(\mathbf{x}),\end{split}$
(83)

for the adjacent case, where $\tau_{a(a+1)}$ denotes the transposition between
$a$ and $a+1$. The difference between the above two equations then gives the
boundary condition

$\psi_{Q}(\mathbf{x}+\mathbf{e_{a}})+\psi_{Q}(\mathbf{x}-\mathbf{e_{a+1}})-\psi_{Q}(\mathbf{x})-\psi_{\tau_{a(a+1)}Q}(\mathbf{x}\tau_{a(a+1)})+2\Delta\delta^{a,a+1}\psi_{Q}(\mathbf{x})=0.$
(84)

Plugging in our Bethe trial wave function for $\phi_{Q}$, we have

$\sum_{P\in S_{n}}[(\mu_{Pa}+\mu_{P(a+1)}^{-1}-1+\Delta
p_{Qa,Q(a+1)})A_{Q,P}-A_{\tau_{a(a+1)}Q,P}]\prod_{i=1}^{n}\mu_{Pi}^{q_{Qi}}=0.$
(85)

Denoting $Pa,P(a+1)$ with $Pa=i$, $P(a+1)=j$, we have

$\begin{split}&\sum_{P\in
S_{n}/Z_{2}}\biggl{(}\bigl{[}\mu_{i}\mu_{j}+1+(2\Delta\delta^{a,a+1}-1)\mu_{j}\bigr{]}A_{Q,P}-\mu_{j}A_{\tau_{a(a+1)}Q,P}\\\
&+\bigl{[}\mu_{i}\mu_{j}+1+(2\Delta\delta^{a,a+1}-1)\mu_{i}\bigr{]}A_{Q,P\sigma_{ij}}-\mu_{i}A_{\tau_{a(a+1)}Q,P\sigma_{ij}}\biggr{)}(\mu_{i}\mu_{j})^{a}\prod_{i\neq
a,a+1}\mu_{Pi}^{x_{Qi}}=0.\end{split}$ (86)

A sufficient condition for this equation to hold for any choice of
$\\{x_{i}\\}$, is that the coefficients in the sum vanish term by term.
Treating now $A_{Q,P}$ as components of $n!$ dimensional column vectors
$\xi_{P}$, this becomes

$\begin{split}&\bigl{(}\mu_{i}\mu_{j}+1+(2\Delta\delta^{a,a+1}-1)\mu_{j}-\pi(\tau_{a(a+1)})\mu_{j}\bigr{)}\xi_{P}\\\
+&\bigl{(}\mu_{i}\mu_{j}+1+(2\Delta\delta^{a,a+1}-1)\mu_{i}-\pi(\tau_{a(a+1)})\mu_{i}\bigr{)}\xi_{P\sigma_{ij}}=0,\end{split}$
(87)

where $\pi$ is the left regular representation of $S_{n}$, which is associated
to the scattering process of a Bethe wave function [82]. From this, we can
solve Yang’s scattering $Y$ matrix, as defined in (54), to be

$Y_{ij}^{a(a+1)}(\Delta)=R_{ij}^{a(a+1)}(\Delta)\mathbb{1}+T_{ij}^{a(a+1)}(\Delta)\pi(\tau_{a(a+1)}),$
(88)

where the reflection coefficient

$R_{ij}^{ab}(\Delta)=-\frac{(\mu_{i}\mu_{j}+1+2\Delta\delta^{ab}\mu_{j})(\mu_{i}\mu_{j}+1+2(\Delta\delta^{ab}-1)\mu_{i})+(\mu_{i}\mu_{j}+1)(\mu_{i}-\mu_{j})}{(\mu_{i}\mu_{j}+1+2\Delta\delta^{ab}\mu_{i})(\mu_{i}\mu_{j}+1+2(\Delta\delta^{ab}-1)\mu_{i})},$
(89)

and transmission coefficient

$T_{ij}^{ab}(\Delta)=-\frac{(\mu_{i}\mu_{j}+1)(\mu_{i}-\mu_{j})}{(\mu_{i}\mu_{j}+1+2\Delta\delta^{ab}\mu_{i})(\mu_{i}\mu_{j}+1+2(\Delta\delta^{ab}-1)\mu_{i})}.$
(90)

##  Appendix D Nested Bethe Ansatz diagonalization for the isotropic
Hamiltonian

The periodic boundary condition in the multiple identical species case relates
wave functions in different sectors of the coordinate space

$\psi_{Q}(x_{1},...,\\{q_{j}=1;c_{j}\\},...,x_{n})=\psi_{\tilde{\tau}Q}(x_{1},...,\\{q_{j}=L+1;c_{j}\\},...,x_{n}),$
(91)

where $\tilde{\tau}$ denotes the translation $(2\ 3\ \cdots\ n\ 1)$. Upon
applying the ansatz wave function, this implies

$A_{Q,P}=A_{\tilde{\tau}Q,\tilde{\tau}P}\mu_{P(1)}^{L},$ (92)

and consequently

$\pi(\tilde{\tau})\xi_{P}=\mu_{P(1)}^{L}\xi_{\tilde{\tau}P},$ (93)

since the left regular representation acts as
$\pi(\tilde{\tau})A_{Q,P}=A_{\tilde{\tau}^{-1}Q,P}$. Taking
$P^{0}=\tau_{12}\cdots\tau_{(j-1)j}$, $\varphi=\xi_{P^{0}}$, and using the
braiding property $S_{kj}\pi(\tau_{ik})=\pi(\tau_{ik})S_{ij}$, the periodic
boundary condition then leads to

$\mu_{j}^{L}\varphi=S_{(j+1)j}S_{(j+2)j}\cdots S_{nj}S_{1j}S_{2j}\cdots
S_{(j-1)j}\varphi,$ (94)

for $j=1,\cdots,n$. Unlike the single identical species case, this is no
longer a scalar equation, but an eigenvalue problem that requires
diagonalization of the r.h.s. To proceed from here, we have to specify an
irreducible representation of the permutation group $\mathfrak{S}_{n}$
corresponding to sector
$\\{c_{1}^{n_{1}},\cdots,c_{m}^{n_{m}},c_{m+1},\cdots,c_{L-n+m}\\}$ of the
Hilbert space. To find the lowest energy of this sector, let’s take the irrep
$R=[m^{n_{m}},(m-1)^{n_{m-1}-n_{m}},...,2^{n_{2}-n_{3}},1^{n_{1}-n_{2}}]$.
Sutherland’s approach of diagonalization of this irrep corresponding to the
multi-row Young tableaux using Bethe-Yang Hypothesis is to first treat all the
other $m-1$ species as the same, except the first one, then all the other
$m-2$ species as the same, except the second one, and so on, so forth. This
way, at each stage of the nesting, we are dealing with an irrep corresponding
to a 2-row Young tableaux, which can be solved by recourse to the algebraic
Bethe Ansatz as used in diagonalizing the spin-half Heisenberg model. Notice,
however, at each stage, while the irrep of the spin part of the wave function
forms a 2-row Young Tableaux, the spatial wave function must form the irrep
corresponding to the 2-column tableaux $[2^{n_{1}},1^{n-n_{1}}]$ at the first
stage, for instance. Yang’s prescription to this difficulty is to consider
instead of the eigenvalue problem of (94), the equivalent eigenvalue problem
of

$\mu_{j}^{L}\varphi=S^{\prime}_{(j+1)j}S^{\prime}_{(j+2)j}\cdots
S^{\prime}_{nj}S^{\prime}_{1j}S^{\prime}_{2j}\cdots
S^{\prime}_{(j-1)j}\varphi,$ (95)

where

$S^{\prime}_{ij}=T_{ij}\mathbb{1}-R_{ij}\tilde{\pi}(\tau_{ij})=\frac{\lambda_{i}-\lambda_{j}+i\pi_{\tilde{R}}(\tau_{ij})}{\lambda_{i}-\lambda_{j}+i}$
(96)

is written in terms of the conjugate representation
$\tilde{R}=[n_{1},n-n_{1}]$. Their equivalence is manifested by the fact that
$\pi_{\tilde{R}}(\tau_{ij})=-\pi_{R}(\tau_{ij})$. The advantage of adopting
such a conjugate representation is that it admits a realization in terms of
the the scattering matrix in the Heisenberg spin-$\frac{1}{2}$ problem of a
length $n$ chain. So we can readily apply the results from algebraic Bethe
Ansatz at each stage to write $\varphi$ as

$\varphi=\sum_{\sigma\in\mathfrak{S}_{n-n_{1}}}A_{\sigma}F(\Lambda_{\sigma
1},y_{1})F(\Lambda_{\sigma 2},y_{2})\cdots
F(\Lambda_{\sigma(n-n_{1})},y_{n-n_{1}}),$ (97)

where $y_{1}<y_{2}<\cdots<y_{n-n_{1}}$ are the coordinates of the $n_{1}$
“down spins”, and

$F(\Lambda,y)=\prod_{j=1}^{y-1}\frac{\lambda_{j}-\Lambda+i/2}{\lambda_{j+1}-\Lambda-i/2}$
(98)

are defined in terms of the set of unequal numbers to be solved from the set
of coupled algebraic equations

$\displaystyle\mu_{j}^{L}=$
$\displaystyle-\prod_{k=1}^{n}\frac{\lambda_{j}-\lambda_{k}-i}{\lambda_{j}-\lambda_{k}+i}\prod_{\alpha=1}^{n-n_{1}}\frac{\lambda_{j}-\Lambda_{\alpha}^{(1)}+i/2}{\lambda_{j}-\Lambda_{\alpha}^{(1)}-i/2},$
$\displaystyle j=1,...,n,$ (99)
$\displaystyle\prod_{j=1}^{n}\frac{\Lambda_{\alpha}^{(1)}-\lambda_{j}-i/2}{\Lambda_{\alpha}^{(1)}-\lambda_{j}+i/2}=$
$\displaystyle-\prod_{\beta=1}^{n-n_{1}}\frac{\Lambda_{\alpha}^{(1)}-\Lambda_{\beta}^{(1)}-i}{\Lambda_{\alpha}^{(1)}-\Lambda_{\beta}^{(1)}+i}\prod_{\gamma=1}^{n-n_{1}-n_{2}}\frac{\Lambda_{\alpha}^{(1)}-\Lambda_{\gamma}^{(2)}+i/2}{\Lambda_{\alpha}^{(1)}-\Lambda_{\gamma}^{(2)}-i/2},$
$\displaystyle\alpha=1,...,n-n_{1},$ (100) $\displaystyle\vdots$
$\displaystyle\prod_{\delta=1}^{n_{m-1}+n_{m}}\frac{\Lambda_{\zeta}^{(m)}-\Lambda_{\delta}^{(m-1)}-i/2}{\Lambda_{\zeta}^{(m)}-\Lambda_{\delta}^{(m-1)}+i/2}=$
$\displaystyle-\prod_{\epsilon=1}^{n_{m}}\frac{\Lambda_{\zeta}^{(m)}-\Lambda_{\epsilon}^{(m)}-i}{\Lambda_{\zeta}^{(m)}-\Lambda_{\epsilon}^{(m)}+i},$
$\displaystyle\zeta=1,...,n_{m}.$ (101)

This set of equations are historically called Lieb-Wu equations [83], the
details of their derivation of these equations can be found in standard texts
such as [82]. Now if there’s a cut-off on the number of components, the ground
state of our Hamiltonian will no longer be a superposition of fully
antisymmetrized spin configurations. In stead, from the reasoning above, it
will be in the sector where each color appear the same number of times along
the chain, say the number of colors $s$ divides $L$, and $l=L/s$. In this
sector, the Young tableaux corresponding to the lowest energy eigenstate in
this sector will be $[l^{s}]$. The Lieb-Wu equations for this irrep become

$\displaystyle L\theta(\lambda_{j})=$ $\displaystyle 2\pi
J_{j}^{(0)}+\sum_{k=1}^{L}\theta\big{(}\frac{1}{2}(\lambda_{k}-\lambda_{j}\big{)}+\sum_{\alpha=1}^{(s-1)l}\theta(\lambda_{j}-\Lambda^{(1)}_{\alpha}),$
$\displaystyle j=1,...,L,$ (102)
$\displaystyle\sum_{j=1}^{L}\theta(\lambda_{j}-\Lambda^{(1)}_{\alpha})=$
$\displaystyle 2\pi
J_{\alpha}^{(1)}+\sum_{\beta=1}^{(s-1)l}\theta\big{(}\frac{1}{2}(\Lambda^{(1)}_{\beta}-\Lambda^{(1)}_{\alpha})\big{)}+\sum_{\gamma=1}^{(s-2)l}\theta(\Lambda^{(1)}_{\alpha}-\Lambda^{(2)}_{\gamma}),$
$\displaystyle\alpha=1,...,(s-1)l,$ (103) $\displaystyle\vdots$
$\displaystyle\sum_{\delta=1}^{2l}\theta(\Lambda_{\delta}^{(s-1)}-\Lambda_{\zeta}^{(s)})=$
$\displaystyle 2\pi
J^{(s)}_{\zeta}+\sum_{\epsilon=1}^{l}\theta\big{(}\frac{1}{2}(\Lambda^{(s)}_{\epsilon}-\Lambda^{(s)}_{\zeta})\big{)},$
$\displaystyle\zeta=1,...,l,$ (104)

where $\theta(\lambda)=2\cot^{-1}(2\lambda)$, and the quantum numbers that
labels the eigenstates, $J$’s, are half-integers.

## References

  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * Gorshkov  _et al._ [2010] A. V. Gorshkov, M. Hermele, V. Gurarie, C. Xu, P. S. Julienne, J. Ye, P. Zoller, E. Demler, M. D. Lukin, and A. M. Rey, Two-orbital s u(n) magnetism with ultracold alkaline-earth atoms, Nature Physics 6, 289 (2010).
  * Taie  _et al._ [2012] S. Taie, R. Yamazaki, S. Sugawa, and Y. Takahashi, An su(6) mott insulator of an atomic fermi gas realized by large-spin pomeranchuk cooling, Nature Physics 8, 825 (2012).
  * Pagano  _et al._ [2014] G. Pagano, M. Mancini, G. Cappellini, P. Lombardi, F. Schäfer, H. Hu, X.-J. Liu, J. Catani, C. Sias, M. Inguscio, and L. Fallani, A one-dimensional liquid of fermions with tunable spin, Nature Physics 10, 198 (2014).
  * Scazza  _et al._ [2014] F. Scazza, C. Hofrichter, M. Höfer, P. C. De Groot, I. Bloch, and S. Fölling, Observation of two-orbital spin-exchange interactions with ultracold su(n)-symmetric fermions, Nature Physics 10, 779 (2014).
  * Zhang  _et al._ [2014] X. Zhang, M. Bishof, S. L. Bromley, C. V. Kraus, M. S. Safronova, P. Zoller, A. M. Rey, and J. Ye, Spectroscopic observation of su(n)-symmetric interactions in sr orbital magnetism, Science 345, 1467 (2014).
  * Yang [1967] C. N. Yang, Some exact results for the many-body problem in one dimension with repulsive delta-function interaction, Phys. Rev. Lett. 19, 1312 (1967).
  * Sutherland [1968] B. Sutherland, Further results for the many-body problem in one dimension, Phys. Rev. Lett. 20, 98 (1968).
  * Sutherland [1975] B. Sutherland, Model for a multicomponent quantum system, Phys. Rev. B 12, 3795 (1975).
  * Babelon  _et al._ [1983] O. Babelon, H. J. de Vega, and C. M. Viallet, Exact excitation spectrum of the zn + 1 ×zn + 1 generalized heisenberg model, Nuclear Physics B 220, 283 (1983).
  * Mussardo  _et al._ [2020] G. Mussardo, A. Trombettoni, and Z. Zhang, Prime suspects in a quantum ladder, Phys. Rev. Lett. 125, 240603 (2020).
  * Kiwata and Akutsu [1994] H. Kiwata and Y. Akutsu, Bethe-ansatz approximation for the s=1 antiferromagnetic spin chain, Journal of the Physical Society of Japan 63, 3598 (1994).
  * Sato [1995] R. Sato, Bethe ansatz calculations for quantum spin chains with partial integrability, Journal of the Physical Society of Japan 64, 2837 (1995).
  * Sato and Akutsu [1996] R. Sato and Y. Akutsu, First-order phase transition in the s=1 quantum spin chain with partial integrability, Journal of the Physical Society of Japan 65, 1885 (1996).
  * Ushveridze [1994] A. G. Ushveridze,  _Quasi-exactly solvable models in quantum mechanics_ (Taylor and Francis Group, 1994).
  * Majumdar and Ghosh [1969] C. K. Majumdar and D. K. Ghosh, On next‐nearest‐neighbor interaction in linear chain. i,  _Journal of Mathematical Physics_ , Journal of Mathematical Physics 10, 1388 (1969).
  * Affleck  _et al._ [1987] I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Rigorous results on valence-bond ground states in antiferromagnets, Phys. Rev. Lett. 59, 799 (1987).
  * Bravyi  _et al._ [2012] S. Bravyi, L. Caha, R. Movassagh, D. Nagaj, and P. W. Shor, Criticality without frustration for quantum spin-1 chains, Phys. Rev. Lett. 109, 207202 (2012).
  * Movassagh and Shor [2016] R. Movassagh and P. W. Shor, Supercritical entanglement in local systems: Counterexample to the area law for quantum matter, Proceedings of the National Academy of Sciences 113, 13278 (2016), https://www.pnas.org/content/113/47/13278.full.pdf .
  * Zhang  _et al._ [2017] Z. Zhang, A. Ahmadain, and I. Klich, Novel quantum phase transition from bounded to extensive entanglement, Proceedings of the National Academy of Sciences 114, 5142 (2017), https://www.pnas.org/content/114/20/5142.full.pdf .
  * Salberger and Korepin [2017] O. Salberger and V. Korepin, Entangled spin chain,  _Reviews in Mathematical Physics_ , Reviews in Mathematical Physics 29, 1750031 (2017).
  * Salberger  _et al._ [2017] O. Salberger, T. Udagawa, Z. Zhang, H. Katsura, I. Klich, and V. Korepin, Deformed fredkin spin chain with extensive entanglement, Journal of Statistical Mechanics: Theory and Experiment 2017, 063103 (2017).
  * Zhang and Klich [2017] Z. Zhang and I. Klich, Entropy, gap and a multi-parameter deformation of the fredkin spin chain, Journal of Physics A: Mathematical and Theoretical 50, 425201 (2017).
  * Kitaev [2003] A. Y. Kitaev, Fault-tolerant quantum computation by anyons, Annals of Physics 303, 2 (2003).
  * Caspers and Magnus [1982] W. J. Caspers and W. Magnus, Some exact excited states in a linear antiferromagnetic spin system, Physics Letters A 88, 103 (1982).
  * Arovas [1989] D. P. Arovas, Two exact excited states for the s = 1 aklt chain, Physics Letters A 137, 431 (1989).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Fendley  _et al._ [2004] P. Fendley, K. Sengupta, and S. Sachdev, Competing density-wave orders in a one-dimensional hard-boson model, Phys. Rev. B 69, 075106 (2004).
  * Lesanovsky and Katsura [2012] I. Lesanovsky and H. Katsura, Interacting fibonacci anyons in a rydberg gas, Phys. Rev. A 86, 041601 (2012).
  * Turner  _et al._ [2018] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Weak ergodicity breaking from quantum many-body scars, Nature Physics 14, 745 (2018).
  * Schecter and Iadecola [2019] M. Schecter and T. Iadecola, Weak ergodicity breaking and quantum many-body scars in spin-1 $xy$ magnets, Phys. Rev. Lett. 123, 147201 (2019).
  * Moudgalya  _et al._ [2018a] S. Moudgalya, S. Rachel, B. A. Bernevig, and N. Regnault, Exact excited states of nonintegrable models, Phys. Rev. B 98, 235155 (2018a).
  * Iadecola and Schecter [2020] T. Iadecola and M. Schecter, Quantum many-body scar states with emergent kinetic constraints and finite-entanglement revivals, Phys. Rev. B 101, 024306 (2020).
  * Chattopadhyay  _et al._ [2020] S. Chattopadhyay, H. Pichler, M. D. Lukin, and W. W. Ho, Quantum many-body scars from virtual entangled pairs, Phys. Rev. B 101, 174308 (2020).
  * Shibata  _et al._ [2020] N. Shibata, N. Yoshioka, and H. Katsura, Onsager’s scars in disordered spin chains, Phys. Rev. Lett. 124, 180604 (2020).
  * Moudgalya  _et al._ [2020] S. Moudgalya, E. O’Brien, B. A. Bernevig, P. Fendley, and N. Regnault, Large classes of quantum scarred hamiltonians from matrix product states, (2020), arXiv:2002.11725 [cond-mat.str-el] .
  * Shiraishi and Mori [2017] N. Shiraishi and T. Mori, Systematic construction of counterexamples to the eigenstate thermalization hypothesis, Phys. Rev. Lett. 119, 030601 (2017).
  * Ok  _et al._ [2019] S. Ok, K. Choo, C. Mudry, C. Castelnovo, C. Chamon, and T. Neupert, Topological many-body scar states in dimensions one, two, and three, Phys. Rev. Research 1, 033144 (2019).
  * Lee  _et al._ [2020] K. Lee, R. Melendrez, A. Pal, and H. J. Changlani, Exact three-colored quantum scars from geometric frustration, Phys. Rev. B 101, 241111 (2020).
  * Surace  _et al._ [2020] F. M. Surace, G. Giudici, and M. Dalmonte, Weak-ergodicity-breaking via lattice supersymmetry, (2020), arXiv:2003.11073 [cond-mat.stat-mech] .
  * Vafek  _et al._ [2017] O. Vafek, N. Regnault, and B. A. Bernevig, Entanglement of Exact Excited Eigenstates of the Hubbard Model in Arbitrary Dimension, SciPost Phys. 3, 043 (2017).
  * Choi  _et al._ [2019] S. Choi, C. J. Turner, H. Pichler, W. W. Ho, A. A. Michailidis, Z. Papić, M. Serbyn, M. D. Lukin, and D. A. Abanin, Emergent su(2) dynamics and perfect quantum many-body scars, Phys. Rev. Lett. 122, 220603 (2019).
  * Schindler  _et al._ [2022] F. Schindler, N. Regnault, and B. A. Bernevig, Exact quantum scars in the chiral nonlinear luttinger liquid, Phys. Rev. B 105, 035146 (2022).
  * Mark  _et al._ [2020] D. K. Mark, C.-J. Lin, and O. I. Motrunich, Unified structure for exact towers of scar states in the affleck-kennedy-lieb-tasaki and other models, Phys. Rev. B 101, 195131 (2020).
  * Bravyi and Hastings [2011] S. Bravyi and M. B. Hastings, A short proof of stability of topological order under local perturbations, Communications in Mathematical Physics 307, 609 (2011).
  * Gosset and Huang [2016] D. Gosset and Y. Huang, Correlation length versus gap in frustration-free systems, Phys. Rev. Lett. 116, 097202 (2016).
  * Gosset and Mozgunov [2016] D. Gosset and E. Mozgunov, Local gap threshold for frustration-free spin systems, Journal of Mathematical Physics 57, 091901 (2016).
  * Lemm and Mozgunov [2019] M. Lemm and E. Mozgunov, Spectral gaps of frustration-free spin systems with boundary, Journal of Mathematical Physics 60, 051901 (2019).
  * Anshu  _et al._ [2020] A. Anshu, I. Arad, and D. Gosset, Entanglement subvolume law for 2d frustration-free spin systems, in  _Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing_ (Association for Computing Machinery, New York, NY, USA, 2020) pp. 868–874.
  * Nachtergaele  _et al._ [2021] B. Nachtergaele, R. Sims, and A. Young, Stability of the bulk gap for frustration-free topologically ordered quantum lattice systems (2021).
  * Arovas [2008] D. P. Arovas, Simplex solid states of $\mathrm{SU}(n)$ quantum antiferromagnets, Phys. Rev. B 77, 104404 (2008).
  * Greiter  _et al._ [2007] M. Greiter, S. Rachel, and D. Schuricht, Exact results for su(3) spin chains: Trimer states, valence bond solids, and their parent hamiltonians, Phys. Rev. B 75, 060401 (2007).
  * Greiter and Rachel [2007] M. Greiter and S. Rachel, Valence bond solids for $\mathrm{SU}(n)$ spin chains: Exact models, spinon confinement, and the haldane gap, Phys. Rev. B 75, 184441 (2007).
  * Laflorencie [2016] N. Laflorencie, Quantum entanglement in condensed matter systems,  _Quantum entanglement in condensed matter systems_ , Physics Reports 646, 1 (2016).
  * Hastings [2007] M. B. Hastings, An area law for one-dimensional quantum systems,  _Journal of Statistical Mechanics: Theory and Experiment_ ,  2007, P08024 (2007).
  * Vitagliano  _et al._ [2010] G. Vitagliano, A. Riera, and J. I. Latorre, Volume-law scaling for the entanglement entropy in spin-1/2 chains,  _New Journal of Physics_ ,  12, 113049 (2010).
  * Ramírez  _et al._ [2014] G. Ramírez, J. Rodríguez-Laguna, and G. Sierra, From conformal to volume law for the entanglement entropy in exponentially deformed critical spin 1/2 chains,  _Journal of Statistical Mechanics: Theory and Experiment_ ,  2014, P10004 (2014).
  * Alba  _et al._ [2009] V. Alba, M. Fagotti, and P. Calabrese, Entanglement entropy of excited states, Journal of Statistical Mechanics: Theory and Experiment 2009, P10020 (2009).
  * Chen  _et al._ [2013] P. Chen, Z. long Xue, I. P. McCulloch, M.-C. Chung, M. Cazalilla, and S.-K. Yip, Entanglement entropy scaling of the XXZ chain, Journal of Statistical Mechanics: Theory and Experiment 2013, P10007 (2013).
  * Hamermesh [1962] M. Hamermesh,  _Group Theory and its Application to Physical Problems_ (Addison-Wesley Publishing Company, INC., 1962).
  * Elliot and Dawber [1979] J. P. Elliot and P. G. Dawber,  _Symmetry in Physics, Volumn 2: Further Applications_ (Macmillan Publishers LTD, 1979).
  * Nataf and Mila [2014] P. Nataf and F. Mila, Exact diagonalization of heisenberg $\mathrm{SU}(n)$ models, Phys. Rev. Lett. 113, 127204 (2014).
  * Nataf and Mila [2016] P. Nataf and F. Mila, Exact diagonalization of heisenberg $\mathrm{SU}(n)$ chains in the fully symmetric and antisymmetric representations, Phys. Rev. B 93, 155134 (2016).
  * Gaudin [1967] M. Gaudin, Un systeme a une dimension de fermions en interaction, Physics Letters A 24, 55 (1967).
  * Gaudin [2014] M. Gaudin,  _The Bethe Wavefunction_ (Cambridge University Press, 2014).
  * Note [1] Although Sutherland was actually solving the ferromagnetic model for fermionic particles, he showed the spectrum is that of the antiferromagnetic bosonic one in reversed ordering.
  * Yang and Yang [1966] C. N. Yang and C. P. Yang, Ground-state energy of a heisenberg-ising lattice, Phys. Rev. 147, 303 (1966).
  * McGuire [1965] J. B. McGuire, Interacting fermions in one dimension. i. repulsive potential,  _Journal of Mathematical Physics_ , Journal of Mathematical Physics 6, 432 (1965).
  * McGuire [1966] J. B. McGuire, Interacting fermions in one dimension. ii. attractive potential,  _Journal of Mathematical Physics_ , Journal of Mathematical Physics 7, 123 (1966).
  * Flicker and Lieb [1967] M. Flicker and E. H. Lieb, Delta-function fermi gas with two-spin deviates, Phys. Rev. 161, 179 (1967).
  * Lieb and Mattis [1962] E. Lieb and D. Mattis, Theory of ferromagnetism and the ordering of electronic energy levels, Phys. Rev. 125, 164 (1962).
  * Lieb [1989] E. H. Lieb, Two theorems on the hubbard model, Phys. Rev. Lett. 62, 1201 (1989).
  * Hakobyan [2004] T. Hakobyan, The ordering of energy levels for su(n) symmetric antiferromagnetic chains, Nuclear Physics B 699, 575 (2004).
  * Sutherland [2004] B. Sutherland,  _Beautiful Models_ (WORLD SCIENTIFIC, 2004).
  * Moudgalya  _et al._ [2018b] S. Moudgalya, N. Regnault, and B. A. Bernevig, Entanglement of exact excited states of affleck-kennedy-lieb-tasaki models: Exact results, many-body scars, and violation of the strong eigenstate thermalization hypothesis, Phys. Rev. B 98, 235156 (2018b).
  * McGuire and Hurst [1972] J. B. McGuire and C. A. Hurst, The scattering of three impenetrable particles in one dimension,  _Journal of Mathematical Physics_ , Journal of Mathematical Physics 13, 1595 (1972).
  * McGuire and Hurst [1988] J. B. McGuire and C. A. Hurst, Three interacting particles in one dimension: An algebraic approach,  _Journal of Mathematical Physics_ , Journal of Mathematical Physics 29, 155 (1988).
  * Bibikov [2016] P. N. Bibikov, Three magnons in an isotropics = 1 ferromagnetic chain as an exactly solvable non-integrable system,  _Journal of Statistical Mechanics: Theory and Experiment_ ,  2016, 033109 (2016).
  * McGuire [1964] J. B. McGuire, Study of exactly soluble one‐dimensional n‐body problems,  _Journal of Mathematical Physics_ , Journal of Mathematical Physics 5, 622 (1964).
  * Lamacraft [2013] A. Lamacraft, Diffractive scattering of three particles in one dimension: A simple result for weak violations of the yang-baxter equation, Phys. Rev. A 87, 012707 (2013).
  * Orbach [1958] R. Orbach, Linear antiferromagnetic chain with anisotropic coupling, Phys. Rev. 112, 309 (1958).
  * Arutyunov [2019] G. Arutyunov, Elements of classical and quantum integrable systems, UNITEXT for Physics 10.1007/978-3-030-24198-8 (2019).
  * Lieb and Wu [1968] E. H. Lieb and F. Y. Wu, Absence of mott transition in an exact solution of the short-range, one-band model in one dimension, Phys. Rev. Lett. 20, 1445 (1968).