††thanks: These two authors contributed equally.††thanks: These two authors
contributed equally.

# Boundary Quantum Phase Transitions in the Spin $\frac{1}{2}$ Heisenberg
Chain with Boundary Magnetic Fields

Parameshwar R. Pasnoori  Department of Physics, University of Maryland,
College Park, Maryland 20742 USA  Laboratory for Physical Sciences, 8050
Greenmead Dr, College Park, Maryland 20740 USA  Junhyun Lee  Department of
Physics and Astronomy, Center for Materials Theory, Rutgers University,
Piscataway, NJ 08854-8019 USA  J. H. Pixley  Department of Physics and
Astronomy, Center for Materials Theory, Rutgers University, Piscataway, NJ
08854-8019 USA  Center for Computational Quantum Physics, Flatiron Institute,
New York, New York 10010, USA  Natan Andrei  Department of Physics and
Astronomy, Center for Materials Theory, Rutgers University, Piscataway, NJ
08854-8019 USA  Patrick Azaria  Laboratoire de Physique Théorique de la
Matière Condensée, Sorbonne Université and CNRS, 4 Place Jussieu, 75252 Paris,
France<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider the spin $\frac{1}{2}$ Heisenberg chain with boundary magnetic
fields and analyze it using a combination of Bethe ansatz and density matrix
renormalization group (DMRG) techniques. We show that the system exhibits
several different ground states which depend on the orientation of the
boundary magnetic fields. When both the boundary fields take equal values
greater than a critical field strength, each edge in the ground state
accumulates a fractional spin which saturates to spin $\frac{1}{4}$, which is
similar to systems exhibiting symmetry protected topological phases (SPT).
Unlike in SPT systems, the fractional boundary spin in the Heisenberg spin
chain is not a genuine quantum number since the variance of the associated
operator does not vanish, this is due to the absence of a bulk gap. The system
exhibits high energy bound states when the boundary fields take values greater
than the critical field. All the excitations in the system can be sorted out
into towers whose number depends on the number of bound states exhibited by
the system. As the boundary fields are varied, in addition to the ground state
phase transition, we find that the system may undergo an eigenstate phase
transition (EPT) where the number of towers of the Hilbert space changes. We
further inquire how the EPT reflects itself on local ground state properties
by computing the magnetization profile $\langle S^{z}_{j}\rangle$ using DMRG.
We identify a clear qualitative change from low edge fields to high edge
fields when crossing the critical field. We though are unable to conclude on
the basis of our data that EPT corresponds to a genuine phase transition in
the ground state.

##  I Introduction

The Heisenberg model is one of the most celebrated models in condensed matter
and statistical physics. It lies at the cornerstone of our understanding of
many physical phenomenon which, besides magnetism, consists of integrability
Bethe (1931); Sklyanin (1988); Yang and Yang (1966a), many body localization
Pal and Huse (2010); Abanin  _et al._ (2019) and out of equilibrium dynamics
Ilievski  _et al._ (2018, 2015); Dupont and Moore (2020). Thanks to available
analytic and numerical methods the model is quite well understood in one
spatial dimension. This is partly due to the fact that the one dimensional
spin $S=1/2$ Heisenberg model, which is also known as the XXX spin chain, is
integrable. At the same time, the model can be probed experimentaly either in
solid-state compounds comprising quasi one-dimensional spin chains in KCuF3
Steiner  _et al._ (1976); Satija  _et al._ (1980); Nagler  _et al._ (1991);
Tennant  _et al._ (1995); Maillet (2007) or more recently in ultra-cold atom
realizations of the spinful Bose-Hubbard model Jepsen  _et al._ (2021, 2022).
Since it was first solved by Bethe Bethe (1931), the spin chain with periodic
boundary conditions has been very well studied. Both the ground state and the
low energy excitations properties are well understood Yang and Yang (1966a, b,
c); Babelon  _et al._ (1983); Boos and Korepin (2001). The system is non
magnetic and supports massless spin $\frac{1}{2}$ excitations named spinons.
Besides this, integral representations of correlation functions have been
obtained Kitanine  _et al._ (2000); Shiroishi and Takahashi (2005); Takahashi
(1999). Spin chains with open boundaries have also been intensely studied
after the Yang-Baxter algebra was generalized to systems with open boundaries
by Sklyannin and Cherednik Sklyanin (1988); Cherednik (1984). For the XXX spin
chain with open boundaries the ground state, bulk excitations and physical
boundary $S-$matrices have been found Grisaru  _et al._ (1995). More generally
the effects of boundary fields Cao  _et al._ (2002); Qiao  _et al._ (2021);
Sun  _et al._ (2019) have also been investigated and spin chains with non
diagonal boundary fields have been solved Wang  _et al._ (2015).

In this work we shall be interested in the XXX spin chain with magnetic fields
at its edges. Although the subject has been studied to some extent, for a
quantum impurity Wang (1997) as well as for a classical one (e.g. with
boundary magnetic fields) we find that some issues remain to be clarified and
explored in the light of new developments related with one dimensional
topological phases and eigenstate phase transitions. To the best of our
knowledge the results that will be presented in this work have not been found
before.

The Hamiltonian of the XXX spin chain with boundary magnetic fields is given
by

$\displaystyle
H=\sum_{j=1}^{N-1}\sum_{\alpha=x,y,z}\sigma^{\alpha}_{j}\cdot\sigma^{\alpha}_{j+1}+h_{L}\;\sigma^{z}_{1}+h_{R}\;\sigma^{z}_{N},$
(1)

where $\sigma^{\alpha}_{j}$ are the Pauli matrices acting on the spin space at
site $j$, $h_{L}$ and $h_{R}$ are boundary magnetic fields acting at sites
$j=1$ and $j=N$ respectively, and $N$ is the number of sites. The boundary
magnetic fields break the $SU(2)$ spin symmetry down to the $U(1)$ group of
rotation around the “$z$” axis. There exists no other symmetries except when
$h_{L}=h_{R}$, where the model displays space parity invariance $\mathbb{P}$.
Despite this, the system possesses a useful isometry obtained by
simultaneously flipping all spins as well as reversing the orientation of the
boundary fields:

$\displaystyle
H\rightarrow\prod_{i=1}^{N}\sigma^{x}_{i}H\sigma^{x}_{i},\;\;h_{L}\rightarrow-
h_{L},\;\;h_{R}\rightarrow-h_{R}.$ (2)

The latter isometry is a symmetry of the phase diagram of the model in the
plane $(h_{L},h_{R})$. The Hamiltonian (1) is integrable by the method of the
Algebraic Bethe Ansatz. The related Bethe equations has been first obtained in
Ref. Sklyanin (1988); Grisaru  _et al._ (1995); Skorik and Saleur (1995);
Kapustin and Skorik (1996) where a preliminary analysis was given. In the
present work we shall extend their analysis and present a more thorough
picture of the phase diagram associated with (1), thus providing important new
results that have not been, to the best of our knowledge, present in the
literature. But before going into more details let us first discuss
qualitatively our main results.

We first discuss the ground state properties of the model. We present in the
Fig. (1) the ground state phase diagram in the plane $(h_{L},h_{R})$ and for
an even number of sites. A similar analysis can be made for an odd number of
sites as given in section (II). There are four different possible ground
states when labelled by the conserved total “$z$” component of the spin
operator,

$S^{z}=\frac{1}{2}\sum_{j=1}^{N}\sigma_{j}^{z}.$ (3)

When $h_{L}h_{R}<0$ the ground state is unique and has total spin $S^{z}=0$
whereas, contrarily to what was found in Ref.Kapustin and Skorik (1996), in
the quadrants $h_{L}h_{R}>0$ we find that the ground state is doubly
degenerated; each one having spin $S^{z}=0$ and $S^{z}=+1$ or $S^{z}=-1$
depending on whether $h_{L(R)}$ is negative or positive. In the later cases,
the degeneracy is found to be the consequence of the existence of two static
spinons (with infinite rapidity) with spins $\pm\frac{1}{2}$ in the ground
state. To get a better understanding of the ground state structure we have
performed extensive density matrix renormalization group (DMRG) calculations
and calculated the magnetization profile $\langle S^{z}_{j}\rangle$ in the
ground state. Overall we find that the edge magnetic fields induce a spin
polarization close to the two edges which extends into the bulk in a power law
fashion. Furthermore, we find that the corresponding spin accumulations at the
edges are fractional and take the values $\pm\frac{1}{4}$ (opposite to the
orientation of the edge field) at large fields $|h|>2$. The situation at hand
is similar to what happens in gapless symmetry protected topological (SPT)
superconductors Keselman and Berg (2015); Keselman  _et al._ (2018); Pasnoori
_et al._ (2020, 2021) where the edge states Hilbert space is exhausted by
spin-$\frac{1}{4}$ operators. However in the present case, due to the
existence of massless spinon bulk excitations, such an operator do not
represent a genuine fractional sharp quantum observable since, as we show in
the section (IV), its variance is not zero.

Figure 1: The figure shows various ground states occurring for different
values of the boundary magnetic fields. The ground state is unique in the
second and fourth quadrants. In the first and third quadrants, the ground
state contains a spinon with infinite rapidity and whose spin is oriented
either in the positive or negative $z$ direction, resulting in a two fold
degenerate ground state.

The second topic we shall discuss in this work concerns the structure of the
excited states and its relation to the existence of bound states localized at
the edges. One of the hallmarks of the boundary physics induced by the edge
magnetic fields is the existence of boundary bound states localized close to
the edges where the fields are applied. As previously found in Refs.Skorik and
Saleur (1995); Kapustin and Skorik (1996), when the magnitude of fields are
large enough, i.e: when $|h_{L(R)}|\geq 2$, the system hosts bound states with
energy

$\displaystyle m_{L,R}=2\pi/\sin(\frac{\pi}{h_{L,R}}),$ (4)

which carry a spin $\frac{1}{2}$, whose spin orientation is along the boundary
fields at each edge. Contrarily to the zero energy edge states in SPT massless
superconductors, the bound states in the XXX model are high-energy states
whose energies are always above the one spinon branch of massless bulk
excitations, and as we shall see, their existence has important consequences
on the structure of the Hilbert space. We present in Fig. 2 the bound state
phase diagram of the model (1). In each quadrant the different phases are
sorted out as a function of the number of bound states: the $A$ sub-phases
support two bound states (one at each edge), the $B$ sub-phases support one
bound state at either the right or the left edge, whereas in the $C$ sub-
phases there are no bound states. When compared to the ground state phase
diagram we see that each quadrant is split into three different sub-phases
named $A$, $B$ or $C$. Although in all these sub-phases the ground state has
the same total spin $S^{z}$, they differ by the structure of the high energy
states. We show that each bound state generates a whole tower of excited
states that can be built upon it. Hence, the Hilbert space is comprised of a
certain number of towers which depends on the number of bound states exhibited
by the system. In the case of the gapped regime of the spin $\frac{1}{2}$ XXZ
spin chain with open boundaries a similar structure of the Hilbert space is
found, where it was demonstrated in Fendley (2016) that the Hilbert space is
comprised of two towers of degenerate eigenstates which leads to the emergence
of a strong zero energy Majorana operator (which commute with the Hamiltonian)
which map these pairs of states.

In each of the $A,B,C$ sub-phases, the direct sum of these towers span the
complete Hilbert space. When crossing the boundaries between any two of these
sub-phases, since the number of bound states exhibited by the system and hence
the number of towers of the excited states changes, an Eigenstate Phase
Transition occurs which involves a full reorganization of the Hilbert space. A
similar phenomenon is observed in gapless SPT superconductors Pasnoori  _et
al._ (2021).

In summary, we find that similar to the systems exhibiting SPT Pasnoori  _et
al._ (2021) the XXX spin chain exhibits several phases as a function of edge
magnetic fields. We find that the total spin of the ground state is not enough
to completely characterize these phases which differ also by the structure of
the Hilbert space. The later is linked to the number of bound states at the
edges which generate towers of excited states that together span the Hilbert
space. As a consequence, on top of the phase transition corresponding to a
change in the ground state, there exists eigenstate phase transitions
involving the change in the number of towers of the Hilbert space.

The paper is organized as follows. We present our results obtained from the
Bethe ansatz for the ground state and the excited states in the sections II
and III respectively. Section IV is dedicated to the DMRG analysis of the
ground state properties. We finally discuss our results in the section V.

##  II  Ground state phase diagram

As said above, the Hamiltonian (1) is integrable by the method of the
Algebraic Bethe ansatz Grisaru  _et al._ (1995); Wang  _et al._ (2015) for
arbitrary values of the boundary fields. Its ground-state as well as
excitations are obtained from the Bethe equations

$\displaystyle\left(\frac{\lambda_{j}-i/2}{\lambda_{j}+i/2}\right)^{2N}\left(\frac{\lambda_{j}+i(\frac{1}{2}-p_{L})}{\lambda_{j}-i(\frac{1}{2}-p_{L})}\right)\left(\frac{\lambda_{j}+i(\frac{1}{2}-p_{R})}{\lambda_{j}-i(\frac{1}{2}-p_{R})}\right)$
$\displaystyle=\prod_{j\neq
k=1}^{M}\left(\frac{\lambda_{j}-\lambda_{k}-i}{\lambda_{j}-\lambda_{k}+i}\right)\left(\frac{\lambda_{j}+\lambda_{k}-i}{\lambda_{j}+\lambda_{k}+i}\right),$
(5)

where we have introduced $p_{L/R}=1/h_{L/R}$ as the boundary parameters. The
eigenstates of the Hamiltonian are labelled by $M\in\mathbb{N}$ Bethe roots
$\lambda_{j=1,...,M}$ which are solutions of Eq.(II) and have energy

$\displaystyle
E=-\sum_{j=1}^{M}\frac{2}{\lambda_{j}^{2}+\frac{1}{4}}+N-1+h_{L}+h_{R}.$ (6)

The corresponding total spin $S^{z}$ of a state is related to the integer $M$
through the relation (7)

$\displaystyle S^{z}=\pm\left(\frac{N}{2}-M\right).$ (7)

Where $\pm$ corresponds to reference state with all spin up and down
respectively. We have obtained from (II) the ground state phase diagram as a
function of the boundary magnetic fields $h_{L(R)}$ and for both an even and
an odd number of sites. Before going into more details let us first review
briefly the situation at zero fields. In this case the ground state depends on
the parity of $N$ as follows: For even $N$ it is non degenerated and has total
spin $S^{z}=0$ whereas for odd $N$ it is twofold degenerated, each ground
state having total spin $S^{z}=\pm\frac{1}{2}$. The latter degeneracy is due
to the presence of spin-$\pm\frac{1}{2}$ spinons, which have energy

$\displaystyle E_{\theta}=\frac{2\pi}{\cosh(\pi\theta)},$ (8)

and is zero in the limit of infinite rapidity $\theta\rightarrow\infty$. Let
us see now how this scheme is modified by the presence of non zero boundary
fields. The situation further depends on the parity of $N$.

Figure 2: The figure shows various phases occurring for different values of
the boundary magnetic fields. The $\it{A}$ sub-phases exhibit two boundary
bound-states, one at each edge. In $\it{B}$ sub-phases there exists one
boundary bound state at either the left or the right edges. In $\it{C}$ sub-
phases boundary bound-states do not exist.

###  II.1 Ground state for odd number of sites

Since the total number of sites is odd the spins of the ground states have to
be half integers. The phase diagram can be broadly divided into four quadrants
based on the direction of the boundary magnetic fields as shown in the Fig.
(3). In the upper right quadrant, when both the magnetic fields point towards
the positive $z$ direction, and independently of their magnitudes
$|h_{L(R)}|$, the ground state is unique and has a total spin
$S^{z}=-\frac{1}{2}$. In the ground state, that we shall label

$|-\frac{1}{2}\rangle,$ (9)

the total magnetization is due to a static spin configuration which account
for the total spin $-\frac{1}{2}$. In the lower right quadrant, in which
$h_{L}>0$ and $h_{R}<0$, the ground state is doubly degenerated and carry
total spins $S^{z}=\pm\frac{1}{2}$

$|-\frac{1}{2}\rangle\;{\rm and}\;|+\frac{1}{2}\rangle$ (10)

In contrast with the previous case, the spins of the ground states here is due
to the presence of spin-$\pm\frac{1}{2}$ spinons with infinite rapidity
$\theta\rightarrow\infty$ (8). The situation in the two other quadrants, i.e.
the lower left and upper left ones, can be obtained by using the isometry (2)
and reversing the total spin quantum number $S^{z}\rightarrow-S^{z}$. The
ground states are then found to be

$|+\frac{1}{2}\rangle\;{\rm and}\;|-\frac{1}{2}\rangle,$ (11)

in the upper left quadrant, i.e. when $h_{L}<0$ and $h_{R}>0$ and

$|+\frac{1}{2}\rangle,$ (12)

in the lower left quadrant when $h_{L}<0$ and $h_{R}<0$.

###  II.2 Ground state for even number of sites

In this case the spins of the ground states are always integers. As shown in
the Fig.(1) in the upper right quadrant, i.e. when both $h_{L(R)}>0$, the
ground state is doubly degenerated and have total spins $S^{z}=0,-1$
represented by

$|0\rangle\,{\rm and}\;|-1\rangle.$ (13)

The double degeneracy of the ground state is due to the presence of
spin-$\pm\frac{1}{2}$ spinons with infinite rapidity. These spinons have to be
added on top of a background static spin configuration contributing to a total
spin $-\frac{1}{2}$ in such a way that the total spins of the ground states
are integers $S^{z}=0,-1$. Using the isometry (2) we deduce immediately that
in the lower left quadrant the ground states are given by

$|0\rangle\,{\rm and}\;|+1\rangle,$ (14)

and also contain spin-$\pm\frac{1}{2}$ spinons with infinite rapidity on top
of a background spin $\frac{1}{2}$ configuration. This is to be true
independently of the magnitudes of the fields $|h_{L(R)}]$. When the boundary
fields point towards opposite direction, like in the two upper left and lower
right quadrants the ground state is unique with total spin $S^{z}=0$,

