What is an integrable quench?

Inspired by classical results in integrable boundary quantum field theory, we propose a definition of integrable initial states for quantum quenches in lattice models. They are defined as the states which are annihilated by all local conserved charges that are odd under space reflection. We show that this class includes the states which can be related to integrable boundary conditions in an appropriate rotated channel, in loose analogy with the picture in quantum field theory. Furthermore, we provide an efficient method to test integrability of given initial states. We revisit the recent literature of global quenches in several models and show that, in all of the cases where closed-form analytical results could be obtained, the initial state is integrable according to our definition. In the prototypical example of the XXZ spin-s chains we show that integrable states include two-site product states but also larger families of matrix product states with arbitrary bond dimension. We argue that our results could be practically useful for the study of quantum quenches in generic integrable models.


Introduction
The notion of solvability associated to integrable models usually refers to the possibility of diagonalizing analytically the corresponding Hamiltonian [1]. Although this is a remarkable property, comparison with experiments often requires to go beyond the simple computation of the system's eigenspectrum and to provide predictions for more general physical observables. Indeed, a tremendous effort has been made in the past decade to compute several nontrivial ground-state and thermal properties in different models, significantly improving our understanding of the underlying mathematical structures [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16].
More recently, integrable models offered new challenges to the community, as an increasing attention was devoted to the study of quantum quenches [17]. In this framework, one is interested in the unitary time evolution from an assigned initial state. It is evident that this problem is generally more complicated than those arising at thermal equilibrium, mainly because the initial state itself might be chosen arbitrarily, with an exponentially large number of degrees of freedom. This freedom hardly combines with the intrinsic rigidity of integrability, and one might legitimately wonder whether quench problems are within the reach of integrability-based techniques at all [18]. On the other hand, nearly ideal, out-ofequilibrium integrable models can now be realized experimentally in cold-atomic laboratories [19][20][21][22], elevating the relevance of these questions beyond a purely academic curiosity (see also [23] for a volume of reviews on integrable models out of equilibrium).
Among the most elementary quantities appearing in the study of quantum quenches are the overlaps between the initial state and the eigenstates of the model. These are needed as intermediate building blocks for several calculations both at finite and infinite times. Ever since the early stages of the literature on quantum quenches, several works tackled the problem of their computation, with the first studies dating back to almost ten years ago [24,25]. However, until recently the problem appeared too hard to be attacked, strongly limiting our possibility of analytical calculations.
Paralleling these studies, analytical results for the post-quench time evolution were also obtained in cases where the overlaps are not known. In the works [59,60] an exact computation of the Loschmidt echo was performed in the XXZ chain, starting from arbitrary families of two-site product states. These results could be derived by means of a construction relating the initial state to an integrable boundary condition in an appropriate rotated channel, where the space and time directions are exchanged. In fact, the connection between quantum quenches in one spatial dimension and boundary integrable quantum field theory (QFT) in the rotated channel has been known and exploited for a long time both in the conformal [17,[61][62][63][64] and massive cases [48,[65][66][67][68][69].
This connection was in particular beautifully illustrated in the classical work of Ghoshal and Zamolodchikov [65]. Here, integrability of the boundary field theory was defined by the existence of an infinite number of conserved operators (or charges) which persist after the addition of a boundary term to the bulk Hamiltonian. Remarkably, the conditions on this boundary term to preserve integrability were explicitly translated into a constraint for the boundary state in the corresponding rotated picture: the latter has to be annihilated by an appropriately chosen (infinite) subset of the bulk conserved charges.
Inspired by these classical works, we propose a definition of integrable states for quantum quenches in lattice integrable systems. We identify integrable states as those which are annihilated by all local conserved charges of the Hamiltonian that are odd under space reflection. This definition naturally extends to those models defined on the continuum that are scaling limit of lattice ones. We introduce an efficient method to test integrability of given initial states. By means of the latter, we show that in all of the cases recently considered where closed-form analytical results could be obtained, the initial state is integrable according to our definition.
Going further, we prove that integrable states include the ones which can be related to integrable boundaries in the rotated channel, where the space and time directions are exchanged [60]. This result, which completes the analogy with the picture in QFT, is highly non-trivial, as the lattice model does not display Lorentz invariance. We devote a detailed analysis to the prototypical case of XXZ spin-s chains and show that integrable states include two-site product states but also larger families of matrix product states (MPSs) with arbitrary bond-dimension. A particular case of the latter are the states studied in [35,39,40].
It is important to stress that our findings should not be interpreted as a no-go theorem for obtaining exact results for non-integrable initial states. For example, an exact overlap formula was computed in [25] for the special case of the domain-wall state, which strictly speaking is not integrable. Further analytical results were obtained recently for quenches from the same state, for example regarding the computation of the return probability [70], and in the context of spin transport [71]. However, these problems are inherently inhomogeneous in space, and in the present work we focus on the homogeneous global quenches.
Regarding global problems, relevant examples concern the computation of the postquench rapidity distribution functions of the quasi-particles in XXZ spin-s chains. Indeed, the latter can be computed exactly in many (non-integrable) cases [72]. This has been a recent important achievement, based on a mature understanding of the physical relevance of quasi-local charges in integrable models [73][74][75][76][77][78][79][80][81][82] and mainly motivated by studies regarding the validity of the so-called Generalized Gibbs Ensemble [83][84][85][86][87][88][89][90][91][92][93][94][95][96][97]. However, even in this problem further simplifications occur for integrable states, allowing one to reach closed-form analytical formulas, via the so-called Y -system [72,94]. In any case, our work provides a unified point of view for the many exact results that appeared in the past years. Furthermore, it also provides a useful starting point for the study of models where quench problems have not been investigated: if exact results are to be derived, one should first look at the integrable initial states, for which they are expected.
The organization of this article is as follows. In Sec. 2 we review the classical work of Ghoshal and Zamolodchikov and present the general setting. Integrable states are defined in Sec. 3, where we also discuss some of their properties. In Sec. 4 we show that states related to integrable boundaries in the rotated channel are integrable according to our definition, while in Sec. 5 we present a general construction of integrable matrix product states with arbitrary bond dimension. Sec. 6 is devoted to a critical review of several models considered in the recent literature. Finally, our conclusions are reported in Sec. 7.

Boundary states in integrable quantum field theory
A source of inspiration for the present paper is the classical work of Ghoshal and Zamolodchikov [65], where integrable QFTs in the presence of boundaries are studied. It is natural to start our discussion by briefly reviewing the aspects of that work that are directly relevant for us. Although our focus will be on lattice integrable models, this will allow us to introduce some constructions by analogy with the picture in boundary QFT.
We start with an Euclidean field theory defined on a semi-infinite plane x ∈ (−∞, 0), y ∈ (−∞, +∞). In the so-called Lagrangian approach it can be defined by introducing the corresponding action, which generally reads as Here ϕ(x, y) represent a set of local bulk fields, while ϕ B (y) = ϕ B (x, y)| x=0 are the boundary fields. Analogously, a(ϕ, ∂ µ ϕ) and b (ϕ B , dϕ B /dy) are respectively the bulk and boundary action densities. We recall that from the action (1) one can introduce a Hamiltonian picture in Euclidean time, which is closely related to a 1+1 Lorentzian QFT through the Wick rotation. This procedure of analytic continuation between real and imaginary times is routinely performed in order to apply QFT results to equilibrium statistical mechanics [98]. As pointed out in [65], there are two natural ways to introduce the Hamiltonian picture. This is pictorially illustrated in Fig. 1. First, one can identify the coordinate y in (1) to be the Euclidean time direction. At each time, one can associate a physical Hilbert space to the semi-infinite line x ∈ (−∞, 0). The Hamiltonian is written as whereĥ(x) is the bulk Hamiltonian density, while θ B is an additional boundary term. Alternatively, one can introduce an Hamiltonian picture in the rotated channel, by identifying the time direction with the coordinate x. In this case, the Hilbert space at fixed Figure 1. Pictorial representation of a 2D Euclidean field theory in the presence of boundaries. There are two natural ways to introduce the Hamiltonian picture. In the first one, displayed in the sub-figure (a), the Euclidean time direction is chosen to be parallel to the boundary. The physical Hilbert space is associated with a semi-infinite spatial line at fixed time (dashed black line). In the second way, displayed in the sub-figure (b), the time direction is chosen to be orthogonal to the boundary. The latter plays the role of an initial condition, and is identified with a boundary initial state |B .
time corresponds to the infinite line y ∈ (−∞, ∞) and the Hamiltonian coincides with the bulk one without boundary terms The boundary at x = 0 now plays the role of an initial condition in Euclidean time, and can be identified with an initial boundary state |B , in analogy with the classical constructions in conformal field theory [99]. Consider now the bulk action density obtained from (1) by removing the boundary. If the latter is integrable, there exist an infinite set of integrals of motions where T s+1 , Θ s−1 (T s+1 ,Θ s−1 ) are local field of positive (negative) spin s + 1 and s − 1, satisfying appropriate continuity equations [100]. The spin index s takes value in an infinite subset S of the positive integer numbers and it describes the transformation properties under the Lorentz boosts. Following [65], a boundary field theory is defined to be integrable if infinitely many conservation laws survive the addition of the boundary term. More precisely, in this case the boundary Hamiltonian (2) has an infinite number of integral of motions of the form The index s now takes value in an infinite subset S B ⊂ S: the boundary integral of motions are obtained by selecting a subset of the bulk ones, and adding an appropriate boundary term to them. This is completely analogous to the situation encountered in lattice models with open boundary conditions [101,102]. Remarkably, the existence of an infinite number of conservation laws for (2) can be related to a condition on the boundary state |B in the rotated picture. If the Hamiltonian (2) is integrable then [65] The relevance of (6) for our purposes is that it represents a condition of integrability which can be tested by relying uniquely on the knowledge of the initial boundary state and the bulk conservation laws of the theory. This provides the main source of inspiration for our work. The application of the ideas of [65] to lattice models is not evident. A one dimensional quantum model is always equivalent to a two-dimensional one of classical statistical physics, and similar to the QFT setting it is always possible to define a "rotated channel". However, typically there is no Euclidean invariance on the 2D lattice, and the Hamiltonians (or transfer matrices) of the rotated channel differ from the original physical ones. Moreover, even if the identification of the initial states with the integrable boundary conditions of the rotated channel can be carried out, this construction is non-trivial and in any case model dependent.
On the other hand, a definition analogous to (6) can be straightforwardly introduced in a general way in lattice models, where conservation laws are well known in terms of explicit operators on the physical Hilbert space. In the next subsection we review some basic facts on lattice integrable models which are necessary in order to introduce and discuss our definition of integrable states. The latter will be finally presented in Sec. 3.

Lattice integrable models
We consider a generic one-dimensional model defined on the Hilbert space H = h 1 ⊗. . .⊗h L .
Here h j is a local d-dimensional Hilbert space associated with the site j, while L is the spatial length of the system. The Hamiltonian is indicated as whereĥ j is a local operator. In the following, we will always assume periodic boundary conditions. In an integrable model there exists an infinite number of conserved operators commuting with the Hamiltonian (7). Furthermore, these can be written as a sum over the chain of local densities, and are therefore called local charges. The latter can be obtained by a standard construction within the so-called algebraic Bethe ansatz [2,103], which we now briefly sketch.
One of the main objects of this formalism is the R-matrix R 1,2 (u). The latter is an operator acting on the tensor product of two local spaces h 1 ⊗ h 2 (possibly of different dimension), where u is an arbitrary complex number, the spectral parameter. The R-matrix has to satisfy a set of non-linear relations known as Yang-Baxter equations From the R-matrix, another fundamental object can be constructed, namely the transfer matrix where the trace is over an auxiliary space h 0 and where we introduced the monodromy matrix By means of (8), one can show that transfer matrices with different spectral parameters commute Using (11), an infinite set of commuting operators is readily obtained as Importantly, λ * can be chosen in such a way that the operators Q n are written as a sum over the chain of local operators (or densities). For example, within our conventions, the charges Q n are such that in the Heisenberg chains the corresponding densities span n sites. Then, one can define an integrable Hamiltonian as By construction, H L is of the form (7), and has an infinite number of local charges Q n . For the lattice models considered in this work, local conserved charges can be divided into two subsets: the even ones, Q 2n , and odd ones, Q 2n+1 with n ≥ 1. These sets display different behavior under spatial reflection, namely ‡ where Π is the reflection operator Here we introduced the notation where |k j are the basis vectors of the local space h j , with k = 1, . . . d.
We are interested in global quantum quenches where the Hamiltonian driving the time evolution is integrable. For each model, we will consider several families of initial states and in order to allow for a general discussion, we will represent them as matrix product states (MPSs) [105]. They are a class of states which display a number of important properties: they have exponentially decaying correlations and finite entanglement between two semi-infinite subsystems. Furthermore, it is known that MPSs can approximate ground states of gapped ‡ One of the simplest ways to prove the reflection properties is by using the so-called boost operator B, which connects the charges through the formal relation Q k+1 = [B, Q k ] [104]. The boost operator is manifestly odd under space reflection, and this guarantees the alternating signs appearing in (14). Note that the overall normalization of the transfer matrix affects the parity properties of the charges: a rapidity dependent overall factor introduces constant additive terms. However, in the models considered in this work there is always a natural normalization in which the reflection properties (14) hold.
local Hamiltonians with arbitrary precision [106]. This provides a physical motivation to consider them as initial states.
We recall that a generic (periodic) MPS can be defined as Here d is the dimension of the physical spaces h j , while A where d j are arbitrary positive integer numbers, called bond dimensions. The trace is over the Hilbert space h 0 with dimension d 0 . In a finite chain, every vector of the Hilbert space H can be represented in the form (17) [107]. It is common practice, however, to refer to a state as MPS if the bond dimensions in the corresponding representation (17) do not grow with the system size L.
Finally, it is useful to introduce the following definition: we say that a state |Ψ 0 is p-periodic if p is the smallest positive integer such that where U is the shift operator along the chain In order to ensure a proper thermodynamic limit, we will restrict to initial states that are pperiodic, with p arbitrary but finite (and not increasing with L). The constraints imposed so far (namely, finite bond dimensions and p-periodicity) are extremely loose, and allow one to consider a very large family of initial states. Integrable states will be defined in the following as a small subset of the latter.

Defining integrable states
We can now introduce our definition of integrable states, guided by the analogy with the picture in QFT outlined in Sec. 2.1. As we already stressed, identification of initial states and boundary conditions in an appropriate rotated channel requires preliminary work in lattice models, and the construction can be model dependent. However, it is possible to introduce an immediate and general definition of integrable states in terms of annihilation of a subset of the local conserved charges {Q n } ∞ n=1 , similarly to Eq. (6). We propose the following definition: an initial state is integrable if it is annihilated by all local charges of the model that are odd under space reflection. In the lattice models that we consider these coincide with the set {Q 2n+1 }, with n ≥ 1, cf. Eq. (14). Therefore we require in any finite volume L where the charge Q 2k+1 is well defined. Typically Q n is a sum along the chain of local operators spanning n sites, so that (20) is meaningful if (2k + 1) ≤ L. We stress that, even though the definition (20) is inspired directly by [65], the analogy with QFT is a loose one and its usefulness should therefore be appreciated a posteriori. In particular, Eq. (20) seems to hold for all the initial states for which closed-form analytical results could be obtained, at least in the models considered in this work. Furthermore, we will see in Sec. 4 that the class of states satisfying (20) include all of those which can be related to integrable boundaries in the rotated channel. These facts, together with additional considerations presented in the following, constitute strong evidence that (20) provides a meaningful and useful definition.
Based on the experience with the known cases we expect that in new models and new quench situations analytic solutions can be found in the integrable cases. Therefore, it is important to develop tools to efficiently test the integrability of an initial state. In the following we provide such methods, which are connected directly to the definition (20), and are independent of the knowledge of the overlaps or other composite objects such as the Loschmidt echo.
We consider a model defined by a transfer matrix τ (u), with the so-called regularity condition where U is the translation operator (19). Furthermore, we introduce the constants of proportionality in (12) as where α n are chosen such that the charges Q n+1 are Hermitian. From this definition, one can write down the following formal representation Integrability of a p-periodic MPS |Ψ 0 with finite bond dimension is equivalent to requiring where we used that the charges are Hermitian operators. In order to test (24), we introduce the quantities where p is the periodicity of |Ψ 0 . The motivation to introduce a product of two transfer matrices is to cancel those parts of the transfer matrix which are even and therefore irrelevant to the integrability condition. Our statement is that (24) holds if and only if in a neighborhood of u = 0 In the thermodynamic limit this leads to the condition for some K > 0. Indeed, if the initial state is annihilated by all the odd charges in a given finite volume L, then the action of τ (u)τ (−u) reduces to U 2 + O(u L ) and we obtain (27) and hence (28). The proof of the other direction of the statement is also easy, but for the sake of clarity it is reported in Appendix A. Both functions G(u) andG(u) can be efficiently computed on MPSs of finite bond dimension by standard techniques, as we will also show explicitly in Sec. 3.2. As a consequence, Eqs. (27) and (28) provide an efficient test for integrability of given initial states.
An alternative definition of integrability can be given by requiring where Π is the reflection operator (15). First, note that (29) directly implies two-site shift invariance: where we used τ (0) = U . Next, the annihilation by the odd charges follows from (29) simply by Taylor expanding τ (u) at u = 0. Since (29) implies two-site shift invariance, it is a stronger condition than (20). In fact, based on the analytic properties of the transfer matrix it can be argued that (29) follows from (20) and two-site invariance. However, the question whether the latter property is actually a consequence of the annihilation by all odd charges is an open problem.

Transfer matrix evaluation of the integrability condition
We consider a translational invariant MPS defined as and address the computation of the functions G(u) andG(u) introduced in the previous section. It is convenient to employ a graphical notation close to the one routinely used in the literature of tensor networks. This will help us to reduce the level of technicality of our discussions. First, we represent the R-matrix as so that for the transfer matrix (9) one has the graphical representation The horizontal line denotes the auxiliary space h 0 . Trace over h 0 is denoted with dashes, while vertical lines are are in 1-to-1 correspondence with the sites along the chain.
Using this notation the functions G(u) andG(u) in Eq. (25) and (26) can be represented as partition functions of appropriate 2D statistical physical models. This is pictorially depicted in Fig. 2, where filled boxes represent the matrices appearing in the MPS, whereas the internal lines stand for the insertion of the transfer matrices τ (u) and τ (−u).
The partition function of the model displayed in Fig. 2 can be evaluated by building an alternative transfer matrix (generally called the Quantum Transfer Matrix) which acts in the horizontal direction. The calculations are standard, and analogous to those reported, for example, in Ref. [60]. In particular, it is easy to obtain where N = Ψ 0 |Ψ 0 whilẽ and where A 1 and A 2 are d-dimensional spaces on which the matrices A (j) act. HereĀ (i) indicates the complex conjugated of A (i) . An analogous computation can be carried out for G(u) and also for MPSs that are p-periodic.

The integrability conditions (27) translate into a constraint for the eigenvalues {Λ
Comparing with the special point u = 0 we obtain the norm as

It can be seen that (36) is satisfied when for
Condition (37) implies strict annihilation in finite volume. If we only require asymptotic annihilation in the thermodynamic limit then it is enough for (37) to hold for the eigenvalues with maximal magnitude. However, we suggest to use exact annihilation for the definition of integrability, and we will show that all the previous cases studied in the literature fit this definition. It is important that for a given MPS the matrix (35) has finite dimension (the latter growing with the bond dimension of the MPS). Accordingly its eigenvalues can be investigated either analytically or numerically. Therefore, it is immediate to test integrability of given MPSs. In Appendix B we report two examples where the integrability condition is tested both for an integrable and a non-integrable initial state.
It is a relevant question to find all integrable MPSs for a given model. One strategy would be to write down a set of conditions for the matrices A (j) appearing in (31) such that the eigenvalues of the Quantum Transfer Matrix (35) have the necessary properties. However, this approach appears to be extremely complicated and not viable. In general, we have not succeeded in solving this problem. Nevertheless, in the following section we present a general method to construct integrable MPSs of arbitrary bond dimension. We will show later that the integrable MPSs studied in the recent literature all fit into our framework. The question whether all integrable MPSs can be generated by our method is left for future work.

The pair structure
The integrability condition (20) has immediate consequences for the overlaps between integrable states and the eigenstates of the Hamiltonian (which in the following will be also called Bethe states). For a generic model, it is well known that the latter can be parametrized by sets of quasi-momenta or rapidities {λ j } N j=1 , where N is the number of quasi-particles associated to the eigenstate. Any initial state can be written as where the sum is over all the sets of rapidities, and where c {λ j } j are the corresponding overlaps with the initial state. Given a local charge Q r , its action on Bethe states is where q r (λ) is some known function. As a consequence, an integrable state can have a nonvanishing overlap c {λ j } j only with the Bethe states |{λ j } j such that for all n such that Q 2n+1 exists in the chain of length L. In an interacting model this is a very strong constraint for the rapidities, because they have to satisfy a set of additional quantization conditions known as Bethe equations. Accordingly, only special configurations are consistent with (40).
It follows from the space-reflection properties that the functions q 2n+1 (λ) are odd with respect to λ. Depending on the specific model chosen, there can be a finite set S λ of rapidity values such that q 2n+1 (λ) = 0 if λ ∈ S λ . Accordingly, the constraint (40) is obviously satisfied by states corresponding to a set where R is some positive integer number. Eq. (41) encodes the so-called pair structure. More precisely, we say that an initial state has the pair structure if it has non-vanishing overlap only with Bethe states corresponding to sets of rapidities of the form (41). The presence of the pair structure was already observed in the seminal works [28,29,34,46,47] in the Lieb-Liniger and XXZ models, and has been repeatedly encountered in the recent literature for different initial states and systems. In the context of boundary QFT, it leads to the specific "squeezed" form of the boundary states [65], see also [31,38,42,67,[108][109][110][111][112][113]. Furthermore, the pair structure has immediate consequences for the entropy of the steady state arising after a quench [57,58,114], in particular for its relation with the so-called diagonal entropy [114][115][116][117][118][119][120]. Therefore, it is an important question whether the pair structure follows in general from (40). In passing, we note that in the context of QFT it has been argued that for interaction quenches near criticality the pair structure does not occur. Namely the initial state (the ground state of a critical Hamiltonian) does not consist uniquely of pairs of particles with opposite momenta [18], see also the related work [121].
In a generic interacting model the Bethe equations are algebraically independent of (40) and a fine tuning of the couplings and/or an interplay with the volume parameter is needed to find exceptions to the pair structure. Such fine tuned examples can be found in the XXZ model at the special points ∆ = cos(pπ/q), including the free fermionic point ∆ = 0, which will be treated in more detail in section 6.2. In the case of the isotropic Heisenberg chain a rigorous proof of the pair structure can be given using a result of [122]. Since the argument is very simple, we present it here for completeness.
For an arbitrary Bethe state we have the relation Π|{λ} N ∼ |{−λ} N . Taking scalar product of the two sides of the equation (29) with the Bethe state |{λ} N we see that the overlaps can be non-zero only if for all u, where τ (u, {λ} N ) is the eigenvalue of the transfer matrix corresponding to |{λ} N . It was proven in [122] that the spectrum of the transfer matrix τ (u) is simple: if two different eigenstates share the same eigenvalue then they belong to the same SU (2) multiplet. This implies that |{λ} N , and Π|{λ} N are related by an SU (2) rotation. Since Π does not change the S z eigenvalues, we conclude that Π|{λ} N ∼ |{λ} N , and the pair structure holds.
To conclude this section, we stress that while the pair structure might or might not hold for a specific model (and has to be investigated separately), the conditions (20) and (40) are general unifying properties of integrable initial states.

Relation with integrable boundaries
In the previous sections we have introduced the definition of integrable states in terms of the bulk conserved charges of the Hamiltonian. This definition was inspired by the picture in QFT, where integrable boundaries are directly related to integrable initial states. It is then a natural question to ask, whether such a relation holds also in lattice models, and whether our definition (20) is compatible with it.
In this section we prove one direction of this relation: we present a method to relate integrable boundary conditions to initial states, and show that states obtained in this way indeed satisfy the condition (20). This construction naturally produces local two-site product states, which are presented in Sec 4.1. In Section 5 we show how integrable MPSs can also be taken into account in this framework.

The general construction
The construction to relate integrable initial states to integrable boundaries relies on the path integral evaluation of the partition function In QFT the exchange of time and space directions is straightforward due to the Euclidean invariance of the path integral. However, this is less immediate in lattice models, where space is discrete and time is continuous. The standard method to circumvent this problem is to introduce a discretization in the time direction and then to develop a lattice path integral for the resulting partition function. This is achieved by employing a Trotter decomposition In integrable models the Hamiltonian H can be related to the transfer matrices, and this makes it possible to introduce the lattice path integrals. For the sake of clarity here we focus on the XXZ spin-1/2 model, defined by the Hamiltonian where σ α j are the Pauli matrices, and periodic boundary conditions are assumed, σ α L+1 = σ α 1 . In this case, the R-matrix is where η = arccosh(∆). The normalization of the R-matrix is such that the corresponding transfer matrix both satisfies the regularity condition (21) and has charges Q n with the correct even/odd behavior (14). Moreover, it satisfies the so-called crossing relation Here R t 0 0,1 (u) denotes partial transposition in the space h 0 , V 0 is a gauge matrix acting on h 0 while γ(u) is a function. In particular, for the R-matrix (46) one has γ(u) = sinh(u)/ sinh(u + η) and V 0 = σ y . While here we focus on the XXZ spin-1/2 chain, the construction described in this section can be carried out straightforwardly for all the integrable models whose R-matrix satisfies a crossing relation, and is thus very general. An example is given by higher spin versions of the Hamiltonian (45), as we will see in the following. On the other hand, generalization of this construction to models for which there is no relation of the type (47) (such as the SU (3)-invariant spin chain) is not evident and needs further research.
We follow the derivation of [60], to which we refer for all the necessary technical details, providing here only the main ideas. For simplicity, we start by considering an initial state of the form while MPSs with arbitrary bond dimension will be considered in the next section. It is important to recall that, similarly to (2), one can consider the Hamiltonian (45) where and where ξ j are arbitrary inhomogeneities. The 2×2 matrices K ± 0 (u) are boundary operators acting on the auxiliary space h 0 , which in this case read K ± (u) = K(u ± η/2, ξ ± , κ ± , τ ± ), with Here ξ, κ and τ are arbitrary parameters. The K-matrix has to satisfy the so-called reflection (or boundary Yang-Baxter) equations which guarantee the commutativity of the two-row transfer matrices (49).
We have now all the elements to evaluate the partition function (43) from the Trotter decomposition (44). In the particular case of the XXZ Heisenberg chain, one can write [123,124] where τ (λ) is the usual periodic transfer matrix and w * = (w sinh η)/2. Eqs. (44) and (54) can be interpreted as follows: the continuous time evolution can be approximated by a discrete one where the elementary time step is obtained from the application of a two-row transfer matrix. In order to parallel the discussion of Sec. 2.1, we consider Euclidean time.
From the Trotter-Suzuki decomposition (44) and (54) it is evident that the computation of (43) reduces to a classical partition function in two dimensions. This is illustrated in Fig. 3, where we made use of the graphical notation introduced in Sec. 3.2. From Fig. 3 the analogy with the field theory case displayed in Fig. 1 becomes evident. Indeed, in the two-dimensional lattice the Euclidean time direction can now be chosen parallel to the boundaries, as shown in Fig. 4. Accordingly, the partition function can be thought of as generated by iterative application of an open transfer matrix in the rotated channel, which implements the discrete time evolution.
The question now is: for which initial states the open transfer matrix generating the time evolution in the rotated channel is integrable? (cf. Fig. 4). The answer to this question can be found as follows. By working out the algebraic steps encoded in Figs. 3 and 4 (see Ref. [60]), the computation of (43) reduces to the problem of finding the leading eigenvalue of the operator in u = 0. Here we have defined A pictorial representation of the operator (55) is given as follows : Then, we find that the initial state is related to integrable boundaries if T (u) is proportional to an appropriate integrable open transfer matrix τ B (u). Going further, one can easily infer that this can happen only if the initial state is written in terms of the matrix elements of K(u) in Eq. (52). Indeed, following the steps sketched above (and reported in full detail in Ref. [60]), one can explicitly write down a parametrization of the family of initial states related to integrable boundaries. For example, in the XXZ spin-1/2 chain, this reads where the functions k ij (u) are defined in (52). Explicitly Eq. (58) is the final result of this construction. Namely, we have identified a family of initial states which play the role of the boundary states in QFT, cf. Sec. 2.1. In the following, we will sometimes refer to them as lattice boundary states. Once again, the derivation presented in this section is completely general, provided that the R-matrix satisfies the crossing relation (47) for an appropriate matrix V . Repeating the steps above, one can write down an expression for the boundary states in a generic model, which reads, up to a global numerical factor, where we employed the notation (16) for the basis vectors of the local Hilbert spaces h j . Here [K(−η/2)V ] i,j are the matrix elements of the product K(−η/2)V , while d is the dimension of the spaces h j . It is straightforward to check that (59) reduces to (57) for the spin-1/2 chain (where V = σ y ), while an explicit expression for the spin-1 case will be provided in the following.
In the next section we will prove that boundary states are integrable according to our definition, as they satisfy (20). Incidentally, we note that in the special case of the spin-1/2 chain any two-site state can be parametrized as in (58) (where one also allows for the parameters to go to infinity). We will see that this is not usually the case for arbitrary models, where only a subset of two-site product states are integrable.

Integrability from reflection equations
In the previous section we have seen that a two-site product state (48) is related to integrable boundaries in the rotated channel if its building block |ψ 0 is written in terms of the elements of the reflection matrix K(u), see. Eq. (57). In this section we show that for these states the  condition (20) follows from general properties of integrability. As a consequence, we establish a direct connection between integrable boundaries in the rotated picture and integrability in terms of annihilation by bulk conserved charges. In turn, this unveils a direct connection with the pair structure discussed in Sec. 3.3.
For the sake of clarity we once again detail the case of the XXZ spin-1/2 chain. However, it will be clear from our discussion that our treatment is in fact much more general, and it holds straightforwardly in all the cases where the R-matrix satisfies the crossing relation (47) (see, for instance, Sec. 6.3 where the case of spin-s chains is worked out).
From the above definition it is clear that |Φ 0 (−η/2) = |Ψ 0 . Furthermore, it follows from the reflection equations (53) thať Figure 5. Pictorial representation of the reflection equations (62).The orientation of the arrows reflects the fact that (62) is written in terms of boundary states, rather than boundary Kmatrices.
Eq. (62) is a crucial relation, which is pictorially represented in Fig. 5. Even though we are now focusing on the spin-1/2 chain, it is in fact quite general: it is simply a rewriting of the boundary reflection equations in terms of states, rather than boundary K-matrices. As one should expect, for a generic model the state (61) will be replaced by a different one, related to the corresponding K-matrix. An explicit example will be given in Sec. 6.3 for the case of higher spin chains.
Consider now a chain of length L with two auxiliary spaces h 0 and h L+1 , where L is an even number. From repeated use of Eq. (62), as pictorially represented in Fig. 6, one can prově Setting v = 0 we geť From the definition ofŘ i,j , the l.h.s. of (65) can be rewritten as where we introduced |W i,j = P i,j |φ(−η/2 + u) i,j . Analogously, the r.h.s. of (65) yields We now make use of the identities P L+1,L . . . P 1,2 P 0,1 = U −1 , where U is the shift operator (19). Plugging these into (66) and (67) we finally obtain From this equation it is now easy to conclude the proof, by writing down its components. In particular, it is shown in Appendix C that (70) implies which is exactly (29), a pictorial representation of which is given on figure 7. The proof presented in this section has far reaching consequences. Most prominently, it directly relates boundary states on the lattice with the pair structure frequently encountered in the recent literature of quantum quenches. In particular, this also provides us with a direct relation between the presence of latter and the validity of the so-called Y -system. Since this is a rather technical point, we consign its discussion to Appendix D.
Our proof relied on a direct application of the boundary Yang-Baxter relations, which leads to the annihilation by the odd charges. However, it is possible to formulate an alternative proof which derives the eigenvalue condition (37). The idea is to introduce the Quantum Transfer Matrices for the inhomogeneous states (60), and to use their commutativity and  (29) for a generic integrable initial state |Ψ 0 . It is proven in the main text that for states generated from boundary integrability this relation is an immediate consequence of (64), which is depicted in Fig. 6. certain simple properties at degenerate points. For the sake of brevity we omit the details of this second proof.

Constructing integrable matrix product states
In this section we address a systematic construction of integrable MPSs of arbitrary bond dimension. The main idea is to obtain new integrable MPSs starting from the known boundary two-site product states. Once again, for the sake of clarity we focus on the XXZ spin-1/2 chain, but it will be apparent that our construction is in fact more general.
As a first example consider the state where |Ψ 0 is a boundary two-site product state of the form (48) and τ (u) is the fundamental transfer matrix (9). By the proof presented in the previous section |Ψ 0 is integrable. Then, it follows straightforwardly from the commutativity of the transfer matrices that (72) is also integrable, as it is annihilated by all the odd charges. This simple observation is at the basis of our construction. More generally, further integrable states can be constructed using the so-called fused (higher-spin) transfer matrices {τ (d) (u)} ∞ d=1 , where we used the convention τ (1) (u) = τ (u). These operators have a similar matrix product form and can be written as Here L (1,d) j,0 (u) are the fused Lax operators, which are matrices acting on the tensor product of the local spaces h j C 2 and the auxiliary space h 0 C d+1 . The trace is taken over h where 1 is the identity operator on the space h 1 ⊗h 2 , σ α are the Pauli matrices while S α are the operators corresponding to the standard (d + 1)-dimensional representation of SU (2). In the anisotropic case, ∆ = 1, Eq. (74) has to be replaced by an appropriate deformed expression, involving the generators of the quantum group U q (sl 2 ) [125]. The operators (73) are called fused transfer matrices, as they can be obtained from (9) by an appropriate procedure named fusion [126].
Importantly, all the transfer matrices (73) commute, so that we can immediately construct an infinite family of integrable MPSs. Consider where |Ψ 0 is an arbitrary boundary two-site product state of the form (48). It follows immediately, that (76) is integrable, as where we used that the fused transfer matrices commute with the local charges (which follows from (22) and (75)). The state (76) can be cast into the canonical MPS form where H D is an auxiliary space of dimension D = n j=1 (d j + 1). Here A and their explicit form can be derived from the knowledge of the operators L (1,d j ) (u j ). Note that |χ in (78) is in general two-site invariant, analogously to the state |Ψ 0 .
We note that the MPSs (76) could be interpreted as lattice versions of the "smeared boundary states" introduced in the study of quantum quenches in the context of conformal field theories [17,63,64] Importantly, the MPSs (76) can also be related to integrable boundaries in the rotated channel. Consider for example the state where |Ψ 0 is defined in (48), with its building block |ψ 0 satisfying (57). Then, in analogy with (33) one has the following pictorial representation Here we indicated the auxiliary row (corresponding to a Hilbert space of dimension d + 1) with a thick red line, to distinguish it from the 2-dimensional representation appearing in the usual transfer matrix (9). One can consequently repeat the steps outlined in Sec. 4: in this case the time evolution in the rotated channel is represented pictorially in Fig. 8. We see that application of τ (d) (w) only results in the insertion of a line in the two-dimensional partition function with respect to the situation displayed in Fig. 4. In particular, using the properties of fused transfer matrices one can see that the open transfer matrix appearing in this construction is still integrable. Finally, it is clear that the same happens by applying a From (76), we see that the set of integrable states is infinite, and their construction involves a large number of free parameters. Even fixing the number n of transfer matrices, their bond dimensions d j and spectral parameters can still be chosen arbitrarily. It is an important question whether this family exhausts all possibilities of integrable MPSs with finite bond dimension. While we can not give a definite answer to this question, it is remarkable that all known cases indeed fit into this framework, including the MPSs studied in [35,39]. This will be detailed in the following section.
An alternative way to understand our construction is to interpret the MPSs as non-scalar solutions to the boundary Yang-Baxter equations. Then the task of finding all integrable MPSs can be split into two parts: determining whether the boundary Yang-Baxter equations are necessary for the integrability condition to hold, and finding all finite dimensional solutions in terms of local Lax operators acting on local K-matrices. We leave this problem for future research.
To conclude this section we note that the overlaps between integrable MPSs of the form (76) and the eigenstates of the Hamiltonian are immediately obtained if the overlaps with the state |Ψ 0 are known. This follows from the fact that the transfer matrices act diagonally on the Bethe states and the eigenvalues are known from the algebraic Bethe ansatz, see for example [127] for explicit formulas. Employing the same notation of Sec. 3.3, we have where τ (dr) (w r |{λ j }) is the eigenvalue of τ (dr) (w r ) corresponding to the eigenstate |{λ j } .

Integrable quenches: analysis of specific models
In this last section we review several recent studies of quantum quenches in different models, where closed-form analytical results could be obtained. We show that in all of these cases the initial states are integrable according to our definition. These include the MPSs constructed in the works [35,39], which are shown to fit into the framework of the previous sections. We also present new results by producing concrete formulas for the integrable two-site states of the spin-1 XXZ chain.

The XXZ spin-1/2 chain
We begin our analysis with the prototypical case of the XXZ spin-1/2 chain (45), where quench problems have been extensively investigated in the past few years [46,47,59,60,80,86,89,92,94,114,[128][129][130][131]. In particular, a large number of analytical results have been obtained. Exact overlap formulas between assigned initial states and the eigenstates of (45) have been derived for special classes of two-site states [26][27][28][32][33][34] and, in the isotropic case, for more general matrix product states [35,39,40]. Closed-form results for the long-time limit of local observables were derived in [46,47], by means of the Quench Action method [44,45], while exact computations for the Loschmidt echo have been reported in [59,60] for arbitrary twosite product states. Finally, an exact result for the time evolution of entanglement entropies has been obtained in [132]. From the previous section, it is now clear that all two-site states considered in these works are integrable, as they are boundary states. These include, as a special case, the Néel and the so-called dimer states studied in [46,47], but also the tilted Néel and ferromagnet states considered in [94]. Note that even though these states are very simple, the computation of their overlaps is extremely difficult: in fact, the latter are still unknown in the case of tilted Néel and ferromagnet states. Nevertheless, it follows from our derivation that for generic ∆ the pair structure holds for all local two-site states.
Among the most interesting recent developments was the discovery of exact overlap formulas for MPSs in the isotropic case ∆ = 1 [35,39,40]. The overlaps were shown to have the same structure as in the case of the Néel state: they included the same Gaudin-like determinants, and only the pre-factors were different. Here we show how these MPSs can be embedded into our framework of integrable initial states. In particular, we argue that they can be obtained by the action of (fundamental or fused) transfer matrices on simple two-site states; in the first few examples we explicitly calculate the corresponding dimensions and spectral parameters.
In [39] the following family of states was considered where Here S (k) α are the (k + 1) × (k + 1) matrices corresponding to the standard representation of SU (2) generators, and the trace in (82) is over the associated (k + 1)-dimensional space.
In order to test integrability of these MPSs we numerically constructed the corresponding QTMs (35) in the first few cases. By evaluating the eigenvalue condition (37) we confirmed integrability of these states. In fact, this also directly follows from the corresponding overlap formulas which were computed in [39], explicitly displaying the pair structure.
The question of how these states fit into our framework of boundary integrability is not immediately clear from their MPS representation. Below we show by explicit calculations that the two simplest vectors |χ are the translationally invariant components of simple two-site product states; we believe that this is a new result. Going further, the structure of the next few states could be investigated using a recursion relation derived in [39], which expresses |χ . This relation explicitly involves fundamental transfer matrices of the form (9) with an auxiliary space of dimension 2 and is such that higher MPSs are expressed as sums of lower ones. Our goal is to relate these to integrable boundaries, namely we intend to express higher MPSs in the product form (76). This could be achieved by a careful analysis of the recursion formulas of [39] and by the so-called T -system, a set of functional relations for different fused transfer matrices [127]. However, the full implementation of this program is beyond the scope of this paper, and we content ourselves with deriving explicit formulas for |χ  (2), so we are allowed to perform global rotations of the physical space (if necessary, this rotation can easily be restored at the end of our calculations). We consider the operator where W j is a local rotation at the site j, described by the matrix After applying the operator (85), the state (82) is rewritten, up to an irrelevant global numerical phase, as where nowZ For k = 1, it is immediate to simplify this expression. In this case (S ± ) 2 = 0, and we get immediately where U is the shift operator (19), while |N is the Néel state, We see that, for k = 1, (87) is nothing but the zero-momentum component of the Néel state. Restoring the rotation by (85) we see that in this simplest case the MPS (82) is the zeromomentum component of the tilted Néel state in the y-direction: A decomposition similar to (90) can be performed for generic k, using standard techniques in the literature of MPSs. In particular, one can prove that (87) where Here we have defined where we introduced the projectors P 1 and P 2 , which are diagonal matrices with elements One can now check for each value of k that |Φ (k) is of the form (76) for an appropriate choice of the two-site state |ψ 0 . We have done this explicitly up to k = 5. In particular, we obtained where |N is the Néel state (91) while U is the shift operator. After rotating back with the inverse of (85), we obtain the following list for the original MPSs: This representation is a new result of our work. We conjecture that all |χ can be written in a similar product form using higher spin fused transfer matrices. The proof of this conjecture and the derivation of explicit formulas is beyond the scope of this work.
As a final remark on the spin-1/2 XXZ model, we note that additional integrable MPSs can be constructed for special values of the anisotropy ∆, in the regime ∆ < 1. These correspond to the so-called "root of unity points", where additional representations of the underlying quantum group, and hence of transfer matrices, exist [125]. MPSs obtained using these additional transfer matrices can be incorporated into our discussion, because they still commute with the transfer matrix (9).

The XX model
At the non-interacting point ∆ = 0 the XXZ Hamiltonian (45) deserves special attention. In this case, one recovers the so-called XX chain, which has served as a prototypical benchmark for countless studies on non-equilibrium dynamics in isolated many-body systems. In particular, analytic results for global quantum quenches have been presented in [37,128,129], where the dynamics arising from the Néel state was mainly addressed. This model presents an important example where integrable states do not display the pair structure discussed in Sec 3.3. In particular, the special structure of conserved charges lead to a different constraint for the rapidities of the eigenstates annihilated by the odd ones. This is briefly discussed in the following, while we refer to [37] for a systematic treatment of quantum quenches from the Néel state.
We recall that the XX model can be studied by introducing fermionic operators through the Jordan-Wigner transformation which satisfy the relations: The Hamiltonian is then written as H = p ε pc † pc p , wherẽ e −ipj c j , and ε p = cos(p). Alternatively, the model can be studied by considering the limit ∆ → 0 (namely, η → iπ/2) of the Bethe ansatz solution of (45) for generic ∆. In the latter language the quasi-particles are associated to the rapidities λ, which are related to the lattice momentum through It follows from the Bethe equations that the Bethe rapidities are either real or they have imaginary part equal to π/2, and the corresponding one-particle energy eigenvalues are e(λ) = − 1 cosh(2λ) .
As in the case of generic ∆, a set of conserved charges Q j can be generated via the transfer matrix construction explained in Sec. 2.2. Alternatively, a set of charges can be constructed using the fermionic operators, such that their commutativity follows from (108). The relation between the two sets of charges is non-trivial and has been studied in detail in [104], see also [78,133]. From the fermionic point of view, one possible choice for the charges is simply the set of occupation numbersñ k for the Fourier modes. These operators are inherently non-local, but they form a complete commuting family. An alternative choice is to consider the local operators

and the Hermitian combinations
It follows directly from the commutation relations (108) that every I j is conserved. The charges I + j are even and the I − j are odd under space reflection. On the other hand, it was shown in [104] that in the XX model all charges can be expressed using the operators e αβ n = L j=1 e αβ n,j , e αβ n,j = σ α j σ z j+1 . . . σ z j+n−2 σ β j+n−1 .
In particular, there are two families of charges defined as H + n = e xx n + e yy n , n even , e xy n − e yx n , n odd , and H − n = e xy n − e yx n , n even , e xx n + e yy n , n odd .
The canonical charges Q j obtained from the transfer matrix are linear combinations of the first family. Explicit formulas can be found in [104]: in the simplest example Q 3 = −H + 3 , for which the one-particle eigenvalues read Using the Jordan-Wigner transformation it can be seen that the families {H ± n } and {I ± n } contain the same operators, but with alternating identification: In the present case, integrability of the initial state is equivalent to annihilation by the charges H + n for all odd n. In the following, we show that this does not imply the pair structure. Consider the rapidity transformation It is straightforward to verify that this implies Analogously, one can see that the eigenvalues of all odd charges are invariant under (117).
Taking an arbitrary eigenstate, this transformation can be performed for every rapidity individually, so that a large number of states can be generated which share the same eigenvalues for all odd Q j . Eigenstates corresponding to sets of rapidities (41) are obviously annihilated by the odd charges, but so do all the other states obtained after repeated application of (117). The rapidities of these eigenstates do not uniquely display pairs of opposite rapidities. Since integrable states will in general overlap with all eigenstates annihilated by odd conserved charges, it follows from this discussion that they will not posses the pair structure. This is explicitly shown in the case of the Néel state in [37], to which we refer for further detail. As a final remark, we note that the XX chain is closely related to the q → ∞ limit of the q-boson model studied in [134,135]. However, the two models are connected by a highly non-local transformation, which alters the locality of the charges. A detailed analysis of the latter for the q-boson model is out of the scope of the present paper.

XXZ chains with higher spin
Analytic results for quantum quenches in higher-spin chains were presented in [94]. In particular, closed-form characterizations of post-quench steady states were obtained for the spin-1 chain known as the Zamolodchikov-Fateev model [136]. The Hamiltonian reads where the indices a, b in the second sum take the values x, y, z and where periodic boundary conditions are assumed, S α N +1 = S α 1 . The coefficients A ab are defined by A ab (η) = A ba (η) and while η plays the role of the anisotropy parameter along the z-direction. Finally the operators S α j are given by the standard three-dimensional representation of the SU (2) generators Two particular initial states were considered in [94], namely where |j , j = 1, 2, 3 represent the basis vectors of the local Hilbert space h C 3 . For these states analytical formulas were obtained by means of a Y -system, cf. Appendix D.
Not surprisingly, one can show that these states are integrable according to our definition. In particular, they belong to the class of boundary states, as we detail in the following.
It is easy to check that the states (124) and (125) belong to this class. Then, from the proof of Sec. 4, which also applies for the spin-1 case, we immediately obtain that they satisfy the condition (20), and hence they are integrable. Note that the integrability of the class (139) can also be verified by constructing the Quantum Transfer Matrix (35), and numerically evaluating the corresponding eigenvalues in a neighborhood of u = 0. Contrary to the spin-1/2 formula (61), it is not true that all two-site states can be parametrized as in (140). In other words, boundary states for the spin-1 model are only a subset of the two-site states: one can explicitly check that for arbitrary choices of the latter (28) does not hold. This is actually true for all the spin-s generalizations of the XXZ Hamiltonian (45) with s ≥ 1. These models can be obtained by the well-known fusion procedure [126], from which also the corresponding K-matrices can be built [138]. As the spin s increases, the dimension of the local Hilbert spaces h j will also increase. On the other hand, there are always just 3 free parameters of the fused K-matrices, and they are not enough to parametrize all the states in the tensor product h j ⊗ h j+1 .

The SU (3)-invariant chain
Finally, we touch upon higher rank generalization of the SU (2) chain, and focus on the SU (3)-invariant Lai-Sutherland Hamiltonian [139] Here the local Hilbert space is h j C 3 , while S α j are given once again by the standard three-dimensional representation of the SU (2) generators (123). The analytical description of this model is significantly more involved as it requires a nested Bethe ansatz treatment [5]. Nevertheless, many elements of the algebraic construction discussed in section 2.2 are also valid in this case. The R-matrix of the model is where P 1,2 is the permutation matrix (63). The transfer matrix can be simply obtained by (9). Recently, a remarkable overlap formula was conjectured in [41] for a particular matrix product state, with a form which is reminiscent of the one in SU (2) chains. The state is where the trace is over the auxiliary space h 0 C 2 , and where σ α 0 are Pauli matrices acting on h 0 . Note that the auxiliary space has a different dimension from the physical spaces h j C 3 . This state is translational invariant, and the integrability conditions can easily be tested. We constructed the corresponding QTM and verified the eigenvalue condition (37) numerically, demonstrating that the state (143) is indeed integrable. Note that this also follows from the conjectured form of the overlap with the Bethe states [41], from which the pair structure is evident.
An important question is whether this state can be understood in terms of integrable boundary conditions in the rotated channel. On the one hand, in the SU (3) case the study of integrable boundary transfer matrices is more complicated [140]. On the other hand, in the argument of Sec. 4.2 we explicitly used the crossing relation (47) of the R-matrix, which is no longer true for the SU (3) model [141]. Accordingly, care has to be taken to generalize those constructions to this case. We hope to return to these topics in a future work.

Conclusions
In this work we have proposed a definition of integrable states for quantum quenches in lattice integrable systems, which is directly inspired by the classical work on boundary quantum field theory of Ghoshal and Zamolodchikov [65]. We have identified integrable states as those which are annihilated by all the odd local conserved charges of the Hamiltonian. We have proven that these include the states which can be related to integrable boundaries in an appropriate rotated channel, in loose analogy with QFT. Furthermore, we have shown that in all of the known cases where closed-form analytical results could be obtained, the initial state is integrable according to our definition.
In the prototypical case of XXZ spin-s chains we have shown that integrable states include two-site product states together with larger families of MPSs. In the spin-1/2 chain all two-site states are integrable, whereas for higher spin this is true only for a subset of them. We have characterized this subset in terms of the fused K-matrices.
One of the properties of integrable quenches seemed to be the pair structure for the overlaps, because this was observed in simple cases in the spin-1/2 XXZ chain and also the Lieb-Liniger model [28,29,34,46,47]. The pair structure has important consequences for the entropy of the steady state arising after a quantum quench [57,58,114], therefore it is important to clarify its relation with the integrability of the initial state. We have argued that the pair structure indeed follows from integrability for generic values of the coupling constants (and could actually be proven for the isotropic Heisenberg chain). Together with our proof of integrability this constitutes a general confirmation of the pair structure for a wide variety of states, including already known cases and new states where the actual overlaps are not yet known. On the other hand, we have also discussed possible exceptions to the pair structure, such as the XX model detailed in Sec. 6.2. Nevertheless, the integrability condition (20) can be always introduced without modifications.
Remarkably, in almost all of the cases encountered, MPSs annihilated by the odd charges could be understood as boundary states in the rotated channel. It is an important open question whether this family exhausts all integrable MPSs with finite bond dimensions. In the SU (3)invariant model the MPS introduced in [41] and studied here in 6.4 is integrable according to our definition, but its interpretation in terms of boundary integrability is not known yet. We hope to return to this question in future work.
As a final remark, we stress that our results should not be interpreted as a no-go theorem for obtaining exact results for non-integrable initial states. However, our work provides a unified point of view for the many exact results that appeared in the past years. Furthermore, if exact results for time evolution are to be derived (be it results for correlators or the Loschmidt echo or other quantities), one should first look at the integrable initial states, regardless of the model considered. The test of integrability is straightforward, therefore our framework gives an extremely useful starting point for the study of models where quantum quenches have not yet been investigated. An example is the case of SU (N )-invariant spin chains, where only a small number of results are currently available [41,55].

Acknowledgments
We are very grateful to Pasquale Calabrese and Gábor Takács for inspiring discussions and useful comments. B. P. is grateful to Pasquale Calabrese and SISSA for their hospitality. E. V. acknowledges support by the ERC under Starting Grant 279391 EDEQS. B. P. acknowledges support from the "Premium" Postdoctoral Program of the Hungarian Academy of Sciences, and the K2016 grant no. 119204 of the research agency NKFIH.

Appendix A. Integrability condition for p-periodic states
In this appendix we show that the validity of (27) implies the annihilation of all odd conserved charges, namely (20). For simplicity, we assume that the initial state is two-site shift invariant with Ψ 0 |Ψ 0 = 1, as an analogous derivation holds in the general case.
First, it follows from (23) that one can write down the following formal expansion Note that here we used U 2 |Ψ 0 = |Ψ 0 and that for parity invariant states G(u) =G(u). At the first orders in u 2 we have If G(u) ≡ 1 then we have immediately that Ψ 0 |Q 3 |Ψ 0 = 0. However, at the next order we have a sum of two terms. If the state is parity invariant, then Ψ 0 |Q 5 |Ψ 0 = 0 and we also have If the state is not parity invariant, then we use the information coming fromG(u). Using that the latter is also identically 1, we get the vanishing of the two terms separately. If (A.3) holds then Q 3 |Ψ 0 = 0, namely this means that Q 3 can be neglected from the Taylor series completely. Going to the next order, we get the vanishing of the mean value of Q 7 , and similarly we can prove Q 5 |Ψ 0 = 0. Proceeding iteratively, the annihilation by all odd charges is proven.

Appendix B. Numerical study of the Quantum Transfer Matrix for integrable and non-integrable MPSs
In this appendix we numerically study the matrixT (u) introduced in (35) and compute its eigenvalues for both an integrable and a non-integrable case.  Figure B1. Magnitude of the eigenvalues of the operatorT (u) corresponding to the (translationally invariant) dimer state (B.1). We see that the eigenvalues which are non-zero for u = 0 remain constant for a non-vanishing neighborhood of u = 0. In fact, they are completely u-independent, but level crossings can occur. The eigenvalues with the second largest magnitude describe the exponentially decaying overlap between the original and the one-site shifted dimer states. For the integrable case, we consider the zero-momentum dimer state defined as It is easy to check that this state is generated by a translationally invariant MPS of the form (31), with As a non-integrable state, we choose the translationally invariant component of a four-site domain-wall state, namely This state was already considered in [94] (and mentioned in [72]), where it was shown that the corresponding rapidity distribution functions did not satisfy the Y -system, cf. Appendix D. Therefore it is natural to expect that this state is not integrable according to our definition. The state in (B.3) can still be written as a translational invariant MPS of the form (31), where now We computed the eigenvalues of the QTMT (u) for these two MPSs. For concreteness, we focused on the XXX chain, with the R-matrix normalization given by R(u) = (u + iP )/(u + i) (here P is the permutation operator defined in (63)). The QTM is symmetric with respect to the sign of u and in the following we restricted to positive values of u. The numerical results for the magnitude of the eigenvalues are shown in Figs. B1 and B2. Note that a single curve here corresponds to at least two different eigenvalues due to the sign difference. Further degeneracies can be present, but we can omit to specify them as they are irrelevant to the test of integrability.
In the case of the dimer state, we see from

Appendix C. Integrability of boundary states: technical details
In this appendix we provide further technical details on the proof presented in Sec. 4.2. We start by showing that (62) is a simple rewriting of the reflection equations (53).

Appendix D. The Y -system
In this appendix we briefly touch upon another property of the boundary states, which is related to the so-called Y -system, an ubiquitous structure of integrability [142]. The latter is a system of equations for a set of functions in the complex plane. In the case of the XXZ model, the Y -system takes the form where for generic values of ∆, one has j = 1, 2, . . . ∞. In the framework of quantum quenches, the Y -system first appeared in the study of quenches from the Néel state [46]. In this case the Y -functions are obtained starting from the Bethe ansatz rapidity distributions of the corresponding long-time steady state. Subsequently, the same relations were found to hold more generally for all two-site product states in the XXZ spin-1/2 model [72,94,96,97], and for some initial states in the spin-1 chain [94]. An explanation for the Y -system in this context was found in Ref. [60]. Using the identification of two-site states with boundary states reviewed in Sec. 4, the Y -system emerged from the fusion properties of the corresponding boundary transfer matrices. From the practical point of view, the existence of a Y -system represents a major computational advance, allowing for a closed-form analytical characterization of the rapidity distribution functions of the postquench steady state [60,72,94].
So far, it was not clear why the existence of the Y -system should be expected to imply the presence of the pair structure discussed in Sec. 3.3. The results of Sec. 4 gives us further understanding of this point. Explicitly, we have proven that boundary states, which were shown to be characterized by a Y -system in [60], satisfy the condition (20). In turn, the latter implies the pair structure if no fine tuning of the couplings is made. Hence, in this case the Y -system and the pair structure are seen to have the same origin, rooted once again in integrability.
We remark that if the initial state is such that the overlaps can be factorized algebraically, i.e. they can be written in the form of such that v(λ) is a meromorphic function and C({λ j } N/2 ) is O(L 0 ) in the thermodynamic limit, then the resulting Y -functions (as obtained within the Quench Action formalism [44]) necessarily satisfy the Y -system equations. This can be proven using a simple analytic manipulation of the resulting TBA equations, in analogy with the methods presented in an early work on Y -systems by Al. B. Zamolodchikov [143]. The same statement can be proven even for non-integrable quenches, where the pair structure does not hold: the Y -system would still hold if the overlaps can be factorized as It follows that if the Y -system does not hold, then algebraic factorization of the overlaps is not possible. This argument suggests that the specific form (D.2) is a further characteristic property of integrable initial states.

Appendix E. Generalities on matrix product states
In this appendix we provide some technical details on MPSs which are useful for our discussion in Sec. 6.1. In particular, we show that the MPS (87) can always be decomposed as in Eq. (93). We start by noting a special property of the state (87). Consider the operator is the complex conjugate ofZ (j) , and whereZ (1) andZ (2) are given in (88) and (89). The eigenvalues of N form pairs with the same magnitude and different sign. This can be seen by performing a similarity transformation using the matrix C (k) , which represents a rotation of π/2 around the z-axis. We have and C (k) ⊗ C (k) N (C (k) ) −1 ⊗ (C (k) ) −1 = −N .
It can be checked that the leading eigenvalues have no further degeneracies. From general theorems regarding MPSs [105], this implies that there exist two projectors P 1,2 acting in auxiliary space such that P 1 + P 2 = 1 and P 1Z (j) =Z (j) P 2 , P 2Z (j) =Z (j) P 1 , j = 1, 2 . (E.2) In our case, the operators P 1 and P 2 are defined by the matrices in (99). Following [105], one can simply compute Using now P 1,2 = P 2 1,2 and (E.2), it is straightforward to recover (93).