The Folded Spin-1/2 XXZ Model: I. Diagonalisation, Jamming, and Ground State Properties

We study an effective Hamiltonian generating time evolution of states on intermediate time scales in the strong-coupling limit of the spin-1/2 XXZ model. To leading order, it describes an integrable model with local interactions. We solve it completely by means of a coordinate Bethe Ansatz that manifestly breaks the translational symmetry. We demonstrate the existence of exponentially many jammed states and estimate their stability under the leading correction to the effective Hamiltonian. Some ground state properties of the model are discussed.


Introduction
Thermalisation in isolated quantum systems has received enormous attention in the last few decades (see review articles [1][2][3] and references therein), partly due to the advancement of the experimental techniques that allow one to manipulate a large number of quantum particles [4,5]. A particularly important question in this regard has been of how the conserved quantities affect the dynamics of a system. In integrable models this was clarified only recently with the introduction of so-called Generalised Gibbs Ensembles (GGE) [6][7][8], which incorporate the constraints of infinitely many conserved charges [9,10] and describe the expectation values of the local observables at asymptotically large times [11]. In the attempt to apply a similar theory in the presence of inhomogeneities, a generalised hydrodynamic description (GHD) was developed [12,13]. To leading order, which is the one better understood so far, the theory describes the time evolution of a class of inhomogeneous states by means of a "classical" continuity equation 1 , which, together with Bethe Ansatz, predicts the expectation values of the local observables on Euler scales. On the other hand, the dynamics of isolated integrable systems on intermediate time scales remains an intriguing open problem. Strong coupling limits in this respect play an important role. In particular, unusually slow symmetry restoration typically indicates that the system can be seen as a perturbation of an effective one, in which the symmetry is never restored. Knowledge of the effective system and its symmetries is thus necessary to understand pre-relaxation behaviour [15]. Strong coupling limits of quantum models have been thoroughly investigated at equilibrium, especially after the seminal work by MacDonald, Girvin, and Yoshioka [16], which proposed an asymptotic perturbation theory for a class of systems that includes the Hubbard model. Later, an alternative approach to calculate the thermal correlation functions in strong coupling limits was developed; it is based on investigating the spectrum of constrained Hamiltonians describing particles with a finite radius and acting on the reduced Hilbert space [17][18][19][20]. Strong coupling descriptions of time evolution are, however, still under investigation. Based on the results of MacDonald et. al., we consider an asymptotic formulation of quantum mechanics that attaches the fast oscillatory part of the dynamics to the operators, letting the state evolve gently in time under an effective Hamiltonian. We term this formulation "asymptotic folded picture", to emphasise that the spectrum of the effective "folded Hamiltonian" is the same as that of the original Hamiltonian, modulo a typical energy proportional to the largest coupling constant. In the folded picture localised operators are asymptotically mapped to quasilocalised ones, hence the effective Hamiltonian captures the essential features of the original dynamics.
As an example, we consider the large anisotropy limit of the spin-1 2 XXZ chain and present a Bethe Ansatz solution of the model described by the corresponding local folded Hamiltonian. Since the XXZ model is integrable for any value of the anisotropy parameter, one could address this problem by applying the Schrieffer-Wolff transformation of Ref. [16] to the eigenstates of the original model -see also Ref. [21]. This approach, however, conceals the extra symmetries that generally emerge in strong coupling limits. For example, at large anisotropy the time evolution generated by the effective Hamiltonian of the XXZ model does not restore one-site shift invariance if that is broken in the initial state [15]. Being translationally invariant, the standard Bethe Ansatz basis is then arguably inappropriate for investigations into intermediate time scales. We propose to overcome this problem by a change of basis that explicitly breaks one-site shift translational symmetry. We stress that the choice of a basis is important because the theories describing relaxation are usually presented in a way that is, in fact, basis dependent; hidden symmetries can spoil the predictions, despite the theories being correct -see, e.g., Ref. [15]. The dissimilarity between our solution and the standard one is manifested in the striking fact that the total spin along the anisotropy axis, which counts the number of rapidities in the standard Bethe Ansatz solution, is not even diagonal in the basis that we consider.
We will show that the folded XXZ Hamiltonian exhibits a richer structure than the XXZ Hamiltonian, reflecting the non-abelian integrability of the folded model. In particular, the Hamiltonian can be diagonalised by introducing two species of impenetrable particles that interact via a diagonal scattering matrix. This implies the conservation of the particle configuration, which turns out to have visible effects. The Hilbert space divides into three exponentially large sectors: (a), a fully interacting one, with a mixed population of both types of particles, (b), two noninteracting sectors, where only a single particle species is present, and, (c), a sector of jammed inhomogeneous eigenstates in which the particles are densely packed and can not move. Remarkably, a subset of the latter states is stable under the perturbation by the leading correction to the strong coupling limit.

Summary
Section 2. We work out an iterative solution of the strong coupling expansion firstly proposed in Ref. [16] and define the "folded picture" and the "folded Hamiltonian". We also exhibit the first orders of the formal asymptotic expansion of the folded Hamiltonian and the corresponding unitary transformation.
Section 3. We focus on the folded Hamiltonian of the XXZ model. We identify a duality transformation that compresses the range of the operator in the limit of infinite anisotropy.
Section 4. We define two species of particles that sit on macrosites consisting of two neighbouring sites and show that their configuration is a topological invariant. We exhibit an exponentially large number of jammed product states with zero energy and assess their stability under the leading strong-coupling correction.
Section 5. The (dual) folded XXZ Hamiltonian is diagonalised via a coordinate Bethe Ansatz based on the two species of particles. Section 6. We identify the ground state of the folded XXZ Hamiltonian and compute its energy. At a given energy density we also compute the maximal number of consecutive particles of the same species, i.e., the size of a domain.
2 Dynamics in strong coupling limits: the folded picture A large class of quantum lattice models are described by Hamiltonians H(κ) that can be decomposed as [16,[23][24][25][26][27][28][29][30]] where q is a finite integer, κ the inverse coupling constant (a perturbative parameter in the strong coupling limit), while operators H F , H I , and F m , for m = 1, 2 . . . , q, satisfy Constant J has the units of energy and typically corresponds to an overall factor in the Hamiltonian H(κ). Algebraic relations (2) imply that, if φ⟩ and φ ′ ⟩ are simultaneous eigenstates of H I and H F , either ⟨φ F m φ ′ ⟩ = 0 or ⟨φ H I φ⟩−⟨φ ′ H I φ ′ ⟩ = mJ. Thus, the eigenstates of H(κ) can be chosen to belong to sectors where the spectrum of H I is equally spaced, the eigenvalue differences being multiples of J. This is the basic reason why, in all the relevant examples, H I turns out to be a very simple operator, for example, the projection of the total spin in a given direction. For the sake of simplicity, we will assume that there is a single sector, i.e., H I has equally spaced eigenvalues and can be interpreted as a "number operator". The class of models that allow for decomposition (1) includes both integrable and generic systems. Examples with q = 1 are the Hubbard model [16], the Heisenberg XYZ model, the Ising model in transverse and longitudinal magnetic field, the Kitaev model [24], the ANNNI model, and the Bloch-Stark MBL Hamiltonian [25,26]. While the focus of this paper is on the Heisenberg XXZ model, in which q = 1, this section aims to present the framework for a generic q.
To the best of our knowledge, Ref. [16] was the first to observe that, in the strong coupling limit κ ≪ 1 with q = 1, algebra (2) can be exploited within a perturbation theory to construct a weak coupling model with the spectrum of the original Hamiltonian. In this section we propose a formulation of quantum mechanics that exploits the results of Ref. [16] to investigate time evolution on an intermediate time scale. Specifically, we define the asymptotic folded picture as the formulation in which operators O and state Ψ(t)⟩ time evolve according to where z t = e iJt κ and the following conditions hold: 3. B z (κ) is (assumed to be) holomorphic in z, in an open annulus enclosing the complex unit circle S 1 = {z ∈ C; z = 1}, and its Laurent series is Within this picture, the part of the dynamics that oscillates rapidly with a period 2πκ J is applied to the operators, whereas the states time evolve under a weak-coupling Hamiltonian that conserves the strong coupling term H I . Since H I is generally trivial, this allows for a non-perturbative definition of a reference state (for example, the ground state of H I ), over which we can define and organise elementary excitations. Using algebra (2) and the conditions satisfied by B z (κ), we see that formulation (3) of the asymptotic folded picture is equivalent to expressing the time evolution operator as Hence, we can identify H F (κ) + κ −1 H I with the effective Hamiltonian equivalent to H(κ), and e iκB 1 (κ) with the corresponding unitary transformation, both considered in Ref. [16]. As an interesting fact we observe that the unitary transformation (5) is equivalent to the end and subject to the initial condition Φ 0 = H(κ). Roughly speaking, operator B s (κ) acts as a generator of the flow equations which have been introduced by Wegner [22] as a way to suppress the off-diagonal terms of a Hamiltonian by means of a continuous unitary transformation. In our case, however, the ending point H F (κ) + κ −1 H I of the flow is far from being a diagonal operator. In fact, it describes a generic interacting model. Both the folded Hamiltonian H F (κ) and the unitary transformation e iκB 1 (κ) can be worked out asymptotically. In Appendix A we report the recurrence relation for the coefficients of the asymptotic expansion as well as the lowest orders of the expansion for q = 2 (for q = 1 see also Ref. [27]). Here we will focus on the zeroth order and the leading correction, which yield for the folded Hamiltonian (the κ 2 correction is reported in Eq. (156) of Appendix A), and for the unitary transformation. Finally, we note that [H I , H F (κ)] = 0 and Eq. (5) imply that e −iκB 1 (κ) H F (κ)e iκB 1 (κ) is a conservation law for H(κ); by checking the first orders of its asymptotic expansion using Eqs. (9) and (10), we realise that it is exactly the charge perturbatively constructed in Ref. [31].

