††thanks:  * These authors contributed equally to this work.

# Exact Thermal Eigenstates of Nonintegrable Spin Chains at Infinite
Temperature

Yuuya Chiba*<EMAIL_ADDRESS>Department of Basic Science, The
University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan  Yasushi
Yoneta*<EMAIL_ADDRESS>Center for Quantum Computing, RIKEN, 2-1
Hirosawa, Wako, Saitama 351-0198, Japan

###### Abstract

The eigenstate thermalization hypothesis (ETH) plays a major role in
explaining thermalization of isolated quantum many-body systems. However,
there has been no proof of the ETH in realistic systems due to the difficulty
in the theoretical treatment of thermal energy eigenstates of nonintegrable
systems. Here, we write down analytically, for the first time, thermal
eigenstates of nonintegrable spin chains. We consider a class of theoretically
tractable volume-law states, which we call entangled antipodal pair (EAP)
states. These states are thermal, in the most strict sense that they are
indistinguishable from the Gibbs state with respect to all local observables,
with infinite temperature. We then identify Hamiltonians having the EAP state
as an eigenstate and rigorously show that some of these Hamiltonians are
nonintegrable. Furthermore, a thermal pure state at an arbitrary temperature
is obtained by the imaginary time evolution of an EAP state. Our results offer
a potential avenue for providing a provable example of the ETH.

Introduction.— Understanding the mechanism of thermalization in quantum many-
body systems has been a pivotal issue in statistical physics [1, 2, 3].
Notably, the eigenstate thermalization hypothesis (ETH) [4, 5, 6] has served
as a cornerstone in this field. It posits that all the energy eigenstates of
quantum many-body systems exhibit thermal properties, thereby giving a
plausible explanation of thermalization.

While the ETH is anticipated to hold in most nonintegrable systems, the
verification of whether this hypothesis holds in realistic many-body systems
relies on numerical calculations, and a theoretical verification has remained
elusive [7, 8, 9, 10]. Thus, a significant challenge lies in theoretically
addressing the nature of energy eigenstates, particularly in nonintegrable
systems. However, it has not been clear whether thermal eigenstates of
nonintegrable systems can be treated theoretically. This stems from the
difficulty of writing down quantum states whose entanglement entropy obeys a
volume law.

One approach to treat quantum many-body states theoretically is to use
variational wave functions. Particularly for states that contain a small
amount of entanglement, they can be represented via tensor network states such
as a matrix product state (MPS) [11, 12]. Tensor network states are highly
tractable, making them not only practical but also significantly contributing
to theoretical advancements. Indeed, by utilizing the MPS, it has been
successful to exactly describe finite-energy-density low-entangled (thus
nonthermal) eigenstates even for nonintegrable systems [13, 14, 15], which are
examples of many-body scars [16, 17, 18, 19, 20]. However, there has been a
lack of variational wave functions suitable for theoretical analysis of
volume-law states, which is one of the reasons why thermal eigenstates have
not yet been obtained. Hence, there is a craving for a class of volume-law
states amenable to the theoretical treatment [21, 22, 23, 24, 25].

In this Letter, we provide, for the first time, pairs of a nonintegrable
Hamiltonian and its thermal eigenstate at infinite temperature. We consider a
class of volume-law states, which we call the entangled antipodal pair (EAP)
states, that are amenable to theoretical calculations. Then we fully
characterize Hamiltonians having the EAP state as an eigenstate. It is
rigorously shown that some of these Hamiltonians are nonintegrable. In
addition, by evolving an EAP state in imaginary time, we construct a thermal
pure state at arbitrary temperature, which is locally indistinguishable from
the Gibbs state.

Figure 1: Schematic diagram depicting an entangled antipodal pair state
$\ket{\mathrm{EAP}}$. The antipodal pairs of spins linked by dotted lines are
in the Bell states $\ket{\Phi_{pq}}$. For any subsystem with diameter smaller
than or equal to half of the size of the entire system $N$, the reduced
density matrix coincides with the maximally mixed state, i.e., the Gibbs state
$\hat{\rho}^{\mathrm{can}}$ at the inverse temperature $\beta=0$.

Entangled antipodal pair state.— We consider quantum spin-$1/2$ systems on the
one-dimensional lattice $\Lambda=\\{1,2,\cdots,N\\}$ with periodic boundary
conditions. We assume that the number of lattice sites $N$ is even. Let
$\hat{\sigma}_{j}^{\mu}(\mu=x,y,z)$ be the Pauli matrices acting on the $j$-th
site, and $\ket{0}_{j}$ and $\ket{1}_{j}$ be the eigenvectors of
$\hat{\sigma}_{j}^{z}$.

Our goal is to obtain pairs of a nonintegrable Hamiltonian and its thermal
eigenstate. To achieve this, we adopt the following strategy. First, we
introduce a class of volume-law states that are theoretically tractable. Next,
for each of these volume law states, we search for Hamiltonians that have it
as a thermal eigenstate. Finally, we prove that some of these Hamiltonians are
nonintegrable. Following this strategy, as the first step, we introduce the
entangled antipodal pair (EAP) state as 111Some special EAP states appear in
the study of integrable systems [38, 39, 40].

$\displaystyle\ket{\mathrm{EAP}}=\bigotimes_{j=1}^{N/2}\ket{\Phi_{p_{j}q_{j}}}_{j,j+N/2}.$
(1)

Here $\ket{\Phi_{p_{j}q_{j}}}_{j,j+N/2}(p_{j},q_{j}=0,1)$ are the Bell states
between sites $j$ and $j+N/2$ defined by

$\displaystyle\ket{\Phi_{pq}}_{j,j^{\prime}}=\frac{\ket{0}_{j}\ket{p}_{j^{\prime}}+(-1)^{q}\ket{1}_{j}\ket{{\bar{p}}}_{j^{\prime}}}{\sqrt{2}}\quad(p,q=0,1),$
(2)

where $\bar{p}$ represents the negation of $p$. It is straightforward to check
that

$\displaystyle\hat{\sigma}_{j}^{\mu}\hat{\sigma}_{j+N/2}^{\mu}\ket{\Phi_{p_{j}q_{j}}}_{j,j+N/2}=\omega_{j}^{\mu}\ket{\Phi_{p_{j}q_{j}}}_{j,j+N/2},$
(3)

where $\omega_{j}^{z}=(-1)^{p_{j}}$, $\omega_{j}^{x}=(-1)^{q_{j}}$, and
$\omega_{j}^{y}=-\omega_{j}^{x}\omega_{j}^{z}$. Hence, the EAP state is
uniquely characterized by $(\omega_{j}^{x})_{j\in\Lambda}$ and
$(\omega_{j}^{z})_{j\in\Lambda}$. As a consequence of Eq. (3), the action of
$\hat{\sigma}_{j}^{\mu}$ on the EAP state is equivalent to the action of
$\hat{\sigma}_{j+N/2}^{\mu}$, except for the factor of $\omega_{j}^{\mu}$,
i.e.,

$\displaystyle\hat{\sigma}_{j}^{\mu}\ket{\mathrm{EAP}}=\omega_{j}^{\mu}\hat{\sigma}_{j+N/2}^{\mu}\ket{\mathrm{EAP}}.$
(4)

This plays a key role in the proof of our main results.

Characteristically, the reduced density matrix of the EAP state for any
subsystem $X(\subset\Lambda)$ with diameter $D(X)=\max_{j,j^{\prime}\in
X}|j-j^{\prime}|$ of $N/2$ or less equals the maximally mixed state
$\propto\hat{I}_{X}$. Therefore, since the maximally mixed state coincides
with the Gibbs state $\hat{\rho}^{\mathrm{can}}\propto e^{-\beta\hat{H}}$ at
the inverse temperature $\beta=0$, it follows that the EAP state cannot be
distinguished from the thermal equilibrium state at $\beta=0$ by any local
measurement. Moreover, this implies that the bipartite entanglement between
the subsystem consisting of contiguous spins and its complement is maximal.
That is, the entanglement entropy of the EAP state obeys a volume low, with a
coefficient reaching the maximum value $\log 2$. This should be contrasted
with the rainbow state [21, 22, 23, 24], which is a product of the Bell states
between sites $j$ and $N-j+1$. For a subsystem $\\{1,2,\cdots,\ell\\}$, the
entanglement entropy of the rainbow state is $\ell\log 2$ and follows a volume
law. However, for a subsystem
$\\{N/2-\ell/2+1,N/2-\ell/2+2,\cdots,N/2+\ell/2\\}$, it strictly equals $0$
and follows an area law. As can be seen from this, the rainbow state is an
athermal state that can be distinguished from the Gibbs state with an
appropriate local observable.

EAP state as an energy eigenstate.— As explained above, in the following, we
search for Hamiltonians that have the EAP state as a thermal eigenstate. To
this end, we consider the most general form of the Hamiltonian written as

$\displaystyle\hat{H}$
$\displaystyle=\sum_{X(\subset\Lambda)}\sum_{\vec{\mu}\in\\{x,y,z\\}^{X}}J_{X}^{\vec{\mu}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}.$ (5)

In the first sum, $X$ represents a subset of $\Lambda$. In the second sum,
$\vec{\mu}$ represents a combination of $\mu_{j}=x,y,z$ for $j\in X$. Here the
origin of the energy is taken such that $\hat{H}$ is traceless. Now we impose
a very mild condition on the locality of interactions: $J_{X}^{\vec{\mu}}=0$
for any subset $X$ with $D(X)\geq N/4$. (For instance, a nearest neighbor
interacting system, where $J_{X}^{\vec{\mu}}=0$ for $D(X)\geq 3$, obviously
satisfies the above condition for large $N$.) We are interested in the
Hamiltonian where some EAP state $\ket{\mathrm{EAP}}$ is an energy eigenstate.
Such Hamiltonians can be completely characterized by the following theorem
222See Supplemental Material for proofs of theorems, additional numerical
results, and discussions..:

######  Theorem 1.

The following three statements are equivalent:

  1. (i)

An EAP state $\ket{\mathrm{EAP}}$ is an eigenstate of $\hat{H}$.

  2. (ii)

An EAP state $\ket{\mathrm{EAP}}$ is an eigenstate of $\hat{H}$ with the
eigenvalue $0$.

  3. (iii)

For all $X\subset\Lambda$ and $\vec{\mu}\in\\{x,y,z\\}^{X}$,

$\displaystyle J_{Y}^{\vec{\nu}}=-J_{X}^{\vec{\mu}}\prod_{j\in
X}\omega_{j}^{\mu_{j}},$ (6)

where $Y=X+N/2$ and $\nu_{j}=\mu_{j-N/2}$ for $j\in Y$.

Note that Eq. (6) has many solutions. In fact, for any given EAP state, one
can construct many trivial Hamiltonians that consist of only a few local
terms, e.g., $\hat{H}=J_{X}^{\vec{\mu}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}+J_{Y}^{\vec{\nu}}\bigotimes_{j\in
Y}\hat{\sigma}_{j}^{\nu_{j}}$ with $J_{X}^{\vec{\mu}}$ and $J_{Y}^{\vec{\nu}}$
satisfying Eq. (6). However, in most cases, such trivial Hamiltonians will be
outside the scope of statistical mechanics. Therefore, in the following, we
focus on cases where the Hamiltonian is translation invariant and show that
only several EAP states are allowed as solutions of Eq. (6).

Translation-invariant nonintegrable Hamiltonians.— Let us consider
translation-invariant Hamiltonians with nearest neighbor interactions. They
can include only $9$ interaction terms, $J^{\mu\nu}$, and $3$ magnetic field
terms, $h^{\mu}$. Solving Eq. (6) in Theorem 1, we can find all solutions of
$J^{\mu\nu}$, $h^{\mu}$ and
$(\omega_{j}^{x},\omega_{j}^{z},\omega_{j}^{y}=-\omega_{j}^{x}\omega_{j}^{z})_{j\in\Lambda}$
[27]. Interestingly, there are some nontrivial solutions whose
$(\omega_{j}^{x},\omega_{j}^{z},\omega_{j}^{y})_{j\in\Lambda}$ are not
invariant by single-site shift but invariant by $n$-site shift for $n>1$. To
express such solutions efficiently, we introduce the following notation for
EAP states that are invariant by $n$-site shift 333Note that, for some of
states (7), $n$-site shift can cause a phase change, such as
$\hat{\mathcal{T}}^{3}\ket{3;10,11,01}=-\ket{3;10,11,01}$, where
$\hat{\mathcal{T}}$ is the one-site translation operator.: When $N/2$ is a
multiple of $n$,

$\displaystyle\ket{n;p_{1}^{*}q_{1}^{*},\cdots,p_{n}^{*}q_{n}^{*}}$ (7)

