Decay of long-lived oscillations after quantum quenches in gapped interacting quantum systems

The presence of long-lived oscillations in the expectation values of local observables after quantum quenches has recently attracted considerable attention in relation to weak ergodicity breaking. Here we focus on an alternative mechanism that gives rise to such oscillations in a class of systems that support kinematically protected gapped excitations at zero temperature. An open question in this context is whether such oscillations will ultimately decay. We resolve this issue by considering spin models that can be mapped to systems of weakly interacting fermions, which in turn are amenable to an analysis by standard methods based on the BBGKY hierarchy. We find that there is a time scale beyond which the oscillations start to decay, which grows as the strength of the quench is made small.


INTRODUCTION
A key question in non-equilibrium many-body quantum dynamics is to understand how ergodicity can be broken and thermalization avoided.The best known examples are many-body localization (MBL) [1] and quantum integrable systems [1,2].The issue of ergodicity in generic systems is addressed by the Eigenstate Thermalization Hypothesis [3][4][5] (ETH), which provides a description of the structure of the matrix elements of observables in energy eigenstates.In the literature, violations of ergodicity are often referred to as either 'strong' or 'weak' depending on whether the fraction of states failing to satisfy ETH remains non-zero or vanishes in the thermodynamic limit, respectively [6].Weak violations of ergodicity have recently attracted a great deal of attention in the context of so-called quantum many-body scars (QMBS) [7][8][9][10][11].This notion was first introduced to explain the unexpected dynamics of a Rydberg atom quantum simulator [12] where initializing the system in a particular initial state resulted in long-lived oscillations in the time-evolution of observables.
A seemingly related phenomenon has been observed in quenches in certain spin chains, referred to as 'weak thermalization' in [13].Various authors [13][14][15][16][17][18][19][20][21][22][23] have noted that for the Ising chain in a tilted field following a quench, there are long-lived oscillations at frequencies corresponding to the masses of 'meson' bound states.Similar oscillations were observed in the Potts model [24].These studies have been unable to simulate late enough times to conclude if the meson oscillations decayed at late times.In the paramagnetic regime [14] found oscillations beginning to damp following a quantum quench in the scaling limit.Theoretical arguments for the eventual decay have been given in [16].In the latter work it is moreover argued that such oscillations should be generic to quantum systems with a quasiparticle gap and isolated bands such as produced by bound states.
The fate of oscillations at late times is however controversial: oscillations in the post-quench dynamics of observables in quantum field theories were predicted by Delfino and collaborators in [25,26] and subsequently in [27,28] it was argued that they persist for arbitrarily long times.
The Ising chain in a tilted field is unusual in the sense that the perturbation that leads to the formation of meson bound states is non-local with respect to the fermionic elementary excitations of the transverse-field Ising chain.This precludes the analysis of particle decay by standard perturbative approaches.In light of this fact it is important to identify models that exhibit the same phenomenology, but can be studied by such methods.One such example was recently reported by us in the axial next-nearest neighbour Ising model (ANNNI) [29], which has a non-confining (contact) potential for the elementary fermion excitations.In this work we give two quench setups that we consider to be particularly simple examples of such behaviour.The first is a quantum quench in integer spin antiferromagnetic chains, which possess a gapped single particle (magnon) mode according to Haldane's conjecture [30,31].This model very cleanly demonstrates that the phenomenology is due neither to confinement, nor indeed to bound states, but simply due to the system possessing a kinematically protected gapped quasiparticle excitation.However, as integer spin antiferromagnets are strongly interacting systems we are restricted to purely numerical investigations of the non-equilibrium dynamics of this model by matrix product state methods.Secondly, we present and analyze in detail a novel example of these long-lived oscillations in a model with weak interactions that has the simplifying feature of exhibiting a U(1) symmetry associated with particle number conservation: a dimerized XXZ chain in a staggered magnetic field.In the scaling limit, the low-lying excitations of this model can be understood in terms of solitons, anti-solitons and a bound state known as a 'breather' which can give rise to longlived oscillations after quantum quenches.
Many other mechanisms for producing long-lived oscillations are possible in quantum systems that should be arXiv:2307.04466v2[cond-mat.stat-mech]13 Aug 2024 differentiated from the case we discuss.We have already mentioned exact quantum many-body scars, which can cause infinite lifetime oscillations if the initial state has large overlap with scar states lying in the middle of the spectrum.There are also Bloch oscillations of domains in the tilted field Ising model at low domain-wall density and related models of Rydberg atoms [32,33].Long-lived oscillations can also occur in the electric field strength in lattice gauge theories that are related to quantum manybody scars [34,35].Finally, in the presence of a U(1) symmetry undamped oscillations can occur when one considers observables that connect neighbouring charge sectors and applies a Zeeman field that splits all sectors by the same energy difference [36,37].
The organization of this paper is as follows: in Section II we describe the mechanism that can give rise to long-lived oscillations in models with kinematically protected single-particle excitations before giving a simple numerical case study of such behaviour in the spin-1 bilinear biquadratic (BLBQ) chain in Section II A. The spin-1 chain provides a clear example of the phenomenology we are interested in, however it is always strongly interacting.In contrast, the ANNNI and related Ising models can be mapped to fermionic chains for which the interaction can be considered perturbatively.However the lack of a U(1) symmetry (and in the case of a tilted field, long-ranged interactions between fermions) significantly complicates the application of standard methods based on the BBGKY hierarchy [38][39][40].In Section II B we address this problem by introducing a dimerized XXZ chain in a staggered magnetic field, which exhibits a U(1) symmetry as well as long-lived oscillations of observables after quantum quenches.In Section III we present two different approximations based on the BBGKY hierarchy -a self-consistent time-dependent mean-field theory (SCTDMFT) and the Second Born approximation [41,42] -to study the quench dynamics of local observables.