Integrability
One of the interesting aspects of the asymptotic expansion (9) is that the leading order H F of the folded Hamiltonian H F (κ) can be integrable even if the full H(κ) pertains to a generic system. Moreover, since H F (κ) commutes with H I , which is generally a trivial operator, a simple reference state can always be defined and the asymptotic model is likely to be solvable by coordinate Bethe Ansatz. Remarkably, the leading order of the folded Hamiltonian has usually more symmetries than the original model, starting from the U (1) symmetry that is associated with H I by construction. In addition, every partial symmetry acting like UH(κ)U † = H(−κ) becomes an exact symmetry for H F .
Especially interesting phenomenology arises when the folded Hamiltonian is asymptotically non-abelian integrable, that is, when it possesses a non-abelian set of quasilocal conservation laws. The practical effect of non-abelian integrability is that states sharing a complete set of integrals of motion remain locally distinguishable even at late times [15]. This is manifested in a physically relevant degeneracy of stationary states even at extensively high energies, in turn implying that different integrable models share the leading order of the folded Hamiltonian. In this situation, constructing the folded Hamiltonian as a limit of an integrable system is likely to hide the emergent symmetries. It is arguably more effective to consider the folded Hamiltonian as a different model, redefining, if advantageous, reference state and excitations. This is the point of view that we embrace in the next part of the paper, where we work out a Bethe Ansatz solution to the asymptotics of the folded XXZ Hamiltonian.
On the other hand, even if the original model is integrable for any value of κ, it is not reasonable to expect the folded Hamiltonian truncated at a given order to remain integrable. However, if the truncation still has more symmetries than the original Hamiltonian, we envisage the possibility to find distinct integrable systems sharing the first orders of the expansion. In that case, the dynamics of states in the strong coupling regime can exhibit a cascade of pre-relaxation plateaux.

Examples
Generic model with integrable H F (0). An example of a generic system with an asymptotically integrable folded Hamiltonian is the XYZ spin-1 2 chain in an external magnetic field. It is described by the Hamiltonian In the strong coupling limit h → ∞, the folded picture is characterised by the XXZ Hamiltonian which is integrable. Identifying h with κ −1 we find q = 1 and so that the leading correction reads To the best of our knowledge this perturbation breaks integrability, however, we have not investigated whether integrability can be restored by the addition of the subleading terms.
Generic model with non-abelian integrable H F (0). An example of a generic model with a non-abelian integrable strong coupling limit is given by the quantum Ising model in longitudinal and transverse fields, described by the Hamiltonian In the limit h → ∞ we identify κ = h −1 , q = 2, and The expansion in h −1 of the folded Hamiltonian reads where K n,m = σ x n σ x m + σ y n σ y m . At the leading order H F describes the XX model (first term), which is noninteracting and known to be non-abelian integrable [15]. Quite interestingly, the folded Hamiltonian remains noninteracting even at O(h −1 ) (second term). The O(h −1 ) correction, however, breaks the non-abelian integrability. Interactions appear only at O(h −2 ) (second line).
Integrable model with quasilocal H F (κ). Strong coupling limits in noninteracting models generally result in quasilocal folded Hamiltonians. We mention, for instance, the quantum Ising model in a transverse field which corresponds to setting g = 0 in the previous example. Since this Hamiltonian can be represented by a quadratic operator of spinless fermions, it is easy to find a unitary transformation mapping H into a U (1) symmetric operator that commutes with ∑ σ z -see Appendix B. Specifically, we have with , and {a , a n } = 2δ n 1.
A posteriori we see that the expansion in h −1 converges for h −1 < h c = 1, where h c is the magnetic field corresponding to the Ising phase transition. In particular, H F (h −1 ) and B 1 (h −1 ) are quasilocal with typical range (log h) −1 .
Interacting H F (0) that breaks one-site shift invariance. In this class we identify the XYZ spin-1 2 chain described by the Hamiltonian (11) with h = 0. In the strong coupling limit ∆ → ∞, the folded picture is characterised by the Hamiltonian The reader can easily check that H F commutes with which breaks one-site shift invariance and causes the system to retain the memory of one-site shift asymmetry. In the rest of the paper we will focus on the special case γ = 0, which exhibits even more symmetries that give rise to exceptional degeneracies and jamming. We leave the study of Hamiltonian (22) with γ ≠ 0 to future investigations.
3 The folded XXZ Hamiltonian in the strong coupling limit The folded Hamiltonians listed in the previous section are all asymptotically integrable, and, except for the last one, have leading terms H F that describe well-known models. Here, we focus on the folded Hamiltonian of the XXZ model at large anisotropy, which turns out to be a meagrely studied integrable system with striking features. The XXZ spin-1 2 chain is described by the Hamiltonian where periodic boundary conditions σ α L+1 = σ α 1 are understood. It is arguably the simplest model capturing the antiferromagnetic properties of crystals with effective one-dimensional magnetic behaviour, like the KCuF 3 compound [32]. From a theoretical point of view, it is integrable for any value of the anisotropy ∆. The ground state phase diagram exhibits antiferromagnetism for ∆ > 1, ferromagnetism for ∆ < 1, and a critical phase for ∆ ∈ [−1, 1]. We are interested in the limit of large anisotropy. In the folded picture we identify κ with (4∆) −1 , and the operators H I , F 1 ≡ F, and H F are given by To the best of our knowledge, until now the folded Hamiltonian has been considered only in the combination H F +4∆H I in the limit of large ∆. That is, indeed, the Hamiltonian unitarily equivalent to the XXZ model in the limit of large anisotropy. Since H I commutes with H F , the two operators are diagonal in the same basis. As discussed in Ref. [21], H F could thus be diagonalised by transforming the eigenstates of the XXZ Hamiltonian and expanding the corresponding eigenvalues in the limit of large ∆. The drawback to this approach is that the solution is obtained as a limit of a model with less symmetries, forcing one to deal with a redundantly complicated structure (analogously to solving the quantum harmonic oscillator as a limit of the anharmonic one). In order to circumvent these complications, we lift the folded XXZ Hamiltonian into an independent model, providing an ab initio solution that will allow us to unveil some of the properties that make it essentially different from the XXZ model.

