Entanglement spreading in non-equilibrium integrable systems

These are lecture notes for a short course given at the Les Houches Summer School on"Integrability in Atomic and Condensed Matter Physics", in summer 2018. Here, I pedagogically discuss recent advances in the study of the entanglement spreading during the non-equilibrium dynamics of isolated integrable quantum systems. I first introduce the idea that the stationary thermodynamic entropy is the entanglement accumulated during the non-equilibrium dynamics and then use such an idea to provide quantitive predictions for the time evolution of the entanglement entropy in arbitrary integrable models, regardless of the interaction strength.


Introduction
Starting from the mid-noughties, the physics community witnessed an incredibly large theoretical and experimental activity aimed to understand the non-equilibrium dynamics of isolated many-body quantum systems. The most studied protocol is certainly that of a quantum quench [1,2] in which an extended quantum system evolves with a Hamiltonian H after having being prepared at time t = 0 in a non-equilibrium state |Ψ 0 , i.e. [H, |Ψ 0 Ψ 0 |] = 0] (|Ψ 0 can also be thought as the ground state of another Hamiltonian H 0 and hence the name quench). At time t, the time evolved state is simply where we work in units of = 1. A main question is whether for large times these many-body quantum systems can attain a stationary state and how this is compatible with the unitary time evolution of quantum mechanics. If a steady state is eventually reached (in some sense to be specified later), it is then natural to ask under what conditions the stationary properties are the same as in a statistical ensemble. This is the problem of thermalisation of an isolated quantum system, a research subject that has been initiated in 1929 by one of the fathers of quantum mechanics, John Von Neumann, [3]. However, only in the last fifteen years the topic came to a new and active life, partially because of the pioneering experimental works with cold atoms and ions which can probe closed quantum systems for time scales large enough to access the relaxation and thermalisation, see, e.g., the experiments in Refs. [4][5][6][7][8][9][10][11][12][13]. Nowadays, there are countless theoretical and experimental studies showing that for large times and in the thermodynamic limit, many observables relax to stationary values, as reported in some of the excellent reviews on the subject [14][15][16][17][18][19]. In some cases (to be better discussed in the following), these stationary values coincide with those in a thermal ensemble or suitable generalisations, despite the fact that the dynamics governing the evolution is unitary and the initial state is pure. In these lecture notes, I focus (in an introductory and elementary fashion) on the entanglement spreading after a quench. The interested reader can find excellent presentations of many other aspects of the problem in the aforementioned reviews [14][15][16][17][18][19]. Furthermore, I will not make any introduction to integrability techniques in and out of equilibrium because they are the subject of other lectures in the 2018 Les Housches school [20][21][22][23].
These lecture notes are organised as follows. In Sec. 2 it is shown how the reduced density matrix naturally encodes the concept of local relaxation to a stationary state. In Sec. 3 the entanglement entropy is defined and its role for the non-equilibrium dynamics is highlighted. In Sec. 4 we introduce the quasiparticle picture for the spreading of entanglement which is after applied to free fermionic systems (Sec. 5) and interacting integrable models (Sec. 6); in particular in Sec. 6.4 we briefly discuss some recent results within the entanglement dynamics of integrable systems.

