Prethermalization from a low-density Holstein-Primakoff expansion

We consider the non-equilibrium dynamics arising after a quench of the transverse magnetic field of a quantum Ising chain, together with the sudden switch-on of a long-range interaction term. The dynamics after the quantum quench is mapped onto a fully-connected model of hard-core bosons, after a suitable combination of a Holstein-Primakoff transformation and of a low-density expansion in the quasi-particles injected by the quench. This mapping holds for a broad class of initial states and for quenches which do not cross the critical point of the transverse field Ising model. We then study the algebraic relaxation in time of a number of observables towards a metastable, pre-thermal state, which becomes the asymptotic steady state in the thermodynamic limit.


I. INTRODUCTION
The constant progress in manipulating cold atomic gases has provided insight into the non-equilibrium dynamics of isolated, interacting quantum many-body systems. The experimental observation of the absence of relaxation in an (almost) integrable one-dimensional Bose gas 1 and the appearance of an intermediate, metastable regime in the dynamics of a non-integrable system on time-scales much shorter than those required for its equilibration 2-7 call for a better understanding of the mechanisms underlying quantum relaxation.
The inherent unitarity of the evolution makes the emergence of relaxation, thermalization and, more generally, decoherence a subtle issue, since they cannot really occur in isolated systems 8 : in fact, a unitary dynamics cannot turn an initial pure state into a thermal distribution, which is a proper statistical mixture. Consequently, one may look for signs of thermalization in a suitable set of local observables 9 , whose expectations, in the long-time limit, are equivalently captured by the thermal distribution -its temperature being determined by the average energy of the initial state [10][11][12] .
In order to probe the relaxation of quantum manybody systems, the protocol known as quantum quench has provided a convenient conceptual framework, widely employed in recent investigations 13,14 . It consists in suddenly changing a parameter of the many-body Hamiltonian (magnetic field, interaction strength, or coupling among spins, etc) on time scales short enough to leave the system initially frozen in the ground state of the prequench Hamiltonian. This state will then evolve according to the post-quench Hamiltonian: since in general it does no longer constitute an energy eigenstate, in the case of a global quench it will have a finite energy density above its ground state, and hence a finite density of quasi-particles, whenever they are defined.
A richer scenario is expected to emerge when the quench involves two parameters with competing effects, i.e., one preserving and one spoiling the integrable nature of the quantum system. For instance, if a weak integrability-breaking interaction is switched on, while simultaneously quenching a control parameter which in the absence of interactions would preserve the integrability of the model, the intermediate time evolution is expected to be still approximately determined by the integrable (non-interacting) part of the Hamiltonian 44 . As a consequence, a two-stage relaxation emerges: First, the system approaches the GGE corresponding to the integrable part of its Hamiltonian, such that this GGE approximately describes a metastable state at intermediate times after the quench. Afterwards, on time scales which depend on the strength of the interaction, the effect of quasi-particle scattering becomes relevant and drives the system towards an eventual, bona fide thermalization. This mechanism is usually referred to as prethermalization (or as prerelaxation when, instead, the perturbation breaks non-abelian integrability into integrability 45 ), and, although this notion was firstly employed in the context of high-energy quantum field theories 46 , during the last few years it has been extended to the domain of condensed matter physics. Signatures of prethermalization have been reported in long-range interacting quantum simulators 47 , in systems of spinless fermions with tunable interactions [48][49][50] , in two-dimensional spinless Fermi gases 51 , in metastable superfluids 52 , in the three-dimensional isotropic Heisenberg model 53 , for interaction quenches in the fermionic Hubbard model [54][55][56][57] , during aging in interacting φ 4 -theories (either isolated or in contact with a bath) [58][59][60][61] , in noisy Ising models 62,63 , and in interacting Luttinger liquids [64][65][66] .
In this work, we investigate in more detail the prethermalization dynamics of the long-range model we originally introduced in Ref. 67. Starting from the quench of the quantum Ising chain, we introduce a long-range spin-spin interaction which breaks many, but not all, of the original conservation laws, as we detail in Sec. II. We first show in Sec. III that an exact mapping exists to a model of hard-core bosons on a fully-connected lattice. As long as the quasi-particle density generated by the quantum quench remains sufficiently small, one can think of the hard-core constraint as being substantially ineffective, and thus treat the bosons as ordinary ones. Formally, this corresponds to a lowest-order truncation of the Holstein-Primakoff transformation 68 . This approximation -which holds for small quenches up to very large times (cfr with Fig. 4), renders the theory effectively noninteracting and allows us to map the non-equilibrium dynamics of the original model onto the relaxation of an integrable one. We should notice that a similar low-density expansion in the quasi-particles produced by the quench has already been employed for studying the relaxation dynamics of isolated interacting quantum systems in different contexts, and including both integrable 11,12 and non-integrable systems 23,69 .
We then proceed to solve numerically the dynamics of the approximately-equivalent bosonic model, highlighting the presence of plateaus in the relaxation of some physically relevant observables, which are typically approached algebraically in time; the main numerical results are reported in Sec. IV. In Sec. V (which extends the results of 67 ) we discuss in more detail the range of applicability of our approach and offer some numerical evidence that it provides reasonably accurate results even for quenches at the critical point or across it, in spite of the fact that they generate highly-populated modes and hence violate the low-density condition.
Contrary to several previous studies on prethermalization plateaux 48,49,51,[55][56][57]65 , our analysis does not rely on a straightforward perturbative expansion in the parameter controlling the interaction among quasi-particles of the pre-quench integrable model and offers a new perspective into the study of non-equilibrium dynamics of interacting systems, which is in general a formidable analytic challenge.
The model presented here is exactly solvable in the thermodynamic limit, since the long-range nature of the interaction suppresses fluctuations and makes mean-field theory exact (see Section II). This fact has been employed in Ref. 70 in order to solve, inter alia, the time evolution of this model through a mean-field mapping.
Despite the diversity of approaches developed in order to tackle the quench dynamics of this system, we think it is worth reporting here in full detail our approach, since it suggests that in certain cases pre-thermalization in interacting systems can be understood as the dynamics of an approximate integrable model, which emerges at intermediate time scales and which is not perturbatively connected to the pre-quench model.

