Exact solution for the quench dynamics of a nested integrable system

TL;DR

This paper develops an exact method to analyze the non-equilibrium quench dynamics of nested integrable systems, specifically the spin-1 Lai-Sutherland model, revealing detailed steady-state and entanglement propagation characteristics.

Contribution

It extends the Quench Action approach to nested integrable models, providing explicit solutions for the post-quench steady state and entanglement dynamics in the spin-1 Lai-Sutherland model.

Findings

Characterized the post-quench steady state for the model.

Predicted entanglement and mutual information propagation patterns.

Abstract

Integrable models provide an exact description for a wide variety of physical phenomena. For example nested integrable systems contain different species of interacting particles with a rich phenomenology in their collective behavior, which is the origin of the unconventional phenomenon of spin-charge separation. So far, however, most of the theoretical work in the study of non-equilibrium dynamics of integrable systems has been focusing on models with an elementary (i.e., not nested) Bethe ansatz. In this work we explicitly investigate quantum quenches in nested integrable systems, by generalizing the application of the Quench Action approach. Specifically, we consider the spin-$1$ Lai-Sutherland model, described, in the thermodynamic limit, by the theory of two different species of Bethe-ansatz particles, each one forming an infinite number of bound states. We focus on the situation where the quench dynamics starts from a simple matrix product state for which the overlaps with the eigenstates of the Hamiltonian are known. We fully characterize the post-quench steady state and perform several consistency checks for the validity of our results. Finally, we provide predictions for the propagation of entanglement and mutual information after the quench, which can be used as signature of the quasi-particle content of the model.

Figures

Equations

$$H_L=\sum _{j=1}^L\left[{\mathbf{s}}_j\cdot {\mathbf{s}}_{j+1}+{\left({\mathbf{s}}_j\cdot {\mathbf{s}}_{j+1}\right)}^2\right]-2L,$$

$$s^{x}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c c c}0&1&0\\ 1&0&1\\ 0&1&0\end{array}\right),\quad s^{y}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c c c}0&-i&0\\ i&0&-i\\ 0&i&0\end{array}\right),\quad s^{z}=\left(\begin{array}[]{c c c}1&0&0\\ 0&0&0\\ 0&0&-1\end{array}\right)\,.$$

$$\ket{\Uparrow}=\left(\begin{array}[]{c}1\\ 0\\ 0\end{array}\right),\quad\ket{0}=\left(\begin{array}[]{c}0\\ 1\\ 0\end{array}\right),\quad\ket{\Downarrow}=\left(\begin{array}[]{c}0\\ 0\\ 1\end{array}\right)\,,$$

$$|e_1⟩=|\Uparrow ⟩,|e_2⟩=|0⟩,|e_3⟩=|\Downarrow ⟩,$$

$$|e_{{\alpha }_1}{e}_{{\alpha }_2}\dots e_{{\alpha }_L}⟩=|e_{{\alpha }_1}⟩\otimes |e_{{\alpha }_2}⟩\otimes \dots \otimes |e_{{\alpha }_L}⟩.$$

$$H_L^{\mathrm{B}}=\sum _{j=1}^L\left[{\mathbf{s}}_j\cdot {\mathbf{s}}_{j+1}-{\left({\mathbf{s}}_j\cdot {\mathbf{s}}_{j+1}\right)}^2\right],$$

$${\mathcal{N}}_1$$

$${\mathcal{N}}_2$$

$$E_j^i\equiv |e_j⟩⟨e_i|.$$

$$|{\mathbf{k}}_N,{\mathbf{\lambda }}_N⟩$$

$$|\Omega ⟩=|\Uparrow \Uparrow \dots \Uparrow ⟩,$$

$$⟨\mathbf{m}⟩{\mathbf{k}}_{\mathcal{P}},\mathbf{\lambda }$$

$$F_{\mathbf{k}}(\lambda ,s)$$

$$A(\lambda )$$

$$L\ge N,N\ge 2M,$$

$${\left(\frac{k_j+i/2}{k_j-i/2}\right)}^L$$

$$1$$

$$E=-\sum _{j=1}^N\frac{1}{k_j^2+1/4},P=\left[\sum _{j=1}^{N}i\ln \left[\frac{k_j+i/2}{k_j-i/2}\right]\right]\text{mod}2\pi .$$

$${\mathcal{N}}_1|{\mathbf{k}}_N,{\mathbf{\lambda }}_M⟩$$

$${\mathcal{N}}_2|{\mathbf{k}}_N,{\mathbf{\lambda }}_M⟩$$

$$k_{\alpha }^{n,\ell }$$

$${\lambda }_{\alpha }^{n,\ell }$$

$$z_n^{(1)}(k_{\alpha }^n)=\frac{2\pi }{L}{I}_{\alpha }^n,$$

$$z_n^{(1)}(\lambda )=p_n(\lambda )-\frac{1}{L}\sum _{m=1}^{+\infty }\sum _{\beta =1}^{M_m^{(1)}}{\Xi }_{nm}(\lambda -k_{\beta }^m)+\frac{1}{L}\sum _{m=1}^{+\infty }\sum _{\beta =1}^{M_m^{(2)}}{\Theta }_{nm}(\lambda -{\lambda }_{\beta }^m),$$

$$z_n^{(2)}(\lambda )=\frac{1}{L}\sum _{m=1}^{+\infty }\sum _{\beta =1}^{M_m^{(1)}}{\Theta }_{nm}(\lambda -k_{\beta }^m)-\frac{1}{L}\sum _{m=1}^{+\infty }\sum _{\beta =1}^{M_m^{(2)}}{\Xi }_{nm}(\lambda -{\lambda }_{\beta }^m).$$

$$p_n(\lambda )$$

$${\Xi }_{n,m}(\lambda )$$

$${\Theta }_{n,m}(\lambda )$$

$$|I_{\alpha }^n|$$

$$|J_{\alpha }^n|$$

$$E=\sum _{m=1}^{+\infty }\sum _{\beta =1}^{M_m^{(1)}}{\epsilon }_m(k_{\beta }^m),P=\left[\sum _{m=1}^{+\infty }\sum _{\beta =1}^{M_m^{(1)}}\left(p_m(k_{\beta }^m)-\pi \right)\right]\text{mod}2\pi ,$$

$${\epsilon }_n(k)=-\frac{n}{k^2+n^2/4}.$$

$$N_1=\sum _{n=1}^{+\infty }n{M}_n^{(1)},N_2=\sum _{n=1}^{+\infty }n{M}_n^{(2)}.$$

$${\rho }_{t,n}^{(r)}(k)={\rho }_n^{(r)}(k)+{\rho }_{h,n}^{(r)}(k),r=1,2,n=1,\dots ,+\infty ,$$

$${\vartheta }_n^{(r)}$$

Full text

I Introduction

Integrable models provide a simplified conceptual framework where the overwhelmingly complex collective behavior of many-body quantum systems can be investigated. Although special by definition, they provide an exact characterization for a wide variety of phenomena. During the last decade, the interest in integrable models has been reinvigorated by the experimental possibility of engineering quantum many-body systems with ultracold atoms bdz-08 ; ccgo-11 ; gbl-13 ; PolkonikovRMP11 with the unique feature of experimental control over interaction parameters, dimensionality, and isolation. These experiments opened the way to the realization of several one-dimensional systems which are described, to a very good approximation, by integrable models HLFS07 ; kww-06 ; gklk-12 ; fse-13 ; lgkr-13 ; glms-14 ; langen-15 . At the same time, on the theoretical side, the study of integrable models out of equilibrium PolkonikovRMP11 ; GE15 ; cem-16 has brought, among other results, to a better understanding of relaxation mechanisms cem-16 ; ef-16 ; VR:review , of transport dynamics BD:review ; BCDF:transport ; CADY:hydro ; VM:review and of entanglement propagation cc-05 ; CC:review ; ac-16 ; kauf ; d-17 . Concerning the long time dynamics instead, while generic systems locally relax to thermal states D91 ; S94 ; R08 ; RS12 , integrable systems admit a larger class of possible stationary states, described by the generalized Gibbs ensemble (GGE) VR:review ; RDYO07 ; ef-16 . A non-trivial transition between these two ensembles has been observed in the late-time dynamics of weakly non-integrable models MK:prethermalization ; RoschPRL08 ; KollarPRB11 ; worm13 ; MarcuzziPRL13 ; EsslerPRB14 ; fagotti-14 ; konik14 ; BF15 ; CTGM:pret ; knap15 ; SmacchiaPRB15 ; FC15 ; BEGR:PRL ; MenegozJStatMech15 ; BEGR:long ; NIC14 ; pret ; af-17 ; mmgs-16 ; cgdm-17 ; fnr-17 : this phenomenon has been named prethermalization bbw-04 .

