Spreading of entanglement and correlations after a quench with intertwined quasiparticles

We extend the semiclassical picture for the spreading of entanglement and correlations to quantum quenches with several species of quasiparticles that have non-trivial pair correlations in momentum space. These pair correlations are, for example, relevant in inhomogeneous lattice models with a periodically-modulated Hamiltonian parameter. We provide explicit predictions for the spreading of the entanglement entropy in the space-time scaling limit. We also predict the time evolution of one- and two-point functions of the order parameter for quenches within the ordered phase. We test all our predictions against exact numerical results for quenches in the Ising chain with a modulated transverse field and we find perfect agreement.


I. INTRODUCTION
During the last decade, the study of the non-equilibrium dynamics after a quantum quench (i.e. after an abrupt change of a parameter in a quantum Hamiltonian) has been the subject of intense theoretical and experimental investigations, see e.g. Refs. 1-4 as reviews on the subject. One of the main issues concerned the nature of the stationary state that describes local properties of the system. Nowadays we have a rather clear understanding of this stationary state: a generic system for long time attains a thermal state 5-10 while an integrable model relaxes to a generalised Gibbs ensemble [11][12][13][14][15][16][17] (also many-body localised systems have very peculiar non-equilibrium features [18][19][20]. Conversely, the approach to the stationary state and the exact time evolution of physical observables remain less generally understood problems, in spite of a very intense activity. While some approximate numerical and analytical methods to tackle the problem in interacting integrable models exist (see, e.g., Refs. [21][22][23], exact analytical results are scarse even for free systems: only few first principle calculations have been worked out up to a final analytic form [24][25][26][27][28][29][30][31][32] . In this respect, the quasiparticle picture 33,34 proved to be an extremely valuable tool. Although it is not an ab-initio technique, it provides a qualitative and quantitive understanding of the time evolution of some observables under specific conditions. It has been introduced to explain the entanglement evolution in the scaling limit after a quantum quench 33 and originally tested against the exact results in conformal field theories [33][34][35] , in free models 25,33,34,[36][37][38][39][40][41][42][43][44][45] , and against many numerical simulations [46][47][48][49][50][51][52] . Only very recently these concepts have been used to quantitatively predict the time evolution of the entanglement entropy in generic interacting integrable systems [53][54][55] . In the field theoretical context, it has also been shown that the quasiparticle picture can be used to understand the time evolution of the one-and two-point functions of the order parameter 56 (more generically of correlations of primary operators whose expectation value is non-vanishing in the initial states 56 ). Some results in the very few analytically treatable free models are compatible with this picture 26 , but the general regime of applicability of these ideas to generic operators is not clear [57][58][59] , even within the realm of quadratic models.
A central object of this paper is the entanglement entropy 60 where ρ A ≡ TrĀ|ψ ψ| is the reduced density matrix of a subsystem A (havingĀ as complement) of a system in a pure state |ψ . Its time evolution plays a crucial role in the understanding of the non-equilibrium dynamics of isolated quantum systems. Indeed, the growth of the entanglement entropy in time has been related to the efficiency of tensor network algorithms 61-65 such as the time dependent density matrix renormalisation group. Furthermore, the extensive value (in subsystem size) reached by the entanglement entropy at long time has been understood as the thermodynamic entropy of the ensemble describing stationary local properties of the system [53][54][55][66][67][68][69][70][71] . Within the quasiparticle description, the initial state is regarded as a source of entangled quasiparticles which ballistically propagate across the system and carry entanglement. The structure of the pre-quench state in terms of the post-quench excitations is essential in dragging quantitative predictions for the spreading of entanglement. For example, in Ref. 25 the XY spin-chain has been considered and quenched in the magnetic field h t<0 → h, starting from the ground state at h t<0 . The XY model is diagonalised in terms of spinless fermions through a Jordan-Wigner transformation: in the thermodynamic limit N → ∞ the momentum is continuous, thus the fermions obey {η(k), η † (q)} = δ(k − q), and The prequench ground state, identified with the vacuum of the prequench modes |0 h t<0 , is readily written in terms of the post quench vacuum |0 h in the form of a squeezed state with K a non trivial function of h and h t<0 , its specific form being irrelevant for our purposes. Squeezed states are common in quenches in free theories, due to the fact that the pre and post quench modes are usually connected through a Bogoliubov rotation. The form of Eq. (4) is rather appealing: the quench modes are excited in pairs of opposite momenta, distinct pairs being created independently. In this case, the entanglement growth is well described by the following quasiparticle picture 33 .
The initial state is regarded as a source of quasiparticles, homogeneously distributed in space. After the quench, each particle ballistically propagates with velocity v(k) = ∂ k E(k). Pairs originating at different positions or with different momentum are unentangled, only particles belonging to the same pair of momentum (k, −k) and originating in the same position are entangled. Given the partition of the system A ∪Ā and considering a time t, only pairs such that one quasiparticle belongs to A and the other toĀ contribute to the entanglement, their contribution being additive. In the case where A is chosen to be an interval of length , the entanglement entropy, in the space-time scaling limit t, → ∞ with t/ fixed, has the following scaling form where s(k) is the contribution to the entanglement associated with each pair. In the standard homogeneous situation, the weight s(k) can be fixed by requiring that for t → ∞ the entanglement entropy density matches the thermodynamic one of the post quench steady state 53-55 being n(k) = |K(k)| 2 /(1 + |K(k)| 2 ) the density of excitations of momentum k, i.e. η † (k)η(q) = δ(k − q)n(k). The generalisation to several, uncorrelated, particle species is obvious, one has just to sum over the contributions of each independent species. This is a simple further step which allows us to describe quenches in truly interacting models where several particle species may be present. Indeed, in Ref. 53 quenches in the XXZ spin-chain have been considered, and the quasiparticle picture (5) (through a suitable dressing of the velocity v(k) and the weight s(k)) has been found to be correct. Interestingly, also Eq. (6) has a very simple semiclassical interpretation. Indeed, focussing on a given momentum k, s(k) is just the entropy (i.e. the logarithm of the number of equivalent microstates) of a mode which is occupied with probability n(k) and empty with probability 1 − n(k).
The main physical feature for the entanglement evolution captured by Eq. (5) is the so-called light-cone spreading, i.e. a linear increase at short time followed by saturation to an extensive value in . Indeed, if a maximum velocity v max for the quasiparticles exists, then since |v(k)| ≤ v max , the second term vanishes when 2v max t < , and the first integral is over all positive momenta, so that S A (t) is strictly proportional to t. On the other hand as t → ∞, the first term is negligible and S A (∞) is proportional to . Actually, analytical, numerical and experimental results 48,69,[71][72][73][74][75][76] suggest that such a light-cone spreading is more generically valid that what suggested by the quasiparticle picture.
At this point, it must be clear that the physical assumptions behind validity of Eq. (5) for the entanglement evolution is that quasiparticles must be produced uniformly in space and in uncorrelated pairs of opposite momenta. Given the large success of Eq. (5) in describing the entanglement evolution in the scaling regime, a lot of recent activity has been devoted to understand how Eq. (5) gets modified when some of these assumptions are weakened. For example, the quasiparticle picture has been extended to large-scale inhomogeneous setups 45 in which the initial state is regarded as a non uniform source of quasiparticles which ballistically propagate for t > 0. While the final expression (5) is slightly modified in order to keep in account the initial inhomogeneity, the core of the result is still Eq. (6), where n(k) is promoted to have a weak spatial dependence and the semiclassical interpretation of the entanglement weight still holds true.
Another setup in which the semiclassical interpretation can be applied, but which lays outside the usual framework of uncorrelated pairs, has been studied in Ref. 42 in a free-fermion model. In that case, the homogeneous initial state was populated with excitations of n different species, with a constraint on the sum of the excitation densities n i=1 n i (k) = 1 which ultimately introduces non trivial correlations. However, this constraint is classical in nature and it does not spoil the semiclassical interpretation of the entanglement entropy. As a matter of facts, the constraint changes the form of Eq. (6), which nevertheless can still be viewed as the entropy of fermions obeying the extra condition. In particular, the entanglement growth is fully determined in terms of the excitation densities {n i (k)} n i=1 and of the velocities of each species {v i (k)} n i=1 with no other information required. In this work we investigate those situations where the initial state is populated by several quasiparticle species, which are non trivially quantum-correlated. The presence of true quantum correlation among the excitations forces us to dismiss the simple entanglement weight (6) together with its classical interpretation. However, despite this lack of classicality, the quasiparticle paradigm will still hold true and the entanglement growth (in the scaling limit) is fully determined in terms of ballistically-propagating localised excitations.
In particular, we are interested in free Hamiltonians possessing n species of excitations (assumed to be fermionic for concreteness) Each mode has its own group velocity v i (k) = ∂ k E i (k) and we assume the existence of a single Brillouin zone [−B, B].
The initial state |Ψ is taken to be a non trivial generalisation of the single species squeezed state (4), i.e.
where |0 is the vacuum η i (k) |0 = 0. This class of states is gaussian, i.e. the knowledge of all the correlation functions can be reduced, by mean of a repeated use of the Wick Theorem, to the two-point correlators where the correlation matrices C (1) (k), C (2) (k) have dimension n × n and are functions of M(k), the specific relation being inessential for our purposes. The contact point with the previous literature can be made in the case where C (1) (k) and C (2) (k) are diagonal on the particle species. States in the form (8) are not rare, making our generalisation more than a mere academic question. For example, we can revert to suitable inhomogeneous quenches in free models in order to realise states such as (8), as we explain thereafter.
Our quasiparticle ansatz will be formulated in full generality without looking at a precise model, however we ultimately rely on the periodically-modulated inhomogeneous Ising chain as a convenient benchmark where we require h j to be periodic with period n The periodicity of the magnetic field effectively splits the original lattice into n sublattices coupled in a non trivial way, each one with lattice spacing n so that the Hamiltonian (10) can be then diagonalised in terms of n particle species {η i (k)} n i=1 . Starting in the ground state for a given set of magnetic fields and quenching towards different {h j } n j=1 creates initial states in the form (8). Of course, in the scaling region where the quasiparticle description holds true, the lengthscale of the inhomogeneity is negligible: such a quench can be regarded as being homogeneous with several species of quasiparticles.
This work is organised as it follows. in Section II we present a general discussion of our quasiparticle ansatz for states in the form (8). In Section III we benchmark our predictions in the inhomogeneous Ising model, providing the details of its solution. In Section IV we provide a quasiparticle description for the time evolution of the correlators of the order parameter in the inhomogeneous Ising model, thus generalising the results of Ref. 26. In Section V we gather our conclusions. Two appendices support the main text with some technical details.

II. THE QUASIPARTICLE PREDICTION FOR THE ENTANGLEMENT ENTROPY SPREADING
In this section we present our result for the time evolution of the entanglement entropy. Of course, being a quasiparticle prediction, its derivation is not rigorous, but relies on reasonable physical arguments and on the experience gained from the existing literature. Our arguments are somehow related to those of Ref. 45, where the standard quasiparticle picture for a single species and pair excitations has been extended to weakly inhomogeneous non-equilibrium protocols. For definiteness, we focus on a a bipartition A ∪Ā where A is an interval of length . Following the standard quasiparticle picture, we regard the initial state as a source of excitations such that i) quasiparticles generated at different spatial points are not entangled, ii) quasiparticles associated with pairs of different momentum are not entangled. As usual 33,45 , we further assume the contribution to the entanglement of unentangled pairs to be additive Here, s A (x, k, t) is the contribution to the entanglement at time t given by the quasiparticles originated at t = 0 in position x and with momentum ±k (the momentum integration in (12) runs on positive values in order to avoid double counting). The difficulty, as well as the main result of our investigation, is finding the correct ansatz for s A (x, k, t): we propose to construct s A (x, k, t) out of suitable finite-dimensional ancillary Hilbert spaces and partitions thereof.
To each position x and momentum k we associate a Hilbert space constructed as a Fock space starting from a vacuum such a Hilbert space has dimension 2 2n . We recall that n is the number of different species of quasiparticles. We have in mind the following, suggestive, correspondence In order to make clearer and precise such a statement, we consider a state in the ancillary Hilbert space encoded in a density matrix ρ x,k such that: i) it is gaussian in the fermions f x,k i (i.e. the Wick Theorem holds true); ii) its correlators are the same of the corresponding modes. More specifically, let us organise the η i (k) modes and the fermions f x,k i in single vectors as We then consider the correlation functions Γ(k)Γ † (q) and F x,k F † x,k ρ x,k , where the first expectation value is taken with respect to the state in Eq. (8), while the second on the ancillary Hilbert space on the density matrix ρ x,k which is defined in such a way to satisfy (we recall that momenta are positive) We finally impose f x,k i ρ x,k = 0. Given that the density matrix ρ x,k is Gaussian, having established its one-and two-point functions completely fixes the density matrix itself. The matrix C(k) has dimension 4n × 4n and may be written in terms of the correlation matrices C (1) (k) and C (2) (k) in (9) as propagates with its own velocity v k i . At a given time t, the initial set of fermions is divided into two subsets B andB, where in B appear those fermions which are carried in A by the ballistic evolution.
In particular, thanks to the block structure with M and N being 2n × 2n matrices (and M = M † ), it is always possible to define a density matrix ρ x,k such that Eq. (15) holds true. This picture semiclassically describes the initial conditions. Now we consider the time evolution: to each ancillary fermion we associate a velocity through the correspondence (13) Then, at a given time t, we introduce a bipartition of the ancillary Hilbert space B ∪B accordingly to the following rule (see also Fig. 1) We set s A (x, k, t) as the entanglement entropy of such a bipartition, i.e. we construct the reduced density matrix ρ x,k B = TrB[ρ x,k ] and pose Notice that, without explicitly computing ρ x,k B , we can take advantage of the gaussianity of the reduced density matrix and express the Von Neumann entropy in terms of the correlation matrix [77][78][79] . In particular, let C B (k) be the correlation matrix extracted from C(k) in (16) retaining only those degrees of freedom in the B subspace, then it holds The equivalence between Eqs. (20) and (21) is discussed in Appendix A. Notice that the traces in Eqs. (20) and (21) are on very different spaces. In analogy to Eq. (5), we can explicitly perform the integration over x in the case when A is an interval of length . For simplicity, we assume that the quasiparticles are ordered in such a way that v k i > v k j if i < j. This does not imply a loss of generality, since when it is not the case we can always, at fixed momentum, reorder the quasiparticles in such a way this requirement holds true, at the price that the needed reordering is momentum-dependent. Under these assumptions we have where we conventionally set v k 0 = ∞ and v k 2n+1 = −∞. The set of indexes B j,i that must be extracted from the two-point correlation matrix is B j,i = (j, j + 1, ..., i) ∪ (j + n, j + n + 1, ..., i + n) .
The prediction to the entanglement growth provided by our ansatz quantitatively differs from the case where correlations are ignored (i.e., Eq. (5) extended to several species). This is clearly shown in Fig. 2, where we provide a few explicit examples anticipating our analysis of the Ising model of Section III. Notice that although the two curves for correlated and uncorrelated pairs are quantitatively different, the light-cone spreading of entanglement still occurs even in the presence of correlations, as manifested by an initial linear increase followed by saturation to an extensive value in (which must be the same in the two cases). Yet, the growth rate of the entanglement for 2v max t < is different.
A very important physical feature of the new prediction (22) is that it cannot be rewritten only in terms of the mode populations n i (k) and velocities v i (k), but the correlations in the initial state must be taken into account. Thus, contrarily to the standard uncorrelated case 53 , the knowledge of the stationary state is not enough to fix the entire time dependence of the entanglement entropy.