II. THE MODEL AND THE QUENCH PROTOCOL
The interacting spin model under study, with Hamiltonian H = H 0 + V , is a quantum Ising chain of N sites with periodic boundary conditions, where is the standard short-range Ising model, with J the nearest-neighbor coupling and Jg the transverse magnetic field, while σ x,y,z i denote the standard spin-1 2 operators acting on site i, which can be represented as Pauli matrices (see Eq. (A2)). This H 0 is perturbed by a longrange interaction of the form where M z ≡ i σ z i is the global transverse magnetization and M z its time average calculated according to Note that, throughout the paper, we set = 1 and measure energies and inverse times in units of J, thereby fixing J = 1.
In order to drive this system out of equilibrium, we consider here a composite quench protocol: We "prepare" the system in the ground state |ψ 0 (g 0 ) of the pre-quench Hamiltonian H 0 (g 0 ). At t = 0, we vary the magnetic field (g 0 → g) on a time scale so short that the state is effectively left unaffected, simultaneously switching on the interaction term V . In other words, Accordingly, the dynamics for t ≥ 0 takes place under the action of the post-quench Hamiltonian, i.e., The Hamiltonian H 0 (g) is known to be integrable [71][72][73][74] and can be mapped via a sequence of Jordan-Wigner 75 , Fourier and Bogoliubov transformation to a free model of fermionic quasi-particles γ k , γ † k , k being the quasimomentum or reciprocal wave vector (for a brief summary of this approach we refer to Appendix A). The integrability of H 0 is however spoiled in the presence of the many-body interaction V which generate scattering among the quasi-particles and allows a redistribution of momentum and energy among the different mode populationsn k = γ † k γ k . Note that switching on the interaction V at time t = 0 without the quench g 0 → g, would be sufficient to make the ground state |ψ 0 (g 0 ) of the pre-quench Hamiltonian a non-equilibrium one for H; however, the quench g 0 → g in the magnetic field has the effect of providing a nonvanishing quasi-particle density from the very beginning of the evolution; in particular, introducing the shorthand |0 = |ψ 0 (g 0 ) , one obtains where θ k (g) and θ k (g 0 ) are the Bogoliubov angles corresponding to H 0 (g) and H 0 (g 0 ), respectively, and can be determined from Eq. (A9), whilen k =n k (g) is the aforementioned mode population of the post-quench Ising Hamiltonian H 0 (g). Hereafter, all quantities are calculated in terms of the post-quench fermionic basis, unless otherwise stated (e.g., θ k is shorthand for θ k (g)). Some examples of post-quench populationsn k (t → 0 + ) are displayed in Fig. 1 for various values of the initial and final transverse fields g 0 and g. Before proceeding, we comment on the subtraction of M z in Eq. (2), which is the long-time average value of the global transverse magnetization in a quenched noninteracting Ising model (λ = 0). The reason for this choice (further discussed in Appendix B) is grounded in the non-equilibrium dynamics of the connected cor- of the global transverse magnetization M z for which (see also Refs. 76 and 77) This fact indicates the presence of information which is never really lost, as measurements of this observable separated by an arbitrary time τ are still correlated. The "cluster property" can be restored also for this non-local observable with the subtraction of M z , which practically corresponds to subtract ab initio from the interaction term V a term which is extensive in the system size (∼ O(N )) and quadratic in the Ising integrals of motion n k .

A. Remaining conserved quantities
The integrability of the Ising model is typically expressed in terms of the conservation of the quasi-particle densitiesn k . This breaks down with the introduction of the interaction V . However, one can still identify an extensive number of conserved quantities. Their presence is more easily highlighted in the fermionic representation of the model (see Appendix A), in which H 0 reads where are Nambu spinors constructed with the fermionic annihilation (creation) operators γ k (γ † k ), while is the dispersion relation of the quasi-particles. The sum runs over the discrete values k = 2π(n−1/2)/N with n = 1 . . . N/2, while σ x,y,z denote the usual Pauli matrices reported for convenience in Eq. (A2). In the same basis, the interaction term V is in the form In the notation of Eq. (11) the long-range nature of V becomes apparent, as the terms ψ † k ψ k ψ † q ψ q clearly connect every possible pair of momenta k, q. The total Hamiltonian H is thus which proves that, for any finite N , most populations are not conserved. The only exceptions are the two extremal cases k = 0 and k = π, which however are present only in the thermodynamic limit N → ∞ and can therefore be safely disregarded. We notice now that both the Bogoliubov angle θ k and the two-particle operators γ † k γ † −k and γ −k γ k appearing on the r.h.s. of Eq. (13) are odd under the inversion of the momentum k → −k, i.e., γ k γ −k = −γ −k γ k . This implies [H,n k ] = [H,n −k ], and therefore commutes with the total post-quench Hamiltonian H (similar operators depending on fermionic pairs at op- Post-quench population of the Ising quasi-particle modesn k (cfr. Eq. (A14)) at time t = 0 + as a function of the final transverse field g and the momentum k for four different choices of the pre-quench field, from left to right g0 = 8, 2, 1 and 0.5. The dashed lines correspond to g = g0, i.e., the no-quench case, in which the initial fermionic populations vanish. Ising modes become significantly populated only for quenches which cross from the ferromagnetic (g < 1) to the paramagnetic (g > 1) phase or viceversa, or start or end at the critical point g = 1. Note that these plots are invariant under the exchange g0 ↔ g, sincen k has the same property (see Eq. (6)).
posite momenta have been discussed also in Ref. 78).