One of the most important conceptual and technical advances for the description of the many-body dynamics of interacting integrable models has been the introduction of the Quench Action approach caux-16 ; ce-13 . Its main idea is that at long times the local properties of the system can be captured, in the thermodynamic limit, by a representative eigenstate of the Hamiltonian. In this formalism, the time evolution of local observables is written as a sum of contributions coming from excitations over the representative eigenstate. While the ideas upon which the method is built are simple, its technical implementation to get quantitative predictions might be very hard depending both on the initial state and the model under investigation.

The introduction of the Quench Action approach comes as the result of a long and intense theoretical research activity in the field of non-equilibrium dynamics. At an early age of the field, analytical calculations were mostly performed in free theories Caza06 ; BaSc08 ; CDEO08 ; scc-09 ; CaEF11 ; CoSC13 ; fe-13b ; KoCC14 ; sc-14 ; se-12 ; grd-10 ; mc-12 , where the technical obstacles are less severe. While some results on interacting models are known by other means cc-06 ; fcc-09 ; cg-11 ; fm-10 ; pozsgay-11 ; ia-12 ; mussardo-13 ; la-14 ; d-14 ; alba-2015 ; cubero-16 , it is only with the Quench Action that they could be systematically obtained in many integrable models BeSE14 ; BePC16 , with a great amount of work focusing on the prototypical examples of the XXZ Heisenberg chain pozsgay-13 ; wdbf-14 ; PMWK14 ; ac-16qa ; ppv-17 and the Lieb-Liniger Bose gas dwbc-14 ; dc-14 ; pce-16 ; NaPC15 ; Bucc16 ; npg-17 . These models represent the simplest interacting integrable systems that can be solved by Bethe ansatz. They display richer physics than free models: it was only thanks to the quench action solutions of the XXZ spin-chain wdbf-14 ; PMWK14 that it has been possible to discover that the GGE built with known (ultra)local charges POZS2-13 ; fe-13 ; fcec-14 is not describing correctly POZS2-14 ; ga-14 the steady state and that new families of quasi-local conservation laws should be included in the GGE iqdb-16 ; IDWC15 ; iqc-16 ; PoVW17 . In turn, finding a complete set of charges became a parallel research field prosen-11 ; ip-12 ; prosen-14 ; zmp-16 ; imp-15 ; ppsa-14 ; fagotti-16 ; emp-15 ; impz-16 ; dlsb-15 ; doyon-15 ; bs-16 ; vecu-16 ; cardy-16 ; pv-16 also because when this is known, the stationary state can be built circumventing the quench action solution, as e.g. done in pvc-16 ; PVCR17 .

As we are now entering in a more mature age of the field, it is natural to wonder whether it is actually possible to extend the analytical techniques developed so far to investigate more complex systems and in particular those which can be solved by nested Bethe ansatz. In contrast to the XXZ chain or the Lieb-Liniger gas, these models contain multiple species of interacting particles. The class of nested systems includes the celebrated Hubbard model efgk-05 and multi-species point-wise interacting Fermi and Bose gases yang-67 . The latter are of particular relevance for cold-atomic experimental realizations pmcl-14 and allow us to study spectacular effects such as, e.g., spin-charge separation.

In this paper we study quantum quenches to a prototypical nested model and explicitly provide non-trivial predictions for entanglement propagation after the quench. We consider the spin-$1$ Lai-Sutherland model lai-74 ; sutherland-75 , whose Bethe ansatz solution displays two distinct species of particles, each one forming an infinite family of bound states. We focus on a simple initial matrix product state which allows us to generalize the application of the Quench Action approach to nested models and deal with the corresponding technical difficulties. The overlaps between this matrix product state and the Bethe eigenstates have been reported in dkm-16 , but in a completely different context. We fully characterize the post-quench steady state and perform several checks for the validity of our results. As a numerically measurable prediction, we provide explicit results for the time evolution of entanglement entropy and for the mutual information exploiting a conjecture for entanglement dynamics ac-16 which so far has been carefully tested only for the XXZ spin chains. The entanglement dynamics turn out to provide a signature of different species of particles and can be used for quench spectroscopy similarly to what proposed in kctc-17 ; rmc-16 .

The organization of this paper is as follows. In Section II we introduce the Lai-Sutherland model and review its nested Bethe ansatz solution. The initial state is presented in Section III, together with a brief introduction to the Quench Action approach, while Section IV is devoted to the analysis and solution of the equations characterizing the post-quench steady state. In Section V we investigate the propagation of entanglement entropy after the quench. Finally, our conclusions are reported in Section VI. Five appendices contain the most technical aspects of our calculations.

II The model and the Bethe ansatz solution

II.1 The Hamiltonian

We consider the spin-$1$ Lai-Sutherland model lai-74 ; sutherland-75 , described by the Hamiltonian

[TABLE]

which acts on the Hilbert space $\mathcal{H}_{L}=h_{1}\otimes\ldots\otimes h_{L}$. Here $h_{j}\simeq\mathbb{C}^{3}$ is the local (physical) Hilbert space associated with site $j$. The spin-$1$ operators $s^{a}_{j}$ are given by the standard three-dimensional representation of the $SU(2)$ generators, explicitly

[TABLE]

In the following, we also define the local spin-$1$ basis as

[TABLE]

and will use the labeling

[TABLE]

and

[TABLE]

The Hamiltonian (1) is expressed in terms of the $SU(2)$ spin-$1$ operators (2) and is manifestly invariant under action of $SU(2)$. In fact, it is invariant under the action of the larger group $SU(3)$ as it is apparent by considering its algebraic Bethe ansatz construction kr-81 , as briefly discussed in Appendix A.

The Lai-Sutherland model should not be confused with the different spin-$1$ integrable chain

[TABLE]

which is the Hamiltonian of the $SU(2)$-invariant Babujian-Takhtajan model takhtajan-82 ; babujian-82 . This model can be analyzed by means of the so called fusion procedure krs-81 , starting from the Bethe ansatz solution of the spin-$1/2$ Heisenberg chain; quantum quenches in the model defined by (6) have been considered in pvc-16 . To clarify the difference between the two Hamiltonians, it is useful to introduce the operators

[TABLE]

where we defined

[TABLE]

When applied to a state, the operator $\mathcal{N}_{1}$ counts the number of spins which are either $\ket{0}$ or $\ket{\Downarrow}$ while $\mathcal{N}_{2}$ counts the number of spins $\ket{\Downarrow}$. It is straightforward to see that these operators are mutually commuting and moreover commute with the Hamiltonian (1), which follows directly from the $SU(3)$ invariance of $H_{L}$. On the other hand, these operators do not separately commute with the Hamiltonian $H_{L}^{\rm B}$ (cf. (6)): only their sum does.

This seemingly innocent difference has drastic consequences on the physics of the two models. While the quasi-particle content and the structure of elementary excitations of the theory defined by (6) are analogous to that of the spin-$1/2$ case, the one of the theory described by (1) is completely different: two different species of elementary excitations emerge.

II.2 The nested Bethe anstatz solution

The Hamiltonian (1) is diagonalized by nested Bethe ansatz. Here we briefly sketch the main aspects relevant to our work while we refer to the specialized literature for a systematic treatment efgk-05 ; kr-81 ; johannesson-86 ; johannesson2-86 .

The starting point is to construct a Bethe-state on the chain of length $L$. In the model at hand, this state is parametrized by two sets of complex parameters called rapidities, ${\boldsymbol{k}}_{N}=\{k_{j}\}_{j=1}^{N}$ and ${\boldsymbol{\lambda}}_{M}=\{\lambda_{j}\}_{j=1}^{M}$, as follows:

[TABLE]

Here we defined the reference state

[TABLE]

together with the functions

[TABLE]

The numbers of the two kinds of rapidities, $N$ and $M$, must satisfy

[TABLE]

where $L$ is the size of the chain.

In complete analogy with the well-known Bethe ansatz solution of the spin-$1/2$ Heisenberg chain kbi-93 , one can show that the state (10) is an eigenstate of the Hamiltonian (1) provided that the rapidities satisfy a set of non-linear quantization conditions, known as nested Bethe equations, which read

[TABLE]

The energy and the momentum of the eigenstate $\ket{{\boldsymbol{k}}_{N},{\boldsymbol{\lambda}}_{M}}$ are given by

[TABLE]

Note that since the Hamiltonian (1) is integrable there exist infinitely many local conserved operators (or charges) beyond energy and momentum, which can be constructed by standard techniques kbi-93 . In analogy with (18), the expectation value of higher conserved charges on a Bethe state are immediately obtained once the sets of rapidities characterising the state are known. This is discussed in more detail in Appendix A.