Stationary state and reduced density matrix
The reduced density matrix is the main conceptual tool to understand how and in which sense for large times after the quench an isolated quantum system can be described by a mixed state such as the thermal one. Let us consider a non-equilibrium many-body quantum system (in arbitrary dimension). Since the time evolution is unitary, the entire system A ∪Ā is in a pure state at any time (cf. |Ψ(t) in Eq. (1)). Let us consider a spatial bipartition of the system into two complementary parts denoted as A andĀ. Denoting with ρ(t) = |Ψ(t) Ψ(t)| the density matrix of the entire system, the reduced density matrix is defined by tracing out the degrees of freedom inĀ as ρ A (t) = TrĀ ρ(t) .
The reduced density matrix ρ A (t) generically corresponds to a mixed state with non-zero entropy, even if ρ(t) is a projector on a pure state. Its time dependent von Neumann entropy is called entanglement entropy and it is the main quantity of interest of these lectures. Some of its features will be discussed in the following section.
A crucial observation is that the physics of the subsystem A is fully encoded in the reduced density matrix ρ A (t), in the sense that ρ A (t) is enough to determine all the correlation functions local within A. In fact, the expectation value of a product of local operators i O(x i ) with x i ∈ A is given by The equilibration of a closed quantum system to a statistical ensemble starts from the concept of reduced density matrix. Indeed, we will say that, following a quantum quench, an isolated infinite system relaxes to a stationary state, if for all finite subsystems A, the limit of the reduced density matrix ρ A (t) for infinite time exists, i.e. if it exists It is very important to stress that Eq. (5) implies a very precise order of limits; since the infinite time limit is taken for an infinite system, it means that the thermodynamic limit must be taken before the infinite time one; the two limits do not commute and phenomena like quantum recurrences and revivals prevent relaxation for finite systems (anyhow timeaveraged quantities could still attain values described by a statistical ensemble). Another important observation is that although Eq. (5) is apparently written only for a subsystem A, it is actually a statement for the entire system. In fact, the subsystem A is finite, but it is placed in an arbitrary position and it has an arbitrary (finite) dimension. Furthermore, the limit of a very large subsystem A can also be taken, but only after the infinite time limit. Once again the two limits do not commute and their order is important. Summarising, there are three possible limits involved in the definition of the stationary state after a quantum quench; these limits do not commute and only one precise order leads to a consistent definition of equilibration of an isolated quantum system. We are now ready to understand in which sense ρ A (∞) may correspond to a statistical ensemble. A first guess would be that ρ A (∞) is itself an ensemble density matrix (e.g. thermal).
However, this definition would not be satisfactory because we should first properly consider boundary effects; moreover it would be valid only for thermodynamically large subsystems. We take here a different route following Refs. [24][25][26][27][28]. Let us consider a statistical ensemble with density matrix ρ E for the entire system. We can construct the reduced density matrix of a subsystem A as We say that the stationary state is described by the statistical ensemble ρ E if, for any finite subsystem A, it holds This implies that arbitrary local multi-point correlation functions within subsystem A, like those in Eq. (4), may be evaluated as averages with the density matrix ρ E . This definition should not suggest that ρ E is the density matrix of the whole system that would be a nonsense because the the former is a mixed state and the latter a pure one.
In these lectures, we are interested only into two statistical ensembles, namely the thermal (Gibbs) ensemble and the generalised Gibbs one. We say that a non-equilibrium quantum system thermalises after a quantum quench when ρ E is the Gibbs distribution with Z = Tre −βH . The inverse temperature β = 1/T is not a free parameter: it is fixed by the conservation of energy. In fact, the initial and the stationary values of the Hamiltonian are equal, i.e. Tr[Hρ E ] = Ψ 0 |H|Ψ 0 .
This equation can be solved for β, fixing the temperature in the stationary state. Once again, thermalisation leads to the remarkable consequence that all local observables will attain thermal expectations, but some non-local quantities will remain non-thermal for arbitrary large times. Generically, all non-integrable systems should relax to a thermal state, as supported by theoretical arguments such as the eigenstate thermalisation hypothesis [29][30][31][32], by a large number of simulations (see, e.g., [33][34][35][36][37]), and by some cold atom experiments [4,5,9,11]. However, there are some exceptional cases in which chaotic systems fail to thermalise like many-body localised ones [38,39], or those in the presence of quantum scars [40][41][42][43], or when elementary excitations are confined [44][45][46][47][48][49]. The dynamics and the relaxation of integrable models are very different from chaotic ones because of the constraints imposed by the conservation laws. Integrable models have, by definition, an infinite number of integrals of motion in involution, i.e. [I n , I m ] = 0 (usually one of the I m is the Hamiltonian). Consequently, rather than a thermal ensemble, the system for large time is expected to be described by a generalised Gibbs ensemble (GGE) [50] with density matrix Here the operators I n form a complete set (in some sense to be specified) of integrals of motion and Z is the normalisation constant Z = Tr e − n λnIn ensuring Trρ GGE = 1. As the inverse temperature for the Gibbs ensemble, the Lagrange multipliers {λ n } are not free, but are fixed by the conservation of {I n }, i.e. they are determined by the (infinite) set of equations In the above introduction to the GGE, we did not specify which conserved charges should enter in the GGE density matrix (10). One could be naively tempted to require that all linearly independent operators commuting with the Hamiltonian should be considered in the GGE, regardless of their structure; this is what one would do in a classical integrable system to fix the orbit in phase space. In this respect, the situation is rather different between classical and quantum mechanics. Indeed, any generic quantum model has too many integrals of motion, independently of its integrability. For example, all the projectors on the eigenstates O n = |E n E n |, are conserved quantities for all Hamiltonians since H = n E n |E n E n |. For a model with N degrees of freedom, the number of these charges is exponentially large in N , instead of being linear, as one would expect from the classical analogue. All these integrals of motion cannot constrain the local dynamics and enter in the GGE, otherwise no system will ever thermalise and all quantum models would be, in some weird sense, integrable. The solution of this apparent paradox is that, as long as we are interested in the expectation values of local observables, only integrals of motion with some locality or extensivity properties must be included in the GGE [27,28,51]. In the spirit of Neither theorem of quantum field theory, an integral of motion is local if it can be written as an integral (sum in the case of a lattice model) of a given local density. However, it has been recently shown that also a more complicated class of integrals of motion, known as quasi-local [52], have the right physical features to be included in the GGE [53,54]. The discussion of the structure of these new conserved charges is far beyond the goal of these lectures. Our main message here is that we nowadays have a very clear picture of which operators form a complete set to specify a well defined GGE in all integrable models, free and interacting.
We conclude this section by mentioning what happens for finite systems, also, but not only, to describe cold atomic experiments with only a few hundred constituents. When there is a maximum velocity of propagation of information v M (in a sense which will become clearer later), as long as we consider times such that v M t L, with L the linear size of the system, all measurements would provide the same outcome as in an infinite system (away from the boundaries). Thus, a subsystem of linear size can show stationary values as long as L is large enough to guarantee the existence of the time window v M t L.

