Correlation and entanglement spreading in nested spin chains

TL;DR

This paper numerically investigates correlation and entanglement spreading in nested integrable spin chains, specifically the Lai-Sutherland model, validating theoretical predictions about quasiparticle dynamics and light-cone velocities after quantum quenches.

Contribution

It provides the first direct numerical validation of Bethe ansatz and quasiparticle predictions in nested integrable systems like the Lai-Sutherland model.

Findings

Confirmed the quasiparticle-based formula for entanglement growth.

Validated Bethe ansatz predictions for correlation light-cone velocity.

Demonstrated the effectiveness of DMRG and exact diagonalization in nested systems.

Abstract

The past few years have witnessed the development of a comprehensive theory to describe integrable systems out of equilibrium, in which the Bethe ansatz formalism has been tailored to address specific problems arising in this context. While most of the work initially focused on the study of prototypical models such as the well-known Heisenberg chain, many theoretical results have been recently extended to a class of more complicated $\it{nested}$ integrable systems, displaying different species of quasiparticles. Still, in the simplest context of quantum quenches, the vast majority of theoretical predictions have been numerically verified only in systems with an elementary Bethe ansatz description. In this work, we fill this gap and present a direct numerical test of some results presented in the recent literature for nested systems, focusing in particular on the Lai-Sutherland model. Using time-dependent density matrix renormalization group and exact diagonalization methods, we compute the spreading of both correlation functions and entanglement entropy after a quench from a simple class of product initial states. This allows us to test the validity of the nested version of a conjectured formula, based on the quasiparticle picture, for the growth of the entanglement entropy, and the Bethe ansatz predictions for the "light-cone" velocity of correlation functions.

Equations

$$H_L=\sum _{j=1}^{L-1}{\text{s}}_j.{\text{s}}_{j+1}+({\text{s}}_j.{\text{s}}_{j+1}{)}^2-2L,$$

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

$${\left(\frac{k_j+i/2}{k_j-i/2}\right)}^L=\prod _{\genfrac{}{}{0ex}{}{p=1}{p\ne j}}^N\frac{k_j-k_p+i}{k_j-k_p-i}\prod _{\ell =1}^M\frac{{\lambda }_{\ell }-k_j+i/2}{{\lambda }_{\ell }-k_j-i/2},j=1,\dots ,N,$$

$$1=\prod _{j=1}^N\frac{k_j-{\lambda }_{\ell }-i/2}{k_j-{\lambda }_{\ell }+i/2}\prod _{\genfrac{}{}{0ex}{}{m=1}{m\ne \ell }}^M\frac{{\lambda }_{\ell }-{\lambda }_m-i}{{\lambda }_{\ell }-{\lambda }_m+i},\ell =1,\dots ,M.$$

$$E=-\sum _{j=1}^{N}1/(k_j^2+1/4).$$

$${\rho }_{t,n}^{(1)}(\lambda )$$

$${\rho }_{t,n}^{(2)}(\lambda )$$

$$a_{n,m}(\lambda )$$

$$b_{n,m}(\lambda )$$

$${\rho }_{t,n}^{(2)}(\lambda )v_n^{(2)}(\lambda )$$

$${\rho }_{t,n}^{(1)}(\lambda )v_n^{(1)}(\lambda )$$

$$|{\Psi }_0⟩=|{\psi }_0{⟩}_{1,2}\otimes |{\psi }_0{⟩}_{3,4}\otimes ....\otimes |{\psi }_0{⟩}_{L-1,L},$$

$$|{\psi }_0⟩=c_{11}|1,1⟩+c_{22}|2,2⟩+c_{33}|3,3⟩+c_{12}(|1,2⟩+|2,1⟩)+c_{13}(|1,3⟩+|3,1⟩)+c_{23}(|2,3⟩+|3,2⟩).$$

$$|{\psi }_0⟩=c_{11}|1,1⟩+c_{22}|2,2⟩+c_{33}|3,3⟩,$$

$$S^{\alpha }(t)=\frac{1}{1-\alpha }\ln \text{mbox}Tr{\rho }_A^{\alpha }(t),$$

$$S_l(t)\propto 2t{\int }_{2|v|t<l}d\lambda |v(\lambda )|s(\lambda )+l{\int }_{2|v|t>l}d\lambda s(\lambda ),$$

$$\underset{\begin{array}{c}t,l\to \infty \\ l/t\mathrm{fixed}\end{array}}{\lim }{S}_l(t)/l=(S_l^{(1)}+S_l^{(2)})/l,$$

$$S_l^{(r)}/l=\sum _{n=1}^{\infty }\int \mathrm{d}\lambda s_n^{(r)}(\lambda )\left\{2\frac{t}{l}|v_n^{(r)}(\lambda )|{\theta }_{\mathrm{H}}\left(\frac{l}{t}-2|v_n^{(r)}(\lambda )|\right)+{\theta }_{\mathrm{H}}\left(2|v_n^{(r)}(\lambda )|-\frac{l}{t}\right)\right\},$$