denotes the EAP state characterized by $(p_{j},q_{j})_{j\in\Lambda}$
satisfying $p_{mn+j}=p_{j}^{*}$ and $q_{mn+j}=q_{j}^{*}$ for
$m=0,1,\cdots,N/2n-1$ and $j=1,2,\cdots,n$. For example, when $N/2$ is even,
$\ket{2;10,11}=\bigotimes_{j=1}^{N/4}\ket{\Phi_{10}}_{2j-1,2j-1+N/2}\bigotimes_{j=1}^{N/4}\ket{\Phi_{11}}_{2j,2j+N/2}$.
Using this notation, we can list all nontrivial solutions as follows:

######  Theorem 2.

By excluding noninteracting Hamiltonians 444There are only two solutions whose
Hamiltonians are noninteracting [27]., the solution of Eq. (6) whose
Hamiltonian is translation invariant and includes only nearest neighbor
interactions and magnetic fields is restricted to the following:

  1. 1.

The EAP state $\ket{1;00}$ is an eigenstate of

$\displaystyle\hat{H}_{1}=\sum_{j=1}^{N}\bigl{(}$ $\displaystyle
J^{xy}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{y}+J^{yx}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{x}$
$\displaystyle+J^{yz}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{z}+J^{zy}\hat{\sigma}_{j}^{z}\hat{\sigma}_{j+1}^{y}+h^{y}\hat{\sigma}_{j}^{y}\bigr{)}$
(8)

for arbitrary values of $J^{xy},J^{yx},J^{yz},J^{zy},h^{y}$.

  2. 2.

When $N/2$ is a multiple of three, the EAP state $\ket{3;10,11,01}$ is an
eigenstate of

$\displaystyle\hat{H}_{2}=\sum_{j=1}^{N}\bigl{(}J^{xy}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{y}+J^{yz}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{z}\bigr{)}$
(9)

for arbitrary values of $J^{xy},J^{yz}$.

  3. 3.

When $N/2$ is a multiple of four, the EAP state $\ket{4;10,11,00,01}$ is an
eigenstate of

$\displaystyle\hat{H}_{3}=\sum_{j=1}^{N}\bigl{(}J^{xx}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{x}+J^{yz}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{z}\bigr{)}$
(10)

for arbitrary values of $J^{xx},J^{yz}$.

and their equivalents obtained by appropriate permutations of directions of
the Pauli matrices.

As shown in this theorem, there exist not so many solutions in the case of
translation-invariant Hamiltonians, and this implies that our condition Eq.
(6) is strong enough. Considering the importance of these models, we call the
models described by $\hat{H}_{1}$, $\hat{H}_{2}$ and $\hat{H}_{3}$ Models 1, 2
and 3, respectively. Note that, since Model 2 is included in Model 1, the EAP
state $\ket{1;00}$ is also an eigenstate of $\hat{H}_{2}$.

Furthermore, we can obtain the following theorem regarding the
nonintegrability of Models 2 and 3 [27]:

######  Theorem 3.

For Model 2 with $J^{xy},J^{yz}\neq 0$ and for Model 3 with $J^{xx},J^{yz}\neq
0$, there exists no local conserved quantity other than a linear combination
of the identity and the Hamiltonian (i.e., trivial one).

Because it is known that integrable systems have many [$\mathcal{O}(N)$ number
of] nontrivial local conserved quantities [30, 31, 32], the above theorem
implies that these models are indeed nonintegrable. In addition, since Model 2
(which is $J^{yx}=J^{zy}=h^{y}=0$ case of Model 1) is nonintegrable for any
nonzero $J^{xy},J^{yz}$, Model 1 will be nonintegrable, at least for
nonaccidental values of the parameters [27].

Combining Theorems 2 and 3 with the fact that EAP states are thermal as
explained above, we first obtain, to our knowledge, an analytic expression of
a thermal energy eigenstate of a nonintegrable system. This highly contrasts
with the recent progress made by H. Tasaki [33] in attempts to prove the ETH.
He proved, only with respect to special observables, that all energy
eigenstates of a certain noninteracting (integrable) system are
indistinguishable from the thermal state. However, problems in nonintegrable
systems remained elusive in his work. In addition, the restriction on
observables is essential in his result because the integrable system has many
local observables by which energy eigenstates can be distinguished from the
thermal state. By contrast, our approach is to construct, in a nonintegrable
system, an energy eigenstate that is indistinguishable from a thermal state
with respect to any local observables.

Finite-temperature state.— So far, we have only considered states at
$\beta=0$. Can we describe thermal states at $\beta\neq 0$ using the EAP
state? Here, we provide a method for constructing finite-temperature states
using the EAP state. It should be noted that Hamiltonians considered here are
not restricted to those having the EAP state as an eigenstate, which we have
considered thus far.

Consider a system with translationally invariant short-range interactions.
Suppose that $\hat{H}$ is a real matrix with respect to the products of
eigenstates of $\hat{\sigma}_{j}^{z}$, $\\{\ket{00\cdots 0},\ket{10\cdots
0},\ket{01\cdots 0},\cdots\\}$. This class includes many well studied models
in statistical mechanics, such as the Ising model and the Heisenberg model.

Here, we utilize the EAP state characterized by $p_{j}=q_{j}=0$, i.e.,
$\ket{1;00}$ in the notation introduced in Eq. (7). Then, as will be discussed
below, the imaginary time evolution of the EAP state

$\displaystyle\ket{\beta}\propto e^{-\frac{1}{4}\beta\hat{H}}\ket{1;00}$ (11)

is locally indistinguishable from the Gibbs state $\hat{\rho}^{\mathrm{can}}$
at the inverse temperature $\beta$ with an exponentially small error, although
it is not an eigenstate of $\hat{H}$ 555Some readers may wonder, if the EAP
state $\ket{1;00}$ is an energy eigenstate, then it is invariant under the
imaginary time evolution, and hence $\ket{\beta}$ cannot describe a finite
temperature state. However, by using Theorem 1, it immediately follows that
$\ket{1;00}$ cannot be an energy eigenstate provided that $\hat{H}$ is
translation invariant (for single-site translations) and is a real matrix..

Let $\hat{O}$ be a local observable of interest. Since both $\ket{\beta}$ and
$\hat{\rho}^{\mathrm{can}}$ are translation invariant, without loss of
generality, we can assume that $\hat{O}$ has support around $j=N/4$. We then
introduce an approximation $\ket{\tilde{\beta}}$ for $\ket{\beta}$ as

$\displaystyle\ket{\tilde{\beta}}\propto
e^{-\frac{1}{4}\beta[\hat{H}-\hat{H}_{\mathrm{int}}]}\ket{1;00}.$ (12)

Here $\hat{H}_{\mathrm{int}}$ represents the interaction between the left half
$\\{1,2,\cdots,N/2\\}$ and the right half $\\{N/2+1,N/2+2,\cdots,N\\}$ of the
whole system, defined as a sum of interaction terms when $\hat{H}$ is
decomposed into a linear combination of the Pauli strings as in Eq. (5). The
imaginary time evolution is expected to be stable against local perturbations.
Indeed, the Gibbs state, which is proportional to the imaginary time evolution
operator, is not affected by perturbations on infinitely distant points
because systems under consideration are one dimensional and do not exhibit the
first-order phase transition [27]. Therefore, since $\hat{H}_{\mathrm{int}}$
is a sum of local observables defined around $j=1$ or $N/2$ at a distance of
$\mathcal{O}(N)$ from the support of $\hat{O}$, it is expected that
$\ket{\tilde{\beta}}$ approximate $\ket{\beta}$ around $j=N/4$ in the limit of
$N\to\infty$, i.e.,

$\displaystyle\lim_{N\to\infty}\braket{\tilde{\beta}}{\hat{O}}{\tilde{\beta}}=\lim_{N\to\infty}\braket{{\beta}}{\hat{O}}{{\beta}}.$
(13)

Under this assumption, we have the following [27]:

######  Theorem 4.

Suppose that $\hat{H}$ is short ranged and translation invariant and is a real
matrix in the basis formed by products of $\ket{0}_{j}$ and $\ket{1}_{j}$. Let
$\hat{O}$ be an arbitrary local observable. Then, if Eq. (13) is satisfied, it
holds that

$\displaystyle\lim_{N\to\infty}\braket{\beta}{\hat{O}}{\beta}=\lim_{N\to\infty}\mathrm{Tr}[\hat{\rho}^{\mathrm{can}}\hat{O}].$
(14)

In other words, $\ket{\beta}$ is a thermal pure state at finite temperature
while it is not an energy eigenstate. This is true for various systems where
the Hamiltonian is a real matrix in the spin basis.

We finally confirm the validity of Eq. (14) by numerically testing our
prediction on the transverse field Ising chain, defined by the Hamiltonian

$\displaystyle\hat{H}=-\sum_{j=1}^{N}\hat{\sigma}_{j}^{z}\hat{\sigma}_{j+1}^{z}-\sum_{j=1}^{N}\hat{\sigma}_{j}^{x}.$
(15)

We plot in Fig. 2 the $N$ dependence of the difference in the expectation
value of the transverse magnetization
$\hat{m}^{x}=\frac{1}{N}\sum_{j=1}^{N}\hat{\sigma}_{j}^{x}$ between the Gibbs
state in the thermodynamic limit and the imaginary time evolved EAP state
$\ket{\beta}$. We can confirm that the expectation value in $\ket{\beta}$
converges to the correct value. Furthermore, the convergence is exponentially
fast. Thus, by utilizing the EAP state evolved in imaginary time, one can
calculate the thermal equilibrium values of local observables.

Here we discuss the relation with the cTPQ state [35], which is also a type of
thermal pure states. The cTPQ state is the imaginary time evolution of a Haar
random state. The Haar random state is (almost) maximally entangled and is
locally indistinguishable from the maximally mixed state, as is the EAP state.
However, while constructing the TPQ state requires an imaginary time evolution
for $\beta/2$, constructing $\ket{\beta}$ requires only half of that,
$\beta/4$. In addition, unlike the cTPQ state, $\ket{\beta}$ does not entail
statistical uncertainties. This means that generating a single $\ket{\beta}$
suffices, without incurring any sampling costs.

Figure 2: $N$ dependence of the difference in the expectation value of the
transverse magnetization between the Gibbs state in the thermodynamic limit
and the imaginary time evolved EAP state for the transverse field Ising model.
We set the transverse field to $1$ and the inverse temperature to $\beta=1$.

Discussion.— We have studied a class of volume-law states, which we call
entangled antipodal pair (EAP) states, and thoroughly characterized these
states by providing the necessary and sufficient conditions that Hamiltonians
must satisfy in order to have an EAP state as an eigenstate. Moreover, we have
rigorously shown that some of such Hamiltonians are nonintegrable. Our EAP
states are indistinguishable from the Gibbs state at infinite temperature. In
other words, we have written down, for the first time, analytic expressions
for thermal energy eigenstates of nonintegrable many-body systems. We have
also devised a method for constructing finite-temperature thermal states using
an EAP state.

Some readers may be concerned that the energy eigenspace containing the EAP
state is excessively large because, with an exponentially large number of
states, it can be not so difficult to obtain a highly entangled state as their
linear combination, even when each is hardly entangled. However, numerical
results show that in Model 1, the degeneracy of the EAP state is only two in
the momentum sector 666Furthermore, by extending the locality of interactions
from two to three, we have also obtained an example of a pair of the EAP
eigenstate and the nonintegrable Hamiltonian that does not have degeneracy in
the corresponding momentum sector [27]. In addition, we can also confirm that
all states in the degenerate eigenspace obey a volume law with the maximal
coefficient $\log 2$. Hence, the nontriviality of our findings is not
diminished by the degeneracy.

Our result suggests that theoretical analysis may be feasible even for thermal
eigenstates of nonintegrable systems and may pave the way for giving a
provable example of the ETH. Since EAP states themselves are at infinite
temperature, constructing finite temperature thermal eigenstates remains
challenging. However, as we have shown, one can obtain thermal states at
arbitrary temperature by applying an imaginary time evolution, a low-
complexity operation that can be well approximated via a matrix product
operator with small bond dimension [37], to the EAP state. This implies that
the expressivity of states derived from EAP states is quite high. Thus, EAP
states would be one of the most promising starting points for constructing
thermal eigenstates at finite temperature. Furthermore, it will provide
theoretical methodologies not only for thermalization, but also for various
areas of quantum statistical mechanics that use thermal pure states, such as
the formulation of statistical mechanics and finite temperature simulations.

Acknowledgement.— We are grateful to N. Shiraishi, A. Shimizu, H. Katsura, K.
Fujii, K. Mizuta, H. Kono, M. Yamaguchi, and Z. Wei for useful discussions. We
would like to especially thank H. K. Park for careful reading and valuable
comments on a draft of this paper. YY is supported by the Special Postdoctoral
Researchers Program at RIKEN. YC is supperted by Japan Society for the
Promotion of Science KAKENHI Grant No. JP21J14313.