Entanglement entropy in many-body quantum systems
In order to understand the connection between entanglement and the equilibration of isolated quantum systems, we should first briefly discuss the bipartite entanglement of many-body systems (see e.g. the reviews [55][56][57][58]). As we did in the previous section, let us consider an extended quantum system in a pure state |Ψ and take a bipartition into two complementary parts A andĀ. Such spatial bipartition induces a bipartition of the Hilbert space as H = H A ⊗ HĀ. We can understand the amount of entanglement shared between these two parts thanks to Schmidt decomposition. It states that for an arbitrary pure state |Ψ and for an arbitrary bipartition, there exist two bases |w A α of H A and |wĀ α of HĀ such that |Ψ can be written as The Schmidt eigenvalues λ α quantify the non-separability of the state, i.e. the entanglement.
If there is only one non-vanishing λ α = 1, the state is separable, i.e. it is unentangled. Conversely, the entanglement gets larger when more λ α are non-zero and get similar values.
Schmidt eigenvalues and eigenvectors allow us to write the reduced density matrix ρ A = TrĀ|Ψ Ψ| as and similarly for ρĀ with |wĀ α replacing |w A α . A proper measure of the entanglement between A andĀ is the von Neumann entropy of ρ A or ρĀ which is known as entanglement entropy. Obviously many other functions of the Schmidt eigenvalues are proper measures of entanglement. For example, all the Rényi entropies quantify the entanglement for any n > 0. These Rényi entropies have many important physical properties. First, the limit for n → 1 provides the von Neumann entropy and, for this reason, they are the core of the replica trick for entanglement [59]. Then, for integer n ≥ 2, they are the only quantities that are measurable in cold-atom and ion-trap experiments [11][12][13][60][61][62][63]. Finally their knowledge for arbitrary integer n provides the entire spectrum of ρ A [64], known as entanglement spectrum [65]. Rigorously speaking entanglement and Rényi entropies are good entanglement measures in the sense that they are entanglement monotones [66]. While these lectures are not the right forum to explain what an entanglement monotone is (the interested reader can check, e.g., the aforementioned [66]), we want to grasp some physical intuition about the physical meaning of the entanglement entropy. To this aim, let us consider the following simple two-spin state with α ∈ [0, π/2]. It is a product state for α = 0 and α = π/2 and we expect that the entanglement should increase with α up to a maximum at α = π/4 (the singlet state). The reduced density matrix of one of the two 1/2 spins is with entanglement entropy which has all the expected properties and takes the maximum value log 2 on the singlet state.
Let us now consider a many-body system formed by many spins 1/2 on a lattice and a state which is a collection of singlets between different pairs of spins at arbitrary distances (incidentally these states have important physical applications in disordered systems [67]). All singlets within A orĀ do not contribute to the entanglement entropy S A . Each shared singlets instead counts for a log 2 bit of entanglement. Hence, the total entanglement entropy is S A = n A:Ā log 2 with n A:Ā being the number of singlets shared between the two parts. As a consequence, the entanglement entropy measures all these quantum correlations between spins that can be very far apart.
Let us now move back to non-equilibrium quantum systems and see what entanglement can teach us. The stationary value of the entanglement entropy S A (∞) = −Trρ A (∞) log ρ A (∞) for a thermodynamically large subsystem A is simply deduced from the reasoning in the previous section. Indeed, we have established that a system relaxes for large times to a statistical ensemble ρ E when, for any finite subsystem A, the reduced density matrix ρ A,E (cf. Eq. (6)) equals the infinite time limit ρ A (∞) (cf. Eq. (5)). This implies that the stationary entanglement entropy must equal S A,E = −Trρ A,E log ρ A,E . For a large subsystem with volume V A , S A,E scales like V A because the entropy is an extensive thermodynamic quantity. Hence, S A,E equals the density of thermodynamic entropy S E = −Trρ E log ρ E times the volume of A. Given that S A,E = S A (∞), the stationary entanglement entropy has the same density as the thermodynamic entropy. In conclusion, we have just proved the following chain of identities From the identification of the asymptotic entanglement entropy with the thermodynamic one we infer that the non-zero thermodynamic entropy of the statistical ensemble is the entanglement accumulated during the time by any large subsystem. We stress that this equality is true only for the extensive leading term of the entropies, as in Eq. (19); subleading terms are generically different. The equality of the extensive parts of the two entropies has been verified analytically for non-interacting many-body systems [68][69][70][71] and numerically for some interacting cases [72][73][74].