III. MAPPING TO HARD-CORE BOSONS
Thanks to the set of N/2 mutually commuting constants of motion I k (with I −k = −I k ), the spin chain described by H can be exactly mapped onto a quadratic (yet non-diagonal) Hamiltonian of hard-core bosons, as we now proceed to show, completing the brief discussion of Ref. 67.
We start by analysing the structure of the Hilbert space in the presence of these constraints: first of all, each I k involves a pair of modes with opposite momenta k and −k. Consequently, it acts non-trivially only on the fourdimensional subspace spanned by the vectors |∅ k , |k = γ † k |∅ k , |−k = γ † −k |∅ k and |k, −k = γ † k γ † −k |∅ k , where the vacuum |∅ k is annihilated by both γ k and γ −k . From the definition (14) it is not too difficult to see that which identifies the possible eigenvalues {−1, 0, 1} of each I k . Note that 0 is doubly degenerate, while ±1 have no degeneracy. Since the I k 's commute, one can independently fix their eigenvalues for each k > 0, e.g.
indicates a list of possible eigenvalues of I k , ordered with increasing allowed values of k = (2n + 1)π/N on the lattice, starting from n = 0. Each choice uniquely identifies a subspace (sector) of the Hilbert space which is orthogonal to the others and is not dynamically connected to them (i.e., the Hamiltonian H acquires a blockdiagonal structure). For instance, take a chain of N = 4 spins; this has four quasi-momenta ±π/4 and ±3π/4 and, consequently, two conserved quantities I π/4 and I 3π/4 . The Hilbert space fragments into 9 dynamicallydisconnected sectors such as {1, 1} -corresponding to the sole vector |π/4 ⊗ |3π/4 -or {0, −1} which is spanned, instead, by the two vectors ∅ π/4 ⊗ |−3π/4 and |π/4, −π/4 ⊗ |−3π/4 . Labels such as (16) contain N/2 distinct eigenvalues, hence the total number of sectors is 3 N/2 . As mentioned above, the "0"s are doubly-degenerate, implying that the actual dimension of a sector is 2 N0 , with N 0 the total number of "0"s appearing in the corresponding label. Considering still the N = 4 example reported above, the sector {1, 1} has no zeroes, so N 0 = 0 and is one-dimensional, while {0, −1} has instead N 0 = 1 and the sector is two-dimensional. Recalling that the Ising Z 2 symmetry σ x → − σ x in the spin formulation corresponds to the preservation of the parity (−1) kn k of the fermion number, we also see that the number N/2 − N 0 of ±1s in a label determines the Z 2 -parity of the sector, as it counts the number of unpaired quasi-particles present. For example, for N = 8, the sector {0, 0, 1, −1} has N 0 = 2 and is thus four-dimensional, while N/2 − N 0 = 2 is even, implying that all its vectors transform trivially under the . By construction, each sector also carries definite momenta kI k corresponding to each (k, −k) subspace.
We shall now look for operators which leave all sectors invariant; by the properties stated above, those operators must preserve both parity (they can only change the number of fermions by an even amount) and momentum components (they must carry zero net momentum in each (k, −k) subspace). For any fixed choice of k, it is sufficient to consider only combinations of the basic creation and annihilation operators γ ±k , as all the others do not act on the corresponding subspace. The only (non-trivial) choices satisfying all constraints are the quadratic operatorŝ which represent the populations and the creation and annihilation of pairs with zero net momentum, respectively, and the quartic onen All other possible operators can be re-expressed in terms of these ones by making use of the canonical anticommutation relations of the fermionic operators, see (A10). Products of the operators above at equal or different momenta (e.g.,n k b † q b p ) will leave all sectors invariant. Every operator which commutes with all I k 's can be expressed in terms of (products of) the ones in Eqs. (17) and (18) and the Hamiltonian H constitutes no exception. As a matter of fact, it can be written in terms of the sole "pair" operators as with This implies that the relevant dynamics is described by the interaction of pairs of quasi-particles with zero net momentum, rather than of single fermionic modes, and that we can therefore reformulate the problem in terms of these new fundamental modes. In order to do so, we shall first investigate their nature. Since they obey they behave almost, but not exactly, as hard-core bosons. In fact, this would require the last anticommutator above to be 1. On the other hand, by noticing that both b k and b † k act as the null operator in a sector with , we can effectively expunge them from H . This operation leaves behind only those pair operators corresponding to momenta q for which I q = 0; within the corresponding eigensector they then satisfy the hard-core constraint. Thereby, in a sector characterized by having N/2 − N 0 unpaired quasiparticles, the projected Hamiltonian effectively describes a fully-connected model of hard-core bosons on a lattice with N 0 sites. The corresponding basis can be obtained by setting, for every involved k, the correspondence |∅ k → |0 k , |k, −k → |1 k , where 0 and 1 stand for the boson being absent or present, respectively. This reinterpretation clarifies the effect of introducing V in the Hamiltonian of the Ising model: with N/2 conserved quantities still remaining, in fact, we cannot expect that the resulting model is completely nonintegrable and, indeed, we identify sectors in which it is trivially solvable, which are the ones almost completely lacking pairs (i.e., those whose labels display just a few 0s). For example, the 2 N/2 totally-unpaired sectors collectively represent the zero-energy eigenspace of the Hamiltonian H and coincide with the corresponding one of H 0 ; furthermore, each of the (N/2) × 2 N/2−1 sectors having a single pair is two-dimensional and the corresponding reduced Hamiltonian is already cast in the diagonal form in the basis {|0 k , |1 k } introduced above. This is due to the presence of an additional symmetry, namely the conservation of the parity which further splits each sector in two halves of equal dimension. From a physical point of view, this is associated with the fact that H in Eq. (19) either does not affect their total number, k b † k b k , or it simultaneously creates or destroys two pairs. Although the structure of the states space becomes progressively more complicated as N 0 grows, it is clear that the model cannot really display non-integrable features as long as the dynamics remains confined in the small-N 0 sectors.
The situation is reversed for N 0 ≈ N/2 1; although the corresponding eigensectors are exponentially smaller than the global Hilbert space (whose dimension is 2 N ), their dimensions are still exponentially large in the number of sites, as expected for a truly many-body problem. Note that, in spite of the fact that the Hamiltonian (19) is quadratic in the pair operators, it does not define a free theory, due to the hard-core nature of the bosons; indeed, trying to diagonalise it by applying a generic Bogoliubov rotation intertwines the mixed commutation/anticommutation relations in Eq. (21), implying the impossibility of interpreting b k and b † k as annihilation/creation operators of any kind of particles. This relates to the fact that hardcore bosons are intrinsically interacting particles, as they can be thought as ordinary bosons subject to infinite "onsite" interparticle repulsion.
Within our setting, the dynamics always starts from the totally-paired sector N 0 = N/2, independently of the values of the quench parameters g 0 , g and λ. The prequench state |0 , in fact, can be represented as 72,76 where t k = tan (θ k (g) − θ k (g 0 )) and the effect of the operator on the r.h.s. is to generate pairs on the vacuum of the post-quench Hamiltonian. Thereby, this whole class of initial states constitutes a suitable choice for highlighting the effects of V on the dynamics, as it never involves the highly-unpaired sectors.
A. Low-density approximation For N 1, the interacting problem is hard to solve; in this case, instead of employing the usual perturbative expansion in the interaction strength, one can adopt a different approach, which makes use of the quadratic structure of the Hamiltonian. The point is that the hard-core constraint is expected to become effective only when the occupation of a given mode approaches 1; as long as the quasi-particle densities remain much smaller than that, they behave approximately as standard bosons. This condition must be first satisfied by the initial state. By inspecting Fig. 1 one sees that some modes get occupied to a relevant degree only in the case of quenches that connect different phases (g > 1 and g 0 < 1 or vice-versa) or that start (g 0 ≈ 1) or end (g ≈ 1) close to the critical point. Accordingly, we expect this approximation to be accurate for quenches within the same phase.
From a formal point of view, the pair operators b k and b † k introduced in Eq. (17) can be expressed in terms of standard bosonic ones a k and a † k by means of a Holstein-Primakoff transformation 68 for low densities, one can then think of expanding the square roots as power series of their arguments, with where a † k a k is the bosonic number operator. Hence, as long as the average and the fluctuations of a † k a k remain small, one can conveniently truncate the expansions above. Further details are provided in App. C. What makes this approach particularly convenient is that, by expanding at the lowest order b † k ≈ a † k and b k ≈ a k , one obtains in each sector a quadratic Hamiltonian which can now be diagonalised by a Bogoliubov rotation. Denoting by K S the set of paired momenta present in a given sector S (see also note 79 ), the expression of the reduced Hamiltonian acting on S is Note that the hard-core constraint b 2 k = b †2 k = 0 is reflected in the vanishing diagonal part of the matrix α kq (see Eq. (20a)), which prevents terms such as a 2 k from appearing in H S .
As we show in the following, for most choices of the initial state -as long as g is not too close to the critical point g c = 1 -the approximation above successfully describes the early stages of the dynamics. At longer times, however, the quadratic approximation, encoded in the Hamiltonian (28), is spoiled because higher-order terms in the expansion of the Holstein-Primakoff representation (27) become important, introducing interactions among the bosons.
IV. PRE-THERMALIZED REGIME With the notation of the previous Section, the ground state |ψ 0 (g 0 ) of H 0 (g 0 ) lies in the totally-paired sector of H, corresponding to the string of eigenvalues {000 . . . 0}. Thereby, the dynamics of any (combination) of the invariant quantities listed in Eqs. (17) and (18) can be entirely determined within this subspace. Of course, working with a free bosonic Hamiltonian such as Eq. (28) comes with the substantial advantage of reducing the many-body problem to a one-body one; in other words, it is not necessary to diagonalise the full operator on the whole eigensector, but it is sufficient to determine the one-particle spectrum by applying an appropriate Bogoliubov transformation which casts H in the diagonal form where {E q } q are the energies of the bosonic modes (obtained through exact numerical diagonalization) and C is an incosenquential constant. This problem amounts to the diagonalization of a N × N matrix, and is thus of polynomial complexity in N . As a consequence, one can perform a numerical analysis up to rather large system sizes. Details about the diagonalization procedure and the numerical computation of the relevant observables are provided in Appendix D. By truncating the power series (27) derived from the Holstein-Primakoff representation to its leading terms, we are actually neglecting any form of interaction, as we end up with an Hamiltonian that is non-interacting in each sector. On the other hand, taking into account next-to-leading orders clearly produces the appearance of higher-than-quadratic terms in Eq. (28), which therefore give rise to quasi-particle interactions. In light of this fact, we must conclude that there can be no sign of thermalization as long as the low-density approximation captures the physics of the system. In other words, if the system eventually leaves the prethermal regime highlighted here, it must do so on time scales longer than the regime of validity of the truncation which, as shown below, may encompass quite long time intervals. Prethermal features may only be observed within this regime, as it represents a stage at which interactions among quasiparticles have not yet become dominant, leaving the system observables to relax on the quasi-stationary values associated to pre-thermalization. The typical prethermal behavior and time scales are illustrated in Fig. 2, where representative plots are reported of the temporal evolution of the central mode population n π/2 for various values of the system size N . Quite clearly, after a short initial transient (of duration t ∼ 10, independently of N ), well-defined plateau values are attained in all cases, which define an intermediate, metastable stage of the evolution. As mentioned above, however, the dynamics of observables such as the occupations n k can be decomposed in a finite collection of modes oscillating with some specific frequencies; accordingly, at finite size, the destructive interference which gives rise to the aforementioned plateaus cannot last indefinitely and, indeed, oscillations reappear at a recurrence time t R N/2. The pre-thermal behavior highlighted above for single populations is also reflected in the evolution of extensive observables such as the total number N p of quasiparticles and the Ising energy which represent two examples of a wider class of quantities which can be generically expressed as corresponding to C k = 1 and C k = k , respectively. The evolution of N p and E I has been investigated in Ref. 67, highlighting the presence of similar prethermal plateaus. The oscillations preceding these plateaus decay algebraically as t −α , with α 3. The fact that the same exponent α appears in both quantities is likely to be related to the fact that they belong to the same class defined by Eq. (33). Accordingly, we can reasonably expect this exponent to characterize the whole set of observables O(t) (see Eq. (33)), possibly apart from very specific choices of the coefficients C k . In Fig. 3 we display N p and E I for various values of N and λ, showing that the typical amplitude of the oscillations around the corresponding quasi-stationary values (E I,qs and N k,qs ) scales as λN as long as λ is sufficiently small (i.e., λ 0.5) and N sufficiently large (N 40). Note, however, that the recurrence time t R explicitly depends on N , thus curves with different N can be collapsed one onto the other, as in Fig. 3, only until the first recurrence appears. On the other hand, this proves that the exponent α of the algebraic decay depends on neither N nor λ.