Note added: After the completion of this work, we became aware of a recent
related work by A. N. Ivanov and O. I. Motrunich (arXiv:2403.05515) which
constructs an exact volume-law eigenstate of the PXP and related models.

## References

  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Adv. Phys. 65, 239 (2016).
  * Gogolin and Eisert [2016] C. Gogolin and J. Eisert, Rep. Prog. Phys. 79, 056001 (2016).
  * Mori  _et al._ [2018] T. Mori, T. N. Ikeda, E. Kaminishi, and M. Ueda, J. Phys. B 51, 112001 (2018).
  * von Neumann [1929] J. von Neumann, Z. Phys. 57, 30 (1929).
  * Deutsch [1991] J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Phys. Rev. E 50, 888 (1994).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
  * Kim  _et al._ [2014] H. Kim, T. N. Ikeda, and D. A. Huse, Phys. Rev. E 90, 052105 (2014).
  * Beugeling  _et al._ [2014] W. Beugeling, R. Moessner, and M. Haque, Phys. Rev. E 89, 042112 (2014).
  * Steinigeweg  _et al._ [2014] R. Steinigeweg, A. Khodja, H. Niemeyer, C. Gogolin, and J. Gemmer, Phys. Rev. Lett. 112, 130403 (2014).
  * Verstraete  _et al._ [2008] F. Verstraete, V. Murg, and J. Cirac, Adv. Phys. 57, 143 (2008).
  * SchollwÃ¶ck [2011] U. SchollwÃ¶ck, Ann. Phys. 326, 96 (2011).
  * Moudgalya  _et al._ [2018a] S. Moudgalya, S. Rachel, B. A. Bernevig, and N. Regnault, Phys. Rev. B 98, 235155 (2018a).
  * Moudgalya  _et al._ [2018b] S. Moudgalya, N. Regnault, and B. A. Bernevig, Phys. Rev. B 98, 235156 (2018b).
  * Moudgalya  _et al._ [2020] S. Moudgalya, E. O’Brien, B. A. Bernevig, P. Fendley, and N. Regnault, Phys. Rev. B 102, 085120 (2020).
  * Shiraishi and Mori [2017] N. Shiraishi and T. Mori, Phys. Rev. Lett. 119, 030601 (2017).
  * Mori and Shiraishi [2017] T. Mori and N. Shiraishi, Phys. Rev. E 96, 022153 (2017).
  * 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. Vuletic, and M. D. Lukin, Nature (London) 551, 579 (2017).
  * Turner  _et al._ [2018a] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Nat. Phys. 14, 745 (2018a).
  * Turner  _et al._ [2018b] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Phys. Rev. B 94, 155134 (2018b).
  * Vitagliano  _et al._ [2010] G. Vitagliano, A. Riera, and J. I. Latorre, New J. Phys. 12, 113049 (2010).
  * RamÃrez  _et al._ [2014] G. RamÃrez, J. RodrÃguez-Laguna, and G. Sierra, J. Stat. Mech. 2014, P10004 (2014).
  * RamÃrez  _et al._ [2015] G. RamÃrez, J. RodrÃguez-Laguna, and G. Sierra, J. Stat. Mech. 2015, P06002 (2015).
  * Langlett  _et al._ [2022] C. M. Langlett, Z.-C. Yang, J. Wildeboer, A. V. Gorshkov, T. Iadecola, and S. Xu, Phys. Rev. B 105, L060301 (2022).
  * Bettaque and Swingle [2023] V. Bettaque and B. Swingle, NoRA: A Tensor Network Ansatz for Volume-Law Entangled Equilibrium States of Highly Connected Hamiltonians (2023).
  * Note [1] Some special EAP states appear in the study of integrable systems [38, 39, 40].
  * Note [2] See Supplemental Material for proofs of theorems, additional numerical results, and discussions.
  * Note [3] Note that, for some of states (7), $n$-site shift can cause a phase change, such as $\hat{\mathcal{T}}^{3}\mathinner{|{3;10,11,01}\rangle}=-\mathinner{|{3;10,11,01}\rangle}$, where $\hat{\mathcal{T}}$ is the one-site translation operator.
  * Note [4] There are only two solutions whose Hamiltonians are noninteracting [27].
  * Baxter [1982] R. J. Baxter,  _Exactly Solved Models in Statistical Mechanics_ (Academic Press, San Diego, 1982).
  * Korepin  _et al._ [1993] V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin,  _Quantum Inverse Scattering Method and Correlation Functions_ (Cambridge University Press, Cambridge, 1993).
  * Takahashi [1999] M. Takahashi,  _Thermodynamics of One-Dimensional Solvable Models_ (Cambridge University Press, Cambridge, 1999).
  * Tasaki [2024] H. Tasaki, Macroscopic Irreversibility in Quantum Systems: ETH and Equilibration in a Free Fermion Chain (2024).
  * Note [5] Some readers may wonder, if the EAP state $\mathinner{|{1;00}\rangle}$ is an energy eigenstate, then it is invariant under the imaginary time evolution, and hence $\mathinner{|{\beta}\rangle}$ cannot describe a finite temperature state. However, by using Theorem 1, it immediately follows that $\mathinner{|{1;00}\rangle}$ cannot be an energy eigenstate provided that $\hat{H}$ is translation invariant (for single-site translations) and is a real matrix.
  * Sugiura and Shimizu [2013] S. Sugiura and A. Shimizu, Phys. Rev. Lett. 111, 010401 (2013).
  * Note [6] Furthermore, by extending the locality of interactions from two to three, we have also obtained an example of a pair of the EAP eigenstate and the nonintegrable Hamiltonian that does not have degeneracy in the corresponding momentum sector [27].
  * Kuwahara  _et al._ [2021] T. Kuwahara, Á. M. Alhambra, and A. Anshu, Physical Review X 11 (2021).
  * Caetano and Komatsu [2022] J. Caetano and S. Komatsu, J. Statist. Phys. 187, 30 (2022).
  * Ekman [2022] C. Ekman, Crosscap states in the XXX spin-1/2 spin chain (2022).
  * Udupa  _et al._ [2023] A. Udupa, S. Sur, S. Nandy, A. Sen, and D. Sen, Phys. Rev. B 108, 214430 (2023).
  * Chiba [2024] Y. Chiba, Phys. Rev. B 109, 035123 (2024).
  * Shiraishi [2019] N. Shiraishi, Europhys. Lett. 128, 17002 (2019).
  * Park and Lee [2024] H. K. Park and S. Lee, Proof of the nonintegrability of PXP model and general spin-$1/2$ systems (2024).
  * Mehta [2004] M. L. Mehta,  _Random Matrices_ , 3rd ed. (Academic Press, New York, 2004).
  * Atas  _et al._ [2013] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Phys. Rev. Lett. 110, 084101 (2013).
  * Note [7] To remove the possibility of accidental discrete symmetries, we avoid eigenspace of special momentum such as $0$ and $\pi$.
  * Note [8] We can always take $J^{zy}=0$ by an appropriate rotation around $y$ axis.
  * Araki [1969] H. Araki, Commun. Math. Phys. 14, 120 (1969).
  * Araki [1975] H. Araki, Commun. Math. Phys. 44, 1 (1975).
  * Note [9] The results of exact diagonalization also show that, at least for $N=8,16,20$, degeneracy at $E=0$ in the whole Hilbert space remains only two. This means that all eigenstates of the Hamiltonian (S221) with eigenenergy $E=0$ are given by (linear combinations of) the EAP states $\mathinner{|{2;11,10}\rangle},\mathinner{|{2;10,11}\rangle}$.

Supplemental Material for
“Exact Thermal Eigenstates of Nonintegrable Spin Chains at Infinite
Temperature”

##  I Proof of Theorem 1

For the proof of Theorem 1, the following lemma is crucial:

######  Lemma 1.

Let $X$ and $Y$ be subsets of $\Lambda$ satisfying $D(X)<N/4$ and $D(Y)<N/4$.
For any EAP state $\ket{\mathrm{EAP}}$, we have

$\displaystyle\braket{\mathrm{EAP}}{\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}\bigotimes_{k\in
Y}\hat{\sigma}_{k}^{\nu_{k}}}{\mathrm{EAP}}=\begin{cases}1&\text{when $Y=X$
and $\nu_{k}=\mu_{k}$ for all $k\in Y$}\\\ \displaystyle\prod_{j\in
X}\omega_{j}^{\mu_{j}}&\text{when $Y=X+N/2$ and $\nu_{k}=\mu_{k-N/2}$ for all
$k\in Y$}\\\ 0&\text{otherwise}\end{cases}.$ (S1)

###### Proof.

We divide the lattice $\Lambda$ into four equally sized parts,
$\Lambda_{a}:=\\{aN/4+1,aN/4+2,...,(a+1)N/4\\}$ $(a=0,1,2,3)$. Since
$D(X)<N/4$, without loss of generality, we can take $X\subset\Lambda_{0}$.
Since $D(Y)<N/4$, $Y$ cannot have any intersection with both $\Lambda_{0}$ and
$\Lambda_{2}$, or with both $\Lambda_{1}$ and $\Lambda_{3}$.

Suppose that $Y\cap\Lambda_{1}\neq\emptyset$ and let $k^{*}\in
Y\cap\Lambda_{1}$. Because no Pauli operator acts on its antipodal site
$k^{*}+N/2\in\Lambda_{3}$ in the left hand side of Eq. (S1), we have

$\displaystyle\braket{\mathrm{EAP}}{\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}\bigotimes_{k\in
Y}\hat{\sigma}_{k}^{\nu_{k}}}{\mathrm{EAP}}=\braket{\mathrm{EAP}}{\left(\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}\right)\left(\bigotimes_{k\in
Y\setminus\\{k^{*}\\}}\hat{\sigma}_{k}^{\nu_{k}}\right)\hat{\sigma}_{k^{*}}^{\nu_{k^{*}}}}{\mathrm{EAP}}=0.$
(S2)

If $Y\cap\Lambda_{3}\neq\emptyset$, we can obtain the same result in almost
the same manner. Therefore, in the following, we only need to consider the two
cases $Y\subset\Lambda_{0}$ and $Y\subset\Lambda_{2}$.

Next we consider the case of $Y\subset\Lambda_{0}$. In the left hand side of
Eq. (S1), at most two Pauli operators act on each site $j\in\Lambda_{0}$, and
no Pauli operator acts on its antipodal site $j+N/2\in\Lambda_{2}$. Therefore,
unless $Y=X$ and $\nu_{k}=\mu_{k}$ for all $k\in Y$, the left hand side of Eq.
(S1) becomes zero. On the other hand, if $Y=X$ and $\nu_{k}=\mu_{k}$ for all
$k\in Y$, we obviously have

$\displaystyle\braket{\mathrm{EAP}}{\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}\bigotimes_{k\in
Y}\hat{\sigma}_{k}^{\nu_{k}}}{\mathrm{EAP}}=\braket{\mathrm{EAP}}{\mathrm{EAP}}=1.$
(S3)

Finally we consider the case of $Y\subset\Lambda_{2}$. In the left hand side
of Eq. (S1), at most one Pauli operator acts on each site $j\in\Lambda_{0}$,
and at most one Pauli operator acts on its antipodal site
$j+N/2\in\Lambda_{2}$. Using Eq. (4) of the main text, we have

$\displaystyle\braket{\mathrm{EAP}}{\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}\bigotimes_{k\in
Y}\hat{\sigma}_{k}^{\nu_{k}}}{\mathrm{EAP}}=\braket{\mathrm{EAP}}{\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}\bigotimes_{k\in
Y}\hat{\sigma}_{k-N/2}^{\nu_{k}}}{\mathrm{EAP}}\prod_{k\in
Y}\omega_{k}^{\nu_{k}}.$ (S4)

Applying the arguments of the previous paragraph, we can obtain the following:
Unless $Y-N/2=X$ and $\nu_{k}=\mu_{k-N/2}$ for all $k\in Y$, the left hand
side of Eq. (S1) becomes zero. On the other hand, if $Y-N/2=X$ and
$\nu_{k}=\mu_{k-N/2}$ for all $k\in Y$, we have

$\displaystyle\braket{\mathrm{EAP}}{\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}\bigotimes_{k\in
Y}\hat{\sigma}_{k}^{\nu_{k}}}{\mathrm{EAP}}=\braket{\mathrm{EAP}}{\mathrm{EAP}}\prod_{k\in
Y}\omega_{k}^{\nu_{k}}=\prod_{j\in X}\omega_{j}^{\mu_{j}}.$ (S5)

∎

###### Proof of Theorem 1.

In order to show that three statements are equivalent, we need to show that
(i) $\Rightarrow$ (ii), (ii) $\Rightarrow$ (iii), and (iii) $\Rightarrow$ (i).