A. Consistency checks
Given that our prediction (21) for the time evolution of the entanglement entropy is after all just a well thought conjecture, its validity must be ultimately tested against ab-initio calculations (either exact or numerical simulations). However, we can perform some non-trivial consistency checks comparing with the already well-established literature. First, we remark that, since the initial state (8) is pure, the density matrix associated with the ancillary Hilbert space ρ x,k corresponds to a pure state as well, i.e. ρ x,k = |Ψ x,k Ψ x,k |, for a certain state |Ψ x,k . This automatically guarantees the following properties.
1. The entanglement entropy S A (t) must be symmetric under exchange A ↔Ā. In Eq. (20) this follows from the facts that exchanging A withĀ is equivalent to B ↔B and that for any set of pairs with weight s A (x, k, t) Eq. (20) is symmetric under this operation.
2. Consider a given set of pairs and their weight s A (x, k, t): if all the particles at time t belong to the same set (either A orĀ), we must have s A (x, k, t) = 0. This is immediately guaranteed by the fact that Tr ρ x,k log ρ x,k = 0, being ρ x,k associated with a pure state.
3. If only one particle species is present, then we must recover the standard quasiparticle prediction 33 . In the single species case, C (1) (k) and C (2) (k) are simple numbers, moreover C (1) (k) is, by definition, the density of excitations C (1) (k) = n(k). Thus the block matrix C(k) (16) specialised to a single species case reads When the quasiparticle of momentum k is in A, while the companion at momentum −k belongs toĀ, the reduced correlation matrix is Thus, from Eq. (21) we get which coincides with the single species weight Eq. (6). Then, in the case where we are interested in a single interval, Eq. (5) is readily recovered from Eq. (22). 4. In the case of several particle species, but uncorrelated (i.e. C (1) (k) and C (2) (k) are diagonal), the generalisation of Eq. (5) to many species is readily obtained, as a straightforward extension of the single-species case.