$|0\rangle,$ (15)

and do not contains spinons.

Figure 3: Ground state phase diagram as a function of the boundary fields for
an odd number of sites. The ground state is unique in the first and third
quadrants. In the second and fourth quadrants, the ground state contains a
spinon with infinite rapidity and whose spin is oriented either in the
positive or negative $z$ direction, resulting in a two fold degenerate ground
state.

As we have seen the ground state phase diagram exhibits four distinct phases
depending solely on the orientations of the boundary fields. In each of the
four quadrants defined by the sign of $h_{L}$ and $h_{R}$ the ground state
degeneracy depends on the parity of $N$. It is two-fold degenerate when
$h_{L}h_{R}<0$ for $N$ odd and when $h_{L}h_{R}>0$ for $N$ even. In all other
cases the ground state is non degenerate contrary to what was found in
Kapustin and Skorik (1996). Overall, our understanding of the spin quantum
numbers in the different phases relies on a static background spin
distribution on top of which spins $\pm\frac{1}{2}$ spinons may or may not be
added. Independently of the parity of the number of sites $N$, the background
spin distribution contributes to a total spin
$S^{z}_{B}=-\operatorname{sgn}(h_{L})\;\frac{1}{2}$ when $h_{L}h_{R}>0$
whereas $S^{z}_{B}=0$ in the opposite case when $h_{L}h_{R}<0$. Such a
background spin structure is due to the presence of the boundary fields
$h_{L},h_{R}$ which are expected to induce a spin accumulation close to the
edges. We shall return to this point in the section (IV) when we shall study
the ground state properties in more detail. For the time being we shall argue
that the phase structure induced by the presence of the boundary fields is
much richer than the one we have just presented. When considering the whole
structure of the Hilbert space, which calls for a detailed description of the
excited states, we shall show that each of the four quadrants, $(h_{L}\gtrless
0,h_{R}\gtrless 0)$, splits into four distinct sub-phases where the excited
states organizes into different towers. This is the consequence of the well-
known fact that the edge fields, when large enough, induce boundary bound
states that are exponentially localized close to the edges.

##  III Excited states

Before going into the description of the full solution of (LABEL:hamiltonian),
let us review quickly how, in the case of a periodic chain without fields, the
structure of the excitation spectrum depends on the parity of the number of
sites $N$. When $N$ is even the ground-state is a singlet and the excitations
are obtained by adding an even number of propagating holes or spinons. The
spinons carry spins $\pm\frac{1}{2}$ and have energy (8). Since physical
excitations correspond to flipping a certain number of spins in the chain they
carry integer spins $S^{z}\in\mathbb{Z}$ and therefore the spinons always come
in pairs in an even chain. For the spin chain with an odd number of sites, the
ground state is two-fold degenerate with total spins $S^{z}=\pm\frac{1}{2}$.
Each of the ground-states contains a spinon, with zero energy in the
thermodynamical limit, and rapidity $\theta\rightarrow\infty$. In contrast
with the even chain case, single spinon excitations with a finite rapidity,
i.e: $\theta\neq\infty$, and energy (8) are allowed when $N$ is odd. All other
excitations are then obtained by adding an even number of spinons to the above
states. Hence the total number of spinons in the odd spin chain is always an
odd integer.

This scheme is to be modified in an open chain with boundary fields for which
the ground-state and, more importantly, the very structure of the Hilbert
space of excitations strongly depend on the boundary fields $h_{L/R}$. As we
shall see, when the boundary fields are strong enough (i.e. when
$|h_{L/R}|\geq 2$) their main effect is to stabilize bound-states which are
localized at either the left or the right edge. These bound-states have a
finite, i.e. non zero, energy above the ground-state,

$m_{L,R}=\frac{2\pi}{\sin(\frac{\pi}{h_{L,R}})},$ (16)

and carry a spin $\frac{1}{2}$ which points towards the direction of the
boundary field at each edge. The bound states at the left and the right edges
are exponentially localized as $\sim e^{-\kappa_{L}x}$ and $\sim
e^{-\kappa_{R}(N-x)}$ respectively (see Appendix B), where

$\kappa_{j}=\log\left(\frac{h_{j}+1}{h_{j}}\right),\;\;\;j=L,R.$ (17)

When the bound states exist, they generate independent towers of excited
states on top of the ground-state one. All these towers of states eventually
span the whole Hilbert space.

We distinguish between three regions, or sub-phases, $A,B$ and $C$ depending
on the number of localized bound-states in the spectrum. In the region $A$
both boundary field strengths exceed a critical value $|h_{L/R}|\geq 2$ and
there exists two boundary bound-states localized at both ends of the chain.
Depending on the relative orientations of the fields $h_{L/R}$ with respect to
the ”z” axis we further distinguish between four sub-phases $A_{j=(1,2,3,4)}$.
In the region $B$ only one boundary field strength exceeds the critical value
and there exists a single bound-state which is localized at either the left or
the right edge. Taking into account the orientations of the fields we end up
with eight sub-phases $B_{j}$ $j=(1...8)$. Finally in region $C$, both
$|h_{L/R}|<2$ and there are no localized bound-states; the four sub-phases
$C_{j=(1,2,3,4)}$ account for all possible boundary fields orientations. The
phase diagram is depicted in Fig.(2).

In the following we shall present our results for the ground-states as well as
the Hilbert space structures in each phase. Since, as with the PBC case
discussed above, the spectral properties are very sensitive to the evenness of
the number of sites $N$, we shall discuss separately both even and odd chains.

###  III.1 $\it{A}$ sub-phases

We start with the $A$ sub-phases where two boundary bound-states are
stabilized. The four $A_{j=(1,2,3,4)}$ sub-phases correspond to the domains of
boundary fields $(h_{L}\geq 2,h_{R}\geq 2),(h_{L}\leq-2,h_{R}\geq
2),(h_{L}\leq-2,h_{R}\leq-2)$ and $(h_{L}\geq 2,h_{R}\leq-2)$ respectively. In
the following we shall distinguish between odd end even chains and discuss
separately the sub-phases $A_{j=(1,3)}$ and $A_{j=(2,4)}$.

####  III.1.1 Odd number of sites

##### The $A_{1}$ and $A_{3}$ sub-phases.

In these cases both boundary magnetic fields point towards the same direction:
along the positive $z$ axis for the $A_{1}$ sub-phase and negative $z$ axis
for the $A_{3}$ sub-phase. Both cases are related by the isometry (2).
Qualitatively speaking, in the sub-phases $A_{1,3}$ and for $N$ odd, the
boundary magnetic fields are not frustrating in the sense that in the Ising
limit of (LABEL:hamiltonian) the ground-state would exhibit perfect
antiferromagnetic order.

In the $A_{1}$ sub-phase we find that the ground-state is unique and has a
total spin $S^{z}=-\frac{1}{2}$. We accordingly label the ground-state in this
phase by

$|-\frac{1}{2}\rangle,$ (18)

and denote by $E_{0}$ its energy. The expression of $E_{0}$ as a function of
$h_{L,R}$ is given in the Appendix (see Eq.(59)). We notice that due to the
presence of the boundary fields the spin $-\frac{1}{2}$ of the ground-state is
not carried by a spinon in contrast with the periodic chain with $N$ odd. It
is rather the consequence of a static spin density distribution. We shall
discuss this topic in more detail in the next section. Similarly to the case
of periodic boundary conditions, one can build up excitations in the bulk on
top of this ground state by adding an arbitrary even number of spinons, bulk
strings and quartets Destri and Lowenstein” (1982). These bulk excitations
built on top of the state $|-\frac{1}{2}\rangle$ form a tower of excited
states that we shall denote the ground-state tower.

As said above in the $A$ sub-phases there exists two boundary bound-state
solutions exponentially localized at either the left or the right edge. In the
language of the Bethe ansatz they correspond to purely imaginary solutions of
(II) (see Appendix A). These bound-states carry a spin $\frac{1}{2}$, whose
spin orientation is along the boundary fields at each edge, and have an energy
(4) Since the bound-states carry a spin half, in order to add a bound-state to
the ground-state one also needs to add a spinon. This spinon may have spin
$+\frac{1}{2}$ or $-\frac{1}{2}$ and an arbitrary rapidity $\theta$. The
energy cost in the process is

$E_{0}+m_{L,R}+E_{\theta},$ (19)

and is minimal when $\theta\rightarrow\infty$. The corresponding states

$|\pm\frac{1}{2}\rangle_{L}\;{\rm and}\;|\pm\frac{1}{2}\rangle_{R}\;,$ (20)

have total spins $S^{z}=\pm\frac{1}{2}$ and energies $E_{0}+m_{L}$ and
$E_{0}+m_{R}$. The lowest excited states above (20) consist of spinon branches
with energies given by (19) and $\theta\neq\infty$. On top of these, the
states (20) generate, each, a tower of excited states obtained by adding an
arbitrary even number of spinons, bulk strings and quartets. In both the left
and right towers, built upon (20), a localized bound-state at the left and the
right edge is present and the number of spinon excitations is always odd.

On top of the above three towers there exists a fourth one which correspond to
states which host two bound-states. The state with the lowest energy in this
tower is obtained by adding a localized bound-state at the left and the right
edges to the ground-state (18). Since in the process the total spin of the
state is shifted by $1$, no spinon is required. The resulting state

$|+\frac{1}{2}\rangle_{LR}\;,$ (21)

which has a total spin $S^{z}=\frac{1}{2}$ and an energy $E_{0}+m_{L}+m_{R}$,
generates a tower of excited states that comprises an arbitrary even number of
spinons, bulk strings and quartets. The number of spinon states in the whole
tower is always even. We thus see that, in the $A_{1}$ sub-phase, the whole
Hilbert space can be split into four towers generated by the states (18, 20)
and (21) as illustrated in figure 4. On top of the ground-state tower which
governs the low-energy physics, the remaining three towers contain at least
one bound-state at the edges and are high-energy states. In particular, we
notice that in the $A_{1}$ sub-phase there exists excitations which contain a
single spinon, and although the system is massless, their minimum energy is
greater than the boundary gap $m_{L}$ or $m_{R}$. These four towers can be
labelled by the bound state parities

$\mathcal{P}_{L,R}=(-1)^{\mathcal{N}_{L,R}},$ (22)

where $\mathcal{N}_{L,R}$ correspond to number of bound states at the left and
right edges respectively.

The situation in the $A_{3}$ sub-phase can be described in the very same way
as above. Using the isometry (2), we can obtain all the states in the sub-
phase $A_{3}$ starting from the states in the sub-phase $A_{1}$ by reversing
the sign of the total spin $S^{z}$ of the states. Hence, we obtain four towers
of states in the sub-phase $A_{3}$ generated by the states
$|+\frac{1}{2}\rangle$, $|\pm\frac{1}{2}\rangle_{L,R}$ and
$|-\frac{1}{2}\rangle_{LR}$ at energies $E_{0},E_{0}+m_{L,R}$ and
$E_{0}+m_{L}+m_{R}$.

Figure 4: The figure illustrates the structure of the Hilbert space in the $A$
sub-phases. In each $A_{j}$ ($j=1,2,3,4$) sub-phase, there exist four towers
of excited states, each labelled by the bound state parities
[$\mathcal{P}_{L}$, $\mathcal{P}_{R}$]. The lowest energy states in each tower
have energies close to the dashed lines corresponding to the energies $E_{0}$,
$E_{0}+m_{L}$, $E_{0}+m_{R}$ and $E_{0}+m_{L}+m_{R}$.

##### The $A_{2}$ and $A_{4}$ sub-phases.

In these cases the boundary fields are frustrating for $N$ odd in the sense
discussed above. As we shall see in these sub-phases the Hilbert space is also
split into four towers of states corresponding to the presence of boundary
bound-states. However, since the boundary magnetic fields at the two edges
point toward opposite directions, the nature of these towers differ from the
ones described above. Consider for instance the $A_{2}$ sub-phase in which the
left boundary field points towards the negative $z$ axis while the one at the
right boundary points in the opposite direction. In this case we find that the
ground-state is two-fold degenerated, each one containing a spinon (but no
bound-state) with spin $\pm\frac{1}{2}$ and rapidity
$\theta\rightarrow\infty$. These two states, i.e:

$|\pm\frac{1}{2}\rangle,$ (23)

have total spin $S^{z}=\pm\frac{1}{2}$ corresponding to the spin of the
spinon, and generate a tower of excited state. It is obtained by adding an
arbitrary even number of spinons, bulk strings and quartets on top of the two
spin $\pm\frac{1}{2}$ massless spinons branches with spectrum (8) and rapidity
$\theta\neq\infty$. In contrast with the $A_{1,3}$ sub-phases the ground-state
tower contains an odd number of spinons.

Just as in the sub-phase $A_{1}$, there exists two boundary bound-state
solutions one at each edge. The bound state’s spin is always oriented along
the boundary magnetic field. Hence, in the sub-phase $A_{2}$ the bound-state
localized at the left edge has spin $-\frac{1}{2}$ whereas the bound-state
localized at the right edge has spin $+\frac{1}{2}$. We find that in order to
add the bound-state at the left edge with spin $-\frac{1}{2}$ one has to
remove the $-\frac{1}{2}$ spinon at $\theta=\infty$ in the
$|-\frac{1}{2}\rangle$ ground-state (23). The resulting state has total spin
$S^{z}=-\frac{1}{2}$ and energy $E_{0}+m_{L}$. Similarly adding a spin
$+\frac{1}{2}$ bound-state at the right edge requires to remove the spin
$\frac{1}{2}$ spinon from the ground-state $|+\frac{1}{2}\rangle$ (23). The
resulting state has total spin $S^{z}=+\frac{1}{2}$ and energy $E_{0}+m_{R}$.
The two states with a bound-state at either the left or right edge

$|-\frac{1}{2}\rangle_{L}\;{\rm and}\;|+\frac{1}{2}\rangle_{R},$ (24)

generate, each, a tower of excited states upon adding an arbitrary even number
of spinons, bulk strings and quartets. In these two towers the number of
spinons in every state is always even.

Finally, the fourth tower is obtained by adding a bound-state at each edge to
the two ground-states (23). The total spin of the resulting state does not
change since the two, left and right, bound-states have opposite spins. We
obtain the states

$|\pm\frac{1}{2}\rangle_{LR},$ (25)

which have an energy $E_{0}+m_{L}+m_{R}$ and generate a tower of excited
states. It is obtained by adding even number of spinons, bulk strings and
quartets. In this tower the number of spinons is always odd.

Using the symmetry 2, we can obtain all the states in the sub-phase $A_{4}$
from the states in the sub-phase $A_{2}$ by reversing their spins. The Hilbert
space in the sub-phase $A_{4}$ can be similarly sorted out in terms of four
towers of states built upon the states $|\pm\frac{1}{2}\rangle$,
$|+\frac{1}{2}\rangle_{L},|-\frac{1}{2}\rangle_{R}$ and
$|\pm\frac{1}{2}\rangle_{LR}$ with energies $E_{0},E_{0}+m_{L,R}$ and
$E_{0}+m_{L}+m_{R}$.

Table 1: Energies and bound state parities of the ground state and the lowest energy states corresponding to each tower in all the $A$ sub-phases for odd number of sites is shown below. The subscripts $L,R$ denote the location of the bound states at the left or the right boundary. Phase |  State |  Energy-$E_{0}$ |  $\mathcal{P}_{L}$ |  $\mathcal{P}_{R}$
---|---|---|---|---
| $|-\frac{1}{2}\rangle$ |  $0$ (g.s) | 1 | 1
| $|\pm\frac{1}{2}\rangle_{R}$ | $m_{R}$ | 1 | -1
$A_{1}$ | $|\pm\frac{1}{2}\rangle_{L}$ | $m_{L}$ | -1 | 1
| $|\frac{1}{2}\rangle_{L,R}$ | $m_{L}+m_{R}$ | 1 | 1
| $|\frac{1}{2}\rangle$ |  $0$ (g.s) | 1 | 1
| $|\pm\frac{1}{2}\rangle_{R}$ | $m_{R}$ | 1 | -1
$A_{3}$ | $|\pm\frac{1}{2}\rangle_{L}$ | $m_{L}$ | -1 | 1
| $|-\frac{1}{2}\rangle_{L,R}$ | $m_{L}+m_{R}$ | 1 | 1
| $|-\frac{1}{2}\rangle_{L}$ | $m_{L}$ | -1 | 1
| $|\pm\frac{1}{2}\rangle$ |  $0$ (g.s) | 1 | 1
$A_{2}$ | $|\pm\frac{1}{2}\rangle_{L,R}$ | $m_{L}+m_{R}$ | -1 | -1
| $|\frac{1}{2}\rangle_{R}$ | $m_{R}$ | 1 | -1
| $|-\frac{1}{2}\rangle_{R}$ | $m_{R}$ | 1 | -1
| $|\pm\frac{1}{2}\rangle$ |  $0$ (g.s) | 1 | 1
$A_{4}$ | $|\pm\frac{1}{2}\rangle_{,L,R}$ | $m_{L}+m_{R}$ | -1 | -1
| $|\frac{1}{2}\rangle_{L}$ | $m_{L}$ | -1 | 1

####  III.1.2 Even number of sites

When the number of sites is even the frustrating effect of the magnetic fields
is reversed as compared to the $N$ odd case. The boundary fields are
frustrating in sub-phases $A_{1,3}$ while non-frustrating in the sub-phases
$A_{2,4}$.

##### The $A_{1}$ and $A_{3}$ sub-phases.