(i) $\Rightarrow$ (ii): We assume that an EAP state $\ket{\mathrm{EAP}}$ is an
eigenstate of $\hat{H}$ with eigenvalue $\lambda$, i.e.,

$\displaystyle\hat{H}\ket{\mathrm{EAP}}=\lambda\ket{\mathrm{EAP}}.$ (S6)

Then, taking the inner product with $\ket{\mathrm{EAP}}$ and substituting Eq.
(5), we have

$\displaystyle\lambda=\braket{\mathrm{EAP}}{\hat{H}}{\mathrm{EAP}}=\sum_{X(\subset\Lambda)}\sum_{\vec{\mu}\in\\{x,y,z\\}^{X}}J_{X}^{\vec{\mu}}\braket{\mathrm{EAP}}{\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}}{\mathrm{EAP}}=0,$ (S7)

where the last equality follows from Lemma 1.

(ii) $\Rightarrow$ (iii): We assume that an EAP state $\ket{\mathrm{EAP}}$ is
an eigenstate of $\hat{H}$ with eigenvalue $0$, i.e.,

$\displaystyle\hat{H}\ket{\mathrm{EAP}}=0.$ (S8)

Then, substituting Eq. (5), we have

$\displaystyle\sum_{Y(\subset\Lambda)}\sum_{\vec{\nu}\in\\{x,y,z\\}^{Y}}J_{Y}^{\vec{\nu}}\bigotimes_{k\in
Y}\hat{\sigma}_{k}^{\nu_{k}}\ket{\mathrm{EAP}}=0.$ (S9)

Taking the inner product with $\displaystyle\bra{\mathrm{EAP}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}$, we get

$\displaystyle J_{X}^{\vec{\mu}}+J_{Y}^{\vec{\nu}}\prod_{j\in
X}\omega_{j}^{\mu_{j}}=0,$ (S10)

where $Y=X+N/2$ and $\nu_{k}=\mu_{k-N/2}$ for $k\in Y$. Since
$\omega_{j}^{\mu_{j}}=\pm 1$, this implies Eq. (6).

(iii) $\Rightarrow$ (i): We assume that Eq. (6) holds. Using Eqs. (4) and (5),
we have

$\displaystyle\hat{H}\ket{\mathrm{EAP}}$
$\displaystyle=\sum_{X(\subset\Lambda)}\sum_{\vec{\mu}\in\\{x,y,z\\}^{X}}J_{X}^{\vec{\mu}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}\ket{\mathrm{EAP}}$ (S11)
$\displaystyle=\sum_{X(\subset\Lambda)}\sum_{\vec{\mu}\in\\{x,y,z\\}^{X}}J_{X}^{\vec{\mu}}\prod_{j\in
X}\omega_{j}^{\mu_{j}}\bigotimes_{j\in
X}\hat{\sigma}_{j+N/2}^{\mu_{j}}\ket{\mathrm{EAP}}.$ (S12)

Then, substituting Eq. (6), we obtain

$\displaystyle\hat{H}\ket{\mathrm{EAP}}$
$\displaystyle=-\sum_{X(\subset\Lambda)}\sum_{\vec{\mu}\in\\{x,y,z\\}^{X}}J_{X+N/2}^{\vec{\mu}}\bigotimes_{j\in
X}\hat{\sigma}_{j+N/2}^{\mu_{j}}\ket{\mathrm{EAP}}=-\hat{H}\ket{\mathrm{EAP}},$
(S13)

which is equivalent to statement (ii). It obviously implies statement (i). ∎

Remark: The above proof also shows that, if we only need to show (iii)
$\Rightarrow$ (ii), we can relax the condition on the locality of interactions
in the main text, “$J_{X}^{\vec{\mu}}=0$ for any subset $X$ with $D(X)\geq
N/4$” to “$J_{X}^{\vec{\mu}}=0$ for any subset $X$ with $D(X)\geq N/2$.”

##  II List of translation-invariant and nearest-neighbor-interacting
Hamiltonians

We provide a list of translation-invariant (for single-site translations) and
nearest-neighbor-interacting Hamiltonians having an EAP state as an
eigenstate. First, we exclude the case of free spins, as it is trivial. Next,
for some cases where only one of $J^{\mu\nu}$ is non-zero, we find that
$2^{N/2}$ EAP states are degenerate, so we also exclude such cases. Then pairs
of the EAP state and the Hamiltonian are limited to the following five types
(and their equivalents obtained by appropriate permutations of directions of
the Pauli matrices). It can be readily confirmed through direct calculations
that the pairs of the EAP state and the Hamiltonian listed below satisfy Eq.
(6).

###  II.1 Case where the EAP state is invariant under $1$-site translation

The EAP state $\ket{1;00}$ is an eigenstate of the Hamiltonian defined by

$\displaystyle\hat{H}=\sum_{j=1}^{N}\left(J^{xy}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{y}+J^{yx}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{x}+J^{yz}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{z}+J^{zy}\hat{\sigma}_{j}^{z}\hat{\sigma}_{j+1}^{y}+h^{y}\hat{\sigma}_{j}^{y}\right),$
(S14)

for arbitrary values of $J^{xy},J^{yx},J^{yz},J^{zy},h^{y}$. This model is
Model 1 in Theorem 2.

###  II.2 Case where the EAP state is invariant under $2$-site translation

Suppose that $N/2$ is a multiple of two.

The EAP state $\ket{2;10,01}$ is an eigenstate of the Hamiltonian defined by

$\displaystyle\hat{H}=\sum_{j=1}^{N}\left(J^{xx}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{x}+J^{zz}\hat{\sigma}_{j}^{z}\hat{\sigma}_{j+1}^{z}\right),$
(S15)

for arbitrary values of $J^{xx},J^{zz}$. This model can be mapped onto free
fermions via the Jordan-Wigner transformation.

In addition, the EAP state $\ket{2;10,11}$ is an eigenstate of the Hamiltonian
defined by

$\displaystyle\hat{H}=\sum_{j=1}^{N}\left(J^{xx}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{x}+J^{yy}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{y}+J^{xy}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{y}+J^{yx}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{x}+h^{z}\hat{\sigma}_{j}^{z}\right),$
(S16)

for arbitrary values of $J^{xx},J^{yy},J^{xy},J^{yx},h^{z}$. This model can
also be mapped onto free fermions via the Jordan-Wigner transformation.

###  II.3 Case where the EAP state is invariant under $3$-site translation

Suppose that $N/2$ is a multiple of three.

The EAP state $\ket{3;10,11,01}$ is an eigenstate of the Hamiltonian defined
by

$\displaystyle\hat{H}=\sum_{j=1}^{N}\left(J^{xy}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{y}+J^{yz}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{z}\right),$
(S17)

for arbitrary values of $J^{xy},J^{yz}$. This model is Model 2 in Theorem 2.
As shown in Theorem 3, this model is nonintegrable.

###  II.4 Case where the EAP state is invariant under $4$-site translation

Suppose that $N/2$ is a multiple of four.

The EAP state $\ket{4;10,11,00,01}$ is an eigenstate of the Hamiltonian
defined by

$\displaystyle\hat{H}=\sum_{j=1}^{N}\left(J^{xx}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{x}+J^{yz}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{z}\right)$
(S18)

for arbitrary values of $J^{xx},J^{yz}$. This model is Model 3 in Theorem 2.
As shown in Theorem 3, this model is nonintegrable.

##  III Proof of Theorem 3

First we define a $k$-local conserved quantity (which is the same as one given
in Ref. [41]) by the operator $\hat{Q}$ that commutes with the Hamiltonian

$\displaystyle[\hat{Q},\hat{H}]=0$ (S19)

and can be written as

$\displaystyle\hat{Q}=\sum_{\ell=1}^{k}\sum_{\bm{A}^{\ell}}q_{\bm{A}^{\ell}_{j}}^{(\ell)}\hat{\bm{A}}^{\ell}_{j}+q_{I}^{(0)}\hat{I}.$
(S20)

Here $\bm{A}^{\ell}$ represents a sequence of symbols,
$A^{1},A^{2},...,A^{\ell}$ satisfying

$\displaystyle A^{1},A^{\ell}\in\\{X,Y,Z\\}$ (S21) $\displaystyle
A^{2},...,A^{\ell-1}\in\\{X,Y,Z,I\\}$ (S22)

and $\hat{\bm{A}}^{\ell}_{j}$ represents the product of the corresponding
Pauli operators on the sites $\\{j,j+1,...,j+\ell-1\\}$:

$\displaystyle\hat{\bm{A}}^{\ell}_{j}=\hat{A}^{1}_{j}\hat{A}^{2}_{j+1}...\hat{A}^{\ell}_{j+\ell-1}.$
(S23)

In Eq. (S20), $q_{\bm{A}^{\ell}_{j}}^{(\ell)}\in\mathbb{R}$ are the expansion
coefficients. [We add the superscript $(\ell)$ in order to emphasize its
value.] The crucial point of Eq. (S20) is that $\hat{Q}$ does not include
$\hat{\bm{A}}^{\ell}_{j}$ with $\ell>k$.

Now we give the precise expression of Theorem 3 of the main text, which is
represented by the following two theorems:

######  Theorem 3.A.

In Model 2 with $J^{xy},J^{yz}\neq 0$, and for $k\leq N/2$, there is no
$k$-local conserved quantity that is linearly independent of the Hamiltonian
and the identity.

######  Theorem 3.B.

In Model 3 with $J^{xx},J^{yz}\neq 0$, and for $k\leq N/2$, there is no
$k$-local conserved quantity that is linearly independent of the Hamiltonian
and the identity.

In the remaining of this section, we prove these theorems by adapting the
theoretical approach to prove the absence of local conserved quantities, which
was introduced by N. Shiraishi [42]. There are only a few examples of such
proofs [42, 41, 43]. This approach starts from solving Eq. (S19) with respect
to the coefficients with largest locality, $q_{\bm{A}^{k}_{j}}^{(k)}$, and
showing that $q_{\bm{A}^{k}_{j}}^{(k)}=0$. When solving Eq. (S19), we need to
calculate many commutators such as

$\displaystyle[q_{X_{j}IX}^{(3)}\hat{X}_{j}\hat{I}_{j+1}\hat{X}_{j+2},\sum_{j^{\prime}}J^{yz}\hat{Y}_{j^{\prime}}\hat{Z}_{j^{\prime}+1}]$
$\displaystyle=q_{X_{j}IX}^{(3)}J^{yz}[\hat{X}_{j}\hat{I}_{j+1}\hat{X}_{j+2},\hat{Y}_{j-1}\hat{Z}_{j}+\hat{Y}_{j}\hat{Z}_{j+1}+\hat{Y}_{j+1}\hat{Z}_{j+2}+\hat{Y}_{j+2}\hat{Z}_{j+3}]$
(S24)
$\displaystyle=2iq_{X_{j}IX}^{(3)}J^{yz}\bigl{(}-\hat{Y}_{j-1}\hat{Y}_{j}\hat{I}_{j+1}\hat{X}_{j+2}+\hat{Z}_{j}\hat{Z}_{j+1}\hat{X}_{j+2}-\hat{X}_{j}\hat{Y}_{j+1}\hat{Y}_{j+2}+\hat{X}_{j}\hat{I}_{j+1}\hat{Z}_{j+2}\hat{Z}_{j+3}\bigr{)}.$
(S25)

For simplicity of notation, we write $q_{X_{j}IX}^{(3)}$ in place of
$q_{X_{j}I_{j+1}X_{j+2}}^{(3)}$. In order to express such calculations
efficiently, we use the following diagrammatic notation:

$\displaystyle\begin{array}[]{rlllllll}&&X&I&X\\\ &Y&Z&&\\\
\hline\cr-&Y_{j-1}&Y&I&X\end{array}\quad\begin{array}[]{rlllllll}&X&I&X\\\
&Y&Z&\\\ \hline\cr&Z_{j}&Z&X\end{array}\quad\begin{array}[]{rlllllll}&X&I&X\\\
&&Y&Z\\\
\hline\cr-&X_{j}&Y&Y\end{array}\quad\begin{array}[]{rlllllll}&X&I&X&\\\
&&&Y&Z\\\ \hline\cr&X_{j}&I&Z&Z\end{array}.$ (S38)

These four diagrams correspond to the four terms in Eq. (S25). In each
diagram, the first row represents the term from $\hat{Q}$, the second row the
term from $\hat{H}$, and the third row the result of the commutator. For
simplicity of notation, we add the site index only for the leftmost operators
in the third row. In addition, we call the first row of the diagram
“$\ell$-local input”, and the third row of the diagram “$\ell$-local output”,
when they consist of $X$, $Y$, $Z$, and $I$ on $\ell$ consecutive sites.

###  III.1 Proof of Theorem 3.A