Symmetries and short-range conservation laws
The folded Hamiltonian inherits the following symmetries of the XXZ model: one-site shift invariance, reflection symmetry, spin flip symmetry ([∏ L =1 σ α , H] = 0, for α = x, y, z) and conservation of the total magnetisation along the z-axis By construction, it also commutes with H I (25). The first unusual symmetry of H F is associated with the conservation of the staggered version of H I , which breaks one-site shift invariance. This symmetry is, in fact, broken at a higher level.
In particular, Ref. [15] showed that the model has noninteracting sectors (characterised by σ z 2 −1 σ z 2 = −1) that break translational invariance and can be effectively described by a noninteracting Hamiltonian with a non-abelian set of local conservation laws. The origin of these sectors will be clarified in Section 4.1.
The model described by the Hamiltonian (27) has other short-range shift-invariant conservation laws. One can identify them by imposing the vanishing of the commutator between a generic local translationally invariant operator and H F . We denote them by where q ± n, is the local density acting on sites , . . . , + 2n + 1±1 2 , n ∈ Z > ≡ N, while the sign specifies whether the charge is even (+) or odd (−) under reflection symmetry. The densities of the first few conservation laws read where K n,m is defined below Eq. (17) and D n,m = σ x n σ y m − σ y n σ x m is a Dzyaloshinskii-Moriya term. Note that Q + 1 = H F is exactly the folded Hamiltonian. Interestingly, the folded Hamiltonian H F and other computed short-range charges Q ± n commute with the transfer matrix of the four-vertex model, obtained from the XXZ model's transfer matrix in the ∆ → ∞ (Ising) limit [18,19]. The four-vertex model transfer matrix enables the algebraic Bethe Ansatz technique, which produces a basis of eigenstates corresponding to the nonzero energies of the following two equivalent projected Hamiltonians: Here, H XX = J ∑ L =1 K , +1 is the Hamiltonian of the XX model, while the projector prevents neighbouring sites to be occupied simultaneously, thus inducing a hard-core repulsion [18][19][20]. Because of the latter, such projected Hamiltonians can only reproduce the lowenergy sector of H F . 3 The eigenstates that are destroyed by the projector P unfortunately turn out to be inaccessible by the algebraic Bethe Ansatz that diagonalises the four-vertex transfer matrix. Nevertheless, our numerical analysis is consistent with the existence of an infinite sequence of local charges, resembling the family of local conservation laws of the XXZ Hamiltonian. The latter set is known not to characterise the excited states completely: the remaining conservation laws are quasilocal and are, within the framework of algebraic Bethe Ansatz, generated by the transfer matrices associated with higher-spin representations of the Lax operator [33]. We will come back to the completeness of the local charges of the folded XXZ Hamiltonian after having worked out a coordinate Bethe Ansatz that produces a complete basis of eigenstates.

Duality transformation
In this section we exhibit and work out a duality transformation that, up to boundary terms, maps H I into J 2 S z and reduces the range of the local density of H F from four to three neighbouring sites. The transformation consists of two steps 4 : 1. rotation around the y-axis 2. left antiperiodic pseudo-site shift of the Jordan-Wigner Majorana fermions (21) We will callÕ the image of operator O under this transformation. In particular, it maps the folded Hamiltonian with periodic boundary conditions (σ α L+n = σ α n , for α = x, y, z) intõ where we introduced the notation with Π z = ∏ L =1 σ z . PBC and aPBC refer to periodic and antiperiodic boundary conditions, respectively (σ x,y L+n = ±σ x,y n ). The operators F and H I are instead mapped intõ Regarding the other local charges, we conjecture that the transformation O ↦Õ reduces the ranges of local densities q + n, by one site and leaves the ranges of q − n, intact. Thus,q ± n, both act on the same number of sites. Numerical calculation confirms our conjecture for charges with ranges of up to six sites. We also note that 2 J H − I , given in Eq. (29), is mapped into the staggered spin along the z-axis To diagonalise the dual folded HamiltonianH F , it is now enough to consider the operator 4 Combining the two steps, the transformation reads which commutes with the total magnetisation S z along the z-axis 5 . Specifically, we have the following correspondence between the eigenstates ofH F and those ofH η F : Considering this mapping in the opposite direction we note that, in both cases, only the eigenstates ofH η F that belong to the sector Π z = 1 are needed to diagonaliseH F . In particular, for η = 1 this is due to the additional action of σ x L that flips the spin on the last site of the chain and therefore maps the state from one sector to the other. In light of correspondence (41) we will restrict our considerations toH η F in the rest of this work.

Dual XXZ Hamiltonian
For the sake of completeness, we also report the dual XXZ Hamiltonian, which reads This operator commutes with Π z , indeed each of its terms flips two next-nearest neighbour spins simultaneously. On the other hand, the total magnetisation in the direction of the anisotropy is mapped into the operator which anti-commutes with Π z . As a consequence, the eigenstates ofH in the sectors Π z = ±1 are also eigenstates ofS z only if they are in the kernel of the latter operator. Otherwise, they are degenerate between the two sectors. The same holds for the eigenstates ofH F , and this is arguably the biggest difference between the standard basis of eigenstates and the basis that is implicitly chosen when diagonalising the folded Hamiltonian throughH η F (or the XXZ Hamiltonian through Eq. (42)).
The action ofH η F in the standard spin basis Since the total magnetisation along the z-axis S z = 1 2 ∑ σ z commutes withH η F , there are bases in which the eigenstates ofH η F are also eigenstates of S z . We denote by s⟩ = s 1 , . . . , s L ⟩ the eigenstate of the projections of all local spins on the z-axis, so that σ z s⟩ = s s⟩, for = 1, 2, . . . , L. Understanding the action ofH η F on such states will later be necessary for the identification of the Bethe states that diagonalise the model and the conservation laws. We haveH whereH η F, denotes the local density ofH η F , acting on sites , + 1 and + 2, with boundary condition σ x L+n = (−1) η σ x n . The local action ofH η F amounts tõ where ↓ stands for s = −1 and ↑ for s = 1. The other combinations of three neighbouring spins are destroyed. In words, two oppositely aligned next-nearest neighbour spins are exchanged if the spin between them points downwards.

Invariants and jamming 4.1 Elementary particles
Given that the local densityH η F, ofH η F is supported on three neighbouring sites labeled by , +1, +2 (see, for example, Eq. (40)), the interaction becomes of the nearest-neighbour type if we define macrosites that consist of two neighbouring lattice sites. We use the following notation: (46) Clearly there is a local ambiguity in the position of the macrosite, which could be defined to start either from an odd or from an even site. If the chain has an even number L of sites, this ambiguity is irrelevant, as it corresponds to fixing a convention. If instead L is odd, the two options are linked by the boundary. We will see later how this can be taken into account; for the moment the reader can assume that the chain has an even number of sites. The macrosite will be denoted by a primed index ′ and consists of the neighbouring sites 2 ′ − 1 and 2 ′ .
The dynamics generated byH η F can be summarised by the following actions of local densities: The remaining (orthogonal) cases of occupations of the macrosites ′ and ′ + 1 are sent to zero by the local densitiesH η F,2 ′ −1 andH η F,2 ′ . In particular, we see that the state with all spins down (all macrosites are occupied by ∅) is an eigenstate ofH η F with zero eigenvalue. The first four rules in Eq. (47) then suggest to identify such a state with the vacuum of two species of particles, ⊕ and ⊖, that move freely until another particle is met; their interaction is then encoded in the remaining four rules of Eq. (47). Within this picture, the state ⊖⊕ can be interpreted as a particle ⊖ and a particle ⊕ sharing the same macrosite. Note that a particle of one type cannot jump across a particle of the other, thus the domains of particles ⊕ and ⊖ do not change under time evolution. Going back to the folded Hamiltonian, we note that the vacuum state φ⟩ is associated with an eigenstate ofH F only if the size L of the system is even (condition Π z φ⟩ = 1 should be fulfilled -see Eq. (41)).