In the sub-phase $A_{1}$ we find that the ground-state is two-fold
degenerated. It does not contain bound-states but does contain a spinon with
rapidity $\theta\rightarrow\infty$ and spins $\pm\frac{1}{2}$. Despite this,
since $N$ is even, the total spins of the two degenerate ground-states have to
be integers. Indeed, as it comes out from our exact solution the two ground-
states have total spins $S^{z}=0$ and $S^{z}=-1$. Our interpretation of this
fact is that the two ground-states with spin $S^{z}=0$ and $S^{z}=-1$ contain
a spin $+\frac{1}{2}$ and a spin $-\frac{1}{2}$ spinon respectively on top of
a static background spin $-\frac{1}{2}$ distribution corresponding to the
ground state in $A_{1}$ sub-phase when $N$ is odd. In the following we denote
these two ground-states by

$|0\rangle\;{\rm and}\;|-1\rangle.$ (26)

The ground-state tower of excitated states comprises spin $\pm\frac{1}{2}$
massless spinon states with energy $E_{0}+m_{R}+E_{\theta}$ and finite
rapidity $\theta\neq\infty$. The rest of the tower is then obtained by adding
an arbitrary even number of spinons, bulk strings and quartets. In this tower
the number of spinon states is always odd.

Starting from one of the two ground-states (26), one may add a bound-state at
either the left or the right edge. To this end one needs to remove the spin
$\pm\frac{1}{2}$ spinon. The resulting total spin is then the sum of the
bound-state spin $+\frac{1}{2}$ with that of the static background spin
$-\frac{1}{2}$ distribution mentioned above. As a result, we end up with two
states of total spin $S^{z}=0$. The corresponding states with the bound-state
at the left or the right edge are denoted

$|0\rangle_{L}\;{\rm and}\;|0\rangle_{R},$ (27)

and have energies $E_{0}+m_{L}$ and $E_{0}+m_{R}$. Each of these two states
generates a tower of excited states. In these towers the number of spinon
states is always even.

The fourth tower is obtained from the ground-states (26) by adding a bound-
state at each edge. Since the change of total spin is $1$ there is no need to
add or remove a spinon. In the process we obtain two degenerate states, with
total spins $S^{z}=1$ and $S^{z}=0$ and energy $E_{0}+m_{L}+m_{R}$,

$|1\rangle_{LR}\;{\rm and}\;|0\rangle_{LR},$ (28)

that host spin $\pm\frac{1}{2}$ spinons with infinite rapidity as in the
ground-states. The fourth tower of excited states comprises, as in the ground-
state tower, spin $\pm\frac{1}{2}$ spinon states. These states have energy
$E_{0}+m_{L}+m_{R}+E_{\theta}$ and are gapped high energy states. The
remaining states of this towers are then built up by adding an even number of
spinons, bulk strings and quartets. The number of spinon states is always odd.

Similar to the odd number of sites case, using the symmetry (2), we can obtain
all the states in the phase $A_{3}$ starting from the states in the phase
$A_{1}$ described above.

##### The $A_{2}$ and $A_{4}$ sub-phases.

In the sub-phase $A_{2}$ we find that the ground-state is non-degenerated

$|0\rangle,$ (29)

and has total spin $S^{z}=0$ with energy $E_{0}$. Starting from this ground
state we can add a bound-state at the left edge with spin $-\frac{1}{2}$. As
already emphasized one also needs to add a spinon, with infinite rapidity and
zero energy, for the total spin shift to be an integer. Depending on the
spinon spin, which can be either $\pm\frac{1}{2}$, one ends up with two states

$|-1\rangle_{L},|0\rangle_{L},$ (30)

which have total spins $S^{z}=-1$ and $S^{z}=0$ and energy $E_{0}+m_{L}$. One
may repeat the same line of arguments with the right edge paying attention
that the bound-state’s spin in this case is $+\frac{1}{2}$. The resulting two
states

$|+1\rangle_{R},|0\rangle_{R},$ (31)

hosting a bound-state at the right edge have total spins $S^{z}=1$ and
$S^{z}=0$ and energy $E_{0}+m_{R}$. Each left and right states (30) and (31)
generate two towers of excited states that comprise spin $\pm\frac{1}{2}$
spinons with energies $E_{0}+m_{L,R}+E(\theta)$. The rest of the towers are
obtained by adding even number of spinons, bulk strings and quartets.

The forth tower is obtained from the ground-state (29) by adding a bound-state
with spin $-\frac{1}{2}$ at the left edge and spin $+\frac{1}{2}$ at the right
edge. No spinons are needed in the process and one ends up with a single state
‘

$|0\rangle_{LR},$ (32)

with total spin $S^{z}=0$ and energy $E_{0}+m_{R}+m_{L}$. The latter state
generates also a tower of states with even number of spinons, bulk strings and
quartets.

Using the symmetry (2), similar to the odd number of sites case, we can obtain
all the states in the sub-phase $A_{4}$ starting from the states in the sub-
phase $A_{2}$ described above. The ground state and the lowest energy state
corresponding to each tower in all the $A$ sub-phases for odd and even number
of sites chain are summarized in the tables 1 and 2 respectively.

Table 2: Energies and bound state parities of the ground state and the lowest energy states corresponding to each tower in all the $A$ sub-phases for even number of sites is shown below. Phase |  State |  Energy-$E_{0}$ |  $\mathcal{P}_{L}$ |  $\mathcal{P}_{R}$
---|---|---|---|---
| $|-1\rangle,\;|0\rangle$ | 0 (g.s) | 1 | 1
| $|0\rangle_{R}$ | $m_{R}$ | 1 | -1
$A_{1}$ | $|0\rangle_{L}$ | $m_{L}$ | -1 | 1
| $|1\rangle_{L,R},\;|0\rangle_{L,R}$ | $m_{L}+m_{R}$ | 1 | 1
| $|1\rangle,\;|0\rangle$ | 0 (g.s) | 1 | 1
| $|0\rangle_{R}$ | $m_{R}$ | 1 | -1
$A_{3}$ | $|0\rangle_{L}$ | $m_{L}$ | -1 | 1
| $|-1\rangle_{L,R},\;|0\rangle_{L,R}$ | $m_{L}+m_{R}$ | 1 | 1
| $|0\rangle$ | 0 (g.s) | -1 | -1
| $|-1\rangle_{L},\;|0\rangle_{L}$ | $m_{L}$ | 1 | -1
$A_{2}$ | $|1\rangle_{R},\;|0\rangle_{R}$ | $m_{R}$ | -1 | 1
| $|0\rangle_{L,R}$ | $m_{L}+m_{R}$ | 1 | 1
| $|0\rangle$ | 0 (g.s) | -1 | -1
| $|1\rangle_{L},\;|0\rangle_{L}$ | $m_{L}$ | 1 | -1
$A_{4}$ | $|-1\rangle_{R},\;|0\rangle_{R}$ | $m_{R}$ | -1 | 1
| $|0\rangle_{L,R}$ | $m_{L}+m_{R}$ | 1 | 1

###  III.2 $\it{B}$ phases

####  III.2.1 Odd number of sites

In the $B_{1}$ sub-phase, the ground state has total spin $S^{z}=-\frac{1}{2}$
which corresponds to a static spin distribution and is represented by

$|-\frac{1}{2}\rangle$ (33)

Unlike in the $A$ phases, there exists only a single boundary bound state
solution corresponding to the bound state at the left edge. Starting from the
ground state, this bound state can be added (which has spin
$S^{z}=\frac{1}{2}$) by adding a spinon whose spin orientation can be either
in the positive or negative $z$ direction resulting in the state with total
spin $S^{z}=\pm\frac{1}{2}$ respectively. This state has energy
$E_{0}+m_{L}+E_{\theta}$, and hence has the lowest energy in the limit
$\theta\rightarrow\infty$. It is represented by

$|\pm\frac{1}{2}\rangle_{L}$ (34)

In the sub-phase $B_{2}$, the state which does not contain a bound state at
either edge contains a spinon whose spin orientation is either in the positive
or negative $z$ direction. The energy of this state is $E_{0}+E_{\theta}$ and
thus forms a continuous branch parameterized by $\theta$. The ground state is
obtained in the limit $\theta\rightarrow\infty$ and is represented by

$|\pm\frac{1}{2}\rangle$ (35)

Starting from this ground state one can add a bound state at the left edge
(which has spin $S^{z}=-\frac{1}{2}$) by removing the existing spinon. The
resulting state has energy $E_{0}+m_{L}$ with total spin $S^{z}=-\frac{1}{2}$,
and is represented by

$|-\frac{1}{2}\rangle_{L}.$ (36)

By using the transformation $L\rightarrow R$, the states in the phases $B_{8}$
and $B_{7}$ can be obtained by starting with the states in the phases $B_{1}$
and $B_{2}$ respectively. By using the transformation 2, the states in the
phases $B_{5},B_{6},B_{3}$ and $B_{4}$ can be obtained from the states in the
phases $B_{1},B_{2},B_{7}$ and $B_{8}$ respectively.

Table 3: Energies and local bound state parities of the ground state and the lowest energy states corresponding to each tower in all the $B$ phases for odd number of sites is shown below. Phase |  State |  Energy-$E_{0}$ |  $\mathcal{P}_{L}$ |  $\mathcal{P}_{R}$
---|---|---|---|---
$B_{1}$ | $|-\frac{1}{2}\rangle$ | 0 (g.s) | 1 | 1
| $|\pm\frac{1}{2}\rangle_{L}$ | $m_{L}$ | -1 | 1
$B_{8}$ | $|-\frac{1}{2}\rangle$ | 0 (g.s) | 1 | 1
| $|\pm\frac{1}{2}\rangle_{R}$ | $m_{R}$ | 1 | -1
$B_{2}$ | $|-\frac{1}{2}\rangle_{L}$ | $m_{L}$ | -1 | 1
| $|\pm\frac{1}{2}\rangle$ | 0 (g.s) | 1 | -1
$B_{7}$ | $|\pm\frac{1}{2}\rangle$ | 0 (g.s) | 1 | 1
| $|-\frac{1}{2}\rangle_{R}$ | $m_{R}$ | 1 | -1
$B_{4}$ | $|\frac{1}{2}\rangle$ | 0 (g.s) | 1 | 1
| $|\pm\frac{1}{2}\rangle_{R}$ | $m_{R}$ | 1 | -1
$B_{5}$ | $|\frac{1}{2}\rangle_{L}$ | 0 (g.s) | 1 | 1
| $|\pm\frac{1}{2}\rangle_{L}$ | $m_{L}$ | -1 | 1
$B_{3}$ | $|\pm\frac{1}{2}\rangle$ | 0 (g.s) | 1 | 1
| $|\frac{1}{2}\rangle_{R}$ | $m_{R}$ | 1 | -1
$B_{6}$ | $|\pm\frac{1}{2}\rangle$ | 0 (g.s) | 1 | 1
| $|\frac{1}{2}\rangle_{L}$ | $m_{L}$ | -1 | 1

####  III.2.2 Even number of sites

In the phase $B_{1}$, the state with no bound states at both the edges is two
fold degenerate. It contains a spinon on top of the static spin distribution
of the ground state in the phase $B_{1}$ corresponding to odd number of sites
case. The spin orientation of the spinon can be either in the positive or
negative $z$ direction which results in a doubly degenerate state with total
spin $S^{z}=0,-1$. This state has energy $E_{0}+E_{\theta}$ and thus forms a
continuous branch which is parameterized by $\theta$. The ground state is
obtained in the limit $\theta\rightarrow\infty$ and is represented by

$|0\rangle,|-1\rangle.$ (37)

We can add the bound state at the left edge (with spin $S^{z}=\frac{1}{2}$) to
the ground state by removing the existing spinon. This results in a state

$|0\rangle_{L}$ (38)

with total spin $S^{z}=0$ with energy $E_{0}+m_{L}$.

In the phase $B_{2}$, the state which does not contain bound state at either
edge has total spin $S^{z}=0$ and has energy $E_{0}$. It is represented by

$|0\rangle.$ (39)

We can add the bound state at the left edge (with spin $S^{z}=-\frac{1}{2}$)
by adding a spinon with spin oriented either in the positive or negative
z-direction and hence resulting in a doubly degenerate state with total spin
$S^{z}=-1,0$. This state has energy $E_{0}+E_{\theta}+m_{L}$, and hence the
lowest energy of this state corresponds to the limit $\theta\rightarrow\infty$
and is represented by

$|0\rangle_{L},|-1\rangle_{L}$ (40)

Similar to the odd number of sites case, the states in the phases $B_{8}$ and
$B_{7}$ can be obtained by starting with the states in $B_{1}$ and $B_{2}$
respectively, by making the transformation $L\rightarrow R$. By using the
transformation 2, the states in the phases $B_{5},B_{6},B_{3}$ and $B_{4}$ can
be obtained from the states in the phases $B_{1},B_{2},B_{7}$ and $B_{8}$
respectively.

Table 4: Energies and local bound state parities of the ground state and the lowest energy states corresponding to each tower in all the $B$ phases for even number of sites is shown below. Phase |  State |  Energy-$E_{0}$ |  $\mathcal{P}_{L}$ |  $\mathcal{P}_{R}$
---|---|---|---|---
$B_{1}$ |  $|-1\rangle$, $|0\rangle$ | 0 (g.s) | 1 | 1
| $|0\rangle_{L}$ | $m_{L}$ | -1 | 1
$B_{8}$ |  $|-1\rangle$, $|0\rangle$ | 0 (g.s) | 1 | 1
| $|0\rangle_{R}$ | $m_{R}$ | 1 | -1
$B_{2}$ |  $|-1\rangle_{L}$, $|0\rangle_{L}$ | $m_{L}$ | -1 | 1
| $|0\rangle$ | 0 (g.s) | 1 | 1
$B_{7}$ |  $|-1\rangle_{R}$, $|0\rangle_{R}$ | $m_{R}$ | 1 | -1
| $|0\rangle$ | 0 (g.s) | 1 | 1
$B_{4}$ |  $|1\rangle$, $|0\rangle$ | 0 (g.s) | 1 | 1
| $|0\rangle_{R}$ | $m_{R}$ | 1 | -1
$B_{5}$ |  $|1\rangle$, $|0\rangle$ | 0 (g.s) | 1 | 1
| $|0\rangle_{L}$ | $m_{L}$ | -1 | 1
$B_{3}$ |  $|1\rangle_{R}$, $|0\rangle_{R}$ | $m_{R}$ | 1 | -1
| $|0\rangle$ | 0 (g.s) | 1 | 1
$B_{6}$ |  $|1\rangle_{L}$, $|0\rangle_{L}$ | $m_{L}$ | -1 | 1
| $|0\rangle$ | 0 (g.s) | 1 | 1

Unlike in the $A$ sub-phases where there exists bound states at both the
edges, we have seen that in $B$ sub-phases there exists only one bound state
at either the left or the right edge. Similar to the $A$ sub-phases,
excitations can be built up starting from the ground state and from the state
containing a bound state either at the left or the right edge by adding even
number of spinons, strings and quartets. This leads to the Hilbert space in
each $B$ sub-phase consisting of only two towers. For example, in the phase
$B_{1}$, the two towers have the bound state parities
$\mathcal{P}_{L}=1,\mathcal{P}_{R}=1$ and
$\mathcal{P}_{L}=-1,\mathcal{P}_{R}=1$, whereas in the $B_{8}$ phase they
correspond to $\mathcal{P}_{L}=1,\mathcal{P}_{R}=1$ and
$\mathcal{P}_{L}=1,\mathcal{P}_{R}=-1$. The ground states and the lowest
energy states corresponding to the two towers in all the $B$ sub-phases are
summarized in the tables 3, 4.

###  III.3 $\it{C}$ sub-phases

##### Odd number of sites

In the sub-phases $\it{C_{1}},\it{C_{3}}$, the ground state is represented by

$|\mp\frac{1}{2}\rangle$ (41)

and have total spin $S^{z}=\mp\frac{1}{2}$ respectively, which corresponds to
a static spin distribution. In the sub-phases $\it{C_{2}},\it{C_{4}}$ the
lowest energy state contains a spinon with spin pointing either in the
positive or negative z-direction resulting in a two fold degenerate state
parameterized by rapidity $\theta$. The ground state is obtained in the limit
$\theta\rightarrow\infty$. The spin orientation of the spinon dictates the
total spin $S^{z}=\pm\frac{1}{2}$ of the state. They are represented by

$|\pm\frac{1}{2}\rangle.$ (42)

##### Even number of sites

In the sub-phase $\it{C_{1}}$, the lowest energy state contains a spinon with
rapidity $\theta$ with spin oriented either in the positive or negative
z-direction on top of the static spin distribution of the ground state in the
sub-phase $C_{1}$ corresponding to odd number of sites case. This state is two
fold degenerate and is parameterized by rapidity $\theta$. The ground state is
obtained in the limit $\theta\rightarrow\infty$ and is represented by

$|0\rangle,|-1\rangle$ (43)

with total spin $S^{z}=0$, $S^{z}=-1$ corresponding to the spin orientation of
the spinon which is along the positive and negative $z$ direction
respectively. Similarly, in the sub-phase $\it{C_{3}}$, the lowest energy
state contains a spinon with spin pointed either in the positive or negative
z-direction with rapidity $\theta$ on top of the static spin distribution of
the ground state in the phase $C_{3}$ corresponding to odd number of sites
case. It is two fold degenerate and is parameterized by rapidity $\theta$. The
ground state is obtained in the limit $\theta\rightarrow\infty$ and is
represented by

$|0\rangle,|1\rangle$ (44)

with total spin $S^{z}=0$, $S^{z}=1$ corresponding to the spin orientation of
the spinon which is along the negative and positive $z$ direction
respectively. In the sub-phases $C_{2},C_{4}$, the ground state has total spin
$S^{z}=0$ and is represented by

$|0\rangle.$ (45)

Similar to the $A$ and $B$ sub-phases, in each $C$ subphase, excitations can
be built on top of the ground state by adding even number of spinons, strings
and quartets generating a single tower of excited states which can be labelled
by $\mathcal{P}_{L}=1,\mathcal{P}_{R}=1$.

###  III.4 Eigenstate phase transition