The Bethe states $|\boldsymbol{k}_{N},\boldsymbol{\lambda}_{M}\rangle$ are common eigenstates of the Hamiltonian and of the operators $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ introduced in (7), (8). In particular, one has

[TABLE]

The physical interpretation for the state (10) is straightforward in the case where all the rapidities $\boldsymbol{k}_{N},\boldsymbol{\lambda}_{M}$ are real. In this case $\boldsymbol{k}_{N}$ and $\boldsymbol{\lambda}_{M}$ can be thought as the rapidities of two different species of quasi-particles created on a vacuum represented by the reference state (11); we will call these two species of quasi-particles “bare quasi-particles”. The state (10) is then nothing but a scattering state of bare quasi-particles ZZ79 . Bare quasi-particles of the first species contribute to the energy and momentum of the state, while those of the second species do not (cf. Eqs. (18)). Note that the two species of bare quasi-particles do not directly correspond to the two spin-flips $\ket{0}$ and $\ket{\Downarrow}$. Pictorially one could imagine that $\ket{0}$ is a bare quasi-particle of the first species and $\ket{\Downarrow}$ splits into a bare quasi-particle of the first and one of the second species.

II.3 String hypothesis and thermodynamic description

The thermodynamics of integrable models is naturally described within the well-known thermodynamic Bethe ansatz formalism takahashi-99 . Within this framework, the quasi-particle content of the model emerges in analogy with the case of non-interacting spin chains.

The thermodynamic Bethe ansatz for the nested system of interest in this work has been widely studied in the literature, see, e.g., Refs. johannesson-86 ; johannesson2-86 ; afl-83 ; jls-89 ; mntt-93 ; dn-98 . Here we review the aspects which are relevant for our work. As usual, the starting point is provided by the so called string hypothesis, according to which, for large $L$, both sets of rapidities arrange themselves in the complex plane forming patterns called strings. In the present case the parametrization of the strings reads

[TABLE]

Here the numbers $n=1,2,\ldots,+\infty$ are labeling the string types, the real numbers $k^{n}_{\alpha}$, $\lambda^{n}_{\alpha}$ are the string centers and $\{M^{(1)}_{n},M^{(2)}_{n}\}$ are respectively the number of strings of the first and of the second species. Note that the strings here do not couple the two different species of rapidities and for each species the strings have the the well-known structure encountered in the case of $XXX$ spin-$1/2$ chain takahashi-99 . Physically, the different strings are interpreted as bound states of the bare quasi-particles (see e.g. Chapter IV in Ref. efgk-05 ).

Under the string hypothesis, the Bethe equations (17) can be turned into equations for the string centers; it is convenient to consider the logarithmic form of these equations which reads as johannesson2-86

[TABLE]

Here $\{I^{n}_{\alpha}\},\{J^{n}_{\alpha}\}$ are integers or semi-integers depending on $\{M^{(1)}_{n},M^{(2)}_{n}\}$ and we introduced the counting functions

[TABLE]

Here

[TABLE]

The counting functions (24) and (25) are assumed to be monotonic for any solution of the equations (23). The range of $\{I_{\alpha}^{n},J_{\beta}^{m}\}$ can be obtained taking the $\lambda\rightarrow\infty$ limit of the counting functions (and imposing $\{I_{\alpha}^{n},J_{\beta}^{m}\}$ to be integers or semi-integers according to the values of $\{M^{(1)}_{n},M^{(2)}_{n}\}$). The result reads as dn-98

[TABLE]

where $t_{nm}\equiv 2\min(n,m)-\delta_{n,m}$. The energy and the momentum of a state described by the string centers $\{k_{\alpha}^{n}\},\{\lambda^{n}_{\alpha}\}$ are given by

[TABLE]

where we introduced the bare energies johannesson-86

[TABLE]

Finally, the eigenvalues of $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ read as

[TABLE]

In the thermodynamic limit the string centers become continuous variables on the real line. Accordingly, we introduce the rapidity distribution functions $\rho^{(1)}_{n}(k)$ and $\rho^{(2)}_{n}(k)$ which generalize the momentum distribution functions for free systems. Analogously, we also introduce the hole distribution functions $\rho^{(1)}_{h,n}(k)$ and $\rho^{(2)}_{h,n}(k)$, which corresponds to the distribution of holes (i.e. values of the rapidity for which there is no particle) in free Fermi gases at finite temperature. As customary, we define

[TABLE]

as well as the functions

[TABLE]

Particle and hole distribution functions of the two species can not be chosen arbitrarily, but are constrained by a set of linear integral equations which are nothing but the thermodynamic version of the Bethe equations (17). They are usually called Bethe-Takahashi equations and are derived in complete analogy with the XXX spin-$1/2$ case takahashi-99 . They read mntt-93

[TABLE]

Here we defined the convolution between two functions as

[TABLE]

and

[TABLE]

where

[TABLE]

Following takahashi-99 , the Bethe-Takahashi equations (37), (38) can also be cast in partially decoupled form which is more convenient for numerical analysis. We refer to Appendix B for the derivation, while here we only report the final result

[TABLE]

which uses the conventions

[TABLE]

and the function

[TABLE]

The set of rapidity distribution functions completely characterizes the thermodynamic properties of a given macrostate. In particular, the density of the quasi-particles of the species $(1)$ and $(2)$ can be computed as

[TABLE]

For later convenience, it is also useful to introduce the density of particles forming bound states as

[TABLE]

Analogously, by means of the string hypothesis, one can obtain the density of energy from (31) as

[TABLE]

III The quench protocol

We are now interested in the standard quench dynamics where the system is prepared in an initial state $\ket{\Psi_{0}}$, which is not an eigenstate of the Hamiltonian, and then it is let evolve for $t>0$ with the Lai-Sutherland Hamiltonian (1). The main goal of our analysis is not to reconstruct the entire time evolution, which is currently out of reach, but just to have an exact characterization of the stationary state.

As we have mentioned in the introduction, to tackle this problem, we have two possible strategies. Either we rely on the knowledge of the expectation value of a complete set of quasi-local charges or we construct the stationary state from the Quench Action. Concerning the former approach, while we expect that quasi-local charges might be successfully employed also in the study of quenches to the $SU(3)$ chain (1), a systematic analysis of this problem (i.e. a systematic characterization and classification of the quasi-local charges) in the case of nested systems has not yet been carried out. As a consequence, we can only rely on the Quench Action approach ce-13 , which can be implemented quite generally for any Bethe ansatz integrable model, but is limited to those quenches for which the overlap between the initial state and the Bethe eigenstates are known exactly. The computation of the overlaps has in fact turned out to be a very hard technical problem, and no general scheme to solve it has yet been developed kp-12 ; pozsgay-14 ; bdwc-14 ; pc-14 ; fz-16 ; msca-16 ; hst-17 ; cd-14 . In the following we then restrict ourselves to consider a special initial state for which these overlaps are known. Indeed, an exact formula has been conjectured and tested in dkm-16 in the context of the AdS/CFT correspondence. While the special initial state that we consider is admittedly artificial, we will see that it has a very simple structure. This, in turn, makes it possible to test our theoretical predictions by means of efficient numerical methods.

Before presenting the specific initial state considered and analysing the corresponding overlap formula, let us briefly sketch the main aspects of the Quench Action method which are relevant for our purposes.

III.1 The Quench Action approach

Despite its introduction being still quite recent ce-13 , the Quench Action approach has already been employed for many models and quench settings. It is not our intention here to provide a detailed review of the method, so we simply sketch the key formulas in order to set the notations. A pedagogical introduction can be found in caux-16 .

In the equilibrium case, at finite temperature, a well-established result of the thermodynamic Bethe ansatz is that the thermodynamic properties of an integrable system are encoded in the rapidity distribution functions of a representative eigenstate. The latter are determined as the solution of a set of non-linear integral equations which are derived by minimizing the thermal free energy functional takahashi-99 . Remarkably, Ref. ce-13 showed that the same idea can be applied to study the stationary state describing local observables after a quantum quench. In this case, the thermal free energy is replaced by another functional called Quench Action. In complete analogy to the thermal case, the saddle-point rapidity distributions fully characterize the thermodynamic properties of the post-quench steady state. Relevant examples of such properties are the local correlations and the velocities of the quasi-particle excitations.

The Quench Action functional explicitly depends on the initial state $|\Psi_{0}\rangle$ and on some symmetry properties of the latter. In the case considered in this work it reads

[TABLE]