II. OSCILLATIONS AT "EARLY" TIMES
The physics underlying the oscillations explored within this paper arises from two key requirements, these are: 1.The existence of a kinematically protected mode, i.e. a quasiparticle with a gapped dispersion ϵ(q) ≥ ∆ ex such that there is some region in the energymomentum plane that has no other energy eigenstates.
2. Initial density matrices ρ(t = 0) such that the energy density deposited into the system by the quench (relative to the post-quench ground state energy E GS ) is small compared to the spectral gap of the post-quench Hamiltonian If the above requirements are met the physics can be viewed in terms of a dilute gas of long-lived particles, whose scattering is to a good approximation purely elastic, cf.[43][44][45][46][47][48][49][50].The possible emergence of long-lived oscillations can then be understood by considering the linear response regime of ground state quenches.This is equivalent to the approach of Refs [25][26][27][28].To that end we consider an initial state |ψ 0 ⟩ that is the ground state of a Hamiltonian H 0 , and time-evolve with H = H 0 +λV .Here V is a global operator that is assumed to be translationally invariant (as is H 0 ).By construction both H 0 and H feature kinematically protected gapped singleparticle excitations.Linear response theory then gives The response function χ OV (t, t ′ ) can be expressed in terms of a Lehmann representation using the eigenstates of H 0 , which gives (3) We now consider perturbations V and operators O that have non-vanishing matrix elements between the ground state and the kinematically protected gapped excitation.As by construction we are dealing with a ground state calculation the contribution of this excited state will provide the dominant contribution to the linear response function for large t where F O (k) = ⟨ψ 0 |O|k⟩ is the matrix element of the operator O between the ground state of H 0 and the single quasiparticle excitation of H 0 with momentum k and dispersion ε(k).As V is by assumption a global, translationally invariant operator the only non-zero matrix element is with the zero-momentum single-particle state which establishes that in linear response theory we obtain persistent oscillations with frequency ε(0), cf.Refs [25][26][27][28].If instead either V or the initial state has invariance only under translations by two sites, the matrix element will select out momenta k = 0 and k = π.In this case oscillations will occur at two frequencies, ε(0), ε(π), as long as the kinematically protected mode exists at these momenta.The approach of [25][26][27][28] can be straightforwardly modified by expanding the initial state in terms of the eigenstates of the post-quench Hamiltonian using perturbation theory.This gives Here ϵ(k) is the dispersion of the kinematically protected mode in H (rather than H 0 ) and F O (k) = ⟨0|O| k⟩ is the matrix element of the operator O between the ground state of H and the single quasiparticle excitation of H with momentum k.This way of approaching the problem is important for some of the cases considered below, in which H has a kinematically protected single-particle mode, but H 0 does not.In these cases the small perturbation λV leads to the creation of a bound state, which is a non-perturbative effect.In these cases it turns out that the perturbation theory around H as sketched above gives a (qualitatively) correct description of the observed dynamics.
The question we want to address is what happens outside the linear response regime.Linear response theory is usually expected to describe the short-time regime for very small but finite values of λ, but fail at late times.The question of its regime of applicability is related to the properties of nonlinear response functions, which have recently been analyzed in the class of systems discussed here [51] and shown to acquire late-time divergences in some cases.
The linear response viewpoint summarized above obscures the fact that the quench deposits a finite energy density into the system.A complementary viewpoint on long-lived oscillations is obtained by employing a spectral representation in terms of energy eigenstates of the post-quench Hamiltonian Tr Oρ(t) = m,n e i(Em−En)t ⟨m|O|n⟩⟨n|ρ(t = 0)|m⟩ .(7) Assuming that the operator O connects states with quasiparticle numbers that differ by one, oscillations with frequency ∆ ex may ensue for the following reason.In the gas phase energy eigenstates can be viewed as scattering states of the stable quasiparticles and adding a single quasiparticle with momentum q to an energy eigenstate (approximately) leads to an energy eigenstate that differs in energy and momentum by ϵ(q).If ρ is translationally invariant, then the only nonzero contributions to the sum occur at q = 0.As this process works for all energy eigenstates at the (low) energy density of interest (which is set by the factor ⟨n|ρ(t = 0)|m⟩), one may expect long-lived oscillations in the expectation value (7) to occur.
We stress that the oscillations produced by the mechanism outlined above are not a finite-size effect but can persist in the thermodynamic limit.For all examples discussed below we have verified that varying the system size does not affect the amplitude or frequency of the oscillations observed.This is very different to the oscillations reported in [52], which indeed are finite-size effects.

