An integrable spin chain with Hilbert space fragmentation and solvable real time dynamics

We revisit the so-called folded XXZ model, which was treated earlier by two independent research groups. We argue that this spin-1/2 chain is one of the simplest quantum integrable models, yet it has quite remarkable physical properties. The particles have constant scattering lengths, which leads to a simple treatment of the exact spectrum and the dynamics of the system. The Hilbert space of the model is fragmented, leading to exponentially large degeneracies in the spectrum, such that the exponent depends on the particle content of a given state. We provide an alternative derivation of the Hamiltonian and the conserved charges of the model, including a new interpretation of the so-called"dual model"considered earlier. We also construct a non-local map that connects the model with the Maassarani-Mathieu spin chain, also known as the SU(3) XX model. We consider the exact solution of the model with periodic and open boundary conditions, and also derive multiple descriptions of the exact thermodynamics of the model. We consider quantum quenches of different types. In one class of problems the dynamics can be treated relatively easily: we compute an example for the real time dependence of a local observable. In another class of quenches the degeneracies of the model lead to the breakdown of equilibration, and we argue that they can lead to persistent oscillations. We also discuss connections with the $T\bar T$- and hard rod deformations known from Quantum Field Theories.


I. INTRODUCTION
Quantum integrable models are special many body systems that allow for an exact solution, at least for certain observables and in certain situations. They possess a large number of extra conservation laws, which constrain their dynamics, eventually allowing their solvability [1]. Consequently, they have been used fruitfully in the past to compute ground-state or thermal properties in various situations, with applications ranging from condensed matter to high-energy physics [2][3][4].
More recently, new challenges have been raised by the study of the non-equilibrium dynamics for such systems. A first set of questions is concerned with the relaxation dynamics of physical observables, for instance after a quantum quench [5]. While it is now well-established that integrable models do not thermalize and instead equilibrate at late times to steady states described by the Generalized Gibbs Ensemble (VEG) [6], the exact computation of the finite time dynamics using the traditional methods of integrability has remained a very difficult task, and to this date only few results are available, associated with some very specific observables [7] or particular models (see below). Another question is that of transport properties, which can be described by the recent theory of Generalized Hydrodynamics (GHD) [8,9]. Even though GHD is largely successful, it relies on some assumptions which have not yet been proven in general [10,11]. Some elements of the theory can be considered proven, for example key statements about the mean values of current operators (see the review [12]), but it would be desirable to rigorously prove more aspects of the theory. It was shown in the remarkable work [13] that certain statements of GGE and GHD can be checked in the Lieb-Liniger model in a large coupling expansion; up to date this result is one of the most convincing analytic checks of GGE and GHD (for closely related works see [10,11,[14][15][16]). Nevertheless there remains a need for simple toy models, which have genuine interactions in them, and which can lead to exact proofs of the GHD predictions.
Therefore, the question that motivates the present work is the following: What are the simplest integrable models? By this we mean models that have genuine interactions (in contrast with those related with free theories as for example the famous XX chain), but which are nevertheless simple enough to bypass some of the difficulties presented by generic integrable models. We should note that the more conventional problem of equilibrium physics can also benefit from simple toy models, because the computation of correlation functions is a notoriously difficult problem in Bethe Ansatz solvable models (see for example [17][18][19][20][21]).
To address this question, let us discuss some simple models that already appeared in the literature.
One of the most important examples is the hard rod gas in the continuum; for the classical hard rod gas see [22][23][24][25][26][27], for the quantum case see [28][29][30]. In these models the fundamental particles have finite and fixed widths, and the speed of wave propagation depends on this width arXiv:2105.02252v1 [cond-mat.stat-mech] 5 May 2021 through the modification of the free space between the particles. Alternatively, the scattering in these models is such that it leads to a fixed displacement given by the hard rod length. In the classical model the emergence of hydrodynamics could be proven for a wide range of initial conditions in [25].
Another important example is the Rule 54 model, which is sometimes claimed to be the simplest interacting integrable model. It is a classical deterministic cellular automaton originally developed in [31], which has been an object of interest in the last couple of years (see the review [32]). There are two fundamental particles in the model: the left-and right-movers, and their scattering leads to a constant displacement of 1 lattice unit. Despite its simplicity, the model shows many distinguishing features of generic interacting many-body systems, such as relaxation at late times and coexistence of ballistic and diffusive transport [33]. The exact real time evolution in this system was computed for special initial states in the recent work [34], which quite remarkably does not use the standard methods of integrability. Instead, it is built on methods developed for "dual unitary quantum gate" models, see [35][36][37]. The dual unitary models are not integrable in the traditional sense, nevertheless they allow for exact solutions [36,37].
A less well known model which nevertheless belongs in this list is the so-called phase model, which arises as the q → ∞ limit of the q-bison lattice model. The model was first studied in [38][39][40], and later also in [41][42][43][44][45][46] with the focus being on equilibrium correlation functions and relations with dominatrix and the theory of symmetric functions. Real time dynamics in the model was first investigated in [47], where it was rigorously shown for a certain quench that the system equilibrates to the GGE prediction. This was achieved by computing the exact time dependence of a local observable and comparing its asymptotic value to the GGE average; to our best knowledge this was the first example for such an exact computation in a genuinely interacting case. Afterwards further quench problems were considered in [48], together with a special non-local connection with the XX chain. Although it is not explicitly stated in [47,48], the scattering of fundamental bosons in this model is such that the trajectories of the particles get displaced by one lattice site.
A common property of these simple models is that they can be considered as deformations of free models, although the deformation is highly non-local. In the case of the hard-rod gas the deformation in question is in the class of the famous TT -deformations, as discussed in [28,29]. For the q-boson the deformation is given by the non-local mapping to the XX chain. In the case of the Rule 54 model such an explicit deformation has not yet been worked out, but connections with the TTdeformations were already pointed out in [49].
More recently another relatively simple model was studied in [50] and later in [51,52], where it was called the "folded XXZ model"; the first appearance of the Hamilto-nian was apparently in [53], where it was obtained as an effective Hamiltonian. In [51][52][53] the model was derived by considering the large ∆ behaviour of the famous XXZ chain (for closely related works see [54,55]). The model has a four site Hamiltonian and a very special dynamics: it has a sector which is equivalent to the so-called constrained XXZ model at the free fermion point [56][57][58], where the particles have genuine interactions originating from a hard rod constraint. However, the full Hilbert space of the model is much larger due to the presence of an additional type of particle (which can be called domain wall, or DW). The DW's are not dynamical: in the absence of particles they lead to frozen configurations and exponentially degenerate energy levels; this phenomenon was interpreted as "Hilbert space fragmentation" in [50] and it was also discussed in [51]. It is also important that the DW's interact with the particles and thus they affect the dynamics of the model. A particle can be considered as a bound state of two DW's, thus becoming dynamical; this phenomenon bears some similarities with the so-called fractonic excitations [59,60].
The exact coordinate Bethe Ansatz solution of the model was presented in [51,52] and connections were found to other existing models in the literature. However, the algebraic origin of the model and its conserved charges was not clarified completely, and a number of interesting features of the model were not yet explored.
In this paper we contribute by a completely independent derivation of the "folded XXZ model" and its Bethe Ansatz solution [61]. Also, we point out new connections to existing models in the literature. We explain that the model is in a special class of systems, which have constant scattering lengths. This class includes the hard rod gas, the Rule 54 model, and the phase model mentioned above. We also discuss relations with the TT deformation, and in an independent line of computation we present an exact solution of a quantum quench, mirroring some of the computations of [47,48]. The results show that this special integrable model sits at the intersection of three very active yet seemingly distant fields: out-of-equilibrium dynamics, TT deformation and fracton systems. This makes the model truly unique and highly interesting. The exact solution of the model might shed lights on the possible exciting connections between these research areas.

II. THE MODEL
We consider a spin-1/2 chain with periodic boundary conditions. The Hamiltonian is [50][51][52] H = Q 4 + hQ 1 + µQ 2 . (1) Here Q 1 , Q 2 and Q 4 are mutually commuting extensive operators that are given below, and h and µ are to be understood as a magnetic field and a chemical potential.
The kinematical part of the Hamiltonian is This is a 4-site operator, which generates a spin exchange between neighbouring sites, controlled by the state of two further neighbours. To be precise, the exchange between sites j + 1 and j + 2 has amplitude -1/2 if the state of the sites j and j + 3 is the same, and zero amplitude if it is different. The numerical pre-factor is added for later convenience.
The remaining charges are Q 1 can be interpreted as particle number and Q 2 as a domain wall number. It can be checked that the three operators are mutually commutative. The Hamiltonian (1) first appeared in [53] as an effective Hamiltonian. Afterwards it was studied independently in [50], although that work only treated it with open boundary conditions. This case will be studied separately in Section V. Afterwards the model also appeared in [51,52], where an exact solution was given. To be precise, the papers [51,52] treated an equivalent dual model, which is discussed later in our Section VI. These papers focused on the periodic case and they did not work out the solution for h = 0. There is considerable overlap between the results of [51,52] and our work. We choose to present a complete treatment of the model from our point of view, meanwhile also explaining what is new and what was already given in [50] and/or [51,52].

A. Relation to the constrained XXZ model
There are various ways of writing the Hamiltonian, and there are multiple connections to known integrable models in the literature. One of these connections is with the constrained XXZ model, which describes a spin chain with XXZ type interaction, where particles have a finite "width" l ∈ Z + . The model was treated in a number of works [54][55][56][57][58], and its Hamiltonian is as follows. Let us choose a convention that the down spins are interpreted as particles. Then the model is given by where P l is a projector onto the states of the Hilbert space where there are at least l up spins between two down spins. In other words, the linear space selected by P l describes particles that have a "width" l + 1.
It is known that the constrained XXZ model is integrable for every l and ∆, its Bethe Ansatz solution can be found in [56,57]. Furthermore, an algebraic treatment of its integrability properties was given in [58], where it was related to a vertex model with long range interactions. This implies the existence of an infinite family of commuting conserved charges for the model.
To establish a relation with our model we write where E ± j are projection operators to the up/down spins on site j. Then we can see that The two operators Q α 4 are related by spin reflection. It is then easy to see that the projected operator P 1 Q + 4 P 1 is equivalent to (4) with l = 1 and ∆ = 0. The spin reflected version obtained from Q − 4 is equivalent to a spin reflected constrained XXZ model. Thus Q 4 can be considered the spin reflection invariant version of (4) without any constraints.
The connection between the two models was already noted in [50] and in [51,52]. The work [50] actually considered a perturbation of (1) by an interaction term, such that after the projection the model becomes equivalent to the constrained XXZ chain with finite ∆. It was shown in [50] that this full model with the interaction term is not integrable, and only the constrained sector is Bethe Ansatz solvable. However, it was not understood in [50] that for the special point of ∆ = 0 (corresponding to our Q 4 ) the full model is integrable and Bethe Ansatz solvable; this point was correctly given in [51,52].

III. DYNAMICS AND BETHE ANSATZ
Let us investigate the dynamics of the model. The hopping term Q 4 describes spin exchange between neighbouring sites, and this can be interpreted as particle hopping.
There are two ways to identify fundamental particles in the model, by choosing two different reference states. We can choose the reference state |∅ consisting of all spins up, or the state ∅ with all spins down. Single particle spin wave excitations can be constructed above either vacuum state.
Focusing on the state |∅ we can define an N -particle basis as where we apply the restriction 1 ≤ x 1 < · · · < x N ≤ L to avoid double counting.
A. Single particle states Single particle states with lattice momentum p are given simply as The associated energy is found to be E = e(p), e(p) = − cos(p).
It is useful to define a semi-classical (bare) speed, which describes the propagation of wave packets. It is given by the well known expression v(p) = de(p) dp = sin(p).
In a finite volume L the quantization condition is simply B. Two-particle scattering states Let us then focus on the scattering of two spin waves. It is clear from the Hamiltonian that as the two incoming particles approach each other, they can not occupy neighbouring sites. The reason is that any amplitude which would bring two particles to neighbouring sites is forbidden by the form of Q 4 . Nevertheless there will be an interaction between the two particles: the exact wave function for the scattering of particles with momenta p 1 and p 2 is found to be with It can be checked that this is an exact eigenstate with the energy given by This is close to a free fermionic wave function, but there are extra phases in the two terms, which can not be transformed away. In fact we can read off the two-body scattering factor The physical meaning of the scattering phase can be read off directly the wave function (13). If we identify the first term as an incoming wave and the second term as the outgoing wave (corresponding to v 1 > v 2 ), then we see that both particles suffer a displacement of ±1 sites.
This displacement is most easily understood in a semiclassical picture, by constructing wave packets and looking at their peaks before and after the scattering. We observe that the particle coming from the left (right) is moved by 1 site to the right (left). This is consistent with attractive interactions in a semi-classical or classical picture.
It is worthwhile to recall that generally the scattering displacements can be computed from the derivatives of the scattering phase δ(p 1 , p 2 ) with respect to the momenta [62]. Generally such a scattering also results in the distortion of the wave packet. In contrast, now the function δ(p 1 , p 2 ) is linear in both variables, and we obtain a complete and exact displacement of the wave packets.
Note that the scattering phase is such that the neighbouring positions x 2 = x 1 + 1 are automatically discarded, they receive zero amplitude. Therefore it is not necessary to impose the constraint by hand as in (4), in this process it emerges dynamically, and the resulting Smatrix phase is identical to the one of the constrained XXZ model with ∆ = 0 and = 1.
In a finite volume L the periodicity condition for the wave function results in the Bethe Ansatz equations Using the concrete form of the scattering phase, and denoting e ip1 e ip2 = e iP we get the equivalent set of equations We see that the total momentum is quantized as usually, but for the single particle momenta the conditions are different from a free theory. The volume appearing in the conditions is changed by −2, and there appears a twist which depends on the total momentum. These two changes signal that there is indeed true interaction in the model, even though the wave function appears almost free. C. N -particle scattering states The previous wave function can be written as a determinant It turns out that this structure can be generalized to higher particle states. The N -particle sector is found to be integrable, with elastic factorized scattering, with the two-body S-matrix given by (15). This will be proven in Section IV. As a result, the N -particle wave function is written as a Vandermonde-like determinant where C is a matrix of size N × N with elements This is the same wave function as given in [54] and related works. The energy of such a state is given by We note again that the wave function is such that neighbouring particle positions x j+1 = x j +1 are automatically forbidden, and this condition naturally emerges from the dynamics of the model. In finite volume the Bethe equations can be written as Introducing again the total momentum we get the equivalent set We see again that the quantization conditions are almost free. Every momentum p j is quantized with a simple relation where a modified volume L − N appears, together with a twist depending on the particle number and the overall momentum. Each quantization condition can then be solved almost independently, with the only requirement that the overall constraint (23) is satisfied. The appearance of a modified volume signals a relation to the TT or hard rod deformations. We discuss this connection in more detail in Section XIII.

D. Domain walls
So far we investigated scattering states consisting of separate particles. In these states neighbouring positions are forbidden. However, configurations with neighbouring down spins are not excluded from the Hilbert space, thus we need to investigate them separately.
It turns out that in this model blocks of spins with the same orientation are completely frozen if the block length is at least 2. The simplest example is probably the case of two down spins on neighbour sites, for example Direct calculation shows that this state is an eigenstate of Q 4 with eigenvalue 0. The stability of this state originates from the control on the particle jumps in Q 4 : neither down spin is allowed to hop, because the control spins result in zero amplitude. Similarly, the same eigenvalue zero is obtained for any state which has n consecutive down spins embedded in the vacuum state.
To understand the situation we introduce the concept of the Domain Wall (DW): We call a DW the boundary between two regions of the chain with different spin orientation, such that each region has at least length 2. With this definition the state (25) can be considered as having two DW's at a distance 2 from each other. In this picture a single particle can be interpreted as a bound state of two DW's, such that the DW's are at neighbouring positions. However, it is best to distinguish this situation and call the particle simply a "particle". The frozen blocks of spins can be regarded as bound states of particles.
In a finite volume the number of the domain walls is always even, but in infinite volume we can also have an odd number of them, when the asymptotic reference states on the left and the right are different.
The states of the computational basis with an arbitrary number of DW's at arbitrary positions are eigenstates of Q 4 with eigenvalue 0, given that the state does not include any particles. This results in a huge degeneracy of the reference state, which scales exponentially with the volume (see [50][51][52] and our discussion in Sec. VIII).
In order to fully understand the dynamics of the model we need to consider the scattering of particles with domain walls. The simplest situation is the case with one particle and one DW. In this situation there is an incoming particle which meets a standing domain wall. Let us situate the domain wall initially at position 2. Then the incoming wave is written as Note that the incoming particle is not allowed to hop to the site x = 1, because this process is forbidden by the control in the Hamiltonian. However, if the particle occupies x = 0, then the up spin at site x = 1 can start to propagate into the vacuum formed by the down spins on the right. Thus a hole will propagate as an outgoing wave. We can find the total wave function to be where it is now understood thatx is missing in the list, that is the position of the hole. Direct computation shows that this wave function is an exact eigenstate with energy E = e(p). The interpretation of this wave function is the following: As a result of the scattering the domain wall gets displaced by 2 sites, and the particle also obtains a displacement of 1 site forwards, see figure 1. Interestingly, the domain wall does not contribute to the energy, but it has an effect on the dynamics due to the scattering displacement. The scattering event is shown schematically in Fig. 1.
We can also compute the situation with one particle and two domain walls, which is equivalent to having a Figure 1. Schematic representation of the scattering of a particle with a domain wall. When the particle arrives at position 0, the jump to site 1 is forbidden by the control factor in (2). However, in this case the up spin at site 1 becomes mobile and it will start to propagate to the right. Eventually the domain wall gets replaced by 2 sites to the left, and the trajectory of the particle receives a displacement of 1 site to the right.
particle and a bound state with a given number of down spins. For example if the bound state is of length two, we obtain the exact wave function The interpretation of the wave function is the following: originally the bound state of two DW's occupies positions 2 and 3, and there is an incoming wave. Afterwards the bound state occupies positions 0 and 1 and the outgoing wave also suffers a displacement of 2 in the forward direction. This is shown in figure 2. Notice that the hard core property is still satisfied, although it is not enforced, it is completely dynamical even in this case. This state also has energy E = e(p), the two domain walls do not contribute.
The states (27) and (28) only exists in infinite volume, because the periodicity conditions can not be satisfied by them. However, the state (28) can be modified to fit into a finite volume situation with periodic boundaries. To this order we introduce a Fourier transform over the position of the block of 2 spins and write |Ψ = y>x+1 e i(px+k(y−2)) |x, y, y + 1 + x>y+2 e i(p(x−2)+ky) |y, y + 1, x . (29) We can extract from this wave function the scattering phase The state (29) has energy E = e(p) and it is periodic if p and k satisfy the equations Substituting (30) into the quantization conditions and introducing once again the total momentum P = p + k we obtain The domain walls do not contribute to the energy, but they modify the propagation of the free particles, and also the quantization conditions. The effect of the interaction is a change in the apparent length of the system by -3 and the appearance of the total momentum as a twist for the quantization conditions. It turns out that this picture can be generalized to many body states: the model is integrable, and scattering of the particles among each other and on the domain walls leads to the same phases that we computed in this Section. In particular the scattering phase between a particle and a frozen block of spins of length n does not depend on n.
The complete integrability of the model is most easily understood by considering a special connection to the XXZ spin chain, which we discuss in the next Section. The degeneracies resulting from the presence of the domain walls in generic states is discussed later in Section VIII.

IV. RELATION TO THE XXZ SPIN CHAIN
Here we show that our model can be derived as a special limit of the XXZ spin chain. Our procedure is different from the one in [51,52], even though we also consider the large anisotropy limit.
The XXZ chain is given by the Hamiltonian where ∆ is the anisotropy parameter. It is known that the XXZ chain is integrable, and it possesses a set of conserved charges that will be denoted asQ α . They satisfy The chargeQ α has an operator density that spans α sites; as usuallyQ 2 can be identified with the Hamiltonian, and Q 1 can be chosen as the global S z or alternatively as Q 1 written in (3). The chargesQ α can be obtained either from a transfer matrix, or using the so-called boost operator [63][64][65]; we describe here the latter method. The boost operator is defined as the formal expression whereh j,j+1 is the Hamiltonian density acting on sites j, j + 1. Then the charges are constructed through the recursive relationsQ Even though the r.h.s. is just a formal expression, the actual commutation relations result in a well defined extensive and local charge at each new step. This relation and the resulting charges were studied in detail in [65,66], where concrete formulas were given in a number of different cases. Furthermore, explicit formulas for Q α with arbitrary α were derived in the recent works [67,68]. We obtain our model by performing a special ∆ → ∞ limit on the set of the conserved chargesQ α . It can be seen from the boost relation that in the XXZ model the operatorsQ α are polynomials in ∆. Our main idea is to select the leading terms in ∆ using a recursive procedure. Then the commutativity (34) will then ensure that our new charges Q α also commute.
We start with Q 1 which does not depend ∆, thus it will stay constant during the limit. The next charge is Q 2 , which is the Hamiltonian (33). We select the leading piece in ∆, which is (apart from the additive normalization, and stripping away the factor of ∆) equal to Q 2 as given in (3). It is clear that [Q 1 , Q 2 ] = 0. The charge Q 2 is not dynamical: the kinematical piece of the orig-inalQ 2 is scaled to zero, and only the "classical" part remains, which describes the classical Ising model. The disappearance of the dynamical terms in Q 2 prompts us to turn to the higher charges of the XXZ chain.
Therefore we take the explicit representation ofQ 3 and Q 4 found in [65,67,68] and we select the leading terms in ∆. It turns out that in the normalization of [65] the maximum power of ∆ inQ 3 is linear, whereas it is quadratic inQ 4 . Selecting these terms and stripping away the factors of ∆ we obtain Q 4 and a new charge It is clear from the construction that all four charges commute with each other, because their commutativity is simply the leading order term from the relation (34). And Q 3 and Q 4 are already dynamical, even after the ∆ → ∞ limit has been taken. The numerical pre-factors in Q 3 and Q 4 are added only for convenience, such that they have simple one-particle eigenvalues. Note that Q 3 is actually a three-site operator, it is just written conveniently in a four site representation as above.
Continuing to even higher charges we observe that more steps are needed. For example we find that the leading term inQ 5 is of order ∆ 2 , and it is proportional to Q 3 as given above. Therefore a new charge can be obtained only if we subtract fromQ 5 the appropriate multiple ofQ 3 and then consider the leading term of the remainder. This procedure leads to a new operator Q 5 which commutes with all previous charges and which is given (after stripping away an irrelevant numerical prefactor) as Q 5 appears to be a 6 site operator, but this is only apparent: after expanding the products it can be seen that every single term spans a maximum of 5 sites only. It follows from the construction that Q 5 commutes with all previous charges of our model. The charge Q 5 is identical to the operator Q − 2 of [51], see formula (31) in that paper.
We can continue this recursive procedure to obtain the higher charges Q α with α ≥ 6. The idea is always to focus on the leading terms in ∆, and to subtract contributions that appeared in the earlier charges. We conjecture that this procedure leads to an infinite set of linearly independent charges Q α .
At present we do not have closed form results for the higher charges, and we do not have a direct transfer matrix construction for the new charges either. However, the existence of the charges Q α with α = 1 . . . 5 is already enough to claim the integrability of the model. In fact, for integrability it is enough to have 2 independent dynamical conserved charges [69], which in the present case are Q 3 and Q 4 . The results of the works [67,68] could be used to give explicit formulas for all the charges of our model, but this is beyond the scope of the present work.
We note a curious property of our construction: the boost relation (36) does not survive the ∆ → ∞ limit. This is most easily seen from the relation between Q 2 and Q 3 : whereas the originalQ 3 is obtained fromQ 2 using the boost, this is not true after the scaling. The new charge Q 2 is not dynamical, and it does not give Q 3 using the relation (36). However, we observe that the next relation is kept intact, which means that with B given by the boosted Q 2 . We comment further on this issue in the Conclusions.

A. Bethe Ansatz for the XXZ chain
Here we present a summary of the Bethe Ansatz solution of the XXZ chain, which is well known [70]. Our goal is to consider the ∆ → ∞ limit of the full solution, this is presented below. Throughout this Section we will use the conventional parametrization It is known that the N -particle eigenstates of the XXZ model can be characterized by Bethe rapidities λ j , j = 1 . . . N , that describe the one-particle momenta within the interacting states. The spin waves can form bound states, which are described by the so-called string solutions. The string hypothesis states that in the thermodynamic limit almost all eigenstates can be described by Bethe roots that organize themselves into the string patterns.
In our normalization λ ∈ [−π/2, π/2] and the n-string propagation factors are where we defined The scattering of an n-string and an m-string is described by Here we omitted the dependence on λ in the notations. We denote the string centers for the n-strings as λ j,n . Then the approximate Bethe equations for the strings are The energy of such an eigenstate is wherẽ In a normalization dictated by the boost operator construction (see for example Section 3 of [71]) the eigenvalues ofQ 4 areQ or with a more explicit formulã B. The special limit of the Bethe Ansatz Here we consider the ∆ → ∞ limit of the previous formulas. This particular limit of the XXZ chain was already studied in a series of works [72,73]. Also, it was known that in this limit the constrained XXZ chain can also be obtained [54,55]. However, the discussion of the dynamics under Q 4 is new.
First let us discuss the interpretation of the strings. It is known that in the ∆ → ∞ limit the bound states become more and more bound [74]. Eventually in the strict ∆ → ∞ limit an n-string will describe n down spins placed beside each other. Thus it is natural in this special limit that the string solutions will correspond to the bound states in our model. Now we compute this correspondence, and the associated scattering phases.
We put forward that one limitation of this Bethe Ansatz picture is that it does not mirror the spin reflection symmetry of our model. As a consequence, the simple scattering state (27) of a particle and a DW is a very complicated object in this picture: the sea of down spins on the right would be described by an infinitely large string. In Section VI we present an alternative representation of the same model, where a single DW is treated as a separate particle. This picture will avoid some of the complications of the present description. Now we compute the special limit of the eigenvalue functions. First of all we find that The explicit form of Q 2 given in (3) implies the limit which together with (45) gives We see that the particles and the bound states give an equal contribution to Q 2 , which does not depend on the momentum. This is clear from the actual form of Q 2 , which shows sensitivity only to the domain walls.
Regarding the scaling of theQ 4 eigenvalues we find the large η behaviour where we identified 2λ = p.
Combining with (47) and a re-scaling ofQ 4 according to we get the scaling leading to This coincides with our real space computations in Section III.
Let us now consider the Bethe equations. Note that the limit of the functions Θ n is simply where we identified p = 2λ.
For the limit of the S-matrix factor between 1-strings we find simply: This agrees with the scattering phase found earlier in (15). For the limit of the S-matrix factor between 1-strings and n-strings for n > 1 we find: This is the same factor as found earlier in (30). We stress that this factor does not depend on n: the scattering phase is independent of the length of the bound state, it only depends on the number of domain walls crossed, which is two in this case. We also compute the S-matrix between an n-string and an m-string with n, m > 1. We find Then we obtain the final Bethe equations For particles (or equivalently for 1-strings) the Bethe equations take a particularly simple form. Using the above results we get This is then rewritten as where we introduced the total string number (without the particles) and We see that for the quantization condition of these particles the various strings are grouped together, such that the information on the string length disappears. Only the total number of the strings and their total momentum enters the equations. This is perfect agreement with other descriptions of the model which will be discussed below.
The energy of the states is We observe that degeneracies between different string lengths remain if h = 0. This is discussed in more detail below in VIII.

V. BOUNDARY CONDITIONS
In this Section we investigate the model with open boundaries; this is the case that was originally considered in [50]. We will show that the boundary problem is even simpler than the periodic one: it gives simpler Bethe equations, we can easily prove completeness of the Bethe Ansatz, and the study of the thermodynamics is also easier (see Section X).
This simplicity can be explained by a standard semiclassical picture. We consider the motion of the particles in a finite volume, both in the periodic and in the boundary case. In the periodic case the particles move along the chain and eventually travel around the full volume. In this process every particle hits every DW one time, therefore the domain walls are displaced. This process is then repeated over time, thus the DW's can not have fixed positions. In contrast, in the boundary case the particles travel back and forth between the left and right boundaries. In this process every particle hits every DW one time from the left and one time from the right, and the displacements of the DW's add up to zero. In other words the DW's stay in their original positions as the particles finish a full cycle of traversing the available volume. This means that in the boundary case we can label the states with the DW positions and the particle momenta. This phenomenon shows that the open boundary condition is "more compatible" with the bulk scattering than the periodic one.
In the boundary setting the Hamiltonian of a chain of length L + 2 is (67) This Hamiltonian can be obtained from a double row transfer matrix (using the identity as the so-called Kmatrices) of the XXZ spin chain in the same way as for the periodic case (see section IV). Using the same argument as section IV we can convince ourselves that the Hamiltonian is integrable. We will not treat this procedure in detail here, and for the algebraic framework of the boundary XXZ chain we refer the reader to [75].
The following charges are also conserved for open boundaries In addition, specific to the boundary case, we have the following conserved charges We see that the first and the last sites are not dynamical and the Hilbert space is decomposed as Each subspace has dimension 2 L . In the following we focus on the subspace H ↑↑ . The other sectors can be analyzed similarly.

A. Bethe Ansatz
In this subsection we present the characterization of the spectrum of H ↑↑ . Let us start with the one-magnon state. The Bethe Ansatz eigenstate reads We can calculate the energy, the reflection factor R(p) and the Bethe Ansatz equation in the usual way, which leads to Solutions to be Bethe equation are The number of solutions agrees with the dimension of the one-magnon subspace of H ↑↑ . We can generalize this Ansatz to general states. We learned from the periodic case that the spectrum contains magnons and DW's. The width of the DW's have to be at least two since the one distance means there is a hole excitation which will be transformed to a magnon at some point.
There is now a key difference as opposed to the periodic case; above we have already touched upon this issue when discussing the semi-classical picture. In the boundary setting we can have bound states with fixed positions, up to the scattering displacements with the particles. In contrast, in the periodic case we needed a Fourier transform over the bound state positions such as in (29), otherwise there was no way to satisfy the periodicity conditions. At the heart of this problem lies the periodic property itself, which makes it impossible to identify which particle is to the left or to the right as compared to any other particle. In the boundary setting this is clearly possible, and we can assign well defined positions to the DW's, which are changed only as a result of the scattering with the particles. Now we present the construction of the Bethe states. We present the final results without proofs; the form of the wave functions follow from the previous results given above. A long and detailed proof would not contribute considerably to the understanding of the model. However, to avoid simple mistakes we checked the wave functions below in small volumes and found that they indeed produce the eigenstates of the boundary model.
The wave functions of the N particle states without DW's are The summation over the signs s j correspond to whether the particle is moving from the left boundary to the right one or vice versa.
Going further we can also insert bound states into the chain, which will be displaced as particles scatter on them. Such states will thus be given by a complicated summation over the various possibilities for the positions of the particles. In any case the states can be labeled by the positions of the DW's as the particles occupy the leftmost possible positions. Thus we can label the states with domain walls at positions ( where y i − x i ≥ 2, x i+1 − y i ≥ 2, and the dots stand for other components of the wave function. To describe the full wave function with N particles and M DW's it is best to consider the motion of N free fermions in a box with effective lengthL = L + 1 − N − 2M . It is convenient to use effective positions 1 ≤ x 1 < · · · < x N ≤L and real positions 1 ≤x 1 < · · · <x N ≤ L. The effective position is the position of the particle in the auxiliary free fermion picture and the real position is its actual position. Ignoring the DWs we saw that the connection between real and effective positions isx k = x + k − 1 (see (76)) which is a manifestation of the hard rod property (see also Section XIII).
We also saw that the DWs are displaced by the scattering. Let (ỹ + j ,ỹ + j ) be the position of the DW when the effective positions are {x k }. Figure 3 shows the connections between these positions if we have one particle and one bound state. We can see that the real position is increased by one when the particle reaches the left of the domain wall and one more when it reaches the right. We can also see thatỹ − andỹ + are decreased by two when the particle reaches the left and right respectively.
Having more particles and bound states we have to count how many particles reached a DW to get the real positions of given effective positions. If the kth particle reaches the left of jth DW then the k + 1, . . . , N particles have to already reached it therefore its position when the kth particle reaches it isỹ − l = y − l − 2(N − k). The kth particle can reach left of jth DW only when it already went through 1, . . . , j − 1 bound states therefore its real position isx k = x k + k − 1 + 2(j − 1) therefore the kth particle pasted the left of the jth DW if x k −y − j +2N −k+ 2j −2 ≥ 0. In an analogous way the kth particle is pasted the right of the jth DW if x k − y + j + 2N − k + 2j − 1 ≥ 0. In summary the full wave function can be written as (Θ is the unit-step function) and We can now easily calculate the Bethe equations for the particles. Let us pick up a particle and move it through the chain. It will hit every other particle and bound state two times. The single time scattering phases are −e −i(pj ∓p k ) and −e −2ipj for magnons and bound states. We can see that the momenta p k drop out from the Bethe equations since the products of the scattering phases are e −i(pj −p k ) e −i(pj +p k ) = e −i2pj , and therefore the Bethe equations are where N and M are the number of the magnons and the bound states. We can see that the Bethe equations are decoupled and they are the same as the one magnon Bethe equation with lengthL = L + 1 − N − 2M . Therefore we obtained an effectively free theory with a modified volume. We can check that the above defined states span the entire Hilbert space H ↑↑ . Fixing the number of particles and bound states it is obvious that we can place the bound states in several ways. The number of these possibilities is The effective length for the magnons isL = L + 1 − N − 2M , therefore the number of the solutions of the Bethe equations is From (85) and (86) we can count number of states we described above (87) We can convince ourselves that We can see that we created as many states as the dimension of the Hilbert space. We checked numerically that the free Bethe equations above with the bound state degeneracies as given above produce the complete spectrum of the open spin chain up to L = 8.

VI. BOND-SITE TRANSFORMATION
The goal of this Section is to build a different representation of the same model, such that the standalone Domain Walls can be interpreted as particles (of a new particle species). We perform a non-local transformation: we put dynamical variables on the bonds between sites, and build a dynamical model for the bonds. It turns out that our original charges Q 4 and Q 2 can be represented by local operators after the transformation. The advantage of such a representation is that it leads to a simpler Bethe Ansatz description, with only two particle types. This Bethe Ansatz naturally respects the spin reflection invariance of the original model. However, the advantages come at a cost: the original charge Q 1 becomes a non-local operator in the new basis.
We put variables onto the bonds (links) between the sites, and we perform a change of basis from the old computational basis to the new one. We have two options for each link, which we denote by • and •. We put • if the two neighbours are of the same spin, and we put • if the two neighbours are different. This transformation is completely invariant with respect to a global spin flip, which is an invariance of the operators Q 2 and Q 4 . This transformation is identical to the one used in [51,52] to derive their dual Hamiltonian; this is most easily seen from footnote 4 on page 11. of [51], where the transformation rule for the σ z operators is given. The interpretation as a bond model is new.
We can define this bond model on a periodic lattice, but for simplicity we first consider the boundary setting. We consider the sector of the original model where the spin at site j = 1 is in the up spin position. Then we construct the bond basis for the L = L − 1 bonds, without any restriction on the last site. This way the Hilbert space will have a dimension of 2 L .
Let us now give the operator representation of the charges in the bond basis. In terms of spins, we can interpret |• as the up spin, and |• (the particle) as the down spin, and below we will use the Pauli matrices in this new basis accordingly.
Under this transformation the charges Q 2,3,4 become the local operators where P • j is the projector onto the state |• on site j. This form of the charges is obtained by direct computation. They were already given in [51,52] using a transformation on the level of the operators.
Interestingly, we can build the non-Hermitian combinations Q B 4 ± iQ B 3 , that describe propagation terms towards the left or to the right only.
A disadvantage of the bond picture is that Q 1 given in (3) is not a local operator anymore. Instead, it is given by the highly non-local expression As noted earlier, we can also define the bond model with the operators above, assuming periodic boundary conditions. In this case the bond model is not completely equivalent to the original spin chain: for example, it allows the presence of an odd number of domain walls, which were forbidden by definition in the old chain. Nevertheless, the sectors with an even number of DW's are equivalent to those of the original chain. In these sectors every state of the bond model actually describes two different states from the original chain, which are related by spin reflection.
Sectors of the bond model with an odd number of DW's can be accommodated in the original spin chain if twisted boundary conditions are chosen, with a twist given by spin reflection. Such a model is also interesting on its own right [76,77], but we do not discuss it here.

A. Dynamics and Bethe Ansatz in the bond picture
In the bond picture a single • represents an original DW, while two bullets on neighbouring sites represent an original particle. Correspondingly, the kinematical terms in Q B 3 and Q B 4 move the double bullets, but they leave the single • invariant. To be more precise, the only non-vanishing matrix elements of Q B 3 and Q B 4 are those corresponding to the moves We can regard the DW as the fundamental excitation, and the original particle (the doublet ••) as a bound state of two DW's. The phenomenon that a single excitation is stable but a bound state of two excitations is dynamical is similar to the situation in fracton models (see the reviews [59,60]). Let us now construct the Bethe Ansatz wave functions in the bond model. We set h = 0 (thus we discard the non-local charge Q B 1 ) and consider the local Hamiltonian with periodic boundary conditions. The construction below is not rigorously proven: we construct the Bethe Ansatz by assuming factorized scattering and using the S-matrix factors extracted from the two body problem. We put forward that the Bethe Ansatz is not unique: different constructions lead to different choices for the basis vectors in the highly degenerate subspaces. We discuss this issue at the end of this Section. We have two types of excitations in the model: the single • which is a DW, and the •• which is a particle. We will use the notations p and DW .
Correspondingly we introduce local creation operators Note that we have automatic exclusions: Let us consider a state with N particles and M DW's; the set of the particle and DW momenta will be denoted as p N and k M . The wave function is most easily written down by merging these sets. Therefore we introduce a set of momenta q N and a set of particle types a N with N = N + M . We assume that there are no coinciding rapidities.
Note that the S-matrix is such that for DW's the occupation of neighbouring sites is forbidden. This ensures that we do not mistake two domain walls close-by with a particle. Also, if we have two particles at positions x 1 < x 2 , then x 2 = x 1 +1 is forbidden by the action of the creation operators, but the next possibility x 2 = x 1 + 2 is allowed. Furthermore, if we have a DW at x 1 and a particle at x 2 then x 2 = x 1 + 1 is allowed, but the other ordering x 1 = x 2 + 1 is forbidden. Thus, in this wave function a sequence of |• • • embedded in the vacuum is interpreted as a DW and a particle from the left to the right. It is merely a choice which follows from our definition of the creation operators.
The energy of this state is The sum runs over the particles only. The domain walls only contribute to the energy through the chemical potential µ.
The Bethe equations for the p and k variables are Substituting the factors we get where we defined The energy is carried only by the particles, and the effect of the domain walls is only a change in the available volume. Correspondingly, the actual values of the k j variables do not matter for the particle momenta p j , which is influenced only by the sum K. Thus the distribution of K among the variables k j only contributes to the degeneracies of the energy levels. This can be seen more explicitly by taking the product of the second set of equations, which leads to the following coupled equations for the variables p j and K: Explicit solutions to the equations above are found as follows. The overall momenta are expressed as and where I and J are arbitrary integers. Afterwards the particle momenta p j are expressed as where quantum numbers I j are given by I = N s=1 I s : 0 < I 1 < I 2 · · · < I N < L.
Let us compare these equations with (59) that was derived using the strings of the XXZ model. We can see in (59) that the contribution to the volume change is the same for every string, thus for every bound state. In the bond picture this means that the volume change does not depend on the relative position of the DW's, just their total number. This is consistent with the equations above. However, the "twist" felt by the particles depends on the momentum of the DW's.
This basis is certainly different from the one obtained by the string picture in the original model. The key difference is that here the DW's are allowed to move independently, whereas in the string picture the two DW's of the bound state are always at a fixed distance. Both pictures describe highly degenerate energy levels, and such differences only amount to a free choice of the basis.
In fact, it is easy to see that the Bethe states obtained in this bond picture can not be identical to those of the original model. This happens because in the bond picture the two DW's associated to an original string solution can not have the same momenta, therefore they can not move together. Furthermore, the wave function (95) naturally gives a linear combination of states where the "string lengths" vary. States with fixed string lengths are not reproduced by this formula. However, we stress that the two descriptions merely correspond to two different choices for the basis vectors in a highly degenerate eigenspace.
It is worth mentioning that the Bethe Ansatz solution presented above is completely different from the one of [51,52], and it is not straightforward to find a dictionary between the two solutions, and most probably they lead to different basis vectors in the highly degenerate eigenspaces. We checked by numerical computations that both constructions correctly reproduce the spectrum up to L = 8.

VII. NON-LOCAL MAPPING TO THE SU (3) XX MODEL
Here we show that the bond model can be mapped to the Maassarani-Mathieu (MM) spin chain, which is also known as the SU (N ) XX model [78]; our mapping concerns the SU (3)-related case. The mapping is similar in spirit to the non-local mapping between the phase model and the XX model which first appeared in [46] and which was treated detail in [48]. The mapping for the present model is new, and it is different from the connection with the Bariev model [79] found in [51].
We start with the bond model with length L in the boundary setting. This model will be mapped to a spin chain with local dimension 3, with local basis states denoted as |1 , |2 and |3 . As before, we construct the mapping on the level of the basis states. We will see that the mapping is volume changing: different states are mapped to states of the new model with different lengths.
The rules for the mapping in the computational basis are as follows. We represent each basis state as a sequence of •'s and •'s. We proceed from the left to the right of this sequence and we translate it to a new sequence consisting of the numbers 1, 2, 3. If at a given position we encounter a • then we add a 1 to the new sequence. Then we move further to the next entry. If we encounter a • then we also need to check the next entry: a pair •• is then mapped to 2, whereas the pair •• is mapped to 3. This rule is then applied as we proceed along the chain.
This mapping works flawlessly on a half-infinite chain, but it runs into a problem on a finite chain if the last entry is a single • for which no rule is specified. Here we are not concerned with the boundary conditions, we are focusing on the mapping on the bulk, thus we discard this problem. It is possible that an exact mapping with open boundary conditions could be constructed, but we leave this problem to a future work.
The mapping backwards is more easily summarized as Let us discuss the changes in the length. If the numbers of the components 1, 2 and 3 in the new sequence are given by N 1,2,3 , with the total length being then the length of the original chain is We see that different basis states of our model are mapped to MM chains with different length. It is useful to note that the original particle and DW numbers are thus the volume of the MM chain for a given N, M reads We focus on the Hamiltonian in the bond model given by Q B 4 in (89), which results in the transition rules (91). It can be checked that after the basis transformation this translates into the transition matrix elements We also write down the Hamiltonian that encodes these transitions. Let E a,b j with a, b = 1 . . . 3 be the elementary matrices that contain a single 1 in row a and column b, acting on the local Hilbert space on site j of the new model. Then the MM Hamiltonian is [78] (111) This is now written down for periodic boundary conditions and we included a factor of -1/2 to match our previous normalizations. We note that this model is also equivalent to the infinite interaction limit of the Hubbard model (also known as the t-0 model), see for example [80,81].
It is useful to consider the specific sector of the bond model, which includes only particles. This sector is equivalent to the constrained XXZ model treated in Sec. II A. Applying the non-local mapping to this sector we see that only the local states 1 and 2 are populated. In this sector the MM chain is equivalent to the standard XX model. In [48] it was shown that the so-called phase model is equivalent to the XX model by a similar non-local transformation. We have thus obtained three different models that are equivalent to each other via the non-local mapping.
We stress that the equivalence of these models is established only for the bulk, because generally the boundary conditions spoil the mapping. This situation is analogous to the one treated in [48]. As mentioned above, an exact mapping could be found perhaps with open boundary conditions, but we do not pursue this direction here.
We now briefly discuss the spectrum of the MM Hamiltonian (111) with periodic boundary conditions, and the construction of its eigenstates through Bethe Ansatz. An algebraic treatment was presented in [78], and here we discuss the coordinate Bethe Ansatz solution. In fact we will see that there are two different constructions, which are related to each other by a particle-hole transformation.
The first version of the coordinate Bethe Ansatz is close in spirit to the usual nested Bethe Ansatz, where we start with a chosen reference state and build excitations above it, such that the orientation of the excitations will involve the "nesting", or second Bethe Ansatz. We start with the sector where only the local states |1 and |2 are populated; in this sector the MM model is equivalent to the standard XX model. The eigenstates in this sector are therefore constructed as excitations on top of the pseudo-vacuum |vac = |2 . . . 2 . They are parameterized by a set of momenta k = {k 1 , . . . k N }, and take the form where we have used the notation |1 x1 . . . 1 x N = E 2,1 j1 . . . E 2,1 jn |vac . Periodicity of the wave function then imposes the quantization condition e ikj L = (−1) N −1 for each j, and for a given value of N the momenta k 1 , . . . k N may take any of the L distinct values compatible with this condition, provided that they are all distinct. The corresponding energies are given by Turning to more general sectors, that is, with N 3 = 0, we introduce the states |k α1,...,α N , with α j ∈ {1, 3}. These have the same form as (112), but where the sequence 1 x1 . . . 1 x N is replaced by the corresponding sequence of 1s and 3s (we then have N = N 1 +N 3 ). Since the order of 1s and 3s along the chain is not modified by the transition matrix elements (110), the states |k α1,...,α N are locally eigenstates of the MM Hamiltonian, with energies again given by (113). However, on a chain of finite size requiring periodicity of the wave function mixes states with sequences α 1 , . . . , α N related through cyclic permutations. We therefore construct eigenstates as linear combinations of the form N n=1 e iqn k U n (α1,...,α N ) , where U cyclically permutes the indices α i . The pseudomomentum q is quantized through e iqN = 1, but depending on the sequence α 1 , . . . , α N only a subset of the solutions for q might give rise to a non-vanishing wave function. In other terms, constructing the linear combinations (114) amounts to diagonalizing the one-site translation operator on an auxiliary spin-1/2 chain of N sites, in the sector with N 1 spins ↑ and N 3 spins ↓. Imposing the periodicity of the wave function, we now have the following quantization conditions  In summary, the eigenstates in a sector of given N 1 , N 2 , N 3 are obtained by first, finding the allowed values of q by diagonalizing an auxiliary spin-1/2 problem, and, second, solving the quantization condition (115) for the momenta k 1 , . . . k N . This is a particularly simple form of nested Bethe Ansatz, as the quantization condition of the auxiliary momenta q does not depend on the values of the momenta k 1 , . . . k N . We checked against exact diagonalization for finite size systems that this construction indeed reproduces the entire spectrum of the periodic MM model. Now we turn to the second construction for the coordinate Bethe Ansatz, which is in fact closer in spirit to the bond model presented in the previous section. Rather than using a unique pseudo-vacuum, now we start with a collection of pseudo-vacua |(α 1 , . . . α N ) made of arbitrary sequences of 1s and 3s, which are all zero energy eigenstates of the MM Hamiltonian. On top of a pseudovacuum |(α 1 , . . . α N ) , we create excitations by inserting local states 2 in between the 1s and 3s. The resulting Bethe wave functions take the form of the usual XX wave functions on a chain of L = N 2 + N sites, analogous to (112) but where the locally excited sites 1 xi are now replaced by 2 xi , and the states on other sites are distributed according to the sequence (α 1 , . . . α N ). As in the previous construction such wave functions are locally eigenstates of the Hamiltonian, but periodicity now mixes excitations on top of different pseudo-vacua, namely related by cyclic permutations. Eigenstates of (111) are therefore obtained as linear combinations of excitations over cyclically permuted pseudo-vacua, involving as before a pseudo-momentum q with is a N th root of unity. The resulting quantization conditions for the momenta k j of the "2" particles take an analog form to (115), more precisely the two sets are mapped onto each other by a particlehole transformation of the underlying XX model. Note the key difference that now the local states 2 are the excitations, whereas they formed the reference state in the previous construction.
It is useful to compare (115) to the Bethe equations obtained in the bond model, see eq. (101). We observe some similarities, for example the apparent volume for the p-variables is the same. However, the twists appearing in those equations are different. This is a consequence of the fact that the non-local mapping is not compatible with the periodic boundary conditions in these models.

VIII. DEGENERACIES AND HILBERT SPACE FRAGMENTATION
Here we discuss two closely related features of the models: Hilbert space fragmentation and a large number of degeneracies present in the spectrum. These questions were already discussed in [50][51][52], here we summarize the key statements and present a complementary view of the matter.
The expression "Hilbert space fragmentation" means that there are a large number of sectors in the Hilbert space such that the Hamiltonian does not have transition matrix elements between the sectors. It is also required that the sectors should be constructed using relatively simple rules, before actually solving the full dynamics in the model. In a typical case there is an exponentially large number of disconnected sectors in the Hilbert space. Fragmentation is known to happen in models with fractonic excitations [59,60].
One of the mechanisms for the fragmentation is the presence of conserved charges [82,83]. If there is a U (1) conserved charge of the model, then its eigenvalues already split the Hilbert space into various sectors, however, such a splitting is very common and natural, and it is usually not called "fragmentation". In contrast, fragmentation can happen in the presence of two conserved charges with local densities, such that the charges are not dynamical and their densities can be diagonalized simultaneously. Examples with a charge and a dipole conservation were discussed in [82,83].
In our case we have two such charges Q 1 and Q 2 , which commute with the dynamical charge Q 4 . Accordingly, the eigenvectors organize themselves into sectors corresponding to the eigenvalues Λ 1,2 of the charges Q 1,2 Furthermore, we observe that even the sectors with a given eigenvalue pair (Λ 1 , Λ 2 ) further split into subsectors, whose number grows exponentially with the volume. These sectors correspond to the presence of the domain walls placed at various distances from each other.
An intuitive way to understand the phenomenon is to first consider the completely frozen states which have an eigenvalue 0 under Q 4 . As explained above, such states can be created by placing an arbitrary number of domain walls on top of a reference state with minimum distances of at least 2. Then the particle excitations are created above such a frozen state, such that the Hilbert space remains fragmented into these sectors.
Despite the simplicity of this picture, it is not completely precise in a finite volume situation with periodic boundaries: Even though the domain walls themselves are frozen, they are displaced once a scattering with a particle occurs. Therefore, we need to impose the proper periodicity conditions also for the movement of the DW's, resulting in the full set of Bethe equations. In contrast, the wave function given in Section V shows that in the boundary case the fragmentation can be understood using this intuitive picture, because it is relatively easy to treat the displacements of the bound states.
Let us now also discuss the pattern of degeneracies in the model. We study the most general Bethe Ansatz equations (61). First we discuss the degeneracies at h = 0. The energy is carried only by the particles, thus it depends only on the set p N1 , and the global quantum numbers N 1 and N s . The Bethe equations for p N1 do not depend on the distribution of string lengths; they are sensitive only to the total string number N s and the total string momentum. Thus we observe a large degeneracy, which is expected to be exponentially growing with the system size. In the bond picture this degeneracy results from an arbitrary placement of the DW's, which does not affect the quantization conditions for the particles.
The exponential growth of the degeneracies strongly depends on the particle content of a state. Each domain wall decreases the space available for the propagation of the particles, thus a given particle number N 1 also constrains the possibilities for the DW's; this was explicitly demonstrated in the boundary case in Section V above. For the reference state the exponential growth was computed in [50] and also in [51], and it was found that the degeneracy behaves in large volumes as α L , where α = (1 + √ 5)/2 is the golden ratio. In contrast, the ground state has a finite degeneracy even in the thermodynamic limit (see below).
Let us also discuss the splitting of the energy levels as a magnetic field is switched on. The additional term in (66) implies that most of the degeneracies are split due to the various distributions of the total string number N s into the different strings. However, an exponential amount of degeneracy still remains, corresponding to the various ways of obtaining the same N s and the same total magnetization, and also to the various ways of solving the Bethe equations for the strings.

IX. GROUND STATE
Here we discuss the nature of the ground state of the original Hamiltonian given by (1) and we compute the ground state energy density. We consider the cases h = µ = 0 and also the situations when either or both are switched on.
For the case of h = µ = 0 our results agree with those of [51] and in the earlier work [56]. Simple arguments and numerical checks show that the ground state is populated by particles only, but their overall density (or equivalently, the Fermi boundary) is a non-trivial quantity. The ground state is doubly degenerate, corresponding to the spin flip invariance, and here we treat the state that has positive overall magnetization (N < L/2).
The Bethe equations are: Taking the logarithm, we get where I j ∈ Z. Since in the ground state P = N k=1 p k = 0, this simplifies to HereĨ j = I j + (N − 1)/2 is an integer/ half integer for odd/even N . In the ground state theĨ j -s are distributed symmetrically around zero, from − N −1 2 to N −1 2 . The energy of the state is In the thermodynamic limit (L, N → ∞, N/L = n is fixed) the energy density becomes: where we used that The particle density that minimizes is found from which leads to This equation can be solved numerically and gives n ≈ 0.3008. The energy density is ≈ −0.2172. If we include bound states, their only effect for the particles is decreasing the effective length by 2M , where M is the number of bound states (with arbitrary length). Therefore the energy density in a state with particle density n and bound state density M/L = m is: This expression has its minimum at n = n 0 ≈ 0.3008 and m = 0. Since the Q 2 term gives the same energy to a particle and a bound state and the Q 1 term energetically prefers particles over bound states, the situation will not change even if we turn on non-zero h and µ. The ground state does not contain bound states. Let us now investigate the case with h = 0 but still µ = 0. We can restrict ourselves to h > 0 by spin reflection invariance. In the thermodynamic limit the energy density of the system in a state characterized by a particle density n is The value of n that minimizes this energy density is given by the relation The expression on the r.h.s. has a maximum at n = 0 with a value of 1. This means that if h < 1, the ground state is characterized by a finite particle density n 0 , given by (126). On the other hand, if h > 1 than the reference state with all spins up becomes the ground state. By expanding (126) around 0, we can calculate how n goes to 0, as h approaches h c = 1. We find the scaling Now we consider the case of h = 0, µ = 0. It is still enough to consider h > 0, but µ can be both positive and negative. The energy density in a state characterized by particle density n is By taking the derivative and rearranging the equation we get The expression on the r.h.s has a maximum at n = 0 with a value of 1, and has a minimum at n = 1/2 with a value of -2. Therefore, if h + 2µ > 1, then the reference state with all spins up becomes the ground state. On the other hand, if h+2µ < −2, the ground state will be given by n = 1/2 which corresponds to the doubly degenerate Néel and anti-Néel states. Between these two regions the ground state is characterized by a finite n, given by the following constraint: Here we investigate the Gibbs states of the model, and compute the free energy and particle densities as a function of the temperature and chemical potential. We present three different computations using the various formulations of the model, all leading to the same free energy density.

A. Thermodynamics in the boundary case
Here compute the thermodynamic limit of the boundary chain discussed in V. In a certain sense this is the simplest case which leads to an exact determination of the free energy density, without recourse to additional assumptions.
For simplicity we focus on the H ↑↑ sector of the open chain. We saw that the spectrum can be built from particles and bound states. Let us define the particle densitỹ ρ(p) for which (L − N − 2M )ρ(p)∆p gives the number of particles in the interval ∆p. Note that his is an unconventional definition because it uses the modified (apparent) volume. Using this density the number of particles can be written as (p)dp.
where the energy term is π 0 e(p)ρ(p)dp (133) with e(p) given by (9). The Yang-Yang entropy is and the entropy S DW comes from the degeneracy (85) Let us introduce a Lagrange multiplier for the particle number (131) Taking the functional derivative with respect toρ we obtain that therefore we obtain the usual free fermion density where the Lagrange multiplier acts as a chemical potential: The advantage of the multiplier is that we can change the functionals to simple functions The particle numbers n, m and a Lagrange multiplier are given by ∂f ∂m : 0 =f 0 (α, n, m) − nα where f 0 (α, n, m) = E(α, n, m) − T S Y Y (α, n, m).
We can simplify the expression of f 0 as × π 0 dp π log 1 + e −β(e(p)−α) + αn. (144) Therefore the solutions of the following equation system gives the finite temperature ground state (n, m, α) π 0 dp π 1 e β(e(p)−α) and π 0 dp π log 1 + e −β(e(p)−α) We solved this equation system numerically and the result is diaplayed on Figure 4. In the T → 0 limit it agrees with the ground state result given in Section IX. Substituting (145-147) into (136) we obtain that the free energy can be written as If someone is interested only on the free energy and not on the values of n and m then the equation system simplifies to a single equation. Let us use the following variable where we used (146). The equation (147) can be rewritten as π 0 dp π log x + e −βe(p) = 2 log(x − 1).
The free energy can be also expressed by x as  The degeneracy of states with these quantum numbers is therefore in the thermodynamic limit the DW entropy is The free energy is (p)dp − n + h(n + 2m 2 + 3m 3 + 4m 4 + . . . ). (156) We can see that theρ and the α dependent parts are the same as before with the same constraint. For the derivatives with respect to n, m k we obtain the following equations ∂f ∂m k : f 0 (α, n, m) = nα Making proper subtractions we obtain therefore The number of all DWs can be expressed as therefore we can expressed all m k and y as At this point we only have five parameters (α, n, m, q, y) which satisfy the following system of equations π 0 dp π 1 e β(e(p)−α) π 0 dp π log(1 + e −β(e(p)−α) ) Let check the h → 0 limit. In this limit the equations (168,169) simplify as (170) Substituting (166) and (167) we obtain the h = 0 equations (146),(147) exactly. Figure 5 shows the numerical results for various h. The blue, orange and green curves are the values of n, m and (L/2 − S z )/L w.r.t. the temperature. The graph on the left shows the h = 0 case which was already plotted on Figure 4, but now we are able to calculate the total spin for the thermal states. We can see that the total spin is zero for every temperature. In the second plot we can see the h = 0.5 case. The low temperature limit agrees with ground state analysis i.e. the ground state is characterized by a finite particle density n 0 ≥ 0. The third picture shows the h = 1.5 case. Now we can see that the n and the total spin goes to 0 and L/2 for zero temperature as we expected from the ground state analysis. We can also see that the DW density m goes to zero at T → 0 and the T → ∞ limit agrees for all h.
The equations are simplified again if we use only the variable x. Now we define it as Since the equation (167) can be written as π 0 dp π log(x + e −βe(p) ) = 2 log F h (x), where The variable q can be expressed from (168) which reads as (176) Substituting (174) we can obtain that The free energy can be also expressed with x as

B. Thermodynamics in the periodic case
An alternative way is to consider the periodic model and to derive the Thermodynamic Bethe Ansatz (TBA) equations based on the Bethe equations (61). Then we obtain a set of non-linear integral equations, as usually in the TBA framework [70]. However, the present case is different from the situation in the XXZ chain, because the scattering kernels are constants. This leads to rather simple equations as opposed to the full XXZ problem.
In the thermodynamic limit we consider again the Bethe root and hole densities ρ n (p), ρ h,n (p), but now with the usual normalization such that in volume L the number of particles/holes for n-strings between p and p + dp is Lρ n (p)dp and Lρ n,h (p)dp. The densities satisfy ρ k (p) + ρ h,k (p) = δ k,1 + π −π dp 4π (ρ h,k−1 (p ) + ρ h,k+1 (p )) . (179) This equation can be derived from (61) using the same steps that lead to similar decoupled equations in the XXZ chain for ∆ > 1, see [70]. We see that the total densities are constants, because the r.h.s. of these equations does not carry a p-dependence.
In the finite temperature situation we introduce the so-called Y -functions as Y n = ρ h,n /ρ n . They satisfy the TBA equations log Y n (p) = β cos(p)δ n,1 It is easy to see that only Y 1 (p) is a non-trivial function of p all the rest of the Y -functions are constants. Written more explicitly, we have the set of equations 2 log Y 1 (p) = 2β cos(p) + log(1 + Y 2 ), The TBA equations can be solved when combining with the asymptotics of Y -functions in the large n limit. For h > 0 the asymptotics is given by This implies that for n → ∞, Y n grows exponentially in n as Y n ∼ e βhn . At h = 0 the asymptotics of the Y n for large n grows polynomially in n. Finally, the free energy is given by We now discuss the solution of TBA. Let us denote The TBA equation can be written as We first consider the second equation in (185), which is a second order difference equation. The general solution to this equation is given by where a 0 and y 0 are two constants. Using the asymptotics (182), we can fix a 0 = βh/2. To fix y 0 , we consider (184), which can be written as From this equation, we can solve y 0 in terms of α. Plugging back to the first equation of (185), we obtain the following equation for the variable x = e −βα 2 log F h (x) = π −π dp 2π log 1 + x e β cos(p ) , where This result is the same as in (177). The TBA equation (188) coincides precisely with the one (173) derived in the previous section. Although (173) is derived by considering open boundary condition, the differences due to boundary conditions are negligible in the thermodynamic limit. Thus we obtained the same TBA equation as expected.

C. Thermodynamics in the bond picture
We can also compute the thermodynamics using the Bethe equations (98). A direct application of the TBA framework would lead to two coupled integral equations, corresponding to the two particle species. However, instead of copying the standard formulas of TBA we can also derive the result using a simple reasoning and formulas known from free theories.
Our first argument is that in the direct equations (98) we can assume that the total momentum of the particles and the DW's is zero; the thermodynamical computations are not sensitive to this assumption. Then we are faced with quantization conditions for two free models, such that the energy is computed only from the particles but it is insensitive to the DW's. Our strategy is similar as in the previous subsection: we compute the free energy for arbitrary overall densities n andm of the particles and DW's and afterwards we perform a saddle point analysis. For a direct comparison with the previous formulas we will usem = 2m, because in the earlier computations M and m denoted the total number and density of the bound states, whereas herem stands for the total density of the DW's. Each bound state has two domain walls at the two ends, thus the factor of 2 in the relation.
The free energy density will be a sum of the particle and DW contributions: The computation of the free energy associated to the DW's is rather simple: as there is no energy, the free energy contribution comes only from the entropy, which is computed easily. We find The free energy of the particles at fixed particle number n is computed using a Lagrange multiplier α: π −π dp 2π log(1 + e −β(e(p)−α) ) + βαn.
(192) The value of α is then found from leading to the condition This equation is the same as (145) after the appropriate re-normalizationm = 2m.

XI. SOLVABLE QUENCH DYNAMICS
The model given by (1) is special because in certain cases the real time dynamics can be solved exactly, leading to exact formulas for the time dependence of local observables. The techniques to be used are completely analogous to those used in [47,48], and they were originally developed in [43,45,46,72,73].
Let us focus on quench problems, where we prepare the system in an initial state |Ψ 0 and let it evolve with the model Hamiltonian. The simplest problems are those which involve only the particle excitations and not the Domain Walls. Such a situation arises if |Ψ 0 has strictly zero overlap with states including DW's. It is relatively easy to construct such initial states: the only requirement is that in the computational basis there can not be two or more spin excitations on neighbouring positions. For such initial states the overlaps can be computed as a sum over determinants using the exact wave function (19). The simplest cases are those when the initial state is an element of the computational basis, such that the overlaps are single determinants given by (19). We put forward that interesting physics arises in those cases where the bound states are also present; this is discussed in the next Section.
The simplest state with zero overlap with bound states is probably the Néel state where there is a particle excitation at every second site. Using the notation (7) the state is given by However, in this model the Néel state is actually an eigenstate with eigenvalue 0. The next simplest case is the period 3 state |Ψ 0 = |3, 6, 9, . . . .
Here we have N particles in a volume L = 3N . This state is not an eigenstate and the overlaps with the Bethe states are given by This follows from (19) after the substitution x k = 3k. We will see that the treatment of the quench from this state is very much analogous to the quench problem treated in Section 6 of [47]. For simplicity we restrict ourselves to the zero momentum sector, therefore we consider the initial state where U is the one-site cyclic shift operator. In this case the overlaps can be expressed simply as This overlap will be evaluated in the following.
In the zero momentum sector the Bethe equations are where we used L = 3N and we assumed an even N for simplicity. Solutions are given by The zero momentum Bethe states are given by the subsets satisfying the constraint j e ipj = 1. The numbers ω k can be paired such that It follows from formula (199) that the overlap is nonvanishing only if exactly one rapidity is chosen from each pair. Therefore, the states with non-vanishing overlap are given by with the constraint that the total momentum is zero. We have We assumed that N is even and the zero momentum constraint implies that there are a total number of 2 N −1 states with non-vanishing overlap. It can then be seen that the overlaps (199) are all equal and it was shown in [47] that their value is Our goal is to compute the time evolution of local observables O using the spectral expansion We focus on simple local observables, with our main candidates being the operators measuring the emptiness formation probability (EFP). We define the -site local operator E (x) positioned at x as We also define their space average: The operatorsĒ 1 andĒ 2 are conserved, because they are linear combinations of the charges Q 1 , Q 2 and the identity. Thus the first member of the series with nontrivial time evolution isĒ 3 .
The matrix element of a single operator p|E (x)|k for arbitrary was computed in [72,73] p|E (x)|k = where The norm of the Bethe states is given by (see [54]) The matrix elements of T can then be rewritten as It is easy to see that for diagonal part of T the rank is given by the number of coinciding elements in the sets p, k. For = 1 the non-diagonal part is equal to zero, thus the only non-zero contribution to the sum (204) is the case p = k and we obtain trivial time dependence. In the case = 2 the non-diagonal part has rank 1, thus it is possible to have a non-zero contribution to (204) in a case the diagonal part of matrix T has the rank at least N − 1. This happens if there are N − 1 rapidities p j coinciding with some of rapidities from set k. However, due to the P = 0 condition it is not possible to have a case where the two sets differ only by one rapidity. Thus, for = 2 case the selection rule for the form factors is again p = k and we have only the static contribution. These findings are consistent with the fact that E 1 and E 2 are related to the conserved charges Q 1 and Q 2 .
Finally, for = 3 the non-diagonal part of T has rank 2, thus in order to have a non-zero determinant the rank of the diagonal part should be at least N − 2. Then there is a time dependence of E 3 (x) expectation value if two rapidities are different in the sets p and k.
The further computation of dynamics for the case = 3 is absolutely similar to the one performed in [47] for the q-boson model. Thus we omit simple details and give here the final result where c a = π(2a − 1)/(2N ), a = 1, . . . , N . In the thermodynamic limit the last expression can be presented as The only difference with the corresponding formula of [47] (see eq. (6.16) of that work) is the appearance of the extra factor of 1/3, which can be traced back to the difference in the norm (209) and an overall factor of 3 originating in the definition (198) of the initial state. The integrals above actually describe the Bessel functions of the first kind, so we can write (214) The asymptotic behaviour for large t is We compared (212) to results from exact diagonalization and found complete agreement.

XII. BREAKDOWN OF THE GGE: PERSISTENT OSCILLATIONS
The large number of degeneracies in the spectrum has an interesting consequence: it can lead to the appearance of persistent oscillations in certain quench problems. Such persistent oscillations can appear in SU (2)symmetric spin chains if a magnetic field is applied, or in more complicated situations [84] where the degeneracies have a somewhat unexpected algebraic origin [84][85][86]. We will argue below that the present model also supports persistent oscillations. However, in contrast to the solvable quench treated in the previous Section, the treatment of the oscillations is much more involved from a theoretical point of view. At present we do not have an exact derivation of the oscillations, therefore we present their existence as a conjecture, which we support by numerical evidence.
In order to understand the effect we briefly discuss the process of equilibration in integrable models. Let us again assume that we release the model from an initial state |Ψ 0 and we investigate the real time dynamics of local observables O(t). Inserting complete sets of states |a and |b in finite volume we get In the long time limit (or for long time averages) it is usually assumed that this double sum can be constrained to the diagonal contributions. The idea is that in the generic situation there are not too many degeneracies in the spectrum, the energy levels E a can be considered sufficiently independent, and for large enough times the phases of the off-diagonal contributions average out to zero. This leads to the Diagonal Ensemble (DE) From this expression it is then argued, that the mean values are given by the Generalized Gibbs Ensemble (GGE), which incorporates all conserved charges of the model [6]. The derivation of the GGE is possible if the socalled Generalized Eigenstate Thermalization Hypothesis (GETH) is satisfied. If these assumptions hold then the system will equilibrate and the steady states will not support oscillations in any local observable. The situation is different in the presence of a large number of degeneracies. Let us now denote by E a the energy eigenvalues, and by a further discrete index j the states |a, j in the degenerate eigenspaces. Then we have the double sum lim t→∞ O(t) = a j,k Ψ 0 |a, j a, j|O|a, k a, k|Ψ 0 .
(218) In such a situation the long time limit is given not only by the mean values, but rather by a complicated sum involving many off-diagonal matrix elements. Generally we can expect the breakdown of the GGE in such a case. Furthermore, if the degeneracies can be lifted such that the energy differences remain commensurable, then we can observe persistent oscillations.
Let us consider concrete examples in our model. We perform quenches with the Hamiltonian (1), where we set µ = 0 for simplicity, but we consider the h = 0 and h = 0 cases as well. The larges number of degeneracies appear at h = 0, while some of these are broken up with fixed differences if h is switched on.
For the initial state our main candidate will be the ferromagnetic state in the x direction, which is given by This is a good choice because it excites various combinations of particles and domain walls. A key property is that the state breaks the U (1)-symmetry of the Hamiltonian, and after the quench the system will be populated by states with different values of the S z global charge. We choose local operators which also break this U (1)symmetry. We look specifically at These operators change the global S z charge. Therefore their mean values in the DE (or in the GGE) are identically zero.
We conjecture that in our model these observables either tend to a finite mean value (for h = 0) or they display oscillations (for h = 0). Both behaviour signals the breakdown of the GGE and the presence of the offdiagonal contributions in (218). The difference between the two behaviours is explained easily: if h is switched on then some of degeneracies are split with fixed amounts according to the global charge S z . The operator O = D k changes S z by a fixed amount, thus every non-zero contribution in the double sum (218) will receive the same phase e ihkt .

A. Numerical results
We support our conjecture with numerical evidence, obtained from iTEBD [87,88] simulations. We used the example code in [89] as a starting point and modified it to our purposes to simulate real-time evolution governed by a four-site Hamiltonian: the state of the system is represented as a four-site translational invariant matrix product state (MPS) [90] with four sets of matrices |Ψ = ...,j k ,j k+1 ,j k+2 ,j k+3 ,...
. . . × |. . . , j k , j k+1 , j k+2 , j k+3 , . . . , where both the Γ j k and Λ k (k = 0, 1, 2, 3) are χ × χ matrices, where χ is the so-called bond dimension. The latter are diagonal and they contain the singular values corresponding to the bi-partition of the system at the given bond. The algorithm uses a first order Suzuki-Trotter decomposition for the time evolution operator. We initialize our system in the state |Ψ 0 given by (219) (since it is a simple product state, it can be easily written in a form like (221)) and evolve it with a Trotter step of δt = 0.01. Unfortunately the entanglement entropy grows rapidly in time, which limits the applicability of the method to short times. To check the validity of our results we use several different maximal bond dimensions χ max . The results are shown in Figure 6. The curves confirm the theoretical expectations in the investigated time frame.

XIII. RELATION WITH TT -DEFORMATIONS
Here we point out an interesting relation with the socalled TT -deformations of QFT's. Such deformations can be introduced for any QFT, and for the integrable cases it is known that they preserve the integrability of the model [91][92][93]. The TT -deformation modifies the scattering matrix of the model, which is well understood and which we discuss below. Analogues of the TT -deformation for spin chains were discussed in [94,95] pointing out that these transformations are essentially the same as a certain type of long-range deformation studied earlier in the context of the AdS/CFT correspondence [96].
A TT -like deformation can be introduced for any pair of extensive conserved quantities Q α and Q β . To first order the deformation consists in modifying the Hamiltonian H of a given model as H = H + κ (J α (x)q β (x) − q α (y)J β (x + 1)) + O(κ 2 ), (222) where q α (x) are the charge densities and J α (x) are the current operators describing the flow of the charges [94,95]. Here we gave the formula for the spin chain situation, but an analogous expression gives the corresponding perturbation in the QFT case as well [91][92][93].
If S(p, k) = e iδ(p,k) is the two-particle scattering matrix of the model, then under the perturbation above it gets deformed as δ (p, k) = δ(p, k) where h α,β are the one-particle eigenvalue functions of the charges Q α and Q β .
In QFT the actual TT -deformation corresponds to choosing the energy and the momentum as two charges. An other important case is the so-called hard rod deformation, which corresponds to choosing the particle number and the momentum as the two charges [28,29]. In this case the modification of the S-matrix is simply S(p, k) → S(p, k) exp {iκ(p − k)} .
It is important that the parameter κ can be varied continuously, and usually it is assumed to be small compared to some characteristic length scale.
In contrast with the QFT situation it is not possible to construct the TT and hard rod deformations on the lattice. On a technical level this happens because on the spin chain there is no momentum operator which would be an extensive local charge. On a physical level this absence comes simply from the discrete structure of the chain.
Quite interestingly, the S-matrices that we encountered in this paper have the structure of (224) with an integer κ with the undeformed S-matrix being ±1. For example the S-matrix factor (15) of the particles in our model is of this form, thus it can be considered a hard rod deformation of the free S-matrix of the XX model with κ = −1. There are other similar examples for this in the literature. For example the S-matrix of the phase model (see [47,48]) is given by S(p, k) = −e i(p−k) .
This is the hard rod deformation of the XX model Smatrix S(p, k) = −1 with deformation parameter κ = +1. Similarly, the S-matrix of the constrained XXZ model found in [56][57][58] can be considered a hard rod deformation of the scattering phases of the XXZ model.  Figure 6. (a) The real part of the expectation value of D1 for h = 0 with different maximal bond dimensions. On the inset the same curves are magnified for t > 1.6. The data seems to confirm a non-vanishing asymptotic value, but we believe that this evidence is not decisive. The observation of longer times would be desirable, which is not possible with our numerical implementation. (b) The real part of the expectation value of D1 for different values of h. The maximal bond dimension is χmax = 1000, the Trotter time step is δt = 0.01. We observe the predicted oscillation of the expectation value with a frequency directly given by h.
In all of the cases mentioned above there is a non-local transformation which actually connects the deformed model with its undeformed parent model. Then the additional phase in (224) can be interpreted as the additional effect of certain displacements dictated by the non-local transformations. These displacements can be interpreted as the particles having a finite width. Thus we can interpret these models and their non-local transformations as concrete and explicit examples of the hard rod deformation discussed [28,29].
It is known that in QFT the actual TT -deformation can be also be interpreted by non-local and fielddependent transformations, and that the particles can be seen as having acquired a (momentum dependent) finite width [28,29,97,98]. The examples discussed so far demonstrate such a relation in lattice systems.
Regarding the full spectrum of the present model the situation is somewhat more complicated than in the constrained XXZ model or the phase model. Now the ultimate non-local transformation that maps the interacting model to an (almost) free model concerns the Maassarani-Mathieu chain, which is a spin chain with a different local Hilbert space. Thus we have found a new type of non-local map, which nevertheless belongs to the class of hard rod deformations.

XIV. CONCLUSIONS
In this work we studied the so-called folded XXZ model which appeared in the earlier works [50][51][52]. The main new results of our work are a) the direct connection with the charges of the XXZ model; b) the various forms of the Bethe Ansatz, both with periodic and open boundary conditions; c) the three different descriptions of the thermodynamics; d) the connection with the Maassarani-Mathieu chain; e) the exact solution for a specific quantum quench problem; f) a conjecture for the presence of persistent oscillations, supported by numerical proof; g) making a connection to the TT and hard rod deformations, and also to similar models in the literature.
There are a number of interesting open questions. First of all, it would be interesting to understand the direct algebraic origin of the present model, and how it fits into the Quantum Inverse Scattering Approach [99]. Our derivation is capable of producing the charges and the eigenstates of the model, nevertheless a direct treatment with the Algebraic Bethe Ansatz is not possible. The reason for this is that the Lax operators which are used in the construction of the transfer matrix of the XXZ model do not survive the ∆ → ∞ limit, or at least this does not happen in a straightforward way. In fact we can prove that the model can not be constructed using the standard way, where the Lax operator is chosen as the Rmatrix in the fundamental representation: in such cases the charges can be constructed using the boost operator, whereas in our model Q 2 is not dynamical and a direct application of the boost will not produce the next charge Q 3 ; this was discussed in Section IV. Thus the model requires a different Lax operator construction, and a direct ∆ → ∞ of the Lax operators of the XXZ chain will not work.
It would be interesting to consider other models with constant scattering lengths, and hard rod-like S-matrix factors, and to formulate the common algebraic origin of these models. So far the known models in this family are the phase model (the q → ∞ limit of the q-boson model), the Rule 54 model, and the present one, the folded XXZ model. It would be desirable to develop a unified algebraic framework for such spin chain models. As a particular case of this problem, it would be important to uncover relations between the Rule 54 model and the quantum chains that we treated.
The exact results in these relatively simple models could be used to confirm the predictions of Generalized Hydrodynamics (GHD) in various situations. Such a highly non-trivial check was performed in the recent work [100] in the case of the Rule 54 model. We believe that the present spin chain and the related models are perfect candidates for such checks: they are quantum models that have genuine interactions, nevertheless they are simple enough that even the spectral sums can be evaluated in certain situations. This presents unique opportunities for exact computations.