After this rather lenghtly, but complete, description of the excited states
let us now summarize our results. As we saw there exists a critical value of
the edge fields $h_{c}$, $|h_{c}|=2$, at each edge associated with the
existence of an edge bound state. When $|h_{i=(L,R)}|>2$ a localized bound
state is stabilized close to the corresponding edge $i=(L,R)$. The three types
of phases $A_{j}$, $B_{j}$ and $C_{j}$ distinguish themselves by the number of
bound states they support, i.e: two, one and zero. Independently of the parity
of $N$ we showed that in the $A$-type phases the Hilbert space splits into
four towers of excited states while there exists two towers in the $B$-type
phases and only one tower in the $C$-type phases. When compared to the ground
state phase diagrams (see Figs.(1,3)) each quadrant splits into one $C_{j}$
sub-phase, two $B_{j}$ sub-phases and one $A_{j}$ sub-phase as displayed in
the Fig. 2. At this point a natural question arises: what is the nature of the
transition that occurs as one moves from an $A_{j}$ sub-phase to a $B_{j}$
sub-phase or from a $B_{j}$ sub-phase to a $C_{j}$ sub-phase by varying the
edge fields.

Without loss of generality let us fix on quadrant with $h_{L}>0$ and
$h_{R}>0$. Consider first the situation where both $h_{L,(R)}>2$, that is one
sits in the $A_{1}$ sub-phase. Then let the left boundary magnetic field
$h_{L}$ be fixed while the right boundary fields $h_{R}$ is decreased. As
$h_{R}$ is decreased below the critical value $h_{c}=2$, we move into the sub-
phase ${B}_{1}$. The two states which contain the bound state at the right
edge no longer exist. On the boundary between the $A_{1}$ and $B_{1}$ sub-
phases, the energy of the bound state and energy of the spinon with zero
rapidity coincide $m_{R}\sim 2\pi=E_{\theta\rightarrow 0}$. Hence it is
natural to interpret that the bound state at the right edge leaks into the
bulk by taking the form of a spinon with rapidity $\theta\sim 0$. Similarly,
moving from ${A}_{1}$ to ${B}_{8}$ (see Fig. 2), the bound state corresponding
to left boundary leaks into the bulk. Similarly, moving from ${B}_{1}$ to
${C}_{1}$, the value of the left boundary field takes values lesser than
critical value, and hence the bound state present at the left edge leaks into
the bulk in a similar way, resulting in $C_{1}$ having no bound states at
either edge. The same phenomena of bound states leaking into the bulk occurs
as one moves from any $A$ sub-phase into the respective $B$ and $C$ sub-
phases.

Figure 5: The Hilbert space in the $A$ phases is comprised of four towers
whereas it is comprised of two towers in the $B$ phases and a single tower in
the $C$ phases. Figure illustrates the eigenstate phase transitions (EPT) that
occur between $A$, $B$ and $C$ type phases, where the number of towers of the
Hilbert space changes.

More importantly, associated with the appearance or disappearance of localized
bound states is the fact that when one goes from any sub-phase to another, the
whole structure of the Hilbert space changes. The excited states organize
themselves into towers whose number is different in the $A$, $B$ or $C$ type
phases. We saw that the towers are labelled by additional quantum numbers
which are the bound state parities $\mathcal{P}_{L,R}$ (see Eq.(22)). The four
towers in $A$-type phases are labeled by
$(\mathcal{P}_{L},\mathcal{P}_{R})=(\pm 1,\pm 1)$, the two towers in the
$B$-type phases by $(\mathcal{P}_{L},\mathcal{P}_{R})=(\pm 1,+1)$ and
$(\mathcal{P}_{L},\mathcal{P}_{R})=(+1,\pm 1)$ and the unique tower of the
$C$-type phases by $(\mathcal{P}_{L},\mathcal{P}_{R})=(+1,+1)$. It is
interesting to notice that one may have also labeled these towers in the
$A$-type phases by the spinon parity $\mathcal{P}_{s}=(-1)^{n}$, where $n$
denotes the number of spinons in a given eigenstate of the Hamiltonian. The
relation between $\mathcal{P}_{L},\mathcal{P}_{R}$ and $\mathcal{P}_{s}$
depends on the phase and the parity of the number of sites $N$ as follows

$\displaystyle\mathcal{P}_{L}\mathcal{P}_{R}(-1)^{n}(-1)^{N}=(-1)^{k},$ (46)

where $k$ labels the different sub-phases $A_{k}$ as given in Fig.(2). As when
crossing from an $A_{j}$ sub-phase to either a $B_{j}$ or $C_{j}$ sub-phase
the structure of the Hilbert space changes, and we coin the corresponding
phase transition a (boundary) eigenstate phase transition. Such transitions
might be probed through dynamical properties at infinite temperatures that
involve operators localized close enough to the boundaries. We shall elaborate
on this topic in a forthcomming work. At present we shall content ourselves,
in the next section, with the simpler question of how this transition reflects
itself in the ground state properties of local observables.

##  IV Ground state magnetization profile and spin accumulation

Figure 6: The magnetization profile of the ground state in (a) $A1$, (b) $A2$,
(c) $A3$, and (d) $A4$ phase, respectively. The boundary fields are $h_{L}=\pm
4$, $h_{R}=\pm 4$ on a chain of length $N=100$, and the arrow indicates the
direction of the boundary fields.

To this end we shall be interested in the behavior of the local magnetization
profile induced by the edge magnetic fields in the different sub-phases. We
shall use the DMRG method, which is ideal for one-dimensional systems, to
calculate the ground state of Eq. (1) on finite size systems. In particular,
we considered system sizes up to $N=1600$ sites, where a maximum of $1000$
states are kept to keep the truncation error below $10^{-12}$. The DMRG
calculations in this paper are performed using the ITensor Library Fishman
_et al._ (2020). Once the ground state is obtained we compute the spin
expectation value $S^{z}(x_{i})\equiv\langle\sigma_{i}^{z}/2\rangle$ as a
function of position $x_{i}$ for various edge magnetic fields $h_{L},h_{R}$ in
the different phases of the problem.

###  IV.1 Magnetization Profile

From the magnetization obtained from the DMRG calculations, we use the
following ansatz for the magnetization near the boundary:

$\displaystyle
S^{z}(x_{i})=(-1)^{i}\left(A+\frac{B}{\sqrt{x}_{i}}+Ce^{-x_{i}/\xi}\right)+\frac{D}{x_{i}}.$
(47)

We have introduced a constant staggered magnetization ($A$) that vanishes in
the thermodynamic limit as well as the alternating $1/\sqrt{x}$ and uniform
$1/x$ terms which account, to leading order, for the gapless bulk. In the
bosonization language they correspond to the staggered and uniform component
of the magnetization in the long distance limit. In addition to the above
terms, we have also included a term which goes like $\sim e^{-x_{i}/\xi}$ to
account for any exponentially localized spin accumulation.

Overall we find excellent agreement between our DMRG results and its fit (47)
for several values of $h_{L},h_{R}$. We show in Fig. 6, as an example, the
magnetization profile $S^{z}(x_{i})$ deep in the four $A$-phases for boundary
fields $|h_{L}|=|h_{R}|=4$ and system size $N=100$. As expected, the
magnetizations at the boundaries are all opposite to the boundary field
directions and one clearly observes a spin accumulation close to the edges. We
notice though that, since the bulk is gapless, the spin accumulation is not
expected to be sharply localized at the edge. This can be seen in the
magnetization profile which exhibit asymptotic power-law antiferromagnetic
decay sufficiently far away from the edges. Notice that, since in the $A_{1}$
and $A_{3}$ phases the ground state contains a spinon in the bulk for even
chains, a node in the bulk antiferromagnetic configuration [Fig. 6(a)(c)] is
clearly seen. Alternatively, for odd chain lengths the $A_{1}$ and $A_{3}$
phases have no spinons whereas the $A_{2}$ and $A_{4}$ phases do.

Figure 7: (a) Magnetization (dots) and its fit to Eq. (47) (dashed lines) for
the first 20 sites of $N=1000$ chain. (b) The fit-parameters $C$ and $D$ from
Eq. (47), as a function of $h$. Both parameters vanishes at the critical
$h_{c}=2$. The shaded region is where the fitting to Eq. (47) numerically
fails as it is close to $h_{c}=2$. (See Appendix C for more details)

Figure 8: (a) The magnetization for the $N=1000$ chain with critical boundary
fields ($h=h_{c}=2$). The magnetization of a critical chain decays as $\sim
1/\sqrt{x_{i}}$ as can be seen from the good agreement with the
$|S^{z}(x_{i})|=0.413/\sqrt{x_{i}}+0.00385$ obtained from fitting the data.
(b) The fit-parameters $A$ and $B$ from Eq. (47), as a function of boundary
field $h$. (c) The length scale $\xi$ of the exponentially localized boundary
spin, also from Eq. (47), as a function of $h$. The critical point is
$h_{c}=2$ where the length scale diverges.

We now discuss in more detail about our results for the fit in the particular
case $h_{L}=h_{R}=h$, where the system exhibits a $\mathbb{Z}_{2}$ space
parity symmetry. As one varies $h$, this allows us to study the spin
magnetization profile when going from the $C_{1}$ phase to the $A_{1}$ phase
(see Fig 2). We show for example in Fig. 7(a) the magnetization profile for
different values of $h$ in the $A_{1}$ phase and the critical point. Fitting
the DMRG data with the form (47), we can extract the parameters $A$, $B$, $C$
and $D$ as well as the length scale $\xi$ as a function of $h$. The
coefficients $C$ of the exponential term and $D$ of the uniform component
$1/x$ are shown in Fig. 7(b). The constant term $A$ and the coefficient $B$ of
the staggered magnetization component $1/\sqrt{x}$ are shown in Fig. 8(b).
Finally $\xi$ is shown in Fig. 8(c).

From our data we first observe that when $h=h_{c}=2$ the magnetization profile
takes a particular simple form in the thermodynamical limit

$S^{z}(x_{i})\simeq(-1)^{i}0.413/\sqrt{x_{i}},$ (48)

as both $C$ and $D$ are zero (even at a finite size) and $A$ goes to zero as
$N\rightarrow\infty$. We display in Fig. 8(a) our best fit for the
magnetization when $h=h_{c}=2$ and $N=1000$ which shows an almost perfect
$\sim 1/\sqrt{x}$ behavior asymptotically. It is remarkable that for this
value of the edge field the bulk uniform component $1/x$ disappears from the
magnetization profile.

This case $h=h_{c}=2$ seems to play a special role in the magnetization
profile. Indeed we find that both the coefficient $C$ and $D$ change sign when
going accross the $h=h_{c}=2$ point where they vanish. The change in the sign
of $C$ means that the exponential term enhances (diminishes) its contribution
to $S^{z}(x_{i})$ when $h>h_{c}$ ($h<h_{c}$). On the other hand the change of
sign of $D$ can be interpreted as a $\pi$ phase-shift of the uniform component
term in Eq. (47). At the transition, both contributions vanish as we find
$C,D=0$. Another important feature is that coefficients $A$ and $B$ saturate
as one increases $h$ above $h_{c}=2$ as seen in Fig. 8(b). This means that at
magnetic fields $h$ larger than $h_{c}=2$ the constant contribution as well as
the staggered component of the magnetization are insensitive to the edge
magnetic field. As these are the dominant contributions for large $x$, this
means that the magnetizations far from the edges are essentially insensitive
to the edge fields when $h>h_{c}$, in contrast to low fields $h<h_{c}$ where
both $A$ and $B$ wildly varies.

In the light of the discussion given in the precedent section it seems that at
the critical point where the eigenstate phase transition occurs (here between
the $A_{1}$ and $C_{1}$ phases on the $\mathbb{Z}_{2}$ symmetric line
$h_{L}=h_{R}$) a qualitative change also occurs in the magnetization profile
in the ground state. Whether this change corresponds to a genuine critical
point for the ground state properties (reflecting itself into a singular
behavior of the magnetization profile at $h_{c}$) is a non trivial issue.
Indeed, although we find that the length scale $\xi$ diverges as $h$
approaches $h_{c}$, at the same time $C,D\rightarrow 0$ and the numerical
fitting overfits the data in the vicinity of the critical point ($h=h_{c}=2$).
We mark this region with a grey shade in Fig. 7 (b) and Fig. 8 (c) (for more
detail see Appendix C). We therefore find it difficult to conclude that these
coefficients, or length scale $\xi$, serve as order parameters for a genuine
ground state phase transition.

###  IV.2 Boundary spin accumulation

Figure 9: (a) The spin accumulation $S^{z}_{L}(\alpha,N)$ for various system
sizes when $h_{L}=-h_{R}=4$. The inset is a blowup near $\alpha=0$. (b) The
difference between the spin accumulation for system size $N$ and that of the
system of $N=300$. The inset is the same plot in linear scale, and the gray
dashed lines are the conjectured $(1/4)e^{-300\alpha}$. (c) The boundary spin
accumulation obtained from the linear fit of (b) to $Ae^{-300\alpha}$ as a
function of boundary field $h_{L}$.

Due to the edge magnetic fields we naturally expect that some amount of spin
is getting accumulated close to the edges. To calculate the spin accumulation
associated with the edges of the system, we use the following definition of
the spin accumulation on the left boundary $S^{z}_{L}$ as Jackiw  _et al._
(1983):

$\displaystyle S^{z}_{L}=\lim_{\alpha\rightarrow
0}\lim_{N\rightarrow\infty}S^{z}_{L}(\alpha,N),$ (49)

where,

$\displaystyle S^{z}_{L}(\alpha,N)=\sum_{i}e^{-x_{i}\alpha}S^{z}(x_{i}).$ (50)

Note the order of the limit is relevant and it is important to take the
thermodynamic limit first. Otherwise, the $\alpha\rightarrow 0$ limit removes
the cutoff and the result becomes merely the total $S^{z}$ of the system.
Also, $\alpha\gg 1/L$ should be satisfied for the cutoff to be meaningful. In
the following we shall compute $S^{z}_{L}$ in both the $C_{2}$ and $A_{2}$
phases where $h_{L}=-h_{R}=-h$. The reason for this is to keep the odd parity
of the system which simplifies the finite-size conjectures which will follow.

To obtain $S^{z}_{L}$ systematically, we infer
$\lim_{N\rightarrow\infty}S^{z}_{L}(\alpha,N)\equiv S^{z}_{L}(\alpha)$ from
the finite size calculations. Fig. 9(a) shows the $S^{z}_{L}(\alpha,N)$ as a
function of $\alpha$ for different system sizes, when $h=4$. One important
observation is that, as the system size grows $S^{z}_{L}(\alpha,N)$ converges
to the $S^{z}_{L}(\alpha)$ curve. Therefore, each $S^{z}_{L}(\alpha,N)$ is
converged to $S^{z}_{L}(\alpha)$ for $\alpha$ larger than a certain value of
$\alpha_{N}$, which $\lim_{N\rightarrow\infty}\alpha_{N}=0$. We conjecture the
leading difference of the finite $S^{z}_{L}(\alpha,N)$ and infinite
$S^{z}_{L}(\alpha)$ as:

$\displaystyle
S^{z}_{L}(\alpha)=S^{z}_{L}(\alpha,N)-\frac{1}{4}e^{-N\alpha}+\cdots,$ (51)

for $h\geq h_{c}$. This equation is suggestive that it consists of the
proposed spin (1/4) and the only length scale ($N$), and directly implies the
fractional 1/4-spin. Taking the $\alpha\rightarrow 0$ limit to Eq. (51), the
first term vanishes because $S^{z}_{L}(0,N)=\sum_{i}S^{z}(x_{i})$ for any
finite $N$ and the total $S^{z}$ is zero for the even chain in $A_{2}$ phase.
The remaining terms give the fractionalized $S^{z}_{L}=\pm 1/4$ per Eq. (49)
where the sign depends on the direction of the boundary field.

To see how this work let us give an example of our finite size scaling
procedure in the particular case $h=4$. In Fig. 9(b) we plot
$S^{z}_{L}(\alpha,300)-S^{z}_{L}(\alpha,N)$ for various $N$ for the same
parameters in Fig. 9(a). As we expect from our conjecture [Eq. (51)], the
difference converges to $\frac{1}{4}e^{-N\alpha}$ (dashed line) for large $N$.
To quantify the numerical value of $S_{L}^{z}$, we find the best exponential
fit to the plots similar to Fig. 9(b) for different $h$ values, and obtain the
spin accumulation $S_{L}^{z}(h)$ as the overall coefficient.

Our final result for the spin accumulation $S_{L}^{z}$ as a function of $h$ is
shown in Fig. 9(c). Our results for $S^{z}_{L}$ are consistent with
$S^{z}_{L}=1/4$ for $h\geq h_{c}=2$ and decreases for smaller values of $h$.
We thus find that at large fields $h\geq h_{c}=2$ a fractional quarter spin is
likely to be accumulated at the edge in the $A_{2}$ phase. The situation at
hand is similar to what happens in topological one dimensional gapless Spin
Triplet Superconductors (STS) where there also a fractional spin-$1/4$ is
getting localized at the edge. However in the present case we do not expect
$S^{z}_{L}$ to be a sharp quantum observable in contrast with the STS case
where eigenvalues of $S^{z}_{L}$ label the different edge states of the
system. The reason for this stems from the absence of a gap in the bulk of the
Heisenberg chain in contrast with STS. The best way to check this is to
compute the variance of the operator $S^{z}_{L}$ [Eq. (49)]

$\displaystyle\sigma^{2}_{S}$ $\displaystyle=$
$\displaystyle\lim_{\alpha\rightarrow
0}\lim_{N\rightarrow\infty}\sigma^{2}_{S}(\alpha,N),$
$\displaystyle\sigma^{2}_{S}(\alpha,N)$ $\displaystyle=$
$\displaystyle\langle(S^{z}_{L}(\alpha,N))^{2}\rangle-\langle
S^{z}_{L}(\alpha,N)\rangle^{2}.$ (52)