This subsection proves Theorem 3.A. Throughout this subsection, we consider
Model 2 and assume $J^{xy},J^{yz}\neq 0$ and $k\leq N/2$. (The reason for the
assumption $k\leq N/2$ is the same as one discussed in Sec. VI A of Ref.
[41].)

###### Proof.

The proof of Theorem 3.A is divided into three parts. The first part
investigates the coefficients with largest locality,
$q_{\bm{A}^{k}_{j}}^{(k)}$. For the coefficients of the form
$q_{Z_{j}A^{2}...A^{k-1}X}^{(k)}$, we have

$\displaystyle\begin{array}[]{rlllllll}&Z&A^{2}&...&A^{k-1}&X&\\\
&&&...&&Y&Z\\\ \hline\cr&Z_{j}&A^{2}&...&A^{k-1}&Z&Z\end{array}.$ (S42)

Because a $(k+1)$-local output can be obtained only when the Hamiltonian term
is applied to the edges of $k$-local inputs, there are at most two $k$-local
inputs that contribute to one $(k+1)$-local output. However, since the left
end of the output of Eq. (S42) is $Z$, the other contribution does not exist.
Furthermore, from Eq. (S19), the sum of all contribution to the output
$Z_{j}A^{2}...A^{k-1}ZZ$ must vanish, and therefore we have

$\displaystyle J^{yz}q_{Z_{j}A^{2}...A^{k-1}X}^{(k)}=0.$ (S43)

In a similar manner, we can obtain the following lemma:

######  Lemma 2.

For $2\leq k\leq N/2$, the solution of Eq. (S19) satisfies

$\displaystyle q_{Z_{j}A^{2}...A^{k-1}A^{k}}^{(k)}$ $\displaystyle=0$ (S44)
$\displaystyle q_{A^{1}_{j}A^{2}...A^{k-1}X}^{(k)}$ $\displaystyle=0$ (S45)

for all $j\in\Lambda$. Here the symbols $A^{1},A^{2},...,A^{k-1},A^{k}$ that
are not specified can be any symbols satisfying Eqs. (S21) and (S22).

Next we examine the coefficients of the form
$q_{X_{j}A^{2}...A^{k-1}Y}^{(k)}$. They have the following contributions

$\displaystyle\begin{array}[]{rlllllll}&X&A^{2}&A^{3}&...&A^{k-1}&Y&\\\
&&&&...&&X&Y\\\ \hline\cr-&X_{j}&A^{2}&A^{3}&...&A^{k-1}&Z&Y\end{array}.$
(S49)

When $A^{2}=I,Y$, this $(k+1)$-local output does not have the other
contribution. When $A^{2}=X$, it has the other contribution,

$\displaystyle\begin{array}[]{rlllllll}&&Z&A^{3}&...&A^{k-1}&Z&Y\\\
&X&Y&&...&&&\\\ \hline\cr-&X_{j}&X&A^{3}&...&A^{k-1}&Z&Y\end{array},$ (S53)

which however vanishes from Lemma 2. In a similar manner, we can obtain the
following lemma:

######  Lemma 3.

For $3\leq k\leq N/2$, the solution of Eq. (S19) satisfies

$\displaystyle q_{X_{j}A^{2}...A^{k-1}A^{k}}^{(k)}$
$\displaystyle=0\quad\text{for }A^{2}=I,Y,X$ (S54) $\displaystyle
q_{Y_{j}A^{2}...A^{k-1}A^{k}}^{(k)}$ $\displaystyle=0\quad\text{for
}A^{2}=I,Z$ (S55) $\displaystyle q_{A^{1}_{j}A^{2}...A^{k-1}Z}^{(k)}$
$\displaystyle=0\quad\text{for }A^{k-1}=I,Y,Z$ (S56) $\displaystyle
q_{A^{1}_{j}A^{2}...A^{k-1}Y}^{(k)}$ $\displaystyle=0\quad\text{for
}A^{k-1}=I,X$ (S57)

for all $j\in\Lambda$. Here the symbols $A^{1},A^{2},...,A^{k-1},A^{k}$ that
are not specified can be any symbols satisfying Eqs. (S21) and (S22).

Furthermore, we can obtain relation between two of the remaining $k$-local
inputs as in

$\displaystyle\begin{array}[]{rlllllll}&X&Z&A^{3}&...&A^{k-1}&Y&\\\
&&&&...&&X&Y\\\
\hline\cr-&X_{j}&Z&A^{3}&...&A^{k-1}&Z&Y\end{array}\quad\begin{array}[]{rlllllll}&&X&A^{3}&...&A^{k-1}&Z&Y\\\
&X&Y&&...&&&\\\ \hline\cr&X_{j}&Z&A^{3}&...&A^{k-1}&Z&Y\end{array},$ (S64)

which results in

$\displaystyle-J^{xy}q_{X_{j}ZA^{3}...A^{k-1}Y}^{(k)}+J^{xy}q_{X_{j+1}A^{3}...A^{k-1}ZY}^{(k)}=0.$
(S65)

If $A^{3}$ is not $Z$, we have $q_{X_{j+1}A^{3}...A^{k-1}ZY}^{(k)}=0$ from
Lemma 3, and hence $q_{X_{j}ZA^{3}...A^{k-1}Y}^{(k)}=0$. By using such a
relation, we can shift the symbols $A^{3},...,A^{k-1}$ to the left and we can
determine these symbols. As a result, we can obtain the following proposition:

######  Proposition 1.

For any $j\in\Lambda$ and for any $\bm{A}^{k}$, the solution of Eq. (S19)
satisfies

$\displaystyle q_{\bm{A}^{k}_{j}}^{(k)}=0,$ (S66)

except for

$\displaystyle q_{X_{j}(Z)^{k-2}Y}^{(k)},\quad q_{Y_{j}(X)^{k-2}Z}^{(k)},\quad
q_{Y_{j}(X)^{n}Y(Z)^{k-n-3}Y}^{(k)}\quad\text{with }n=0,1,...,k-3.$ (S67)

In addition, these remaining coefficients are independent of the site $j$ and
satisfy

$\displaystyle\frac{q_{Y_{j}(X)^{n}Y(Z)^{k-n-3}Y}^{(k)}}{(J^{xy})^{k-n-2}(J^{yz})^{n+1}}=-\frac{q_{Y_{j}(X)^{k-2}Z}^{(k)}}{(J^{yz})^{k-1}}=-\frac{q_{X_{j}(Z)^{k-2}Y}^{(k)}}{(J^{xy})^{k-1}}=-\frac{q_{X_{1}(Z)^{k-2}Y}^{(k)}}{(J^{xy})^{k-1}}\quad\text{for
}n=0,1,...,k-3.$ (S68)

Here we used a shorthand notation of a sequence of symbols

$\displaystyle(A)^{m}:=\underbrace{AA...A}_{m\,\mathrm{times}}\quad\text{ for
}A=X,Y,Z,I.$ (S69)

Therefore we only need to show that one of these remaining coefficients is
zero.

As the second part of the proof, we examine the coefficient
$q_{X_{j}(Z)^{k-2}Y}^{(k)}(=q_{X_{1}(Z)^{k-2}Y}^{(k)})$. We consider the
contribution from $q_{X_{j}(Z)^{k-2}Y}^{(k)}$ to a $k$-local output which can
also include the contribution from $(k-1)$-local inputs. For instance,

$\displaystyle\begin{array}[]{rlllllll}&X&(Z)^{k-4}&Z&Z&Y\\\ &&&Y&Z&\\\
\hline\cr-&X_{j}&(Z)^{k-4}&X&I&Y\end{array}\quad\begin{array}[]{rlllllll}&&X&(Z)^{k-5}&X&I&Y\\\
&X&Y&&&&\\\ \hline\cr&X_{j}&Z&(Z)^{k-5}&X&I&Y\end{array}$ (S76)

are the only contribution to the $k$-local output $X_{j}(Z)^{k-4}XIY$, and
hence we have

$\displaystyle
J^{xy}q_{X_{j+1}(Z)^{k-5}XIY}^{(k-1)}=J^{yz}q_{X_{1}(Z)^{k-2}Y}^{(k)}.$ (S77)

For coefficients of $(k-1)$-local inputs, we obtain the following: All
contributions to $k$-local output $X_{j}(Z)^{k-n-5}XI(Z)^{n+1}Y$ (for
$n=0,...,k-6$) are given by

$\displaystyle\begin{array}[]{rllllllllll}&X&(Z)^{k-n-5}&X&I&(Z)^{n}&Y&\\\
&&&&&&X&Y\\\
\hline\cr-&X_{j}&(Z)^{k-n-5}&X&I&(Z)^{n}&Z&Y\end{array}\quad\begin{array}[]{rllllllllll}&&X&(Z)^{k-n-6}&X&I&(Z)^{n+1}&Y\\\
&X&Y&&&&&\\\ \hline\cr&X_{j}&Z&(Z)^{k-n-5}&X&I&(Z)^{n+1}&Y\end{array}$ (S84)
$\displaystyle\begin{array}[]{rllllllllll}&X&(Z)^{k-n-5}&Z&Z&(Z)^{n+1}&Y\\\
&&&Y&Z&&\\\
\hline\cr-&X_{j}&(Z)^{k-n-5}&X&I&(Z)^{n+1}&Y\end{array}\quad\begin{array}[]{rllllllllll}&X&(Z)^{k-n-5}&X&X&X&(Z)^{n}&Y\\\
&&&&X&Y&&\\\ \hline\cr&X_{j}&(Z)^{k-n-5}&X&I&Z&(Z)^{n}&Y\end{array},$ (S91)

which result in

$\displaystyle
J^{xy}\bigl{(}q_{X_{j+1}(Z)^{k-n-6}XI(Z)^{n+1}Y}^{(k-1)}-q_{X_{j}(Z)^{k-n-5}XI(Z)^{n}Y}^{(k-1)}\bigr{)}=J^{yz}q_{X_{1}(Z)^{k-2}Y}^{(k)}\quad\text{for
all }n=0,...,k-6.$ (S92)

For the coefficient $q_{X_{j}XI(Z)^{k-5}Y}^{(k-1)}$, which appears in $n=k-6$
case of the above equation, we can obtain another relation

$\displaystyle
J^{xy}\bigl{(}-q_{Z_{j+1}I(Z)^{k-4}Y}^{(k-1)}-q_{X_{j}XI(Z)^{k-5}Y}^{(k-1)}\bigr{)}=J^{yz}q_{X_{1}(Z)^{k-2}Y}^{(k)},$
(S93)

in a similar manner. We can also obtain the following relation for the
coefficient $q_{Z_{j}I(Z)^{k-4}Y}^{(k-1)}$,

$\displaystyle
J^{xy}q_{Z_{j}I(Z)^{k-4}Y}^{(k-1)}=2J^{yz}q_{X_{1}(Z)^{k-2}Y}^{(k)},$ (S94)

by considering the contributions to $Z_{j}I(Z)^{k-3}Y$,

$\displaystyle\begin{array}[]{rllllllllll}&Z&I&(Z)^{k-4}&Y&\\\ &&&&X&Y\\\
\hline\cr-&Z_{j}&I&(Z)^{k-4}&Z&Y\end{array}\quad\begin{array}[]{rllllllllll}&Z&X&X&(Z)^{k-4}&Y\\\
&&X&Y&&\\\
\hline\cr&Z_{j}&I&Z&(Z)^{k-4}&Y\end{array}\quad\begin{array}[]{rllllllllll}&Y&Y&(Z)^{k-3}&Y\\\
&X&Y&&\\\
\hline\cr-&Z_{j}&I&(Z)^{k-3}&Y\end{array}\quad\begin{array}[]{rllllllllll}&X&Z&(Z)^{k-3}&Y\\\
&Y&Z&&\\\ \hline\cr&Z_{j}&I&(Z)^{k-3}&Y\end{array}.$ (S107)

Because the sum of the left-hand sides of Eqs. (S77), (S92)–(S94) (by choosing
the site $j$ appropriately) become zero, we have

$\displaystyle 0$
$\displaystyle=J^{xy}q_{X_{1}(Z)^{k-5}XIY}^{(k-1)}+\sum_{n=0}^{k-6}J^{xy}\bigl{(}q_{X_{n+2}(Z)^{k-n-6}XI(Z)^{n+1}Y}^{(k-1)}-q_{X_{n+1}(Z)^{k-n-5}XI(Z)^{n}Y}^{(k-1)}\bigr{)}$
$\displaystyle\hskip
10.0pt+J^{xy}\bigl{(}-q_{Z_{k-3}I(Z)^{k-4}Y}^{(k-1)}-q_{X_{k-4}XI(Z)^{k-5}Y}^{(k-1)}\bigr{)}+J^{xy}q_{Z_{k-3}I(Z)^{k-4}Y}^{(k-1)}$
(S108) $\displaystyle=(k-1)J^{yz}q_{X_{1}(Z)^{k-2}Y}^{(k)}.$ (S109)

