The Folded Spin-1/2 XXZ Model: II. Thermodynamics and Hydrodynamics with a Minimal Set of Charges

We study the (dual) folded spin-1/2 XXZ model in the thermodynamic limit. We focus, in particular, on a class of local macrostates that includes Gibbs ensembles. We develop a thermodynamic Bethe Ansatz description and work out generalised hydrodynamics at the leading order. Remarkably, in the ballistic scaling limit the junction of two local macrostates results in a discontinuity in the profile of essentially any local observable.


Introduction
This is the second part of our investigation into the large-anisotropy limit of the Heisenberg spin-1 2 XXZ model. In the first part [1] we defined the "folded picture" as an asymptotic formulation of quantum mechanics of many-body systems that are described by Hamiltonians with a large coupling constant. In this picture the fast oscillatory part of the dynamics is attached to the operators, whereas the state evolves slowly in time through an effective "folded Hamiltonian". As an example we considered the folded XXZ model, which reveals a structure otherwise hidden in the standard Bethe Ansatz solution [2] of the anisotropic Heisenberg magnet. Our solution of the folded XXZ model is peculiar as it diagonalises the folded Hamiltonian, but not the total spin S z in the direction of the anisotropy. This is unveiled by a duality transformation that maps the folded Hamiltonian into a block-diagonal local operator and S z into a block off-diagonal pseudolocal one.
The present manuscript is a detailed account of the thermodynamic properties of the dual folded XXZ model, both in and out of equilibrium. Up to boundary terms, the Hamiltonian that we consider reads as J ∑ (σ x −1 σ x +1 + σ y −1 σ y +1 ) 1−σ z 2 , where σ α are the Pauli matrices. It describes a special point of the two-component Bariev model [3], which is solvable with a nested Bethe Ansatz technique. In order to exploit the symmetries of that special point, we will however use the results of the non-nested Bethe Ansatz method of Ref. [1]. We first develop a thermodynamic Bethe Ansatz that describes the infinite chain in thermal states as well as in generalised Gibbs ensembles [4] constructed with the local conservation laws and a special pseudolocal charge. The predictions of the thermodynamic Bethe Ansatz are checked against numerical data from DMRG algorithms.
The macrostates represented by generalised Gibbs ensembles are the key ingredients of the generalised hydrodynamic theory. Introduced in Refs [5,6], and recently experimentally corroborated in cold-atom setups [7], this mesoscopic description of the dynamics in integrable systems allows one to treat the evolution of inhomogeneous states on large space-time scales, and also in presence of inhomogeneous space-and time-dependent interactions [8,9]. Generalised hydrodynamics is based on the assumption that the so-called root densities, which characterise the local properties of stationary states, can be promoted into functions of space and time. The validity of such an assumption is hard to establish rigorously [10,11], and it is typically justified by comparing the predictions with state-of-the-art numerical simulations. Provided that the assumption holds, the state can be coarse-grained into fluid cells locally equivalent to macrostates, which can be accessed through the thermodynamic Bethe Ansatz.
Within this framework, we investigate the transport phenomena that emerge after two grand canonical ensembles are joined together and left to evolve in time under the folded Hamiltonian. Remarkably, due to the rich particle content of the folded XXZ model, a local macrostate is not completely determined by a root density, and an additional independent (pseudolocal) thermodynamic variable is necessary for describing the junction of two local macrostates. We show that this special feature results in discontinuous ballistic-scale profiles of all local charges, independently of how the initial grand canonical ensembles are prepared.
In the rest of the introduction we recapitulate the results of Ref. [1] necessary to understand the present work.

The asymptotic folded picture
When a Hamiltonian has a large coupling constant, it is possible to derive asymptotic expansions of the time evolution operator in the strong coupling limit and at fixed time. Starting from Ref. [12], explicit results have been obtained for Hamiltonians H(κ) of the form where q is a finite integer and In Ref. [1] we have recast the strong coupling expansion κ → 0 of such models into a formulation in which operators and state time evolve according to where z t = e iκ −1 Jt . An asymptotic expansion of H F (κ) and B z (κ) in the limit of small κ is reported in Ref. [1]. As long as Jt ≪ κ −1 , the folded picture can be used at the leading order, where the state time evolves with H F (0) ≡ H F . This was used for example in Ref. [13] to explain the very slow restoration of one-site shift invariance in XY and XXZ models. Following this perspective, the present paper investigates time evolution under H(κ) in the asymptotic limit 1 ≪ Jt ≪ κ −1 ≪ L → ∞. In fact, we will assume the relation between H(κ) and H F as understood, and interpret H F as the Hamiltonian of a new model, which we refer to as "folded" model. We are using this nomenclature because, in the limit κ → 0, the spectrum of H(κ) turns out to be equal to the spectrum of H F , modulo a typical energy, inversely proportional to κ.

The folded XXZ model and its dual
The XXZ spin-1 2 chain is described by the Hamiltonian We are interested in the limit of large anisotropy ∆. In the folded picture we identify κ with (4∆) −1 , and the operator H F is given by In the following we provide a brief overview of the diagonalisation of H F that we worked out in the first part of our work [1]. The preliminary step is a duality transformation mapping H F into a local operator with a density that has a range of three sites. The transformation reads where L is the chain's length. Under the transformation (6), the folded Hamiltonian with periodic boundary conditions (σ α L+n = σ α n , for α ∈ {x, y, z}) is mapped into the following operatorH where Π z = ∏ L =1 σ z , and The eigenstates ofH η F with Π z = 1 are mapped into eigenstates ofH F either trivially, or by the unitary transformation σ x L , depending on whether η = 0 or η = 1. It is therefore convenient to focus onH η F . As already mentioned, the latter Hamiltonian describes a strong repulsion limit of the two-component Bariev model [3]. In our previous work [1] we circumvented the more general solution (based on nested Bethe Ansatz) by exploiting the additional symmetries emerging in such a strong repulsive regime. This allowed us to solve the Bethe equations in terms of closed-form relations between momenta and quantum numbers, which will be summarised in the next section.

Coordinate Bethe Ansatz
A basis of eigenstates ofH η F can be constructed within a coordinate Bethe Ansatz starting from the reference state vac⟩ = ↓↓ . . . ↓⟩ .
For the sake of simplicity we assume that L is even; the reader can find some information about the odd case in Ref. [1]. If we indicate by 2 ′ j − b j the positions of the spins up in the state, the set B N = {(b 1 , . . . , b N )} c of cyclic permutations of the sequence (b 1 , . . . , b N ) is preserved by the Hamiltonian. Here, N is the number of spins up. In the following we will refer to the information encoded in B N as that pertaining to the "configuration space". The eigenstates are characterised by N momenta (rapidities) {p}, each one associated with a spin up in an even (b = 0) or odd (b = 1) position. Interactions are characterised by the two-body scattering matrix The energy of the state p 1 , . . . , p N ⟩ B N , specified by the momenta and the configuration, is given by The momenta p 1 , . . . , p N solve Bethe equations that can be integrated explicitly, as shown in Ref. [1]. The solution reads where the total momentum and the shift of quantum numbers are defined as respectively. The integer quantum numbers I satisfy N g = min{m b n+m = b n , ∀n} denoting the size of the unit cell in the sequence (b 1 , . . . , b N ), and M = ∑ N j=1 b j (1 − b j+1 ), with b n+1 = b 1 the number of subsequences (1, 0). We point out that some configurations of integers within the domain specified by Eq. (14) give rise to the same set of momenta. This subtlety, however, does not affect the thermodynamic limit, so we will not be more specific about it; the interested reader can find more details in Ref. [1]. Since the momenta are invariant under the transformation I → I − n, ϕ → ϕ + 2π n, with n ∈ Z, the parameter ϕ can be defined modulo 2π.
An alternative parametrisation of the solution to the Bethe equations is through the rational quantum numbers J = I + ϕ 2π , in terms of which we have We emphasise that the rational quantum numbers lie in the affine lattice Z + ϕ 2π , whose shift, with respect to integers, depends on the state itself. Importantly, for large N almost all the states have g = 1 [1]. In light of this, we will refer to them as generic states.