At infinite time and choosing
A to be a finite interval, we must recover the Von Neumann entropy constructed on the late time steady state, i.e. it must hold true where This property is simply proven: for a very large time, looking at a set of pairs originated in (x, k), at most one quasiparticle belongs to A 80 . The distance between two quasiparticles associated with fermions f x,k i and f x,k j is t|v k i − v k j |, thus if t > /|v k i − v k j | they cannot both belong to A (under the assumption of absence of velocity degeneracies). Therefore, if only the fermion f x,k i belongs to A, constructing the reduced correlation matrix we find exactly Eq. (25) with n(k) → n i (k). Thus, s A (x, k, t) reduces to Eq. (28). Considering the spatial integration in Eq. (12) we simply get a prefactor and Eq. (27) immediately follows.

III. THE INHOMOGENEOUS ISING MODEL
We now discuss the solution of the inhomogeneous Ising model (10) which provides a benchmark for our ansatz. We introduce fermionic degrees of freedom {d j , d † j } = δ j,j through a Jordan Wigner transformation In fermionic variables, the Ising Hamiltonian (10) can be rewritten as Hereafter, we are interested in the thermodynamic limit N → ∞ so that the boundary terms do not play any role and can be discarded. In order to diagonalise the Hamiltonian, it is convenient to move to Fourier space and define where {α k , α † q } = δ(k − q). In Fourier space, the Hamiltonian reads where we exploited the periodicity of h j and defined In the Hamiltonian (32), the periodic field h j couples the modes accordingly to the "roots of unity rule". Now, we introduce several fermionic species splitting the Brillouin zone; in the momentum basis, we define β i (k) fermions as The momentum k of the β i (k) fermions runs on the reduced Brillouin zone k ∈ [0, 2π/n) and they satisfy standard anticommutation rules {β i (k), β † j (q)} = δ i,j δ(k − q). In terms of these fermions, the Hamiltonian becomes We get a more compact notation organising the β i (k) fermions in an unique vector as, and writing the Hamiltonian as The matrix H h (k) can be written as where the matrix elements of the n × n blocks are The desired modes are identified through the diagonalisation of H h (k). Indeed, given the block form of H h (k), it always exists an unitary transformation U h (k) such that Here E h (k) are positive defined diagonal matrices, which are the energies of the modes The modes η i (k) are then identified as the solution of the linear equation where G h (k) is defined as With this definition, we made the Brillouin zone symmetric around zero.
Having diagonalised the Hamiltonian for arbitrary magnetic field, we can now consider a quench changing the magnetic field from {h t<0 j } n j=1 to {h j } n j=1 . Then the pre and post quench modes are connected through a proper Bogoliubov transformation. In particular, from Eq. (42) we have We finally need to show that the initial state (i.e. the ground state {h t<0 j } n j=1 ), when expressed in the post quench modes for a magnetic field {h j } n j=1 , is of the form in Eq. (8). The initial state is the vacuum for the prequench modes, therefore, using Eq. (44), we can readily write the set of equations where for compactness we set It is then straightforward to realise that |0 h t<0 can be written in the form Eq. (8), i.e.
provided the matrix M i,j (k) satisfies the equation We can then conclude that quenches in the inhomogeneous Ising spin chain (10) fall within the framework of our ansatz for the entanglement spreading, which can now be tested. As we saw, finding the energies of the modes E i (k) and the t = 0 correlations of the post quench modes boils down to diagonalising finite-dimensional matrices; even though pushing further the analytical calculations can be cumbersome (especially if several species are involved), this last step can be quickly carried out numerically.
In Fig. 3 we test the ansatz against direct exact numerical calculations in the Ising model for various choices of the pre and post quench magnetic fields. Numerical calculations have been carried out on a lattice of 1200 sites with periodic boundary conditions, by mean of a direct solution of the free fermion model. We consider time t and subsystem sizes such that the finite size of the entire system does not play a role. We focus on a bipartition where A is a finite interval of length . The figure shows that as becomes larger, the numerical results clearly approach the quasiparticles ansatz. In order to provide a stronger evidence for the correctness of our conjecture we provide an extrapolation to → ∞. In this perspective, we assume a regular expansion in powers of −1 where s QP A is the quasiparticle prediction. Performing a fit of our data with the above form at order O( −1 ), we obtain the extrapolation that are represented as crosses in the figure. It is evident that these extrapolations perfectly match the quasiparticle ansatz for all the considered quenches.