Thus we obtain the following proposition:

######  Proposition 2.

For $3\leq k\leq N/2$, the solution of Eq. (S19) satisfies

$\displaystyle q_{X_{1}(Z)^{k-2}Y}^{(k)}=0.$ (S110)

By combining Propositions 1 and 2, we have

$\displaystyle q_{\bm{A}^{k}_{j}}^{(k)}=0\quad\text{for all }j\text{ and
}\bm{A}^{k}.$ (S111)

This means that $\hat{Q}$ is a $(k-1)$-local conserved quantity. Applying the
same argument to $k-1$, $k-2$,…, and $3$-local conserved quantity, we can show
that any $k$-local conserved quantity with $k\leq N$ have to be a $2$-local
conserved quantity.

As the third part of the proof, we analyze the coefficients
$q_{\bm{A}^{\ell}_{j}}^{(\ell)}$ with $\ell\leq k$ in the case of $k\leq 2$.
From Lemma 2, we only need to consider the coefficients of the form
$q_{X_{j}Y}^{(2)},q_{X_{j}Z}^{(2)},q_{Y_{j}Y}^{(2)},q_{Y_{j}Z}^{(2)}$ and
$q_{A_{j}}^{(1)}$. Furthermore, the coefficient $q_{X_{j}Z}^{(2)}$ vanishes
because

$\displaystyle\begin{array}[]{rlllllll}&X&Z&\\\ &&X&Y\\\
\hline\cr&X_{j}&Y&Y\end{array}$ (S115)

is the only contribution to the $3$-local output $X_{j}YY$. The coefficient
$q_{Y_{j}Y}^{(2)}$ vanishes by a similar reason. In addition, because
$2$-local output comes only from $1$-local input, we can easily show that
$q_{X_{j}}^{(1)}=q_{Y_{j}}^{(1)}=q_{Z_{j}}^{(1)}=0$. For the remaining
coefficients $q_{X_{j}Y}^{(2)}$ and $q_{Y_{j}Z}^{(2)}$, we can easily show
that they are independent of the site $j$ and are related to each other by

$\displaystyle J^{xy}q_{Y_{j}Z}^{(2)}-J^{yz}q_{X_{j+1}Y}^{(2)}=0,$ (S116)

which results in the following proposition:

######  Proposition 3.

Any $2$-local conserved quantity $\hat{Q}$ can be written as

$\displaystyle\hat{Q}=a\hat{H}+b\hat{I},$ (S117)

with arbitrary constants $a,b\in\mathbb{R}$.

From Eq. (S111) and Proposition 3, we obtain Theorem 3.A. ∎

###  III.2 Proof of Theorem 3.B

This subsection proves Theorem 3.B. Throughout this subsection, we consider
Model 3 and assume $J^{xx},J^{yz}\neq 0$ and $k\leq N/2$. (The reason for the
assumption $k\leq N/2$ is the same as one discussed in Sec. VI A of Ref.
[41].)

###### Proof.

The proof of Theorem 3.A is divided into three parts. The first part
investigates the coefficients with largest locality,
$q_{\bm{A}^{k}_{j}}^{(k)}$. In a manner similar to the proof of Lemma 2, we
can show the following lemma:

######  Lemma 4.

For $2\leq k\leq N/2$, the solution of Eq. (S19) satisfies

$\displaystyle q_{Z_{j}A^{2}...A^{k-1}A^{k}}^{(k)}$ $\displaystyle=0$ (S118)
$\displaystyle q_{A^{1}_{j}A^{2}...A^{k-1}Y}^{(k)}$ $\displaystyle=0$ (S119)

for all $j\in\Lambda$. Here the symbols $A^{1},A^{2},...,A^{k-1},A^{k}$ that
are not specified can be any symbols satisfying Eqs. (S21) and (S22).

Furthermore, we can show the following lemma in a manner similar to the proof
of Lemma 3:

######  Lemma 5.

For $3\leq k\leq N/2$, the solution of Eq. (S19) satisfies

$\displaystyle q_{X_{j}A^{2}...A^{k-1}A^{k}}^{(k)}$
$\displaystyle=0\quad\text{for }A^{2}=I,Y,X$ (S120) $\displaystyle
q_{Y_{j}A^{2}...A^{k-1}A^{k}}^{(k)}$ $\displaystyle=0\quad\text{for
}A^{2}=I,Z$ (S121) $\displaystyle q_{A^{1}_{j}A^{2}...A^{k-1}X}^{(k)}$
$\displaystyle=0\quad\text{for }A^{k-1}=I,X,Z$ (S122) $\displaystyle
q_{A^{1}_{j}A^{2}...A^{k-1}Z}^{(k)}$ $\displaystyle=0\quad\text{for
}A^{k-1}=I,Y$ (S123)

for all $j\in\Lambda$. Here the symbols $A^{1},A^{2},...,A^{k-1},A^{k}$ that
are not specified can be any symbols satisfying Eqs. (S21) and (S22).

By shifting the symbols $A^{3},...,A^{k-1}$ to the left as we did to obtain
Proposition 1, many coefficients can be shown to be zero. To explain the
result, we introduce a version of “doubling product.” It was originally
introduced by N. Shiraishi [42]. Our version is modified for analyzing Model 3
as follows: We call a sequence of the Pauli operators doubling product, if it
can by written as $\overline{A^{1}_{j}A^{2}...A^{n}}$, where

$\displaystyle\overline{A^{1}_{j}}$
$\displaystyle:=\begin{cases}X_{j}X_{j+1}\quad&\text{when }A^{1}=X\\\
Y_{j}Z_{j+1}\quad&\text{when }A^{1}=Y\end{cases}$ (S124)
$\displaystyle\overline{A^{1}_{j}A^{2}...A^{n}A^{n+1}}$
$\displaystyle:=\begin{cases}c\overline{A^{1}_{j}A^{2}...A^{n}}\times
X_{j+n}X_{j+n+1}\quad&\text{when }A^{n+1}=X,\ A^{n}\neq X\\\
c\overline{A^{1}_{j}A^{2}...A^{n}}\times Y_{j+n}Z_{j+n+1}\quad&\text{when
}A^{n+1}=Y\end{cases}.$ (S125)

Here $c$ is chosen from $\\{\pm i\\}$ to make its coefficient $1$. In
addition, we introduce $J_{\overline{A^{1}...A^{n}}}$ by

$\displaystyle J_{\overline{A^{1}}}$
$\displaystyle:=\begin{cases}J^{xx}\quad\text{when }A^{1}=X,\\\
J^{yz}\quad\text{when }A^{1}=Y\end{cases}$ (S126) $\displaystyle
J_{\overline{A^{1}...A^{n}A^{n+1}}}$
$\displaystyle:=\begin{cases}J_{\overline{A^{1}...A^{n}}}\times
J^{xx}\quad&\text{when }A^{n}=Y,\ A^{n+1}=X\\\
J_{\overline{A^{1}...A^{n}}}\times J^{yz}\quad&\text{when }A^{n}=X,\
A^{n+1}=Y\\\ -J_{\overline{A^{1}...A^{n}}}\times J^{yz}\quad&\text{when
}A^{n}=Y,\ A^{n+1}=Y\end{cases}.$ (S127)

Then we can obtain the following proposition:

######  Proposition 4.

For any $j\in\Lambda$ and for any $\bm{A}^{k}$ other than doubling product,
the solution of Eq. (S19) satisfies

$\displaystyle q_{\bm{A}^{k}_{j}}^{(k)}=0.$ (S128)

For the case where $\bm{A}^{k}$ is given by a doubling product Eqs. (S124) and
(S125), these remaining coefficients are independent of the site $j$ and are
related to each other by

$\displaystyle\frac{q_{\overline{A^{1}_{j}...A^{k-1}}}^{(k)}}{J_{\overline{A^{1}...A^{k-1}}}}=\frac{q_{Y_{j}(X)^{k-2}Z}^{(k)}}{J^{yz}(-J^{yz})^{k-2}}=\frac{q_{Y_{1}(X)^{k-2}Z}^{(k)}}{J^{yz}(-J^{yz})^{k-2}}$
(S129)

for any $j\in\Lambda$.

Therefore we only need to show that one of these remaining coefficients is
zero.

As the second part of the proof, we examine the coefficient
$q_{Y_{1}(X)^{k-2}Z}^{(k)}$. We consider the contribution from the $k$-local
input $Y_{j}(X)^{k-2}Z$ to a $k$-local output which can also include the
contribution from $(k-1)$-local inputs. For instance,

$\displaystyle\begin{array}[]{rlllllll}&Y&X&(X)^{k-3}&Z\\\ &X&X&&\\\
\hline\cr-&Z_{j}&I&(X)^{k-3}&Z\end{array}\quad\begin{array}[]{rlllllll}&X&Z&(X)^{k-3}&Z\\\
&Y&Z&&\\\
\hline\cr&Z_{j}&I&(X)^{k-3}&Z\end{array}\quad\begin{array}[]{rlllllll}&Z&Y&Y&(X)^{k-4}&Z\\\
&&Y&Z&&\\\
\hline\cr&Z_{j}&I&X&(X)^{k-4}&Z\end{array}\quad\begin{array}[]{rlllllll}&Z&I&(X)^{k-4}&Z&\\\
&&&&Y&Z\\\ \hline\cr-&Z_{j}&I&(X)^{k-4}&X&Z\end{array}$ (S142)

are the only contribution to the $k$-local output $Z_{j}I(X)^{k-3}Z$. Note
that the contribution from the third diagram vanishes because of Lemma 4, and
the contribution from the second diagram satisfies

$\displaystyle
q_{X_{j}Z(X)^{k-3}Z}^{(k)}=-\frac{J^{xx}}{J^{yz}}q_{Y_{1}(X)^{k-2}Z}^{(k)}$
(S143)

because $XZ(X)^{k-3}Z$ can be written as a doubling product
$\overline{X(Y)^{k-2}}$ and Eq. (S129) in Proposition 4 is applicable. Hence
we have

$\displaystyle
J^{yz}q_{Z_{j}I(X)^{k-4}Z}^{(k-1)}=-2J^{xx}q_{Y_{1}(X)^{k-2}Z}^{(k)}.$ (S144)

In a similar manner, we can obtain

$\displaystyle
J^{yz}q_{X_{j}YI(X)^{k-5}Z}^{(k-1)}-J^{xx}q_{Z_{j+1}I(X)^{k-4}Z}^{(k-1)}$
$\displaystyle=-\frac{(J^{xx})^{2}}{J^{yz}}q_{Y_{1}(X)^{k-2}Z}^{(k)}$ (S145)
$\displaystyle-J^{yz}q_{Y_{j}YYI(X)^{k-6}Z}^{(k-1)}-J^{yz}q_{X_{j+1}YI(X)^{k-5}Z}^{(k-1)}$
$\displaystyle=-\frac{(J^{xx})^{2}}{J^{yz}}q_{Y_{1}(X)^{k-2}Z}^{(k)}$ (S146)
$\displaystyle-J^{yz}q_{Y_{j}(X)^{n+1}YYI(X)^{k-7-n}Z}^{(k-1)}+J^{yz}q_{Y_{j+1}(X)^{n}YYI(X)^{k-6-n}Z}^{(k-1)}$
$\displaystyle=-\frac{(J^{xx})^{2}}{J^{yz}}q_{Y_{1}(X)^{k-2}Z}^{(k)}\quad\text{for
}n=0,...,k-7$ (S147) $\displaystyle J^{yz}q_{Y_{j}(X)^{k-6}YYIZ}^{(k-1)}$
$\displaystyle=-\frac{(J^{xx})^{2}}{J^{yz}}q_{Y_{1}(X)^{k-2}Z}^{(k)}$ (S148)

from the diagrams