Conservation laws
In the first part of our work [1] we have shown that, in finite chains, the momenta of the Bethe Ansatz characterise the eigenvalues of an "extended" 1 family of local conservation laws. Among them, we have identified two charges that are diagonal in the standard σ z basis: the latter operator's eigenvalue being equal to the parameter M in the Bethe equations (12). The remaining charges of the family are not diagonal and can be organized into two sequences: The local density q ± n, can be defined so as to be supported on 2n + 1 neighbouring sites and, by convention, Q + 1 =H η F . The expectation value of any charge Q ± n in the Bethe state can be written as where the sum runs over the momenta. Here, q(p) denotes the one-particle eigenvalue of the charge, for example, q + 1 (p) = 4J cos p = E(p) is the one-particle energy. In Eq. (17) we subtracted the vacuum expectation value of the charge on the right-hand side, since the reference state vac⟩ has no momenta, i.e., the left-hand side of the equation is zero in that case.
In general, the single-particle eigenvalues of the charges in the integrable hierarchy read and notably span the Fourier basis in the space of functions of the momentum. As an example, we state the local densities of the first few conservation laws, Q ± 1 and Q + 2 : where K n,m = σ x n σ x m + σ y n σ y m , D n,m = σ x n σ y m − σ y n σ x m . For any fixed n, the conserved charges Q ± n are connected by the adjoint action of In integrable spin chains one usually refers to such an operator as a "boost" or "ladder" operator if it allows the reconstruction of the entire hierarchy of local conserved quantities from a single charge. In our case it is evidently not so. Moreover, there seems to be no local boost operator of this form, by means of which one could translate between conserved charges with different n. This is consistent with the fact that the energy current (its explicit form is provided in Section 4.2) is not conserved [14], as well as with the fact that the R-matrix that constitutes the Algebraic Bethe Ansatz for the two-component Bariev model is not of the difference form [15].
We point out that the family of local charges with local densities (19) is not complete, and we can even exhibit a local charge that does not belong to the family: the staggered spin along the z-axis We will denote its expectation value by S z st = L 2 m z st and refer to it as staggered magnetisation. The staggered spin along the z-axis is one of the diagonal operators that commute with the Hamiltonian and span the configuration space. The common property of such diagonal charges is that their eigenvalues are functionals of the configuration B N . Summary Section 2 considers the thermodynamic limit of the Bethe equations and develops a (generalised) thermodynamic Bethe Ansatz description. Thermal states are analysed in detail, including low-temperature and high-temperature asymptotic expansions.
Section 3 proposes a definition of elementary quasiparticles (excitations) and works out their dressed charges.
Section 4 develops generalised hydrodynamics (GHD) in the dual folded Hamiltonian for a class of states that are characterised by a minimal set of charges.
Section 5 uses GHD to investigate time evolution after two (local) macrostates have been joined together. It is shown that the profiles of essentially all local observables remain discontinuous in the ballistic scaling limit.

Thermodynamic Bethe Ansatz
The thermodynamic Bethe Ansatz (TBA) was originally developed by Yang and Yang for the one-dimensional Bose gas with Dirac-delta repulsive interactions [16]. In this section we will follow their approach with few minimal changes necessary for accommodating the TBA in the richer structure of the folded XXZ Hamiltonian. We note that an alternative, albeit more complicated, approach towards TBA exists for our model. It can be obtained as a particular limit of the nested Bethe Ansatz solution to the thermodynamics of the more general multi-component Bariev model [17].

Root densities and expectation values of charges
We start by noting that the solution (12) to the Bethe equations can be recast into the standard form Here we have introduced the monotonous counting function h(p): ∂ p h(p) = 1 + µ ξ > 0, where we denoted µ = M N and ξ = L 2N . Whenever evaluated in the momentum that solves the Bethe equations, its value falls into the affine lattice Z + g N I 0 + η+g−1 2 populated by quantum numbers J . Bethe equations (21), in particular the prefactor L 2 in front of p , manifest that our momentum generates translations for two sites on the spin chain.
The counting function associates certain elements of the affine lattice Z + g N I 0 + η+g−1 2 to momenta that form a particular solution of Bethe equations. Specifically, a Bethe state with quantum numbers {J 1 , J 2 , . . . , J N } contains momenta {p 1 , p 2 , . . . , p N } that solve L 4π h(p ) = J .
If, for instance, J 2 + 1 is not among the quantum numbers {J 1 , J 2 , . . . , J N }, the momentum p, for which L 4π h(p) = J 2 + 1, does not belong to the set {p 1 , p 2 , . . . , p N } of momenta in the Bethe state. We refer to the vacant quantum number J 2 + 1 as a hole, in the sense that adding that quantum number corresponds to an elementary excitation -see Section 3. In the thermodynamic limit (TD) N, M, L → ∞, with fixed ratios µ = M N and ξ = L 2N , the momenta characterising the excited states become densely distributed. Then, L 4π dh(p) yields the number of vacancies (both, particles and holes) in the infinitesimal interval [p, p + dp) ⊂ [−π, π). Defining the root density ρ and the density of holes ρ h as L 2 ρ(p)dp = number of particles with momentum in [p, p + dp) , L 2 ρ h (p)dp = number of holes with momentum in [p, p + dp) , we obtain dh(p) = 2π[ρ(p) + ρ h (p)]dp. The derivative of the counting function is thus proportional to the total density of vacancies which is notably independent of the momentum. From Eqs (17) and (22) it then follows which is the standard way to express the charge densities as functionals of the root density in the thermodynamic limit.
We are now in a position to recast Eq. (23) in the form of a standard TBA equation. Indeed, the macrosite density of particles is given by whence we have -cf. Eq. (23): It is also customary to define the filling function as the density of occupied vacancies, namely, n(p) = ρ(p) ρ t , which clearly satisfies 0 ≤ n(p) ≤ 1. We then have the analogous equation

Macrostates and Yang-Yang entropy
In the thermodynamic limit different states can share the same local properties and are usually said to be "locally equivalent". The concept of state is then replaced by the concept of macrostate, which represents the set of all locally indistinguishable states. It is usually understood that macrostates are stationary, and this will be assumed in this section.
We consider first the macrostates characterised by the strictly local integrals of motion; in order to distinguish them from more general macrostates, we will refer to them as "local macrostates". A local macrostate corresponds to a set of microstates (representative states) that, by construction, are simultaneous eigenstates of all quasilocal conservation laws [18]. The expectation value of a quasilocal charge takes, however, the most probable value at fixed local integrals of motion. In particular, as long as the Hamiltonian is local, the equilibrium properties at finite temperature are described by local macrostates. The situation is more complicated when the system is out of equilibrium. For example, local macrostates in the XXZ model were originally conjectured to describe the late-time properties after quantum quenches [19,20]. The failure of such an assumption was reported in Refs [21,22], and Ref. [23] later pointed out that the model possesses additional quasilocal charges (i.e., conserved operators with exponentially localised densities) that constrain the dynamics to a greater extent. It soon became evident that the state at late times becomes equivalent to a "quasilocal macrostate", characterised by all quasilocal integrals of motion [24]. In our specific case the stationary properties in the particle space are completely determined by the local charges and a single quasilocal one, M. The rest of the charges, diagonal and generally quasilocal, characterise solely the configuration space, i.e., their eigenvalues are functionals of the configuration B N of a Bethe state.