where we have indicated with compact notation ${\boldsymbol{\rho}}$ the sets $\{\rho^{(1)}_{n}\}_{n=1}^{\infty}$, $\{\rho^{(2)}_{n}\}_{n=1}^{\infty}$. Here $S_{\rm YY}[{\boldsymbol{\rho}}]$ is the so called Yang-Yang entropy johannesson-86

[TABLE]

while we will comment in the next subsection on the factor $1/2$ in front of $S_{\rm YY}[{\boldsymbol{\rho}}]$ appearing in (51). The first term in (51) is defined in terms of the thermodynamically leading part of the overlap between the initial state $|\Psi_{0}\rangle$ and the Bethe state corresponding to $|{\boldsymbol{\rho}}\rangle$ as follows

[TABLE]

This expression has to be interpreted as the thermodynamic limit of the logarithm of the overlap between $|\Psi_{0}\rangle$ and a Bethe state which corresponds to the functions $\{\rho^{(1)}_{n}\}_{n=1}^{\infty}$ and $\{\rho^{(2)}_{n}\}_{n=1}^{\infty}$. The post-quench steady state is then described by the distributions $\bar{\rho}^{(r)}_{n}(\lambda)$ which are the saddle-point of the Quench Action (51), namely the solution of the system of equations

[TABLE]

In the next subsection, we will explicitly write the Quench Action for the initial state studied in this work, and derive the corresponding saddle-point equations.

III.2 The initial state and the saddle-point equations

We consider the following initial state for our quench problem

[TABLE]

Here we used the notations (4) and (5), while $\sigma^{\alpha}$ are the Pauli matrices

[TABLE]

so that the trace in (55) is over the auxiliary space $h_{0}\simeq\mathbb{C}^{2}$. Finally, the normalization $\mathcal{N}=3^{L}+3(-1)^{L}$ is chosen such that $\braket{\Psi_{0}}{\Psi_{0}}=1$.

The initial state (55) is a matrix product state with local dimension $2$ which satisfies cluster decomposition. The simple structure of $\ket{\Psi_{0}}$ allows for the investigation of its time evolution by means of efficient numerical methods such as DMRG and iTEBD simulations wf-04 ; dksv-04 ; schollwock-05 . Consider a chain of length $L$ and a Bethe state characterized by the rapidities $\{k_{j}\}_{j=1}^{N},\{\mu_{j}\}_{j=1}^{M}$ such that $L$, $N$ and $M$ are even numbers and

[TABLE]

Then, its overlap with the initial state (55) is given by dkm-16

[TABLE]

Here $G_{\pm}$ are Gaudin-like matrices defined by

[TABLE]

where

[TABLE]

with the additional definitions

[TABLE]

An analogous formula exists for the case where $N$ is even, while $L$ and $M$ are odd, and the sets of rapidity distributions still satisfy (57), (58). Conversely, for Bethe states not satisfying (57), (58) the overlap is zero dkm-16 . Formula (59) was conjectured in dkm-16 based on an analogy with the case of the $XXX$ spin-$1/2$ chain, where a similar state can be constructed and the corresponding overlaps computed LeKZ15 ; BLKZ16 ; fz-16 .

It is not simple to extract from (59) the thermodynamically leading part of the overlap. In fact, due to divergences arising in the matrices $G_{\pm}$, one should take into account finite-size deviations from perfect strings (21), (22). Note however that the situation is completely analogous to the one encountered in other models displaying bound states. This is, for instance, the case for quenches from the Néel state to the $XXZ$ spin-$1/2$ chain wdbf-14 ; PMWK14 or from non-interacting to attractive one-dimensional Bose gases pce-16 . Following these works, it can be argued that the ratio of the determinants in (59) only gives a sub-leading contribution in the thermodynamic limit and can thus be neglected. Given the similarity of the argument, we do not report it here, and refer to wdbf-14 ; PMWK14 ; pce-16 for more details.

Dropping the ratio of the determinants, it is straightforward to take the thermodynamic limit of (59). Since the calculations are analogous to the ones performed in wdbf-14 ; pce-16 , here we only report the final result, which reads

[TABLE]

where we defined

[TABLE]

From (51), we see that we have now all the elements necessary to explicitly write down the saddle-point equations (54). Note that the factor $1/2$ in front of the Yang-Yang entropy (52) comes from the fact that the initial state $|\Psi_{0}\rangle$ has non-vanishing overlaps only with parity-symmetric Bethe states dwbc-14 . Putting all together, we can explicitly perform the functional derivative in (54) and obtain the desired saddle-point equations. By exploiting the Bethe-Takahashi equations (37) and (38) we finally obtain

[TABLE]

where the functions $a_{n,m}$ and $b_{n,m}$ are defined in (40) and (41) respectively.

Once again, Eqs. (72) and (73) can be cast in a partially decoupled form which is more suitable for numerical evaluation, i.e.,

[TABLE]

where $s(\lambda)$ is defined in (46), and we introduced

[TABLE]

with the convention

[TABLE]

Interestingly, Eqs. (72) and (73) can be cast in yet another form, which is still partially decoupled. We report it in Appendix C, together with the derivation of (74) and (75).

One can immediately see that (72) are (73) are symmetric under exchanging the particle species $(1)\leftrightarrow(2)$. This observation allows us to simplify the saddle-point equations as follows. Suppose we find a solution $\eta^{(1)}_{n}(\lambda)=\Theta_{n}(\lambda)$, $\eta_{n}^{(2)}(\lambda)=\Xi_{n}(\lambda)$. If $\Theta_{n}(\lambda)\neq\Xi_{n}(\lambda)$, then we find another solution $\eta^{(1)}_{n}(\lambda)=\Xi_{n}(\lambda)$, $\eta^{(2)}_{n}(\lambda)=\Theta_{n}(\lambda)$. We rule out this possibility by assuming uniqueness of the solution of (72) and (73). Hence, we conclude that the two sets of functions coincide, namely

[TABLE]

As a consequence, one can write a unique set of non-linear integral equations for $\eta_{n}(\lambda)$. From (43) and (44) it follows

[TABLE]

The corresponding root densities $\rho^{(1)}_{n}(\lambda)$ and $\rho^{(2)}_{n}(\lambda)$ are found by solving the Bethe-Takahashi equations (37) – (38). Note that, even if $\eta^{(1)}_{n}(\lambda)=\eta^{(2)}_{n}(\lambda)$, the root densities will generically be different due to the asymmetric form of (37) – (38).

The single set of decoupled equations (79) easily allows us to understand the asymptotic behavior of $\eta_{n}(\lambda)$ for large $\lambda$ and $n$. Defining

[TABLE]

and taking the limit $|\lambda|\to\infty$ in (79), one obtains the set of algebraic equations

[TABLE]

with $\eta_{0,\infty}=0$. It is straightforward to verify that the following ansatz satisfy (81)

[TABLE]

As we discuss in the following, Eq. (82) is recovered by our numerical solution of (79).

IV The post-quench steady state

Our strategy to numerically determine the saddle-point root densities is straightforward. First, we solve Eq. (79) for $\eta^{(r)}_{n}(\lambda)$ and then we find $\rho^{(r)}_{n}(\lambda)$, $\rho^{(r)}_{h,n}(\lambda)$ by solving the partially decoupled Bethe equations (44) – (45). For the sake of presentation, details on our numerical solution of (79) and (44) – (45) are postponed to Sec. IV.2, while here we present and discuss the final result.

The rapidity distribution functions $\rho_{n}^{(1)}(\lambda)$, $\rho_{n}^{(2)}(\lambda)$ characterizing the post-quench steady state are displayed in Fig. 1. Note that we have rescaled the rapidity distributions corresponding to bound particles, as they are significantly smaller than those of unbound particles. The knowledge of these rapidity distributions in principle allows us to compute the long-time limit of any local observable after the quench.

The bound-state content of the post-quench steady state is displayed in Fig. 2. The density of unbound particles is the prominent one, even though finite densities of $n$-particle bound-states are non-negligible for small $n$. Also note that the sequence of densities is not monotonic in $n$, but displays an even/odd effect.

The post-quench steady state lies in the same magnetization sector of the ground-state of the model: they both have $D^{(1)}=2/3$ and $D^{(2)}=1/3$ johannesson-86 . The ground-state, however, displays absence of bound-states so that ${\rho^{(1)}_{n}(\lambda)\equiv\rho^{(2)}_{n}(\lambda)\equiv 0}$ for $n\geq 2$. A comparison between the rapidity distributions $\rho^{(1)}_{1}(\lambda)$, $\rho^{(2)}_{1}(\lambda)$ of the ground-state and the steady state is displayed in Fig. 3. We see that even though higher bound-states have small densities in the steady state they significantly influence the rapidity distribution $\rho^{(2)}_{1}(\lambda)$ of the second species of particles.