$\displaystyle\begin{array}[]{rlllllll}&&Z&I&(X)^{k-4}&Z\\\ &X&X&&&\\\
\hline\cr-&X_{j}&Y&I&(X)^{k-4}&Z\end{array}\quad\begin{array}[]{rlllllll}&X&Y&I&(X)^{k-5}&Z&\\\
&&&&&Y&Z\\\
\hline\cr-&X_{j}&Y&I&(X)^{k-5}&X&Z\end{array}\quad\begin{array}[]{rlllllll}&X&Z&X&(X)^{k-4}&Z\\\
&&X&X&&\\\
\hline\cr&X_{j}&Y&I&(X)^{k-4}&Z\end{array}\quad\begin{array}[]{rlllllll}&X&Y&Y&Y&(X)^{k-5}&Z&\\\
&&&Y&Z&&\\\ \hline\cr-&X_{j}&Y&I&X&(X)^{k-5}&Z\end{array},$ (S161)
$\displaystyle\begin{array}[]{rlllllll}&&X&Y&I&(X)^{k-5}&Z\\\ &Y&Z&&&\\\
\hline\cr-&Y_{j}&Y&Y&I&(X)^{k-5}&Z\end{array}\quad\begin{array}[]{rlllllll}&Y&Y&Y&I&(X)^{k-6}&Z&\\\
&&&&&&Y&Z\\\
\hline\cr-&Y_{j}&Y&Y&I&(X)^{k-6}&X&Z\end{array}\quad\begin{array}[]{rlllllll}&Y&Y&Z&X&(X)^{k-5}&Z\\\
&&&X&X&&\\\
\hline\cr&Y_{j}&Y&Y&I&(X)^{k-5}&Z\end{array}\quad\begin{array}[]{rlllllll}&Y&Y&Y&Y&Y&(X)^{k-6}&Z\\\
&&&&Y&Z&&\\\ \hline\cr&Y_{j}&Y&Y&I&X&(X)^{k-6}&Z\end{array},$ (S174)
$\displaystyle\begin{array}[]{rlllllllllll}&&Y&(X)^{n}&Y&Y&I&(X)^{k-6-n}&Z\\\
&Y&Z&&&&&&\\\
\hline\cr&Y_{j}&X&(X)^{n}&Y&Y&I&(X)^{k-6-n}&Z\end{array}\quad\begin{array}[]{rlllllllllll}&Y&(X)^{n+1}&Y&Y&I&(X)^{k-7-n}&Z&\\\
&&&&&&&Y&Z\\\ \hline\cr&Y_{j}&(X)^{n+1}&Y&Y&I&(X)^{k-7-n}&X&Z\end{array}$
(S181)
$\displaystyle\begin{array}[]{rlllllllllll}&Y&(X)^{n+1}&Y&Z&X&(X)^{k-6-n}&Z\\\
&&&&X&X&&\\\
\hline\cr&Y_{j}&(X)^{n+1}&Y&Y&I&(X)^{k-6-n}&Z\end{array}\quad\begin{array}[]{rlllllllllll}&Y&(X)^{n+1}&Y&Y&Y&Y&(X)^{k-7-n}&Z\\\
&&&&&Y&Z&&\\\ \hline\cr&Y_{j}&(X)^{n+1}&Y&Y&I&X&(X)^{k-7-n}&Z\end{array},$
(S188) $\displaystyle\begin{array}[]{rlllllll}&&Y&(X)^{k-6}&Y&Y&I&Z\\\
&Y&Z&&&&&\\\
\hline\cr&Y_{j}&X&(X)^{k-6}&Y&Y&I&Z\end{array}\quad\begin{array}[]{rlllllll}&Y&(X)^{k-5}&Y&Z&X&Z\\\
&&&&X&X&\\\
\hline\cr&Y_{j}&(X)^{k-5}&Y&Y&I&Z\end{array}\quad\begin{array}[]{rlllllll}&Y&(X)^{k-6}&Y&Y&X&Y\\\
&&&&&X&X\\\ \hline\cr-&Y_{j}&(X)^{k-6}&Y&Y&I&Z\end{array},$ (S198)

respectively. Because the sum of the left-hand sides of Eq. (S144) $\times
J^{xx}/J^{yz}$ and of Eqs. (S145)–(S148) (by choosing the site $j$
appropriately) become zero, we have

$\displaystyle 0$
$\displaystyle=J^{yz}q_{Y_{1}(X)^{k-6}YYIZ}^{(k-1)}+\sum_{n=0}^{k-7}J^{yz}\bigl{(}-q_{Y_{n+1}(X)^{k-n-6}YYI(X)^{n}Z}^{(k-1)}+q_{Y_{n+2}(X)^{k-n-7}YYI(X)^{n+1}Z}^{(k-1)}\bigr{)}$
$\displaystyle\hskip
10.0pt+J^{yz}\bigl{(}-q_{Y_{k-5}YYI(X)^{k-6}Z}^{(k-1)}-q_{X_{k-4}YI(X)^{k-5}Z}^{(k-1)}\bigr{)}+\bigl{(}J^{yz}q_{X_{k-4}YI(X)^{k-5}Z}^{(k-1)}-J^{xx}q_{Z_{k-3}I(X)^{k-4}Z}^{(k-1)}\bigr{)}$
$\displaystyle\hskip 10.0pt+J^{xx}q_{Z_{k-3}I(X)^{k-4}Z}^{(k-1)}$ (S199)
$\displaystyle=-(k-1)\frac{(J^{xx})^{2}}{J^{yz}}q_{Y_{1}(X)^{k-2}Z}^{(k)}.$
(S200)

Thus we obtain the following proposition:

######  Proposition 5.

For $3\leq k\leq N/2$, the solution of Eq. (S19) satisfies

$\displaystyle q_{Y_{1}(X)^{k-2}Z}^{(k)}=0.$ (S201)

By combining Propositions 1 and 2, we have

$\displaystyle q_{\bm{A}^{k}_{j}}^{(k)}=0\quad\text{for all }j\text{ and
}\bm{A}^{k}.$ (S202)

This means that $\hat{Q}$ is a $(k-1)$-local conserved quantity. Applying the
same argument to $k-1$, $k-2$,…, and $3$-local conserved quantity, we can show
that any $k$-local conserved quantity with $k\leq N$ have to be a $2$-local
conserved quantity.

As the third part of the proof, it is straightforward to show the following
proposition:

######  Proposition 6.

Any $2$-local conserved quantity $\hat{Q}$ can be written as

$\displaystyle\hat{Q}=a\hat{H}+b\hat{I},$ (S203)

with arbitrary constants $a,b\in\mathbb{R}$.

From Eq. (S202) and Proposition 6, we obtain Theorem 3.B. ∎

##  IV Level spacing statistics of Model 1

In the paragraph below Theorem 3 of the main text, we explained that Model 1
is expected to be nonintegrable. To confirm this expectation, we numerically
investigate the level spacing statistics [44, 45] of Model 1 by exact
diagonalization. Figure S1 plots the distribution of the ratio of consecutive
level spacings [45]

$\displaystyle
r=\min\Bigl{\\{}\frac{E_{n+2}^{q}-E_{n+1}^{q}}{E_{n+1}^{q}-E_{n}^{q}},\frac{E_{n+1}^{q}-E_{n}^{q}}{E_{n+2}^{q}-E_{n+1}^{q}}\Bigr{\\}},$
(S204)

constructed from eigenenergy $E_{n}^{q}$ in the eigenspace of translation with
momentum $q=2\pi/N$ 777To remove the possibility of accidental discrete
symmetries, we avoid eigenspace of special momentum such as $0$ and $\pi$.
(sorted in descending order). We set the parameters as
$J^{xy}=e,J^{yx}=1,J^{yz}=\pi,J^{zy}=0,h^{y}=\ln 7$ 888We can always take
$J^{zy}=0$ by an appropriate rotation around $y$ axis.. The plot is well
described by the Gaussian unitary ensemble distribution [45] (dashed line) and
well separated from the Poisson distribution [45] (dotted line). This shows
that Model 1 (with above mentioned parameters) has no nontrivial local
conserved quantity, implying nonintegrability of the model.

Figure S1: Distribution of the ratio of consecutive level spacings of Model 1
in Theorem 2 in the main text. We use eigenenergies in the subspace of
momentum $2\pi/N$.

##  V Proof of Theorem 4

###### Proof.

Let $\vec{\sigma}=(\sigma_{1},\sigma_{2},\cdots,\sigma_{N})$ be a bit string
of length $N$. Then, using $\vec{\sigma}$, we define the the computational
basis as

$\displaystyle\ket{\vec{\sigma}}$
$\displaystyle=\ket{\sigma_{1}}\otimes\ket{\sigma_{2}}\otimes\cdots\otimes\ket{\sigma_{N}}.$
(S205)

As a preparation, we clarify the properties of the Hamiltonian satisfying the
assumptions of the theorem. Since Pauli strings form an orthogonal basis of
operators on the whole Hilbert space, the Hamiltonian $\hat{H}$ can be
uniquely expressed as a linear combination of them:

$\displaystyle\hat{H}$
$\displaystyle=\sum_{X(\subset\Lambda)}\sum_{\vec{\mu}\in\\{x,y,z\\}^{X}}J_{X}^{\vec{\mu}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}.$ (S206)

With this notation, $\hat{H}_{\mathrm{int}}$ in the main text is written as

$\displaystyle\hat{H}_{\mathrm{int}}$ $\displaystyle=\sum_{\text{$X\cap
L\neq\emptyset$ and $X\cap
R\neq\emptyset$}}\sum_{\vec{\mu}\in\\{x,y,z\\}^{X}}J_{X}^{\vec{\mu}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}},$ (S207)

where $L=\\{1,2,\cdots,N/2\\}$ and $R=\\{N/2+1,N/2+2,\cdots,N\\}$. Since,
$\hat{H}$ is translation invariant by assumption, we have

$\displaystyle\hat{H}-\hat{H}_{\mathrm{int}}=\hat{H}_{N/2,\mathrm{OBC}}\otimes\hat{I}_{R}+\hat{I}_{L}\otimes\hat{H}_{N/2,\mathrm{OBC}}.$
(S208)

Here, $\hat{H}_{N/2,\mathrm{OBC}}$ is the Hamiltonian for the same system of
length $N/2$, but with open boundary conditions rather than periodic boundary
conditions:

$\displaystyle\hat{H}_{N/2,\mathrm{OBC}}$ $\displaystyle=\sum_{X\subset
L}\sum_{\vec{\mu}\in\\{x,y,z\\}^{X}}J_{X}^{\vec{\mu}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}.$ (S209)

Under the complex conjugation with respect to the computational basis, the
Pauli string behaves as

$\displaystyle\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}\longmapsto(-1)^{P_{X}^{\vec{\mu}}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}},$ (S210)

where $P_{X}^{\vec{\mu}}=|\\{\mu_{j}|j\in X,\mu_{j}=y\\}|$ is the number of
Pauli matrices along the $y$-direction, $\hat{\sigma}_{j}^{y}$, within the
Pauli string. Thus, the complex conjugation transforms the Hamiltonian as

$\displaystyle\hat{H}$
$\displaystyle=\sum_{X(\subset\Lambda)}\sum_{\begin{subarray}{c}\vec{\mu}\in\\{x,y,z\\}^{X}\\\
\text{s.t. $P_{X}^{\vec{\mu}}$ is
even}\end{subarray}}J_{X}^{\vec{\mu}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}+\sum_{X(\subset\Lambda)}\sum_{\begin{subarray}{c}\vec{\mu}\in\\{x,y,z\\}^{X}\\\
\text{s.t. $P_{X}^{\vec{\mu}}$ is
odd}\end{subarray}}J_{X}^{\vec{\mu}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}$ $\displaystyle\longmapsto\qquad\hat{H}^{*}$
$\displaystyle=\sum_{X(\subset\Lambda)}\sum_{\begin{subarray}{c}\vec{\mu}\in\\{x,y,z\\}^{X}\\\
\text{s.t. $P_{X}^{\vec{\mu}}$ is
even}\end{subarray}}J_{X}^{\vec{\mu}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}-\sum_{X(\subset\Lambda)}\sum_{\begin{subarray}{c}\vec{\mu}\in\\{x,y,z\\}^{X}\\\
\text{s.t. $P_{X}^{\vec{\mu}}$ is
odd}\end{subarray}}J_{X}^{\vec{\mu}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}.$ (S211)

Since the expansion in terms of Pauli strings is unique, for $\hat{H}$ to be a
real matrix in the computational basis (i.e., $\hat{H}=\hat{H}^{*}$),
$J_{X}^{\vec{\mu}}$ must be zero when $P_{X}^{\vec{\mu}}$ is odd. Hence, we
obtain

$\displaystyle\hat{H}_{N/2,\mathrm{OBC}}=\sum_{X\subset
L}\sum_{\begin{subarray}{c}\vec{\mu}\in\\{x,y,z\\}^{X}\\\ \text{s.t.
$P_{X}^{\vec{\mu}}$ is even}\end{subarray}}J_{X}^{\vec{\mu}}\bigotimes_{j\in
X}\hat{\sigma}_{j}^{\mu_{j}}.$ (S212)

Therefore, $\hat{H}_{N/2,\mathrm{OBC}}$ is also a real matrix in the
computational basis.

We now proceed to prove Eq. (14). The EAP state $\ket{1;00}$ can be expanded
in the computational basis for subsystems $L$ and $R$ as

$\displaystyle\ket{1;00}\propto\sum_{\vec{\sigma}}\ket{\vec{\sigma}}_{L}\otimes\ket{\vec{\sigma}}_{R}.$
(S213)

Thus, using Eq. (S208), we have