Local macrostates
Except for the staggered magnetisation ⟨S z st ⟩ ≡ S z st = L 2 m z st , all local integrals of motion that we identified are completely determined by the momenta (rapidities), which depend on the configuration through M . An important representation of a local macrostate is the (generalised) canonical state, which is described by a density matrix of the form Such a state minimises the "generalised free energy", namely, the functional 2 under a generic trace-preserving variation of ρ. A local macrostate of the folded XXZ model is characterised by three thermodynamic quantities: µ, m z st , and the root density ρ(p). In turn, the variation of ρ is realised by an arbitrary variation of µ, m z st , and ρ(p) in the functional f loc [µ, m z st , ρ], which represents the thermodynamic limit of f loc [ρ]. The thermodynamic limit of the second term on the right hand side of Eq. (29) has already been worked out -it is an expectation value of the form (24). The first term, on the other hand, is the entropy of the state per macrosite and is proportional to the logarithm of the number of microstates associated with the particular macrostate [16]. In our specific case, it consists of two terms: the entropy in the configuration space and the entropy in the particle space; both are computed below.
Entropy in the configuration space. macrosite is therefore where H(p) = −p log p − (1 − p) log(1 − p) is the binary entropy function. An analogous calculation gives the entropy at fixed staggered magnetisation S z st = L 2 m z st [1], which reduces to the previous expression for m z st = 0.
Entropy in the particle space. In the thermodynamic limit the distribution of momenta is encoded in the root density ρ(p). Following Ref. [16], the entropy associated with the number of microstates with the same root density is where index p in the product runs over the momenta representing the centres of N cells of width dp; they are coarse-grained momenta in the interval [−π, π).
Yang-Yang entropy. We call Yang-Yang entropy the thermodynamic limit of 2 L log Ω, where Ω is the size of the space associated with the macrostate. In our specific case, it is given by the sum of Eqs (31) and (32), namely, Here, ξ and ρ t must be interpreted as the functionals of ρ(p) and µ shown in Eqs (25) and (26), respectively.
We can now work out the variation of f loc [µ, m z st , ρ]. To that aim, we represent the local macrostate as where Q can be any linear combination of the local charges forming the integrable hierarchy described in Section 1.4, namely, S z or Q ± n . According to Eq. (24), the expectation value of the charge Q can be written as where q(p) is its single-particle eigenvalue. By definition, the staggered magnetisation per unit macrosite is equal to m z st , the generalised free energy thus reads We can readily find the minimum of the functional f loc [µ, m z st , ρ]. It is reached in the state with the staggered magnetisation and the filling function where µ is implicitly defined as the solution to the equation Alternatively, we can interpret these equations as the statement that a local macrostate with given local integrals of motion satisfies π −π dp 2π The latter equation can be used to express µ as a functional of ρ(p) and m z st ; this is the point of view used when describing the late-time dynamics after quantum quenches. There, the initial state fixes the integrals of motion, i.e., in our specific case, ρ(p) and m z st . In Section 5.2 we will use Eq. (39) to show that the locally quasistationary state emerging at late times after joining two thermal states is not locally equivalent to a local macrostate. This eventually necessitates the generalisation of local macrostates to quasilocal ones, which we discuss in Section 2.2.2.
Range of µ in local macrostates. The height µ[m z st , ρ] of the surface characterised by Eq. (40) extends over a restricted interval of µ. In particular, we have where we used the Jensen's inequality ∫ . We note that the inequality in Eq. (41) is saturated for ρ(p) = 1 2πξ , which is a physical value for the root density. Thus, for generic ρ(p), the bound on µ that originates in Eq. (41) can not be improved.
The parameter ξ is also constrained: we observe m z Regions of allowed µ and ξ, for different values of m z st , are plotted in Fig 1. For less generic states it is instead useful to start from Eq. (39) and use the Jensen's inequality for the convex function f (x) = log(1 + e w e −x ). In that case one has Figure 1: Bounds on µ and ξ, imposed by inequalities (41) and (42), for different values of staggered magnetisation m z st .

A class of quasilocal macrostates
Equation (40) where Q is, again, any combination of local charges Q ± n and S z . The generalised free energy is now given by The minimum of the functional f q-loc [µ, m z st , ρ] is the state in which staggered magnetisation and filling function are again given by Eqs (37) and (38), respectively, while µ is now implicitly defined as the solution to the equation Contrary to Eq. (39), the additional parameter χ allows us to vary µ (almost) independently of the other integrals of motion. We point out that Eq. (45) does not describe the most general quasilocal macrostate, as there are still infinitely many quasilocal charges in the configuration space [1]. It will nevertheless suffice to describe the late-time behaviour after the junction of two local macrostates (e.g., thermal states). More generally, indicating with C a generic quasilocal charge in the configuration space, the most general quasilocal macrostate can be represented as follows: Like the staggered spin along the z-axis, the charge C affects the entropy in the configuration space, which, for given C, becomes a functional of µ, ξ, and 2 L ⟨C⟩. This, in turn, moves the minimum of the generalised free energy without however affecting the general functional form of n(p), given on the left-hand side of Eq. (38). 3

Thermal states
The Gibbs canonical ensemble at temperature T = 1 β is described by the density matrix ρ = e −βH tr[e −βH ], and it corresponds to setting h = 0 and q(p) = βE(p) in Eq. (34). Here E(p) = 4J cos p denotes the single-particle eigenvalue of the energy. The TBA equations can be solved numerically to obtain, for example, the macrosite energy density as a function of the inverse temperature β. In the latter case, comparison with the DMRG-based numerical simulation confirms the TBA prediction, as evident in Fig. 2.
Consistently with the results of Ref. [1], we find that µ tends to 1 2 as the temperature drops to zero. In addition, Fig. 3 clearly shows that the limit is approached exponentially fast in β. It is then convenient to parametrise µ as Plugging this Ansatz into Eq. (39), splitting the integration domain, and changing the integration variables into linear functions of cos p results in the identity In Appendix A it is shown that the first two terms on the left-hand side of Eq. (50) asymptotically behave as  where ζ is the Riemann zeta function and we performed a Sommerfeld expansion of the integrand. On the other hand, we can analytically integrate the third term on the left-hand side of Eq. (50) to obtain the asymptotic expansion which is effective as long as j ≪ Jβ log(Jβ).
Ground state. In the limit T → 0 we can approximate n(p) = [1 + e 4Jβ(cos p−cos k(β)) ] −1 by a characteristic function where k F ≡ lim β→∞ k(β) is the Fermi momentum. In addition, in this limit Eq. (52) reduces to a simple transcendental equation The solution to Eq. (54) exists and is such that cos k F is nonzero: consistently with the numerical observation and with Ansatz (49), we find that the ground state has µ = 1 2. The energy per unit macrosite is now computed by enforcing the replacement (53) in the root density ρ(p) = ρ t n(p), where ρ t is reported in Eq. (27). We finally obtain where Eq. (54) has been used in the last equality. The numeric value of the ground state energy per unit macrosite reads −8J cos k F ≈ −1.73787J and perfectly agrees both with numerical investigations and with the thermodynamic limit of the analytic result that was obtained in our first work [1] focused on finite chains. We can therefore conclude that the thermodynamic ground state matches the actual ground state. Finally, we report the number of particles per unit macrosite, ξ −1 , which is given by Low temperature. For large but finite β the right-hand side of Eq. (52) can not be neglected. Instead, the leading finite-temperature corrections are obtained by setting j = 2 in the latter equation, yielding By expanding k(β) = k F + ∑ n>1 c n (4βJ) −n , for some real coefficients c n , we can then obtain cos k(β) as a function of Fermi momentum k F : Using this in Eq. (49) we see that µ remains exponentially close to 1 2. As shown in Appendix A, ξ(β) and the energy E(β) per unit macrosite can be obtained by carrying out expansions analogous to the one in Eq. (51). In particular we have and the first terms reproducing the thermodynamic limits of the ground state values 2N GS L − 1 and 2E GS L, respectively.
We also report the first orders of the asymptotic low-temperature expansion of the specific heat, which is shown in Fig. 4. We point out that the sub-leading correction, i.e., O((Jβ) −4 ) in the energy and in the magnetisation, is practically irrelevant, somehow manifesting the asymptotic character of the expansion.
High temperature. In the limit of infinite temperature β = 0 one immediately obtains µ = 1 3, n(p) = 3 4, and hence ρ(p) = 1 2π . This is consistent with all traceless local operators having zero expectation value in the infinite-temperature state. An example is given by the conserved charges with single-particle expectation values given in Eq. (18).
At small, but finite β, µ remains close to 1 3, and we are about to show that the deviation scales as β 2 . It is then convenient to parametrise µ as where z(β) is to be computed. Plugging this into Eq. (39) yields π −π dp 2π which can easily be expanded about Jβ = 0. In particular we have confirming the quadratic scaling of the deviation of µ from the infinite-temperature value 1 3. Since, in the limit, the root density remains smooth, we can directly expand it for small Jβ, and find The asymptotic values of the energy per unit macrosite and the specific heat now readily follow (see Fig. 4): (65) Figure 4: The specific heat as a function of the temperature (black solid line). The solid blue and orange lines represent the first orders of the low-temperature (given in Eq. (60)) and the high-temperature (given in Eq. (65)) asymptotic expansions, respectively.
Analogously, the magnetisation per macrosite is given by