$$s_n^{(r)}(\lambda )=\left({\rho }_n^{(r)}(\lambda )+{\rho }_{h,n}^{(r)}(\lambda )\right)\ln \left({\rho }_n^{(r)}(\lambda )+{\rho }_{h,n}^{(r)}(\lambda )\right)-{\rho }_n^{(r)}(\lambda )\ln {\rho }_n^{(r)}(\lambda )-{\rho }_{h,n}^{(r)}(\lambda )\ln {\rho }_{h,n}^{(r)}(\lambda ),$$

$$I_{A:B}(t)=$$

$$+\left(d+2l-2\left|v_n^{(r)}(\lambda )\right|t\right){\chi }_{[d+l,d+2l]}\left(2\left|v_n^{(r)}(\lambda )\right|t\right)]s_n^{(r)}(\lambda ),$$

$$⟨{\Psi }_0⟩O_{L/2}(t){\Psi }_0=⟨{\Psi }_0⟩O_{L/2+1}(t){\Psi }_0.$$

$$⟨{\Psi }_0⟩O_j(t){\Psi }_0=\frac{1}{L}\sum _{k=1}^L⟨{\Psi }_0⟩O_k(t){\Psi }_0.$$

$$\frac{S_l}{l}=s_{\infty }+\frac{a}{l}+\frac{b}{l^2},$$

Full text

I Introduction

After more than a decade of intense investigations, the theory of many-body integrable systems out of equilibrium is arguably entering a more mature stage, as several experimentally feasible protocols are now falling under the reach of our theoretical understanding CaEM16 ; KiWW06 ; LaGS16 ; LEGR15 ; SBDD19 ; FSEH13 . In this respect, the study of the quantum quench cc-06 , in which a well-defined initial state is left to evolve unitarily according to some known Hamiltonian, has played a major role for the development of the theory. Indeed, besides offering a conceptually simple framework to investigate questions of fundamental interest, the need to overcome the significant computational challenges of this problem led to the introduction of powerful methods EsFa16 ; Caux16 ; IMPZ16 , which now represent significant milestones in the field.

One of the most important ones is arguably the Quench Action approach introduced in Ref. CaEs13 , and first applied in Ref. DWBC14 to a genuinely interacting model (see Ref. Caux16 for a review). This method provided an intuitive picture to describe the local properties of the system at large times after a quench: according to the Quench Action approach, the unitary dynamics brings the system to a stationary state whose local properties can be computed from the knowledge of a single representative eigenstate of the Hamiltonian. Importantly, the latter can be characterized in terms of the corresponding quasimomentum (or rapidity) distribution functions of the stable quasiparticles; in fact, it was later realized BWFD14 ; PMWK14 ; IDWC15 ; PoVW17 ; IlQC17 ; IQDB16 ; PiVC16 ; PVCR17 that this characterization is equivalent to that of the generalized Gibbs ensemble (GGE) RDYO07 ; ViRi16 ; CaEF11 ; FaEs13 ; Pozs13 , a statistical ensemble which generalizes the Gibbs density matrix by taking into account all the local and quasilocal conservation laws of the Hamiltonian Pros11 ; IlMP15 ; PiVe16 ; DeCD17 ; PrIl13 .

The description of the post-quench steady state in terms of the quasiparticle rapidity distribution functions is particularly convenient, as they are directly related to several aspects of the post-quench dynamics. For instance, the “light-cone” velocity of correlation functions was argued to coincide with that of the fastest quasiparticle excitation built on top of the stationary state BoEL14 , while exact formulas have been derived to compute the post-quench steady value of local correlators as a function of these rapidity distributions Pozs11 ; MePo14 ; Pozs17 ; BaPi18 ; BePC16 ; NeSm13 . Furthermore, as a notable result, it was shown in Refs. AlCa17 ; AlCa18 that they also give us access to the exact evolution of the entanglement after the quench (in the scaling limit of asymptotically large times and spin blocks), as the knowledge of the occupation numbers and dispersion relations of the quasiparticles complements the standard picture of Ref. CaCa05 introduced in the context of conformal field theory.

Unfortunately, the actual computation of the rapidity distributions of the quasiparticles remains in general a difficult task. Indeed, the prescriptions of Ref. CaEs13 to compute them require, as an input, the overlaps between the initial state and the eigenstates of the post-quench Hamiltonian, which have proven to be very hard to obtain KoPo12 ; Pozs14 ; cd-14 ; PiCa14 ; Broc14 ; LeKZ15 ; LeKM16 ; HoST16 ; BrSt17 ; Pozs18 ; LeKL18 ; HoKT19 . In order to get around this problem, two different approaches were later introduced: the first one is the “string-charge duality” method IQDB16 ; IDWC15 , based on the idea that the post-quench representative eigenstate could be uniquely fixed from the expectation value of a complete set of quasi-local conserved operators on the initial state. While this method allows us to study quenches from rather general initial states, it is limited to the models for which the latter are known; in particular, so far it could be applied only to XXZ Heisenberg chains IQDB16 ; PVCR17 ; PiVC16 . The second method is the Quantum Transfer Matrix approach introduced in Refs. PiPV17 ; PiPV17_II (see also PiPV18_losch ), which derives the post-quench rapidity distribution functions based on some underlying mathematical structures of integrability. This approach is, by construction, limited to a special class of “integrable” initial states PiPV17_II ; PoPV18 , whose definition is inspired by classical works in integrable quantum field theory Ghos94 (see also Refs. Delf14 ; Schu15 ); still it was immediately realized it could be applied quite generally to different models, relying exclusively on defining properties of integrability PiPV17_II .