V. RANGE OF VALIDITY OF THE TRUNCATION
We address here the range of validity of the truncated Holstein-Primakoff expansion and its dependence on the quench parameters. First, we assess the accuracy of the approximation on the basis of the comparison between the bosonic formalism introduced above (see Eq. (29)) and the full diagonalization of the fermionic Hamiltonian (19). Although the latter is still restricted to the totally-paired sector, it represents an exponentiallycomplex problem, which we could diagonalize only up to N = 20 sites. Still, this allows us to verify that there is a wide range of values of the quench parameters for which the approximation and the original model are quantitatively almost indistinguishable. This is demonstrated, e.g., by Fig. 4, where we report the evolution of the Ising energies E  includes time scales much longer than those typically required for pre-thermal features to emerge.
As expected, as the post-quench field g approaches the critical value g c = 1, the approximation becomes less accurate. This agrees with the intuitive picture that the closer g is to g c , the larger the initial populations n k (t = 0) are. As a matter of fact, we have numerically verified that for g = 1.01 and all the other parameters fixed as above, the agreement between the two curves in Fig. 4 remains within 2% up to t 10 3 , which implies that, as long as the interaction λ is small, the low-density approximation features a considerably wide range of applicability. In order to probe the effects of the interaction, we have performed an analogous comparison upon varying λ. We verified that the relative discrepancy between the kinetic energy computed in the bosonic and fermionic represenatations, increases upon increasing the interaction λ. Still, for intra-phase quenches this relative discrepancy remains small also for very long times: for instance, for g 0 = 8, g = 3.5, and λ = 0.9 we have verified that | | keeps within 0.2% up to the time scales t ∼ 10 4 probed in Fig. 4: A numerical analysis carried out for smaller system sizes also suggests that the accuracy improves upon increasing N . Quenches within the ferromagnetic phase also display low discrepancies for very long times, whereas for quenches across the phases these become important on much shorter time-scales (e.g., for g 0 = 3.5, g = 0.8 and λ = 0.1 one finds that | | is larger than 1% already for t 10 3 ).
Due to the interacting nature of the model in Eq. (12), it is not possible to employ the exact diagonalization for a direct comparison much beyond the system size considered in Fig. 4. However, an indirect measure of the accuracy of the approximation can also be extracted from the   (12), whereas red ones refer to the bosonic approximation, Eq. (30). This shows that at very long times (see the time scales on the horizontal axis) the latter starts to deviate from the actual evolution, signalling a breakdown of the low-density approximation, upon which it is based. However, this breakdown occurs at increasingly longer times as the post-quench value of g departs from the critical point in the paramagnetic phase, as exemplified by the curves in the inset, which display the same two quantities but for g = 3.5, which are practically indistinguishable. effective integrable bosonic description, Eq. (30). As we prove in Appendix C, as long as the bosonic occupations N k ≡ a † k a k satisfy, in the course of the quench dynamics, the inequality (where we averaged over the pre-quench state) the effects of the terms neglected in the truncation can be expected to be small. When this condition is not satisfied, the approximate dynamics departs from the real one, making the approximation no longer accurate. In Fig. 5(a) we show the evolution of both N k (left scale, black lines) and N 2 k (right scale, red lines) for a chain of N = 50 spins up to time t = 50, corresponding to approximately twice the recurrence time t R ≈ 25 (the curves in Fig. 5 are displayed for various values of wave-vector k from top to bottom). This shows that the two quantities are still closely related within the prethermal regime and that the free-boson approximation is still able to capture the relaxation reasonably well, even for highly-populated modes. Figure 5(b) shows the corresponding relative error  for three different values of the momentum k. Interestingly enough, lower-populated modes display smaller errors compared to higher-populated ones. The worst case among the ones displayed in panel (b) corresponds to momentum k = π/50 within this time window, and it anyhow shows a discrepancy smaller than 1%. This constitutes thereby an upper bound on the error committed by employing this approximation in these parameter ranges. As expected, the maximal error becomes much smaller for intra-phase quenches (within the same time range as the one in Fig. (5(b)), it remains < 10 −6 with g = 3.5 and N = 30, data not shown) and increases when crossing the critical point (oscillating around 1% for g = 0.5 and N = 30). Furthermore, the relative er-ror increases for quenches across the two phases and it displays a weaker increase upon increasing N .