In the next subsection, we discuss on the computation of the local conserved operators both on the initial and the post-quench steady state. This will be crucial in order to test the validity of Eq. (79) and the accuracy of our numerical solution. Next, we provide further details on the numerical scheme employed to solve the saddle-point integral equations.

IV.1 The local conserved charges

Since our system is integrable there exists an infinite number of local and quasi-local conserved operators, or charges, commuting with the Hamiltonian. Let us focus only on the local charges and indicate them as $\{Q_{n}\}_{n}$, where $Q_{2}=H$ by convention. Since these operators are conserved, their expectation values on the initial state and on the long-time stationary state have to be equal. This fact provides the basis for the main test of the validity of our results.

In Appendix A we derive the following expression for the expectation value of a given charge on a Bethe state

[TABLE]

where

[TABLE]

We indicated with $\ket{\{\rho^{(r)}_{n}\}}$ a Bethe state which corresponds to the rapidity distributions $\{\rho^{(r)}_{n}\}$ in the thermodynamic limit. Note that the second species of particles does not contribute to the value of any of the local conserved charges. Equations (83), (84) and (85) immediately allow us to numerically compute, after integration, the value of the charges on the post-quench steady state.

In order to compute the expectation value of the local charges on the initial state, we exploit the method outlined in fe-13 for the case of the XXZ spin-$1/2$ chain. As discussed in Appendix D, we find that the same method can straightforwardly be applied also in our case. First, we define the generating function $\Omega_{\Psi_{0}}(\lambda)$ such that

[TABLE]

Then, we show (cf. Appendix D) that the generating function has the following simple expression

[TABLE]

From (87) we immediately obtain

[TABLE]

together with the explicit expression of the even charges. As an example we report the first few charges, i.e.,

[TABLE]

We can also readily write down two additional local conserved charges which are independent from the operators $Q_{n}$. These are $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ defined in (7) and (8). The expression for their expectation value on Bethe states is given in (47) and (48). In addition, by exploiting the simple matrix product form of the initial state, it is easy to compute that

[TABLE]

As we will discuss in the next subsection, all these values are correctly recovered by our numerical solution for the steady state rapidity distributions (79).

IV.2 The numerical solution

The partially decoupled version of TBA integral equations is in general more convenient from the numerical point of view. Nevertheless, as a general rule we have always solved both the coupled and partially decoupled versions for all the integral equations that we have considered. The agreement between the two results gives us a useful check for our numerical methods.

All the systems considered in this work feature an infinite number of equations, therefore, to provide a numerical solution one needs to “truncate” them retaining only a finite number $N_{\rm eq}$ of equations. The solution to the truncated system, indicated by $X_{n}(\lambda)$, is then an approximation to the exact result and becomes exact in the limit $N_{\rm eq}\to\infty$. The partially decoupled version of the integral equations allows us to consider larger $N_{\rm eq}$ and usually allows to reach better accuracy. The agreement between the numerical solution of the coupled and partially decoupled equations is almost perfect for the functions $X_{n}(\lambda)$ with $n<N_{\rm eq}$ and $n$ not too close to $N_{\rm eq}$, while small discrepancies arise for $X_{n}(\lambda)$ with $n\lesssim N_{\rm eq}$ due to the effect of truncation. In the following we briefly comment on our numerical solution of the partially decoupled equations, which yielded our most accurate results.

First, we solved the saddle point equations (79) for $\eta_{n}^{(r)}(\lambda)$ using successive iterations. The truncation procedure followed here is analogous to that of Refs. wdbf-14 ; PMWK14 . Since the driving terms $d_{n}(\lambda)$ in (79) are all equal for $n$ of a given parity, it is natural to expect the same even/odd effect to be present for the solutions $\eta_{n}(\lambda)$. In particular, one expects that for large $n$ it holds $\eta_{n}/\eta_{n+2}\sim 1$. More precisely, based on the asymptotic behavior (82), we truncated the equations (79) using the condition

[TABLE]

We checked that for fixed $n$ and $\lambda$ the value $\eta_{n}(\lambda)$ converges to a well defined number increasing $N_{\rm eq}$. In addition, the solution obtained using (94) is consistent with the one of the coupled version of the integral equations (79).

The second step of the numerical solution is to apply the same iterative technique to the partially decoupled Bethe-Takahashi equations (43)-(44) determining $\rho_{n}^{(r)}(\lambda)$. The truncation in this case is performed in the standard fashion

[TABLE]

Note that, while the functions $\eta_{n}^{(r)}(\lambda)$ and $\rho_{n}^{(r)}(\lambda)$ live on the real line, their numerical realization is taken on an evenly distributed mesh of $N_{\mathrm{points}}$ points on a finite interval $[-\Lambda,\Lambda]$. Extrapolation of $\eta_{n}^{(r)}(\lambda)$ and $\rho_{n}^{(r)}(\lambda)$ in $N_{\mathrm{points}}\rightarrow\infty$ is necessary to obtain an accurate solution, while the other two parameters $\Lambda$ and $N_{\rm eq}$ are fixed to reasonably large values. The results shown in the figures are computed using $\Lambda=50$ and $N_{\mathrm{eq}}=40$.

We tested our numerical solution by comparing the theoretical expectation values of the conserved charges (89)-(91) to the corresponding numerical results obtained via (83). After the extrapolation of $\rho_{n}^{(r)}(\lambda)$ to $N_{\mathrm{points}}\rightarrow\infty$, the numerical values of the charges have less than $1\%$ error. The expectation values of $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ are slightly less accurate as they receive significant contributions from higher strings whose densities are affected by the truncation in the number of equations. Choosing the parameters as specified above, the errors in the expectation values of the particle numbers are $3.6\%$ for the first species and $6.2\%$ for the second. As an additional test of the numerical solution, we checked that the quench action (51) evaluated at the saddle point is indeed zero, i.e., the two terms appearing in the difference (51) are equal. After the extrapolation of the densities to $N_{\mathrm{points}}\rightarrow\infty$ we find that the difference between the overlap and the entropy term in (51) is less than $1\%$ of each term.

V Entanglement Dynamics and Elementary Excitations

In this section we exploit the knowledge of the post-quench stationary state, determined in the above section, to investigate the finite-time dynamics of the system after the quench. We focus on the time evolution of the entanglement entropy after the quench. The amount of entanglement between a subsystem $A$ and the rest of the system $\bar{A}$, is measured by the the entanglement entropy $S_{A}(t)$ defined as

[TABLE]

where $\rho_{A}(t)$ is the time evolving density matrix of the system reduced to the subsystem $A$, i.e. $\rho_{A}\equiv{\rm Tr}_{\bar{A}}|\Psi(t)\rangle\langle\Psi(t)|$. This entanglement entropy is known to give very important information about the system, both in and out of equilibrium, see, e.g., the reviews entanglement:review ; ECP:review ; CCD:review ; L:review .

A convenient way to describe the evolution of entanglement is by means of the quasi-particle picture originally proposed in Ref. cc-05 ; cc-06 . In essence, one postulates that the quench creates pairs of correlated quasi-particles in any spacial point of the system. Only pairs created at the same point are correlated and carry entanglement through the system. For $t>0$, the quasi-particles move ballistically in opposite directions and, as a consequence of momentum conservation, the two correlated quasi-particles have opposite velocities $\pm v(\lambda)$, where $\lambda$ is the rapidity parametrising the dispersion relation. When moving through the system, the quasi-particles correlate regions which were initially uncorrelated as pictorially shown in Fig. 4. The entanglement entropy $S_{A}$ is then a weighted (by a function $s(\lambda)$) measure of the number of pairs with one quasi-particle in $A$ and the other in $\bar{A}$. It has been shown in many non-interacting models NR:entharmonicoscillators ; FC:ent ; CTC:entanglement ; EP:entperiodicquench ; chmm-16 , that the predictions of the this quasi-particle picture become exact in the space-time scaling limit of large times and subsystem sizes. Furthermore the qualitative picture for the entanglement entropy evolution has been shown to be correct even in numerical simulations of many interacting integrable and non-integrable models, as e.g. in Refs. dmcf-06 ; lk-08 ; kh-13 ; cce-15 ; nahum-17 . The same picture also provides the entanglement dynamics in local quenches cc-07l ; ep-07 ; sd-11 ; B:localquench and in inhomogeneous situations dsvc17 ; a-17 .