Given the significant technical challenges, these important developments were initially restricted to classes of models solvable through an elementary (as opposed to nested) Bethe ansatz (with a notable exception in the context of relativistic quantum field theory BeSE14 ). Still, nested integrable systems are particularly interesting as they contain different species of interacting particles efgk-05 , and include experimentally relevant examples such as, most prominently, multi-component Fermi and Bose gases, as already realized in several pioneering cold atomic experiments BlDZ08 ; guan2013 ; PMCL14 ; exp2 . As a recent development, a series of works was able to extend several of the aforementioned results to the case of nested systems, focusing on the simplest example of the Lai-Sutherland chain lai-74 ; sutherland-75 , a special case of the Perk-Schultz model PeSc81 ; Schu83 . In particular, the first application of the Quench Action method in this system appeared in Ref. MBPC17 , where some of us analyzed the quench from a very special initial state, whose overlaps with the eigenstates of the Hamiltonian were computed in Ref. LeKM16 . Later, the Quantum Transfer Matrix approach was generalized in Refs. PVCP18_I ; PVCP18_II to the nested models, finally allowing us to study the post-quench dynamics for an entire class of integrable initial states. Note that the special state considered in Ref. LeKM16 was also shown to be a particular example in this class PiPV17_II ; PoPV18 .

The theory that has been developing with the introduction of the Quench Action approach relies on several assumptions, although physically very sound, and extensive numerical tests have been necessary to establish its validity beyond reasonable doubt. In particular, the predictions of the theory for the asymptotic values of local correlations have been tested against tDMRG WhFe04 ; DKSV04 and iTEBD Vida07 computations FCEC14 ; PVCR17 ; PMWK14 , as well as by means of exact diagonalization techniques and numerical-linked cluster expansions BWFD14 ; PVCR17 . Analogously, the evolution of the entanglement has been tested extensively against tDMRG simulations in Refs. AlCa17 ; AlCa18 , as well as other entanglement measures AlCal-b ; AlCal-c . However, all of these works focused on the prototypical XXZ Heisenberg chain, while no direct verification was ever performed for the case of nested systems (with the only exception of Fermi-Hubbard models with infinite onsite repulsion BeTC17 ; BeTC18 ; ZhVR18 ; ZhVR19 , which is however simpler, since it can be mapped to a free model). Evidently, it is very important to corroborate at least some of the recent results obtained in nested systems. This is the aim of this work, where we focus on the Lai-Sutherland chain, and compute numerically different quantities after the quench from the family of integrable initial states of Refs. PVCP18_I ; PVCP18_II . This allows us to test directly several analytic predictions formulated in the recent literature for a truly interacting nested intergrable model, thus corroborating the existing theoretical framework for the study of quantum quenches in nested models.

The paper is organized as follows. In Sec. II we introduce the Lai-Sutherland model, and discuss the class of initial states which will be considered in this work. Our numerical results for the evolution of the entanglement are reported in Sec. III, together with a comparison against analytic predictions available in the literature. Sec. IV is devoted to the analysis of the correlation functions after the quench, while our conclusions are consigned to Sec. V. Finally, technical details on our numerical computations are reported in Appendix A.

II The model and the initial states

We consider the Lai-Sutherland model lai-74 ; sutherland-75 , which is described by the Hamiltonian

[TABLE]

where $\textbf{s}_{j}=(s^{x}_{j},s^{y}_{j},s^{z}_{j})$. The spin-$1$ operators $s^{a}_{j}$ act on the local Hilbert space $h_{j}\simeq\mathbb{C}^{3}$, and are given by the standard three-dimensional representation of the $SU(2)$ generators; in particular, introducing the local basis $\ket{1}=(1,0,0)$, $\ket{2}=(0,1,0)$, $\ket{3}=(0,0,1)$, they read

[TABLE]

The calculations presented in this work will be exclusively numerical. However, since they will be compared to analytic predictions, it is useful to recall a few basic aspects of the theoretical tools used to analyze the Hamiltonian (1), which provides one of the simplest models which can be solved via the nested Bethe ansatz method efgk-05 ; kr-81 ; johannesson-86 ; johannesson2-86 . Specifically, by means of the latter it is possible to show that the eigenstates of (1) are parametrized by two distinct sets of quasi-momenta (or rapidities) $\{k_{j}\}_{j=1}^{N}$, $\{\lambda_{j}\}_{j=1}^{M}$, satisfying an appropriate set of quantization conditions (the nested Bethe equations)