IV. TIME EVOLUTION OF THE ORDER PARAMETER CORRELATIONS
Although the main focus of this work is the entanglement entropy, the quasiparticle picture may provide useful information even for other quantities, primarily the order parameter σ x j and its two-point correlation 56 . For the quantum quench in the homogeneous Ising model, it has been shown that the quasiparticle prediction is qualitatively and quantitatively correct only for quenches within the ordered phase 26,38 . For quenches from the ordered phase to the paramagnetic one, the quasiparticle picture can be heuristically adapted to provide correct results 26 . Instead in the case of quenches starting from the paramagnetic phase, it is still not known whether it is possible to use these ideas to have an exact ansatz; anyhow a discussion of this issue is beyond the scope of this paper, see Ref. 26. Consequently, in this section we limit ourselves to extend the quasiparticle picture for the order parameter correlations to initial states with correlated quasiparticles for quenches within the ferromagnetic phase and to test the prediction in the inhomogeneous Ising chain (10).
Correlation functions of the order parameter are non local objects when expressed in the fermionic basis Although the computation of the correlator (and its time evolution) ultimately boils down to an extensive use of the Wick theorem and evaluating determinants, the large number of the involved degrees of freedom makes the calculation very complicated. In Ref. 26 a first-principle calculation was carried out in the homogeneous Ising model, resulting in a scaling behaviour that can be interpreted a posteriori in terms of quasiparticles. While in principle possible, we do not try to generalise the complicated methods of Ref. 26, but, inspired by the resulting expression, we directly attempt a quasiparticle ansatz which is then numerically verified. A posteriori, we will see that our ansatz fails to describe those situations where σ x j = 0 on the initial state, as expected on the basis of the results for the homogeneous case 26 . Inspired by Ref. 26, we conjecture the following ansatz for the logarithm of the two-point correlator where A is the interval of extrema j and j , the quantity p A (x, k, t) is discussed hereafter. An explicit space integration may eventually lead to a formula similar to the entropy one in Eq. (22). Our ansatz, relays on the following assumptions, similar to those used in the entanglement entropy case: i) quasiparticles generated at different spatial points or belonging to pairs of different momenta are uncorrelated; ii) the large distance behaviour of the correlator Eq. (50) is ultimately determined by the string (In passing: most likely ii) is the hypothesis that is failing for quenches from the disordered phase.) Consider then a semiclassical computation of e −iπ j l=j+1 d † l d l : since pairs originated at different position or having different momentum are uncorrelated, the expectation value should factorise in the contribution of each set of pairs. Equivalently, the logarithm must be additive, justifying the form of Eq. (51). Now, in order to find the proper ansatz for p A (x, k, t), we recognise that the string (52) simply counts the parity of the number of fermions within the interval (j, j ). Therefore, using the same notation of Section II for the auxiliary fermions f x,k j , we take as an ansatz where in the product only those fermions which semiclassically belong to the interval x + v k i t ∈ A must be considered. The expectation value is taken on the auxiliary Hilbert space as in Section II.
The time evolution of the order parameter itself may be accessed from Eq. (51) using the cluster decomposition principle, obtaining Under the further assumption that σ x j is translational invariant in the scaling limit (as it is the case here), we therefore obtain a quasiparticle prediction for the one point function. This observation, besides providing an additional result, also helps a posteriori to understand the regime of applicability of the quasiparticle ansatz. Indeed, from Eq. (53) |p A (x, k, t)| ≤ 1 (it is a product of terms that are all smaller than 1). Therefore, the r.h.s. of Eq. (51) is surely negative (or at most zero), implying an exponentially decaying | σ x j |, which cannot be correct if the order parameter is zero in the initial state.
We close this section mentioning that by the repeated use of the Wick Theorem, p A (x, k, t) can be efficiently formulated in terms of a determinant of a correlation function of the fermions. Assume that the fermions f x,k ij for a set of indexes i j , j ∈ {1, ..., n } are those belonging to A. Then, we consider a 2n × 2n antisymmetric matrix A in the inhomogeneous Ising model. We test the time evolution of the two-point correlator of the order parameter against the quasiparticle ansatz for increasing separation . It is evident that increasing the numerical data quickly approach the quasiparticle ansatz. The initial exponential decay is due to the evolution of the one-point function of the order parameter since, on that time scale, we have σ x 1 σ x σ x 1 σ x . At late times, saturation to the steady state value is observed.
defined as In terms of the matrix A, it holds (see Appendix B) We tested our ansatz against exact numerical calculations for quenches in the inhomogeneous Ising chain. We found that for all quenches within the ferromagnetic phase, our quasiparticle prediction perfectly reproduces the numerical data in the space-time scaling limit t, → ∞ with t/ fixed. Two representative cases are shown in Fig. 4. It is evident that by increasing the numerical data approach the prediction, although finite size correction are clearly visible at small , but these are much smaller than those for the entanglement entropy. It is clear from the result that even for the two-point function of the order parameter, the "light-cone" spreading of correlations persists in the presence of correlated quasi-particles. We also checked that for other quenches (i.e. from and to the paramagnetic phase) the conjecture does not work, as expected. Finally, to be exhaustive, we remind the reader that the quasiparticle prediction is correct for quenches to the critical point from the ordered phase, but not for quenches originating from the critical point 26 .

