Entanglement spreading and quasiparticle picture beyond the pair structure

The quasi-particle picture is a powerful tool to understand the entanglement spreading in many-body quantum systems after a quench. As an input, the structure of the excitations' pattern of the initial state must be provided, the common choice being pairwise-created excitations. However, several cases exile this simple assumption. In this work, we investigate weakly-interacting to free quenches in one dimension. This results in a far richer excitations' pattern where multiplets with a larger number of particles are excited. We generalize the quasi-particle ansatz to such a wide class of initial states, providing a small-coupling expansion of the Renyi entropies. Our results are in perfect agreement with iTEBD numerical simulations.


Introduction
One of the most fundamental questions in physics is how collective statistical features emerge from a microscopic deterministic time evolution, both in the case where the model at hand is classical or quantum. In particular, recent times witnessed an outburst of interest in the dynamics of closed quantum many-body systems: the rise of new experimental techniques on the one hand [1], as well as the sharpening of our understanding on the theoretical side, laid the foundations of a fruitful synergy among experimentalists and theoreticians. Indeed, the possibility of realizing almost isolated quantum systems and tuning their effective dimensionality, together with a high-precision manipulation of their interactions, gave theoreticians the perfect arena to study their favorite models. Among the various out-of-equilibrium protocols, a prominent role is covered by the quantum quench [2]: in its simplicity, it provides a neat and clear framework to ask several fundamental questions. Focusing on homogeneous sudden quenches, one imagines the system to be prepared in a given steady state of an initial Hamiltonian H 0 : common choices for the initial conditions are the ground state of H 0 or, more generally, other physically-motivated states or density matrices (e.g. thermal ensembles). Then, at time t = 0 the system is suddenly brought out-of-equilibrium by means of a sudden change of the Hamiltonian H 0 → H, which is then kept constant: after a long time, despite the non-dissipative (and thus unitary) dynamics, the system is expected to locally relax. The study of the final steady state, and how it is approached, have been at the center of an intensive research (see Refs. [3][4][5][6] for a review).
Among the various achievements, the late time relaxation of these models has been shown to be described by a Generalized Gibbs Ensemble [28] rather than the common thermal one (see Ref. [5] and reference therein). Very recently, weakly inhomogeneous and timedependent protocols have been made accessible through the Generalized Hydrodynamics [29][30][31][32]. Furthermore, spreading of both correlations and entanglement has been framed in a simple, yet powerful, quasi-particle picture [33]. In particular, how to generalise this ansatz is at the focus of our investigation.
To be concrete, let us consider a many-body quantum system, e.g. a spin chain, initialized in a certain pure state |ψ . After dividing the system in two parts A ∪ B, we may wonder how much the two are entangled: this can be measured by the von Neumann entanglement entropy [34][35][36] S A (t) = −Tr AρA (t) logρ A (t) , where we defined the reduced density matrix of the subsystem A,ρ A (t) = Tr B |ψ(t) ψ(t)|, obtained by tracing out the degrees of freedom within the B part of the time-evolved state |ψ(t) . Together with the von Neumann entropy, other common quantifiers of entanglement are the Renyi entropies, defined for arbitrary real indexes N as S (N ) Above, N can be any real number, but later on we restrict to the case of integers N . In general, entanglement entropies are not easy to compute: apart from free systems [37], progresses can be made in critical systems taking advantage of conformal invariance [38,39] with extensions to out-of-equilibrium setups [33] and recently to inhomogeneous cases [40][41][42][43][44]. In particular, bringing a low-entanglement state out-of-equilibrium will result in the entanglement growing with time: a powerful method to describe (the scaling part of) the entanglement spreading is provided by the quasi-particle picture. The scaling part of the entanglement growth is defined as it follows. Let the bipartition A∪B be identified by a set of coordinates is analogously defined). The scaling part of the entanglement is defined as S (N ) Despite not being an ab-initio calculation, the quasi-particle picture is able to exactly reproduce several features of the entanglement growth after a quantum quench. Originally introduced in critical systems [33,45,46], its applicability has been promptly extended first to generic non-interacting models [47][48][49][50][51][52][53][54], then to interacting integrable ones [55][56][57][58][59][60] (with recent developments in non-homogeneous systems [61][62][63]); finally, it has provided important insight even into not strictly-integrable systems [64,65].
The idea is based on the following semiclassical argument: let us imagine to perform a homogeneous quantum quench. At time t = 0 + the system is brought out-of-equilibrium and, from the perspective of the post-quench Hamiltonian, the system lays in a highly excited state. Semiclassically, one can depict such a state as a gas of excitations, which then start spreading across the system. If the post-quench model is a CFT or an integrable model, these excitations are stable and undergo a ballistic motion. If instead the state is evolved with a non-integrable dynamics, the quasi-particle are not stable anymore and acquire a finite lifetime. For the sake of simplicity, we consider only the case of stable excitations. The traveling particles are ultimately responsible of the entanglement (and correlation) spreading. In this regard, the understanding of the excitations' pattern is of utmost importance in order to drag quantitative predictions.
In the original formulation [33], quasi-particles are assumed to be produced in pairs of opposite momenta: this is the simplest assumption which respects the translational invariance of the protocol, and it is the common situation in most of free-to-free quenches and in exactly solvable quenches in integrable models [66]. This is, for example, the case for the protocol studied in Refs. [67,68].
For the sake of clarity (and later comparison with previous results), we briefly review free quenches in the Ising spin chain, whose Hamiltonian is given bŷ where σ x,y,z are the standard Pauli matrices and the system is assumed infinitely large. By mean of a Jordan-Wigner transformation, the Hamiltonian (4) is readily mapped into a freefermion model which is then diagonalized in the Fourier space with modes {η(k), We will thoroughly analyze the solution of the Ising model in Section 2. Initializing the system in the ground state of the Hamiltonian for a certain magnetic fieldh and quenchingh → h, one can express the pre-quench state in terms of the post-quench modes in the form of a squeezed state where N is a normalization factor, and K 2 (k, −k) is a non-trivial function ofh and h. The state |0 is the ground state (or vacuum) of the post quench HamiltonianĤ I , i.e. η(k)|0 = 0. As the form in Eq. (6) suggests, excitations are indeed independently created in pairs, their density being fixed by the mode density η † (k)η(q) = δ(k − q)n(k) Within the quasi-particle picture, backed by the assumption of the pair structure, the initial state is regarded as a source of paired excitations: pairs created at different spatial points or belonging to different momenta are uncorrelated and disentangled. Only excitations belonging to the same pair are responsible for the entanglement spreading, the total effect being additive on the different pairs. More precisely, at time t > 0 particles ballistically travel with velocity v(k) = ∂ k ω(k) and only pairs such that one excitation lays in A and the other in B contribute to the entanglement. The Renyi entropies are thus S (N ) where χ A (x) is the characteristic function of the subsystem A, namely χ A (x) = 1 if x ∈ A, otherwise zero. The factor s (N ) (k) is 1 As powerful as the quasi-particle picture is, the pair-structure is a very restrictive constraint: even though free-to-free quenches generally display this pattern, this is not true in more general setups. For example, preparing the system in the ground state of an interacting Hamiltonian (where the Wick theorem does not hold) and quenching to the free dynamics, the state cannot be in the form given by Eq. (6), which does satisfy the Wick theorem.
Hence, one is naturally led to wonder whether the quasi-particle approach can be generalized, dropping the pair-assumption in favor of a more general structure: a first step in this direction has been pursued in Refs. [69,70].
In particular, Ref. [69] considers quenches in a free model starting from a suitablyengineered initial state, where excitations were created in multiplets which shared classical correlation among the components. As such, the quasi-particle ansatz could be generalized in terms of purely-classical concepts. On the other hand, Ref. [70] considered quenches in free lattice models whose couplings are periodically modulated in space with period T , leading to a quasi-particle picture with excitations generated in pairs of momenta (k 1 , k 2 ), with k 1 + k 2 = 0 mod 2π/T , but with non-trivial quantum correlation among the different pairs. In both cases, the Wick theorem held on the initial states.
In this paper, we consider the very generic situation where the system is initialized in the ground state of a weakly-interacting Hamiltonian and then quenched to the free point: this generates a complicated pattern of excitations which exile the pair structure, as well as the studies of Refs. [69,70]. In particular, we consider initial states in the following "generalized" squeezed form [71], where multiplets (or clusters) of more than two particles are produced and the Wick theorem does not apply anymore. Above, we are assuming an underlying lattice, as such the wavevectors k are confined to a Brillouin zone [−π, π], but the extension to continuum models is straightforward.
We provide a generalization of the quasi-particle picture for the Renyi entropies of integer index N to these states with a multiplet structure in form of a systematic expansion valid for small quenches (i.e. K 2n are supposed to be small), which perfectly reproduces exact numerical simulations.
The paper is organized as follows: in Section 2 we provide a specific quench protocol which has Eq. (10) as an initial state. The quench consists in turning off a small local interaction in a spin chain, letting then the system to evolve with the Ising Hamiltonian (4). We show how to perturbatively compute the amplitudes K 2n in terms of the quench's parameter. In Section 3 we discuss our generalized quasi-particle picture and test it with iTEBD simulations, finding excellent agreement; whilst, in contrast, the naive application of the method assuming a pair structure gives incorrect results. Section 4 gathers our conclusions, and some technical details are left in the Appendices.

The model
Albeit our generalized quasi-particle picture can be applied to any initial state in the form Eq. (10), for the sake of clarity it is useful to refer to a simple model which displays such a structure. Hence, we consider an infinitely long spin chain with the following Hamiltonian where C 1 , C 2 , C 3 , C 4 are tunable parameters and γ controls the interaction strength: for the time being, we consider arbitrary values for C i . Later on, we tune these constants in order to enhance the role of the multiparticle clusters in Eq. (10) when compared with the pairs. We assume the system to be initialized in the ground state for γ = 0, then at t = 0 we switch the interaction off γ = 0 and let the system evolve with the Ising Hamiltonian. For C 1 = C 4 = 0, the protocol reduces to a transverse-field quench in the Ising chain: this free-to-free quench has been already studied in Ref. [68]. Hence, the pair structure holds in this case. On the contrary, if C 1 and C 4 are not zero, the pre-quench Hamiltonian is truly interacting: the ground state does not satisfy the Wick theorem any longer, and thus the pair structure for the post quench excitations is inevitably spoiled. Before of discussing the structure of the interacting ground state, we briefly review the diagonalization of the Ising Hamiltonian (4) which governs the post-quench dynamics. First, the spin variables can be written in terms of standard fermionic operators through the following Jordan-Wigner (JW) transformation with σ ± = (σ x ± iσ y )/2. The fermions satisfy standard anticommutation rules {c j , c † j } = δ j,j . After the JW transformation, the Ising Hamiltonian is rewritten in the following form (we drop additive constants which do not affect the dynamics) In terms of the fermions, the Ising Hamiltonian is quadratic and promptly diagonalized in the Fourier space, once the proper fermionic modes η(k) are defined where the Bogoliubov angle Ω k is and ω(k) = (cos k − h) 2 + sin 2 k are the single-particle energies of the modes entering in the diagonal Ising Hamiltonian (5). Since ω(k) > 0, the ground state can be identified with the vacuum of the theory |0 , namely the state annihilated by all the modes η(k)|0 = 0. We now move to discuss the pre-quench Hamiltonian (11). We use the JW transformation to recast it in terms of the fermionic degrees of freedom, and then rewrite it in terms of the diagonal mode-operators of the non-interacting Ising Hamiltonian. This tedious, albeit simple, calculation results in the following general form where, the summation is over all the positive integers n 1 and n 2 , and the coefficients H n 1 ,n 2 are complicated functions of the pre-quench parameters C i and of the magnetic field h. The δ−function in the momentum space is a consequence of the translational symmetry of the problem. The coefficients H n 1 ,n 2 must satisfy certain constraints discussed hereafter. From the anti-commutation relations of the fermions we have that H n 1 ,n 2 (k 1 , ..., k n 1 |q 1 , ..., q n 2 ) is an antisymmetric function of the k and q momenta separately. The hermicity of the Hamiltonian (16) implies H n 1 ,n 2 (k 1 , ..., k n 1 |q 1 , ..., q n 2 ) = H * n 2 ,n 1 (q n 2 , ..., q 1 |k n 1 , ..., k 1 ) .
Moreover, for parity reasons, the coefficient H n 1 ,n 2 vanishes whenever n 1 +n 2 is odd. The task of computing the coefficients H n 1 ,n 2 in terms of the interactions in Eq. (11) is straightforward. For example, in the simplest case of a quench among two Ising spin chains, namely posing C 1 = C 3 = C 4 = 0 and C 2 = 1, one obtains and all the other coefficients are zero. In the more general case with C 1 , C 3 and C 4 not zero, also the coefficients H n 1 ,n 2 with n 1 + n 2 = 4 are present, while all the terms with n 1 + n 2 > 4 vanish.
We will explicitly write the expression for H n 1 ,n 2 for a specific choice of the parameters C i after having discussed the structure of the ground state in the forthcoming section.

Multiplets structure of the pre-quench state
There are several reasons that point to the structure in Eq. (10) for the pre-quench state in terms of the post-quench excitations, all of them extensively discussed in Ref. [71]. Such a state satisfies, at any order in a small K expansion, basic thermodynamic properties that ground states of a local Hamiltonian must have; namely, the extensive behavior of thermodynamic quantities and the cluster property for local observables [71]. In addition, the form in Eq. (10) can be recovered from standard perturbation theory. We thus assume the ground state to be in this form and discuss how to determine the coefficients K n starting from the Hamiltonian in Eq. (16): we follow the same method presented in Ref. [71], namely we impose the eigenvalue equationĤ|ψ = E G |ψ , together with the ansatz in Eq. (10). In this respect, one should notice that applying an annihilation operator η(p) to the state (10) is equivalent to take a functional derivative of the exponential, namely Above, while taking the functional derivatives, the anticommutation relations of the fermions impose that the field η † must be regarded as a Grassmanian variable, leading to extra signs. For example, one has Therefore, replacing the annihilation operators in Eq. (16) with the functional derivatives, and applying the Hamiltonian on |ψ , one reaches the following set of equations where I 2n [K] are non-trivial functions of K 2n which couple multiplets with different indexes and vanish if one sets K = 0. The computation of I 2n can be framed into a Feynman diagrams picture [71] which, due to its technical aspect, we leave to Appendix A. The condition (21) can be satisfied if and only if the infinite set of equations is imposed. In general, even in the case where a finite number of H n 1 ,n 2 appears in the Hamiltonian (16), these equations cannot be closed with a finite number of K 2n wavefunctions. The free-to-free case is exceptional and the equations can be closed with only K 2 being not zero (as it should be): in Appendix A we analyze this case and obtain the exact K 2 amplitude, which can be computed by mean of other methods, thus providing a consistency check.
Here, we rather approach the problem in the small interaction limit: within this assumption, Eq. (22) is extremely useful in developing a γ−power expansion for the K−amplitudes by means of an iterative solution.
This perturbative expansion differs from the standard ground-state perturbation theory [72]. Even retaining only up to a finite order in the γ−expansion of the functions K 2n , contributes to any order in γ in the expansion of the state |ψ are present, because of the series expansion of the exponential in Eq. (10). In other words, one can look at the γ−expansion of the functions K 2n as a partial re-summation of the standard perturbation theory on the non-interacting ground state. As already mentioned, among the advantages of the present approach there is the fact that at any perturbative order the state retains some features that a proper ground state should have [71]; namely the extensive behavior of thermodynamics quantities, and the cluster property of the correlation functions. From Eq. (22), one gets Hence, the coefficients H 2n,0 determine, at first order in γ, the multiplets' population. As we previously commented, setting C 1 = C 3 = C 4 = 0 in the Hamiltonian (11) results in a squeezed state (6): consistently, H 2n =2 = 0 and Eq. (23) predicts a non-zero amplitude K 2 , while all other terms vanish. On the other hand, the coefficients C i can be suitably engineered to suppress the pair production and enhance other multiplets.
In particular, feeding the mode decomposition (14) into the general Hamiltonian (11) (after it has been expressed in the fermionic basis through the Jordan-Wigner transformation) and rearraging the whole expression in the form (16), one can impose H 2,0 = 0 choosing Above, ... 0 is the expectation values of the fermionic operators on the non-interacting ground state, namely In this case, all the H 2n =4,0 coefficients vanish, while Therefore, this quench creates only multiplets of 4 particles (at first order in γ). Above, the sum is over all the possible permutations P of four elements.
3 The generalized quasi-particle picture for the entanglement growth After having discussed the structure of the pre-quench state in terms of the post-quench excitations, we can now generalize the quasi-particle picture to include the presence of multiplets beyond the simpler pairs' case. It should be stressed that the quasi-particle picture is an ansatz: its formulation is built on reasonable assumptions and on the experience gained from the previous literature, which also provides an important consistency check. However, its correctness must be ultimately confirmed by numerical simulations, which turn out to be in perfect agreement with the proposed ansatz.
Similarly to the case with paired excitations, we look at the initial state as a homogeneous source of quasi-particles (but generalizations to inhomogeneous cases are straightforward [73]), which organize in multiplets as Eq. (10) suggests. In this respect, following Refs. [70,73], we introduce a set of auxiliary Hilbert spaces H x , where x labels the position (in the perspective of a coarse grain analysis, we can replace the underlying lattice with a continuum). Similarly, we introduce auxiliary fermionic mode operators {η x (k), η † x (q)} = δ(k − q) and a local vacuum |0 x , then we promote Eq. (10) to a state |ψ x in the local Hilbert space H x Above, the amplitudes K 2n are those of the homogeneous model. The quasi-particle ansatz states that excitations created at different points in space are uncorrelated and disentangled, as such they contribute additively to the entanglement: using of the auxiliary Hilbert spaces, we write S (N ) where s

(N )
Ax,t is the entanglement density contribution coming from each |ψ x state. We are considering bipartitions of the system A ∪ B in the real space, but the local Hilbert spaces are constructed in terms of creation-annihilation operators with the momentum k as quantum number. The correct dictionary to translate the bipartion A ∪ B into a proper bipartition in the momentum space of the local auxiliary Hilbert space relies on the ballistic propagation of the quasi-particles and B x,t is analogously defined. In Fig. 1 we provide a semiclassical representation of Eq. (32). The contribution to the entanglement entropy of each spatial point is defined tracing on the local degrees of freedom. In this respect, we define the reduced density matrix and then the entanglement entropy. While doing so, a subtlety arises: indeed, the entanglement entropy coming from Eq. (33) is formally divergent (as we discuss in Appendix B) and its computation requires a proper finite-volume regularization with playing the role of a large (but finite) volume. Then, one defines the entanglement density s In order for the quasi-particle picture to be predictive, one is left with the highly challenging task of computing s (N ) Ax,t . In general, the structure of Eq. (30) when other multiplets besides pairs are present is too much complicated to drag an exact result. However, in the limit of small quenches (i.e. K 2n small), one can revert to a perturbative approach and expand the Renyi entropies for integer indexes N in powers of the amplitudes K 2n . We present the full technique, based on Feynman diagrams, in Appendix B: hereafter, we provide and discuss the first non-trivial order and compare with ab-initio numerical simulations of the model (11). For s  Figure 1: Pictorial representation of the bipartition induced in the auxiliary Hilbert space H x . On the horizontal axis we consider the momenta within the first Brillouin zone, on the vertical axis we consider the function x + tv(k). We focus on two different times. We consider a bipartition in the real space A ∪ B with A an interval in the vertical axis, which is then prolongated in the whole red-shaded area. The domain A x,t (blue-shaded area) in the momentum space is determined by those momenta such that x + tv(k) lays within the interval A.
As a simple consistency check, one can consider the case where only pairs are present and compare with the usual quasi-particle ansatz (8), using the mode-density (7) for small K 2 : of course, the two calculations are in complete agreement. As we have already stressed, Eq. (36) holds for integers N , but one could wonder if the result can be analytically continued to real indexes and, through the limit N → 1, reaching the von Neumann entanglement entropy. Unfortunately, this is not the case: dropping further corrections in K and promoting N to be a real variable, the limit N → 1 of Eq. (36) diverges. On the other hand, lim N →1 s

(N )
Ax,t is expected to exists and to give the quasi-particle picture for the von Neumann entropy. This discrepancy is due to the fact that the operation of truncating the K−expansion does not commute with the analytic continuation: there could be terms, which we are neglecting in the Renyi entropies, such that after a proper analytic continuation become of the same order of, or even more relevant than, the terms written explicitly in Eq. (36). Now, we specialize the discussion to the model of Section 2 and compare the quasiparticle prediction with iTEBD simulations. For general values of the couplings C i in Eq. (11), both pairs and quartuplets are generated at the first order in the interaction (23): higher clusters of particles are created as well, but at further orders in the perturbation theory. Hence, Eq. (36) predicts an entanglement growth ∝ γ 2 (with further corrections O(γ 3 )) where pairs and quartuplets contribute additively. In order to enhance the role of the quartuplets, we now consider the choice of parameters Eq. (24) which, at the first order in γ, does not excite pairs.
We first focus on a bipartition of the system where A = [0, ∞]: in this case, the integrals and summations appearing in Eq. (36) can be greatly simplified. Indeed, considering a multiplet of quasiparticles {k i } 4 i=1 which ballistically propagates for a time t, it can contribute to the entanglement if and only if the excitations were created at position x ∈ (t min k i v(k i ), t max k i v(k i )). This causes the entanglement to linearly grow with time. More precisely, one gets the following scaling form with F (h) = 1 2π In Fig. 2 we provide a numerical check of the scaling form for different Renyi entropies, together with a comparison of the slope F (h) computed as per Eq. (38) with the same quantity extracted from iTEBD results. In order to stress the importance of considering the multiplets, we plot the slope F (h) one would find expanding Eq. (9) and Eq. (8) at the same order in γ 2 , using for n(k) the mode density computed on the prequench state: The two curves are clearly apart and the numerical data are in perfect agreement with the multiplet result. In Fig. 3, instead, we consider a finite interval A = [0, L]: in this case, the quasi-particle picture can be cast in a scaling form similar to Eq. (37) Eq. (41) is exact: checking its consistency with Eq. (40) using Eq. (39) is a straightforward calculation.
Accessing the scaling part of the entanglement entropies for finite intervals by mean of iTEBD methods is difficult: in Fig. 3 we provide a plot of F(tL −1 , h) for different values of the magnetic field. Again, on top of the correct multiplet result, we plot what one would get employing the mode density n(k) into the ansatz based on the pair structure. The curves are clearly far apart, confirming once again the importance of considering the multiparticle structure in the quasi-particle ansatz Eqs. (8)(9).
It should be stressed that, for the current choice of the C i parameters, the quasiparticle ansatz based on the pair structure is not justified at all, since only quartuplets are present. For a generic choice of the couplings C i , also paired excitations will be present and the standard ansatz is expected to improve. However, for the sake of dealing with more compact expressions, we consider only the choice Eqs. (24) (25).

Conclusions
In this paper we investigated the applicability of the quasi-particle ansatz for the entanglement growth in those frameworks exiling the commonly-assumed pairwise structure of the excitations' pattern. Indeed, there are several quenches exhibiting a far richer structure than the simplest paired one: this is verified, for example, in any quench from the ground state of an interacting model (on which the Wick theorem does not hold) to a free one, where the paired-excitation pattern implies the Wick theorem. As a specific example, we considered a sudden homogeneous quench from a weakly interacting spin chain to a free one, unveiling the rich structure of the initial state in terms of the post-quench modes, where excitations' multiplets beyond simple pairs are present. We generalize the quasi-particle picture to this wide class of states, capturing the growth (in the scaling limit) of the Renyi entropies of integer index N . Albeit the complicated structure of the initial state makes difficult to exactly determine the input of the quasi-particle ansatz, we provided a systematic expansion for small interactions which is in excellent agreement with exact iTEBD simulations.
Our findings confirm the wide applicability of the quasi-particle ansatz to our understanding of entanglement spreading, pointing out at the same time the urge of going beyond the simple pair-structure.
Several interesting developments are left for future investigations. First of all, the quasiparticle picture has been proven to be a precious tool in investigating the correlation functions as well [68]: it would be interesting to explore the consequences of higher multiplets in this analysis.
While we accessed the Renyi entropies for integers N , the von Neumann entanglement entropy still remains an elusive quantity, the reason being the impossibility of performing an analytic continuation N → 1 term by term in the small coupling expansion. In this perspective, partial resummations of the perturbative series could lead to meaningful analytic continuations: a better understanding of the pattern of the diagrammatic representation will be surely helpful in this perspective.
Very recently, an ab-initio technique to address entanglement entropies, based on formfactors of twist fields, has been proposed in Ref. [74] and successfully applied to the Ising field theory, starting from an initial state with pair-structure. It would be of great interest applying those techniques to the more general form of states considered in this paper.
Lastly, a surely compelling question concerns quenches towards truly interacting integrable models: as clearly shown by the form-factor perturbation theory of the pre-quench state [75], quenching towards an interacting integrable model will generally result in a much more complicated structure than pairwise excitations. The lack of a pair-structure in generic interaction quenches in integrable models have been pointed out in Ref. [76], with the further development of a perturbative study in terms of the pre-quench dynamics [76][77][78]. However, the question whether states in the form Eq. (10), together with the entire quasi-particle picture, can be generalized to truly interacting integrable case is still an open problem, which we leave for future investigations.

Acknowledgements
We thank P. Calabrese and V. Alba for useful comments on the manuscript.

Funding information This work has been supported by the European Research Council under the ERC Advanced grant 743032 DYNAMINT.
A The diagrammatic expansion for the pre-quench state In this appendix we discuss a systematic description of Eq. (22) in terms of Feynman diagrams, whose iterative solution ultimately leads to a power-expansion of the K wavefunction in terms of the small parameter γ. This representation was already discussed in Ref. [71] for the bosonic case and can be promptly generalized to the fermionic one. First of all, we represent the wavefunction K 2n with a dot with 2n departing arrows, each of them representing a momentum of the wavefunction. We use a similar representation for the interaction vertexes H n,m , assigning them empty circles with n departing and m incoming arrows, representing the momenta in the vertex.
The value of I 2n in Eq. (22) is then computed according with the rules below. As a guide for better understanding the Feynman rules, in Fig. 4 we provide a specific example of a possible diagram and its value.
1. Draw a single vertex associated with H n ,m and a few vertexes associated with K 2n .
Then, contract the incoming m arrows of the vertex H n ,m with the arrows departing from the K−vertexes. All the in-going arrows must be contracted, the resulting diagram must not have disconnected parts and there must be exactly n out-going legs.
2. At any vertex one must impose the conservation of momenta, according with the direction of the arrows. A global Dirac−δ for the conservation of the n−outgoing momenta always appears and must be removed (as it should be clear from Eq. (21), where the global δ−function is explicitly factorized).
3. Extra caution must be payed because of the fermionic nature of the operators, giving the correct sign to the diagram. This can be done as follows: write the vertex H n ,m and, after that, the product of the wavefunctions appearing in the diagram. Then, as per the Wick theorem, add the sign of the permutation needed to reshuffle the position of the momenta in order to contract those that are equal and put the out-going one into the right order.
4. One must integrate over the momenta carried by the incoming legs and divide by a symmetry factor equal to the number of permutations of legs and vertexes that leave the diagram the same.
5. Finally, sum over all the possible diagrams with n−outgoing legs.
In general, if the interaction vertexes H n,n with n + n > 2 are present, the equations (22) cannot be closed with a finite number of wavefunctions K: an important exception is the case when only vertexes n + n = 2 are present. In this case, we are considering a quench among quadratic, i.e. free, models where the squeezed state form Eq. (6) must be recovered. It is instructive to approach the problem within this formalism, providing in the meanwhile a consistency check: let us assume only the interactions H 2,0 , H 1,1 , H 0,2 are present. In Fig.  5 we draw the diagrams relevant for the case n = 1 of Eq. (22), where we already assumed Figure 4: Above, we provide an example of a Feynman diagram contributing to I 6 (k 1 , ..., k 6 ), together with its value. The symmetry factor in front of the integral comes from the possible permutations of legs: more specifically, the factorials 3! and 2! come from the permutation of the triplet and pair of external legs, while an additional term 2! is due to the permutations of the internal legs.
K 2n =2 = 0. These graphs lead to the following non-linear equation for K 2 which has two possible solutions The correct choice is the minus sign, which ensures that in the γ → 0 limit the amplitude K 2 vanishes. Moreover, plugging the two solutions into Eq. (22) and looking at the case n = 0, one can check that choosing the minus sign in Eq. (43) corresponds to the minimum energy E G . Using in the above equation the amplitudes associated with quenches of the magnetic field in the Ising model (18), one promptly recovers the known result obtained by mean of Bogoliubov rotations [67].

B The diagrammatic expansion for the entanglement spreading
In this appendix we present the systematic expansion of the quasi-particle ansatz which has, as a first non trivial order, the result Eq. (36). We will proceed through a direct computation of the reduced density matrix and powers thereof, applying the definition Eq. (35): as we already commented, states in the form (30) seem unfeasible of exact computations, but an expansion valid for small K can be straightforwardly performed through Feynman diagrams, as we are going to discuss. Rather than considering the state Eq. (30), it is more convenient to use the unnormalized version In this way, the norm appearing in Eq. (30) is N = Ψ x |Ψ x . As a warm up, we start presenting the Feynman rules for computing the norm of the state, then generalize them to the partial traces and powers of the reduced density matrices. Expanding the exponential in Eq. (44) and contracting with the bra, one contracts the momenta in the wavefunctions K 2n coming from the ket with the conjugated K * 2n coming from the bra. With the same notation of the previous appendix, we represent the K 2n wavefunctions with vertexes with 2n outgoing arrows. On the other hand, the wavefunctions K * 2n are represented as vertexes with 2n ingoing arrows: performing the contractions, only legs with compatible orientation can be joined. While computing the norm N , all the legs must be contracted: one should then integrate over the momenta associated with internal legs enforcing momentum conservation at each vertex. Then, one must divide for a symmetry factor equal to the number of permutations of legs and vertexes that leave the diagram unchanged. Lastly, the sum over all the possible diagrams must be considered: in this respect, as it is standard in perturbative expansions of partition functions, the sum over the Feynman diagrams can be represented as an exponential of the sum over the connected diagrams N = Feynman diagrams = exp connected Feynman diagrams .
In writing this expression, a subtlety arises. Indeed, in any connected Feynman diagram the global momentum conservation is always satisfied, resulting in a formally divergent contribution δ(0), where the Dirac-δ acts in the momentum space. As it is standard in these computations, one can make sense of this divergence by mean of a finite-volume regularization (34) of the system, which would lead to the replacement δ(0) → /(2π), with the system's size. In this perspective, the logarithm of the norm behaves extensively in the system's size log Ψ x |Ψ x ∝ : this is in agreement with the fact that Ψ x |Ψ x can be interpreted as a partition function [71] and hence its logarithm, namely a free energy, is extensive. We can now generalize the method to compute reduced density matrices, then powers thereof. Again, we work with the unnormalized state Eq. (44) and define the unnormalized reduced density matrixσ Ax,t asσ Note that we clearly have Tr Ax,tσAx,t = N ,therefore the normalized density matrix isρ Ax,t = σ Ax,t /N and the Renyi entropies are s (N ) The Feynman rules to compute Eq. (46) are exactly the same of those for the norm N , with two exceptions: i) we can now have free (namely, uncontracted) external legs representing the degrees of freedom of A x,t (on which we are not tracing yet) and ii) the momenta associated with the contracted legs on which we are integrating are confined to the domain B x,t , rather than to the whole Brillouin zone. Again, the sum over the Feynman diagrams can be exponentiated in terms of the connected diagrams, as in Eq. (45): this exponential of diagrams represent the reduced density matrix Eq. (46). Armed with this machinery, we can now discuss powers of the reduced density matrix. Referring to Fig. 6, we divide the plane into different sectors by mean of vertical dashed lines, each sector represents a reduced density matrix. In particular, while computingσ N Ax,t , one needs N + 1 vertical lines, splitting the plane in N parts). The diagrams drawn in each region belong to the respective density matrix in the product, which is then performed contracting the free legs with the following rules. Legs can be connected only between adiacent regions and outgoing arrows can be contracted only towards the left, while ingoing arrows only towards the right: this because in(out) arrows are respectively associated with bra(ket) states. The internal legs crossing the vertical lines represent contraction among different density matrices: therefore, their momenta are integrated over the A x,t domain. After this operation, all the legs are contracted, with the exceptions of the out-going legs in the leftmost region and ingoing arrows in the rightmost one. Now, we can finally take the trace of the product joining together the legs in the leftmost and rightmost regions, their momenta being still integrated over the A x,t domain. As usual, momentum conservation must be imposed at any vertex and the proper symmetry factors kept in account.
Finally, the sum over all the possible diagrams must be considered. Again, the sum over the Feynman diagrams is equivalent to the exponential of the sum over the connected diagrams.
Tr Ax,tσ N Ax,t = exp connected Feynman diagrams , and, as we have already discussed while considering the norm N , a global momentum conservation results in a divergent term δ(0), which is regularized with a finite volume scheme δ(0) → /(2π). Hence, log Tr Ax,tσ N Ax,t ∝ and the density of entanglement s