In Ref. ac-16 it has been shown that the quasi-particle picture gives, in the space-time scaling limit, exact predictions even for interacting integrable models, provided that an appropriate choice for the functions $v(\lambda)$ and $s(\lambda)$ is made. One has to introduce multiple species of quasi-particles moving at the velocities $v_{n}(\lambda)$, which are the group velocities of elementary excitations over the stationary state described by $\{\bar{\rho}_{n}(\lambda)\}$. In interacting models the velocities $v_{n}(\lambda)$ are generically state-dependent and non-trivially encode the effects of the interactions – they fulfil a set of integral equations depending on $\{\bar{\rho}_{n}(\lambda)\}$. In the non-interacting limit they reduce to the bare velocities $v^{0}_{n}(\lambda)=\varepsilon_{n}^{\prime}(\lambda)/p_{n}^{\prime}(\lambda)$, where $\varepsilon_{n}(\lambda)$ and $p_{n}(\lambda)$ are respectively the bare energy and momentum. The choice of the quasi-particles’ velocities is in agreement with the one found in transport problems BCDF:transport ; CADY:hydro ; BF:defect and it is ultimately related to the fundamental observation BEL:lightcone that the group velocities of the elementary excitations are the relevant velocities for the propagation of information in interacting integrable models. The natural choice for the weighting functions $s_{n}(\lambda)$ is to set them equal to the Yang-Yang entropy density per rapidity and species. This choice guarantees that the extensive parts of entanglement entropy and thermodynamic entropy coincide at infinite times, in agreement with some general expectations santos-2011 ; gurarie-2013 as well as analytical findings in non-interacting models FC:ent ; CKC:TGLL ; KBC:Ising .

In our case, using the quasi-particle interpretation, we can write the entanglement entropy between a subsystem of contiguous spins of length $\ell$ and the rest of the system as

[TABLE]

Here we used that the indices labeling quasi-particles are $n=1,2,\ldots$ and $r=1,2$. The Yang-Yang entropy density $s_{n}^{(r)}(\lambda)$ appearing in (97) is given by

[TABLE]

and velocities $v_{n}^{(1)}(\lambda)$ and $v_{n}^{(2)}(\lambda)$ fulfil the integral equations

[TABLE]

Before discussing the predictions of (97) for the entaglement dynamics we briefly sketch the derivation of Equations (99) – (100), and refer the reader to Appendix E for the detailed calculation.

V.1 Velocities of the Excitations

Elementary excitations over the stationary state $\{\bar{\rho}_{n}^{(r)}(\lambda)\}$ are constructed by considering a microscopic representation of the stationary state, described by a configuration of integers $\{I_{\alpha}^{n},J_{\beta}^{m}\}$ (cf. (23)), and adding or removing one integer. The operation of adding such an integer to those in the sector labelled by $n$ and $r$, which we call excitation of string type $n$ and species $r$, can be interpreted as the addition of a quasi-particle with rapidity $\lambda$, energy $\varepsilon_{n}^{d\,(r)}(\lambda)$ and momentum $p_{n}^{d\,(r)}(\lambda)$. The latter quantities are given by

[TABLE]

The shift functions $F^{11}_{kn}(z|x)$ and $F^{12}_{kn}(z|x)$ appearing here are found by solving a system of integral equations, reported in Equations (176) – (179) of Appendix E. The quantities $\varepsilon_{n}^{d\,(r)}(\lambda)$ and $p_{n}^{d\,(r)}(\lambda)$ are called dressed energy and momentum. They differ from the bare ones $\varepsilon_{n}^{(r)}(\lambda)$ and $p_{n}^{(r)}(\lambda)$ because of the effects of the interaction encoded in $\{F^{1a}_{kn}(z|x)\}_{a=1,2}$. All this has a very simple physical interpretation – adding a new quasi-particle to the system has a feedback on all the others because of the interaction, so its effective dispersion relation changes.

Given the dispersion relation of an excitation, i.e. its dressed energy and momentum, we can find its group velocity from the formula

[TABLE]

Taking the derivative of (101) and (102) and using (176) – (179) we find that

[TABLE]

and that $\{\varepsilon^{d\,(r)\,\prime}_{n}(\lambda)\}_{r=1,2}$ fulfil the following system of integral equations

[TABLE]

Substituting the definition (103) readily gives the system (99) – (100).

For the the purpose of the numerical solution, it is convenient to follow the steps illustrated in Appendix B and write the system (99) – (100) in a partially decoupled form

[TABLE]

These integral equations can be readily solved numerically. In Fig. 5 we report the velocities of the two species of elementary excitations with the first four string types, constructed over the stationary state $\{\bar{\rho}_{n}^{(r)}(\lambda)\}$ (cf. Eq. (54)). From the plot we see that the velocities are odd functions of $\lambda$ with a minimum and a maximum reached for finite values of $\lambda$. The maximal velocity for the propagation of information is given by that of excitations of the first species and string type $1$, $v_{\rm max}=\max_{\lambda}v_{1}^{(1)}(\lambda)$. The maximal velocity for excitations of the second species is given by $v_{\rm max}^{(2)}=\max_{\lambda}v_{1}^{(2)}(\lambda)$ with $v_{\rm max}^{(2)}\approx 0.5\,v_{\rm max}$.

V.2 Entanglement Dynamics

Let us now consider the entanglement dynamics predicted by Eq. (97). In Fig. 6a, we report the time evolution of the entanglement entropy after a quench from the initial state (55). As is customary, we plot $S_{\ell}(t)/\ell$ as a function of the scaling variable $2v_{\rm max}t/\ell$. The plot clearly shows the standard spreading of entanglement entropy cc-05 : there is a linear increase of the entanglement entropy for $t<\ell/2v_{\rm max}$ governed by the fastest quasi-particles, followed by a slow saturation dictated by all the other slower quasi-particles. For the initial state considered, the largest contribution to the entanglement is coming from the fastest quasi-particles: those of species $r=1$ and string type $n=1$ (cf. Fig. 2). This observation is confirmed by the species resolved lines in Fig. 6a, which show that the particles of the first species bring almost twice as much entanglement as those of the second species. A further confirmation comes from the string-type resolved plots in Figs 6b and 6c, which demonstrate that the contribution of bound states is strongly suppressed. A final observation is that the asymptotic value of the entanglement entropy is approximately 0.7 which is smaller than $\log 3\approx 1.1$, the maximal density of entropy per site in the spin-$1$ chain (indeed $\log 3$ corresponds to the density of thermodynamic entropy in the infinite temperature state).

V.3 Mutual Information

The entanglement entropy is not the ideal quantity to highlight the contribution of all the different quasi-particles: the contribution of the quasi-particles bringing more entanglement covers all the others. The contribution of different quasi-particles can be resolved more effectively considering two disjoint intervals. Let us take two spin blocks $A$ and $B$ of length $\ell$, separated by a distance $d$, as depicted in Fig. 4, and focus on the mutual information

[TABLE]

between the two subsystems $A$ and $B$. The mutual information after a quench is also believed to signal the non-integrability and chaotic behavior of a system via the breakdown of the quasi-particle picture abgh-15 ; sm-15 .

In the quasi-particle picture, the mutual information counts all the pairs of quasi-particles with one quasi-particle in $A$ and the other in $B$, as pictorially shown in Fig. 4. Its time evolution is then simply written down generalizing the result of ac-16 to two species of excitations

[TABLE]

where $\chi_{[a,b]}(x)$ is the characteristic function of $[a,b]$, i.e. it is equal to $1$ if $x\in[a,b]$ and equal to [math] otherwise.

The time evolution of the mutual information is reported in Fig. 7, where we plot $I_{A:B}(t)/\ell$ as a function of $2v_{\rm max}t/\ell$ for three different values of the ratio $d/\ell$. We see that the contributions of different quasi-particles are easily detected as they give rise to peaks in $I_{A:B}(t)$ – the peak due to the quasi-particles of species $r$ and string type $n$ is appearing at approximately

[TABLE]

where we introduced $v_{\max,n}^{(r)}=\max_{\lambda}v_{n}^{(r)}(\lambda)$. Once again, the most prominent peak corresponds to the fastest quasi-particles, as they carry the largest amount of correlations. As shown in Figs. 7a –  7c, increasing the ratio $d/\ell$ we can increase $2v_{\rm max}t^{(1)}_{1}/\ell$ and separate the peaks, in this way it is easier to discern the contribution of different quasi-particles.

VI Conclusions