Elementary excitations
We define here the elementary excitations of the model. In Ref. [1] we have classified the Bethe states with a set of N momenta that depend on the particle configuration We distinguish two types of elementary excitations: A: creation or annihilation of a momentum along with a Bethe shift of the remaining momenta; B: global "fractional" shift of the momenta.
A close look at the solution (12) to the Bethe equations reveals the following. Excitations of type A are obtained by adding or removing a quantum number I , and choosing I 0 so as to make the change in ϕ negligible in the thermodynamic limit. These are the standard excitations, analogous to the ones in the XXZ model, which are defined through addition or removal of a Bethe-Takahashi quantum number. On the other hand, the excitation of type B can be interpreted as a global shift of quantum numbers I by a fractional amount. It is achieved by changing ϕ (through I 0 ). Figure 5: Examples of elementary excitations over an eigenstate assuming g = 1 both before and after the excitation. In the example above the hole excitation decreases M by 1, whereas the particle excitation does not change M .
Excitations of type A. In terms of the rational quantum numbers, {J 1 , . . . , J N }, excitations of type A correspond to removing or adding a quantum number at an integer distance from the others; A hole excitation corresponds to removing an integer I or, equivalently, the corresponding rational quantum number J . Since N → N − 1, this is necessarily accompanied by a change in the configuration; we propose to update it by removing an elementary particle -see, e.g., Figure 5. As a consequence, M can either remain unchanged or decrease by 1, i.e. ∆M ∈ {0, −1}. More quantitatively, representing this operation by an operator C(p ; b j ) acting on the eigenstate, we have Here, using Bethe equations (12), we find Analogously, a particle excitation corresponds to adding an integer I . The configuration is updated by adding an elementary particle. As a consequence, N → N + 1 and M can either remain unchanged or increase by 1, i.e., ∆M ∈ {0, 1}. Representing this operation by an operator B(p ; b j ) acting on the eigenstate, we have where δ ⃗ p is still given by Eq. (69). Almost always, both before and after the excitation, g = 1, and thus I 0 can be left unchanged. If, before or after the excitation, g ≠ 1, but still g ≪ N , the change in ϕ can always be compensated up to O(N −1 ) by a change in I 0 . If, instead, after the excitation g ∝ N , then the compensation is not always possible and generally such excitation is composite, consisting of an elementary excitation of type B, followed by an elementary excitation of type A.
Excitation of type B. An excitation of type B does not modify the configuration and corresponds to an extensive change of I 0 -see Figure 5. In terms of the rational quantum numbers {J 1 , . . . , J N } we have If we agree on representing this operation by an operator e i 2ϕ L X acting on the eigenstate, we have 4 We note that e i 4π L X can be written in terms of excitations of type A -cf. Eqs (67) and (72): We warn the reader that we conceived these elementary excitations so as to form a basis of the space of excitations and to preserve the macroscopic properties of the excited states. We did not investigate, however, whether the matrix elements of local observables could be nonzero even between states differing in a macroscopically large number of elementary excitations.

Dressed charges
The changes in the expectation values of the local charges (17) after an elementary excitation are of the order O(1). It is customary to interpret such discrepancies as the charges carried by the excitation. In interacting systems these values depend on the state and are usually called "dressed charges". For example, after a hole excitation of type A, generated by operator C(p ; b j ), the integral of motion with a single-particle eigenvalue q(p) changes as and hence, using Eq. (69), we find Here, we have extended the notation ⟨•⟩, for the average of an operator in the macrostate described by ρ(p), to functions of momenta: The dressed momentum ℘ thus reads Comparing this with Eq. (26), we obtain the standard relation between the dressed momentum derivative and the total root density, i.e., 2πρ t = ∂ k ℘(k; b j ). Analogously, the dressed energy ε ≡ E dr of the excitation corresponds to q(p) = 4J cos p and is given by In the specific case of a thermal state at temperature β −1 the expectation value of any charge that is odd under spatial reflections zeroes, so we have where ⟨•⟩ β signifies that the average is taken in a thermal macrostate. Contrary to the XXZ model, the dressed energy of an elementary excitation does not match β −1 log 1 n(p) −1 , which is unhappily called "dressed energy" as well 6 . In a thermal state, however, their derivatives with respect to the momentum coincide.
For the sake of completeness, we finally report the derivative of the dressed charge with respect to the momentum (rapidity), which plays a key role in the thermodynamic Bethe Ansatz. It reads

Hydrodynamics with a minimal set of charges
This section develops the tools necessary to describe time evolution of inhomogeneous states in the folded XXZ model. A theory, now known as "generalised hydrodynamics", was proposed in the seminal papers [5,6] to describe the late-time behaviour of the expectation values of local observables in integrable systems that evolve in time in the presence of inhomogeneities. Specifically, it was assumed that one can describe the time evolution of a large class of inhomogeneous states by lifting the root densities that characterise stationary states into functions of space and time (in addition to the momentum). In this way, the concept of local equilibrium is generalised to integrable systems, where the stationary states are characterised by infinitely many conservation laws. The structure underlying local equilibrium in noninteracting spin chains was clarified in Ref. [25]. There, it was shown that generalised hydrodynamics can be interpreted as a phasespace description of the dynamics in a subspace of the Hilbert space that is invariant under time evolution. When scrutinised locally, the states belonging to that subspace are quasistationary -for that reason they have been called "locally quasistationary states" [26] -and the generalised hydrodynamic equation of motion is simply the projection of the Schrödinger equation onto the invariant subspace.
In interacting integrable systems the structure is less transparent and the aforementioned invariant subspaces have not yet been identified. The theory is therefore still based on the weaker assumption that, after a long-enough time, the unspecified degrees of freedom that can not be described by space-time dependent root densities become irrelevant. This corresponds to identifying the leading order(s) of an asymptotic expansion in the degree of the inhomogeneity. To that aim, we first need to define the charge densities and the currents of the integrable hierarchy introduced in Section 1.4. We will then develop first-order generalised hydrodynamics in the dual folded XXZ model.