A. Haldane-gap chains
As a first example of a model that exhibits undamped oscillations after quantum quenches in the linear response regime we consider the antiferromagnetic bilinearbiquadratic (BLBQ) chain.This is a family of spin-1 chains described by [53][54][55] For γ = 0 the model reduces to the spin-1 Heisenberg antiferromagnet whilst for γ = 1/3 it is the AKLT chain [56], whose ground state is an exact MPS with bond dimension χ = 2.Both values of γ lie within the gapped 'Haldane gap phase' [31] −1 < γ < 1.At γ = 1 the model is the SU (3) symmetric Lai-Sutherland model which is gapless [57,58].In Fig. 1 we plot the low-energy spectrum obtained by exact diagonalization (ED) on L = 16 sites and periodic boundary conditions for γ = 0.25 which is representative of the Haldane phase, with a gap of ∆ ex (π) ≈ 0.62J and group velocity v ≈ 1.26J.
The ground state is at k = 0 and there is a gap to a triplet band of magnons with energy minimum at k = π.
Consider first global quenches where a finite energy density is generated by quenching the ratio of exchange constants γ → γ ′ .For the Hamiltonian (8) the ground state has momentum k = 0 but the Haldane gap is at k = π, with the magnon mode only persisting in a region of the Brillouin zone around this that does not extend to k = 0.For such a quench both the Hamiltonian and the initial state are translationally invariant, so local operators cannot "access" single magnon excitations in the way Result of a quantum quench keeping γ = 0.25 and quenching the initial staggered field hs : 0.01 → 0. This quench produces an energy per site ϵ Quench ≈ 3.8 × 10 −4 J which corresponds to an equilibrium temperature of T ≈ 0.12J (estimated using ED on 12 sites).TEBD performed using L = 400 and χ = 400.The final time is window is determined by requiring that the results do not change on increasing the bond dimension to χ = 600.
described above because ρ nm = ⟨n|ρ(t = 0)|m⟩ = 0 for energy eigenstates that differ by a single magnon, as they differ in momentum.As a result there are no long-lived oscillations for translationally invariant initial states.On the other hand these considerations suggest a way out: we need to choose an initial state that is invariant only under translations by two sites, which can be achieved simply by choosing the pre-quench Hamiltonian to have an additional staggered magnetic field H pre now has a ground state that is invariant only under translation by two sites and ρ nm ̸ = 0 and so as discussed in the previous section we therefore expect to see oscillations at ϵ(π) (the magnon does not extend to k = 0 so this will be the only frequency present).As a numerical test of these ideas in the BLBQ chain we perform a quench using the ITensor [61] library, which enables us to use DMRG to compute an approximation to the ground state of H(γ i , h s ), and then to time-evolve the state according to H(γ f , 0) using time-evolving block decimation (TEBD).Here and elsewhere, the time window in which we plot TEBD data is determined such that the TEBD results do not change when suitably increasing the bond dimension, in Figs. is to be understood as due to a discrete spin-flip symmetry as follows: the oscillations can only occur when the matrix element ⟨GS|O|QP(k)⟩ ̸ = 0 where |QP(k)⟩ is the quasiparticle at momentum k and |GS⟩ is the ground state of the post-quench Hamiltonian.For the BLBQ chain the ground state is invariant under S z → −S z but the magnon mode is odd under this symmetry.Therefore, the matrix elements are only non-zero when the observable is Z 2 odd, such as S z itself, and conversely the oscillations decay rapidly when the observable is Z 2 even.The oscillations in the magnetization have frequency ω very close to the magnon gap at k = π as expected Fig. 3 shows a somewhat larger quench from h s = 0.05, which still has an energy per site well below the gap.The average inter-particle distance after this quench is approximately From the perspective of a low density gas of quasiparticles, the finite lifetime of the oscillations is caused by scattering events [23].Therefore one would not expect to see appreciable decay at times vt ≲ ℓ.From ED we estimate that the group velocity is roughly v ≈ 1.26J, and so the slight decay by Jt = 50 makes sense within the quasiparticle gas picture.Conversely, for the quench in Fig. 2 the mean free path is ℓ ≈ 1600 and thus the lack of decay is also consistent with this rough estimate.
For even larger quenches than in Fig. 3 we find that, as expected on the basis of our quasiparticle gas picture, the decay becomes more easily visible at short times.The quenches in Haldane-gapped models explored above yield several insights -the first is that for weak quenches oscillatory behaviour results as predicted [25][26][27][28] when there is a quasiparticle mode.This confirms that such oscillations are unrelated to the formation of bound states or of confinement, except inasmuch as they provide a mechanism for kinematically protected quasiparticles to exist.It also highlights the importance of symmetries which can cause the relevant matrix elements to be zero and relaxation to be consequently much faster.Perhaps most importantly we have evidence that the oscillations do decay, which was suggested in previous studies [17,19] but not observed in the models considered therein.Our findings also show that in the case studied above the regime of validity of the perturbative approach used in [25][26][27][28] is limited to short times.