In this paper we studied a quantum quench in a nested Bethe-ansatz solvable model, the spin-$1$ Lai-Sutherland chain. The thermodynamic description of the latter is in terms of two different species of quasi-particles, each one forming an infinite number of bound states. We considered a simple initial matrix product state for which the overlaps with Bethe states are known dkm-16 and determined the corresponding post-quench steady state by means of the Quench Action approach. We tested the validity of our result by checking that the expectation value of all the local conserved charges on the initial and long-time steady state are equal. Finally, we investigated the post-quench entanglement dynamics using a recently proposed conjecture for interacting integrable models  ac-16 . Importantly, due to the simple structure of our initial state, all our predictions for the entanglement dynamics can in principle be tested against efficient numerical methods such as tDMRG or iTEBD, on the same lines of what done in Ref. ac-16 .

Our work raises a number of compelling questions. From the physical point of view, it is natural to wonder how the two different species of excitations affect the dynamics of local correlations after the quench. It is clear from our results that each species carries a non-vanishing contribution on entanglement dynamics. As a consequence, it is not possible to individually probe the dynamics of the single species by computation of entanglement entropy or mutual information. An interesting question then is whether some appropriate local correlators can be found, such that their dynamics is entirely determined by one and only one of the two species. This would result in a mechanism of effective separation in the spreading of correlations after the quench.

From the technical point of view, an interesting question stems from the structure of the overlap formula of our initial state. Since the seminal works kp-12 ; pozsgay-14 ; dwbc-14 , determinant expressions of the form (59) have appeared in several cases in the past few years for a variety of models. It is even more surprising to encounter a formula of this kind in the nested spin chain studied in this work, possibly signalling the presence of an hidden structure common to all Bethe ansatz integrable models. Characterising the states for which such a formula exists is an important open problem, with far reaching consequences, both from the purely theoretical and practical point of view.

Our work also calls for a systematic analysis of the quasi-local charges in the spin-$1$ Lai-Sutherland chain. This is a necessary step towards the explicit construction of the complete GGE in this model. Such a construction would in turn lead to the possibility of considering a wider set of initial states, overcoming the technical difficulty of the computation of the overlaps required by the Quench Action approach. The explicit calculations presented in this work for the initial state (55), provide a strict test for any future construction of the complete GGE.

Finally, it would be highly desirable to find explicit formulae for calculating expectation values of local operators on the stationary state in the spin-$1$ Lai Sutherland chain. Such formulae are currently known only for few non-nested Bethe ansatz integrable models, namely the Lieb-Liniger model KCI:correlatorsLL11 ; P:correlatorsLL11 , the XXZ spin-$1/2$ chain MP:correlatorsXXZ14 ; P:correlatorsXXZ16 , and the sinh-Gordon field theory BePC16 ; NS:sinhgordon . A promising route is to first consider expectation values at finite temperature and then generalize them to arbitrary excited eigenstates. This approach has been found fruitful in the case of the XXZ spin-$1/2$ chain and of the sinh-Gordon field theory.

VII Acknowledgements

We are grateful to Fabian Essler for inspiring discussions. BB and PC acknowledge the financial support by the ERC under Starting Grant 279391 EDEQS.

Appendix A The algebraic Bethe ansatz

In this section we briefly sketch the algebraic Bethe ansatz analysis of the Hamiltonian (1), which also allows for a systematic derivation of higher local conservation laws.

The fundamental object of the algebraic Bethe ansatz is the so called $R$-matrix, which in our case reads

[TABLE]

where $\mathcal{P}_{12}$ is the permutation matrix exchanging the local spaces $h_{1}$, $h_{2}$. The $R$-matrix (112) is invariant under $SU(3)$ and it is a simple exercise to show that it satisfies

[TABLE]

where ${\rm id}$ is the identity operator on $h_{1}\otimes h_{2}$. Next, introducing the Lax operator

[TABLE]

one can define the family of transfer matrices

[TABLE]

which depend on the spectral parameter $\lambda$. Crucially, different transfer matrices commute

[TABLE]

and as a consequence one can define a sequence of commuting operators as

[TABLE]

These are usually called charges and can be written down explicitly by exploiting the properties of the $R$-matrix, as nicely explained in faddeev-96 . However, the actual expression becomes increasingly unwieldy for many practical purposes as $n$ increases.

The importance of this algebraic construction lies in the fact that the first charge is equal to the Hamiltonian (1), namely

[TABLE]

Hence, the transfer matrix (116) can be thought of as the generator of the local conservation laws of the model.

As we have mentioned in Section II.2, the knowledge of the rapidities corresponding to a given eigenstate also allows for a direct calculation of expectation values of local charges. Indeed, the analytic expression of the eigenvalue of the transfer matrix (116) on a Bethe state with rapidities $\{k_{j}\}$, $\{\lambda_{j}\}$ is known kr-81 , and one can explicitly write

[TABLE]

where

[TABLE]

and

[TABLE]

with

[TABLE]

Then, it is straightforward to compute the action of a charge $Q_{n}$ on the state $\ket{\{k_{j}\},\{\lambda_{j}\}}$: from the definition (118), the calculation is reduced to taking derivatives of the expression (121). This calculation is completely analogous to the corresponding one in the spin-$1/2$ XXX chain and reads as

[TABLE]

In particular, note that one can neglect the first term in (121) since

[TABLE]

It is also possible to take the thermodynamic limit and evaluate the expectation value of conserved charges on Bethe states corresponding to the distributions $\rho_{n}^{r}(\lambda)$, $r=1,2$. Starting from the expression at finite size, making use of the string hypothesis and taking the thermodynamic limit, one finally obtains

[TABLE]

where

[TABLE]

while we indicated with $\ket{\{\rho^{(r)}_{n}\}}$ a Bethe state which corresponds to the rapidity distributions $\{\rho^{(r)}_{n}\}$ in the thermodynamic limit. Remarkably, we see from (124) and (128) that the second species of quasi-particles does not contribute to any of the higher local conserved charges.

Appendix B Bethe equations in the thermodynamic limit

In this appendix we derive the partially decoupled form of the Bethe-Takahashi equations (43), (44) following the standard Bethe ansatz literature takahashi-99 . We also introduce a series of identities which are used throughout this work.

We start by explicitly writing down the Fourier transforms

[TABLE]

where $s(\lambda)$ and $a_{n}(\lambda)$ are defined in (46) and (42) respectively. These allow us to show, by the convolution theorem, the identity

[TABLE]

where $\ast$ denotes the convolution as in (39). Eq. (131) is valid for $j\geq 1$ provided that one defines $a_{0}(x)=\delta(x)$, with $\delta(x)$ being the Dirac delta function. Defining now

[TABLE]

one can also easily show

[TABLE]

where $b_{nm}$ is defined in (41). The above identities allow us to rewrite the Bethe-Takahahshi equations (37) – (38) in compact notation as

[TABLE]

Here $C_{r,s,}$ is defined by the components of the $2\times 2$ matrix

[TABLE]

Define now

[TABLE]

where

[TABLE]

while $A_{n,m}$ and $b_{n,m}$ are defined in (132) and (41) respectively. One can straightforwardly verify the following identities

[TABLE]

Plugging (138) into (134) for $n\geq 2$, we obtain

[TABLE]

Using again (134) to remove the infinite sum, we readily obtain

[TABLE]

where we used (131). This is precisely (43). Analogously, using identity (139) for $n=1$, one has

[TABLE]

which, removing the infinite sum using (134) and using (130), is cast into

[TABLE]

which coincides with (44).

Appendix C Partially decoupled form of the saddle-point equations

In this appendix we derive the partially decoupled form of the saddle-point equations (72) – (73).

Using the identities (138) – (139), one can immediately derive from (72) – (73)

[TABLE]

where we used the definition (137), while

[TABLE]

where $g_{n}$ is given by (70) and $g_{0}\equiv 0$. The kernels (146) can be simplified by systematic use of the convolution theorem. Using some known results on the Fourier transform of generalized functions vladimirovbook , we find the following identity

[TABLE]

Here $\gamma$ is the Euler constant $\gamma\simeq 0.577$ and $\mathcal{P}[1/|x|]$ is such that

[TABLE]

with $f(x)$ smooth.

The identity (147) allows us to establish

[TABLE]

where $f_{n}(\lambda)$ is defined in (71). Then, in the case of $n$ odd, we find

[TABLE]

Analogously, if $n$ is even we have

[TABLE]

Finally, using again (147) we obtain

[TABLE]

This proves (76).

For completeness, we now show how to cast (144) – (145) in yet another form which is similar to the integral equations derived in the finite temperature case jls-89 . The derivation closely follows Ref. afl-83 . First, we define

[TABLE]

where $s^{-1}[\mathcal{G}](\lambda)$ is formally defined as

[TABLE]

namely it is the inversion of the integral functional with kernel $s$ defined in (46). In (154) we used the convention

[TABLE]