Noninteracting sectors
Let us consider a state with a single species of particles, say ⊕. The entire dynamics is described by the first two rules in Eq. (47). If L is even, the mapping ⊕ ↦↑, ∅ ↦↓ sendsH η F into the Hamiltonian of the XX model with L 2 spins, This Hamiltonian describes a noninteracting spin chain. For (−1) η = (−1) L 2 , H η XX possesses a non-abelian set of local conservation laws. Since both boundary conditions are compatible with the existence of non-commuting charges, this sector of H F has the phenomenology typical of non-abelian integrable systems, as first discussed in Ref. [15].
For L odd, the hypothesis of having a single species of particles is meaningless: the two species are transformed one into the other when crossing the boundary. The closest situation is having no more than two groups of particles of different species. Such configurations would still form an invariant subspace, like in the even case, but now in presence of interaction. We postpone the analysis of such a case to the next sections.

Topological invariants
If we indicate by ′ j the macrosite of the j-th particle, the actual position of the corresponding spin up is j = 2 ′ j − b j , where b j = 0 for a particle of type ⊕ and b j = 1 for a particle of type ⊖. Assuming L to be even, rules (47) imply that the set B  Figure 2: The same as in Fig. 1 for an odd chain with L = 7 sites; we report all the configurations analogous to those in Fig. 1 (there are more options because the sequences in Fig. 1 are defined up to cyclic permutations). The particle in parentheses is transmuted due to the periodic boundary conditions (dashed vertical line). Note that in odd chains each configuration represents a space that is invariant under a shift by one site. We point out that the states in (a 1 ), (a 2 ), (b 1 ) and (b 2 ) are jammed (indeed the particles can not move), so each one separately generates an invariant subspace.
other -see Fig. 1. This invariance is not present in the original XXZ Hamiltonian, so it is a candidate for generating pre-relaxational behaviour at large anisotropy. Even if L is odd a similar result holds, however, in order to account for the transmutation of particles at the boundary, the configuration of particles must be encoded in a larger set of sequences, namely B (o) Fig. 2. 6 As a result, contrary to the even-size case, the space associated with a configuration is invariant under a translation by one site.
In general, the configuration is a "topological invariant" in the sense that it is independent of the actual positions of the particles and hence does not change when the particles are arbitrarily moved without crossing.
We point out that all functionals of the set B (e) , which has eigenvalues ∑ N j=1 b j , such charges are not local, and it is not evident whether they could be defined as quasilocal [33], i.e., as translationally invariant sums of local densities that are exponentially suppressed in their range 7 . Assuming L to be even, an important example is provided by the following charge whose eigenspaces are characterised by the configurations B These eigenvalues count the number of oriented domain walls separating particles of type ⊖ on the left-hand side from those of type ⊕ on the right-hand side. We point out that M is quasilocal, in particular, it is extensive, in the sense that where ⟨•⟩ denotes the expectation value in some physically relevant state with a finite correlation length, e.g., a thermal state, or a generalised Gibbs ensemble [34]. Such states will be considered in the second part of our work [35], where we elaborate on the thermodynamic limit of the Bethe Ansatz. The importance of M will become evident in Section 5, here we only mention its appearance in a necessary and sufficient condition for the existence of the configuration B (e) N : . This is a consequence of the fact that generally each particle corresponds to a pair ↑↓ or ↓↑, for a total of 2N sites; exceptions, however, can arise on the M oriented domain 6 Note that the configuration (b1, . . . , b N ; 1 − b1, . . . , 1 − b N ) describes the particle content of two side-by-side copies of the spin configuration, i.e., (s1, . . . , s N ; s1, . . . , s N ). 7 The range of a localised operator is the number of sites on which the operator acts nontrivially. By exponential localisation we mean that the local density can be written as a sum of operators that have distinct ranges, and whose fluctuations in a generic macro-state decay exponentially with their range. For example, in equilibrium at infinite temperature the fluctuation of a projector Πn acting on n adjacent sites reads walls, where ⊖ and ⊕ can occupy a single macrosite (we then denote them by ⊖⊕), releasing 2M sites. Incidentally, we note that the average distance between two particles, measured in the number of macrosites, is ξ = L 2N . Since L ≥ 2(N − M ), ξ can not be smaller than 1 − µ, where µ = M N parametrises the average degree of macro-degeneracy of particles. The configurations with minimal average distance between particles, namely, the ones with L = 2(N − M ), are special and will be investigated in detail in the next section.

Jammed states
In this section we study the subspace with L = 2(N −M ), which is the ground-state eigenspace of the quasilocal charge M − S z . Without making assumptions on the behaviour of the states under translations, we exhibit an orthonormal set of stationary states spanning the entire subspace. We nonetheless anticipate that the Bethe Ansatz solution that will be worked out in Section 5 is complete, including also such states in the form of linear combinations invariant under shifts by two sites.
In the states with L = 2(N − M ) a spin down is always surrounded by spins up, and from Eq. (45) it follows that any product state Ψ⟩ with that property belongs to the kernel ofH η F . Such a condition can be expressed as that is to say, Ψ⟩ is in the eigenspace of the ground state of the classical Ising model Each product state corresponds to a configuration where the particles (⊕ and ⊖) are jammed -see Fig. 2; for that reason we call them jammed states. Remarkably, translational invariance can be broken in a chaotic way, as the jammed states consist of strings of spins up with arbitrary length, interspersed with single spins down. As a result, jammed states can play a key role in nonequilibrium time evolution, preventing restoration of one-or two-site shift invariance when the initial state is not symmetric. The connection with the classical Ising model allows us to easily compute the size of the space which is, remarkably, exponentially large. 8 The number of jammed states at fixed number N of spins up is instead given by The first term counts the product states in which the left-most position is occupied by a spin up. Since two neighbouring spins cannot point downwards, spins down should always be paired with a spin up. We therefore have to distribute N "fragments" of consecutive spins, of which L − N are of type ↑↓ and the rest of type ↑. In contrast, the second term in Eq. (56) counts the configurations in which the left-most position is occupied by a spin down and the right-most one by a spin up (the latter can not be down because of periodic boundary conditions). We have then to distribute N − 1 fragments, of which L − N − 1 are of type ↑↓, the rest being spins up. Incidentally, we note that, since every spin down comes in a pair with a spin up, the total magnetisation of a jammed product state is always non-negative, i.e., S z = ⟨S z ⟩ = N − L 2 ≥ 0. Transforming back to the folded Hamiltonian (27), only the jammed states with an even number of spins down, i.e., those belonging to the sector Π z = 1, correspond to eigenstates of H F . Indeed, according to Eq. (41), if s⟩ is a jammed state ofH as one can easily check by applying the inverse of the duality transformation, mappingH F into H F : Finally, a symmetric and antisymmetric combination of the states (57) with η = 0 and η = 1 produces the jammed product (eigen)states of H F , namely ±1, ±s 1 , ±s 1 s 2 , . . . , ±s 1 ⋯s L−1 ⟩ .
These are the product states in which all the domains consisting of spins aligned in the same direction have length larger than 1. This follows from Eq. (53), which prohibits two consecutive s (in the dual basis) from being equal to −1. Interestingly, cluster decomposition properties, which were absent in the state shown in Eq. (57), are restored by mixing jammed states belonging to different sectors of the dual folded HamiltonianH F .