[TABLE]

The energy eigenvalues can then be written as

[TABLE]

Physically, the sets $\{k_{j}\}_{j=1}^{N}$ and $\{\lambda_{j}\}_{j=1}^{M}$ correspond to two distinct types, or species, of quasiparticles. This situation is different from that of ordinary Bethe-ansatz solvable models (e.g. the well-know Heisenberg chain kbi-93 ), where the eigenstates of the Hamiltonian are parametrized by a single set of rapidities. We note that the existence of more than one species of quasiparticles is not just a mathematical fact, but have implications on observable quantities, even in out-of-equilibrium protocols MBPC19 .

The thermodynamic description of the model, defined by $L,N,M\to\infty$ keeping the ratios $D_{1}=N/L$ and $D_{2}=M/L$ constant, is analogous to the case of ordinary Bethe ansatz takahashi_book . In particular, in this limit the rapidities arrange themselves in the complex plane according to specific patterns called strings, which correspond to bound-states of the quasiparticles; in each string the rapidities are parametrized as $k^{n,\ell}_{\alpha}=k^{n}_{\alpha}+i\left[(n+1)/2-\ell\right]+\delta_{1,\alpha}^{n,\ell}$, $\lambda^{n,\ell}_{\alpha}=\lambda^{n}_{\alpha}+i\left[(n+1)/2-\ell\right]+\delta_{2,\alpha}^{n,\ell}$, where $\ell=1,\ldots n$. Here the real numbers $k^{n}_{\alpha}$, $\lambda^{n}_{\alpha}\in(-\infty,+\infty)$ are the string centers, which can be interpreted as the quasimomenta of the bound-states, while $\delta_{r,\alpha}^{n,\ell}$ are negligibly small deviations; finally, $n$ is the length of the string (namely, the number of quasiparticles in the bound-state).

In the thermodynamic limit the string centers for the two species become continuous variables on the real line, distributed according to rapidity distribution functions $\rho^{(1)}_{n}(k)$ and $\rho^{(2)}_{n}(\lambda)$. One also needs to introduce the functions $\rho^{(1)}_{h,n}(k)$ and $\rho^{(2)}_{h,n}(\lambda)$ describing the distribution of “holes”, which are the available rapidities for which there is no quasiparticle. These functions are related by the following thermodynamic version of the Bethe equations (see e.g. Ref. MBPC17 )

[TABLE]

Here we employed the standard definition $\rho^{(r)}_{t,n}(k)=\rho^{(r)}_{n}(k)+\rho^{(r)}_{h,n}(k)$, together with $\left(f\ast g\right)(\lambda)=\int_{-\infty}^{\infty}{\rm d}\mu f(\lambda-\mu)g(\mu)$, and

[TABLE]

where $a_{n}(\lambda)=n/[2\pi(\lambda^{2}+n^{2}/4)]$. Finally, from the quasiparticle distribution functions several quantities can be computed directly. Among these, the (dressed) velocities of the quasiparticles (which are once again parametrized by the rapidities $\lambda$) are obtained through the system of integral equations

[TABLE]

Here, $v^{(1)}_{n}(\lambda)$ and $v^{(2)}_{n}(\lambda)$ are the dressed velocities of the $n$-quasiparticle bound states of rapidity $\lambda$, for the first and second species.

II.1 The initial states

As we have mentioned in Sec. I, until recently no initial state was known in nested systems for which an analytical study of the quench dynamics could be carried out. As a first piece of progress, in Refs. LeKM16 ; LeKL18 a family of matrix product states (MPSs) PVWC06 with increasing bond dimension $\xi$ was found, for which the exact overlaps with the eigenstates of the Lai-Sutherland Hamiltonian could be worked out. While the quench dynamics from the simplest case $\xi=2$ could be analyzed in Ref. MBPC17 by means of the Quench Action method, these states have large initial entanglement, which makes them not convenient for numerical simulations. Subsequent studies clarified the status of the special states found in Refs. LeKM16 ; LeKL18 as belonging to a broader class of “integrable” initial states PiPV17_II ; PoPV18 . Furthermore, in Refs. PVCP18_I ; PVCP18_II additional initial product states in this class were found, for which the quasiparticle rapidity distribution functions could be derived exactly. From the numerical point of view, these states are particularly convenient: on the one hand, they display small entanglement at short times; on the other hand, they are parametrized by continuous variables, which offer the possibility of a more complete analysis.

The class of integrable states found in Refs. PVCP18_I ; PVCP18_II read

[TABLE]

with the two-site block

[TABLE]

Here $c_{ij}$ are arbitrary complex numbers. Note that, as it was shown in Ref. PVCP18_I , each state in this class can be rotated via a global $SU(3)$ transformation to a product state with a two-site block of the form