Equations (144) – (145) can then be rewritten as

[TABLE]

where $\bar{r}$ is defined in (137). A formal solution of (157) can be obtained in Fourier transform. It is given by

[TABLE]

Using the expression for the inverse Fourier transform gr-14

[TABLE]

one can directly take the inverse Fourier transform of (158) – (159). Defining

[TABLE]

and using

[TABLE]

we obtain

[TABLE]

with the convention (156) and where the driving terms are given by

[TABLE]

This integral can be performed explicitly. We use gr-14

[TABLE]

Choosing $b=2\pi$, $a=4\pi/3$, and exploiting the properties of the Fourier transform, one arrives at

[TABLE]

Note that (165) – (166) are exactly of the same form of the decoupled thermal equations reported in jls-89 . As expected from the experience with previous applications of the Quench Action, the only difference is in the driving terms $w_{n}(\lambda)$.

Appendix D Higher local conserved charges

In this appendix we discuss the computation of the initial-state expectation values of local conserved charges. More specifically, we detail the derivation of formula (87) for the corresponding generating function.

As we anticipated in the main text, we exploit the method used in fe-13 for the case of the XXZ spin-$1/2$ chain. The starting point of our derivation is provided by the formal expansion of the transfer matrix (116) in terms of the local charges (118), namely

[TABLE]

Since the the operators $Q_{k}$ are self-adjoint, Eq. (170) suggests that for large $L$ one can write

[TABLE]

in the sense that the power series expansions in $\lambda$ of the two sides of (171) coincide for $\lambda\in\mathbb{R}$. Next, from (86) we observe that the generating function $\Omega_{\Psi_{0}}(\lambda)$ can be defined as

[TABLE]

Using (171) we then obtain

[TABLE]

The computation is then reduced to evaluating the expectation value of the product of two transfer matrices. This can be performed analytically, due to the simple structure of our initial state and exploiting the representation of transfer matrix (116) as a matrix product operator. In fact, the steps involved for the explicit evaluation of the r.h.s. of (173) are completely analogous to those explained in fcec-14 for the case of the spin-$1/2$ chain and will not be reported here. Note in particular that the derivative in (173) can be efficiently performed by means of the Jacobi formula fcec-14 , which makes it possible to reach a closed form analytical result. The intermediate analytical calculations can be easily carried out with the program Mathematica, so that one finally derives Eq. (87).

Appendix E Elementary excitations over a stationary state

Here we consider excitations on a stationary state described by a set of root densities $\{\rho_{n}^{(1)}(k),\rho_{n}^{(2)}(\lambda)\}$. We consider in a large, finite volume $L$ and take the system in the stationary state described by the set of integers $\{I_{\alpha}^{n},J_{\beta}^{m}\}$, corresponding to a set of string centers $\{k^{n}_{\alpha},\lambda_{\beta}^{m}\}$ through Eqs. (23). Using the counting functions (24) – (25) we can establish a one to one correspondence between the integers $\{I_{\alpha}^{n},J_{\beta}^{m}\}$ and the rapidities $\{k^{n}_{\alpha},\lambda_{\beta}^{m}\}$. The integers $\{I_{\alpha}^{n},J_{\beta}^{m}\}$ are chosen such that the distribution of string centers describing the state in the thermodynamic limit are $\{\rho_{n}^{(1)}(k),\rho_{n}^{(2)}(\lambda)\}$.

Excitations on the state are constructed by changing the value of a finite number of integers in $\{I_{\alpha}^{n},J_{\beta}^{m}\}$. The corresponding modified set of rapidities can be described as follows. The rapidities corresponding to the unchanged integers $\{\tilde{k}^{n}_{\alpha},\tilde{\lambda}_{\beta}^{m}\}$; the ones of particle excitations $\{k^{p}_{n},\lambda_{m}^{p}\}$, corresponding to the added integers; those of hole excitations $\{k^{h}_{n},\lambda_{m}^{h}\}$, corresponding to the removed integers. The latter set is composed by fictitious rapidities, cancelling the contribution of those corresponding to the removed integers. The rapidities $\{\tilde{k}^{n}_{\alpha},\tilde{\lambda}_{\beta}^{m}\}$ in the state with $n^{(r)}_{p,m}$ particle excitations and $n^{(r)}_{h,m}$ hole excitations for every species $r=1,2$ and string type $m=1,\ldots$ are related to $\{k^{n}_{\alpha},\lambda_{\beta}^{m}\}$ as

[TABLE]

where we introduced the shift functions $F^{11}_{nm}(x|y),F^{12}_{nm}(x|y),F^{21}_{nm}(x|y)$, $F^{22}_{nm}(x|y)$ and assumed the total number of excitations to be even. The shift functions are defined through the following integral equations

[TABLE]

where the terms $\Xi_{n,m}$ and $\Theta_{n,m}$ are defined in (27) and (28), while the kernels $a_{n,m}$ and $b_{n,m}$ are given in (40) and (41). These equations are derived following the same procedure employed in the study of excitations on thermal states kbi-93 . Using the Bethe equations for the state described by $\{I_{\alpha}^{n},J_{\beta}^{m}\}$ and those for the state with the excitations one finds

[TABLE]

Here $z^{(a)}(x)$ are the counting functions of the state, while $z_{\rm ex}^{(a)}(x)$ are those of the state with the excitations. Taking the thermodynamic limit of (180) and using the definition (174), (175) one readily finds that the shift functions satisfy Eqs. (176) – (179).

One can easily compute the expectation value on excited states of the local conserved charges discussed in appendix D. Using (124), one can introduce for a given local conserved operator $Q$ the corresponding bare charge densities $q_{n}(\lambda)$ by

[TABLE]

Then, the difference in the expectation value of the excited state and the one without excitations is given by

[TABLE]

Here we introduced the dressed charges

[TABLE]

where $q^{\prime}_{k}(z)=({\rm d}/{{\rm d}z})q_{k}(z)$. In particular, dressed energy and momentum are given by

[TABLE]

where $p_{n}(\lambda)$ and $\varepsilon_{n}(\lambda)$ are given in (26) and (32) respectively. We are now in a position to formally solve the equations (176)–(179), finding a closed expression for the dressed charges. In doing so, it is convenient to introduce the following compact notations BCDF:transport . We represent functions $w_{j}(\lambda)$ by a vector $\vec{w}$, such that

[TABLE]

Kernels $A_{jk}(\lambda,\mu)$ which are functions of two rapidities with two string indices are instead represented by an operator $\hat{A}$ acting as follows

[TABLE]

The inverse of $\hat{A}$ is the operator $\hat{A}^{-1}$ satisfying

[TABLE]

In this way, the distribution $\delta(\lambda-\mu)\delta_{ij}$ corresponds to the identity $\hat{1}$. It is also convenient to define diagonal operators $\hat{w}$ associated with a function of a single rapidity

[TABLE]

Rewriting (176)–(179) in compact notations we have

[TABLE]

where $\hat{\vartheta}^{(r)}$ are diagonal operators

[TABLE]

while $\hat{a}$, $\hat{b}$ and $\hat{F}^{rs}$ are non diagonal

[TABLE]

Equations (191) - (194) are coupled in pairs. Let us consider the first pair composed by (191) and (193), from the latter we find

[TABLE]

Plugging into (191) we have

[TABLE]

so that

[TABLE]

Analogously, considering (192) and (194) we have

[TABLE]

In compact notation, the equations for the dressed charges read as

[TABLE]

where $t$ is the transposition operation. For a given kernel $A$, acting as in (188), the transposition is defined by

[TABLE]

Let us consider $a=1$ and plug (199) into (201)

[TABLE]

where in the first step we used that $\hat{\Theta}$ and $\hat{\Xi}$ are antisymmetric, $\hat{a}$ and $\hat{b}$ are symmetric and $\hat{\vartheta}^{(r)}$ diagonal. Let us now consider the equation “in components” $[\cdot]_{j}(\lambda)$ and take the derivative with respect to $\lambda$. Re-expressing then everything in compact notations we have

[TABLE]

where we used

[TABLE]

Analogously we find

[TABLE]

Applying our general results to the case of dressed momenta we have

[TABLE]

where the second equality follows by inverting Eqs. (37) – (38). We can finally write an integral equation for the group velocities of the quasi-particle excitations. These are defined as BCDF:transport ; BEL:lightcone

[TABLE]

In compact notations we then have

[TABLE]

Note now that these equations have the exact same form of those obtained by inverting Eqs. (37) – (38), the driving term $a$ being replaced by $\varepsilon^{\prime}/(2\pi)$. Then, it is straightforward to bring these equations to the following final explicit form

[TABLE]