Jammed states under the leading correction
The leading asymptotic correction to the dual folded Hamiltonian readsH , whereF is given in Eq. (38) (see also Eq. (9)). Explicitly, we havẽ which can be rewritten asH are the dual conserved charges H I and Q + 2 , whose densities are given in Eqs. (25) and (31), respectively.
The leading correction spoils the conservation of the configuration and, in turn, the jamming, as evident in the following actions of its local density: The subset of states Ψ⟩ that remain jammed even when the leading correctionH ′ F is included has an additional property: no substrings of the form ↓↑↓ are present in the spin configuration of Ψ⟩. This can be compactly expressed as follows: Since each spin down has to be surrounded by distinct spins up, the minimal number of spins up is ⌈ 2 3 L⌉; the total magnetisation is therefore bounded from below by S z ≥ ⌈ 1 6 L⌉. We note that the energy of the jammed states is not zero anymore -it is now given by Since all the energies are multiples of 4Jκ, the dynamics in the subspace of the jammed states is periodic with period π (2Jκ). For a given N , and hence energy, the number of jammed product states reads To see this, let us use that a spin down should always be surrounded by spins up, while a spin up can stand alone. The first term accounts for the configurations starting and ending with a spin up. It counts the number of ways in which L − N fragments of type ↑↓↑ can be distributed among the stand-alone spins up, ↑. The total number of fragments ↑↓↑ and ↑ to distribute is 2N − L. Conversely, the second term counts the configurations that either start or end with a spin down. There, one has to count the number of distributions of L − N − 1 fragments of type ↑↓↑ among the stand-alone spins up, fixing ↓↑ on the first two sites and ↑ on the last site (or vice versa, whence the prefactor 2). As in the absence of the leading correction, fixing the boundary spins is a consequence of periodic boundary conditions and of the constraint that prohibits substrings of the form ↓↑↓. The total number of fragments ↑↓↑ and ↑ that have to be distributed is now 2N − L − 1. The total number of jammed states is obtained by summing Eq. (64) over all allowed N : The (leading) asymptotic approximation χ L , where χ < 1 is the real solution to χ 2 (χ − 1) = 1, arises from one of the terms in the first sum. Again, the kernel of the folded Hamiltonian contains an exponentially large number of jammed states, but its size is negligible with respect to the one without the leading correction -cf. Eq. (55). Transforming back to the folded XXZ Hamiltonian, the additional constraint means that the domains consisting of spins aligned in the same direction must contain at least three neighbouring sites.
We point out that, in order to apply this result to the XXZ model in the limit of large anisotropy, a further step should be made: when the leading correction is considered, the second formula of Eq. (3) implies that the actual asymptotically jammed states Ψ(0)⟩ are obtained by acting on the jammed states of H F (κ) with operator e −iκB 1 (κ) = e (F † −F) (4J∆)+O(∆ −2 ) -see Eq. (10). They read where s satisfy

Coordinate Bethe Ansatz
In this section we diagonalise the asymptotic dual folded HamiltonianH F , defined in Section 3.2. Remarkably, the solution that we present is independent of the standard Bethe Ansatz for the XXZ model. This is possible due to three conditions that arise in the strong coupling limit: 1. other reference states can be used as vacua of the Bethe Ansatz description; 2. the particle structure of the XXZ model becomes redundant, in the sense that particles can be reorganised in more efficient ways; 3. the spectrum becomes highly degenerate, and inequivalent bases can be chosen.
The difference with the standard Bethe Ansatz can be recognised a priori. Firstly, the total magnetisation S z , which counts the number of rapidities associated with the standard Bethe Ansatz eigenstates, is not diagonal in our basis of eigenstates of H F . This is because its dualS z does not preserve the sectors that block diagonalise the dual folded Hamiltonian -see Eqs. (37) and (41), and the discussion after Eq. (43). Secondly, since the particle configuration is conserved, the asymptotic dual folded HamiltonianH F can be diagonalised at fixed configuration. This is an additional structure that is not present in the coordinate Bethe Ansatz of the XXZ model, where one only exploits the conservation of the particle number.

Even number of sites
The eigenstates corresponding to the set B in which each term has the same relative spatial distribution of particles, enforced by the in the sum. The coefficients of the terms that describe a particular spatial configuration b of particles read Here, S b 1 ,...,b N p 1 ,...,p N p π(1) ,...,p π(N ) is the multi-particle scattering matrix that describes the permutation (p 1 , . . . , p N ) ↦ (p π(1) , . . . , p π(N ) ) of momenta due to two-particle scattering processes, which leave the relative spatial distribution of particles intact, exchanging only their momenta. In particular, we assume i.e., the amplitude of the process without collisions between particles is 1, and each process is reversible. We remind the reader that the exceptionality of this Ansatz is in preserving the set of momenta under the scattering process. The boundary conditions σ x L+n = (−1) η σ x n imply i.e., a prefactor (−1) η appears each time a particle is moved across the boundary from the first to the last position in the configuration. We then find for a generic permutation π. Note that any ratio Z satisfies Eq. (72), provided that we shift the configuration b in the scattering matrix accordingly (the indices of the momenta should, however, remain the same).
The Bethe equations are now a result of enforcing, (i), the π-independence of Eq. (72), and, (ii), the conservation of the total momentum. The latter comes from imposing on the ratios exhibited in Eq. (72) for the sequence of permutations π(j) = j in the first ratio, π(j) = j + 1 in the second, and so on and so forth. The scattering matrices then cancel due to the reversibility of the scattering process (right-hand side of Eq. (70)), and we obtain  ∶ (j, 1, . . . , N ) ↦ (1, 2, . . . , N ). The dots represent particles of types b 1 , . . . , b N , to which momenta are assigned. The lines connecting them represent the exchange of momenta during collisions. Each crossing of two lines contributes a two-particle scattering matrix to the product that, at the end, represents the entire scattering process.
The condition that Eq. (72) should not depend on the choice of permutation can be used to obtain an additional identity. Choosing, for instance, permutation π ∶ (1, 2, . . . , N ) ↦ (j, 1, . . . , N ) on the one hand and π ′ ∶ (1, 2, . . . , N ) ↦ ( , 1, . . . , N ) on the other, we get where, by virtue of integrability, we expressed the multi-particle scattering matrices from Eq. (72) as products of the two-particle ones -see Fig. 3 for a schematic guide. In expressing the difference p − p j of the momenta, our choice of permutations π and π ′ is arguably the simplest to factorise. For more complicated permutations, the factorisation of the scattering matrix can be determined, for example, by applying the Steinhaus-Johnson-Trotter algorithm. The two-particle scattering matrix is computed in Appendix C and reads By plugging it in Eq. (75) we finally obtain where M is defined in Eq. (51). The constraints on the momenta given by Eqs.  . The angular variable is 2πgI 0 N + π(g − 1) mod 2π; the radius is proportional to the maximal value of I >1 that can be reached for given I 0 (each dashed circle corresponds to an integer from L 2 to L). The arrow points at the circle representing the noninteracting sector. Symbols become larger and more transparent as the number of particles increases, whereas the colour represents the maximal I >1 associated with I 0 = 0, for each given configuration. The black dots identify points that are associated also with jammed configurations -the only ones in which, independently of I 0 , the maximal value of I >0 is L 2 + M − 1.
In summary, any Bethe state is specified by the configuration, b = (b 1 , . . . , b N ), of N particles (equivalently, the set B (79)

Quantum numbers
Remarkably, Bethe equations (74) and (77) can be solved explicitly, allowing us to overcome one of the complications of the Bethe Ansatz description. Specifically, the momenta p 1 , . . . , p N are given by where I ∈ Z and I 0 ∈ {0, 1, . . . , N − 1}. The total momentum P = ∑ N =1 p in turn satisfies Importantly, the phases e ip , = 1, . . . , N , do not change if one shifts I 0 → I 0 − M and one of the other quantum numbers, I n → I n + L 2 + M . Hence, all independent solutions can be found by imposing the folllowing constraints In practice, there are sets of quantum numbers corresponding to zero-norm states. This happens because there are other conditions, hidden in Eq. (78), that relate the coefficients of the state. To see that, let us introduce the cardinality B i.e., the product on its left-hand-side needs to be truncated to avoid repeated sequences b. For an arbitrary subset of N g momenta {p j }, j = 1, . . . , N g, we can now rewrite this condition using Eq. (78), obtaining a new Bethe equation for the total momentum Together with the constraint (77) imposed on the differences of momenta, this equation is finally solved by where 0 ≤ I 0 < N g and 0 ≤ I 1 < I 2 < . . .
Despite the fact that these conditions are sufficient to generate a complete basis of eigenstates with nonzero norm, the same rapidities (modulo 2π) can in fact be generated by different sets of quantum numbers. For example, for L = 6 there are 85 different configurations of quantum numbers satisfying Eq. (88) but only 64 independent states; for L = 4 there are 21 different configurations but only 16 independent states. We conjecture that the independent eigenstates can be obtained by exhausting all choices of the quantum numbers that satisfy the following constraints: where the functional F [I 0 ] and I max 0 are such that Here, d B where each binomial counts, in how many ways n j spins up can be distributed among the odd sites between 2 ′ j and 2 ′ j+1 . The first sum is over all the cyclic permutations that bring the sequence (n 1 , . . . , n N+ ) into an independent one, and thus exhausts all possible distributions of spins up on the odd sites. Our numerical investigations based on chains up to L = 20 sites suggest that, whenever possible, the bounds F [I 0 ] are such that I max The only circumstance where we found the theoretical maximal value N g − 1 unreached is when  Figure 4.
Finally, we note that the quantum number I 0 can be incorporated into the following rational quantum numbers which live in the affine lattice Z + gI 0 N + (η + g − 1) N . In terms of them we have whence the one-to-one correspondence between the N momenta and the N quantum numbers becomes apparent. Rational quantum numbers, however, depend on the state. In this regard, the correspondence between {p } and {J } somewhat differs from the standard scenario that arises in models solvable via standard Bethe Ansatz techniques.

Degeneracy of momenta
To compute the degeneracy, consider first all sequences b of ones and zeros, with M subse- 51)). The number of such sequences is where each term in the sum on the left-hand side represents the number of binary sequences at fixed M , composed of N 1 ones and N − N 1 zeros, and ending with a zero, b N = 0. The sum is over all N 1 compatible with the existence of M subsequences (1, 0), while the prefactor 2 comes from counting also the sequences that end with a one instead of a zero, i.e., b N = 1.
On the other hand, the number of N -digit binary sequences with a unit cell of size N g reads N g d( N g , M g ), and thus for any two positive integers n, m ∈ Z > . Here, µ(n) denotes the Möbius function, which yields the sum of all primitive n-th roots of unity. If n and m are common divisors of N and M , we can restrict the sums on the right-hand side of Eq. (99) to k N, M , whence we recognize the elements of the inverse matrix D −1 as (D −1 ) k,n = µ n k D k,n (all indices now run over the common divisors of N and M ). Finally, plugging this into ⃗ d = G⋅D −1 ⋅ ⃗ v, gives the degeneracies