We show in Fig. 10(a) $\sigma^{2}_{S}(\alpha,N)$ for different system sizes
and $h=4$. We observe the $\sigma^{2}_{S}(\alpha,N)$ converges to
$\sigma^{2}_{S}(\alpha)$ as system size increases. We are thus led to
conjecture the leading finite-size correction to the infinite $N$ limit:

$\displaystyle\sigma^{2}_{S}(\alpha)=\sigma^{2}_{S}(\alpha,N)+\frac{a}{(N\alpha)^{2}+b}+\cdots.$
(53)

This conjecture is based on the empirical observation that the difference of
two curves in Fig. 10(a) follows $1/\alpha^{2}$ for large $\alpha$ and remains
finite at $\alpha=0$, which is qualitatively reminiscent of how the connected
spin-spin correlation function (that the variance is related to) vanishes for
large momentum (represented by $\alpha)$. Fig. 10(b) shows
$\Delta\sigma^{2}_{S}(\alpha,N)\equiv\sigma^{2}_{S}(\alpha,N)-\sigma^{2}_{S}(\alpha,300)$,
together with the fitted equations using Eq. (53). Importantly, we find that
$a$ and $b$ are essentially independent of $N$ (with discrepancies within
$1.61\%$) demonstrating the quality of the conjectured functional form with
only two fit parameters.

We again take the $\alpha\rightarrow 0$ limit of Eq. (53) and obtain
$\sigma^{2}_{S}$. Using the fitted parameters $a$ and $b$ from
$\sigma^{2}_{S}(\alpha,600)-\sigma^{2}_{S}(\alpha,300)$, we get
$\sigma^{2}_{S}=\sigma^{2}_{S}(0,N)+a/b$. Since $\sigma^{2}_{S}$ and $a/b$
does not depend on $N$, the remaining term $\sigma^{2}_{S}(0,N)$ should also
be $N$ independent. We indeed find that $\sigma^{2}_{S}(0,N)$ is nearly zero,
and three orders of magnitude smaller than $a/b$. From this, we plot
$\sigma^{2}_{S}(\approx a/b)$ as a function of $h$ in Fig. 10(c). Although the
variance decrease as $h$ increase and crosses the phase transition, it remains
nonzero. This means that $S^{z}_{L}$ does not represent a sharp quantum
observable, and as such, we cannot, in contrast with topological STS
superconductor; label the states according to the value of the spin
accumulations at the edges.

Figure 10: (a) The variance of $S^{z}_{L}(\alpha,N)$ for various system sizes
when $h_{L}=-h_{R}=4$. (b) The difference between the variance for system size
$N$ and that of the system of $N=300$. The linear behavior for large $\alpha$
in this log-log plot suggests the difference is asymptotically
$\sim\alpha^{-2}$. The solid lines are the result of a fit based on Eq. (53).
The discrepancies at large $\alpha$ are exaggerated due to the log scale. The
inset is the same figure in linear scale which shows good agreement between
data and fit. (c) The variance obtained by the fit parameters of (b). Variance
remains nonzero beyond the phase transition. Inset is the same quantity
plotted in a $y$-axis range down to zero.

Overall we find that the ground state properties, as probed by the
magnetization profile, displays interesting and non trivial behavior as a
function of the edge fields $h$. In particular we observe that above the
critical value $h>h_{c}=2$ a fractional spin $1/4$ is likely to be accumulated
at the edge of the chain. However, in contrast with what happens in
topological spin triplet superconductor, this fractional spin is not a sharp
quantum observable and cannot be used as a genuine quantum number. The reason
for this stems from the absence of a gap in the bulk. Nevertheless, the
behavior of the spin accumulation as well as the change of the magnetization
profile in the ground state might reflect, at least qualitatively, the
eigenstate phase transition between phases $A$ and $C$ discussed in the
previous section.

##  V Discussion

We considered spin $\frac{1}{2}$ Heisenberg chain with boundary magnetic
fields and analyzed it analytically using Bethe ansatz and also numerically
using DMRG. Although the Heisenberg chain has been immensely studied and is
very well understood, the results that we have presented in this work have not
been found before. We find that the system exhibits four different ground
states for four different orientations of the boundary fields. The total spin
$S^{z}$ in each ground state may differ, which depends on the orientation of
the boundary fields and also depends on the evenness or oddness of the number
of sites of the chain. As the orientation of the magnetic fields is changed,
the system undergoes a phase transition where the ground state of the system
changes. The nature of this phase transition is currently unknown to us and
will be analyzed in the future work.

For a given orientation of the boundary fields, the system exhibits a high
energy bound state exponentially localized at an edge when the boundary
magnetic field takes values greater than the critical field $h_{c}$. Every
phase corresponding to a certain ground state can be further divided into four
sub-phases. In one of the sub-phases the system exhibits bound states at both
edges, and in one sub-phase the system exhibits no bound states while in the
remaining two sub-phases the system exhibits a bound state at the left or the
right edges.

Starting from the ground state, one can build up excitations in the bulk by
adding spinons, strings, quartets etc., and one obtains a tower of excited
states. Similarly, starting from either the state which contains one bound
state at the left or the right edge or from the state which contains two bound
states, one can build up excitations in the bulk and one obtains different
towers of excited states. Hence in the region where the system exhibits two
bound states, the Hilbert space is comprised of four towers and in the regions
where the system exhibits one bound state, the Hilbert space is comprised of
two towers and in the regions where there exists no bound states, the Hilbert
space is comprised of a single tower of excited states.

For a particular orientation of the boundary fields where the system exhibits
a certain ground state, as the values of the magnetic fields is changed, the
system undergoes an eigenstate phase transition, where the system may gain or
lose a bound state at a particular edge which results in the change in the
number of towers in the Hilbert space. Across this phase transition line where
the structure of the Hilbert space changes, the total spin $S^{z}$ of the
ground state of the system remains unchanged. To analyze the properties of the
ground state across this phase transition, we chose the $\mathbb{Z}_{2}$
symmetric point where the values of the magnetic fields at the edges take
equal values. By using DMRG we obtained the edge magnetization profile in the
regions where the magnetic fields take values greater than and lesser than the
critical field $h_{c}$. We find that when both the magnetic fields take values
greater than the critical field, the total spin accumulation at each edge
saturates to $\frac{1}{4}$. To check whether this fractional spin is a genuine
quantum observable, we calculated the variance and found that although it
saturates to a small value, it remains non zero, indicating that the
fractional spin $\frac{1}{4}$ is not a genuine quantum number. Nevertheless,
in the region where both the magnetic fields take values greater than the
critical field, there exists non zero probability to observe a non zero spin
$S^{z}$ close to each edge.

Although there exist no genuine spin fractionalization in the ground state,
the structure of the Hilbert space and the eigenstate phase transition the
system exhibits are also found in the gapless superconductors which exhibit
SPT. Even though the Hilbert space is comprised of towers of excited states,
unlike the gapped regime of the XXZ spin $\frac{1}{2}$ chain which exhibits
SSB, degenerate pairing in the spectrum is not apparent due to the gapless
nature of the bulk excitations. Although the model we considered is
integrable, the structure of its Hilbert space and the eigenstate phase
transition it exhibits might provide insight into systems with disorder which
exhibit phenomena such as many body localizaton.

###### Acknowledgements.

J.H.P. thanks Wolfgang Ketterle for useful discussions about realizing
boundary fields in experiment. P.A and P.R.P thank F. H. L. Essler for very
helpful discussions. J.L. and J.H.P. are partially supported by the Air Force
Office of Scientific Research under Grant No. FA9550-20-1-0136 and the Alfred
P. Sloan Foundation through a Sloan Research Fellowship. J.H.P. acknowledges
the Aspen Center for Physics, where some of this work was discussed, which is
supported by National Science Foundation grant PHY-1607611. P.R.P.
acknowledges support from Rutgers HEERF Fellowship during the stay at Rutgers
University, where most part of the work was done.

## References

  * Bethe (1931) H. Bethe, Zeitschrift fÃŒr Physik 71, 205 (1931).
  * Sklyanin (1988) E. K. Sklyanin, Journal of Physics A Mathematical General 21, 2375 (1988).
  * Yang and Yang (1966a) C. N. Yang and C. P. Yang, Phys. Rev. 150, 321 (1966a).
  * Pal and Huse (2010) A. Pal and D. A. Huse, Phys. Rev. B 82, 174411 (2010).
  * Abanin  _et al._ (2019) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Ilievski  _et al._ (2018) E. Ilievski, J. De Nardis, M. Medenjak, and T. c. v. Prosen, Phys. Rev. Lett. 121, 230602 (2018).
  * Ilievski  _et al._ (2015) E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler, and T. Prosen, Phys. Rev. Lett. 115, 157201 (2015).
  * Dupont and Moore (2020) M. Dupont and J. E. Moore, Phys. Rev. B 101, 121106 (2020).
  * Steiner  _et al._ (1976) M. Steiner, J. Villain, and C. Windsor, Advances in Physics 25, 87 (1976), https://doi.org/10.1080/00018737600101372 .
  * Satija  _et al._ (1980) S. K. Satija, J. D. Axe, G. Shirane, H. Yoshizawa, and K. Hirakawa, Phys. Rev. B 21, 2001 (1980).
  * Nagler  _et al._ (1991) S. E. Nagler, D. A. Tennant, R. A. Cowley, T. G. Perring, and S. K. Satija, Phys. Rev. B 44, 12361 (1991).
  * Tennant  _et al._ (1995) D. A. Tennant, R. A. Cowley, S. E. Nagler, and A. M. Tsvelik, Phys. Rev. B 52, 13368 (1995).
  * Maillet (2007) J.-M. Maillet, Prog. Math. Phys. 53, 161 (2007).
  * Jepsen  _et al._ (2021) P. N. Jepsen, W. W. Ho, J. Amato-Grill, I. Dimitrova, E. Demler, and W. Ketterle, Phys. Rev. X 11, 041054 (2021).
  * Jepsen  _et al._ (2022) P. N. Jepsen, Y. K. Lee, H. Lin, I. Dimitrova, Y. Margalit, W. W. Ho, and W. Ketterle, Nature Physics , 1 (2022).
  * Yang and Yang (1966b) C. N. Yang and C. P. Yang, Phys. Rev. 150, 327 (1966b).
  * Yang and Yang (1966c) C. N. Yang and C. P. Yang, Phys. Rev. 151, 258 (1966c).
  * Babelon  _et al._ (1983) O. Babelon, H. de Vega, and C. Viallet, Nuclear Physics B 220, 13 (1983).
  * Boos and Korepin (2001) H. E. Boos and V. E. Korepin, Journal of Physics A: Mathematical and General 34, 5311 (2001).
  * Kitanine  _et al._ (2000) N. Kitanine, J. Maillet, and V. Terras, Nuclear Physics B 567, 554 (2000).
  * Shiroishi and Takahashi (2005) M. Shiroishi and M. Takahashi, Journal of the Physical Society of Japan 74, 47 (2005), https://doi.org/10.1143/JPSJS.74S.47 .
  * Takahashi (1999) M. Takahashi,  _Thermodynamics of One-Dimensional Solvable Models, by Minoru Takahashi_ (Cambridge University Press, 1999, Tokyo, 1999).
  * Cherednik (1984) I. V. Cherednik, Theoretical and Mathematical Physics 61, 977 (1984).
  * Grisaru  _et al._ (1995) M. T. Grisaru, L. Mezincescu, and R. I. Nepomechie, Journal of Physics A: Mathematical and General 28, 1027 (1995).
  * Cao  _et al._ (2002) J. Cao, H.-Q. Lin, K.-j. Shi, and Y. Wang, arXiv e-prints , cond-mat/0212163 (2002), arXiv:cond-mat/0212163 [cond-mat.str-el] .
  * Qiao  _et al._ (2021) Y. Qiao, J. Cao, W.-L. Yang, K. Shi, and Y. Wang, Phys. Rev. B 103, L220401 (2021).
  * Sun  _et al._ (2019) P. Sun, Z.-R. Xin, Y. Qiao, K. Hao, L. Cao, J. Cao, T. Yang, and W.-L. Yang, Journal of Physics A: Mathematical and Theoretical 52, 265201 (2019).
  * Wang  _et al._ (2015) Y. Wang, W.-L. Yang, J. Cao, and K. Shi,  _Off-diagonal Bethe ansatz for exactly solvable models_ (Springer, Berlin, 2015).
  * Wang (1997) Y. Wang, Phys. Rev. B 56, 14045 (1997).
  * Skorik and Saleur (1995) S. Skorik and H. Saleur, Journal of Physics A: Mathematical and General 28, 6605 (1995).
  * Kapustin and Skorik (1996) A. Kapustin and S. Skorik, Journal of Physics A: Mathematical and General 29, 1629 (1996).
  * Keselman and Berg (2015) A. Keselman and E. Berg, Phys. Rev. B 91, 235309 (2015).
  * Keselman  _et al._ (2018) A. Keselman, E. Berg, and P. Azaria, Phys. Rev. B 98, 214501 (2018).
  * Pasnoori  _et al._ (2020) P. R. Pasnoori, N. Andrei, and P. Azaria, Phys. Rev. B 102, 214511 (2020).
  * Pasnoori  _et al._ (2021) P. R. Pasnoori, N. Andrei, and P. Azaria, Phys. Rev. B 104, 134519 (2021).
  * Fendley (2016) P. Fendley, Journal of Physics A: Mathematical and Theoretical 49, 30LT01 (2016).
  * Destri and Lowenstein” (1982) C. Destri and J. Lowenstein”, ”Nuclear Physics B” ”205”, ”369 (”1982”).
  * Fishman  _et al._ (2020) M. Fishman, S. R. White, and E. M. Stoudenmire, “The ITensor software library for tensor network calculations,” (2020), arXiv:2007.14822 .
  * Jackiw  _et al._ (1983) R. Jackiw, A. Kerman, I. Klebanov, and G. Semenoff, Nuclear Physics B 225, 233 (1983).
  * Alcaraz  _et al._ (1987) F. C. Alcaraz, M. N. Barber, M. T. Batchelor, R. J. Baxter, and G. R. W. Quispel, Journal of Physics A: Mathematical and General 20, 6397 (1987).

##  Appendix A Solution of Bethe ansatz equations

In this section we provide a detailed calculation of the ground state and the
boundary excitations for both odd and even number of sites chain in the phases
$A_{1}$, $B_{1}$ and $C_{1}$, and describe how the solution in all the other
phases can be constructed using the solution obtained in these phases.

###  A.1 $A$ phases

Consider the phase $A_{1}$. In this phase both the magnetic fields point in
the positive $z$-direction and take values $h_{L},h_{R}>h_{C}$, ($h_{C}=2$),
which translates to the boundary parameters taking values
$0<p_{L},p_{R}<\frac{1}{2}$. The Bethe equations corresponding to Bethe
reference state with all down spins are Sklyanin (1988); Wang  _et al._ (2015)

$\displaystyle\left(\frac{\lambda_{j}-\frac{i}{2}}{\lambda_{j}+\frac{i}{2}}\right)^{2N}\left(\frac{\lambda_{j}+i(\frac{1}{2}+p_{L})}{\lambda_{j}-i(\frac{1}{2}+p_{L})}\right)\left(\frac{\lambda_{j}+i(\frac{1}{2}+p_{R})}{\lambda_{j}-i(\frac{1}{2}+p_{R})}\right)=\prod_{j\neq
k=1}^{M}\left(\frac{\lambda_{j}-\lambda_{k}-i}{\lambda_{j}-\lambda_{k}+i}\right)\left(\frac{\lambda_{j}+\lambda_{k}-i}{\lambda_{j}+\lambda_{k}+i}\right)$

The Eigenvalues of the Hamiltonian are given by

$\displaystyle
E=-\sum_{j=1}^{M}\frac{2}{\lambda_{j}^{2}+\frac{1}{4}}+N-1-\frac{1}{p_{L}}-\frac{1}{p_{R}}$
(54)

####  A.1.1 $A_{1}$: Odd number of sites

Let us first consider a state with all real Bethe roots. By applying logarithm
to A.1 we obtain

$\displaystyle(2N+1)\Theta(2\lambda_{j})-\Theta\left(\frac{\lambda_{j}}{\frac{1}{2}+p_{L}}\right)-\Theta\left(\frac{\lambda_{j}}{\frac{1}{2}+p_{R}}\right)=\sum_{k\neq
j,\sigma=\pm}\Theta(\lambda_{j}+\sigma\lambda_{k})+\pi\nu(\lambda_{j})$ (55)

where $\Theta(x)=\text{ArcTan(x)}$. We introduced the counting function
$\nu(\lambda)$, where $\nu(\lambda_{j})=I_{j}$. Differentiating 55 and noting
that $2\rho(\lambda)=\frac{d}{d\lambda}\nu(\lambda)$, we obtain

$\displaystyle(2N+1)a_{\frac{1}{2}}(\lambda)-a_{\frac{1}{2}+p_{L}}(\lambda)-a_{\frac{1}{2}+p_{R}}(\lambda)=\pi\delta(\lambda)+2\pi\rho_{|-\frac{1}{2}\rangle}(\lambda)+2\int\rho_{|-\frac{1}{2}\rangle}(\mu)a_{1}(\lambda-\mu)d\mu$
(56)

where a hole at $\lambda=0$ is added as $\lambda=0$ is a trivial solution to
the Bethe equations A.1 and leads to a zero wavefunction Wang  _et al._
(2015). The reason for the subscripts for the density distribution will become
evident soon.

Taking Fourier transform we obtain

$\displaystyle\tilde{\rho}_{|-\frac{1}{2}\rangle}(\omega)=\frac{(2N+1)e^{-\frac{|\omega|}{2}}-e^{-(\frac{1}{2}+p_{L})|\omega|}-e^{-(\frac{1}{2}+p_{R})|\omega|}-1}{2\left(1+e^{-|\omega|}\right)}$
(57)

The total number of Bethe roots is
$M_{|-\frac{1}{2}\rangle}=\tilde{\rho}_{|-\frac{1}{2}\rangle}(0)=\frac{N-1}{2}$,
hence $S^{z}=-\frac{1}{2}$. We represent this state by $|-\frac{1}{2}\rangle$.
Notice that the number of roots is an integer only for a spin chain with odd
number of sites.