[TABLE]

with the additional restriction $|c_{11}|\geq|c_{22}|\geq|c_{33}|$. In the rest of this work, we will study the quench dynamics from the initial states (11), with two-site blocks given in (13).

Before leaving this section, we recall an important property of these states: their overlaps with an eigenstate of the Hamiltonian is non vanishing only if the latter is labeled by sets of rapidities $\{k_{j}\}$, $\{\lambda_{j}\}$ which are parity invariant, namely $\braket{\{k_{j}\},\{\lambda_{j}\}}{\Psi_{0}}\neq 0$ implies $\{k_{j}\}=\{-k_{j}\}$, $\{\lambda_{j}\}=\{-\lambda_{j}\}$. This property has been chosen as the very definition of integrability in Refs. PiPV17_II ; PoPV18 : there, an initial MPS with finite bond dimension was called integrable if this condition is fulfilled. As we will see in the next section, this property has also consequences on the description of the post-quench entanglement dynamics.

III The entanglement dynamics

We start our analysis with the entanglement dynamics after the quench from an integrable initial state. Specifically, we consider a subsystem $A$ of length $l$, and compute numerically the evolution of the entanglement between $A$ and the rest of the system $B$, as given by the Rènyi entropy

[TABLE]

where $\alpha$ is an arbitrary positive real number. Here $\rho_{A}(t)$ is the time-evolved reduced density matrix of the subsystem $A$, i.e. $\rho_{A}=\mbox{Tr}\,_{B}|\Psi(t)\rangle\langle\Psi(t)|$, where $|\Psi(t)\rangle=e^{-iHt}|\Psi_{0}\rangle$. We will be particularly interested in the limit $\alpha\to 1$, which yields the well-known von Neumann entropy $S_{V}(t)=-{\rm tr}\left[\rho_{A}(t)\log\rho_{A}(t)\right]$.

As we have discussed in Sec. I, an important result in the recent literature has been the discovery in Refs. AlCa17 ; AlCa18 of a formula to compute the evolution of the von Neumann entanglement entropy in the scaling limit $t,l\to\infty$, keeping the ratio $l/t$ fixed. Such formula was derived in Ref. AlCa17 based on an intuitive quasiparticle picture, which goes along the following lines. One interprets the quench as a process generating everywhere and homogeneously uncorrelated pairs of entangled quasiparticles with opposite momenta. In this picture, two regions may be entangled only if they share at least a pair of quasiparticles emitted from an arbitrary initial point. Accordingly, the total entanglement entropy between a region $A$ and its complement $B$ is proportional to the number of pairs with one quasiparticle in $A$ and the other in $B$. Assuming there exists only one type of quasiparticles, this semiclassical description immediately gives us the following expression for the von Neumann entanglement entropy AlCa17

[TABLE]

where $s(\lambda)$ is the contribution to the entanglement carried by the pair of quasiparticles with rapidity $\pm\lambda$, while $|v(\lambda)|$ is the absolute value of their (opposite) velocities. In Ref. AlCa17 it was argued that the first could be identified with the Yang-Yang entropy density YaYa69 while the second could be obtained from the dressed velocity of the quasiparticles as computed from the post-quench rapidity distribution functions BoEL14 . A straightforward generalization of this picture to the case where the quasiparticles can form bound states was employed in Ref. AlCa17 to obtain quantitative predictions for the $XXZ$ Heisenberg chain. The same ideas can be also extended to inhomogeneous situations Al-18 ; abf-19 ; bertini-2018a within the generalized hydrodynamics approach CaDY16 ; BCDF16 .

The analytic formula found in Ref. AlCa17 is now widely believed to be exact for the class of integrable initial states: on the one hand, based on their defining property, the application of the quasiparticle picture is particularly intuitive and natural for these states; on the other hand, the formula for the entanglement growth was extensively tested against tDMRG simulations in Refs. AlCa18 ; AlCa17 and very good agreement was always found for integrable quenches. Conversely, it is expected that these results should be modified for initial states with a more complicated structure of the overlaps, as it was found in Refs. BeTC18 ; bc-18 for the case of free models. It is important to stress that a similar formula for the time evolution of Rényi entropies after a quench to an interacting integrable models is not yet known AlCal-b .

In the case of nested systems, it is natural to assume that a similar picture holds for integrable initial states, and that the quasiparticles of each species and their bound states carry an independent contribution to the entanglement. In the case of the Lai-Sutherland model, this leads to the formula MBPC17 ; PVCP18_I

[TABLE]

where

[TABLE]

and where $\theta_{H}(x)$ is the Heaviside Theta function ($\theta_{H}(x)=0$ if $x<0$, $\theta_{H}(x)=1$ otherwise). Here, $r$ and $n$ are the indexes labeling the different species and bound states of the quasiparticles; $s_{n}^{(r)}(\lambda)$ is the Yang-Yang entropy density

[TABLE]