Generic states
In the thermodynamic limit N, M, L → ∞ with the ratios ξ = L 2N and µ = M N fixed, the vast majority of states have g = 1, i.e., the configurations that characterise them have no periodic substructure. This can be inferred from Eq. (100) in the following way. First, we remind the reader that the number of distinct sets of momenta associated with a configuration depends only on N , M , and g. In addition, such a number is maximal when g = 1; indeed, for a larger g the norm of some states becomes zero. The total number of states with given g, N , and M is then proportional to d g,{p } , the proportionality factor being maximal for g = 1.
where we have isolated the first term and bounded the rest from below by replacing k with 2 in each term (in particular, we have replaced µ(k) → −1). Finally, we have used the fact that the number of common divisors is smaller than √ M . Analogously, we can bound d g,{p } , for g > 1, from above as which comes from replacing k with 1 (in particular µ(k) → 1) and g with 2, using again that the number of common divisors is smaller than for some γ > 0 depending on M and N . That is to say, the number of states with g > 1 is exponentially smaller than the number of states with g = 1. In view of this, we refer to the Bethe states with g = 1 as generic states.

Degeneracy at fixed staggered magnetisation
We consider here the effect of fixing the staggered magnetisation S z − , defined as the expectation value of the staggered spin along the z-axis -see Eq. (39). In a Bethe state at fixed configuration, S z − is the difference between the numbers N 0 and N 1 of particles ⊕ and ⊖, respectively: S z − = N 0 − N 1 . The number of binary sequences b at fixed N = N 0 + N 1 , M , and S z − reads where the first (second) term on the left-hand side counts the binary sequences of N 0 zeros and N 1 ones ending with one (zero). Note that, when the staggered magnetisation is not fixed, summing the left-hand side of Eq. (105) over the allowed values of N 1 yields exactly Eq. (96). For large N , the degeneracy at fixed staggered magnetisation can be approximated by the leading-order contribution to Eq. (105) divided by the number N of cyclic permutations of the binary sequence b:

Odd number of sites
By analogy with the even L case, the eigenstates in the sector characterized by the configuration B (o) can be represented by a set of N momenta (rapidities) {p 1 , . . . , p N }. We use the Bethe Ansatz where the ceiling function in the lower bounds of the sums accounts for the possibility of ⊖ and ⊕ sharing the same macrosite, i.e., ⊖⊕. The coefficients are again given in Eq. (69), but now satisfy boundary conditions that follow from σ x L+n = (−1) η σ z n and account for the transmutation of the particle crossing the boundary. Applying the Ansatz to the boundary conditions and removing the dependence on ′ j we then find As a result of the factorisation of the scattering matrices, shown in Fig. 3, and of the explicit form of the two-particle one, (76), this is equivalent to Here M − , defined as is invariant under "anti-cyclic" permutations of the configuration: b N +1 ≡ 1 − b 1 . Since the ratio (110) holds for any momentum p j , we obtain On the other hand, the analogue of Eq. (73) with 2N cyclic permutations instead of N yields the conservation of momentum The energy is again given by Eq.

Quantum numbers
Like in the even case, the Bethe equations can be solved explicitly: where I ∈ Z and I 0 ∈ {0, 1, . . . , 2N − 1}. The phases e ip , for = 1, . . . , N , do not change under the transformation I 0 → I 0 − 1 − 2M − and I n → I n + L+1 2 + M − , for a given n > 0. Therefore, all the independent solutions can be found by imposing To avoid the zero-norm Bethe states, Eq. (113) has to be amended similarly as in the case of even L. Specifically, if the configuration (b; 1 − b) has a periodic substructure with a unit cell of size B (o) N = 2N g (copying the unit cell g-times yields the configuration), the product of ratios in Eq. (73) has to be truncated accordingly: Using the explicit form of the ratio, given in Eq. (110), this becomes for any set of 2N g momenta {p j }, j = 1, . . . , 2N g. Together with Eq. (112), it is solved by Remarkably, calculation of the momenta shows that they do not depend on η. In fact, it seems that the last term in Eq. (118) can safely be ignored.
We warn the reader that, like in the even case, a set of rapidities could be generated by different sets of quantum numbers and, in order to have a one-to-one correspondence between quantum numbers and momenta, the allowed values for the quantum numbers should be appropriately reduced.
We conclude here the discussion about Bethe states; for the sake of simplicity we assume, in the rest of the paper, that the number L of spins is even.