V. CONCLUSIONS
In this manuscript we generalised the semiclassical quasiparticle picture for the spreading of entanglement and correlations to global quantum quenches with multiple species of quasiparticles that show momentum-pair correlations in the initial state. The main new physical result (compared to the standard uncorrelated case) is that the information encoded in the mode populations n i (k) of the single species and their velocities v i (k) are not enough to determine the time evolution of the entanglement entropy and correlations. We show that among the systems displaying this phenomenology for the quasiparticles, a remarkable example is the Ising chain with a periodically-modulated transverse field, which is easily mappable to a free fermionic theory. We then use this model to test our predictions for entanglement and correlations against exact numerical calculations, finding perfect agreement in the scaling regime.
A simple and straightforward generalisation of our results concerns the case when the initial state is inhomogeneous on a large scale so to be describable by the generalised hydrodynamics approach [81][82][83] . For example, we have in mind the joining of two different thermal states or groundstates of different Hamiltonians producing correlated quasiparticles. In this case, it is very simple to merge the results of the present paper with those of Ref. 45. The spatial variation of the entanglement entropy may be properly captured by a term that locally is given by Eq. (21).
Another inhomogeneous setup in which several correlated pairs can be produced is that of moving defects 84 , where an external localised perturbation is dragged at constant velocity in an otherwise homogeneous free lattice system. The moving impurity can be regarded as a source of quasiparticles and the propagation of entanglement could fall within our frame, again supplemented with the generalised hydrodynamics [81][82][83] .
A simple generalisation of our results is the time evolution of Rényi entanglement entropies that for free models just requires a minor variation of the form of the kernel (21), see for example the discussion in 85.
Finally, a more difficult open problem concerns the spreading of entanglement in interacting integrable models. For these models there are almost always multiple species of quasiparticles (they are bound states of the lightest species), but known integrable initial states do not have correlations between them and the general evolution of the entanglement entropy has been understood in 53-55 (for the Rényi entropies see 85-88 while for inhomogeneous systems see 89). This lack of correlations in solvable initial states has been sometimes related to the integrability of the quench problem itself 90,91 . Yet, it would be very interesting to understand whether, at least in in some models, it is possible to have correlated quasiparticles which would require the generalisation of our approach to interacting integrable systems.
The gaussianity of the ensemble and the diagonal correlator, necessarily implies that the density matrix in the new basis is written as Hence the entanglement entropy easily follows To recover Eq. (A2) it is enough to notice that the spectrum of C consists in . Therefore we can equivalently write Eq. (58) can be proven, e.g., following Ref. 26. Hereafter, we consider the general problem of determining where f i are fermionic operators {f i , f † j } = δ i,j and the state is Gaussian in these fields. We first rewrite p as where we naturally introduced the operators b i as It holds true {b i , b j } = 0 for i = j. Let us introduce the following antisymmetric matrix A i,j Therefore, as noticed in Ref. 26, by mean of a simple application of the Wick Theorem it can be shown where Pf(A) is the Pfaffian of the matrix A. Given that |Pf(A)| = det(A), this concludes the proof. Indeed, the definition of the matrix A in Eqs. (55,56,57) of Section IV is nothing else that the analogue of the matrix A.