VI. CONCLUSIONS
In this work we reported on a viable approach to study the dynamics of a quenched interacting quantum spin chain with long-range interactions, which originates from the perturbation of a one-dimensional quantum Ising model.
We employed a composite quench of the Ising transverse field and of the interaction coupling, in order to initiate the non-equilibrium dynamics. While the quench of the transverse field does not affect the conserved quantities of the Ising chain (the occupation number of Bogoliubov quasiparticles), the interaction introduces scattering among them, reducing the integrals of motion of the model to a smaller subset.
For the specific quench protocol and initial states considered here, the study of the non-equilibrium dynamics can be based on a lowest-order Holstein-Primakoff expansion in the density of the quasi-particles injected after the quench (see Eq. (27)). Specifically, for quenches that do not cross the critical point, the density of quasiparticles produced is small enough to allow a mapping of the model and of its intermediate-time dynamics into the effective evolution of an integrable model of fullyconnected bosons (see Eq. (28)). This allows us to extract the algebraic relaxation of a number of observables towards a state which constitutes the asymptotic steady state of the system in the thermodynamics limit: for N → ∞ the model becomes exactly solvable, since the long-range nature of the interaction suppresses fluctuations and makes the mean-field theory exact (in this respect see Ref. 70).
On the other hand, for large but finite system sizes, the steady-state approached by observables is expected to be destabilized in the long-time limit, leading eventually to a thermal steady state, as a result of the energy redistribution due to the scattering among quasi-particles. In Ref. 67 this aspect has been explored via a kinetic equation for the time averages of occupation numbers of the Ising quasi-particles. In addition to that, we think that inspecting, in the future, the impact of higher order terms in the Holstein-Primakoff expansion could provide insight into the relaxation dynamics of this model; they could be included, for instance, by writing down kinetic equations for two-point functions in the bosonic basis, and numerically studying their evolution after the quench. This kind of analysis could provide important information on the time scales of the departure from the pre-thermal state, which is a topic of current interest both for theory 50,64 and experiments 7 . Finally, since the interest towards prethermalization in long-range interacting models has increased in recent years, it would be tempting to apply the method presented here to the investigation of those cases 47,80-82 .

