Entanglement dynamics after quantum quenches in generic integrable systems

The time evolution of the entanglement entropy in non-equilibrium quantum systems provides crucial information about the structure of the time-dependent state. For quantum quench protocols, by combining a quasiparticle picture for the entanglement spreading with the exact knowledge of the stationary state provided by Bethe ansatz, it is possible to obtain an exact and analytic description of the evolution of the entanglement entropy. Here we discuss the application of these ideas to several integrable models. First we show that for non-interacting systems, both bosonic and fermionic, the exact time-dependence of the entanglement entropy can be derived by elementary techniques and without solving the dynamics. We then provide exact results for interacting spin chains that are carefully tested against numerical simulations. Finally, we apply this method to integrable one-dimensional Bose gases (Lieb-Liniger model) both in the attractive and repulsive regimes. We highlight a peculiar behaviour of the entanglement entropy due to the absence of a maximum velocity of excitations.


Introduction
In recent years, understanding the entanglement structure of out-of-equilibrium many-body quantum systems has become an emerging research theme at the crossroad between statistical physics, condensed matter physics, quantum field theory, and quantum information. In one dimension, the growth of entanglement has been related to the capability of a classical computer to simulate non-equilibrium quantum systems with matrix product states (see, e.g., [1][2][3][4][5]). Moreover, the thermodynamic entropy in a stationary state has been interpreted as the asymptotic entanglement of a large subsystem [6][7][8][9][10].
One of the prototype protocols for driving a system out-of-equilibrium is the quantum quench [11][12][13][14][15][16][17][18]: An isolated system is initially prepared at t = 0 in a given pure state |ψ 0 (usually the ground state of a quantum many-body hamiltonian H 0 ) and for t > 0 the unitary dynamics is governed by a hamiltonian H (with [H, H 0 ] = 0 e.g., at t = 0 a parameter of the hamiltonian is suddenly changed). Besides the theoretical interest, in recent years it has become possible to investigate quantum quenches experimentally with cold-atom systems [19][20][21][22][23][24][25][26][27][28][29][30][31]. Since the post-quench dynamics is unitary, the full system never reaches stationary behaviour, which, instead, can arise locally. The central object to define local equilibration is the reduced density matrix. Given a subsystem A of the full system, the reduced density matrix ρ A is defined as where the trace is over the degrees of freedom of the complement B of the subsystem A, and |ψ ≡ e −iHt |ψ 0 is the time-dependent state of the system. For quantum quenches in generic models, the stationary behaviour of local and quasilocal observables is described by the Gibbs (thermal) ensemble [32][33][34][35][36][37][38]. In contrast integrable models possess an extensive number of conserved quantities, besides the hamiltonian, which highly constrain the post-quench dynamics. As a consequence, integrable systems fail to thermalise, meaning that the reduced density matrix for long times is not thermal. Remarkably, a statistical description of local properties of the steady state is possible in terms of a Generalised Gibbs Ensemble (GGE) [12,16,17,, which is obtained by complementing the Gibbs ensemble with all the local and quasilocal conserved quantities [56,67].
The problem of understanding how entanglement spreads after a quench is deeply intertwined with that of equilibration and thermalisation. The standard measure of the entanglement is the entanglement entropy [68] which is defined as the von Neumann entropy of the reduced density matrix (1): The out-of-equilibrium dynamics of the entanglement entropy following quantum quenches has been the focus of intense research during the last decade [6,10,. Remarkably, in recent years it has become possible to measure entanglement and its evolution in cold-atom experiments [9,88,93]. For a wide variety of global quenches, the quasiparticle picture of Ref. [6] provides an understanding of the main qualitative features of the entanglement dynamics. In the quasiparticle picture, the pre-quench initial state is a source of pairs of excitations with opposite momentum that travel ballistically through the system. Let us assume that there is only one type of excitations (quasiparticles) identified by their quasi-momentum λ, and moving with group velocity v(λ). The main assumption of the quasiparticle picture is that excitations that are created far apart from each other are incoherent, whereas those emitted at the same point in space are entangled (more precisely, quasi-particles emitted within the initial correlation length, but this refinement just provides a subleading correction to the result [91] and will be ignored in what follows). As the quasiparticles propagate, larger regions of the system get entangled. At time t the entanglement entropy of a subsystem A is proportional to the total number of quasiparticles that after being emitted from the same point in space are shared between subsystem A and its complement. Specifically, for an interval A of length embedded in an infinite one-dimensional system, by counting the quasiparticles with a given weight s(λ), one obtains [6] S(t) = 2t 2|v(λ)|t< dλv(λ)s(λ) + 2|v(λ)|t> dλs(λ). ( Here the function s(λ) depends on the production rate of quasiparticles with quasimomentum ±λ and on their individual contribution to the entanglement entropy. Formula (3) holds true in the space-time scaling limit t, → ∞ with the ratio t/ fixed. Notice that (3) does not take into account O(1) terms, which are subleading in the scaling limit. When a maximum quasiparticle velocity v M exists, such that |v(λ)| ≤ v M (e.g., as a consequence of the Lieb-Robinson bound [92]), Eq. (3) predicts that for t ≤ /(2v M ), S grows linearly in time because the second term in (3) vanishes. In contrast, for t /(2v M ), only the second term contributes and the entanglement is extensive in the subsystem size, i.e., S ∝ . Eq. (3) describes the light-cone spreading of the entanglement evolution which has been analytically confirmed in few free models [69][70][71][72][73] and also verified in several numerical studies (see e.g. [77][78][79][80][81]). However, in order to give some predictive power to (3), we should have a way to fix the entropy density s(λ) and the velocity of the entangling quasiparticles v(λ). Yet, determining s(λ) ab-initio from the dynamical problem is a formidable task even for free models (see e.g. [69]); furthermore, for interacting integrable models, also the identification of the velocity v(λ) is a non-trivial issue. A major breakthrough in this respect has been achieved in [10] where it has been shown that, at least for certain classes of quenches in integrable models, the function s(λ) can be conjectured from the equivalence between the entanglement and the thermodynamic entropy in the stationary state. The latter can be straightforwardly calculated with equilibrium techniques from the GGE describing the stationary state. In this way Eq. (3) becomes a quantitative analytic conjecture for the entanglement evolution which can be obtained only from the stationary state without solving the many-body dynamics. Suggestively, we can state that the main idea of Ref. [10] is to reconstruct the entanglement evolution going back in time from the stationary state. Physically, Eq. (3) highlights the transformation during the dynamics of the entanglement into the thermodynamic entropy. This transformation happens for non-integrable systems as well, but in that case, the entanglement entropy becomes the thermal entropy [7,9,94].
In a generic interacting integrable model there are several families of quasiparticles. The generalization of (3) is obtained by summing all the contributions of the different species. The final result of Ref. [10] for the entanglement dynamics is where the index n labels the different families of elementary quasiparticles present in a generic integrable model, and λ is their associated momentum label. The sum over the quasiparticle families and momenta reflects the presence in integrable models of well-defined excitations with an infinite lifetime. According to the ideas of Ref. [10] outlined above, in Eq. (4), s n (λ) can be conjectured from the contribution of the individual quasiparticles to the thermodynamic entropy of the GGE describing the steady state. Furthermore, the velocities v n (λ) are assumed to be the group velocities of the low-lying excitations around the steady state. The validity of (4) has been checked numerically for several quenches in the Heisenberg XXZ chain [10]. A generalisation of (4) has been provided to describe the entanglement evolution after inhomogeneous quenches in the XXZ chain [95].
In this work we discuss in detail several applications of (4). We start focusing on free fermionic and free bosonic models for which we provide generic results valid for a wide class of quenches. We show that it is possible to recover, in an elementary manner, the known result for the entanglement dynamics after a generic quench in the transverse field Ising chain [69]. For the bosonic case, the quasiparticle picture provides new exact results for the entanglement dynamics in the harmonic chain (the lattice discretisation of the one-dimensional Klein-Gordon field theory). This result is remarkable also because its ab initio derivation is not available yet, although we are dealing with a free model. Then, we turn to discuss the entanglement dynamics in the anisotropic Heisenberg chain (XXZ chain). We provide several new theoretical predictions, which complement the results already presented in [10]. For instance, we provide exact results for the post-quench dynamics of the mutual information between two intervals starting from several initial states. This is important because the mutual information is a useful tool to probe the validity of the quasiparticle picture, in which welldefined quasiparticles entangle different regions of the system, so that the mutual information exhibits a peak at intermediate times. An alternative picture is the information scrambling scenario [96][97][98][99][100], which should apply to many non-integrable models such as irrational 1+1 conformal field theories. In the scrambling scenario the quasiparticles loose coherence during the dynamics, due to scattering. As a consequence, for large time the mutual information vanishes independently of the separation of the intervals. Conversely, in integrable models, well-defined quasiparticles exist, and the the mutual information in the space-time scaling limit has a peak also at large times for large enough separation of the intervals, ruling out the scrambling scenario. Numerical evidence supporting the validity of the quasiparticle picture for the mutual information has been provided in [10] considering the quench from the Néel state in the XXZ chain. Moreover, in this work we investigate the signatures of composite excitations (multi-particle bound states) in the mutual information dynamics. An interesting result is that the presence of bound states leads to an anomalous decay of the mutual information at late times and, for some quenches, to multi-peak structures (as already highlighted for other models in [101]).
Another main result obtained here is a quasiparticle prediction for the entanglement dynamics in the one-dimensional Bose gas. We focus on the quench from the Bose-Einstein condensate (BEC), considering both the attractive and repulsive Lieb-Liniger model. In both cases, at short times the von Neumann entropy exhibits a non-linear increase with time due to the fact there is no maximum velocity of propagation of excitations. Nevertheless, at long times the entanglement entropy saturates. An important difference between attractive and repulsive interactions, is that while for repulsive interactions only one species of quasiparticles is present, for attractive ones multi-boson bound states appear. Interestingly, for weak interactions, bound states contribute significantly to the entanglement dynamics. Moreover, similar to the XXZ chain, their presence is reflected in a slow vanishing behaviour of the mutual information between two intervals at late times.
The outline of the paper is as follows. Section 2 is devoted to the entanglement dynamics after quantum quenches in free-fermion and free-boson models. In section 3, we detail the approach of [10] for the entanglement dynamics in a generic Bethe ansatz integrable model. In section 4 we provide several results for the entanglement dynamics in the XXZ chain. In section 5 we present the quasiparticle results for the entanglement dynamics after the quench from the Bose-Einstein condensate in the Lieb-Liniger model. In the last section we discuss several points and developments which deserve further investigation.

Entanglement dynamics in free models
In this section we employ the quasiparticle scenario of [10] to derive analytically the entanglement dynamics in free-fermion and free-boson models after rather generic quenches. We test these results against exact analytical and numerical results for the entanglement dynamics after a global quench in the transverse field Ising/XY chain and in the harmonic chain. These models can be mapped onto a system of free fermions and free bosons, respectively. For the Ising model our result agrees with the ab initio derivation in [69], providing a further benchmark of the ideas pursued in this paper and in [10]. For the harmonic chain our results have been anticipated in [102] and appeared, for a similar bosonic model, also in [103].

Models of free fermions
If a translational invariant fermionic model is free, it means that the hamiltonian in momentum space can be mapped into (apart from an unimportant additive constant) where b k are fermionic mode occupation operators satisfying standard anticommutation relations and k is the energy of the mode k (i.e. the dispersion relation). For all these models, the GGE built with local conservation laws is equivalent to the one built with the mode occupation numbersn k = b † k b k since they are linearly related [46,50]. Thus the local properties of the stationary state are captured by the GGE density matrix where Z = Tre − k λ knk ensures the normalisation Trρ GGE = 1. The thermodynamic entropy of the GGE is obtained by elementary methods, leading, in the thermodynamic limit, to where n k ≡ n k GGE = Tr(ρ GGEnk ) and the function H is The interpretation of Eq. (7) is obvious: the mode k is occupied with probability n k and empty with probability 1 − n k . Given thatn k is an integral of motion, one does not need to compute explicitly the GGE (6), but it is sufficient to calculate the expectation values ofn k in the initial state ψ 0 |n k |ψ 0 which, by construction, equals n k = n k GGE . At this point, following [10], we identify the stationary thermodynamic entropy with the density of entanglement entropy to be plugged in Eq. (3), obtaining the general prediction where k = d k /dk is the group velocity of the mode k. This formula is generically valid for arbitrary models of free fermions with the crucial but rather general assumption that the initial state is writable in terms of pairs of quasiparticles. More general and peculiar structures of initial states can be also considered, see [104,105].

Test for the transverse field Ising chain
Eq. (9) can be tested against available exact analytic results for the transverse field Ising chain with hamiltonian where σ x,z j are Pauli matrices and h is the transverse magnetic field. We use periodic boundary conditions in (10).
The hamiltonian (10) is diagonalised by a combination of Jordan-Wigner and Bogoliubov transformations, leading to Eq. (5) where the single-particle energies are We focus on a quench of the magnetic field in which the chain is initially prepared in the ground state of (10) with h 0 and then, at t = 0 the magnetic field is suddenly changed from h 0 to h. As in the general analysis above, the steady-state is determined by the fermionic occupation numbers n k given by [46,106] where ∆ k is the difference of the pre-and post-quench Bogoliubov angles [106] where 0 (p) and (p) stand for pre-and post-quench dispersion relations respectively. The quasiparticle prediction for the entanglement dynamics after the quench is then Eq. (9) with n k in (12). This coincides with the ab initio derivation performed in [69]. The same derivation is valid also for a generic quench in the XY chain reported in [69].

Free bosonic models
For a free bosonic model, the hamiltonian can be written after some suitable transformations as (apart from a unimportant additive constant) with [a k , a † k ] = δ k,k being bosonic mode operators. The stationary values of local observables can be described by a generalised Gibbs ensemble (GGE) constructed from the mode occupation numbersn k = a † k a k with the GGE density matrix where λ k are Lagrange multipliers and Z is a normalisation. The Lagrange multipliers λ k in (15) are fixed by imposing that the expectation value ofn k in the initial state coincides with its GGE average. The initial value n k ≡ ψ 0 |n k |ψ 0 is easily calculated from the initial state. The GGE expectation value ofn k is obtained as with Thus, one has After imposing the conservation ofn k , i.e., that (18) equals n k = ψ 0 |n k |ψ 0 , one obtains λ k as At this point, calculating the thermodynamic entropy is a trivial exercise in statistical physics: Using that Z = k (1 − e −λ k ) −1 (cf. Eq. (17)), we obtain where we used n k = 1/(e λ k − 1), cf. (19). In the thermodynamic limit the sum over the momenta becomes an integral and (22) becomes where in the rightmost side of the equation we introduced the entropy contribution s(k) of the quasiparticle with momentum k as 2πs(k) = (n k + 1) ln(n k + 1) − n k ln n k .
At this point, we are ready to use the fact that the entanglement entropy is the stationary thermodynamic entropy (23) so that the quasiparticle picture for the entanglement evolution (4) gives where the entropy density s(k) is given by (24) and v k = d k /dk.

Tests for the harmonic chain
Here we focus on one of the simplest bosonic models with an exactly solvable non-equilibrium dynamics, i.e., the harmonic chain defined by the hamiltonian with periodic boundary conditions. Eq. (26) defines a chain of N harmonic oscillators with frequency (mass) m and with nearest-neighbour quadratic interactions. Here φ n and π n are the position and the momentum operators of the n-th oscillator, with equal time commutation relations [π m , π n ] = iδ nm [φ n , φ m ] = [π n , π m ] = 0 .
In the context of quench dynamics the harmonic chain was first discussed in [11] to which we refer for a detailed analysis; here we only report the results relevant for our aims. The harmonic chain is easily diagonalised in momentum space where it assumes the standard diagonal form (14) with disperion relation We now consider the quantum quench in which the harmonic chain is initially prepared in the ground-state |ψ 0 of (26) with m = m 0 , and at time t = 0 the mass is quenched to a different value m = m 0 . We use the notation 0 k for the dispersion relation in the initial state and k for the one for t > 0.
In order to give predictive power to Eq. (25) we just need to fix the conserved value n k = ψ 0 |n k |ψ 0 which is obtained by elementary methods [12] Also the group velocity from (28) is The quasiparticle prediction (cf. (25)) for the entanglement dynamics after the mass quench in the harmonic chain is reported in Figure 1 The entanglement entropy exhibits the expected linear behaviour at short times followed by a saturation at asymptotically long times. Clearly, the steady-state value of the entanglement entropy increases with m. In the limit m m 0 , the steady-state entropy at the leading order in 1/m is S ≈ ln m. The crossover time from the linear to the saturation regime increases with m, because the maximum velocity v M decreases upon increasing m, as it is clear from (30).

Numerical checks
We now provide numerical checks of the validity of (25). The entanglement dynamics after a global quench in the harmonic chain has been studied numerically in several papers [71,84,85]. These papers focused on the critical (m → 0) and continuum limit, in which several simplifications occur because there is a single velocity of excitations. The quasi-particle prediction turned out to be correct, but with additive logarithmic corrections due to the presence of a zero mode [85]. In the following we focus on the massive regime that so far received only little attention. For systems of free bosons, at any time after the quench the entanglement entropy of a finite subsystem can be calculated effectively [107,108] from the time-dependent twopoint correlation functions reported in [12]. In Figure 2 we present numerical results for the entanglement entropy S(t) after a mass quench in the harmonic chain. The results are for a chain with L = 1000 sites and subsystems sizes = 10, 20, 100. We numerically checked that for these values of the effect of the finite L is negligible. The two panels (a) and (b) show results for the quenches with m = 2 and m = 3, respectively. The pre-quench value of the mass m 0 = 1 is the same for both quenches. The theoretical prediction obtained using the quasiparticle picture (cf. (25)) is reported in the Figure as dashed-dotted line. For any finite scaling corrections are expected because Eq. (25) holds only in the space-time scaling limit with , t → ∞, at t/ fixed. These corrections are clearly visible in the data. However, they rapidly decrease upon increasing , and the results for = 100 are almost indistinguishable from the thermodynamic limit predictions.

Entanglement dynamics in a generic Bethe ansatz integrable model
In this section, following the ideas of [10], we show how the quasiparticle prediction (4) can be applied to a generic Bethe ansatz integrable model. In order to do so, in the next two subsections we provide explicit conjectures for the values of s n (λ) and v n (λ) to be plugged in (4). As explained in the introduction, s n (λ) can be read off from the thermodynamic entropy in the stationary state that can be worked out in the thermodynamic Bethe ansatz framework. For v n (λ), we will instead use the velocity of low-lying particle-hole excitations built on top of the stationary state. In the following subsections, we will show how to derive these velocities by Bethe ansatz techniques following Ref. [109].

The thermodynamic Bethe ansatz
In a Bethe anstaz integrable model of length L, with N elementary particles, and with periodic boundary conditions, the eigenstates are in one to one correspondence with a set of N complex quasi-momenta λ j (known as rapidities) which satisfy model dependent quantisation condition denoted as Bethe equation. (Here we focus on models with an "elementary" Bethe ansatz; there are models with more than one type of rapidities leading to the so-called nested Bethe ansatz [110]; in that case the modification of (4) is straightforward because one has just to perform a further sum on the types of the rapidities, see [101] for an illustrative example.) The prototype integrable model that we consider here is the XXZ spin-1/2 chain in the regime with ∆ > 1, although the TBA results that we will discuss can be generalized to the case with ∆ < 1 and to other integrable models with minor modifications. In the thermodynamic limit and for a generic translational invariant model, the vast majority of the solutions of the Bethe equations obey the string hypothesis [111]. Specifically, solutions of the Bethe equations form string patterns in the complex plane. Rapidities forming a n-string are parametrised as [111] λ j n,γ = λ n,γ + i where η is an interaction parameter, j = 1, . . . , n labels the different string components, λ n,γ is the "string centre", and δ j n,γ are the string deviations, which for the majority of the eigenstates are δ j n,γ = O(e −L ), implying that they can be neglected in the thermodynamic limit (string hypothesis [111]). Physically, a n-string corresponds to a bound state of n elementary particles. For the XXZ chain with ∆ < 1 the structure of the string solutions is more complicated [111] than (31), although major simplifications occur for ∆ = ∆ k ≡ − cos(π/k) with k = 1, 2, . . . (roots of unity).
Within the framework of the string hypothesis, the string centres λ n,γ are obtained by solving the Bethe-Gaudin-Takahashi (BGT) equations [111] Lπ n (λ n,α ) = 2πI n,α + (n,α) =(m,β) Here I n,α are (integer or half-integer) quantum numbers, π n (x) are model dependent functions for the string momentum. The scattering phases for the bound states Θ n,m (λ) can be written as in terms of a model dependent elementary phase shift θ n (λ). Each different choice of I n,α identifies a different set of solutions of (32), which correspond to a different eigenstate of the considered integrable model. The corresponding eigenstate energy E and total momentum P are obtained by summing over all the BGT rapidities [111] as where n (λ) is the model dependent string energy, while z n (λ n,α ) = 2πI n,α /L so that the total momentum P depends only on the I n,α . In the thermodynamic limit the solutions of the BGT equations (32) become dense on the real axis. The central quantities to describe local properties of the system are then the rapidity densities ρ n (λ) (n labelling different string types) which are formally defined in the thermodynamic limit as To fully specify the thermodynamic state of the system, the densities ρ (h) n (λ) of the n-string holes, i.e., of the unoccupied string centres are also required. Finally, it is also custom [111] to introduce the total densities ρ (t) n (λ). Some TBA relations are written in a more compact form in terms of the ratio that we introduce for future convenience. The ρ (h) n (λ) and ρ n (λ) are obtained via the thermodynamic version of the BGT equations which are obtained from (32) by taking the thermodynamic limit. The symbol f g denotes the convolution between two functions as The functions b n (λ) and a nm (λ) are related to π n (λ) and Θ nm (λ) as In the thermodynamic limit, the expectation values of local conserved quantities are functionals of the densities ρ n (λ); for example the particle and energy densities are The set of rapidity densities ρ ρ ρ ≡ {ρ n } ∞ n=1 defines a thermodynamic macrostate, which encodes all the expectation values of local or quasi-local observables in the thermodynamic limit. A generic thermodynamic macrostate corresponds to an exponentially large (with L) number of microscopic eigenstates of the model, all leading to the same set of rapidity densities in the thermodynamic limit. The total number of possible The Yang-Yang entropy represents the thermodynamic entropy of a given macrostate, as it should be clear from a generalised microcanonical argument. For example, it has been proved that for systems in thermal equilibrium S Y Y coincides with the thermal entropy [111]. Our conjecture for the time evolution of the entanglement starts from the Yang-Yang entropy since we assume that at long times the entanglement entropy is the thermodynamic one. Furthermore, we also assume that the Bethe quasiparticles are the one entangling the system and appearing in (4). Thus it is natural to identify s n (λ) with the integrand in (42), i.e.
Here the three sets of root densities ρ n , ρ (h) n , and ρ (t) n refer to the macrostate that describes the stationary state. This is in principle calculable by Bethe ansatz techniques from the overlaps of the initial state with the Bethe states [113,114] or equivalently from the GGE [56].

Group velocities over a macrostate
Having identified s n (λ) in Eq. (4), the other crucial ingredient for the quasiparticle picture for the entanglement dynamics is the group velocity of the entangling quasiparticles. In the approach of [10] the entangling quasiparticles are identified with the low-lying excitations (particle-hole excitations) around the thermodynamic macrostate describing the steady state.
The low-lying excitations over a given macrostate can be constructed explicitly in the framework of TBA as originally pointed out for the stationary state after a quench in [109] and only briefly summarised in the following. The first step is to choose, among the equivalent eigenstates of the macrostate identified by the densities ρ n , ρ (h) n , one representative microstate at finite, but large, volume L. This corresponds to a particular set of BGT quantum numbers I n,α in (32) chosen in such a way that the resulting rapidities from the BGT equations are a discretisation of the desired macrostate. A particle-hole excitation in each n-string sector is obtained by replacing I n,h → I n,p , where I n,p (I n,h ) is the BGT number of the new added particle (hole). Due to interactions, this local change in quantum numbers implies a rearrangement of all the rapidities. The excess energy of the particle-hole excitation is easily calculated as δE n = e n (λ n,p ) − e n (λ n,h ).
Remarkably, apart from the dressing of the "single-particle" energy e(λ) (44) is the same as for free models. Similarly, the change in the total momentum is obtained from (54) as Finally, the group velocity of the particle-hole excitations is by definition .
Here we used that dz n (λ)/dλ = 2πρ The function e n (λ) is determined by solving an infinite system of Fredholm integral equations of the second kind as Equations (47) are routinely solved numerically by truncating the system, i.e., considering n ≤ n max and checking convergence with varying n max . The method outlined above for calculating the group velocities has been introduced in [109] in order to study velocity of the spreading of correlation after a quench from a thermal state. Very recently it has also been used to study transport properties in integrable models [115][116][117][118]. At this point, we have a Bethe ansatz procedure to calculate the velocities of the entangling quasiparticles and we are ready to use the conjecture (4) to provide quantitative predictions for the entanglement spreading in generic integrable systems.

Entanglement dynamics in Heisenberg spin chains
In this section we focus on the spin-1/2 anisotropic Heisenberg chain (XXZ chain). The goal of this section is to provide a thorough discussion of some results that have been already presented in [10] and to extend them in several directions.
The XXZ chain is defined by the hamiltonian Here S α i are spin-1/2 operators acting at site i of the chain, and ∆ is the anisotropy parameter. Periodic boundary conditions are used in (48).

Bethe ansatz solution of the XXZ chain
In the Bethe ansatz solution of the XXZ chain, the eigenstates of (48) can be labeled by the total number of down spins (particles). Eigenstates in the sector with M particles are in correspondence with a set of M rapidities λ j . The rapidities are obtained by solving a set of non linear algebraic equations (Bethe equations) as [111] sin where η ≡ arccosh(∆). In the thermodynamic limit the vast majority of the solutions of the Bethe equations (52) organise according to the string hypothesis (31). For the XXZ spin-chain, physically, a n-string corresponds to a bound states of n down spins. The BGT equations are given in (32) in which one should identify θ n (λ) = π n (λ) = 2 arctan tan(λ) tanh(nη/2) .
For ∆ > 1, the string centres are in the interval [−π/2, π/2). The eigenstate energy E and total momentum P are given by Eq. (34) with string energy The thermodynamic version of the BGT equations are given by (37). We also consider the XXZ chain in the limit ∆ = 1 (XXX chain). The Bethe ansatz results for the XXX chain can be obtained from those for the XXZ chain by taking an appropriate scaling limit. The first step is to rewrite the formulas for the XXZ chain in terms of the rescaled rapidities µ defined as As η → 0, i.e., for ∆ → 1, the rescaled rapidities µ are now defined in the whole real axis [−∞, ∞]. Moreover, the spacing along the imaginary axis between different string components is i/2. Using (54) one obtains that for the XXX chain n (µ) is

Entanglement dynamics
Let us repeat here the quasiparticle prediction (4) for the entanglement dynamics where the sum is over the quasiparticle families n (strings of different length), v n (λ) is the velocity of the entangling quasiparticles numerically calculated above, and s n (λ) denotes the contribution of each quasiparticle to the Yang-Yang entropy of the steady state in Eq. (43). The exact numerical results obtained using (73) are illustrated in Figure 4. The different panels are for quenches from different initial states in the XXZ chain and several values of ∆. For the quenches from the tilted Néel and ferromagnetic states ϑ is the tilting angle. In all panels the entropy density S/ is plotted versus the rescaled time v M t/ with v M the maximum velocity, which is extracted from the Bethe ansatz. In all panels the expected behaviour with a linear increase at short times followed by an asymptotic saturation is observed. Interestingly, for all the quenches the larger steady-state entanglement is obtained for the smaller ∆. The largest amount of entanglement is produced in the quench from the tilted Néel state (panel (c)). For the Néel quench the entropy vanishes in the limit ∆ → ∞, which follows from the fact that the Néel state is the ground state of the XXZ chain in that limit, whereas it is finite for all other initial states. Finally, as already noticed in [10], for the quench from the tilted ferromagnet (see panel (b)) the linear regime seems to extend for v M t/ > 1. However, the true linear regime extends only up to v M t/ = 1. The behaviour observed in panel (b) is due to the large contributions to the entanglement entropy of slow quasiparticles (see [10]). A lot of important information is extracted by looking at the contribution to the entanglement dynamics of the individual quasiparticles. This is investigated in Figure 5 focusing on the steady-state entropy density S/ (panels (a,c)) and on the slope of the linear growth at short times (b,d) S / (we denote with entanglement production rate the quantity S (t)/ for t < /(2v M ) when it does not depend on time). Both quantities are plotted against the quasiparticles rapidity λ. All the results are for the quench in the XXZ chain with ∆ = 2. Panels (a) and (b) are for the quench from the tilted ferromagnet (with tilting angle ϑ = π/10). Remarkably, the largest contribution to the steady-state entropy and to the entanglement production rate is in the region with large λ, which correspond to slow quasiparticles (see Figure 3). Also, the largest contribution is in the sector with n = 1 (continuous line in the Figure). The contributions of higher strings are negligible (the dotted line in Figure 5 (a) (b) is the contribution of the two-particle bound states). A striking different behaviour is ob-  served for the quench from the Néel state (panels (c) and (d) in the Figure); now the largest contribution to the stationary entanglement and to the entanglement production rate is the region with small (but non-zero) rapidities, corresponding to fast quasiparticles. Similar to the quench from the tilted ferromagnet, the bound state contribution to the entanglement dynamics decays rapidly with their size. Finally, we discuss the dependence of the stationary entropy and of the entanglement production rate on the chain anisotropy ∆. This is shown in Figure 6. Clearly, for the quench from the Néel state the entropy is vanishing in the limit ∆ → ∞, as already discussed. On the other hand, it remains finite for all the other quenches. Moreover, for the quench from the Néel state and the dimer state, both the steady-state entropy and the entanglement production rates exhibit their maximum value for ∆ ≈ 1. In contrast, they vanish in the limit ∆ → 1 for the quench from the tilted ferromagnet. This is expected because at ∆ = 1 the tilted ferromagnet becomes an eigenstate of the XXZ chain for any tilting angle. The large ∆ behaviour can be understood analytically using perturbative methods. Here we discuss the behaviour of the steady-state entropy, although similar results can be derived for the entanglement production rate. It is straightforward to show that for the quench from the Néel state, in the limit ∆ → ∞, the stationary entropy is For the quench from the Majumdar-Ghosh state one has In all panels the different curves correspond bound states (strings) of different size n = 1, 2, 3. Panels (a)(b) show the results for the quench from the tilted ferromagnet with tilting angle ϑ = π/10 and chain anisotropy ∆ = 2. Notice in both cases the peaks at λ ≈ ±π/2, which signal a large contribution of the slow quasiparticle to the entanglement dynamics. Panels (c)(d) show results for the quench from the Néel state in which the largest contributions correspond to λ with small, but non-zero, value. Similar results are obtained for the quench from the dimer state and the tilted Néel. The contribution of the bound states with n > 1 are always much smaller than that for n = 1.
On the other hand, for the quenches from the tilted states the dependence of the root densities ρ n , ρ (h) n on the tilting angle is non-trivial even in the limit ∆ → ∞, implying a non-trivial dependence for the entanglement entropy as well.

Numerical checks
In this section, using tDMRG simulations [131][132][133], we provide numerical evidence supporting our main result (73). Numerical results are presented in Figure 7. Panels (a) and (b) show tDMRG simulations for the quench from the Néel state in the XXZ chain at ∆ = 1 and ∆ = 2, respectively. The results in panel (b) are the same as in [10]. Both panels plot the entropy density S/ , with the size of subsystem A, as a function of the rescaled time v M t/ , where v M is the maximum velocity calculated using the Bethe ansatz (see section 3.2). The continuous curves are tDMRG results for a chain with L = 40 sites and = 5 − 20. The dashed-dotted line is the theoretical result (73) in the scaling limit. For both ∆ = 1 and ∆ = 2 scaling corrections are visible. The diamonds are extrapolations to the thermodynamic limit. These are obtained by fitting the data at fixed v M t/ to  (76). We now turn to discuss further checks of (73) using the infinite Time-Evolving Block Decimation (iTEBD) [134] which works directly in the thermodynamic limit. Our results are discussed in Figure 8 (some results have been already reported in [10]). Different panels in the figure show the entanglement production rate S (t) plotted as a function of time for quenches with different initial states in the XXZ chain. The data shown in Figure 8 are the entanglement entropies for the half-infinite chain. Although no finite-size corrections are expected, finite-time corrections are visible in the Figure. The data exhibit a non-trivial dynamics at short times, often with oscillating behaviour. Interestingly, already at t ≈ 10 for most of the quenches the data exhibit stationary behaviour. The horizontal lines in the Figure mark the quasiparticle prediction The agreement between (77) and the iTEBD data is spectacular for all the quenches. Note that in the vicinity of ∆ = 1 a slower relaxation to the stationary behaviour takes place, especially for the quenches from the Néel state and from the tilted ferromagnet: longer times would be needed in order to provide a more robust check of (77).

Mutual information
In this section we focus on the post-quench dynamics of the mutual information between two blocks. Considering the tripartition A 1 ∪ A 2 ∪ B (with A 1 and A 2 two intervals of equal length and at distance d and B the rest of the chain), the von Neumann mutual information is defined as with S A 1(2) and S A 1 ∪A 2 being the entanglement entropies of A 1(2) and A 1 ∪ A 2 , respectively. Using the quasiparticle picture, it is straightforward to derive a prediction for the mutual information. When only one type of quasiparticles is present with fixed group velocity v (as in a conformal field theory), the prediction for the mutual information is simply obtained by counting the quasiparticles arriving to each interval, obtaining [6] I A 1 :A 2 ∝ −2 max((d + )/2, vt) + max(d/2, vt) + max((d + 2 )/2, vt).
Formula (79) predicts I A 1 :A 2 = 0 for vt ≤ d/2, followed by a linear increase for d/2 < vt ≤ (d + )/2 and a linear decrease up to vt = (d + 2 )/2. The first region corresponds to A 1 and A 2 being entangled with the environment B but not mutually entangled. At time t = d/(2v) quasiparticles originated at the same point in space start to connect A 1 and A 2 . The linear increase up to t = (d + )/(2v) correspond to entangled quasiparticles traveling in the two subsystems. At time t = (d + )/(2v) the entangled quasiparticles start leaving the two subsystems. Finally, at t = (d + 2 )/(2v) there are no entangled quasiparticles connecting A 1 and A 2 and the mutual information vanishes again.
In the presence of different species of quasiparticles with different velocities, one has to integrate (79) over the full quasiparticle content to obtain which is valid for infinite systems. For a finite chain, (80) applies before the revival time.
The exact numerical results for I A 1 :A 2 obtained using (80) for quenches in the XXZ chain are shown in Figure 9. Panel (a) shows results for the quench from the Néel state in the XXZ chain with ∆ = 4. The result for I A 1 :A 2 (full line in the Figure) is for two disjoint intervals of equal length = 10 at distance d = 10. Clearly, one has that for d/(2v M ), with v M ≈ 2 the maximum velocity, the mutual information is zero. A linear behaviour is clearly visible at larger times up to (d + )/(2v M ), where the mutual information reaches a maximum. A linear decrease is subsequently observed. Interestingly, the presence of slow quasiparticles leads to a slow decay of the mutual information at long times, instead of a sudden vanishing behaviour at t = (d + 2 )/(2v M ). A similar slow decay has been numerically observed in free bosonic models [84].
It is also interesting to investigate the effects of the bound states on the mutual information dynamics. The dotted and dashed lines in Figure 9 denote the contributions of the two-particle and three-particle bound states, respectively. Interestingly, the contributions of the bound states rapidly decay with their size. Moreover, the bound-state contributions are shifted at longer times, reflecting their smaller group velocities (see Figure 3). Similar qualitative results are observed for the quench from the dimer state (reported in Figure 9 (b)). Finally, Figure 9 shows results also for the quench from the tilted ferromagnet. The data are for ∆ = 4 and tilting angle ϑ = π/2. The results are for two adjacent equal-length intervals with = 10. In contrast with panels (a) and (b), an additional second peak is observed in the mutual Post-quench dynamics of I A 1 :A 2 for the quench from the tilted ferromagnet in the XXZ chain with ∆ = 4. Here A 1 and A 2 are two equal-length intervals with = 10 at distance d = 0. Note the second peak at t ≈ 30 resulting from the contribution of the two-particle bound states.
information. As it is clear from the Figure, this is due to the contribution of the two-particle bound states (dashed line). This last result suggests that the mutual information, at least in some case, can be used to reveal the bound state content of integrable models. This idea has already been put forward in [101] during the study of quenches in the spin-1 Lai-Sutherland model.

Entanglement dynamics in the Lieb-Liniger model
In this section we provide exact results for the entanglement dynamics after the quench from the Bose-Einstein condensate (BEC) in the Lieb-Liniger model. We discuss both the Lieb-Liniger model with repulsive interactions, as well as with attractive ones. Quantum quenches in the Lieb-Liniger model have been the focus of intensive investigations  and here we will largely use the results from Refs. [141] and [156] for repulsive and attractive cases respectively. We should mention that, in contrast with the XXZ chain, here we cannot provide a numerical check of our preditions. This is due to the fact that as of now for models in the continuum there are no efficient numerical methods, such tDMRG.

Lieb-Liniger model and its Bethe Ansatz solution
The Lieb-Liniger model consists of a system of N interacting bosons on a ring of length L. The model is defined by the hamiltonian where m is the mass of the bosons and c is the interaction strength. In the following we set = 2m = 1. In second quantisation (81) reads with Ψ(x) bosonic fields satisfying the standard commutation relations [Ψ(x), Ψ † (y)] = δ(x − y). In the limit c → ∞, (82) becomes equivalent to a system of hard-core bosons. For any value of c, the Lieb-Liniger model can be solved using Bethe ansatz [162]. In this work we consider both the repulsive regime with c > 0, as well as the attractive one with c < 0. We definec ≡ |c|. We also introduce the dimensionless coupling γ as The Bethe equations for the Lieb-Liniger model are [111,162] 2πI The eigenstates energy E and total momentum P are given as The structure of the solutions of the Bethe equations depends dramatically on the sign of the interactions. Specifically, for c > 0, i.e., for repulsive interactions, only real solutions of (84) are present. Consequently (84) is of the form (32) for the only species of particles after the straightforward identification of the various functions. In the thermodynamic limit the solutions of the BGT equations become dense on the real axis and, since there are no bound states, there is a single particle density ρ and hole density ρ (h) , with ρ (t) = ρ + ρ (h) which is a major simplification compared to the standard case. The Bethe equations for these root densities (84) where the kernel K is given as K(λ) = c/[π(λ 2 + c 2 )] For attractive interactions c < 0, the eigenstates of the model contain non-trivial multiparticle bound states that, as usual, can be understood with the string hypothesis, i.e. they have the form (31) with η =c. The Bethe-Gaudin-Takahashi (BGT) equations for the Lieb-Liniger gas are of the form (32) with π n (λ) = nλ and elementary kernel θ n (λ) given by [111] θ n (λ) = 2 arctan 2λ nc .
For the attractive Lieb-Liniger the energy and momentum in (85) can be rewritten as (34) with In the thermodynamic limit, there are infinite particle densities {ρ n } ∞ n=1 , hole densities {ρ with and a n (λ) = 2 π|c|n

Quench from the Bose condensate
Here we briefly detail the TBA treatment for the quantum quench from the Bose condensate state in the Lieb-Liniger model. In the BEC the bosons are uniformly distributed in the interval [0, L]. The steady state arising at infinite time after the quench is fully described by a particular thermodynamic macrostate.

Repulsive case
The quench action solution for the quench in the repulsive Lieb-Liniger has been provided in [141]. The thermodynamic macrostate describing the post-quench steady-state is identified by the densities ρ(λ), η(λ) [141]: written in terms of the the auxiliary function Here τ = 1/γ and I α (x) are the modified Bessel functions of the first kind. The calculation of the group velocities of the low-lying excitations around the macrostate that describes the steady-state follows the general derivation of section 3.2 with the major simplification of having a single set of rapidities. Figure 10 Figure). At large |λ| the interactions are negligible and the linear behaviour v ∝ 2λ is observed, reflecting the "free" dispersion E = λ 2 and the absence of a maximum velocity. We anticipate that this fact will have striking consequences in the behaviour of the entanglement entropy (see 5.3).

Attractive case
We now consider the quench from the Bose condensate in the attractive gas for which the thermodynamic macrostate describing the steady state is identified by the set of densities {ρ n } ∞ n=1 and {η n } ∞ n=1 . The solution for this problem has been provided in [156]. The densities η n satisfy the recursion relations [156] with x ≡ λ/c and The particle densities ρ n (x) are [156] ρ n (x) = τ 4π The group velocities of the low-lying excitations around the macrostate describing the steady-state can be calculated following the general derivation of Sec. 3.2. Figure 10 (b) shows numerical results for these group velocities v n for the different multi-particle bound states as a function of λ. The results are for fixed γ = −2. As for the repulsive case, at large momenta, the interactions are negligible and the free-like behaviour v ∝ 2λ is found.

Entanglement dynamics in the Lieb-Liniger model
We now turn to discuss the post-quench dynamics of the entanglement entropy for the repulsive Lieb-Liniger model as given by the quasiparticle prediction (4). In the present case, (4) greatly simplifies because of the presence of a single species of quasiparticles and it can be written as where v(λ) is the group velocity of the entangling quasiparticles of the previous section, and s(λ) is the thermodynamic Yang-Yang entropy of the steady state. The dynamics of the entanglement entropy obtained from (97)  Interestingly, for all values of γ the entropy exhibits a non-linear growth with time, even at short times, and it saturates at asymptotically long times. The non linear behaviour at short times is due to the absence of a maximum velocity (see Figure 10). The almost linear behaviour of the entanglement entropy for small values of γ is due to the very low weight of fast quasiparticles, as it should be clear from Fig. 10. Anyhow, at a closer analysis, a strictly linear behaviour never takes place for any value of γ. The maximum stationary entanglement entropy is obtained in the limit γ → ∞, when the system is equivalent to a system of hardcore bosons. In this limit, we find S/ = 2, as already known [74]. In order to understand the saturation behaviour at long times it is useful to investigate the quasiparticle contribution to the steady-state entanglement entropy. This is reported in Figure 11 (b) which shows the entropy density S/ contribution versus the quasiparticle rapidity λ. The different curves correspond to different values of the interaction strength. For all values of γ the quasiparticles contributions decay rapidly as λ → ∞. Upon increasing γ, quasiparticles with larger rapidity contribute more significantly to the steady-state entropy, which is one of the factors explaining why the stationary entropy increases with γ. This is better shown in Figure 12 that reports the entropy density S/ versus γ. The entropy density monotonically increases with γ and it vanishes for γ → 0, i.e. in the absence of a quench. In the limit γ → ∞ the result S/ = 2 [74] for hard-core bosons is recovered. We now turn to discuss the entanglement dynamics after the quench from the Bose condensate in the attractive Lieb-Liniger model. In this case, all the multi-boson bound states contribute to the entanglement which is then described by (4) that we repeat here for convenience: Numerical results for the entanglement evolution obtained using (98) are shown in Figure 13. Panel (a) shows results for S/ for the quench with γ = −2 plotted as a function of the rescaled time t/ . The different curves in the panel are the entanglement entropies in which the different multi-boson bound states up to size n have been taken into account in the sum (98). Only results for n ≤ 5 are shown. We verified that for this value of γ and in the time window reported in the plot, the contributions of the bound states with n > 5 are negligible.
As for the repulsive case (cf. Figure 11) there is no linear increase in the short time regime. The contribution of bound states with different rapidity is investigated in Figure 13 (b) plotting S/ as a function of rapidity λ for different values of n. Interestingly, the maximum contribution of the bound states increases with their size, although the support of S/ as a function of λ shrinks with increasing n. As a consequence of this very peculiar velocity distribution, we have that the larger bound states have a dominant velocity that is smaller and smaller as n increases. Thus their effect will manifest at longer times. This is already clear from the panel (a) in Fig. 13 where we can notice that the contributions with n = 3, 4, 5 have a visible effect some time after the quench. Consequently, we expect that larger bound states can have non-negligible contributions at some larger time not displayed in the figure.
We turn now to discuss the steady-state entropy as a function of the interaction strength. Clearly, the entropy density increases with γ, similar to the repulsive case (see Figure 13). The behaviour in the limit γ → ∞ can be understood analytically. In the limit γ → ∞, one has that the support of the root densities ρ n (λ) and ρ in the limit τ → 0 one has that where x ≡ λ/c. Interestingly, Eq. (99) implies that i.e. that in the limit of infinite attractive interaction the quasiparticles with n = 1 behave as free fermions. For generic n, the total density ρ (t) n is consistent with the ansatz Crucially, from (99)- (102) it is clear that the support of the higher densities ρ n and ρ (h) n for n > 1 shrinks faster, i.e., with a higher power of τ , in the limit τ → 0, implying that the multi boson bound states do not contribute to the leading behaviour of the steady-state entanglement entropy. Also, we should remark that, although the functional form of the densities in the limit γ → ∞ appear to be simple, we were not able to generalize the results (99)-(102) to arbitrary n. We can derive the average energy and particle density using (99)- (102). The boson density in the limit γ → ∞ is determined by the strings with n = 1 and it is given as as it should. Using (99)- (102), it is straightforward to check that the contributions of the bound states are vanishing as ∝ τ n−1 . On the other hand, for the energy density the contribution of each bound state diverges in the limit γ → ∞ as expected because the energy of the post-quench hamiltonian calculated on the BEC state diverges as γ → ∞.
Using (103), (99), and (100) in the definition of the Yang-Yang entropy, the stationary entanglement in the limit τ → 0 is determined by the strings with n = 1, and it is given as Interestingly, Eq. (106) is the same as for the BEC quench in the repulsive Lieb-Liniger [74].

Mutual information
Finally, we investigate the behaviour of the mutual information between two intervals. The quasiparticle formula for the mutual information is the same as that for the XXZ chain (80). Figure 14 shows I A 1 :A 2 for two disjoint intervals with equal length = 10 at distance d = 10 for the Lieb-Liniger gas with γ = −2. The continuous line denotes I A 1 :A 2 while the other curves are the individual contributions of the bound states with n = 1, 2, 3. Interestingly, the mutual information exhibits a peak at short times, which is followed by a quite slow vanishing behaviour as t → ∞. This slow relaxation is due to the significant contributions of the multibosons bound states. For all the bound states, a peak is observed at relatively short times, followed by a vanishing behaviour at long times. However, the position of the peak is shifted to longer times for the larger bound states.

Conclusions
In this paper we provided a thorough analysis of the framework put forward in [10] for the time evolution of the entanglement entropy which combines the quasiparticle picture of [6] with the exact knowledge of the stationary state coming from integrability. This approach is expected to hold in generic one-dimensional integrable systems. Here, we provided predictions, valid under rather general conditions, for arbitrary free systems, both bosonic and fermionic. These results have been tested against exact computations for the Ising and the harmonic chains. We also provided new results for the Heisenberg anisotropic spin chain (XXZ chain), which was the only model analysed in [10]. We finally derived theoretical predictions for the entanglement dynamics in the Lieb-Liniger model which have not been checked against numerical simulations, although it would be very interesting to do so. Specifically, it would be useful to verify the non-linear behaviour of the entropy at short times. A promising direction to perform this check is to extend the framework of continuous matrix product states [163] to simulate non-equilibrium systems. Alternatively, one could study the non-equilibrium dynamics of a very dilute Bose Hubbard model (as done in [164]), but this is computationally demanding.
A crucial observation is that Eq. (4) has been conjectured on the basis that the initial state acts as a source of pair of quasiparticle excitations with opposite momentum. In Bethe ansatz language, this assumption reflects the property that only parity-invariant eigenstates (as defined in [138,141]) have non zero-overlap with the initial state. Recently, there is a broad consensus emerging about the idea that only quenches from these initial states are exactly solvable for genuinely interacting integrable models, as first proposed in the context of quantum field theory [165] and later for lattice integrable models [166]. However, states with non-zero overlap with generic eigenstates do exist and it is fundamental to understand how (4) generalises. In this respect, free models can be a useful playground because they can be solved even relaxing this assumption. Examples of exact results for quenches from non parity invariant states have been provided recently for the Hubbard chain with infinite repulsion [104] (which is mappable to free fermions), and the entanglement dynamics can be described by a suitable generalisation of (3) [105].
A main open problem is the generalisation of the approach of this paper to Rényi entanglement entropies. While in Refs. [167][168][169] it has been shown how to derive analytically the stationary value of these entanglement monotones, a complete quasiparticle description for their full-time evolution is still lacking. On the same line of thoughts, it would be important to provide a semiclassical picture for more complex entanglement measures, such as the negativity [170][171][172], which quantify the entanglement also in mixed states. In this respect, a promising direction is to study the dynamics of the negativity in the harmonic chain, for which exact calculations are possible [84].