B. Dimerized XXZ Model
Generally the quasiparticle gas can consist of several species, for instance in models with "elementary" quasiparticle excitations as well as (multi-particle) bound states.We have already mentioned two examples of this situation -the Ising model with both longitudinal and transverse field [15,17,19,26] and the Ising model with transverse field only, but additional next-nearest neighbour Ising interactions [29].Neither of these lattice models analyzed in the context of persistent oscillations exhibits a U(1) symmetry associated with particle number conservation, which greatly complicates the application of perturbative approaches based on the BBGKY hierarchy or the flow equation approach [41,42,[62][63][64][65][66][67].In order to study the fate of these oscillations at very late times we therefore introduce a spin-1/2 dimerized XXZ model in a staggered field, which can be mapped to a model of spinless interacting fermions with particle number conservation.This enables us to apply the equations of motion techniques developed in [41,42].This model features elementary fermionic excitations as well as bosonic two-particle bound states.Moreover, in the appropriate scaling limit the model reduces to the sine-Gordon quantum field theory in the attractive regime.The Hamiltonian of our model is Here α tunes the degree of dimerization in the XY plane and h s is a staggered applied field.For α = h s = 0 and |∆| < 1 the model reduces to the integrable spin-1/2 XXZ chain in the massless Luttinger liquid phase, whilst non-zero values of α, h s break the integrability and open a gap in the zero magnetization sector.A Jordan-Wigner transformation maps the Hamiltonian (12) to one of interacting spinless fermions with interaction strength ∆.As in the case of the BLBQ chain, discrete symmetries play a role in allowing or disallowing persistent oscillations.When only one of α, h s is present in (12) there is a discrete spatial Z 2 symmetry corresponding to reflection across a bond or site, respectively.We will elaborate the importance of retaining both parameters after presenting data from quenches.
Rotational symmetry about the z axis corresponds to U(1) particle number conservation in the fermionic variables.In the low-energy limit, the Hamiltonian (12) for |∆| < 1 reduces to a sine-Gordon model [68], whose low lying excitations for ∆ > 0 are solitons, anti-solitons and soliton-antisoliton bound states known as 'breathers'.In the lattice model we determine the spectrum of low-lying excitations in the S z Tot = 0 sector by exact diagonalization, cf.Fig. 4. We can see that throughout the Brillouin zone there is a bound state visible below a continuum of states.The bound state is kinematically protected and therefore stable.This establishes that our model fulfils the first of our requirements.
We investigate the time evolution using both TEBD, which is capable of exactly describing the evolution of states with sufficiently low entanglement, and perturbative methods which are valid at small ∆.The latter is needed because following a quench the entanglement entropy generically grows linearly in time [69][70][71].As such the true time evolved state quickly leaves the manifold of states that can be accurately described by matrix-product states with finite bond dimensions.We adopt open boundary conditions when performing TEBD numerics; when working with the equivalent fermionic model it suffices to work in the sector with fixed fermion number and we adopt periodic boundary conditions out of convenience.We also ensure that our system sizes are sufficiently large to rule out finite-size effects such as traversals [2] on the time scales we are interested in.
To ensure a long window of applicability of the perturbative approaches [41,42] we consider quantum quenches from an initial thermal state of the non-interacting model ρ(0, α, h s , β), where As the perturbative approaches expand around thermal states of free Hamiltonians, they are able to describe states with volume law entanglement, unlike MPS methods.Instead, they are limited by their assumption that higher particle cumulants are negligible.This is then time evolved using the Schrödinger equation for H(∆, α ′ , h ′ s ).We focus on expectation values of local observables such as the staggered magnetization within a unit cell and nearest-neighbour spin-bilinears We plot these quantities in Fig. 5 computed using TEBD on L = 400 sites with the ITensor [61] library, and two approximate calculations, self-consistent timedependent mean field theory and the second Born Approximation, both detailed in Section III.Fig. 5 shows the case where the time evolution has a Z 2 symmetry (reflection in a bond as h ′ s = 0) but the initial state does not.In this case, the staggered magnetization initially has a non-zero value and then oscillates about the thermal value of 0. The frequency agrees with the energy difference between the post-quench ground state and the bound state, which have opposite Z 2 parities and are thus connected by the Z 2 odd operator m z s .Conversely, the operator S xx m,m+1 is Z 2 even and thus has decaying expectation value at late times.We note that in the bottom panel of Fig. 5 the SCTDMFT does not appear to accurately capture the early time evolution.This is despite the fact that SCTDMFT is expected to be accurate at early times, as will be explained in Section III A. We have checked that the difference between it and the second Born approximation scales like O(∆ 2 ) and so this deviation indicates that for the observable in question the O(∆) contribution is either very small or absent.

III. DECAY OF OSCILLATIONS IN THE STAGGERED XXZ MODEL AT LATE TIMES
We now turn to the "intermediate" time regime.This is no longer accessible to TEBD for the reasons set out above, but can be studied by appropriate truncation schemes of the BBGKY hierarchy.We first consider the simplest such scheme, a self-consistent time-dependent mean field theory [29,[72][73][74][75][76][77].

A. Self-consistent time-dependent mean-field theory (SCTDMFT)
The correct degrees of freedom for constructing a mean-field approximation for (12) are fermions rather than spins.Applying a Jordan-Wigner transformation to the Hamiltonian gives where c n are spinless fermions and n m = c † m c m .If the fermion parity is odd the c n have periodic (Ramond) boundary conditions whilst if the fermion parity is even they have anti-periodic (Neveu-Schwarz) boundary conditions.We work at half filling such that these conditions correspond to the number of unit cells being odd/even respectively.In SCTDMFT the interaction terms in Eq. (15) are decoupled in a time-dependent way, which corresponds to normal-ordering with respect to the timeevolving state of the system and retaining only terms that are quadratic in fermionic creation/annihilation operators Here the expectation values ⟨.⟩ t in the time-evolving state are determined self-consistently.This produces a meanfield Hamiltonian H MFT (t) with the following key properties: 1. H MFT (t) is quadratic at all times and thus time evolving a Gaussian state with H MFT (t) ensures that Wick's theorem holds at all times.
2. The equations of motion for the two-point functions obtained using H MFT agree with the equations of motion obtained using the full Hamiltonian, under the approximation that Wick's theorem is valid.
Since our initial state is a thermal state of the free fermion Hamiltonian H(0, α, h s ) it is Gaussian and the SCTDMFT is thus accurate at early times by construction.The time-dependent mean-field Hamiltonian takes the form where s now labels the unit cell and a = 0, 1 the sites within it such that the spin labelled (s, a) is at position m = 2s + a in the chain.The constant term E 0 (t) does not have any effect on the equations of motion, but ensures that the expectation value of energy is conserved.The mean-field Hamiltonian contains the effective couplings Note that J 0,1 (t) are generically complex at intermediate times and that h eff is (up to a constant shift and rescaling) equal to the staggered magnetization m s (t).The mean-field Hamiltonian can be block-diagonalized by the canonical transformation where k = 2πn/L for n = 0, . . .L/2 − 1.We find Here A(k, t), B(k, t) and B ′ (k, t) are real functions that can be expressed in terms of J ± (t) = J 0 (t) ± J 1 (t) and h eff (t) as The time evolution of the momentum space two-point functions easily obtained from the Heisenberg equations of motion These equations can alternatively be derived by a first order truncation of the BBGKY hierarchy, cf.Ref.We solve these equations of motion numerically, updating the mean fields J ± (t) and h eff (t) every timestep and plot the resulting staggered magnetization following a quench from an initial thermal state in Fig. 6, which shows clear oscillations that become highly monochromatic and undamped at late times.The choice of thermal state is made to ensure that the system has an energy per site well above the ground state but still small compared to the gap.This ensures that there is a low density of quasiparticles in the system and that they can be treated as a dilute gas.
We now return to the physical origin of the oscillations and their frequency.To that end we have considered ground state quenches at α = 0.4, h s = 0.3, i.e. initial density matrices ρ(0, 0.4, 0.3, ∞), for several ∆.In this case we observe essentially a single oscillation frequency ω B at intermediate and late times.For example, in Fig. 7 we show the evolution of the mean fields.Performing a fast Fourier transform using data from t = 50 up to t = 2000 gives a single sharp peak at the frequency ω B .We compare this to the energy of the first excited state computed by exact diagonalization on system sizes up to L = 30 in Fig. 8.We observe that the oscillation frequency observed in SCTDMFT is in very good agreement with the bound state gap at q = 0 for small interaction strengths ∆ ≲ 0.25J.For small interactions strengths ∆ ≈ 0 the ED results exhibit sizeable finite size effects that are discussed in Appendix A. Using the extrapolation procedure summarized there provides us with the curve labelled 'ED (Extrapolated)' in Fig. 8.
The emergence of persistent oscillations of the expectation values of observables in the SCTDMFT can be understood as follows.The solutions to the self-consistency equations reported above exhibit periodic behaviour with a single frequency.As a result the SCTDMFT is equivalent to a periodically driven system with a Hamiltonian that is quadratic in fermions.It is well known that such systems typically synchronize at late times and physical observables then display persistent oscillations at the driving frequency [78].

B. Second Born approximation
The lack of damping in SCTDMFT is in fact not surprising as the method is perturbative to first order in ∆ (at the level of the equations of motion).In thermal equilibrium we have to evaluate the self-energy to second order in ∆ in order to obtain a non-vanishing imaginary part that signals a finite life-time of the fermions.This suggests that finite lifetime effects in the non-equilibrium setting of interest here can be captured by the 'second Born approximation' [41,42,65].We follow [41] in deriving the equations of motion for fermionic bilinears nµν chosen to diagonalize the quadratic part of the Hamiltonian where we have introduced shorthand notations k The equations of motion to second order in ∆ for n µν = ⟨n µν (k, t)⟩ are obtained by truncating the BBGKY hierarchy as derived in [41,42].The result is V µ1µ2µ3µ (k, q, q, k)e i(ϵµ 1 ν (k)+ϵµ 2 µ 3 (q)t n µ1ν (k, 0)n µ2µ3 (q, 0) −4i∆ V νµ2µ3µ1 (k, q, q, k)e i(ϵµµ 1 (k)+ϵµ 2 µ 3 (q))t n µµ1 (k, 0)n µ2µ3 (q, 0) where the kernel functions L µν , K µν are given by and where The derivation of Eq. ( 25) is summarized in Appendix B. We note that n µν are different from the quantities n ±± (k, t) considered in the SCTDMFT of the previous section.Taking this into account one sees that the SCT-DMFT agrees with (25) up to order O(∆ 1 ) and disagrees with the O(∆ 2 ) terms as expected.Solving Eq. ( 25) requires a runtime of O(L 3 × T ) where T is the simulation time reached.Since simulating up to time T requires a system size at least 2v LR T where v LR is the Lieb-Robinson velocity, investigating up to time T scales as O(T 4 ).
The second Born approximation is premised on the assumption that many-particle cumulants are small at intermediate times, whilst going beyond SCTDMFT by allowing for a non-Gaussian state.A priori, this is an uncontrolled approximation.However, we start in a Gaussian state in which all cumulants vanish and so the approximation must be accurate for early times, and becomes better as the interaction strength ∆ becomes small.Furthermore, it gives rise to a Boltzmann equation at late times [41,42] which is believed to become exact in the "Boltzmann scaling limit".This suggests that the approximation remains accurate on time scales t ∼ ∆ −2 .At intermediate times one expects the second Born approximation to continue to provide useful physical insights even if it may not retain full quantitative accuracy.What we wish to establish in this paper is that whilst the SCTDMFT treatment agrees with the prediction of [25][26][27][28], the leading correction provided by the second Born approximation causes the oscillations to damp.Finally, we note that the second Born approximation is complementary to TEBD in such cases, since the former can reproduce volume law entanglement but cannot fully capture strong interactions, whilst the latter method can only describe states with a finite amount of entanglement related to the bond dimension used but can describe strongly correlated states at sufficiently low entanglement.
In order to clearly exhibit some of the issues associated with the damping of oscillations we first consider ground state quenches, in which we initialize the system in the ground state of H(0, 0.4, 0.3) and time-evolve with H(0.2, 0.4, 0.3) for a system size of L = 300.Fig. 9 shows the time evolution of the staggered magnetization and observe long-lived oscillations with no apparent damping on the long time scales considered.This is however entirely corresponds to a very small quench with quasiparticle density ∼ 1.6 × 10 −3 .Accordingly very large system sizes and late times would be required to observe decay of oscillations.Inset: half chain entanglement entropy, which grows linearly for all times shown.TEBD calculation done with maximum bond dimension χ = 800.
expected as the quench produces a very small energy per site ϵ quench ≈ 0.0019J whilst the gap to create a quasiparticle is ∆ ex ≈ 1.18J.The resulting average interparticle distance is therefore ℓ = ∆ ex /ϵ Quench ≈ 640.That this exceeds the system size simulated means that finite size effects such as traversals will matter long before the many-body effects that would dampen the oscillations.The fact that we are effectively dealing with the linear response regime is also apparent from the fact that TEBD is able to access very large time scales Jt ∼ 100, which means that the volume-law contribution to the entanglement entropy is still negligible.
These considerations show that the energy per site deposited by the quench should be small compared to the bound state energy, however it should not be so small that quasiparticle interactions are negligible on accessible time scales.To overcome this problem we consider larger quenches of the interaction parameter as well as thermal initial states, which provide us with a simple parameter -the pre-quench temperature -to vary the post-quench energy density.For finite pre-quench temperatures we cannot use TEBD since the initial state has volume law entanglement, and so only show the second Born and SCTDMFT results.We compare the results obtained by the second Born approximation to SCTDMFT, which as discussed before exhibits persistent oscillations at a frequency that is very close to the bound state energy gap.In Fig. 10 we show results for quench initial- Staggered magnetization for a quench from the ground state of H(0, 0, 0.23) and time evolved with H(0.2, 0.4, 0.3) using mean-field theory, the second Born approximation and a TEBD calculation with maximum bond dimension χ = 800 and L = 200.ized in the ground state of H(0, 0, 0.23) and time evolved with H(0.2, 0.4, 0.3).Here the post-quench energy per site is ϵ Quench ≈ 0.040 corresponding to a mean-free path ℓ ≈ 29, which is much smaller than our system size of L = 200.We observe that the second Born approximation clearly shows the decay of the oscillatory behaviour of the staggered magnetization.We plot the TEBD results only up to times Jt = 40, where our criterion is agreement of the numerical results for bond dimensions χ = 600 and our maximal bond dimension χ max = 800.The TEBD data show the beginning of a decay, consistent with the second Born results.In Figs 11 and 12 we show the behaviour of the staggered magnetization after quenches from thermal initial states.In Fig. 11  3).L = 448 used for the second Born approximation.The initial state is chosen such that the energy density is approximately the same as in Fig. 10.Inset: successive peak to peak amplitudes of the oscillations in the second Born approximation, with an exponential fit (dashed black line).The grey scatter points are excluded from the fit.we initialize the system in the thermal state of the noninteracting system at inverse temperature Jβ i = 4, which corresponds to the same energy density as in Fig. 10.We again observe decaying oscillations.In the inset in that figure, we estimate the decay time by fitting a decaying exponential A exp(−t/τ ) to the successive peak to peak amplitudes -the resulting decay time for this particular quench is Jτ ∼ 200.Finally, in Fig. 12, we consider a lower temperature Jβ i = 10, which corresponds to energy per site ϵ Quench ≈ 0.0037 and mean free path ℓ ≈ 250.
Here the oscillations are seen to decay very slowly.12. Staggered magnetization for a quench starting in the thermal state ρ(0, 0.4, 0.10, 10.0) and time evolved with H(0.2, 0.4, 0) on a ring with L = 400.Other than the finite pre-quench temperature Jβi = 10.0 this is the same as Fig. 5.

IV. SUMMARY AND CONCLUSIONS
In this work we have carried out a detailed study of a mechanism that gives rise to long-lived oscillations in the expectation values of local observables after quantum quenches.This mechanism is very different from quantum scars and occurs after small quenches in interacting many-particle systems with an excitation gap, which generate a regime that can be understood in terms of a low-density gas of (long-lived) kinematically protected quasiparticles.Long-lived oscillations can then occur in expectation values of observables that have matrix elements between the ground state and excited states that contain a single quasiparticle.
We have presented very strong evidence using a combination of matrix-product state methods and perturbative approaches based on truncations of the BBGKY hierarchy that these oscillations decay at late times in all models we have considered.This is an important difference to models with exact quantum scars [8,10].
Our results show that the linear response theory prediction is upheld only at the level of self-consistent meanfield theory.Going beyond mean-field theory to the second Born approximation provides evidence of damping.For small interaction strengths U the timescale of the decay is therefore generally expected to be O(U −2 ).Whilst it is not impossible that higher order corrections would cause the oscillations to remain at late times, this is highly unlikely as there is no reason to anticipate such an effect.Instead, truncating the BBGKY hierarchy at higher orders should merely modify the lifetime.We presented non-perturbative numerics using TEBD for the oscillations in both the spin-1 chain (Fig. 3) and dimerized XXZ chain (Fig. 10) which indicate our qualitative conclusions that the oscillations damp at intermediate times are robust to including higher orders.
Our results differ from the prediction made in [27,28] that the oscillations have infinite lifetime regardless of the quench strength λ.Whilst those papers are formulated in the continuum, the same arguments made therein lead to oscillations on the lattice that we have shown to decay.
Our perturbative analysis generalizes straightforwardly to other interacting fermion and boson models with interactions that are local with regards to the elementary excitations of the unperturbed theory.We note that the much-studied Ising chain in a tilted field [13, 15-19, 21, 23] does not fall within this class of models.This is because when viewed as a perturbation of the transverse-field Ising chain, which maps onto noninteracting fermions by the Jordan-Wigner transformation, the perturbing operator is not local in terms of these fermions, i.e. involves interaction vertices with arbitrarily large numbers of particles.This precludes employing approaches based on truncating the BBGKY hierarchy for the fermionic degrees of freedom.FIG. 13.Dispersion at ∆ = 0, with the ground state indicated.Filled circles indicate states that are occupied and empty circles ones that are unoccupied, drawn at L = 16 for clarity.Note that k = π/2 is not an allowed single particle state for any finite system size.
converged with respect to system size when estimating the gap for small ∆.In the thermodynamic limit the minimum energy excitation corresponds to a hole in the filled band at k = π/2 and a particle in the empty upper band at k = π/2, with a corresponding gap of E gap (∞) = 2ϵ + (π/2) = 2 α 2 J 2 + h 2 s .However, if the number of unit cells L/2 is even then we must work in the Neveu-Schwarz sector and so have anti-periodic boundary conditions with k = (2n + 1)π/L, which never equals exactly π/2.Likewise for an odd number of unit cells we work in the Ramond sector where k = 2nπ/L ̸ = π/2.For finite L the gap is therefore This functional form indeed provides a good description of our numerical results for the gap deduced from ED on L ∈ {20, 22, 24, 26, 28, 30} sites when ∆ is small.

Appendix B: Equations of motion in Second Born approximation
The method we use for deriving the equations of motion ( 25) is given in more detail in [41,42], here we simply briefly recap the main points for completeness.The first step is to diagonalize the free part of the Hamiltonian using Eqs.(A2)-(A4).
The next step is to rewrite the interaction in this basis.Whilst in the original real-space basis the interaction is the same as that considered in [41,42], we find a different final form since we obtain different b µ (k) of the free part of the Hamiltonian, which agrees with the expression there when setting h s = 0. We antisymmetrize the interaction in the first and second pairs of indices sgn(P )V ′ P (µ) (P (k)) . (B1) Here P is an element of Z 2 × Z 2 where the first Z 2 swaps µ 1 ↔ µ 2 , k 1 ↔ k 2 and the second Z 2 factor acts likewise on 3, 4. We define sgn(P ) as the product of the sign of each permutation.The unsymmetrized interaction components are equal to where we have defined The Heisenberg equations of motion are ∂ nµν (k, t) ∂t = iϵ µν nµν + i∆ Y µ µν (k, q) Âµ (q) ,(B4) where ϵ µν (k) = ϵ µ (k) − ϵ ν (k) and the coefficients of the quartic operators appearing are given explicitly in terms of the interaction potential V µ (k) by Y µ αβ (k, q) = δ βµ4 δ k,q4 V µ1µ2µ3α (q) + δ βµ3 δ k,q3 V µ1µ2αµ4 (q) − δ αµ2 δ k,q2 V µ1βµ3µ4 (q) − δ αµ1 δ k,q1 V βµ2µ3µ4 (q) .(B5) The quartic operators Âµ (k) likewise evolve according to the following Heisenberg equations of motion V γ (q) Âγ (q, t), Âµ (k, t) , where E µ (k) = ϵ µ1 (k 1 )+ϵ µ2 (k 2 )−ϵ µ3 (k 3 )−ϵ µ4 (k 4 ).This equation contains products of up to six fermionic operators on the right hand side, and carrying on in this way generates the BBGKY hierarchy of equations of motion.

FIG. 5 .
FIG.5.Time evolution following a quench from the ground state with ∆ : 0 → 0.2, hs : 0.1 → 0 and αs = 0.4 before and after the quench.Upper: staggered magnetization, showing persistent oscillations.Lower: Time evolution of ⟨S x m S x m+1 ⟩ where m = L/2.Calculations use L = 400 sites and a maximum bond dimension of χ = 1000 for the TEBD.

FIG. 6 .
FIG.6.Mean-field evolution of the staggered magnetization after a quench with initial state ρ(0, 0.4, 0.2, 4.0) and time evolved using the Hamiltonian H(0.1, 0.4, 0.2).The oscillations are undamped at late times in this approximation up to the light-cone time set by the system size (here L = 2000).

FIG. 8 .
FIG.8.Estimation of the bound state mass using ED compared to the persistent frequency extracted from the meanfield evolution of the initial state ρ(0, 0.4, 0.3, ∞) by the Hamiltonian H(∆, 0.4, 0.3).The ED exhibits large finite size effects when ∆ → 0 so we plot the result of extrapolating to L = ∞ by fitting the gap to a power series in 1/L 2 up to L −4 .

FIG. 9 .
FIG.9.Staggered magnetization for a ground state quench from ρ(0, 0.4, 0.3, ∞) and time evolving with H(0.2, 0.4, 0.3).corresponds to a very small quench with quasiparticle density ∼ 1.6 × 10 −3 .Accordingly very large system sizes and late times would be required to observe decay of oscillations.Inset: half chain entanglement entropy, which grows linearly for all times shown.TEBD calculation done with maximum bond dimension χ = 800.
FIG. 10.Staggered magnetization for a quench from the ground state of H(0, 0, 0.23) and time evolved with H(0.2, 0.4, 0.3) using mean-field theory, the second Born approximation and a TEBD calculation with maximum bond dimension χ = 800 and L = 200.

FIG. 11 .
FIG.11.Main: Staggered magnetization for a quench from the thermal state ρ(0, 0.4, 0.3, 4.0) of the non-interacting system and time evolved with H(0.2, 0.4, 0.3).L = 448 used for the second Born approximation.The initial state is chosen such that the energy density is approximately the same as in Fig.10.Inset: successive peak to peak amplitudes of the oscillations in the second Born approximation, with an exponential fit (dashed black line).The grey scatter points are excluded from the fit.
FIG.12.Staggered magnetization for a quench starting in the thermal state ρ(0, 0.4, 0.10, 10.0) and time evolved with H(0.2, 0.4, 0) on a ring with L = 400.Other than the finite pre-quench temperature Jβi = 10.0 this is the same as Fig.5.