ACKNOWLEDGMENTS
We would like to thank E. Canovi for helpful comments on the numerical diagonalization of the model and P. Calabrese, M. Fabrizio and E. Tonni for useful discussions. J. M. acknowledges support from the Alexander von Humboldt foundation.
Appendix A: The quantum Ising chain in a transverse field The quantum Ising chain in a transverse field consists of a one-dimensional lattice of N sites, each accommodating a quantum 1/2-spin. These spins are simultaneously subject to a ferromagnetic nearest-neighbour interaction of strength J > 0, which favours configurations in which they are all aligned along a specific spatial direction (say, x), and an external magnetic field Jg directed orthogonally to it (e.g., towards z), which instead tends to destroy such an ordering. The corresponding Hamiltonian is 71 Periodic boundary conditions are chosen, such that σ x/y/z N +1 ≡ σ x/y/z 1 . Here σ µ i (µ = x, y, , z) denotes the standard spin operators acting on the i-th site. In the eigenbasis of σ z they can be conveniently represented as Pauli matrices: (A2) For simplicity, we set J = 1, which is tantamount to measuring energy and time in units of J and J −1 , respectively. By construction, operators acting on different sites commute: and for generic µ and ν. The Hamiltonian (A1) is invariant under the Z 2 transformation σ x i → − σ x i , σ z i → σ z i , which is implemented via the unitary operator U Z2 = i σ z i . This model undergoes a prototypical quantum phase transition 71 at the critical value g = g c = 1; for g > g c the system is paramagnetic and the longitudinal magnetization σ x i identically vanishes, whereas for g < g c a ferromagnetic ordering ensues which entails a spontaneous breaking of the Z 2 symmetry, i.e., σ x i = 0. The inherent integrability of the Hamiltonian in Eq. (A1) is made apparent after a Jordan-Wigner transformation accompanied by a Bogoliubov rotation 71 . The former reads 75,83 and is expressed in terms of fermionic creation and annihilation operators, c † i and c i respectively, with The Hamiltonian in this new basis reads with N = i c † i c i , and it is translationally-invariant (with anti-periodic and periodic boundary conditions in the even-and odd-parity sectors, respectively) and quadratic. It can be therefore mapped to a free-fermion model via a Bogoliubov rotation 71 (including a Fourier transform) with the new operator basis denoted by γ k , γ † k , and k the quasi-momentum running over discrete values The coefficients u k = cos θ k (g) and v k = sin θ k (g) are defined through the Bogoliubov angle θ k (g) given by 71 Therefore, also the operators γ k (g) and γ † k (g), depend on g and they keep preserved the canonical anticommutation relations We employ here the same convention as that one introduced in the main text, i.e., all quantities refer to the post-quench value g unless otherwise specified. An additional k ↔ −k symmetry of the Hamiltonian makes it possible to express it in the compact form where the summation is restricted to the positive values of k in Eq. (A8) and are Nambu spinors (here denotes transposition), while is the dispersion relation of the quasi-particles. Accordingly, the Hamiltonian (A11) is non-interacting: the various modes are independent and preserved under the unitary evolution, i.e., all the fermionic populationŝ commute with H 0 (g) and therefore they do not evolve. These operators constitute a sufficiently large set of conserved quantities to analytically solve the model. The simplicity of this solution might appear to be conflicting with the complexity of the collective modes driving a phase transition; however, the Jordan-Wigner transformation connecting the free-fermionic and the spin models is both non-linear and non-local and makes the correlation function of the order parameter a highly non-trivial combination of the free-fermion expectation values. The N constraints introduced by the conservation of the populations (A14) are actually much stringent than the original Z 2 symmetry mentioned above, as the latter can be implemented in this picture by the operator which describes the parity of the total number N p of fermions. In other words, U Z2 evaluates to 1 if N p is even and to −1 if odd. Note that the rightmost equality above comes from the fact thatn m k ≡n k for every integer m ≥ 1, due to the anticommutation relations in Eq. (A10).
remainder ∆M z ≡ M z − M z satisfies, in the thermodynamic limit, the cluster property at long times with ∆M z = M z − M z . We now analyze more closely Eq. (B1). We first consider the evolution of M z in the Heisenberg picture, where θ k (g) is the Bogoliubov angle. The time average in Eq.
(3) wipes out the oscillating factors, yielding and, correspondingly, On the other hand, the long-time limit of the the twotimes correlation function -without subtraction of M z , is a non-vanishing constant C (see, for instance Refs. 77 and 84), which is equal also to C = (M z ) 2 − M z 2 . The correlation function of M z after a quantum quench, spoils the time clustering property, because in the initial state there are non-vanishing correlations γ † k (g)γ † −k (g) = 0 and γ −k (g)γ k (g) = 0 between modes of opposite momenta. Remarkably, since the global nature of the operator M z implies that C scales with the system size, i.e., C ∼ O(L), the long-time clustering property is spoiled in the thermodynamic limit L → ∞.
Notice, for comparison, that C is also present in correlation functions of local operators, like the on-site transverse magnetization correlation function σ z i+r σ z i but there the contribution by C turns out to be subleading in the thermodynamic limit, i.e., ∼ O(L −1 ). Accordingly, the subtraction of M z , in the interaction term V removes a sort of undesired memory-preserving behaviour.