Charge densities
A local charge Q that commutes with the shift operator e iP , where P is the operator of the total momentum, can always be represented as a sum of localised operators that differ only in their position: Q = ∑ q , with q +1 = e iP q e −iP . The operators q are referred to as charge densities, but their definition is not at all unique. For example, if q is a charge density, the operator q + z − z −1 is as well, whatever z is. Such an ambiguity can be used to simplify either the definition of the operator or its dynamical equations, the latter point of view being taken, for example, in Ref. [25]. In both cases, it is usually convenient to choose a density with as many symmetries as possible. Simple symmetries that can be easily enforced are generated by operators with equally spaced eigenvalues. Consider, for example, the total magnetisation S z . The eigenvalues of L 2 + S z are integers, thus e iϕ( L 2 +S z ) is 2π-periodic in ϕ. This implies that the averaged densitȳ commutes with S z 7 . Since a rotation does not change the range of q, the local densityq has the same locality properties as q, but is, in addition, U (1)-invariant. As a matter of fact, a similar conclusion can be drawn even if we replace S z with a less trivial charge with equally-spaced spectrum, say A. The averaged charge still commutes with A, but its range might increase. Nevertheless, the range can not become arbitrarily large, provided that the local terms of A decay exponentially with their range. This follows from a result of Ref. [27], which showed that, if q is a local operator with range r, the operator e iϕA qe −iϕA can be approximated, with exponential accuracy, by an operator with range r + 2v LR T , where v LR is the Lieb-Robinson velocity of A. Specifically, the error of this approximation, i.e., the part of operator e iϕA qe −iϕA with a range larger than r + 2v LR T , is exponentially small in v LR (T − ϕ). This bound can be readily adapted to the averaged charge density by replacing ϕ with its maximum, namely 2π.
In our specific case, the invariance of the configuration can be expressed as the conservation of infinitely many quasilocal charges with equally spaced spectrum (including S z , S z st , and M). It is therefore natural to define charge densities that preserve the configuration. Remarkably, the local densities (19) are already of such form. This holds even if one redefines them by adding a local operator z ± n, that preserves the configuration and satisfies ∑ L =1 z ± n, = 0. Concerning the diagonal charges (the ones commuting with any σ z ), a sensible choice is to enforce their densities to be diagonal as well. In particular, we will use s z = σ z 2 and as diagonal local densities of S z and M, respectively. Finally, we mention that the expectation value of m can be rewritten in terms of the emptiness formation probability P (n) = ⟨∏ +n−1 j= 1−σ z j 2 ⟩ as follows:

Currents
A fundamental role in generalised hydrodynamics is played by the continuity equations that relate the local charge densities to the corresponding currents. In a spin chain, they are obtained by applying the Heisenberg equation to the charge restricted to a half-open interval ∆X = a[n − , n + ), where n − and n + are the integers denoting the boundary sites, while a is the lattice spacing. In particular, the equation reads ∂ t ∑ ∈∆X a q (t) = i[H, ∑ ∈∆X a q (t)]; since the commutator acts nontrivially only around the boundaries of the region, it can also be written as If the region ∆X includes a mesoscopically large number of sites ((x + − x − ) a ≫ 1) and the local properties of the state vary in a sufficiently smooth way on a larger scale, the continuity equation can also be expressed in a differential form. Here, the expectation values are interpolated by smooth functions ⟨q ⟩ (t) = a ⟨q⟩ (x, t) and ⟨j ⟩ (t) = a ⟨j⟩ (x, t), with x = a , and, according to the trapezoidal rule and the definition of the derivative, we have Importantly, if the ranges of q and j are both small with respect to ∆X a, one can ignore the inhomogeneity of the state in the region where q and j act nontrivially, so that the state can be effectively replaced by a homogeneous one. In our specific case, there exist local and quasilocal conservation laws that break one-site shift invariance, therefore we can not generally expect the smoothness hypothesis to hold on a mesoscopic scale if is associated with a chain site. The natural choice for the constituents of the unit cell are instead macrosites -the index in Eqs (85) and (86) should then refer to a pair of neighbouring spins. Technically speaking, this is a consequence of the fact that the momentum entering the Bethe equations (21) generates two-site translations. For the sake of clarity, let us introduce the notations q ′ = q 2 ′ −1 + q 2 ′ for the macrosite charge density and  ′ for the macrosite current. For example, using definitions (19), the energy current reads whereas the current of Q − 1 is given by Note that the macrosite current equals the spin-site current evaluated at odd sites  ′ = j 2 ′ −1 . 8 The continuity equation for ∆X = { ′ } is represented schematically in Fig. 6. Importantly, having defined the charge densities so as to preserve the configuration, the corresponding currents will commute with S z , M, as well as with other charges related to the conservation of the configuration.
. . . To the best of our knowledge, the staggered spin along the z-axis is the unique strictly local charge that breaks one-site shift invariance. Setting the local density equal to s z st, ′ = (σ z 2 ′ − σ z 2 ′ −1 ) 2, the staggered spin current reads We now turn to the current of the quasilocal charge M. As evident from Eq. (16), we can choose the quasilocal operator as the charge density.

In a macrostate
In a macrostate characterised by a root density ρ(p), the expectation value of a charge Q is computed according to Eq. (24), which can be written as Here, q(p) is the single-particle value of charge Q and vac⟩ = ↓↓ . . . ↓⟩ is the Bethe Ansatz reference state. On the left-hand side of Eq. (92) we subtract the vacuum expectation value, since the right-hand side is always zero in the vacuum state (the latter contains no momenta, so ρ (vac) (p) = 0). The subtraction is allowed, because the charge density expectation value is defined up to an additive constant.
Remarkably, we can also relate the expectation value of the corresponding current to the root density. To that aim, it is convenient to take a step backward and consider first a finite system with L spins. Reference [28] showed that the expectation value of the current in a Bethe state can be expressed in terms of the Gaudin matrix G, which is a quantity discernible directly from the Bethe equations, as follows: Here, (•) ′ denotes the derivative, and N is fixed because both the charge density and the current commute with the total magnetisation along the z-axis, S z . The elements of the Gaudin matrix are given by Analogously to what was said for N , the partial derivative in Eq. (94) is taken at fixed M . Indeed, both the charge density and the current commute with M as well, and hence the continuity equation (85) can be enforced at fixed M . Comparing Eqs (21), (78), and (94), we see that the Gaudin matrix corresponds to the Jacobian of the transformation p ↦ ℘, up to a multiplicative factor. In the end we find Using E(p) = 4J cos p, the expectation value of the current then reads Relation (96) is expected to hold for the currents of all local charges Q ± n , 9 and we have checked it for J ± 1 (see Eqs (87) and (88)) by means of exact diagonalisation, up to L = 14 sites. In the thermodynamic limit Eq. (96) becomes providing the sought-after relation between the expectation value of the current and the root density. Here, v(p) = ∂ p ε(p) ∂ p ℘(p) is the effective velocity of quasiparticles, explicitly expressible as where we used the dressing equations (78) In the noninteracting sector, where M = 0, J d [Z n ] are equivalent to local operators. This property is, however, immediately lost in the presence of interactions. 10