C Numerical simulation of time evolution
The numerical simulations of the entanglement dynamics has been obtained by using the infinite volume Time-Evolving Block-Decimation (iTEBD) algorithm [79]. The algorithm represents a generic translational invariant state in one dimension as |Ψ = ...,s j ,s j+1 ,...
s j+1 e · · · | . . . , s j , s j+1 , . . . , where s j spans the local spin-1/2 Hilbert space, Γ s o/e are χ × χ matrices associated with the odd/even lattice site; Λ o/e are diagonal χ × χ matrices containing the singular values corresponding to the bipartition of the system at the odd/even bond.
Exploiting the representation (51), the Renyi entropies of the semi-infinite subsystem are given by  where a bipartition at the odd/even bond can be indifferently chosen.
The infinite Matrix Product State (iMPS) representation of the ground state |Ψ GS of the pre-quench XY Z Hamiltonian has been obtained by (imaginary)time-evolving the following initial product state We used a second-order Suzuki-Trotter decomposition of the evolution operator with imaginary time Trotter step τ = 10 −5 . The pre-quench parameters {C 1 , C 2 , C 3 , C 4 } have been tuned accordingly to Eqs. (24) (25), with γ ∈ {0.1, 0.08, 0.06} and transverse field h ∈ {0.1, 0.2, 0.3, . . . , 1.9}. Thanks to the energy gap between the ground state and the the rest of the spectrum, the iMPS ansatz is very accurate (up to machine precision) by retaining an auxiliary dimension χ 0 = 16 or 32, depending on the vicinity of the critical point (i.e. γ = 0, h = 1). Similarly, the post-quench time evolution has been obtained by evolving the corresponding |Ψ GS with the Ising Hamiltonian with γ = 0. For this purpose, a second-order Suzuki-Trotter decomposition of the evolution operator was used again, with real time Trotter step dt = 10 −3 . In order to keep the truncation error as small as possible (≤ 10 −6 ), the auxiliary dimension was allowed to grow up to χ M AX = 1024 which was sufficient to reach a maximum time T = 100. We were able to reach such large times thanks to the relatively small quenches considered. As explained in the main text, the entanglement entropy, for all quenches under investigation, growths linearly with a slope proportional to γ 2 , thus remaining very small for all the duration of the simulations.
The numerical data for the Renyi entropies as a function of the rescaled time γ 2 t showed a linear increase (apart from oscillations) whose slope depends on the particular value of the transverse field h. In particular, a numerical estimation of the entanglement entropy slope ∂ γ 2 t S (N ) [0,∞] has been obtained by performing a linear fit of the iTEBD data in the entire timewindow 0 ≤ t ≤ 100, for the three different values of γ we have considered. Thereafter, the γ → 0 value has been evaluated by linear extrapolation.