Appendix C: The Holstein-Primakoff transformation and its truncation
Hard-core bosons emerge naturally in the context of quantum spin-1/2 chains: it is in fact sufficient to introduce the spin ladder operatorsσ ± i as where σ µ i are the spin operators already encountered at the beginning of Appendix A, to yield the identification where indeed b † i and b i act as creation and annihilation operators of hard-core bosons. The Hilbert space can be accordingly reinterpreted by setting the correspondence |↑ i ↔ |1 , |↓ i ↔ |0 , in analogy with the notation introduced after Eq. (21): |0 denotes absence of the hard-core boson, whereas |1 the presence of a single one. The hard-core constraint can be reinterpreted as emerging from strongly repulsive contact interactions among the quasi-particles.
The Holstein-Primakoff transformation 68 expresses b i and b † i in terms of non-linear functions of the standard (free) bosonic operators a i and a † i . Here we shall focus, for simplicity, on a single mode b, b † , (by dropping the site index) which we recast in the form with a and a † obeying the usual commutation relations [a, a † ] = 1 and N = a † a. The Hilbert space is enlarged accordingly, from the two-dimensional space spanned by |0 and |1 , to the infinite-dimensional one generated by the usual bosonic number basis {|n } n∈N . On the other hand, the latter is split into two sectors which cannot be connected by b and b † , and which respectively include all the "physical" states {|n } n=0,1 and all the "unphysical" ones {|n } n>1 . This represents a relevant aspect, as the anticommutation relations are not correctly reproduced at the operatorial level; however, on both physical states (and, thus, in the whole physical subspace) one finds b, b † ≡ 1, thereby recovering the hard-core nature. When expanding the square roots in Eq. (C3) as power series, and approximating the b operators by truncation at any finite order, the separation between the physical and unphysical subspaces becomes weaker. This fact emerges quite clearly when considering the simplest possible case, i.e., b ( †) ≈ a ( †) ; indeed, within this approximation, we find in fact that b † |1 = |2 , which connects the physical state |1 with the unphysical one |2 . Consequently, the regime of validity of such an approximation is determined by the overlap of the state under study with the physical basis: the more it resembles its projection onto the physical space, the more accurate the result is.
We shall now briefly discuss the implications that the truncation has on the populationsn = b † b. In Sec. IV, we extensively used the approximationn ≈ N , which is valid only in the physical sector: in fact, if we restrict the states to just |0 and |1 , we can see that for every (integer) m = 0. Thus, as long as the system is still lying approximately in the physical space, we can approximate N m N , which rendersn N . Conversely, it is true that if N m N , i.e., for any integer m ≥ 2, then the truncation holds, as we prove below. Consider, in fact, a generic normalised state |ψ = n a n |n ; according to the discussion above, |ψ can be considered "approximately physical" as long as n>1 |a n | 2 1. Calculating the averages in Eq. (C6) on this state one can rewrite this condition as Using the equation above, it is easy to prove that the state, |ψ , which satisfies (C6) is indeed an "almost physical" one. From this discussion, we understand that, although the validity of the truncation may be intuitively expected to rely on a small quasi-particle density as suggested by the name "low-density approximation", it relies instead on higher moments of the quasi-particle densities, as shown by the condition (C6). On the other hand, N > 1 constitutes a clear sign that unphysical states are populated; therefore, it is still physically meaningful to focus on a low-density condition N 1.
We now make use of the fact that the system is quadratic for explicitly determining the temporal evolution of the operators η q and η † q : according to Eq. (30) we have where E q are the energies of the Bogoliubov modes determined numerically (see Eq. (30)). Consequently, the expectations in Eq. (D6) oscillate as where the matrices Z 0,q1q2 = η † q1 η q2 and Z 1,q1q2 = η q1 η q2 contain the information about the initial (t = 0) expectation values of all spin operators that can be written quadratically in terms of bosons a q (or equivalently η q ). We introduce the corresponding matrices where ∆θ k = θ k (g) − θ k (g 0 ) and θ k (g) is determined according to Eq. (A9). In the N × N block representation introduced in Eq. (D2) these matrices can be reorganised as where we used the properties W 0 = W 0 and W 1 = W † 1 which can be easily inferred from their explicit forms (D9). Exploiting the inverse change of basis in Eq. (D4) we finally find which allows an exact numerical calculation of the populations, in the sense described above. Unfortunately, this construction, which relies on the Heisenberg picture, has the disadvantage of being specific to the operator chosen; for example, for a quartic one in the operators a, it would be necessary to calculate every possible entry of the average a ⊗ a ⊗ a ⊗ a , which denotes a 4-tensor of dimension (N/2) 4 . On the other hand, once a specific tensor has been obtained, the corresponding dynamical expectation can be in principle calculated for any choice of the times by applying the formula Here U(t) keeps track of the evolution in the diagonal basis η, η † and can be written as and E is the diagonal matrix defined by 2E = diag E 1 , E 2 , . . . E N/2 , where {E i } i is the bosonic oneparticle spectrum. We emphasize that, for any choice of the time coordinates, the only operation left is the contraction of the indices k i in Eq. (D14), which just involves (N/2) m sums. Despite being of polynomial complexity, this can still pose serious difficulties for the investigation of large systems: for instance, the calculation of N 2 k for Fig. 5(a) is limited to N = 50 by the fact that, differently from the simpler populations n k , it involves a 4-index tensor.