The quasiparticle picture
In this section, we descibe the quasiparticle picture for the entanglement evolution [75] which, as we shall see, is a very powerful framework leading to analytic predictions for the time evolution of the entanglement entropy that are valid for an arbitrary integrable model (when complemented with a solution for the stationary state coming from integrability). This picture is expected to provide exact results in the space-time scaling limit in which t, → ∞, with the ratio t/ fixed and finite. Let us describe how the quasiparticle picture works [18,75]. The initial state |Ψ 0 has an extensive excess of energy compared to the ground state of the Hamiltonian H governing the time evolution, i.e. it has an energy located in the middle of the many-body spectrum. The state |Ψ 0 can be written as a superposition of the eigenstates of H; for an integrable system these eigenstates are multiparticle excitations. Therefore we can interpret the initial state as a source of quasiparticle excitations. We assume that quasiparticles are produced in pairs of opposite momenta. We will discuss when and why this assumption is correct for some explicit cases in the following, see Sec. 6.4.2 (clearly the distribution of the quasiparticles depends on the structure of the overlaps between the initial state and the eigenstates of the post-quench Hamiltonian). The essence of the picture is that particles emitted from different points are unentangled. Conversely, pairs of particles emitted from the same point are entangled and, as they move far apart, they are responsible for the spreading entanglement and correlations throughout the system (see Fig. 1 for an illustration). A particle of momentum p has energy Figure 1: Quasiparticle picture for the spreading of entanglement. The initial state (at time t = 0) acts as a source of pairs of quasiparticles produced homogeneously throughout the system. After being produced, the quasiparticles separate ballistically moving with constant momentum-dependent velocity and spreading the entanglement. the system; we assume that there is no scattering between them and that they have an infinite lifetime (assumptions which are fully justified in integrable models [76]). Thus, a quasiparticle created at the point x with momentum p will be found at position x = x + v p t at time t while its entangled partner will be at The entanglement between A andĀ at time t is related to the pairs quasiparticles that are shared between A andĀ after being emitted together from an arbitrary point x. For fixed momentum p, this is proportional to the length of the interval (or region in more complicated The proportionality constant depends on both the rate of production of pairs of quasiparticles of momentum (p, −p) and their contribution to the entanglement entropy itself. The combined result of these two effects is a function s(p) which depends on the momentum p of each quasiparticle in the pair. This function s(p) encodes all information about the initial state for the entanglement evolution.
Putting together the various pieces, the total entanglement entropy is [75] which is valid for an arbitrary bipartition of the whole system in A andĀ. We can see in this formula all elements we have been discussing: (i) particles are emitted from arbitrary points x (the integral runs over [−∞, ∞]); (ii) they move ballistically as forced by the delta functions constraints over the linear trajectories; (iii) they are forced to arrive one in A the other in A (the domain of integration in x and x ); (iv) finally, we sum over all allowed momenta p (whose domain can depend on the model) with weight s(p). We specialise Eq. (20) to the case where A is a single interval of length . All the integrals over the positions x, x , x in Eq. (20) are easily performed, leading to the main result of the quasiparticle picture [75]   is strongly divergent too, but this is not physical). Consequently, the stationary value of the entanglement entropy is where in the rhs we used that s(p) = s(−p) by construction. At this point, we assume that a maximum speed v M for the propagation of quasiparticles exists. The Lieb-Robinson bound [77] guarantees the existence of this velocity for lattice models with a finite dimensional local Hilbert space (such as spin chains). Also in relativistic field theories, the speed of light is a natural velocity bound. Since |v(p)| ≤ v M , the second integral in Eq. (21) is vanishing as long as t < t * = /(2v M ) (the domain of integration again shrinks to zero). Hence, for For finite t such that t > t * , both integrals in Eq. (21) are non zero. The physical interpretation is that while the fastest quasiparticles (those with velocities close to v M ) reached a saturation value, slower quasiparticles continue arriving at any time so that the entanglement entropy slowly approaches the asymptotic value (22). The typical behaviour of the entanglement entropy resulting from Eq. (21) is the one reported in Fig. 2 where the various panels and curves correspond to the actual theoretical results for an interacting integrable spin chain (the anisotropic Heisenberg model) that we will discuss in the forthcoming sections. The last missing ingredients to make Eq. (21) quantitatively robust are the functions s(p) and v p which should be fixed in terms of the quench parameters. The idea proposed in Ref. [78] (see also [79,80]) is that s(p) can be deduced from the thermodynamic entropy in the stationary state, using the fact that the stationary entanglement entropy has the same density as the thermodynamic one, cf. Eq. (19). To see how this idea works, we will apply it to free fermionic models in the next section and then to generic integrable models in the following one.