A family of local conservation laws
The local conservation laws of integrable systems are usually constructed within the framework of algebraic Bethe Ansatz. Indeed, the range of a charge cannot be easily inferred from coordinate Bethe Ansatz. Nevertheless, we conjecture that the local conservation laws introduced in Section 3.1 have single-particle eigenvalues equal to They notably form a complete basis of square-integrable functions on a circle. We have numerically checked this conjecture up to n = 4 by comparing the spectrum of the local conservation laws obtained by imposing the vanishing of their commutator with H F against the Bethe Ansatz predictions. 10 Let us compute how many charges of this form are independent in a system of finite size L. To that aim, we introduce non-Hermitian charges Z n = Q + n + iQ − n , with eigenvalues The knowledge of the expectation value of a generic conservation law with eigenvalue ∑ N =1 f (p ) is equivalent to the knowledge of the set of momenta {p }. Since the momenta are defined modulo 2π, we can assume f (p) to be a 2π-periodic function and write it in the Fourier space This equation implies that the knowledge of I 0 , P , and ⟨Z n ⟩, with n = 1, . . . , L 2 + M , unambiguously fixes the momenta through a functional that does not depend on L explicitly (the latter appears only in the range of the index of charges). In fact, also the dependence on I 0 and P can be removed by specialising the identity to f (p) = 1: Using this, we can express Eq. (121) in terms of ⟨Z n ⟩ andf n only. Assuming, without loss of generality, that the single-particle eigenvalue f (p) is real, we thus obtain This implies that, for any given M , only O L 2 + M Hermitian charges are actually independent. Note also that this relation involves the eigenstates of the operator M, which we have shown in section 5.1.2 to be completely determined by this family of charges (at finite L).
Since the number of independent conservation laws scales linearly with the system size, the Fourier decomposition in Eq. (121) suggests that, in the thermodynamic limit N, M, L → ∞, at fixed ξ = L 2N and µ = M N , the eigenvalue ∑ N =1 f (p ) is linear in the expectation values ⟨Z n ⟩, for n ∈ Z > ≡ N. Indeed the last term of the expansion vanishes in the thermodynamic limit, provided that the Fourier series of f (p) converges. Additional insight is gained by using Eq. (123) to express ⟨Z L 2 +M ⟩ in terms of ⟨Z L 2 ⟩ and ⟨Z M ⟩, and similarly for charges with index n > L 2 + M . In this way Eq. (121) can be shown to be equivalent to for finite L, N and M . Since M ≤ L 2 and ⟨Z 0 ⟩ = N , this equation reveals that the momenta are completely determined by the expectation values of L 2 non-Hermitian operators or, equivalently, L Hermitian ones. These operators extend the family of local charges to all possible ranges.

Low-energy limits Ground state
While the XXZ model with anisotropy ∆ > 1 is gapped antiferromagnetic, due to the absence of the rapidly oscillating part κ −1 H I , the folded Hamiltonian has a critical ground state. The latter is a Fermi sea with two chiral modes, i.e., it is conformal critical, the central charge being c = 1. In view of this, we identify it with a stationary state in which the neighbouring momenta 11 Note that this relation between expectation values can be lifted into a relation between operators (it holds for all their eigenstates). In particular, we have  are occupied, i.e., quantum numbers I j , for j = 1, . . . , N , are consecutive: I j = I 1 + j − 1 12 . Using this Ansatz, the ground state energy reads Since g, I 0 , and I 1 appear only in the argument of the cosine, we can minimize the energy by maximizing the ratio of sines and minimizing the cosine. At fixed average distance ξ = L 2N between the particles, the maximal value of the ratio of sines is reached at µ = M N = 1 2. At finite L this implies that N is even and g = N 2: the configuration of the ground state is B N = {(0, 1, 0, 1, . . . , 0, 1)} c . The cosine, on the other hand, is minimal when its argument is an odd multiple of π. Note that there are so many integer degrees of freedom that the minimum −1 is always reachable. Defining the variable 13 which takes the value ϕ = 2πL N +L in the ground state, the energy Ansatz reduces to As a function of the parameter ϕ, the energy is minimal at the solutionφ L to the transcendental equation tanφ L = L tan(L −1φ L ). The dependence of the parameterφ L on the system size L can be taken into account perturbatively as whereφ is the solution to tanφ =φ in the interval [π, 2π), which is independent of the system size L. Expansion of the energy in ∆ϕ = ϕ −φ L now yields The final step is to minimise the lattice displacement ∆ϕ, which is nonzero because L+N 2 = πL ϕ is an integer. We find which, plugged into Eq. (130), yields where o(L −1 ) means that the correction approaches zero faster than L −1 -see Fig. 5. The term originating in the displacement ∆ϕ is a clear lattice effect: it arises because L and N are (even) integers. On the other hand, the O(L −1 ) correction that remains when the integer condition is relaxed (i.e., keeping only the first term in the rounded brackets in Eq. (132)) is expected to be universal and describable by the underlying conformal field theory. Specifically, the prediction for periodic boundary conditions is [36,37] where c is the central charge, v F is the Fermi velocity, and ε 0 is the ground state energy density in the infinite volume, which is not universal (note that the volume is L 2 because the macrosites consist of two adjecent spins). Since the theory under investigation has two chiral modes, the central change c is equal to 1. By equating the volume correction that we computed with the conformal field theory prediction we can then infer the Fermi velocity This result is in agreement with the thermodynamic Bethe Ansatz calculation worked out in the second part of this work [35], as well as with the calculation in the constrained model, described by the Hamiltonian (32) -see Ref. [20]. This is another indication that constrained Hamiltonians, in which particles have a finite radius, correspond to low-energy limits of the asymptotically folded Hamiltonians. Figure 5 compares the CFT correction against the full leading one (132) and the exact ground-state energy for small chains up to L = 20.
Finally, we note that the ground state has a finite negative magnetisation 1  Figure 6: Energy per unit length ε = E L and µ = M N of the eigenstates (black symbols) in the low-energy part of the spectrum ofH 0 F , for L = 18. The orange curve is prediction (145), which is valid in the thermodynamic limit L → ∞.

Average size of domains in minimal-energy states
In the previous subsection we have shown that, in the ground state, M = N 2 and hence g = N 2. The Ansatz (126) for the energy is, however, still valid if we impose a different value for µ and minimise the energy. For µ ≠ 0, the result will then be the minimal energy at which the groups of consecutive particles of the same species, i.e., the domains, have an average size of (2µ) −1 . As we will see, the minimal energy increases with the average number of consecutive identical particles. The exact calculation is very similar to the one carried out in the previous subsection. Again, the cosine in the Ansatz (126) can be replaced by −1, so that the energy reads where ϕ is again given in Eq. (127) and depends on the parameter µ. The energy is minimal at ϕ =φ L , satisfying the transcendental equation tanφ L = L 2µ tan L −1 [2µφ L + (1 − 2µ)2π] . As before, we take into account the dependence of the parameterφ L on L perturbatively, whereφ is the solution to 2µ tanφ = 2µφ + (1 − 2µ)2π. It lies in the interval [π, 3π 2) and is independent of the system size L. We can now expand in ∆ϕ = ϕ −φ L : The non-zero lattice displacement ∆ϕ is now obtained by imposing πL 2µϕ+(1−2µ)2π = L 2 + M ∈ N, which results in Finally, the energy is expressed as and using the conformal field theory prediction (133) we now find This again agrees with the thermodynamic Bethe Ansatz calculation worked out in the second part of this work [35]. Differently from the genuine ground state with energy (132), the minimal energy states with energy E min have exponentially large degeneracy; indeed only µ = 1 2 (in the genuine ground state) fixes the configuration unambiguously. A question that now arises is, how the ratio µ min = M N changes as a function of the energy density in the minimal energy state. Here,φ depends on µ min , and the dependence of the latter on ε is what we wish to obtain. Eq. (141), on the one hand, implies where the positive branch of the square root has been chosen because the relevant solution ϕ to the transcendental equation 2µ tanφ = 2µφ + (1 − 2µ)2π falls into the interval [π, 3π 2). On the other hand, the transcendental equation, along with Eq. (141), also yields where, again, the branch of arccos is chosen in such a way thatφ ∈ [π, 3π 2). Noting that ε < 0, we now have We now recognise the left-hand side to be equal to ∫ εµ min 2J 1 dy arccos y; the solution can then be expressed implicitly as where f (x) is a solution to ∫ f (x) 1 dy arccos y = x -see Fig. 6. The Heaviside step function θ comes from the fact that εµ min 2J and πε 2J should both have the sign of ε, and Eq. (144) is consistent with such a constraint only if ε ≤ −2J π.
For the sake of completeness we also report the value of the minimal number M of oriented domain walls in the eigenstates with a given energy density ε:

Conclusion
We have proposed a framework for studying time evolution in quantum many-body systems on intermediate time scales at large coupling constant. As an example, we have investigated the effective Hamiltonian that arises in the limit of large anisotropy ∆ in the Heisenberg spin-1 2 XXZ model. We have developed a coordinate Bethe Ansatz that manifests the symmetries emerging in the strong coupling limit. Differently from the Bethe Ansatz solution of the XXZ model [38], the new Ansatz is based on two species of particles that live on macrosites, comprising two neighbouring spins. As a consequence, the Bethe states explicitly break the translational symmetry of the effective Hamiltonian. The Bethe states with both species of particles form an interacting sector, while those with a single species of particles form the noninteracting one. The latter sector can be mapped into an effective XX model, which is well known for its non-abelian set of conservation laws [15]. How they fit into the model's overall integrable structure is an interesting problem, left to future investigations. We have also identified excited states in which particles populate the lattice in a chaotic way and are packed so closely that they keep fixed positions. These jammed states form an exponentially large subspace of the Hilbert space, which we showed to span the eigenspace of the ground state of a quasilocal charge M − S z . We also found that a subset of such states is stable under the O(∆ −1 ) correction to the leading-order effective Hamiltonian. Whether or not the leading correction preserves integrability is still an open question. All the above properties make the proposed effective model a particularly convenient playground for studying pre-relaxation phenomena.
We have found that the ground-state energy and Fermi velocity agree with those calculated by Alcaraz and Bariev in the constrained XX Hamiltonian, describing finite-radius particles [20]. Their projected Hamiltonian corresponds to the four-vertex model and we speculate that it describes the low-energy sector of our effective Hamiltonian. In fact, Alcaraz and Bariev provided also a Bethe Ansatz for a range of larger particle radii; it would be interesting to see how the other constrained Hamiltonians fit in the Bethe Ansatz that we considered.
Finally, we remark that the simplicity of the Bethe Ansatz equations makes the effective Hamiltonian originating in the large-anisotropy limit of the XXZ model an excellent candidate for advancing the theory of generalised hydrodynamics in interacting systems. To set the stage for it, the second part of this work [35] will develop a thermodynamic Bethe Ansatz and a generalised hydrodynamic theory on the Euler scale.

A Asymptotic expansion
For the sake of completeness, we recover the result of Ref. [16] and its generalisation by showing the asymptotic validity of decomposition (5), which is equivalent to where we have used Eqs. (2) and the fact that B z (κ) is a functional of the form To show the asymptotic validity of Eq. (147), we now define the operator satisfying the differential equation where z t = e iJt κ . Assuming decomposition (147), on the other hand, we have and hence i∂ t V(t) = e iH F t −H F + −κ −1 Jz t ∂ zt e −iκBz t (κ) + e −iκBz t (κ) H F (κ) e iκBz t (κ) e −iH F t V(t).
(152) If decomposition (147) is possible, this expression has to match Eq. (150), whence a gauge transformation follows for any z ∈ C, with z = 1. Here, ∂ z denotes the derivation operator, i.e., both sides of the equation should be read as if applied to some differentiable function of z.
The existence of a solution to this equation is the necessary and sufficient condition for the validity of Eq. (147). Using Eq. (4) and notation A † n (κ) ≡ A −n (κ) ≡ A − n (κ) ≡ A + −n (κ) we can recast Eq. (153) as an infinite system of equations B Folded picture in noninteracting spin-1 2 chains In this appending we consider an example -a noninteracting spin-1 2 chain with H I ∝ ∑ σ z -in which the folded picture can be worked out in a non-perturbative way. Other cases could be addressed by applying first a noninteracting mapping that sends H I into an operator proportional to ∑ σ z . Specifically, we consider Hamiltonians that can be written as H = 1 4 ,n a H n a n , where a are the Majorana fermions defined in Eq. (21). The matrix H is conveniently written as a Fourier transform, which is used to define the so-called "symbol" of the noninteracting operator. Let us show the procedure by considering the Hamiltonian of the Ising model (18), with the matrix elements Its symbol is a 2 × 2 matrix function appearing in the Fourier transform where e iθ(k) = Analogously, the symbol of H I is − J 2 σ y . A very important property is that the symbol of the commutator of two noninteracting operators is the commutator of the corresponding symbols. This can be readily exploited to exhibit a unitary transformation that maps H into an operator commuting with H I .
To that aim, let us consider the noninteracting operator W(h −1 ) represented by From Eq. (159) and the Baker-Campbell-Hausdorf formula it follows that the noninteracting operator e iW(h −1 ) He −iW(h −1 ) has the symbol −ε(k)σ y , which commutes with that of H I . This, in turn, implies We can then identify the folded Hamiltonian with that is to say, its symbol readŝ The unitary transformation is instead given by and, from B zt (h −1 ) = e ihH I t B 1 (h −1 )e −ihH I t , it follows that only A 1 (κ) is different from zero, its symbol readingâ

C Energy and scattering matrix
The bare energy of quasiparticle excitations can be calculated by considering the one-particle sector, N = 1. In addition, if a coordinate Bethe Ansatz can be built from such quasiparticles, their scattering properties are completely determined by the solution in the two-particle sector, N = 2, which is supposed to completely characterise the scattering properties. The other sectors are then solved by a Bethe Ansatz, which is presented in full generality in Section 5.

One-particle sector
Even number of sites. In the one-particle sector there are only two possible configurations, B This is solved by plane waves with energy E {(b)}c (p) = 4J cos p.
As expected by one-site shift invariance, the excitation energy does not depend on the particular configuration. Observe that no eigenstate in this sector is mapped into an eigenstate ofH F : since L is even, the vacuum belongs to the sector Π z = 1, hence the single-particle excitations ofH η F are in the sector Π z = −1, which is not associated with eigenstates ofH F -cf. Eq. (41).
Odd number of sites. Normally, one would look for a solution of the form which has energy E(k) = 4J cos(2k), as inferred from the action ofH η F . It is however more convenient to enforce the description in terms of particles lying on macrosites. From that perspective, we must distinguish between even and odd sites Concerning the particle configuration, we note instead that, up to cyclic permutations, there is a single possibility: {(0; 1)} c ≡ {(0; 1), (1; 0)}. It is now natural to identify the momentum of the particles with p = 2k, so that also the energy keeps the same form as in the even-size case, i.e., E {(0;1)}c (p) = 4J cos p.
Quantisation condition for p is obtained by squaring the one for k, that is, This results in some ambiguity in the relation between k and p, which can be lifted by imposing the boundary condition σ x +L = (−1) η σ x . Indeed, using it in the second sum of Eq. (172) and demanding the uniqueness of the state, we obtain where the quantisation condition (174) was imposed in the second equality. Finally, we have Remarkably, only the eigenstates depend on η (the momenta do not). Since, in this sector, one has Π z = 1, each eigenstate is mapped into one ofH F through Eq. (41). Their eigenstates are simply given by

Two-particle sector
where the relative minus sign is due to the fermionic nature of the scattering process. The excitation energy is the sum of the single-particle excitation energies The case of configuration B Note that the sum in the second term allows for the particles ⊖ on the left-hand side and ⊕ on the right-hand side, to share the same macrosite (we denote it by ⊖⊕). The Ansatz that we impose upon the coefficients is now c b 1 ,b 2 ′ ,n ′ (p 1 , p 2 ) = Z b 1 ,b 2 p 1 ,p 2 e i ′ p 1 +in ′ p 2 + S b 1 ,b 2 p 1 , p 2 p 2 , p 1 e i ′ p 2 +in ′ p 1 , for generic ′ and n ′ . From them we obtain 191) and the Bethe equations e iL(p 1 +p 2 ) = 1, e i L+1 2 (p 1 −p 2 ) = 1.