Using 54, one can find the energy of the states. Equation 54 can be written as

$\displaystyle E=-4\int
a_{\frac{1}{2}}(\lambda)\rho(\lambda)+N-1-\frac{1}{p_{L}}-\frac{1}{p_{R}}$
(58)

Using the density distribution 57 in the above equation we obtain

$\displaystyle
E_{|-\frac{1}{2}\rangle}=E_{0}=-(2N+1)\ln(4)+N-1+\pi+\sum_{i=L,R}\Psi\left(\frac{p_{i}}{2}\right)-\Psi\left(\frac{p_{i}-1}{2}\right)-\frac{1}{p_{i}}$
(59)

Where $\Psi$ is the digamma function. There exists two boundary string
solutions $\lambda_{p_{L}^{\prime}}=\pm
i(\frac{1}{2}+p_{L}),\lambda_{p_{R}}=\pm i(\frac{1}{2}+p_{R})$ to the Bethe
equations corresponding to all spin down Bethe reference state A.1. Adding
$\lambda_{p_{L}}$ we have

$\displaystyle(2N+1)a_{\frac{1}{2}}(\lambda)-a_{\frac{1}{2}+p_{L}}(\lambda)-a_{\frac{1}{2}+p_{R}}(\lambda)-a_{\frac{3}{2}+p_{L}}(\lambda)-a_{\frac{1}{2}-p_{L}}(\lambda)=\pi\delta(\lambda)+2\pi\rho_{|0\rangle_{p_{L}^{\prime}}}(\lambda)+2\int\rho_{|0\rangle_{p_{L}^{\prime}}}(\mu)a_{1}(\lambda-\mu)d\mu$
(60)

Taking Fourier transform we obtain,

$\displaystyle\tilde{\rho}_{|0\rangle_{p_{L}^{\prime}}}(\omega)=\tilde{\rho}_{|-\frac{1}{2}\rangle}(\omega)+\Delta\tilde{\rho}_{p_{L}^{\prime}}(\omega),\;\;\;\Delta\tilde{\rho}_{p_{L}^{\prime}}(\omega)=-\frac{e^{(\frac{3}{2}+p_{L})|\omega|}+e^{(\frac{1}{2}-p_{L})|\omega|}}{2(1+e^{-|\omega|})}$
(61)

The number of real roots is given by
$M_{|0\rangle{p_{L}^{\prime}}}-1=\tilde{\rho}_{|0\rangle_{p_{L}^{\prime}}}(0)$.
From this we obtain the total number of roots
$M_{|0\rangle{p_{L}^{\prime}}}=\frac{N}{2}$, hence $S^{z}=0$. We observe that
the number of roots is an integer only if the number of sites is even. Since
we have a chain with odd number of sites, in order for one to add a boundary
string to the state $|-\frac{1}{2}\rangle$, a propagating hole (spinon) needs
to be added as well. Adding a spinon with rapidity $\theta$, to the state
$|-\frac{1}{2}\rangle$ in addition to the boundary string
$\lambda_{p_{L}^{\prime}}$ we have

$\displaystyle(2N+1)a_{\frac{1}{2}}(\lambda)-a_{\frac{1}{2}+p_{L}}(\lambda)-a_{\frac{1}{2}+p_{R}}(\lambda)-a_{\frac{3}{2}+p_{L}}(\lambda)-a_{\frac{1}{2}-p_{L}}(\lambda)$
(62)
$\displaystyle=\pi\delta(\lambda)+\pi\delta(\lambda-\theta)+\pi\delta(\lambda+\theta)+2\pi\rho_{|-\frac{1}{2}\rangle_{(\theta,p_{L}^{\prime})}}(\lambda)+2\int\rho_{|-\frac{1}{2}\rangle_{(\theta,p_{L}^{\prime})}}(\mu)a_{1}(\lambda-\mu)d\mu$
(63)

Taking Fourier transform we obtain

$\displaystyle\tilde{\rho}_{|-\frac{1}{2}\rangle_{(\theta,p_{L}^{\prime})}}(\omega)=\tilde{\rho}_{|-\frac{1}{2}\rangle}(\omega)+\Delta\tilde{\rho}_{p_{L}^{\prime}}(\omega)+\Delta\tilde{\rho}_{\theta}(\omega),\;\;\;\;\Delta\tilde{\rho}_{\theta}(\omega)=-\frac{\cos(\theta\omega)}{(1+e^{-|\omega|})}$
(64)

The number of real roots is given by
$M_{|0\rangle{(\theta,p_{L}^{\prime})}}-1=\tilde{\rho}_{|0\rangle_{(\theta,p_{L}^{\prime})}}(0)$.
From this we find that the total number of Bethe roots is
$M_{|-\frac{1}{2}\rangle_{(\theta,p_{L}^{\prime})}}=\frac{N-1}{2}$, hence
$S^{z}=-\frac{1}{2}$. We represent this state by
$|-\frac{1}{2}\rangle_{(\theta,L)}$. We can calculate the energy of this state
using 58, we have

$\displaystyle E_{|-\frac{1}{2}\rangle_{(\theta,p_{L}^{\prime})}}=-4\int
a_{\frac{1}{2}}(\lambda)\rho_{|-\frac{1}{2}\rangle_{(\theta,L)}}(\lambda)+N-1-\frac{1}{p_{L}}-\frac{1}{p_{R}}-\frac{2}{\frac{1}{4}+(i(\frac{1}{2}+p_{L}))^{2}}$
(65)

Using 67, in the above equation we obtain

$\displaystyle
E_{|-\frac{1}{2}\rangle_{(\theta,L)}}=E_{0}+\frac{2\pi}{\sin(\pi
p_{L})}+\frac{2\pi}{\cosh(\theta)}.$ (66)

The first term is just the energy of the ground state $|-\frac{1}{2}\rangle$.
The second term is the energy of the bound state at the left edge and the
third term is the energy of the spinon with rapidity $\theta$. The energy of
both the spinon and the bound state are strictly positive with

$0<E_{\theta}<2\pi,\;\;\;\;E_{\text{bound state}}>2\pi$

Similarly we can add the boundary string corresponding to the right boundary
$\lambda_{p_{R}^{\prime}}$ along with a spinon and obtain the state
$|-\frac{1}{2}\rangle_{(\theta,R)}$ with total spin $S^{z}=-\frac{1}{2}$
described by the following density distribution

$\displaystyle\tilde{\rho}_{|-\frac{1}{2}\rangle_{(\theta,p_{R}^{\prime})}}(\omega)=\tilde{\rho}_{|-\frac{1}{2}\rangle}(\omega)+\Delta\tilde{\rho}_{p_{R}^{\prime}}(\omega)+\Delta\tilde{\rho}_{\theta}(\omega),\;\;\;\;\Delta\tilde{\rho}_{\theta}(\omega)=-\frac{\cos(\theta\omega)}{(1+e^{-|\omega|})}$
(67)

with energy

$\displaystyle
E_{|-\frac{1}{2}\rangle_{(\theta,R)}}=E_{0}+\frac{2\pi}{\sin(\pi
p_{R})}+\frac{2\pi}{\cosh(\theta)}.$ (68)

Now consider the Bethe equations corresponding to all spin up reference state
which can be obtained by taking $p_{L}\rightarrow-p_{L}$, $p_{R}\rightarrow-
p_{R}$ Sklyanin (1988) in the equations A.1, 54.

$\displaystyle\left(\frac{\lambda_{j}-i/2}{\lambda_{j}+i/2}\right)^{2N}\left(\frac{\lambda_{j}+i(\frac{1}{2}-p_{L})}{\lambda_{j}-i(\frac{1}{2}-p_{L})}\right)\left(\frac{\lambda_{j}+i(\frac{1}{2}-p_{R})}{\lambda_{j}-i(\frac{1}{2}-p_{R})}\right)=\prod_{j\neq
k=1}^{M}\left(\frac{\lambda_{j}-\lambda_{k}-i}{\lambda_{j}-\lambda_{k}+i}\right)\left(\frac{\lambda_{j}+\lambda_{k}-i}{\lambda_{j}+\lambda_{k}+i}\right)$

The Eigenvalues of the Hamiltonian are given by

$\displaystyle
E=-\sum_{j=1}^{M}\frac{2}{\lambda_{j}^{2}+\frac{1}{4}}+N-1+\frac{1}{p_{L}}+\frac{1}{p_{R}}$
(69)

By applying logarithm and following the same procedure as above, we obtain the
following distribution for the state with all real Bethe roots.

$\displaystyle\tilde{\rho}_{|\frac{1}{2}\rangle}(\omega)=\frac{(2N+1)e^{-\frac{|\omega|}{2}}-e^{-(\frac{1}{2}-p_{L})|\omega|}-e^{-(\frac{1}{2}-p_{R})|\omega|}-1}{2\left(1+e^{-|\omega|}\right)}$
(70)

The total number of roots is given by
$M_{|\frac{1}{2}\rangle}=\tilde{\rho}_{|\frac{1}{2}\rangle}(0)$. Using which
we obtain $M_{|\frac{1}{2}\rangle}=\frac{N-1}{2}$. Using this we obtain
$S^{z}=\frac{1}{2}$. We represent this state by $|\frac{1}{2}\rangle$. By
using A.1.1 and 70 we obtain the following expression for the energy of the
state $|\frac{1}{2}\rangle$

$\displaystyle E_{|\frac{1}{2}\rangle}=\frac{2\pi}{\sin(\pi
p_{L})}+\frac{2\pi}{\sin(\pi p_{R})}+E_{|-\frac{1}{2}\rangle}$ (71)

Hence we find that this state contains bound states at the left and right
edges, hence we represent this state by $|\frac{1}{2}\rangle_{L,R}$.

There exists states $|\frac{1}{2}\rangle_{(\theta,L)}$,
$|\frac{1}{2}\rangle_{(\theta,R)}$ that are degenerate (in thermodynamic
limit) to the states $|-\frac{1}{2}\rangle_{(\theta,L)}$,
$|-\frac{1}{2}\rangle_{(\theta,R)}$ respectively obtained above.

The state $|\frac{1}{2}\rangle_{(\theta,L)}$ contains a bound state at the
left edge and a spinon. This state is obtained by adding the boundary string
$\lambda_{p_{R}}=\pm i(\frac{1}{2}-p_{R})$ which is a solution to A.1.1 and a
spinon with rapidity $\theta$. Following the same procedure as above we obtain
the following density distribution

$\displaystyle\tilde{\rho}_{|\frac{1}{2}\rangle_{(\theta,p_{R})}}(\omega)=\tilde{\rho}_{|\frac{1}{2}\rangle}(\omega)+\Delta\tilde{\rho}_{p_{R}}(\omega)+\Delta\tilde{\rho}_{\theta}(\omega),\;\;\;\;\Delta\tilde{\rho}_{p_{R}}(\omega)=-\frac{e^{(\frac{3}{2}-p_{R})|\omega|}+e^{(\frac{1}{2}+p_{R})|\omega|}}{2(1+e^{-|\omega|})}$
(72)