Quasiparticle picture for free fermionic models
The ab-initio calculation of entanglement entropy is an extremely challenging task. For Gaussian theories (i.e. non-interacting ones) it is possible to relate the entanglement entropy to the two-point correlation functions within the subsystem A both for fermions and bosons [81][82][83][84]. Anyhow, for quench problems, extracting analytic asymptotic results from the correlation matrix technique is a daunting task that has been performed for some quenches in free fermions [68], but not yet for free bosons. We are going to see here that instead the quasiparticle picture provides exact analytic predictions in an elementary way, although not derived directly from first principles.
In this section, we consider an arbitrary model of free fermions. We focus on translational invariant models that can be diagonalised in momentum space k. It then exists a basis in which the Hamiltonian, apart from an unimportant additive constant, can be written as in terms of canonical creation b † k and annihilation b k operators (satisfying {b k , b † k } = δ k,k ). The variables k are single-particle energy levels.
We consider the quantum quench in which the system is prepared in an initial state |Ψ 0 and then is let evolve with the Hamiltonian H. For all these models, the GGE built with local conservation laws is equivalent to the one built with the mode occupation numberŝ n k = b † k b k since they are linearly related [28]. Thus the local properties of the stationary state are captured by the GGE density matrix where Z = Tre − k λ knk (under some some reasonable assumptions on the initial state [26,85,86]). 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. (25) 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 (24), but it is sufficient to calculate the expectation values ofn k in the initial state ψ 0 |n k |ψ 0 which equals, by construction, n k = n k GGE .
At this point, following the ideas of the previous sections (cf. Eq. (19)), we identify the stationary thermodynamic entropy with the density of entanglement entropy to be plugged in Eq. (21), obtaining the general result 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 can be written in terms of pairs of quasiparticles. More general and peculiar structures of initial states can be also considered, see Sec. 6.4. Following the same logic, it is clear that Eq. (27) is also valid for free bosons (i.e. Hamiltonians like (23) with the ladder bosonic operators) with the minor replacement of the function H(n) (26) with [79,80] H bos (n) = −n log n + (1 + n) log(1 + n) .