$\displaystyle\ket{\tilde{\beta}}$
$\displaystyle\propto\sum_{\vec{\sigma}}e^{-\frac{1}{4}\beta\hat{H}_{N/2,\mathrm{OBC}}}\ket{\vec{\sigma}}\otimes
e^{-\frac{1}{4}\beta\hat{H}_{N/2,\mathrm{OBC}}}\ket{\vec{\sigma}}$
$\displaystyle=\sum_{\vec{\sigma}}e^{-\frac{1}{4}\beta\hat{H}_{N/2,\mathrm{OBC}}}\ket{\vec{\sigma}}\otimes\left(\sum_{\vec{\sigma}^{\prime}}{\ket{\vec{\sigma}^{\prime}}\bra{\vec{\sigma}^{\prime}}}\right)e^{-\frac{1}{4}\beta\hat{H}_{N/2,\mathrm{OBC}}}\ket{\vec{\sigma}}$
$\displaystyle=\sum_{\vec{\sigma},\vec{\sigma}^{\prime}}e^{-\frac{1}{4}\beta\hat{H}_{N/2,\mathrm{OBC}}}\ket{\vec{\sigma}}\otimes\braket{\vec{\sigma}^{\prime}}{e^{-\frac{1}{4}\beta\hat{H}_{N/2,\mathrm{OBC}}}}{\vec{\sigma}}\ket{\vec{\sigma}^{\prime}}.$
(S214)

Since $\hat{H}_{N/2,\mathrm{OBC}}$ is a real matrix with respect to the
computational basis, it holds that

$\displaystyle\braket{\vec{\sigma}}{e^{-\frac{1}{4}\beta\hat{H}_{N/2,\mathrm{OBC}}}}{\vec{\sigma}^{\prime}}^{*}=\braket{\vec{\sigma}}{e^{-\frac{1}{4}\beta\hat{H}_{N/2,\mathrm{OBC}}}}{\vec{\sigma}^{\prime}}$
(S215)

for any $\ket{\vec{\sigma}}$ and $\ket{\vec{\sigma}^{\prime}}$. Substituting
this into Eq. (S214), we obtain

$\displaystyle\ket{\tilde{\beta}}$
$\displaystyle\propto\sum_{\vec{\sigma},\vec{\sigma}^{\prime}}e^{-\frac{1}{4}\beta\hat{H}_{N/2,\mathrm{OBC}}}\ket{\vec{\sigma}}\bra{\vec{\sigma}}e^{-\frac{1}{4}\beta\hat{H}_{N/2,\mathrm{OBC}}}\ket{\vec{\sigma}^{\prime}}\otimes\ket{\vec{\sigma}^{\prime}}=\sum_{\vec{\sigma}^{\prime}}e^{-\frac{1}{2}\beta\hat{H}_{N/2,\mathrm{OBC}}}\ket{\vec{\sigma}^{\prime}}\otimes\ket{\vec{\sigma}^{\prime}}.$
(S216)

Therefore, for any observable $\hat{O}$ defined on the subsystem $L$, we get

$\displaystyle\braket{\tilde{\beta}}{\hat{O}}{\tilde{\beta}}=\mathrm{Tr}[\hat{\rho}_{N/2,\mathrm{OBC}}^{\mathrm{can}}\hat{O}],$
(S217)

where $\hat{\rho}_{N/2,\mathrm{OBC}}^{\mathrm{can}}$ is the Gibbs state for
$\hat{H}_{N/2,\mathrm{OBC}}$. Hence the thermodynamic limit yields

$\displaystyle\lim_{N\to\infty}\braket{\tilde{\beta}}{\hat{O}}{\tilde{\beta}}=\lim_{N\to\infty}\mathrm{Tr}[\hat{\rho}_{N/2,\mathrm{OBC}}^{\mathrm{can}}\hat{O}].$
(S218)

In the thermodynamic limit, the Gibbs state converges to the KMS state
regardless of whether periodic or open boundary conditions are imposed. Since
we are now considering a one-dimensional system, there exists a unique KMS
state at finite temperature [48, 49]. Consequently, expectation values of
local observables in the Gibbs state do not depend on boundary conditions in
the thermodynamic limit. Thus, using Eq. (13), we finally obtain

$\displaystyle\lim_{N\to\infty}\braket{\beta}{\hat{O}}{\beta}=\lim_{N\to\infty}\mathrm{Tr}[\hat{\rho}^{\mathrm{can}}\hat{O}].$
(S219)

∎

##  VI Degeneracy

###  VI.1 Degeneracy in Model 1

According to Theorem 2, the EAP state $\ket{1;00}$ is an energy eigenstate
with an eigenvalue $E=0$ of Model 1 for arbitrary parameters. Without loss of
generality, we can take $J^{zy}=0$ by an appropriate rotation around $y$ axis,
so we set the parameters of Model 1 as
$J^{xy}=e,J^{yx}=1,J^{yz}=\pi,J^{zy}=0,h^{y}=\ln 7$. By the exact
diagonalization, we find that $\ket{1;00}$ is doubly degenerate in the zero-
momentum sector for $N=10,12,14,16$. Let us investigate entanglement
properties of states in this eigenspace, which we will write as
$\mathcal{H}_{k=0,E=0}$. Let $\ket{\perp}$ denote the state orthogonal to
$\ket{1;00}$ in $\mathcal{H}_{k=0,E=0}$. All states in $\mathcal{H}_{k=0,E=0}$
can be expressed as a linear combination of $\ket{1;00}$ and $\ket{\perp}$:

$\displaystyle\ket{\theta,\lambda}=\sqrt{1-\lambda}\ket{1;00}+e^{2\pi
i\theta}\sqrt{\lambda}\ket{\perp}\qquad(0\leq\lambda\leq 1,0\leq\theta\leq
1).$ (S220)

First, we investigate the bipartite entanglement in $\ket{\perp}$
(corresponding to the case of $\lambda=1$). We plot in Fig. S3 the
entanglement entropy $S_{A}$ between a subsystem $A=\\{1,2,\cdots,\ell\\}$ of
length $\ell$ and its complement as a function of $\ell$. It can be seen that
$\ket{\perp}$ is almost maximally entangled, but is different from EAP states.

Next, we confirm that all states in the eigenspace $\mathcal{H}_{k=0,E=0}$ are
maximally entangled states. To investigate the coefficient of the volume-law
scaling, we compute the entanglement entropy of $\ket{\theta,\lambda}$ between
a subsystem of length $1$ and its complement for various $\lambda$ and
$\theta$ and show in Fig. S3 the deviation from the coefficient of the
maximally entangled state, $1-S_{A=\\{1\\}}/\log 2$. It can be observed that
for all states in $\mathcal{H}_{k=0,E=0}$, the volume-law coefficients are
significantly close to the maximal coefficient.

Thus, there are not any low entangled states in the eigenspace, and hence the
EAP state is not a superposition of such states.

Figure S2: Entanglement entropy of the EAP state and its orthogonal degenerate
state $\ket{\perp}$ between a subsystem $A=\\{1,2,\cdots,\ell\\}$ of length
$\ell$ and its complement as a function of $\ell$ for Model 1 with $N=10$. We
set the parameters as $J^{xy}=e,J^{yx}=1,J^{yz}=\pi,J^{zy}=0,h^{y}=\ln 7$.

Figure S3: Deviation of the volume-law coefficient of the entanglement entropy
of zero-energy eigenstates $\ket{\theta,\lambda}$ defined by Eq. (S220) in the
zero-momentum sector of Model 1 from that of the maximally entangled state. We
set the parameters as $J^{xy}=e,J^{yx}=1,J^{yz}=\pi,J^{zy}=0,h^{y}=\ln 7$ and
$N=10$.

###  VI.2 Nondegenerate Hamiltonian with next-nearest-neighbor interactions

In this subsection, by extending the Hamiltonian (S16) to the next-nearest-
neighbor interacting one, we provide a Hamiltonian having an EAP state as an
eigenstate that is nondegenerate in the corresponding momentum sector.

Suppose that $\hat{H}$ is translation invariant and satisfies
$J_{X}^{\vec{\mu}}=0$ for any subset $X$ with $D(X)\geq 4$. Then, it can be
characterized by $48$ coupling constants, $J^{\mu\lambda\nu}$, $J^{\mu\nu}$
and $h^{\mu}$ ($\mu,\nu=x,y,z$ and $\lambda=x,y,z,0$, where
$\hat{\sigma}^{0}:=\hat{1}$). From Theorem 1 of the main text, it is
straightforward to show that, when $N/2$ is a multiple of $2$, the EAP states
$\ket{2;11,10}$ and $\ket{2;10,11}$ are eigenstates of $\hat{H}$ if and only
if $\hat{H}$ can be written as

$\displaystyle\hat{H}$
$\displaystyle=\sum_{j=1}^{N}\bigl{(}J^{xzx}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{z}\hat{\sigma}_{j+2}^{x}+J^{yzy}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{z}\hat{\sigma}_{j+2}^{y}+J^{xzy}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{z}\hat{\sigma}_{j+2}^{y}+J^{yzx}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{z}\hat{\sigma}_{j+2}^{x}$
$\displaystyle\hskip
36.0pt+J^{zzz}\hat{\sigma}_{j}^{z}\hat{\sigma}_{j+1}^{z}\hat{\sigma}_{j+2}^{z}+J^{xx}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{x}+J^{yy}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{y}+J^{xy}\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{y}+J^{yx}\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{x}+h^{z}\hat{\sigma}_{j}^{z}\bigr{)}.$
(S221)

Here all parameters are arbitrary, and hence it is an extension of Eq. (S16).

Because the EAP states $\ket{2;11,10}$ and $\ket{2;10,11}$ are related to each
other by translation $\hat{\mathcal{T}}$ as

$\displaystyle\hat{\mathcal{T}}\ket{2;11,10}$ $\displaystyle=+\ket{2;10,11},$
(S222) $\displaystyle\hat{\mathcal{T}}\ket{2;10,11}$
$\displaystyle=-\ket{2;11,10},$ (S223)

their superposition $(\ket{2;11,10}\pm i\ket{2;10,11})/\sqrt{2}$ is included
in the eigenspace of translation with the momentum $k=\pm\pi/2$. Therefore, we
investigate degeneracy of energy eigenvalues in the subspace of $k=\pi/2$. We
set $h^{z}=1$, $J^{xx}=e$, $J^{xy}=\pi$, $J^{yx}=\ln 2$, $J^{yy}=\ln 3$,
$J^{xzx}=\ln 5$, $J^{xzy}=\ln 7$, $J^{yzx}=\ln 11$, $J^{yzy}=\ln 13$, and
$J^{zzz}=\ln 17$. By exact diagonalization, we numerically find that, at least
for $N=8,12,16,20$, the eigenvalue $E=0$ is nondegenerate in the subspace of
$k=\pi/2$.

We also verify the nonintegrability of model (S221). Because of the existence
of $J^{zzz}$, model (S221) is not mapped to a free fermionic (integrable)
system by the Jordan-Wigner transformation. We confirm nonintegrability of
model (S221) by calculating the distribution of the ratio of consecutive level
spacings, as in Fig. S1. Figure S4 plots the distribution constructed from
energy eigenvalues in the eigenspace of translation $\hat{\mathcal{T}}$ and
$\mathbb{Z}_{2}$ transformation $\bigotimes_{j=1}^{N}\hat{\sigma}_{j}^{z}$
with momentum $k=\pi/2$ and parity $\mathcal{P}=\pm 1$. This plot is well
described by the Gaussian unitary ensemble distribution [45] (dashed line) and
well separated from the Poisson distribution [45] (dotted line), indicating
the nonintegrability of the model.

Figure S4: Distribution of the ratio of consecutive level spacings of model
(S221). We use eigenenergies in the subspace of momentum $k=\pi/2$ and parity
$\mathcal{P}=\pm 1$ (regarding rotation by $\pi$ around $z$-axis). The
parameters are given below Eq. (S223).

Note that, because the interference term between $\ket{2;11,10}$ and
$\ket{2;10,11}$ does not affect the expectation values of local observables
whose support size are less than $N/2$, any state described by a linear
combination of $\ket{2;11,10}$ and $\ket{2;10,11}$ is locally
indistinguishable from the maximally mixed state.

Combining all results of this subsection, we can say that the state
$(\ket{2;11,10}\pm i\ket{2;10,11})/\sqrt{2}$ is a thermal eigenstate of the
nonintegrable Hamiltonian (S221), and is nondegenerate in the corresponding
momentum sector 999The results of exact diagonalization also show that, at
least for $N=8,16,20$, degeneracy at $E=0$ in the whole Hilbert space remains
only two. This means that all eigenstates of the Hamiltonian (S221) with
eigenenergy $E=0$ are given by (linear combinations of) the EAP states
$\ket{2;11,10},\ket{2;10,11}$..