The number of roots is given by
$M_{|\frac{1}{2}\rangle_{(\theta,p_{R})}}=\tilde{\rho}_{|\frac{1}{2}\rangle_{(\theta,p_{R})}}(0)$.
Using which we obtain $M_{|\frac{1}{2}\rangle_{(\theta,p_{R}}}=\frac{N-1}{2}$.
Using this we obtain $S^{z}=\frac{1}{2}$. The energy of this state can be
obtained following the same procedure as above, we obtain

$\displaystyle
E_{|\frac{1}{2}\rangle_{(\theta,L)}}=E_{|\frac{1}{2}\rangle}-\frac{2\pi}{\sin(\pi
p_{R})}+\frac{2\pi}{\cosh(\theta)}\equiv
E_{-|\frac{1}{2}\rangle}+\frac{2\pi}{\sin(\pi
p_{L})}+\frac{2\pi}{\cosh(\theta)}.$ (73)

Hence it contains a bound state at the left edge and is degenerate with the
state $|-\frac{1}{2}\rangle_{(\theta,L)}$ obtained previously. We represent
this state by $|\frac{1}{2}\rangle_{(\theta,L)}$. Similarly, by adding the
boundary string $\lambda_{p_{L}}=\pm i(\frac{1}{2}-p_{L})$ and a spinon to the
state $|\frac{1}{2}\rangle$, we obtain the state
$|\frac{1}{2}\rangle_{(\theta,R)}$.

Note that the bound state and the spinon both carry spin $\frac{1}{2}$. When a
bound state at either the left or the right edge is added to the state which
has spin $S^{z}=-\frac{1}{2}$, the bound state’s spin is oriented in the
positive $z$ direction. Now when a spinon is added it’s spin can be oriented
along or opposite to that of the bound state, and hence the final resulting
state $|\pm\frac{1}{2}\rangle_{(\theta,j)}$, where $j=L,R$, has total spin
$S^{z}=\pm\frac{1}{2}$ depending on the spin orientation of the spinon.

From the above analysis we see that for the odd chain, $|-\frac{1}{2}\rangle$
is the ground state. One can add one bound state at either the left edge or
the right edge accompanied by a spinon and one obtains the states
$|\pm\frac{1}{2}\rangle_{(\theta,j)}$, where $j=L,R$. On can also add two
bound states one at each edge and one obtains the state
$|\frac{1}{2}\rangle_{L,R}$. Starting from either of these six states, one can
build up excitations in the bulk by adding even number of spinons and other
type of Bethe roots such as 2-strings and quartets Destri and Lowenstein”
(1982).

Here we summarize the construction of the ground state and the excited states
in the phase $A_{1}$ for odd number of sites. In the region $\it{A}_{1}$
magnetic fields at both the boundaries point in the positive $z$ direction.
The ground state has total spin $S^{z}=-\frac{1}{2}$ with energy $E_{0}$,
exact expression of which is given by 59 and is represented by
$|-\frac{1}{2}\rangle$. It contains $\frac{N-1}{2}$ all real Bethe roots and
is constructed by starting with all spin down reference state. There exists a
state $|\frac{1}{2}\rangle_{L,R}$ with total spin $S^{z}=\frac{1}{2}$ which
contains one exponentially localized spin $S^{z}=\frac{1}{2}$ boundary bound
state at each edge. It contains $\frac{N-1}{2}$ all real Bethe roots and is
constructed by starting with reference state with all spin up. This state has
energy $E_{0}+m_{R}+m_{L}$, where $m_{i}=\frac{2\pi}{\sin(\pi p_{i})}$ is the
energy of the bound state.

There exists a state with a bound state only at the right edge, represented by
$|-\frac{1}{2}\rangle_{(\theta,R)}$. It has total spin $S^{z}=-\frac{1}{2}$
and is obtained from the state $|-\frac{1}{2}\rangle$ by adding imaginary
Bethe root $\lambda_{p_{R^{\prime}}}=\pm i(\frac{1}{2}+p_{R})$, which is
called a boundary string. One also needs to add a spinon with rapidity
$\theta$ in order to include the boundary string. There exists another state
with a bound state at the right edge represented by
$|\frac{1}{2}\rangle_{(\theta,R)}$. This state has total spin
$S^{z}=\frac{1}{2}$, and is obtained from the state
$|\frac{1}{2}\rangle_{L,R}$ by adding the boundary string $\lambda_{p_{L}}=\pm
i(\frac{1}{2}-p_{L})$ and also a spinon with rapidity $\theta$. The two states
$|\frac{1}{2}\rangle_{(\theta,R)},|-\frac{1}{2}\rangle_{(\theta,R)}$ are
degenerate in the thermodynamic limit and have energy $E_{0}+E_{\theta}+E_{R}$
but differ in the spin orientation of the spinon which has spin
$S^{z}=\pm\frac{1}{2}$ respectively. Here
$E_{\theta}=\frac{2\pi}{\cosh(\pi\theta)}$ is the energy of the spinon with
rapidity $\theta$.

There also exists two degenerate states $|\pm\frac{1}{2}\rangle_{(\theta,L)}$,
with energy $E_{0}+m_{L}+E_{\theta}$. The state
$|-\frac{1}{2}\rangle_{(\theta,L)}$ has total spin $S^{z}=-\frac{1}{2}$ and
can be obtained from the state $|-\frac{1}{2}\rangle$ by adding the boundary
string $\lambda_{p_{L^{\prime}}}=\pm i(\frac{1}{2}+p_{L})$ and a spinon with
rapidity $\theta$, whose spin is $S^{z}=-\frac{1}{2}$. The state
$|\frac{1}{2}\rangle_{(\theta,L)}$ has total spin $S^{z}=\frac{1}{2}$, and is
obtained from the state $|\frac{1}{2}\rangle$ by adding the boundary string
$\lambda_{p_{R}}=\pm i(\frac{1}{2}-p_{R})$ and a spinon with rapidity $\theta$
whose spin is $S^{z}=\frac{1}{2}$.

####  A.1.2 $A_{1}$: Even number of sites

Now consider the spin chain with even number of sites. As seen in the previous
section, starting with all spin down reference state and considering a state
with all real roots one obtains the ground state $|-\frac{1}{2}\rangle$. The
number of roots in this state is $M=\frac{N-1}{2}$. N has to be odd in order
for the number of roots to be an integer. To obtain the ground state for even
number of sites, one needs to consider a state with one less Bethe root
compared to the state $|-\frac{1}{2}\rangle$, that is starting with all spin
down reference state we need to include one spinon in addition to all real
Bethe roots. Following the procedure described in the previous section, we
obtain the following distribution

$\displaystyle\tilde{\rho}_{|-1\rangle_{\theta}}(\omega)=\frac{(2N+1)e^{-\frac{|\omega|}{2}}-e^{-(\frac{1}{2}+p_{L})|\omega|}-e^{-(\frac{1}{2}+p_{R})|\omega|}-1}{2\left(1+e^{-|\omega|}\right)}+\Delta\tilde{\rho}_{\theta}(\omega),\;\;\;\;\Delta\tilde{\rho}_{\theta}(\omega)=-\frac{\cos(\theta\omega)}{(1+e^{-|\omega|})}$
(74)

The total number of Bethe roots is $M_{|-1\rangle_{\theta}}=\frac{N-2}{2}$,
hence $S^{z}=-1$. The number of roots is an integer only for a spin chain with
even number of sites as desired. Note that the first term is same as the
density distribution describing the state $|-\frac{1}{2}\rangle$, with $N$ now
being even. The energy of this state can be calculated using 54, we obtain

$\displaystyle E_{|1\rangle_{\theta}}=E_{0}+\frac{2\pi}{\cosh(\theta)}.$ (75)

Hence in $A_{1}$, the lowest energy state for even number of sites chain is
parametrized by the rapidity of the spinon $\theta$. The ground state is
obtained in the limit where $\theta\rightarrow\infty$. Starting with all spin
up reference state and considering a state with all real roots and a spinon we
obtain

$\displaystyle\tilde{\rho}_{|1\rangle_{\theta}}(\omega)=\frac{(2N+1)e^{-\frac{|\omega|}{2}}-e^{-(\frac{1}{2}-p_{L})|\omega|}-e^{-(\frac{1}{2}-p_{R})|\omega|}-1}{2\left(1+e^{-|\omega|}\right)}+\Delta\tilde{\rho}_{\theta}(\omega)$
(76)

The total number of Bethe roots is $M_{|1\rangle_{\theta}}=\frac{N-2}{2}$,
hence $S^{z}=1$. The energy of this state can be calculated using 69, we
obtain

$\displaystyle E_{|1\rangle_{\theta}}=E_{0}+\frac{2\pi}{\sin(\pi
p_{L})}+\frac{2\pi}{\sin(\pi p_{R})}+\frac{2\pi}{\cosh(\theta)}.$ (77)

This state contains bound states at both the edges, and hence we represent
this state by $|1\rangle_{(\theta,L,R)}$. Starting with all spin down
reference state, consider a state with all real roots and the boundary string
$\lambda_{p_{L}^{\prime}}=\pm i(\frac{1}{2}+p_{L})$. Following the similar
procedure described in the previous section we obtain

$\displaystyle\tilde{\rho}_{|0\rangle_{p_{L}^{\prime}}}(\omega)=\frac{(2N+1)e^{-\frac{|\omega|}{2}}-e^{-(\frac{1}{2}+p_{L})|\omega|}-e^{-(\frac{1}{2}+p_{R})|\omega|}-1}{2\left(1+e^{-|\omega|}\right)}+\Delta\tilde{\rho}_{p_{L}^{\prime}}(\omega),\;\;\;\Delta\tilde{\rho}_{p_{L}^{\prime}}(\omega)=-\frac{e^{(\frac{3}{2}+p_{L})|\omega|}+e^{(\frac{1}{2}-p_{L})|\omega|}}{2(1+e^{-|\omega|})}$
(78)

The total number of roots is $M_{|0\rangle_{p_{L}^{\prime}}}=\frac{N}{2}$, and
hence $S^{z}=0$. We represent this state by $|0\rangle_{L}$. The energy of
this state can be calculated using the procedure described in the previous
section. We obtain

$\displaystyle E_{|0\rangle_{L}}=E_{0}+\frac{2\pi}{\sin(\pi p_{L})}$ (79)

Similarly we can add the boundary string $\lambda_{p_{R}^{\prime}}=\pm
i(\frac{1}{2}+p_{R})$ and obtain the state $|0\rangle_{R}$ which has energy
$E_{|0\rangle_{R}}=E_{0}+\frac{2\pi}{\sin(\pi p_{R})}$. These states can also
be obtained by starting with all spin up reference state, in which case,
$|0\rangle_{R}$, $|0\rangle_{L}$ contain all real roots and the boundary
strings $\lambda_{p_{L}}=\pm i(\frac{1}{2}-p_{L})$, $\lambda_{p_{R}}=\pm
i(\frac{1}{2}-p_{R})$ respectively.

Starting with all spin down reference state, one can add both the boundary
strings $\lambda_{p_{R}^{\prime}},\lambda_{p_{L}^{\prime}}$ to the state with
all real roots and one spinon. Following the regular procedure, we obtain the
following density distribution

$\displaystyle\tilde{\rho}_{|0\rangle_{(\theta
p_{L}^{\prime},p_{R}^{\prime})}}(\omega)=\frac{(2N+1)e^{-\frac{|\omega|}{2}}-e^{-(\frac{1}{2}+p_{L})|\omega|}-e^{-(\frac{1}{2}+p_{R})|\omega|}-1}{2\left(1+e^{-|\omega|}\right)}+\Delta\tilde{\rho}_{p_{L}^{\prime}}(\omega)+\Delta\tilde{\rho}_{p_{R}^{\prime}}(\omega)+\Delta\tilde{\rho}_{\theta}(\omega)$
(80)

The number of real roots is given by $M_{|0\rangle_{(\theta
p_{L}^{\prime},p_{R}^{\prime})}}-2=\tilde{\rho}_{|0\rangle_{(\theta
p_{L}^{\prime},p_{R}^{\prime})}}(0)$. From this we obtain that the total
number of roots is $M_{|0\rangle_{(\theta
p_{L}^{\prime},p_{R}^{\prime})}}=\frac{N}{2}$. This results in $S^{z}=0$. We
represent this state by $|0\rangle_{(\theta,L,R)}$. Calculating the energy we
find that this state is degenerate with the state $|1\rangle_{(\theta,L,R)}$.
Similarly, we can add both the boundary strings
$\lambda_{p_{R}},\lambda_{p_{L}}$ to the state $|1\rangle_{\theta}$ and obtain
the state $0\rangle_{\theta}$ which is degenerate with the state
$|-1\rangle_{\theta}$.

Here we summarize the construction of the ground state and the excited states
in the phase $A_{1}$ for even number of sites. In the phase $A_{1}$, the
ground state $|-1\rangle_{\theta}$ has total spin $S^{z}=-1$. It is
constructed by starting with all spin down reference state and contains
$\frac{N-2}{2}$ real roots and a spinon with rapidity $\theta$. In the
thermodynamic limit there exists another state $|0\rangle_{\theta}$ with total
spin $S^{z}=0$ which is degenerate with the ground state. This state is
obtained by starting with all spin up reference state and in addition to
$\frac{N-4}{2}$ real roots it contains two boundary strings $\lambda_{p_{R}}$,
$\lambda_{p_{L}}$ and a spinon with rapidity $\theta$. The spin orientation of
the spinon in the two ground states $|-1\rangle_{\theta}$,
$|0\rangle_{\theta}$ is along the negative and positive $z$ directions
respectively.

There exists a state represented by $|0\rangle_{L}$ which contains a bound
state at the left edge and has total spin $S^{z}=0$. This is constructed from
the reference state with all spin down. In addition to $\frac{N-2}{2}$ real
roots it contains one boundary string $\lambda_{p_{L^{\prime}}}$ and has
energy $E_{0}+m_{L}$. This state can also be constructed from the reference
state with all spin up and is made up of $\frac{N-2}{2}$ real roots and one
boundary string $\lambda_{p_{R}}$. There exists a state $|0\rangle_{R}$ which
contains a bound state at the right edge and has total spin $S^{z}=0$. It is
constructed from the reference state with all spin down and in addition to
$\frac{N-2}{2}$real roots it contains one boundary string
$\lambda_{p_{R^{\prime}}}$ and has energy $E_{0}+m_{R}$. This state can also
be constructed from the reference state with all spin up and in addition to
$\frac{N-2}{2}$ real roots it contains one boundary string $\lambda_{p_{L}}$.

There also exists a state $|0\rangle_{(\theta,L,R)}$ with total spin $S^{z}=0$
which contains one bound state at each edge. It is constructed by starting
with all spin down reference state and contains $\frac{N-4}{2}$ real roots and
two boundary strings $\lambda_{p_{L^{\prime}}}$, $\lambda_{p_{R^{\prime}}}$
and a spinon. The energy of the state is $E_{0}+m_{R}+m_{L}+E_{\theta}$. In
the thermodynamic limit there exists another degenerate state
$|1\rangle_{(\theta,L,R)}$ with total spin $S^{z}=1$ that contains one bound
state at each edge. This state is constructed by starting with all spin up
reference state and contains $\frac{N-2}{2}$ real roots and a spinon. The spin
orientation of the spinon in the two states $|0\rangle_{(\theta,L,R)}$ and
$|1\rangle_{(\theta,L,R)}$ is in the negative and positive $z$ directions
respectively.

####  A.1.3 Phase $A_{2}$

Construction of the states in the phase $A_{2}$ is similar to that in the
$A_{1}$ phase, hence we skip the details of the calculation and provide an
overview of how the ground state and boundary excited states are constructed
in $A_{2}$ phase. The Bethe equations for all down and up reference states are
obtained from A.1, 54 and A.1.1, 69 by making the transformation
$p_{L}\rightarrow-p_{L}$.

Consider the case of odd number of sites. The ground state is two fold
degenerate with energy $E_{0}+E_{\theta}$ and is represented by
$|\pm\frac{1}{2}\rangle_{\theta}$ with total spin $S^{z}=\pm\frac{1}{2}$
respectively. $|-\frac{1}{2}\rangle_{\theta}$ is constructed from all spin
down reference state. It contains the boundary string $\lambda_{p_{L}}$ and a
spinon in addition to $\frac{N-3}{2}$ real roots. The state
$|\frac{1}{2}\rangle_{\theta}$ is constructed by starting with all spin up
reference state. It contains the boundary string $\lambda_{p_{R}}$ and a
spinon in addition to $\frac{N-3}{2}$ real roots. The state with a bound state
at the left edge is represented by $|-\frac{1}{2}\rangle_{L}$ and has total
spin $S^{z}=-\frac{1}{2}$ and energy $E_{0}+m_{L}$. It is constructed by
starting with all spin down reference state and contains $\frac{N-1}{2}$ all
real roots. The state with a bound state at the right edge is represented by
$|\frac{1}{2}\rangle_{R}$ and has total spin $S^{z}=\frac{1}{2}$ and energy
$E_{0}+m_{R}$. It is constructed by starting with all spin up reference state
and contains $\frac{N-1}{2}$ all real roots. The state with bound states at
both the edges is two fold degenerate with energy
$E_{0}+m_{L}+m_{R}+E_{\theta}$ and is represented by
$|\pm\frac{1}{2}\rangle_{(\theta,L,R)}$ with total spin $S^{z}=\pm\frac{1}{2}$
respectively. $|\frac{1}{2}\rangle_{(\theta,L,R)}$ is constructed by starting
with all spin up reference state and contains the boundary string
$\lambda_{p_{L}^{\prime}}$ and a spinon in addition to $\frac{N-3}{2}$ real
roots. $|-\frac{1}{2}\rangle_{(\theta,L,R)}$ is constructed by starting with
all spin down reference state and contains the boundary string
$\lambda_{p_{R}^{\prime}}$ and a spinon in addition to $\frac{N-3}{2}$ real
roots.

Now consider the case of even number of sites. The ground state is represented
by $|0\rangle$ and has energy $E_{0}$ and total spin $S^{z}=0$. It is
constructed by starting with all spin up reference state and contains the
boundary string $\lambda_{p_{R}}$ in addition to $\frac{N-2}{2}$ real roots.
It can also be constructed by starting with all spin down reference state and
contains the boundary string $\lambda_{p_{L}^{\prime}}$ in addition to
$\frac{N-2}{2}$ real roots. The state which contains a bound state at the left
edge is two fold degenerate with energy $E_{0}+m_{L}+E_{\theta}$ and is
represented by $|0\rangle_{(\theta,L)}$ and $|-1\rangle_{(\theta,L)}$ with
total spin $S^{z}=0$ and $S^{z}=-1$ respectively. The state
$|0\rangle_{(\theta,L)}$ is constructed by starting with all spin up reference
state. It contains the boundary strings $\lambda_{p_{R}}$,
$\lambda_{p_{L}^{\prime}}$ and a spinon in addition to $\frac{N-4}{2}$ real
roots. The state $|-1\rangle_{(\theta,L)}$ is constructed by starting with all
spin down reference state and contains a spinon and $\frac{N-2}{2}$ real
roots. The state which contains a bound state at the right edge is two fold
degenerate with energy $E_{0}+m_{R}+E_{\theta}$ and is represented by
$|0\rangle_{(\theta,R)}$ and $|1\rangle_{(\theta,R)}$ with total spin
$S^{z}=0$ and $S^{z}=1$ respectively. The state $|0\rangle_{(\theta,R)}$ is
constructed by starting with all spin down reference state. It contains the
boundary strings $\lambda_{p_{R}^{\prime}}$, $\lambda_{p_{L}}$ and a spinon in
addition to $\frac{N-4}{2}$ real roots. The state $|1\rangle_{(\theta,R)}$ is
constructed by starting with all spin up reference state and contains a spinon
and $\frac{N-2}{2}$ real roots. The state with bound states at both the edges
is represented by $|0\rangle_{L,R}$ and has total spin $S^{z}=0$ and energy
$E_{0}+m_{L}+m_{R}$. It is constructed by starting with all spin up reference
state and contains the boundary string $\lambda_{p_{L}^{\prime}}$ in addition
to $\frac{N-2}{2}$ real roots. It can also be constructed by starting with all
spin down reference state and contains the boundary string
$\lambda_{p_{R}^{\prime}}$ in addition to $\frac{N-2}{2}$ real roots.

####  A.1.4 Phase $A_{3}$, $A_{4}$

In the phases $A_{3}$, $A_{4}$, both the boundary magnetic fields point in the
direction opposite that of the phases $A_{1}$ and $A_{2}$ respectively. Using
the property 2, we can obtain all the states in the phases $A_{3}$ and $A_{4}$
from the states obtained in the phase $A_{1}$ and $A_{2}$ respectively.

In constructing a state in the phase $A_{3}$ or $A_{4}$, we can use the
construction of the respective state in the phase $A_{1}$ or $A_{2}$
respectively, and use the following transformation:

$\displaystyle|\uparrow\uparrow...\uparrow\rangle\leftrightarrow|\downarrow\downarrow...\downarrow\rangle,\hskip
14.22636ptp_{L}\rightarrow-p_{L},\hskip 11.38109ptp_{R}\rightarrow p_{R}$ (81)

where the all spin up and all spin down reference states are interchanged and
the boundary magnetic fields change sign.

###  A.2 $B$ phases

Consider the phase $B_{1}$. In this phase both the magnetic fields point in
the positive $z$\- direction and take the values $0<h_{R}<h_{C}$ and
$h_{L}>h_{C}$, which corresponds to the values $p_{R}>\frac{1}{2}$ and
$0<p_{L}<\frac{1}{2}$.

####  A.2.1 Phase $B_{1}$

Consider the chain with odd number of sites. The Bethe equations corresponding
to all spin down reference state are given by A.1. The Eigenvalues of the
Hamiltonian are given by 54. Ground state contains all real roots and we
obtain the distribution 57 with energy $E_{0}$ and spin $S^{z}=-\frac{1}{2}$.
Now by adding the boundary string $\lambda_{p_{L}^{\prime}}$ along with a
spinon we obtain the state $|-\frac{1}{2}\rangle_{(\theta,L)}$ which contains
a bound state at the left edge. This state has total spin $S^{z}=-\frac{1}{2}$
and energy $E_{0}+m_{L}+E_{\theta}$.

Now consider all spin up reference state. The eigenvalues are given by 69
where as the Bethe equations A.1.1 take the form

$\displaystyle\left(\frac{\lambda_{j}-i/2}{\lambda_{j}+i/2}\right)^{2N}\left(\frac{\lambda_{j}+i(\frac{1}{2}-p_{L})}{\lambda_{j}-i(\frac{1}{2}-p_{L})}\right)\left(\frac{\lambda_{j}-i(p_{R}-\frac{1}{2})}{\lambda_{j}+i(p_{R}-\frac{1}{2})}\right)=\prod_{j\neq
k=1}^{M}\left(\frac{\lambda_{j}-\lambda_{k}-i}{\lambda_{j}-\lambda_{k}+i}\right)\left(\frac{\lambda_{j}+\lambda_{k}-i}{\lambda_{j}+\lambda_{k}+i}\right)$

Considering the state with all real roots and following the same procedure as
in the previous section we obtain

$\displaystyle\tilde{\rho}_{|0\rangle}(\omega)=\frac{(2N+1)e^{-\frac{|\omega|}{2}}-e^{-(\frac{1}{2}-p_{L})|\omega|}+e^{-(p_{R}-\frac{1}{2})|\omega|}-1}{2\left(1+e^{-|\omega|}\right)}$
(82)

The total number of roots given by
$M=\tilde{\rho}_{|0\rangle}(0)=\frac{N}{2}$. Hence we find that the number of
roots is not an integer for odd number of sites. Hence we need to consider the
state with one spinon along with real roots. We obtain

$\displaystyle\tilde{\rho}_{|\frac{1}{2}\rangle_{\theta}}(\omega)=\frac{(2N+1)e^{-\frac{|\omega|}{2}}-e^{-(\frac{1}{2}-p_{L})|\omega|}+e^{-(p_{R}-\frac{1}{2})|\omega|}-1}{2\left(1+e^{-|\omega|}\right)}+\Delta\tilde{\rho}_{\theta}(\omega)$
(83)

The number of roots is given by
$M=\tilde{\rho}_{|\frac{1}{2}\rangle_{\theta}}(0)=\frac{N-1}{2}$, and hence
the number of roots is an integer for odd number of sites and we obtain
$S^{z}=\frac{1}{2}$. The energy of this state can be calculated using 69.
Following the procedure described in the previous section we obtain
$E_{|\frac{1}{2}\rangle_{\theta}}=E_{0}+m_{L}+E_{\theta}$. Hence this state
contains a bound state at the left edge and is degenerate with the state
$|-\frac{1}{2}\rangle_{(\theta,L)}$ obtained above. We represent this state by
$|\frac{1}{2}\rangle_{(\theta,L)}$.

Now consider the chain with even number of sites. The Bethe equation
corresponding to all spin up reference state have the boundary string solution
$\lambda_{p_{L}}$. Consider a state with this boundary string and a spinon in
addition to real roots. By following the usual procedure we obtain the
following distribution

$\displaystyle\tilde{\rho}_{|0\rangle_{(\theta,p_{L})}}(\omega)=\frac{(2N+1)e^{-\frac{|\omega|}{2}}-e^{-(\frac{1}{2}-p_{L})|\omega|}+e^{-(p_{R}-\frac{1}{2})|\omega|}-1}{2\left(1+e^{-|\omega|}\right)}+\Delta\tilde{\rho}_{p_{L}}(\omega)+\Delta\tilde{\rho}_{\theta}(\omega)$
(84)

The number of real roots is given by
$M_{|0\rangle_{(\theta,p_{L})}}-1=\tilde{\rho}_{|0\rangle_{(\theta,p_{L})}}(0)=\frac{N-2}{2}$,
hence we obtain that the total number of roots is
$M_{|0\rangle_{(\theta,p_{L})}}=\frac{N}{2}$ and hence $S^{z}=0$. The energy
can be calculated using 69, we obtain
$E_{|0\rangle_{(\theta,p_{L})}}=E_{0}+E_{\theta}$. Starting with the spin down
reference state, consider the state with one spinon in addition to real roots.
We obtain the following distribution

$\displaystyle\tilde{\rho}_{|-1\rangle_{\theta}}(\omega)=\frac{(2N+1)e^{-\frac{|\omega|}{2}}-e^{-(\frac{1}{2}-p_{L})|\omega|}-e^{-(\frac{1}{2}-p_{R})|\omega|}-1}{2\left(1+e^{-|\omega|}\right)}+\Delta\tilde{\rho}_{\theta}(\omega)$
(85)

The number of roots is given by $M_{|-1\rangle_{\theta}}=\frac{N-2}{2}$, hence
$S^{z}=-1$. The energy of this state can be calculated using 54, and we obtain
$E_{|-1\rangle_{\theta}}=E_{0}+E_{\theta}$. Hence the ground state is two fold
degenerate with energy $E_{0}+E_{\theta}$. It is represented by
$|0\rangle_{\theta}$ and $|-1\rangle_{\theta}$ with total spin $S^{z}=0$ and
$S^{z}=-1$ respectively.

Starting with all spin up reference state, consider the state with all real
roots. We obtain the following distribution

$\displaystyle\tilde{\rho}_{|0\rangle}(\omega)=\frac{(2N+1)e^{-\frac{|\omega|}{2}}-e^{-(\frac{1}{2}-p_{L})|\omega|}+e^{-(p_{R}-\frac{1}{2})|\omega|}-1}{2\left(1+e^{-|\omega|}\right)}$
(86)

The number of roots is given by
$M_{|0\rangle}=\tilde{\rho}_{|0\rangle}(0)=\frac{N}{2}$, and hence it has
total spin $S^{z}=0$. The energy of this state can be calculated using 69, we
obtain $E_{|0\rangle}=E_{0}+m_{L}$. Hence it has a bound state at the left
edge. This state is represented by $|0\rangle_{L}$.

In summary, for the spin chain with odd number of sites, in the phases
$B_{1}$, the ground state has total spin $S^{z}=-\frac{1}{2}$ and is
represented by $|-\frac{1}{2}\rangle$. This state is constructed by starting
with all spin down reference state and contains $\frac{N-1}{2}$ real roots and
has energy $E_{0}$. There exists a state which contains a bound state at the
left edge, which is represented by $|-\frac{1}{2}\rangle_{(\theta,L)}$ and has
total spin $S^{z}=-\frac{1}{2}$ and energy $E_{0}+E_{\theta}+m_{L}$. This
state is obtained from the state $|-\frac{1}{2}\rangle$ by adding the boundary
string $\lambda_{p_{L^{\prime}}}$ and a spinon with rapidity $\theta$. There
exists a degenerate state $|\frac{1}{2}\rangle_{(\theta,L)}$ which contains a
bound state at the left edge and has total spin $S^{z}=\frac{1}{2}$. This
state is obtained by starting with reference state with all spin up. It
contains $\frac{N-1}{2}$ real roots and a spinon with rapidity $\theta$.

For the spin chain with even number of sites, in the phase $B_{1}$, the ground
state is two fold degenerate. The ground state $|-1\rangle_{\theta}$, with
total spin $S^{z}=-1$ is constructed by starting with all spin down reference
state and contains $\frac{N-2}{2}$ real roots and a spinon with rapidity
$\theta$. The ground state $|0\rangle_{\theta}$ with total spin $S^{z}=0$ is
constructed by starting with all spin up reference state. It contains a spinon
and the boundary string $\lambda_{p_{L}}$ in addition to $\frac{N-2}{2}$ real
roots. The spin orientation of the spinon in the states $|-1\rangle_{\theta}$
and $|0\rangle_{\theta}$ is along the negative and positive $z$ directions
respectively. There exists a state $|0\rangle_{L}$ with a bound state at the
left boundary. It has total spin $S^{z}=0$ and is constructed by starting with
all spin up reference state and contains $\frac{N}{2}$ real roots. It can also
be constructed by starting with all spin down reference state and it includes
the boundary string $\lambda_{p_{L^{\prime}}}$ in addition to $\frac{N-2}{2}$
real roots.

####  A.2.2 Phase $B_{2}$

The construction of the state in the phase $B_{2}$ is similar to that in the
phase $B_{1}$, hence we skip the details of the calculation and provide an
overview of how the ground state and boundary excited states are constructed
in the phase $B_{2}$.

Consider the chain with odd number of sites. The ground state is two fold
degenerate with energy $E_{0}+E_{\theta}$. It is represented by
$|\pm\frac{1}{2}\rangle_{\theta}$ with total spin $S^{z}=\pm\frac{1}{2}$. The
state $|\frac{1}{2}\rangle_{\theta}$ is obtained by starting with reference
state with all spin up and it contains a spinon in addition to $\frac{N-2}{2}$
real Bethe roots. The state $|-\frac{1}{2}\rangle_{\theta}$ is obtained by
starting with all spin down reference state and contains the boundary string
$\lambda_{p_{L}}$ and a spinon in addition to real roots. The state with the
bound state at the left edge has energy $E_{0}+m_{L}$ and total spin
$S^{z}=-\frac{1}{2}$. It is represented by $|-\frac{1}{2}\rangle_{L}$ and is
constructed by starting with all spin down reference state and contains all
real roots.

Consider the chain with even number of sites. The ground state $|0\rangle$ has
energy $E_{0}$ and total spin $S^{z}=0$. It is obtained by starting with
reference state with all spin up and it contains $\frac{N}{2}$ all real roots.
It can also be obtained by starting with all spin down reference state and it
contains the boundary string $\lambda_{p_{L}}$ in addition to
$\frac{N-2}{2}$real roots. The state with the bound state at the left edge is
doubly degenerate with energy $E_{0}+m_{L}+E_{\theta}$ and is represented by
$|-1\rangle_{(\theta,L)}$ and $|0\rangle_{(\theta,L)}$ with total spin
$S^{z}=-1$ and $S^{z}=0$ respectively. The state $|-1\rangle_{(\theta,L)}$ is
obtained by starting with all spin down reference state and contains one
spinon in addition to $\frac{N-2}{2}$ real roots. The state
$|0\rangle_{(\theta,L)}$ is obtained by starting with all spin all reference
state. It contains the boundary string $\lambda_{p_{L}^{\prime}}$ and a spinon
in addition to $\frac{N-2}{2}$ real roots.

####  A.2.3 Other $B$ phases

The states in the phases $B_{8}$ and $B_{7}$ can be obtained from the states
in the phases $B_{1}$ and $B_{2}$ respectively by the transformation
$p_{L}\leftrightarrow p_{R}$. The states in the phases $B_{5}$,$B_{6}$,
$B_{3}$ and $B_{4}$ can be obtained from the states in the phases $B_{1}$,
$B_{2}$, $B7$ and $B_{8}$ respectively by the transformation 81.

###  A.3 $C$ phases

####  A.3.1 Odd number of sites

In subregion $\it{C_{1}}$, the ground state is $|-\frac{1}{2}\rangle$ with
total spin $S^{z}=-\frac{1}{2}$. This state is constructed from the reference
state with all down spins and contains $\frac{N-1}{2}$ real roots. In
$\it{C_{3}}$ the ground state is $|\frac{1}{2}\rangle$ with total spin
$S^{z}=\frac{1}{2}$. This state is constructed from the reference state with
all up spins and contains $\frac{N-1}{2}$ real roots. In subregions
$\it{C_{2}},\it{C_{4}}$ the ground state is two fold degenerate and contains a
spinon with rapidity $\theta$. The spin orientation of the spinon dictates the
total spin $S^{z}=\pm\frac{1}{2}$ of the state. They are represented by
$|\pm\frac{1}{2}\rangle_{(\pm,\theta)}$, and contain $\frac{N-1}{2}$ real
roots and a spinon with rapidity $\theta$, and constructed from either all
spin up or down reference states and contain all real roots and a spinon with
rapidity $\theta$.

####  A.3.2 Even number of sites

In subregion $\it{C_{1}}$, the ground state is two fold degenerate in
thermodynamic limit and is represented by $|0\rangle_{\theta}$,
$|-1\rangle_{\theta}$ with total spin $S^{z}=0$, $S^{z}=-1$ respectively. The
state $|0\rangle_{\theta}$ is constructed from the reference state with all up
spin and contains $\frac{N}{2}$ real roots and a spinon with rapidity
$\theta$. the state $|-1\rangle_{\theta}$ is constructed with the reference
state with all down spins and contains $\frac{N-2}{2}$ real roots and a spinon
with rapidity $\theta$. The spin orientation of the spinon is in the negative
and positive $z$ direction in the states $|-1\rangle_{\theta}$ and
$|0\rangle_{\theta}$ respectively. In subregion $\it{C_{3}}$, the ground state
is two fold degenerate and is represented by $|0\rangle_{\theta}$,
$|1\rangle_{\theta}$ with total spin $S^{z}=0$, $S^{z}=1$ respectively.
$|0\rangle_{\theta}$ contains a spinon with rapidity $\theta$, the spin
orientation of which is in the negative $z$ direction. It is constructed from
the reference state with all spin down and contains $\frac{N}{2}$ real roots
and a spinon with rapidity $\theta$. $|1\rangle_{\theta}$ is constructed with
the reference state with all spin up and contains $\frac{N-2}{2}$ real roots
and a spinon with rapidity $\theta$ with spin oriented in the positive $z$
direction. In subregions $C_{2},C_{4}$, the ground state has total spin
$S^{z}=0$ and is represented by $|0\rangle$. It can be constructed from the
reference state with all spin up or all spin down and contains $\frac{N}{2}$
real roots.

##  Appendix B Bound state wavefunction

In this section we provide the bound state wavefunction corresponding to the
boundary string in one particle sector (one flipped spin). Let us consider the
sub-phase $A_{1}$. The Bethe equations corresponding to all spin down
reference state in one particle sector are given by

$\displaystyle\left(\frac{\lambda-\frac{i}{2}}{\lambda+\frac{i}{2}}\right)^{2N}\left(\frac{\lambda+i(\frac{1}{2}+p_{L})}{\lambda-i(\frac{1}{2}+p_{L})}\right)\left(\frac{\lambda+i(\frac{1}{2}+p_{R})}{\lambda-i(\frac{1}{2}+p_{R})}\right)=1$
(87)

The wavefunction is given by Alcaraz  _et al._ (1987); Skorik and Saleur
(1995)

$\displaystyle
f(x)=\left(\frac{\lambda+\frac{i}{2}}{\lambda-\frac{i}{2}}\right)^{N-x}\left(\frac{\lambda-i\left(\frac{1}{2}+p_{R}\right)}{p_{R}\left(\lambda-\frac{i}{2}\right)}\right)-\left(\frac{\lambda-\frac{i}{2}}{\lambda+\frac{i}{2}}\right)^{N-x}\left(\frac{\lambda+i\left(\frac{1}{2}+p_{R}\right)}{p_{R}\left(\lambda+\frac{i}{2}\right)}\right)$
(88)

When $\lambda=\pm i(\frac{1}{2}+p_{R})$, which is the boundary string
associated with the right edge, we readily obtain the wavefunction for the
boundstate localized at the right edge

$\displaystyle
f_{R}(x)=\pm\left(\frac{1+2p_{R}}{p_{R}(1+p_{R})}\right)\left(\frac{1+p_{R}}{p_{R}}\right)^{-(N-x)}$
(89)

To obtain the bound state wavefunction associated with the left edge, we
multiply the wavefunction (88) with a normalization constant

$\displaystyle\mathcal{A}=\left(\frac{p_{R}}{p_{L}}\right)\left(\frac{\lambda+i(\frac{1}{2}+p_{L})}{\lambda+i(\frac{1}{2}+p_{R})}\right)\left(\frac{\lambda-i(\frac{1}{2}+p_{L})}{\lambda-i(\frac{1}{2}+p_{R})}\right)$
(90)

and use the one particle Bethe equation (87). We obtain

$\displaystyle\mathcal{A}f(x)=\left(\frac{\lambda+\frac{i}{2}}{\lambda-\frac{i}{2}}\right)^{-x}\left(\frac{\lambda+i\left(\frac{1}{2}+p_{L}\right)}{p_{L}\left(\lambda-\frac{i}{2}\right)}\right)-\left(\frac{\lambda-\frac{i}{2}}{\lambda+\frac{i}{2}}\right)^{-x}\left(\frac{\lambda-i\left(\frac{1}{2}+p_{L}\right)}{p_{L}\left(\lambda+\frac{i}{2}\right)}\right)$
(91)

When $\lambda=\pm i(\frac{1}{2}+p_{L})$, which is the boundary string
associated with the left edge, we obtain the bound state wavefunction
localized at the left edge

$\displaystyle
f_{L}(x)=\pm\left(\frac{1+2p_{L}}{p_{L}^{2}}\right)\left(\frac{1+p_{L}}{p_{L}}\right)^{-x}$
(92)

Hence, we find that the two boundary string solutions correspond to
exponentially localized bound states $\sim e^{-\kappa_{L}x}$ and $\sim
e^{-\kappa_{R}(N-x)}$, where

$\displaystyle\kappa_{j}=\log\left(\frac{h_{j}+1}{h_{j}}\right),\;\;\;j=L,R.$
(93)

Note that when the magnetic fields at the edges take equal values
$h_{L}=h_{R}=h$ ($p_{L}=p_{R}=p$), the normalization constant $\mathcal{A}=1$.
To obtain the bound state wavefunctions in this case, consider

$\displaystyle f(x)\pm\left(\frac{p}{1+p}\right)f(x).$ (94)

In this limit we only have one boundary string solution $\lambda=\pm
i(\frac{1}{2}+p)$, corresponding to the double pole of the Bethe equations
(87). Using this in (94), we obtain

$\displaystyle
f(x)_{\pm}=-\left(\frac{1+2p}{p(1+p)}\right)\left[\left(\frac{1+p}{p}\right)^{-(N-x)}\pm\left(\frac{1+p}{p}\right)^{-x}\right]$
(95)

Hence, we find that when $h_{L}=h_{R}$, in contrast to Skorik and Saleur
(1995), there exist two bound states solutions simultaneously localized at
both the edges

$\displaystyle f(x)_{\pm}\sim(e^{-\kappa(N-x)}\pm e^{-\kappa
x}),\;\;\;\kappa=\log\left(\frac{h+1}{h}\right).$ (96)

##  Appendix C Details and validity of the fitting

Figure 11: (a) The same figure from Fig. 7(b) together with the smooth
interpolating function through the missing region. (b) Fitting parameters $A$
and $B$ (see Eq. (47)) with the substitution of $C$ and $D$ with the
interpolated function obtained in (a). (c) Magnetization data (dots) and its
fit to Eq. (47) (dashed lines) of the $h$ values of the shaded region using
the parameters obtained in (a,b) for the first 20 sites of the $N=1000$ chain.
The good agreement between the data and fit demonstrates the interpolation in
(a) is valid.

In Fig. 7(b) and Fig. 8(c) we have indicated the region where the fitting to
Eq. (47) numerically fails with a shade. Here, we show that this is merely due
to numerical overfittings, and is not of physical consequence.

First, we interpolate the data in Fig. 7(b) through the missing region. Some
data points in the figure near $h_{c}=2$ are additionally excluded in the
interpolation to make the interpolated function smooth. The results, together
with the original Fig. 7(b) are plotted in Fig. 11(a). Now we take the values
from the interpolated function for $C$ and $D$, and re-fit to Eq. (47) to
obtain $A$ and $B$. This result is shown in Fig. 11(b) and in Fig. 8(b) in the
main text. For the special point $h=2$, $C=D=0$ is used in the fit, and note
that the resulting $A$ and $B$ fit parameters are smooth near $h=2$. The $A$,
$B$ data away from $h=2$ in Fig. 11(b) is identical to that obtained from the
calculation in Fig. 7(b) and Fig. 8(c).

Finally, we use the newly fitted parameters $A$ to $D$ from Fig. 11(a,b) and
compare with the original $S^{z}(x_{i})$ data for the shaded region, which are
$h=1.90,~{}1.95,~{}2.05$, in Fig. 11(c). ($h=2.00$ data is included in Fig.
7(a)) The agreement between the data and fit is excellent in all three values
of $h$. This demonstrates that the fitting failure at the shaded region is
because of overfitting, and the true physics in those parameters connect
smoothly to the behavior outside the region. We also claim that small
deviations from the interpolated function and data (for example, $C$ near
$h=1.8$) are also the result of such numerical issues of overfit, while less
serious than that more closer to $h=2$.