Williamson's theorem
A symmetric, positive-definite, 2n × 2n matrix can be always brought into diagonal form by a symplectic transformation and the corresponding spectrum is positive and doubly-degenerate. This is the statement of Williamson's theorem 85 . The proof is constructive and shows how to translate the problem into one of standard diagonalization; since the algorithm we have employed follows its main steps 85 , we will report it here.
We start by recalling that a 2n × 2n matrix S is said to be symplectic (S ∈ Sp (2n, R)) if SΩS = Ω with Ω = −Ω = 0 1 n −1 n 0 , where 1 n is the n × n identity. As we have mentioned at the beginning of Appendix D, the Bogoliubov rotation (D2) defines in general a complex symplectic matrix, whereas here we assumed that it is real. In order to circumvent this complication, we employ the unitary transformation r = U a, ρ = U η, with U = 1 √ 2 which represents (apart from a multiplicative factor) a transformation from ladder operators a k , a † k to "position" and "momentum" operators Clearly, these new operators obey the canonical bosonic commutation relations [r k , p q ] = [ρ k , π q ] = iδ kq . As a shorthand, we group them in N -component vectors according to the basis: r = (r 1 , . . . , r N/2 , p 1 , . . . , p N/2 ) (starting basis) and ρ = (ρ 1 , . . . , ρ N/2 , π 1 , . . . , π N/2 ) (diagonal basis). These two vectors are therefore mapped one onto the other via the change of basis which is now evidently a real matrix. To see that M is symplectic we simply apply the definition (D16) to its form (D2), which yields where on the second line we have used the identities (D1). This implies that the matrix S preserves Ω = U ΩU = iΩ, which is tantamount to say that it is symplectic as well. We now proceed to show how this matrix can be determined. Using the canonical commutation relations, we rewrite the Hamiltonian H in Eq. (28) as where α and β are the n × n matrices defined in Eq. (20) with n = N/2 (we recall that we have assumed N to be even). The corresponding form in coordinate space (r, p) is Note that (β − α) k1k2 = δ k1k2 k1 (with k given in Eq. (10)), so that half of this matrix is diagonal and displays the unperturbed eigenvalues k ≥ |g − 1|, which makes it positive definite (if not at the critical point). The other half is given by (D24) and we can safely assume that, as long as g is kept far from g c = 1 and λ is not too large, also this part is positive definite and therefore Ξ satisfies all the requirements of the theorem. This implies that both the inverse Ξ −1 and its "square root" Ξ −1/2 exist and are symmetric, positive-definite matrices. We now define K = Ξ −1/2 ΩΞ −1/2 , which is skew-symmetric and invertible due to the properties of Ω (see Eq. (D16)). Accordingly, by the spectral theorem, there exists an orthogonal matrix R ∈ O (2n, R) which performs the block diagonalization where E −1 is a positive-definite, diagonal n × n matrix (which, as we are going to show, coincides with the one appearing in Eq. (D15)). Its positivity is guaranteed by the fact that one can always exchange a negative diagonal entry with its opposite lying in the opposite block −E −1 by exchanging the two vectors identified by the corresponding row and column via an orthogonal transformation. Being positive-definite, its inverse square root E 1/2 exists and we can use it to define the diagonal 2n × 2n matrices (D29) which, applying the transformation U † (·)U to retrieve the representation in terms of particle creation and annihilation operators and denoting with {E q /2} q the spec-trum of E, yields which corresponds to Eq. (30) with C = q>0 (E q /2 − q ).