while $v^{(r)}_{n}(\lambda)$ are the quasiparticle velocities introduced in Eq. (10).

We stress that, even though the above conjecture is natural, it is very important to provide numerical tests of its validity. For example, in nested systems it is known that only the first species of quasiparticles contribute to the expectation value of local conserved quantities MBPC17 ; PVCP18_I ; this is apparent, for instance, for the energy eigenvalues appearing in Eq. (5). Accordingly one could wonder whether the correct generalization of (15) to the Lai-Sutherland model should only involve the first, and not both, species of the quasiparticles.

In order to test the validity of (16) we have computed numerically the evolution of the entanglement entropy $S_{l}(t)$ for different integrable initial states and increasing values of the subsystem size $l$. We have then performed a fit to extrapolate the value of $S_{l}(t)/l$ in the limit $t,l\to\infty$, as explained in detail in Appendix A. The numerical results are then compared against the Bethe ansatz prediction (16); importantly, this can be evaluated explicitly for the integrable initial states using the analytic results of Ref. PVCP18_I for the corresponding post-quench quasiparticle distribution functions. We report in Fig. 1 the final outcome of this analysis, from which we see that the numerical results are in very good agreement with the analytic predictions.

Although the fitting procedure to extrapolate the infinite-time limit bears a numerical error, the method is accurate enough to clearly resolve the entanglement evolution for different initial states, as shown in Fig. 2$(a)$. We have also tested our numerics against the contribution $S_{l}^{(1)}(t)$ of the first quasiparticle species in Eq. (16): this is displayed in Fig. 2$(b)$, where we see that the numerical data are in agreement with the prediction (16), which takes into account the contribution of both species, but not with the individual term $S_{l}^{(1)}$. This rules out the possibility of a picture where only the first species of quasiparticles contributes to the entanglement dynamics. Overall, our analysis fully corroborates the validity of Eq. (16) as the correct generalization of the results of Ref. AlCa17 ; AlCa18 to the case of nested spin chains.

The knowledge of the post-quench rapidity distribution functions corresponding to the integrable initial states also allows us to provide predictions for other quantities related to the entanglement entropy, such as the mutual information between two subsystems $A$ and $B$, which has been recently proposed as a simple measure of scrambling asplund ; AlCa19 . We recall that the mutual information is defined as $I_{A:B}=S_{A}+S_{B}-S_{A\cup B}$, where $S_{A}$, $S_{B}$, and $S_{A\cup B}$ are the von Neumann entanglement entropies of the subsystems $A$, $B$ and of their union, $A\cup B$, respectively. By employing the very same quasiparticle picture discussed above, one can derive the following formula MBPC17