In a locally quasistationary state
Reference [26] named "locally quasistationary state" (LQSS) the inhomogeneous macrostate that captures the expectation values of local operators in an inhomogeneous out-of-equilibrium state at asymptotically large times, after the fastest degrees of freedom have stopped affecting the dynamics of local observables. In particular, at each point in space an LQSS is locally equivalent to some stationary state. A route to a formal definition of locally quasistationary states is based on the identification of the smallest families C of operators, including both the set of conserved charges and their densities, that are closed under time evolution [25]. A locally quasistationary state could than be represented by a density matrix ρ = e W tr[e W ], with W ∈ C. A basic property of the family C is that, for any conserved operator O ∈ C, its charge density is defined in such a way that J d [O] (which generalises functional (99) to the conserved operators belonging to C) is in C as well. Notwithstanding the simplicity of Eq. (99) suggesting the potential to identify a family C of such operators in the folded XXZ model, we leave this problem to future investigation. Here, we stick to the heuristic method that has been used to obtain the first-order generalised hydrodynamics (GHD) equation in interacting integrable systems.
To that aim, we consider the continuity equation (85) in the limit described by Eq. (86): The local equivalence of an LQSS to a stationary state indicates the possibility to lift the root density to a function of space and time. Let us then assume that the expectation values in the 10 Note that, as discussed in the first part of our work [1], for finite L the local charges Zn are not functionally independent. Indeed, we have where Π M is the projector on the subspace with M = M (see Ref. [29] for some observations concerning the independence of integrals of motion in the presence of such a functional dependence). This is a further manifestation of the incompleteness of the charges Q n . LQSS characterised by the density matrix e W tr[e W ], with W ∈ C, are completely determined by a root density depending on space and time. Even if we do not really know how to compute expectation values in such an inhomogeneous macrostate, if the range of the operator is much smaller than the typical length of the inhomogeneity, we can simply compute its expectation value in the macrostate characterised by ρ x,t (p), by treating x and t as external parameters. Under such a low-inhomogeneity assumption, using Eqs (92) and (97), we finally obtain The higher order derivatives signify the corrections due to the modulation of the root density on length scales smaller than the typical size of the inhomogeneity. Since the single-particle values q ± n (p) of the charges Q ± n form a basis in the space of square-integrable functions of the momentum p ∈ [−π, π), we end up with the fundamental equation of generalised hydrodynam- By symmetry, the density of holes ρ h (p) is expected to satisfy the same equation, which then implies where we used that, in our specific case, ρ t = ρ(p) + ρ h (p) does not depend on the rapidity. At first sight the latter equation might look wrong, since it relates the rapidity-independent function ρ t to a function that, instead, seemingly depends on the rapidity. This is only apparent: plugging Eq. (98) into the continuity equation (103) and using 2πρ t = 1 + µ ⟨1⟩ (cf. Eq. (26)), we obtain where we used notation (77). By virtue of Eq. (102), this can be recast into a continuity equation for µ x,t : Equations (102) and (103) can also be used to obtain the dynamical equation satisfied by the filling function n x,t (p) = ρ x,t (p) ρ t;x,t , which is given by The O(∂ 2 x ) corrections become irrelevant in a particular scaling limit, often refered to as "Euler scale" or "ballistic scaling limit". In this limit one can truncate generalised hydrodynamics at the first order. The description of diffusive, sub-diffusive, and Tracy-Widom scaling behaviours requires the development of higher-order generalised hydrodynamics [25,30,31].
Finally, we note 2 L ⟨M⟩ x,t = 2πρ t;x,t − 1, therefore Eq. (104) can be used to infer the expectation value of the corresponding current  ′ [m], given by Eq. (91), in a macro-state: In other words, the diagonal part of the total current reads An attentive reader might have noticed that the picture developed so far does not include the staggered magnetisation. We conjecture that the staggered magnetisation -as well as the other quasilocal charges whose eigenvalues describe the configuration -satisfies the same GHD equation as µ = ⟨m⟩ ⟨1⟩, namely In the next section we study the first-order GHD equation, expressed by Eqs (105) and (106), in the inhomogeneous setting where two local macrostates are joined together.

Junction of two local macrostates
A nonequilibrium setting that has been intensively studied for its potential to elucidate transport properties of both classical and quantum many-body systems is the junction of two pieces of materials kept at different temperatures [32,33,34,35,36,37,38,39]. The system of interest is typically sandwiched between the materials at the junction. In this case the materials act as thermal baths, and such a setting is a paradigm of a non-isolated system. If, however, the two pieces are in direct contact, one can take both of them as constituents of a larger isolated system. In the latter, however big the pieces of materials are, they will experience a nontrivial evolution over time. In some cases the evolution can be captured by assuming local thermalisation, where the state of the system is still characterised by a temperature, but the latter depends on space and time. An important exception are the integrable systems, where local thermal equilibrium states do not form invariant subspaces. In other words, the nonequilibrium state can not be characterised solely by the time evolution of the local temperature. Instead, the full set of generalised chemical potentials (or generalised temperatures) is required for the determination of such an out-of-equilibrium state.
We consider here the junction of two local macrostates in the folded XXZ model, for example, two thermal states. Besides its physical relevance, this scenario carries some mathematical simplifications. Firstly, the boundary conditions, i.e., n(p) and µ in the initial state, are functions of the ray ζ = x t, so the solution to the GHD equations is expected to depend on x and t only through their ratio. Secondly, the corrections to Eqs (105) and (106) are irrelevant in the ballistic (B) scaling limit t → ∞, at fixed ζ = x t, so the GHD equations reduce to where we have defined −π dp n ζ (p)E ′ (p) ∫ π −π dp n ζ (p) . (111) The ballistic scaling limit is supposed to exist whenever the limits exist.

On the solution to the GHD equation
The initial conditions (112) fix the values of n ζ (p) and µ ζ in the limits ζ → ±∞. Thus, in order for the solution to be unique, the equations ζ = v ζ (p) (for each given momentum p) and ζ = V ζ should only have one solution. In this section such uniqueness conditions will just be assumed to be satisfied. The solution to the second of the GHD equations (110), namely, the equation for µ ζ , is a piecewise constant function with one discontinuity (for the uniqueness assumption). It is then convenient to define the auxiliary velocities and the auxiliary GHD equations with the same boundary conditions as in the original problem. The sign in the superscript or subscript refers to the region, with respect to the discontinuity in µ ζ , where the equations or quantities are defined. The full solution to Eq. (110), with a single discontinuity, can be constructed from the solutions to the auxiliary GHD equations (114) by joining n − ζ (p) and n + ζ (p) at the ray of the discontinuity in µ ζ , i.e., the ray that solves both ζ = V − ζ and ζ = V + ζ . The latter two equations should therefore have the same solution. To find it, we note that a change in µ + does not affect V − ζ . This implies that the common solution of equations ζ = V − ζ and ζ = V + ζ is, in fact, independent of µ + . We can then specialise the calculation to the noninteracting case µ + = 0, in which the effective velocity is simply v + (p) = E ′ (p) = −4J sin p, so that n + ζ (p) = n sgn(4J sin p+ζ) (p). The equation ζ = V + ζ can then be rewritten as 0 = f (ζ; n + (p), n − (p)) ∶= π −π dp n sgn(4J sin p+ζ) (p) sin p + ζ 4J .
Since f (ζ; n + (p), n − (p)) is a monotonous function of ζ, and lim ζ→±∞ f (ζ; n + (p), n − (p)) = ±∞, equation f (ζ; n + (p), n − (p)) = 0 has a single solution ζ = ζ d . This simple observation allows us to conclude that, assuming a single discontinuity in µ, the ray ζ d of the discontinuity is the unique zero of f (ζ; n + (p), n − (p)). Again assuming uniqueness, the first of the GHD equations (110) is solved by where n ± ζ (p) are the filling functions on each side of the ray of the discontinuity. They separately solve the first of the auxiliary GHD equations (114) and read In addition to this solution, the following two facts should be kept in mind: 1. The maximal (minimal) ray at which the filling function undergoes a transition from n − (p) to n + (p) corresponds to the right (left) edge of the light cone, emanating at the junction of the two local macrostates. It is given by In the ballistic scaling limit, the expectation value of any local observable outside the light cone is constant and equal to its initial value.
2. At the ray ζ = ζ d of the discontinuity we have where we have used ζ d = V ± ζ d and the definitions (113) of the auxiliary velocities. Since for all momenta p. In particular, for momenta p that satisfy ζ d = E ′ (p) = −4J sin p, the filling function changes its value from n − (p) to n + (p) at the ray ζ d . This will be used in the next section, where we discuss the discontinuities in the expectation values of local charges. Figure 7: When the ray ζ crosses v − ζ (k), the value of the filling function at momentum k jumps from n − (k) to n + (k). When ζ d is crossed, µ − changes discontinuously into µ + , and v − ζ (k), at k satisfying E ′ (k) = ζ d , continuously into v + ζ (k). Sweeping the interval ζ ∈ (−∞, ∞), the lefthand side state n − completely transforms into n + . The data shown corresponds to the junction of two thermal states with β − = 0.25 and β + = 2.75. Clearly visible are the asymmetric light cone between the rays ζ − lc ≈ −3.0140J and ζ + lc ≈ 2.9400J, computed according to formula (118), as well as the ray of discontinuity ζ d ≈ 0.2416J. The sketched filling functions are purely qualitative.
To summarise, the picture that unfolds is as follows (see Fig. 7). When the ray ζ sweeps from ζ = ζ − lc to ζ = ζ d , the filling function transforms from n − (p) to n + (p) at the momenta , the change being, in fact, continuous. By "continuous" we mean that, at the momenta p for which the transition of the filling function from n − (p) to n + (p) occurs at ζ d (these momenta satisfy ζ d = E ′ (p)), the velocities coincide:

Non-ballistic behaviour
Let us assume that the filling functions of the boundary conditions are smooth; in this way any non-analytic behaviour can be traced back to the "domain-wall" structure of the initial state. As it generally happens, the discontinuity in the filling functions that is caused by the discontinuity of the initial state does not directly produce non-analytic behaviour in the profiles of the local observables (except for the neighbourhoods of the light cones). Indeed, if µ + = µ − , the ballistic profiles are smooth inside the light cone 11 . On the other hand, a discontinuity in µ at the initial time immediately spoils the smoothness of the profiles. At first sight, this bears resemblance to the effect of a discontinuous sign of magnetisation in the XXZ model [40]. However, in that case, in the ballistic scaling limit the operators that are even under the change in the sign of S z are not affected by the discontinuity. In contrast, in our case µ + ≠ µ − makes the ballistic profile of any generic local charge and current discontinuous. In other words, there is no global symmetry (independent of the initial state) that protects the smoothness of the profiles. In the following we quantify the discontinuities of charges and currents, which turn out to be deducible a priori. Denoting f ζ − d = lim ζ↗ζ d f ζ and f ζ + d = lim ζ↘ζ d f ζ , the left and right limits of the expectation values of charges and currents at the discontinuity are given by where n ± ζ d (p) and v ± ζ d (p) readily follow from ζ d = V ± ζ d and Eqs (113), (117), and (120): Remarkably, these expressions, along with Eq. (115), are written solely in terms of the boundary conditions. We also note that the existence of two rays ζ ± d , described by the same filling function but different µ, rules out the possibility to interpret the locally quasistationary state as a collection of local macrostates, since in the latter µ is fixed by Eq. (39). In particular, since thermal states at different temperatures have different values of µ (see Fig. 3), this observation applies also to the junction of thermal states.
Remarkably, the ratio of expectation values of any two local charge densities q a and q b is continuous π −π dp n sgn(ζ d +4J sin p) (p)q a (p) ∫ π −π dp n sgn(ζ d +4J sin p) (p)q b (p) .
Although this does not extend to the expectation values of currents, it holds true for the shifted currents [q] + q ⋆ q − 1 , where q a ⋆ q b denotes the charge with a single-particle value q a (p)q b (p). Indeed, we have Similar calculations can be performed for the expectation value of m. Here, the left and right limits are given by while the discontinuity reads π −π dp n sgn(ζ d +4J sin p) (p) (2π − µ + ∫ π −π dp n sgn(ζ d +4J sin p) (p))(2π − µ − ∫ π −π dp n sgn(ζ d +4J sin p) (p)) . (129) Finally, we point out that the GHD equation for µ can be replaced by the additional boundary condition n ζ d (p) = n sgn(ζ d +4J sin p) (p) for the filling function in the two independent regions ζ < ζ d and ζ > ζ d , in which µ is constant.

Entanglement entropy
The expectation values of charges and currents are not the only physical quantities that exhibit non-ballistic behaviour. At large time, the entanglement entropy of a subsystem per unit length is supposed to approach the Yang-Yang entropy density in the emerging locally quasistationary state [41,42]. In our case the Yang-Yang entropy is reported in Eq. (33). 12 For the sake of simplicity we assume that the initial state has zero staggered magnetisation. Then, around the discontinuity, the Yang-Yang entropy reads H(2µ ± ) ∫ π −π dp n sgn(ζ d +4J sin p) (p) + ∫ π −π dp H(n sgn(ζ d +4J sin p) (p)) 2π − µ ± ∫ π −π dp n sgn(ζ d +4J sin p) (p) (130) 12 We are assuming that the locally quasistationary state is locally equivalent to a quasilocal macrostate defined in Section 2.2.2. If that is not the case, our estimation of the Yang-Yang entropy becomes an upper bound.
The resulting discontinuity is computed as −π dp H(n sgn(ζ d +4J sin p) (p)) ∫ π −π dp n sgn(ζ d +4J sin p) (p) (2π − µ + ∫ π −π dp n sgn(ζ d +4J sin p) (p))(2π − µ − ∫ π −π dp n sgn(ζ d +4J sin p) (p)) As a case study we could consider the junction of the ground state, in which µ − = 1 2, n − (p) = θ( p − k F ), with a pure state in which any traceless charge has zero expectation value: the late-time behaviour on the right of the light cone would then be described by the infinite temperature state, where µ + = 1 3, n + (p) = 3 4. Calculation then yields ζ d ≈ −0.272252J, ∫ π −π dp n sgn(ζ d +4J sin p) (p) ≈ 4.16302, and ∫ π −π dp H(n sgn(ζ d +4J sin p) (p)) ≈ 1.72832. The entropy discontinuity is then given by s YY;ζ + d − s YY;ζ − d ≈ 0.482976. To the best of our knowledge, this is the first model where a discontinuity is predicted in the ballistic scaling limit of the entanglement entropy. Figure 8 shows a comparison between the generalised-hydrodynamics predictions and the numerical simulations based on tDMRG algorithms [43,44,45,46]. We simulated the dynamics in a spin chain with open boundary conditions; the boundary effects are mitigated by choosing initial states that are locally stationary. We consider two scenarios: in the first scenario two thermal states are prepared at temperatures T L = 10 and T R = 1 on the left half and on the right half of the spin chain, respectively, where we set J = 1. In particular, denoting H L = ∑ L 2−2 =1 q + 1, and H R = ∑ L−2 =L 2+1 q + 1, , with q + 1, given in Eq. (19), the initial state reads

Comparison with numerical simulations
It is prepared by performing imaginary-time evolution of the infinite-temperature state, using the ancilla tDMRG annealing technique [47,46]. Once thermalised, the two halves of the system are coupled and evolved in real time with the full folded Hamiltonian. The boundary conditions for the GHD continuity equations read where the energy shifts w ± are computed by solving Eqs (38) and (39) at zero staggered magnetisation m z st = 0. In the second scenario, the initial state is prepared so as to have a nonzero staggered magnetisation in the left half of the spin chain. It reads where we have chosen h L = 4. To obtain the new left-hand side boundary condition for the GHD continuity equations, one has to compute the staggered magnetisation per macrosite according to Eq. (37). This only affects w L , through Eqs (38) and (39). For more details on numerical simulations see Appendix B. , the GHD prediction for the macrosite energy density has been divided by 2, in order to obtain the average energy per spin site. In panel (d), the magnetisation on the left-hand side oscillates between the solid lines that border the shaded region (it is staggered as a result of the initial condition). The shaded region expands in time, as the discontinuity propagates towards the right-hand side with velocity computed from Eq. (115). In all plots, the profiles are centered at the geometrical center of the spin chain used in the simulation.