The example of the transverse field Ising chain
Eq. (27) 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. The Hamiltonian (29) is diagonalised by a combination of Jordan-Wigner and Bogoliubov transformations, leading to Eq. (23) with the single-particle energies We focus on a quench of the magnetic field in which the chain is initially prepared in the ground state of (29) with h 0 and then, at t = 0, the magnetic field is suddenly switched 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 [87] where ∆ k is the difference of the pre-and post-quench Bogoliubov angles [87] with 0 k and k the pre-and post-quench energy levels, respectively. The quasiparticle prediction for the entanglement dynamics after the quench is then given by Eq. (27) with n k in Eq. (31). This result coincides with the ab initio derivation performed in [68]. The Ising model is only one of the many quenches in non-interacting theories of bosons and fermions in which the entanglement evolution is quantitatively captured by Eq. (27), as seen numerically in many cases [75,[88][89][90][91][92][93][94][95][96][97][98].

Quasiparticle picture for interacting integrable models
We are finally ready to extend the application of the quasiparticle picture to the entanglement entropy dynamics in interacting integrable models. We exploit the thermodynamic Bethe ansatz (TBA) solution of these models and remand for all the technicalities to other lectures in this school [20][21][22], or to the existing textbooks [99][100][101] on the subject. Here we just summarise the main ingredients we need and then move back to the entanglement dynamics

Thermodynamic Bethe ansatz
In all Bethe ansatz integrable models, energy eigenstates are in one to one correspondence with a set of complex quasi momenta λ j (known as rapidities) which satisfy model dependent non-linear quantisation conditions known as Bethe equations. The solutions of the Bethe equations organise themselves into mutually disjoint patterns in the complex plane called strings [99]. Intuitively, a n-string solution corresponds to a bound state of n elementary particles with n = 1. Each bound state (of n particles) has its own quasi momentum λ (n) α . The Bethe equations induce effective equations for the quantisation of the quasi momenta of the bound states known as Bethe-Takahashi equations [99]. In the thermodynamic limit, the solutions of these equations become dense on the real axis and hence can be described by smooth distribution functions ρ (p) n (λ). One also needs to introduce the hole distribution functions ρ (h) n (λ): they are a generalisation to the interacting case of the hole distributions of an ideal Fermi gas at finite temperature [99][100][101]. Because of the non-trivial (i.e. due to interactions) quantisation conditions, the hole distribution is not simply related to the particle one. Finally, it is also useful to introduce the total density ρ (t) In conclusion, in the thermodynamic limit a macrostate is identified with a set of densities ρ ρ ρ ≡ {ρ The Yang-Yang entropy is the thermodynamic entropy of a given macrostate, as it simply follows from a generalised microcanonical argument [102]. In particular, it has been shown that for in thermal equilibrium it coincides with the thermal entropy [99].

The GGE as a TBA macrostate
The generalised Gibbs ensemble describing the asymptotic long time limit of a system after a quench is one particular TBA macrostate and hence it is fully specified by its rapidities (particle and hole) distribution functions. There are (at least) three effective ways to calculate these distributions (see also the lectures by Fabian Essler [20]). The first one is based on the quench action approach [103,104], a recent framework that led to a very deep understanding and characterisation of the quench dynamics of interacting integrable models. This technique is based on the knowledge of the overlaps between the initial state and Bethe eigenstates. Starting from these, it provides a set of TBA integral equations for the rapidity distributions in the stationary state that can be easily solved numerically and, in a few instances, also analytically. In turns, the developing of such approach also motivated the determination of the exact overlaps in many Bethe ansatz solvable models [105][106][107][108][109][110][111][112][113][114][115][116][117][118]. Based on these overlaps, a lot of exact results for the stationary states have been systematically obtained in integrable models [105,[119][120][121][122][123][124][125][126][127][128][129][130][131][132]. We must mention that only thanks to the quench action solutions of some quenches in the XXZ spin chain [122][123][124][125], it has been discovered that the GGE built with known (ultra)local charges [133][134][135] is insufficient to describe correctly [136,137] the steady state; this result motivated and boosted the discovery of new families of quasi-local conservation laws that must be included in the GGE [53,54,[138][139][140]. This finding is extremely important because when a complete set of charges is known, the stationary state can be built circumventing the knowledge of the overlaps required for quench action solution, as e.g. done in Refs. [141][142][143][144][145][146][147]. The direct construction of the GGE based on all the linear independent quasilocal conserved charges is the second technique to access the asymptotic TBA macrostate. The third technique is based on the quantum transfer matrix formalism [121,148,149], but will not be further discussed here. We finally mention that in the quench action formalism, the time evolution of local observables can be obtained as a sum of contributions coming from excitations over the stationary state [103]. This sum has been explicitly calculated for some non-interacting systems [103,150,151], but, until now, resisted all attempts for an exact computation in interacting models [119,152] and hence it has only been numerically evaluated [153].

The entanglement evolution
As we have seen above, in interacting integrable models there are generically different species of quasiparticles corresponding to the bound states of n elementary ones. According to the standard wisdom (based, e.g., on the S matrix, see [76]), these bound states must be treated as independent quasiparticles. It is then natural to generalise Eq. (21), for the entanglement evolution with only one type of particles, to the independent sum of all of them, resulting in where the sum is over the species of quasiparticles n, v n (λ) is their velocity, and s n (λ) the entropy density in rapidity space (the generalisation of s(p) in Eq. (21)). To give predictive power to Eq. (34), we have to device a framework to determine v n (λ) and s n (λ) in the Bethe ansatz formalism. The first ingredient to use is that in the stationary state the density of thermodynamic entropy (see Eq. (33)) equals that of the entanglement entropy in (34). Since this equality must hold for arbitrary root densities, we can identify s n (λ) with the density of Yang-Yang entropy for the particle n, i.e.
Moreover, the entangling quasiparticles in (34) can be identified with the excitations built on top of the stationary state. Their group velocities v n (λ) depend on the stationary state, because the interactions induce a state-dependent dressing of the excitations. These velocities v n (λ) can be calculated by Bethe ansatz techniques [154], but we do not discuss this problem here (see [79,154] for all technical details).
Eq. (34) complemented by Eq. (35) and by the proper group velocities v n (λ) is the final quasiparticle prediction for the time evolution of the entanglement entropy in a generic integrable model. This prediction is not based on an ab-initio calculation and should be thought as an educated conjecture. It has been explicitly worked out using rapidity distributions of asymptotic macrostates for several models and initial states [78,79,130,149,155]. Some examples for the interacting XXZ spin chains, taken from [79], are shown in Fig. 2. The validity of this conjecture has been tested against numerical simulations (based on tensor network techniques) for a few interacting models. In particular, in Refs. [78,79], the XXZ spin chain for many different initial states and for various values of the interaction parameter ∆ has been considered. The numerical data (after the extrapolation to the thermodynamic limit) are found to be in perfect agreement with the conjecture (34), providing a strong support for its correctness. In Ref. [156], the quasiparticle conjecture (34) has been tested for a spin-1 integrable spin chain, finding again a perfect match. This latter example is particularly relevant because it shows the correctness of Eq. (34) also for integrable models with a nested Bethe ansatz solution.
We conclude the section stressing that Eq. (34) represents a deep conceptual breakthrough because it provides in a single compact formula how the entanglement entropy becomes the thermodynamic entropy for an arbitrary integrable model.

Further developments
In this concluding subsection, we briefly go through several generalisations for the entanglement dynamics based on quasiparticle picture that have been derived starting from Eq. (34). Here, we do not aim to give an exhaustive treatment, but just to provide to the interested reader an idea of the new developments and some open problems.

Rényi entropies
A very interesting issue concerns the time evolution of the Rényi entropies defined in Eq. (15). These quantities are important for a twofold reason: on the one hand, they represent the core of the replica approach to the entanglement entropy itself [59], on the other, they are the quantities that are directly measured in cold atom and ion trap experiments [11][12][13][60][61][62][63].
For non-interacting systems, the generalisation of the formula for the quasiparticle picture is straightforward. Taking free fermions as example, the density of thermodynamic Rényi entropy in momentum space in terms of the mode occupation n k is just [68,157] Consequently, the time evolution of the Rényi entropy is just given by the same formula for von Neumann one, i.e. Eq. (27), in which H(n k ) is replaced by s (α) (n k ). One would then naively expect that something similar works also for interacting integrable models. Unfortunately, this is not the case because it is still not known whether the Rényi analogue of the Yang-Yang entropy (33) exists. In Ref. [157] an alternative approach based on quench action has been taken to directly write the stationary Rényi entropy. First, in quench action approach, the α-moment of ρ A may be written as the path integral [157] Trρ α where E[ρ ρ ρ] stands for the thermodynamic limit of the logarithm of the overlaps, S Y Y [ρ ρ ρ] is the Yang-Yang entropy, accounting for the total degeneracy of the macrostate, and the path integral is over all possible root densities ρ ρ ρ defining the macrostates. The most important aspect of Eq. (37) is that the Rényi index α appears in the exponential term and so it shifts the saddle point of the quench action. There is then a modified quench action with saddle-point equation Finally, the stationary Rényi entropies are the saddle point expectation of this quench action where in the rhs we used the property that S A in a form that closely resembles the replica definition of the entanglement entropy [59].
Eq. (39) is a set of coupled equations for the root densities ρ * α ρ * α ρ * α that can be solved, at least numerically, by standard methods. This analysis has been performed for several quenches in the XXZ spin chain [158,159] and the results have been compared with numerical simulations finding perfect agreement.
The main drawback of this approach is that the stationary Rényi entropy for α = 1 is not written in terms of the root distribution of the stationary state ρ * 1 ρ * 1 ρ * 1 for local observables. Since the entangling quasiparticles are the excitations on top of ρ * 1 ρ * 1 ρ * 1 , to apply the quasiparticle picture we should first rewrite the Rényi entropy in terms of ρ ρ ρ * 1 . Unfortunately, it is still not know how to perform this step. We mention that an alternative promising route to bypass this problem is based on the branch point twist field approach [160,161].

Beyond the pair structure
A crucial assumption to arrive at Eq. (22) for the entanglement evolution is that quasiparticles are produced in uncorrelated pairs of opposite momenta. This assumption is justified by the structure of the overlaps between initial state and Hamiltonian eigenstates found for many quenches both in free [68,87,163,164] and interacting models [105,111,112,117,162]. Indeed, it has been proposed that this pair structure in interacting integrable models is what makes the initial state compatible with integrability [162] and, in some sense, makes the quench itself integrable (see [162] for details). This no-go theorem does not apply to non-interacting theories and indeed, in free fermionic models, it is possible to engineer peculiar initial states such that quasiparticles are produced in multiplets [131,132] or in pairs having non-trivial correlations [165,166]. In all these cases, it is possible to adapt the quasiparticle picture to write exact formulas for the entanglement evolution, but the final results are rather cumbersome and so we remand the interested reader to the original references [131,132,165,166].

Disjoint intervals: Mutual information and entanglement negativity
Let us now consider a tripartition A 1 ∪ A 2 ∪Ā of a many-body system (with A 1 and A 2 two intervals of equal length and at distance d andĀ the rest of the system). We are interested in correlations and entanglement between A 1 and A 2 . A first measure of the total correlations is the mutual information 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 and counting the quasiparticles that at time t are shared between A 1 and A 2 , it is straightforward to derive a prediction for the mutual information which reads [75,78,79] where s n (λ) and v n (λ) have been already defined for the entanglement entropy. An interesting idea put forward in the literature is that one can use this formula to make spectroscopy of the particle content [79,130]. In fact, since the typical velocities of different quasiparticles n are rather different, Eq. (42) implies that the mutual information is formed by a train of peaks in time; these peaks become better and better resolved as d grows compared to which is kept fixed. The mutual information, however, is not a measure of entanglement between A 1 and A 2 . An appropriate measure of entanglement is instead the logarithmic negativity E A 1 :A 2 [167] defined as Here ρ T 2 A is the partial transpose of the reduced density matrix ρ A . The time evolution of the negativity after a quench in an integrable model has been analysed in Refs. [168,169]. To make a long story short, the quasiparticle prediction is the same as Eq. (42) but with s n (λ) replaced by another functional ε(λ) of the root densities. This functional is related to the Rényi-1/2 entropy. Hence, as discussed in Sec. 6.4.1, we know it only for free theories. Exact predictions for free bosons and fermions have been explicitly constructed in Ref. [169] and tested against exact lattice calculations, finding perfect agreement.

Finite systems and revivals
How the quasiparticle picture generalise to a finite system of total length L? Starting from Eq. (21), it is clear that the only change is to impose the periodic trajectories of the quasiparticles which are x ± = [(x ± v p t) mod L]. Using these trajectories, the final result is easily worked out as [170][171][172] where {x} denotes the fractional part of x, e.g., {7.36} = 0.36. This form has been carefully tested for free systems [171] in which it is possible to handle very large sizes. For interacting models, tensor network simulations work well only for relatively small values of L, but still the agreement is satisfactory [171]. We must mention that Eq. (44) also applies to the dynamics of the thermofield double [170,173], a state which is of great relevance also for the physics of black holes [174]. Finally, the structure of the revivals in minimal models of conformal field theories is also known [175].

Towards chaotic systems: scrambling and prethermalisation
What happens when integrability is broken? Can we say something about the time evolution of the entanglement entropy? It has already been found, especially in numerical simulations, that, in a large number of chaotic systems, the growth of the entanglement entropy is always linear followed by a saturation, see e.g. [37,[176][177][178][179][180][181]. This behaviour is the same as the one found in the quasiparticle picture, that, anyhow, cannot be the working principle here because the quasiparticles are unstable or do not exist at all. Recently, an explanation for this entanglement dynamics has arisen by studying random unitary circuits [182], systems in which the dynamics is random in space and time with the only constraint being the locality of interactions. In this picture, the entanglement entropy is given by the surface of the minimal space-time membrane separating the two subsystems. It has been proposed that this picture should describe, at least qualitatively, the entanglement spreading in generic non-integrable systems [183]. Although the prediction for the entanglement entropy of a single interval in an infinite system is the same for both the quasiparticle and the minimal membrane pictures, the two rely on very different physical mechanisms and should provide different results for other entanglement related quantities. In fact, it has been found that the behaviour of the entanglement of disjoint regions [184,185] or that of one interval in finite volume [171,182,186] is qualitatively different. For maximally chaotic systems, the mutual information and the negativity of disjoint intervals are constantly zero and do not exhibit the peak from the quasiparticle picture seen in Eq. (42). The explanation of this behaviour is rather easy: in non-integrable models, the quasiparticles decay and scatter and they cannot spread the mutual entanglement far away. It has been then proposed that the decay of the peak of the mutual information and/or negativity with the separation is a measure of the scrambling of quantum information [184,185], as carefully tested numerically [185]. Remarkably, such a peak and its decay with the distance has been also observed in the analysis of the experimental ion-trap data related to the negativity [187]. Also in the case of a finite size system, the decay and the scattering of the quasiparticles prevent them to turn around the system; consequently the dip in the revival of the entanglement of a single interval predicted by Eq. (44) is washed out [186]. In full analogy with the mutual information, the disappearance of such dip is a quantitive measure of scrambling [171].
A natural question is now what happens to the entanglement dynamics when the integrability is broken only weakly. In this case, one would expect the two different mechanisms underlying the above picture to coexist until the metastable quasiparticles decay. This problem has been addressed in Ref. [188] finding that, for sufficiently small interactions, the entanglement entropy shows the typical prethermalization behaviour [189][190][191][192]: it first approaches a quasi-stationary plateau described by a deformed GGE and then, on a separate timescale, its starts drifting towards its thermal value. A modified quasiparticle picture provides an effec-tive quantitative description of this behaviour: the contribution of each pair of quasiparticles to the entanglement becomes time-dependent and can be obtained by quantum Boltzmann equations [192], see for details [188].

Open systems
So far, we limited our attention to isolated quantum systems, but it is of great importance to understand when and how the quasiparticle picture can be generalised to systems that interact with their surrounding. In this respect, a main step forward has been taken in Ref. [193] (see also [194]), where it was shown that the quasiparticle picture can be adapted to the dynamic of some open quantum systems. In these systems, the spreading of entanglement is still governed by quasiparticles, but the environment introduces incoherent effects on top of it. For free fermions, this approach provided exact formulas for the evolution of the entanglement entropy and the mutual information which have been tested against ab-initio simulations

Inhomogeneous systems and generalised hydrodynamics
The recently introduced generalised hydrodynamics [195,196] (see in particular the lectures by Ben Doyon in this volume [197]) is a new framework that empower us to handle spatially inhomogeneous initial states for arbitrary integrable models (generalising earlier works in the context of conformal field theory [198,199]). For what concerns the entanglement evolution, the attention in the literature focused on the case of the sudden junction of two leads [200][201][202][203] (e.g., at different temperatures, chemical potentials, or just two different states on each side). One of the main results is that while the rate of exchange of entanglement entropy coincides with the thermodynamic one for free systems [201] (in analogy to homogenous cases), this is no longer the case for interacting integrable models [202]. Exact formulas, taking into account the inhomogeneities in space and time (and consequently the curved trajectories of the quasiparticles) can be explicitly written down both for free [201] and interacting [202] systems, but they are too cumbersome to be reported here. We finally stress that such an approach applies to states with locally non-zero Yang-Yang entropy, otherwise the growth of entanglement is sub-extensive and other techniques should be used [204,205].