[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. Here we assumed $A$ and $B$ to be of the same size $l$, and we denoted their distance by $d$, while $s^{(r)}_{n}(\lambda)$ and $v^{(r)}_{n}(\lambda)$ are defined as in Eq. (17). We stress that Eq. (19) is conjectured to be exact in the limit $l,d,t\to\infty$, while keeping the ratios $l/t$ and $d/t$ constant. The mutual information is interesting as its time dependence signals, at large times, the presence of different species and bound states of quasiparticles, manifesting themselves as subsequent peaks of $I_{A:B}(t)$ MBPC17 ; PVCP18_I . Unfortunately, the entanglement entropy of multiple disjoint intervals in tDMRG simulations is computationally more demanding, when compared to the entropy of a single interval. Hence, we are able to provide reliable numerical data only for small system and subsystem sizes. An example is shown in Fig. 2$(c)$, where $I_{A:B}(t)/l$ is plotted as a function of $v_{M}t/l$ for a given initial integrable state (here $v_{M}$ is the maximum velocity of the quasiparticles, as computed via Bethe ansatz). From the figure we see that $I_{A:B}(t)\simeq 0$ up to a time $t\simeq d/2v_{M}$, after which it increases linearly, as expected from Eq. (19); however, the time scales and system sizes that we can reach are too small to be compared with the theoretical prediction (19), and to detect the peaks associated with the different species and the bound states of the quasiparticles MBPC17 ; PVCP18_I .

Finally, with the same numerical procedure we have also computed the evolution of the Rényi entropies (14) for different values of $\alpha$. Since there is not yet a Bethe ansatz prediction for their post-quench dynamics AlCal-b , the corresponding numerical results are reported in Appendix A, to which we refer the interested reader.

IV Evolution of correlation functions

In this section, we study another important aspect of the post-quench dynamics, namely the evolution of local correlation functions. It is useful to start by observing a few properties of the Hamiltonian (1) and of the integrable initial states (11). Focusing on the site $k=L/2$, we notice that the Hamiltonian is invariant under the “bond-inversion symmetry” $P_{B}$ around the “link” between $k$ and $k+1$, acting on local operators as $O_{j}\to O_{L+1-j}$, where $O_{j}$ acts as the identity on $h_{n}$ for $n\neq j$. The same transformation leaves invariant an integrable state $\ket{\Psi_{0}}$ [with two-site block (13)], namely $P_{B}\ket{\Psi_{0}}=\ket{\Psi_{0}}$. Finally, using that $P_{B}O_{L/2}P_{B}^{-1}=O_{L/2+1}$, this implies

[TABLE]

In the case of periodic boundary conditions, one could repeat the same argument for each site in the chain, obtaining

[TABLE]

We stress that (21) holds despite $\ket{\Psi_{0}}$ is not translationally invariant. As a result, due to the $SU(3)$ symmetry of the Hamiltonian (1), the expectation value of point-wise operators $O_{j}$ is conserved after the quench from integrable initial states, as it is easily shown using (21).

According to the discussion above, the simplest non-trivial correlation functions that we can consider are two-site correlators. We start by studying $\braket{S^{z}_{j}(t)S^{z}_{j+1}(t)}$ and $\braket{S^{x}_{j}(t)S^{x}_{j+1}(t)}$, which we have computed numerically for a quench from different initial states; an example is reported in Fig. 3. In each case, the correlators depend on the parity of $j$, due to the two-site structure of the initial state (which breaks translational invariance). From Fig. 3 we see that within the available time scales the two-site correlation functions exhibit large oscillations; furthermore, these appear to be around the same value for even and odd $j$, suggesting the expected restoration of translational symmetry at large times FCEC14 . We have found that oscillations appear to be suppressed in the evolution of four-point functions, as shown in Fig. 3$(c)$. In any case, our numerical data do not allow us to extract the asymptotic value of the local correlation functions. Furthermore, while exact result have been recently derived in the thermal case to compute expectation values of local observables in the Lai-Sutherland model RiKl19 , such results have not yet been generalized to the case of GGE-states, as in the case of $XXZ$ Heisenberg chains MePo14 ; Pozs17 . For these reasons, it appears problematic to test the validity of the analytic results of Ref. PVCP18_I ; PVCP18_II for the post-quench rapidity distribution functions, based on the asymptotics of local correlations.

Luckily, one can perform another test, albeit more indirect, which can be carried out within the time scales accessible to our numerical scheme. This is based on the computation of the “light-cone” velocity of correlations BoEL14 ; FSEH13 , namely the speed with which correlations spread after a quench. First, it is useful to recall that the Bethe ansatz provides us with a prediction for the latter in the scaling limit of large distances and times, based on the knowledge of the post-quench rapidity distribution functions. Let us consider $\braket{O_{j}(t)O_{j+l}(t)}$, where $O_{x}$ is an observable localized at site $x$. According to the quasiparticle picture, correlations between $O_{j}$ and $O_{j+l}$ are caused by pairs of quasiparticles produced in the same spatial point and that arrive at positions $j$ and $j+l$ at the same time. From this picture, we obtain that the time needed for these local operators to be correlated is at least $t=l/(2v_{M})$, where $v_{M}$ is the maximum velocity of the quasiparticles. This can be computed by Eq. (10), once the functions $\rho^{(r)}_{n}(\lambda)$ are known, so that the numerical computation of time-evolved correlators provides a test for the validity of the results of Refs. PVCP18_I ; PVCP18_II for the post-quench rapidity distribution functions.

In order to probe the light-cone velocity of the correlation functions, we compute $\braket{S^{\alpha}_{j}(t)S^{\alpha}_{j+l}(t)}$ for $\alpha=x,z$ and increasing values of $l$. Examples of our results for a particular integrable initial state are reported in the plots of Fig. 4, from which one can immediately see the emergence of a light-cone structure. Quantitative results for the corresponding velocity can be obtained from these plots following the procedure employed in Ref. BoEL14 ; MWNM09 : in particular, for a given correlator $\braket{S^{\alpha}_{j}(t)S^{\alpha}_{j+l}(t)}$ one can identify the first “inflection point” [i. e. the first maximum of the its time-derivative, cf. Fig 5$(c)$] as the arrival time $t_{l}$ of the light-cone and plot the sequence $t_{l}$ as a function of $l$. Then, the velocity of the light-cone can be extracted as the slope of the straight line fitting the set of points $(t_{l},l)$ (which should be independent of the particular correlation function which has been chosen, as we have explicitly verified).

A summary of our results is displayed in Fig. 5. First, in subfigure $(a)$ we report the Bethe ansatz predictions for the maximum quasiparticle velocity corresponding to different initial states. We see that the latter has a weak dependence on the initial parameters $c_{jj}$. Furthermore, states with large velocity have large entanglement and are more difficult to simulate. Hence, the expected maximum velocity of the states that we can study numerically vary in a rather restricted range. In Fig. 5 we display our results for the light-cone velocity for three different initial states, from which the aforementioned weak dependence on the state parameters is apparent. The numerical data are compared against the ansatz $y(t)=2v_{M}t+c$, where we have performed a fit to determine the constant $c$. We see that they are in good quantitative agreement with the analytic predictions, and that these are able to resolve the different initial states, despite the weak dependence on the quench parameters. In conclusions, our results are always seen to corroborate the validity of the Bethe ansatz predictions for the light-cone velocity of local correlations.

V Conclusions

In this work we have presented numerical results for quantum quenches in the nested Lai-Sutherland model from a family of integrable initial states, focusing on the spreading of both correlation functions and entanglement entropy. By means of tDMRG simulations, we tested explicitly the validity of the nested version of the conjectured quasiparticle formula for the growth of the entanglement entropy of Refs. AlCa17 ; AlCa18 , and the Bethe ansatz predictions for the “light-cone” velocity of correlation functions BoEL14 . The present paper complements other studies in the literature, extending to the case of nested models the body of numerical evidence corroborating the existing theory of integrable systems out of equilibrium CaEM16 .

The theoretical predictions tested in this work are based on the analytic results of Ref. PVCP18_I ; PVCP18_II , where the post-quench rapidity distribution functions for the integrable states were computed. While the small initial entanglement of these states allowed us to simulate their dynamics by means of the tDMRG algorithm, we have seen that the correlation functions exhibit large relaxation times, so that it is not possible to extract their asymptotic behavior based on short and intermediate-time simulations. In this respect, it would be highly desirable to generalize the string-charge duality method IQDB16 to nested systems: the latter would allow us to enlarge the class of initial states for which analytic predictions can be made, hopefully including states with faster local relaxation. Another open problem with strong connections to our work is the derivation of analytic formulas to compute GGE correlation functions in nested models, analogously to what has been done in the case of Heisenberg chains MePo14 . These would make it possible to give predictions for the asymptotics of local observables based on the knowledge of the post-quench rapidity distribution functions and allow for a direct comparison with numerical results. We hope that our work will serve as a motivation for future investigations in these directions.

Acknowledgments

L. P. acknowledges support from the Alexander von Humboldt foundation and from the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy – EXC-2111 – 390814868. P. C. acknowledges support from ERC under Consolidator grant number 771536 (NEMO). P. C. and L. P. acknowledge Eric Vernier and Balász Pozsgay for discussions and collaboration on related subjects.

Appendix A Details on the numerical computations

In this appendix we provide details on the numerical methods that we have employed for the computation of the von Neumann and Rényi entanglement entropies, and present additional plots for the latter.

All the data displayed in the main text are obtained using the tDMRG algorithm, as implemented in the iTensor library itensor , where open boundary conditions have always been chosen. The data in Figs. 1 and 2 are obtained by considering a chain of $L=48$ sites, and the subsystem is always chosen at one of its edges. Due to open boundary conditions, the numerical data are compared with the Bethe ansatz formula in Eq. (16) with the substitution $t\to t/2$, which corresponds to the fact that the quasiparticles can only enter from one of the two borders of the subsystem under consideration AlCa17 . We choose the maximum bond dimension $\xi=1500$ and a time step $\Delta t=0.02$ for all simulations. We have checked the robustness of our results by performing simulations for different values of $\Delta t$, $\xi$ and $L$. All the entanglement entropy data reported in the manuscript are for $L=48$, which is the maximum system size we could reach with our computational resources. Furthermore, all of our results have been checked against exact diagonalization calculations for small times.

For the initial states in Fig. 1, the entanglement entropy $S_{l}(t=0)$ is different, depending on whether the subsystem considered contains an odd or an even number of sites. More precisely, it is a non-zero constant for $l$ odd, while it is vanishing for $l$ even. Of course, in both cases one has $\lim_{l\to\infty}S_{l}(t=0)/l=0$. In order to extrapolate the value of $S_{l}/l$ to the thermodynamic limit, we have performed a fit of the data (at fixed $v_{M}t/l$) using the function

[TABLE]

where $s_{\infty}$, $a$ and $b$ are fitting parameters. Given the even/odd effect in the finite-size data, we have performed the fit independently for the even and odd sequences, obtaining different numerical results $s^{(e/o)}_{\infty}$. The (small) differences between the fits for even and odd values of $l$ are attributed to finite size effects; hence, as the best estimate for $s_{\infty}$ we chose the average $(s^{(e)}_{\infty}+s^{(o)}_{\infty})/2$. We note that the numerical difference between $s^{(e)}_{\infty}$ and $s^{(o)}_{\infty}$ was found to decrease as $t/l$ increased for our set of finite-size data, consistently with the expected restoration of translational symmetry at large times.

As we have mentioned in Sec. III, we used the same extrapolation procedure to compute the dynamics of the Rényi entropy of order $\alpha$ in the scaling limit $t,l\to\infty$. An example of our result, corresponding to $\alpha=1/2$, is displayed in Fig. 6. From the figure, we see that, as expected, the qualitative features of the Rényi and von Neumann entropies are similar, at least at the time scales available to our numerical procedure. We chose to report the Rényi entropy of order $1/2$ because it is related to the time evolution of the entanglement negativity ctc-14 ; AlCa19 , a relevant entanglement measure for tripartite systems.