Conclusion
In this work we have considered the dual folded XXZ model in the thermodynamic limit. We have developed the thermodynamic Bethe Ansatz that enables description of generalised Gibbs ensembles constructed from the strictly local charges and a special quasilocal conservation law related to the number of domain walls in a particle configuration. We have studied thermal states in some detail, exhibiting the low-and high-temperature asymptotic expansions of the energy, the magnetisation, and the specific heat. We have then moved our attention to the time evolution of inhomogeneous states prepared so as to be locally equivalent to local macrostates on large space and time scales. Such systems can be successfully investigated within the framework of generalised hydrodynamics. We have adapted this theory to the folded XXZ model and applied it to the crucial setting in which two macrostates (e.g., thermal states) are joined together at time t = 0. We found that the profiles of the expectation values of charge densities and currents are discontinuous in the ballistic scaling limit. To the best of our knowledge, there are no global symmetries that would protect some charges from developing discontinuities, contrary to the analogous scenario in the XXZ model.
Besides solving the first-order GHD equations numerically, we have also obtained analytical predictions for the position of the discontinuity and for the macrostates that describe the expectation values of the local observables in its neighbourhoods. We have tested our theoretical predictions against extensive tDMRG simulations of time evolution. The agreement was always fair, the discrepancies being reasonably explainable as finite-size or finite-time effects, and as numerical inaccuracies.
We warn the reader that, despite being sufficient for our purposes, the thermodynamic Bethe Ansatz that we developed (and, in turn, the generalised hydrodynamic theory) is not complete, as we did not include all quasilocal charges that complete the characterisation of the configuration space. Some of these deficiencies could be filled rather easily. However, the solution would nevertheless remain deficient, as it would not take into account the nonabelian integrability of the folded XXZ model. The interest in the latter structure comes from its striking effects on the dynamics at intermediate times, therefore this aspect deserves further investigations, especially in the presence of interactions. Finally, despite non-ballistic behaviour being manifest in the inhomogeneous setting, we only worked out generalised hydrodynamics at the first order (ballistic scaling limit). In this regard we wonder whether the folded XXZ model is simple enough to make a step forward in the development of higher-order generalised hydrodynamics in interacting integrable systems.
Expressions on the right-hand side are the Lagrange remainders of the Taylor polynomials. For some 0 < ν 2 < 1 we thus have where, in the second equality, we used Eq. (135) to rewrite the integral over the interval [0, 1 − z − z 0 ]. The first two terms can easily be bounded with terms that are exponentially small in β. In particular, we have The last but one term in Eq. (136) is instead a power law correction, namely, where ζ is the Riemann zeta function. For any 0 < z 0 < 1 − z all of the correction terms are bounded by a power law. We therefore find the asymptotic expansion where the integral in the last term of Eq. (136) has been worked out. An analogous calculation allows us to work out the asymptotic expansion of expressions of the type ⟨f (cos p)⟩ (2πρ t ) = ∫ dp 2π n(p)f (cos p), for a given smooth function f . Indeed, we have ⟨f (cos p)⟩ 2πρ t = π −π dp 2π where and the last step is exponentially accurate in β. The lowest-order coefficients of the expansion of g(x) about x = cos k F , for f (x) = 1 and f (x) = x, are shown in Table 1.

B DMRG simulations
Numerical calculations are performed using C++ ITensor library [48]. We use the ancillabased DMRG technique [47,46] to prepare initial state with a given inverse temperature β. We fix the maximal bond dimension to be N max = 600 and a truncation, such that the sum of the discarded Schmidt values is smaller than 10 −13 . The imaginary-time step in the preparation of the initial state and the real-time step in the time evolution of the state are both set to δt = δβ = 0.01J −1 . The code and the animated evolution of the energy profile are freely accessible (see also Fig. 9). 13,14,15 Figure 9: QR codes that link to animated version of the Fig. 8. Alternatively, follow the links in the footnotes below.

B.1 Preparing the initial state
In order to use algorithms based on the matrix product state representation, the density matrix is vectorised as where the infinite temperature state is defined as a product of maximally entangled physical and ancilla spins: Here, indices j and denote the physical and the ancilla spin, respectively. 13 The code: https://github.com/kbidzhiev/Folded_XXZ 14 Evolution of the energy profile: https://github.com/kbidzhiev/Folded_XXZ/raw/master/Pictures/ Energy_profile.gif 15 Evolution of the rescaled energy profile: https://github.com/kbidzhiev/Folded_XXZ/raw/master/ Pictures/Energy_profile_rescaled.gif  . Trotter gates act on the wavefunction and gradually decrease the temperature from T = ∞ (β = 0) to T L (β − ). Panel (b) shows the energy profile for a state with temperature T L (T R ) on the left (right) half of the system. We keep acting with Trotter gates on the right-hand side, until the temperature T R is reached.
In numerical simulations we consider the Hamiltonian (8) with open boundary conditions. It is written as a sum of three parts H = H 1 + H 2 + H 3 , such that local terms in each mutually commute. This allows to approximate the evolution operator by Trotter gates [49,50,51]. In the present work we use second order Trotter-Suzuki approximation e −δt(H 1 +H 2 +H 3 ) = e − δt 2 H 1 e − δt 2 H 2 e −δtH 3 e − δt 2 H 2 e − δt 2 H 1 + O(δt 3 ) , which can be factorised in terms of Trotter gates exp i J 2 δt(σ x σ x +2 + σ y σ y +2 )(1 − σ z +1 ) that act on the physical sites only, the ancilla spins remaining untouched. This helps avoid additional truncation.
Starting from the infinite temperature state ψ(0)⟩, we gradually cool the system by imaginary-time evolution until it reaches a state ψ(β − )⟩, with β − = 1 T L . For a typical example, see panel (a) of Fig. 10. The further cooling is provided by Trotter gates that act only on the right half of the system, until the right subsystem reaches the temperature β + = 1 T R , as shown in panel (b) of Fig. 10. This protocol creates initial correlations between the left and the right subsystems. Alternatively, one could zero the initial correlations between the two parts by acting with the Trotter gates on both subsystems separately, omitting the three central Trotter gates that overlap with the junction between the two subsystems. The latter protocol, however, introduces larger oscillatory finite-size effects. The energy profile that results from this last protocol is presented in panel (a) of Fig. 8.  Figure 11: Panel (a) shows the energy density as a function of Jβ. Discrepancy between the analytical and numerical result for the energy ⟨H total ⟩ is a consequence of the boundary Friedel-like oscillations. To avoid the inaccuracies due to these oscillations, we average the energy in the vicinity of the centre, i.e., on sites −L 4 ≤ i < L 4, as shown in panel (b) (the corresponding average is denoted by ⟨H centre ⟩). The effect of oscillations is larger at low temperatures and reaches the maximum in the ground state -see panel (d). Panel (c) shows the energy density for different system sizes. Already for relatively small systems of ∼ 40 spins the results perfectly agree with the analytical predicitions. Panel (d) shows the energy profile in the ground state (β → ∞). The Friedel-like oscillations reach their maximal value and slowly vanish as ∼ 1 √ x, where x is the distance from the boundary. For small finite sizes the envelope of the oscillations is described by the "chord length" formula. boundary conditions. In the vicinity of the centre and boundaries one can observe Friedellike oscillations, which increase as the temperature drops to zero. The oscillations reach the maximal amplitude in the ground state, as shown in panel (b) of Fig. 11. They decay as the inverse square root of the "chord length", i.e.,

B.2 Energy density
where x denotes the distance from the boundary of the system, as measured in the number of spin sites. The energy front in the left (right) half of the system propagates with velocity ζ ± lc ≈ 3, measured in macrosites, or ≈ 6, measured in terms of spin sites. This gives rise to a "light cone" shown in Figs 7 and 8. The energy smoothly tends to the initial-state values near the light cone edges. The discontinuity in the numerically obtained profile is relatively small on the achievable time scales.