Persistent many-body quantum echoes

We consider quantum many-body systems evolving under a time-independent Hamiltonian $H$ from a nonequilibrium initial state at time $t=0$ towards a close-to-equilibrium state at time $t=\tau$. Subsequently, this state is slightly perturbed and finally propagated for another time period $\tau$ under the inverted Hamiltonian $-H$. The entire procedure may also be viewed as an imperfect time inversion or"echo dynamics". We unravel a remarkable persistence of such dynamics with respect to the observable deviations of the time-dependent expectation values from the equilibrium expectation value: For most perturbations, the deviations in the final state are essentially independent of the inversion time point $\tau$. Our quantitative analytical predictions compare very well with exact numerical results.


I. INTRODUCTION
A trademark of classical chaos [1][2][3] is the sensitivity against small perturbations: Tiny changes in the initial conditions grow exponentially in time. While this effect is readily observable in low-dimensional systems, things are less obvious in the case of many-body dynamics. Indeed, already the mere fact that macroscopic experiments are reproducible suggests that the practically unavoidable small differences in the microscopic initial conditions usually do not result in rapidly growing differences of the actually measurable quantities.
A first key role in this context is played by the following purely geometric effect: Given a system state (classical phase space point) φ and an observable (phase space function) A(φ), small variations δφ are best detectable if they point into the direction n(φ) of the (normalized) gradient of A(φ). In contrast, variations perpendicular to n(φ) are (in leading order) undetectable by means of this observable. Since the subspace orthogonal to n(φ) is of very high dimensionality, the vast majority of all possible variations δφ exhibit very small components parallel to n(φ), and hence the observable changes A(φ + δφ) − A(φ) become unresolvably small.
Focusing on isolated systems with a large (but finite) number of degrees of freedom, yet another key role is played by the fact that they (generically) exhibit equilibration [4][5][6][7], i.e., most initial states φ(0) evolve in time so that the observable A(φ(t)) initially shows a transient relaxation behavior, and later remains very close to a constant value (apart from exceedingly rare recurrences). Moreover, different initial values A(φ(0)) will result in practically identical long-time averages (equilibrium values) of A(φ(t)) under quite weak requirements, namely similar values of the energy and of possibly existing further conserved quantities. In other words, initial differences are largely "forgotten" at later times.
In short, it would be nice to directly "see" (by means of an observable) the exponential growth of small initial perturbations in isolated many-body systems, but this seems not easy to achieve in practice.
A first step consists in evolving φ(t) up to some time point τ , comparable to or larger than the equilibration time of A(φ(t)), and considering the so-obtained state φ(τ ) as a "new" initial condition, which is then evolved under the time-inverted dynamics (so-called echo dynamics [8]). After another time span τ one thus recovers the original initial state φ(0). However, intuitively it now seems reasonable to expect that this inverted evolution will exhibit a more easily observable instability in the sense that small perturbations of φ(τ ) may typically lead to significant deviations of A(φ(2τ )) from the unperturbed "final" value A(φ(0)) of the echo dynamics. Furthermore, one may expect that this instability is closely related to the sensitivity of the original (forwardin-time) solution φ(t) against small initial perturbations. Though all those expectations appear intuitively quite plausible, providing a more tangible justification seems to be a surprisingly difficult task. Fortunately, such a more detailed analytical and numerical confirmation has been provided, for instance, in Refs. [9,10].
A second step in making this exponential instability directly observable is based on the fact that the time inversion is equivalent to a sign change of the Hamiltonian, which in turn can be experimentally realized -at least for spin systems -by inverting the relative orientation of the externally applied (and possibly inhomogeneous) magnetic field in cases with non-interacting spins (so-called Hahn echo [8]), and by means of more sophisticated "tricks" in cases with interacting spins [11][12][13][14][15][16][17][18].
The main topic of our present work is the corresponding echo dynamics within a quantum mechanical instead of a classical approach. Besides its conceptual interest, this is obviously also of relevance with respect to the experimental realizations in terms of spin systems, cold atoms, and quantum simulators [8,[12][13][14][15][16][17][18][19][20][21][22].
The above arguments of why small initial perturbations usually also remain unobservable at all later times (or, equivalently, why macroscopic experiments are reproducible) can be readily adapted to the quantum case. But what about the echo dynamics?
In related previous works, the main focus was on the similarities and differences between classical and quantum spin models [23,24] and on questions of whether and in which sense the aforementioned signatures of classical chaos can indeed be recovered also in quantum systems in the thermodynamic limit [25][26][27]. Notably, the thermodynamic limit was thus formally carried out before letting the echo imperfections go to zero, and it was found that the two limits actually do not commute [27].
In our present approach, the general setup is similar as in those previous works [23][24][25][26][27], whereas the main focus is laid on a somewhat different issue, namely on quantitative analytical approximations for generic, large but finite systems in the presence of sufficiently small (but finite) imperfections. The pertinent key results of our paper consist in the prediction (21) about the detailed behavior of the observable echo dynamics, and in the demonstration that it agrees very well with numerical results. In particular, the analytical result (21) predicts that the observable perturbation of the quantum echo due to tiny imperfections during the reversal procedure does not grow at all upon increasing the inversion time point τ , and this behavior is indeed recovered numerically up to possible transient effects for very short values of τ .
The paper is organized as follows: In Sec. II, we introduce the general setting of the considered echo dynamics and types of imperfections. The main result is then derived and verified numerically in Sec. III, before we combine it with results from [28] for a related scenario in Sec. IV. Finally, we summarize and discuss our findings in Sec. V.

II. ECHOES AND IMPERFECTIONS
In this section, we introduce the general framework of our echo approach and the objectives of the later sections.
To begin with, we consider a finite, isolated many-body quantum system with time-independent Hamiltonian H = n E n |n n| . (1) Before turning to questions of time reversal and its instability against small perturbations, we collect the essential properties of the "clean" (time-independent) dynamics.
If the system is prepared in a (pure or mixed) state with density operator ρ(0) at time t = 0, it then evolves according to the Schrödinger or von Neumann equation, so that the state at any later time t > 0 is ρ(t) := e −iHt ρ(0)e iHt ( = 1). Given some observable (selfadjoint operator) A, its expectation value in an arbitrary state ρ follows as A ρ := Tr{ρA}. For the time evolution of the expectation value in the system state ρ(t) we thus obtain Despite the quasiperiodic nature of this expectation value, a many-body system is known to exhibit equilibration under rather weak conditions [29][30][31], i.e., after the initial transients have died out, it spends most of its time close to the time-independent value A ρeq , where ρ eq is (essentially) the long-time average of ρ(t) [32]. Moreover, generic (nonintegrable) systems also thermalize [4][5][6], meaning that A ρeq can be identified with A ρmc , where ρ mc is the pertinent microcanonical ensemble for the isolated system at hand.
We emphasize that our main result will actually turn out to apply irrespective of whether or not the system exhibits equilibration, let alone thermalization. However, if it does thermalize, then to study the effect of imperfect echo dynamics, the system must start out sufficiently far from equilibrium, i.e., there must exist an initial time interval during which notably deviates from zero. Assuming that the system is in a nonequilibrium state at time t = 0, we then probe how special this situation is by exploring how difficult it is to return to this state by an effective, but possibly imperfect reversal of the dynamics after the system has relaxed for some time. In other words, imagine that the many-body system with Hamiltonian H is prepared in a well-controlled far-from-equilibrium state ρ(0) = ρ T , called target state, such that A(0) = 0. Then, the system evolves according to the Hamiltonian H until time τ , resulting in a return state ρ R := ρ(τ ). From Eq. (2), we understand that evolving the system backward in time with the Hamiltonian H is equivalent to evolving it forward in time with the "backward" Hamiltonian −H. Switching to this backward Hamiltonian at t = τ , the resulting evolution from ρ R to ρ T during the interval [τ, 2τ ] is just the reflection of the initial process during [0, τ ] and constitutes our reference dynamics, i.e., the perfect time reversal with for any t ∈ [0, τ ]. Note that even in cases which are not expected to exhibit thermalization, the difference in (3) will usually be quite appreciable at early times t, but may possibly no longer approach zero for (most) sufficiently late times. In any case, the quantity in (3) will later turn out to play a crucial role within the analytical approach adopted in our present paper (see Sec. III B below).
Finally, two different types of imperfections are often considered [23][24][25][26][27] leading to deviations from the above reference scenario. In the first case [sketched in Eq. (6) below], one slightly perturbs the return state ρ R by evolving it for a short time δ with a "scrambling Hamiltonian" V [33] to obtain a different state ρ ′ R , modeling an imprecise preparation of initial conditions for the timereversed setup. Thereafter, the system is evolved with the "true" backward Hamiltonian −H. (Strictly speaking, the backward evolution thus takes place in the time interval [τ + δ, 2τ + δ], but for simplicity we will sometimes approximately denote it as [τ, 2τ ].) In the second case [cf. Eq. (31)], one starts from the original return state ρ R , but now includes a small perturbation in the backward Hamiltonian H ′ := −H + ǫW , modeling an imperfect implementation of the time-reversed system.
Henceforth, those two types of perturbations will be referred to as imperfect preparation and imperfect reversal, respectively.
The main focus of our present paper will be on perturbations of the first type (imperfect preparation). But for the sake of completeness and comparability, we will also briefly summarize our previous findings from Ref. [28] for perturbations of the second type. Moreover, we will provide the pertinent extension when both types of perturbations are simultaneously present.
In either case, one intuitively expects that a generic imperfection of the time-reversed dynamics typically brings the system closer to equilibrium because it impairs the fine-tuned correlations between state and observable that are vital to obtain a reasonably far-from-equilibrium expectation value in (2). Hence, the imperfect backward dynamics will usually remain closer to the equilibrium state than the forward one, i.e., The chaoticity and irreversibility of the many-body dynamics may thus be gauged in terms of the sensitivity of the deviations between the perfect and scrambled echo dynamics with respect to the two above introduced scrambling parameters δ and ǫ: The more A(τ + t) is attenuated with increasing δ or ǫ compared to A(τ − t), the harder it is to design a reversible process and the more extraordinary or special are the initially probed nonequilibrium states. Accordingly, the ratio A(τ + t)/A(τ − t) for times t ∈ [0, τ ] is at the heart of our present study, most importantly in the region around the revival or echo peak of A(τ + t) at t ≈ τ , where observable deviations from equilibrium will stand out most prominently, cf. (3).
Our main results will be analytical predictions for the ratio A(τ +t)/A(τ −t), based on a few basic properties of the scrambling operators V and W . In the derivation, we assume that these operators are -essentially by definition -uncontrolled (random) and, in particular, independent of the actual system Hamiltonian H, in a sense that will become more clear below.
Finally, we will also verify our predictions by comparison with numerically exact results for a spin-1 2 XXX model.

A. Setup
As outlined above, our echo protocol modeling an imperfect initial state for the backward evolution consists of three phases: After preparing the system in the target state ρ T , it is first evolved with the forward Hamiltonian H for a time τ to reach the return state ρ R . From there, second, the system is subject to a perturbing Hamiltonian V for the short scrambling time δ, yielding the state ρ ′ R . Third, we then evolve it backwards with the Hamiltonian −H for a time τ , ending up in a state ρ ′ T . The entire process can thus be summarized as By analogy with classical systems and their pertinent indicators of chaos [1,2], one may roughly regard this approach as a probe for sensitive dependence on initial conditions since it quantifies how a small difference between the states ρ R and ρ ′ R translates into a difference in the states ρ T and ρ ′ T after time τ . In this sense, the scrambling phase taking ρ R to ρ ′ R need not necessarily be viewed as an actual time evolution. Instead, more abstractly, one may assume a Lie group perspective and regard e iV δ as a "rotation" in Hilbert space by an angle δ in the direction V . When comparing different systems, the operator V should then be properly normalized. Nevertheless, for notational convenience, we will include the scrambling as part of the time evolution in the following.
We denote the time-dependent state of the system in the forward, scrambling, and backward phases by respectively. We also write ρ(t) to refer to the state during the entire process, i.e., ρ(t) An implicit assumption in all what follows is that the considered many-body system under study is finite and exhibits a well-defined macroscopic energy E. In particular, it is understood that the perturbation taking the return state ρ R to ρ ′ R is so small that it does not change this energy, i.e., there should be no macroscopically measurable energy intake during the scrambling phase. Consequently, the state ρ(t) at any time can only significantly populate energy levels within a macroscopically small window [E − ∆E, E]. Therefore, we can and will implicitly restrict summations over energy levels to this window in the following. Given the extremely small level spacing in many-body systems, the number of levels N within such an energy shell is still exponentially large in the system's degrees of freedom f [36], Another temporary assumption we make is that the return state ρ R populates the levels of the scrambling Hamiltonian V approximately uniformly within the relevant energy interval [E − ∆E, E]. This assumption seems reasonable in the absence of any additional knowledge about the experimental imperfections modeled by V . Should such information be available, it is still possible to generalize our present approach along similar lines as in Ref. [37] in order to account for significant nonuniformities in the level populations. We will comment on the necessary modifications at the end of Sec. III D.
Aiming at a prediction for the echo dynamics, we now turn to the third phase in the evolution (6), where the state of the system is [cf. Eqs. (7)] Recalling that the eigenvalues and eigenstates of the Hamiltonian H are E n and |n , cf. Eq. (1), and denoting those of V by E V ν and |ν V , respectively, the matrix elements m|ρ b (t)|n can be written as Here U µk := V µ|k are the transformation matrices between the eigenbases of H and V . For a particular observable A, the time-dependent expectation value is then given by We recall that the Hamiltonian H corresponds to a given many-body quantum system, whereas the perturbation V is supposed to describe an uncontrolled experimental uncertainty in the time-reversal procedure with the result of an imperfectly prepared return state ρ ′ R . In this spirit, we will thus model the scrambling Hamiltonian V by a random Hermitian matrix, conforming with our partial knowledge and partial ignorance about the actual experimental imperfections.
Regarding the statistical properties of V , we adopt the rather general framework from Ref. [38], where eigenvalues and eigenvectors are considered as statistically independent quantities. While the former (i.e., the E V ν 's in (10)) are considered as "given", i.e., as largely arbitrary but fixed (non-random), the latter are generated according to a uniform (unbiased) distribution, i.e., the U µk 's in (10) are sampled via the Haar measure of the appropriate symmetry group, e.g. the circular unitary ensemble (CUE) for U(N ) or the circular orthogonal ensemble (COE) for symmetric unitary matrices. For instance, the well-known standard case of a Gaussian unitary Vensemble (GUE) is recovered by combining the CUE with a semicircular density of states (DOS) for the eigenvalues [39]. In general, i.e., unless there is specific additional knowledge about the structure of typical experimental uncertainty in the time-reversal procedure, it may however be reasonable to also leave room for some different (non-semicircular) DOS, for instance one which exhibits some similarity with the DOS of the system Hamiltonian H.
This approach provides us with a useful prediction for the behavior of an actual system in two steps: First, we calculate how such a perturbation affects the observed dynamics on average. Second, we show that the squared magnitude of fluctuations around the average is inversely proportional to the number of states N in the pertinent energy shell. For large N , this concentration of measure property ensures that a single realization will typically not deviate from the ensemble average in practice. Bearing in mind Eq. (8), we may then conclude that such a modeling should indeed account for most uncontrolled uncertainties in setting up the return state.

B. Results
The average effect of the scrambling with e −iV δ in (9) is obtained by averaging over the corresponding, above specified U ensemble in Eq. (10). Recalling that we keep the spectrum of V fixed, the only random quantities are the four factors of U µk , i.e., the eigenvectors of V in the eigenbasis of H. Denoting this average over the Haardistributed U µk by E[ · · · ], we can consult Ref. [40] to find the general form where the so-called symmetry factors v 1,1 and v 2 are given by v 1,1 = and v 2 ≃ − 1 N 3 in both cases. Exploiting all this in carrying out the average of Eq. (10), one finally obtains In the second line, we retrieve the forward state from (7a). In the third line, we can identify the matrix elements of the microcanonical density operator ρ mc , The remaining terms can be conveniently rewritten by means of the DOS of the scrambling Hamiltonian V in the relevant energy shell, and its Fourier transform, Altogether, Eq. (13) then takes the form (18) Observing that the 1/N term is negligible due to N ≫ 1, and exploiting the definition (3) with ρ(t) as introduced below (7), we then find This is the first key result of this section, relating the relative strength of the ensemble-averaged echo signal under imperfect preparation of the return state to the density of states of the scrambling operator.
To assess deviations between the observed dynamics for one particular realization of V and the above derived average behavior, we consider the variance of A(τ +δ+t). In view of (10) and (11), we thus need averages over products of eight matrix elements U µk . Employing Ref. [40] once again, the result is structurally similar to Eq. (12), i.e., it consists of Kronecker-δ contractions multiplied by symmetry factors such as v 1,1 , etc. Referring to Appendix A for the details, one finally obtains where ∆ A indicates the measurement range of the observable A, i.e., the difference between the largest and smallest possible measurement outcomes (eigenvalues). Similarly as below Eq. (12), the differences between the CUE and the COE only appear in the higher-order corrections (last term in (20)), and similarly as below (18), these corrections will be henceforth neglected.
The result (20) establishes the following concentration of measure property of the considered ensembles of V operators: Roughly speaking, the fluctuations for one particular V typically decrease like 1/ √ N and thus exponentially fast in the system's degrees of freedom, see Eq. (8). More precisely, Chebyshev's inequality implies that the probability for A(τ +δ+t) to differ by more than ∆ A /N 1/3 from the average E[A(τ + δ + t)] at a certain instant in time t is less than 11/N 1/3 . This is our second key result, providing predictive power to our first key result from Eq. (19). Namely, since observable deviations from the average are extremely scarce for systems with a reasonably large number of degrees of freedom, we can conclude that amounts to a very good approximation for nearly all imperfections in the return state preparation which can be modeled by (6). Eq. (21) is our main result regarding the generic echo dynamics under imperfect initial conditions for the backward evolution. As announced below (4), we see that the quantity from (3) indeed plays a crucial role, independently of whether or not the considered system is expected to exhibit thermalization.
The particularly interesting revival or echo peak signal is recovered by choosing t = τ in Eq. (21). Moreover, it can be recast by means of (21) and (3) into the equivalent form of an approximately uniformly spread-out return state ρ R in the eigenbasis of V [41], this result also confirms the finite-size analysis of Ref. [27] within the present typicality approach.
According to (21), we thus expect that the "intensity" of the revival dynamics (left-hand side) is determined by the Fourier transformD V (δ) of the perturbation's density of states evaluated at the scrambling time δ. The longer this scrambling time, the farther the effective return state is carried away from the perfect state ρ R , and the weaker the echo peak signal will typically be. Notably, though, the echo peak signal is predicted to be independent of the propagation time τ , see also (22). In this respect, the quantum echo dynamics is thus of a surprisingly persistent character. This should be contrasted with the classical case where -as one might intuitively expect from chaos theory -imperfectly prepared echo signals decay with τ [9, 10].

C. Example
As an illustrational example, we consider a spin-1 2 XXX chain of length L with Hamiltonian where σ i = (σ x i , σ y i , σ z i ) is a vector collecting the Pauli matrices acting on site i. The initial (target) state ρ T = |ψ ψ| is a Néel state, |ψ = |↓↑↓↑ · · · , and we observe the staggered magnetization  (23) with L = 14 under the "imperfect preparation" protocol (6), starting from a Néel target state ρT, for various scrambling times δ. The time reversal is initiated after time τ = 2 in (a) and τ = 10 in (b). The perturbation Hamiltonian V is of the "spin-glass" form (25) for one realization of the normally distributed random couplings J αβ ij . (For other realizations, we found practically the same results.) Solid lines correspond to the numerical results using exact diagonalization, color-coded as black for t ∈ [0, τ ] (forward evolution), and brown-to-purple as indicated in the legend for t ∈ [τ, 2τ + δ] with the peak height decreasing as δ increases. Dashed lines show the analytical prediction for t ∈ [τ + δ, 2τ + δ] (backward or echo dynamics) according to (21), (26), and (27), complemented by (3)  For the scrambling operator V , we choose with independent, normally distributed couplings J αβ ij of vanishing mean and unit variance. Hence the direction in Hilbert space in which the return state ρ R is perturbed is somewhat random and reflects the erratic character of the imperfections, but these still respect the overall setting and spin structure of the model under study. Note that in the more abstract picture of V as a generator of Hilbert space rotations [see below (6)], V should be normalized by dividing by the number of terms on the right-hand side in order to compare perturbation strengths across different system sizes.
For L ≫ 1, the DOS of V from (25) is known to be a Gaussian with mean zero and variance [42] With (17) and (21), we thus expect that the echo signal is characterized byD In Fig. 1 we show the resulting dynamics under the protocol (6) for two durations τ and several scrambling times δ, and compare it to our analytical result (21). We emphasize that, by exploiting (26) and (27), this theoretical prediction does not involve any fit parameters. Its agreement with the numerics is practically flawless, illustrating the decay of the echo peak as the intensity δ of the perturbation is increased. At the same time, the persistence of the echo signal becomes apparent since the peak heights for fixed δ are identical for both values of τ . As we will discuss at the end of Sec. III D (see also Appendix B), the numerically observed echo signal does show some τ dependence for very short τ 1 (below the relaxation time of the reference system), so that deviations from our prediction occur in this regime.

D. Perspectives
It is instructive to reconsider the first of the three steps in (6) by switching from the so-far adopted Schrödinger picture to the equivalent Heisenberg picture and, moreover, by considering time as evolving backward. In other words, we observe that ρ(τ − t) = U t ρ(τ )U † t with backward-in-time propagator U t := e iHt , and hence for all t ∈ [0, τ ]. Likewise, upon observing that U t from above amounts to the forward-in-time propagator for the inverted Hamiltonian −H, the last step in (6) can be rewritten as A ρ(τ +δ+t) = Tr{ρ(τ + δ)A(t)} = Tr{ρ ′ R A(t)} (30) for all t ∈ [0, τ ]. Finally, one readily verifies that the measurement range ∆ A(t) [see below (20)] and the microcanonical expectation value A(t) ρmc are independent of t. Altogether, one thus can conclude that if the relation (21) is valid for arbitrary observables at one single t value, say t = 0, then it is valid for any given observable at arbitrary t ∈ [0, τ ].
On the one hand, this conclusion is helpful to develop some intuitive understanding of why the ratio on the lefthand side of (21) is predicted not to depend on t (see right-hand side). On the other hand, in order to verify (21) for arbitrary A and t, it is sufficient to solely focus on t = 0. Put differently, the first and the last of the three steps in (6) are actually irrelevant for our present purpose (both time evolutions can be omitted by considering the modified observables {A(t)} τ t=0 instead of the original A). As a consequence, the result (21) is expected to remain valid even for explicitly time-dependent Hamiltonians H(t) in (6).
Another immediate implication is that the system does not need to (approximately) approach a steady state at the end of the first step in (6) (even if τ is very large), i.e., neither equilibration nor thermalization are required. Fig. 2 depicts a particularly interesting example of this kind, which arises when the initial state ρ(0) only populates two levels of the Hamiltonian H, hence the expectation value in (2) exhibits perpetual oscillations in time. According to (21), the echo dynamics thus also exhibits everlasting oscillations in spite of the scrambling effects via the second step in (6). The reason is that also the level populations, corresponding to observables of the form A = |n n|, are only moderately changed according to (21) for not too large scrambling time δ.
Similarly as in Fig. 1, our analytical theory again reproduces the numerically exact results in Fig. 2 remarkably well without any free fit parameter.
Finally, we come back to the assumption from below Eq. (8), stating that the return state ρ R should occupy the energy levels of V roughly uniformly. If this is not the case, it is straightforward to adapt the formalism of Ref. [37] to obtain a refined prediction for the echo dynamics that includes information about occupation imbalances. The essential modification affects the function D V (δ) such that the integrand in (17) receives an additional energy-dependent weight according to the populations of the levels around E. While the definition of the functionD V (δ) appearing on the right-hand side in (21) may thus be modified, the key point is that it still remains independent of t. Moreover, it also remains (approximately) independent of τ , provided τ exceeds the characteristic relaxation time scale of the forward dynamics (governed by H in (6)). For smaller τ values, things may in general become considerably more complicated. Indeed, for the same example as in Fig. 1, we found that the differences between numerics and theory quite notably increase in the regime τ 1, see also Appendix B. Analogous numerical findings for small-tomoderate τ values have also been observed and discussed in Refs. [23,24,27]. One might be tempted to conjecture that the few-body character of the scrambling operator V in those examples invalidates a faithful modeling in terms of the random matrix ensemble introduced below Eq. (11) in this regime. However, we observed no such deviations in the setting of Fig. 2, involving the same H and V as in Fig. 1, but different A and ρ T . Hence the few-body form of V alone does not explain the deviations. Since such small-τ deviations are also absent in a "direct" implementation of the perturbation ensemble from Sec. III A, which we demonstrate by means of example in Appendix B, they hint at (mild) correlations between the observable A and the initial state ρ T on the one hand and the Hamiltonian H and the scrambling operator V on the other hand. Crucially, by further increasing τ in the examples from Fig. 1 or Appendix B, we found, as predicted above, that the agreement between numerics and theory remains unchanged. In particular, the persistence of the echo peak signal is maintained.

IV. IMPERFECT REVERSAL AND COMBINED IMPERFECTIONS
As pointed out in Sec. II, a second basic type of experimental inaccuracy in the considered echo dynamics consists in an imperfect implementation of the time-reversed Hamiltonian [25]. Such an "imperfect reversal" scenario is particularly relevant with regard to experimental applications because (a) it is usually only feasible to invert the dominant part of the respective Hamiltonian, and (b) the accuracy to carry out the elaborate protocols necessary to do so is subject to experimental constraints [11][12][13][14][15][16][17][18]. In fact, inverting only parts of the Hamiltonian can even be desirable as the degree of echo-signal attenuations can reveal information about the "imperfections" and thus the structure of the probed material, the paradigmatic example being magnetic resonance imaging (MRI).
To separate the influences of such an imperfect reversal procedure from the imperfect preparation scenario discussed in Sec. III, we assume for now that the return state is set up perfectly for the time-reversed evolution. However, instead of the clean backward Hamiltonian −H, we consider a perturbed H ′ = −H + ǫW , leading to a total echo process of the form Formally, this protocol may be regarded as a generalization of the Loschmidt echo [43][44][45], since the latter is recovered (as a function of τ ) for the special choice A = ρ T (whose experimental realization in a many-body system may be difficult in general [25][26][27] but not impossible in specific examples [45]). Combinations of the two scenarios from (6) and (31) will be addressed at the end of this section. Similarly as in the previous section, we have in mind a given system with fixed Hamiltonian H whose dynamics is subject to some (partially) uncontrolled or unknown experimental imperfection, described by W in (31). Accordingly, this lack of control is again modeled by choosing an appropriate ensemble of random perturbations W , whose characteristics should nevertheless still emulate the key features of typical imprecisions in many-body systems. Notably, we therefore allow for a sparse distribution of the matrix elements m|W |n and a dependence of their amplitude on the energy difference E m − E n of the involved levels [46][47][48][49][50]. In other words, the random matrices m|W |n may (but need not) exhibit a so-called sparse and/or banded structure [51].
The just-described setup was laid out and investigated in detail in Ref. [28], so that we only summarize the key findings here. Moreover, the conceptual approach is also very similar to the one from the previous Sec. III: The first step consists in averaging the echo dynamics from the protocol (31) over the considered ensemble of perturbations W . Denoting this average by E[ · · · ], too, the so-obtained average suppression of the echo signal then takes the form where α = πσ 2 w D 0 is a constant involving the varianceσ 2 w of the matrix elements m|W |n for small E m − E n . Furthermore, D 0 is the density of states of H in the vicinity of the system energy E (see also Sec. III A).
Similarly as for the result from Sec. III B, this average (32) gains predictive power in systems with many degrees of freedom because the corresponding variance is found to satisfy the bound [28] where c is a constant of the order of 10 3 or smaller, and N w := 2αǫ 2 D 0 quantifies how many unperturbed energy levels the perturbation W couples on average. Crucially, this latter number usually scales exponentially with the system's degrees of freedom [52], so that (33) indeed becomes exceedingly small for typical many-body systems. Altogether, we thus can conclude that in a system with sufficiently many degrees of freedom, nearly all imperfections W in the time-reversed Hamiltonian which can be faithfully modeled by any of the above specified random matrix ensembles lead to an attenuation of the echo signal of the form [28] A(τ + t) Contrary to the imperfect initial conditions [see Eq. (21)], the intensity of the echo signal at t ≈ τ now decreases as the reversal time τ grows. Essentially, this may be understood as a consequence of the fact that the system is continuously exposed to the perturbation W during the entire backward evolution of duration τ . Finally, it is conceivable and in fact practically inevitable that both types of imperfections discussed here actually contribute to the deterioration of an echo signal. Analogously to (6) and (31), such a scenario may be symbolically indicated as Provided that those two types of imperfections have "independent" origins (corresponding to statistically uncorrelated randomization effects of V and W ), the results from Secs. III and IV can be readily combined by multiplying the right-hand sides of Eqs. (21) and (34) to obtain a prediction for the total effect of the form While both types of imperfections thus generically cause a mitigation of the echo signal, only the persistent exposure to perturbations of the imperfect reversal type from this section entails a stronger decay with increasing reversal time τ .

V. SUMMARY AND CONCLUSIONS
At the focus of our present work is the echo dynamics of large but finite many-body quantum systems when weakly perturbed according to the so-called imperfect preparation scheme from (6): An initial state ρ T evolves under the action of a Hamiltonian H during a time period τ , is then slightly perturbed by a scrambling operator V for a short time δ, and is finally propagated by the inverted Hamiltonian for another time period τ . The key question is: How does such an imperfect echo dynamics manifest itself in the time-dependent expectation values of an observable A, and, in particular, how do the deviations from a perfectly unperturbed echo dynamics depend on the propagation time τ , the perturbation time δ, and on the details of the perturbation operator V ? As our main result we predicted that most perturbations V , which exhibit some rather weak general properties, will satisfy in very good approximation the relation (21), connecting the expectation values of A from (3) with the Fourier-transformed density of states of V from (16) and (17). Figs. 1 and 2 illustrate the very good agreement of this prediction with numerically exact results without any free fit parameter.
The most remarkable feature of Eq. (21) is that the right-hand side is independent of t and τ . In particular, for t = τ the deviation between the initial and the final expectation values in (22) is independent of the propagation time τ in (6), see also Fig. 1 for an explicit example. In contrast to the classical case, where echo signals usually attenuate with increasing τ [9,10], the quantum echo dynamics is thus found to be of a remarkably persistent character. Some heuristic arguments in support of this property have been provided in Sec. III D as well as in Ref. [27], but for a fully satisfying explanation, one apparently cannot circumvent the non-trivial calculational details from Sec. III B and Appendix A.
For the sake of completeness, we also summarized in Sec. IV our previously obtained results from [28] for the so-called imperfect reversal scenario (31), and we provided the pertinent generalization for cases which simultaneously exhibit an imperfect preparation and an imperfect reversal, cf. (35).
As discussed at the end of Sec. III D, our present work is complementary to the earlier studies from Refs. [23,24], where the main focus was on small-to-moderate propagation times τ in (6). Moreover, our present approach is also complementary to the one adopted in Refs. [25][26][27], where the main focus was on the echo dynamics in the thermodynamic limit, probing to what extent classical chaos indicators reemerge there [25,26], and how this limit is related to the semiclassical limit → 0 [27]. In passing, we note that also Refs. [23,24] had mainly in mind this limiting case in the theoretical discussions, while the numerical examples were actually closer in spirit to our present setup. Furthermore, it may be recalled that the concept of chaos in classical systems usually does not involve the thermodynamic limit, but rather is mostly considered in systems of finite and possibly even low dimensionality.
While our findings here suggest that one signature of classical chaos in finite-dimensional systems, namely sensitive dependence on initial conditions, may not have a corresponding quantum analog, we finally mention that there are various well-established or recently proposed alternative indicators of chaos (or nonintegrability) in fewand many-body quantum systems, such as level statistics [49], eigenstate thermalization [4], or out-of-timeorder correlators [34,35], to name but a few. Since we are not aware of any tangible mathematical or physical connection with our present approach, we desist from a more detailed compilation regarding the agreement or disagreement with all those alternative indicators for the "chaoticity" of any given specific model.
In this appendix, we derive the bound (20) on the fluctuations of A(τ + δ + t) under imperfect preparation of the return state as studied in Sec. III. We first observe that (A1) Combining Eqs. (10) and (11), we can write (A2) for the time-dependent expectation values during the backward evolution. As derived in the main text, cf. Eq. (18), the average over all perturbations V of this expression is given by from which we can read off E[ A ρ b (t) ] 2 straightforwardly. Squaring (A2) leads to Computing the average thus amounts to averaging over the eight factors of U µk . Following Ref. [40], this average is given by a sum over all possible combinations of pairing up the first and second indices of U and U * factors, i.e., Here Sym(4) denotes the symmetric group of order 4, i.e., the set of all permutations of (1,2,3,4), so that the sum in (A5) consists of (4!) 2 = 576 terms. The symmetry factors v P,P ′ depend only on the cyclic structure of P −1 P ′ . There are five different combinations of cycle lengths in Sym (4), and the leading order as N ≫ 1 of the corresponding factors for CUE or COE matrices is [40] v 1,1,1, As announced in the main text, we restrict to this leading order analysis here. Substituting (A5) into (A4), the sums over µ j and k j , l j factorize, so that we can analyze them individually. The factors involving µ j have the general form with P ∈ Sym(4). The order of these terms in N is again determined by the cyclic structure of the relevant permutation P . The different contributions are summarized in Tab. I. The general form of the factors with k j , l j for a permutation P ′ ∈ Sym(4) is These reduce to the 10 different combinations of state and observable summarized in Tab. II, of which there are two of order N 2 , four of order N , and five of order 1.
Combining the factors v P,P ′ from (A6), F P from (A7) and Tab. I, and G P ′ from (A8) and Tab. II, there are nine terms of order 1 in the average of (A4). These cancel exactly against the nine terms of order 1 obtained by squaring (A3), so that the variance E[( A ρ b (t) ) 2 ] − E[ A ρ b (t) ] 2 vanishes to order 1.
At order N −1 , there are 42 terms in the average of (A4) and six terms in the square of (A3), giving Using the triangle inequality, the inequality |D V (δ)| ≤ 1, and bounding the absolute value of all the different traces ( 1 2 3   over factors of A and ρ by A 2 (operator norm), we then arrive at the estimate As in the main text, we tacitly focus on cases where A models an observable with a finite measurement range ∆ A (difference between largest and smallest eigenvalues of A, see below (20)) and thus with finite A . Without loss of generality, we can and will assume (possibly after adding an irrelevant constant to A) that the largest and smallest eigenvalues of A are equal in modulus and of opposite sign, implying that A = ∆ A /2. Together with (A10), this finally yields Eq. (20) from the main text.  Fig. 1, but now for three rather small reversal times τ = 0.2, 0.5, 1. In other words, the ratio A(2τ + δ)/A(0) of echo and initial peak heights is depicted as a function of the scrambling time δ for the staggered magnetization (24) of a spin-1 2 XXX chain (23) with L = 14, normally distributed J αβ ij in (25), "imperfect preparation" protocol (6), and Néel target state ρT. The colored dots are the numerical solutions for increasing τ from top to bottom, the solid red curve is the analytical prediction from (21), (26), and (27). . Results are shown for various scrambling parameters δ (color-coded as indicated, with the peak height decreasing as δ increases) and reversal times τ = 0.5, 2, 5. Solid lines correspond to the numerics, dashed lines are the theoretical predictions for the backward evolution according to (21) with (B3) and A ρmc = 0. The echo peak height is entirely independent of τ .

Appendix B: Further numerical examples
In this appendix, we elaborate in more detail on the deviations between our theoretical prediction (21) and the spin-model numerics from Sec. III C that occur for τ 1 when starting from a Néel state and observing the staggered magnetization. We emphasize that these deviations do not affect our main conclusion, namely the persistence of many-body quantum echoes at arbitrarily large times: For τ 1 and sufficiently small δ, the echo peaks remain finite and become independent of τ . (For the examples of Figs. 1 and 2, we checked reversal times up to τ = 20 and τ = 36, respectively. ) We also emphasize that there are no such deviations for small τ in the setup from Fig. 2, where the reference Hamiltonian H and scrambling operator V are identical to the ones from Fig. 1 [Eqs. (23) and (25)], but the initial state and observable are superpositions of two energy eigenstates instead.
Coming back to the setup from Fig. 1, we depict in Fig. 3 the echo peak height in the spin model from Sec. III C for reversal times τ = 0.2, 0.5, 1 as a function of the scrambling parameter δ with the solid red line indicating the theoretical prediction from Eq. (21). This figure should be contrasted with the insets of Fig. 1. For fixed δ, the actual (numerical) echo signal is attenuated less than our theory predicts, with the deviations becoming stronger as τ decreases.
One reason for this may be that the return state ρ R for small τ violates the assumption [see below Eq. (8)] of being uniformly spread across the V eigenvectors in the energy shell (see also the discussion at the end of Sec. III D), i.e., the special structure of the Néel initial state is only gradually stirred up during the forward evolution. A related problem can arise when the perturbation V is correlated with H, A, and/or ρ T , so that a reliable modeling in terms of the unbiased ensemble introduced in Sec. III A is not guaranteed. Such issues are expected to be most compromising when V is applied on a state ρ R that is still out of equilibrium (i.e., for small τ ) since this is the regime dominated by correlations.
Indeed, we observe that there are no small-τ deviations when one avoids all these correlations by choosing a sufficiently randomized system. To this end, we consider a Hamiltonian H = n n D |n n| (B1) with a fixed level spacing 1/D and corresponding level density D = 512, in combination with an observable A which is randomly sampled from the Gaussian orthogonal ensemble (GOE). The initial state ρ T = |ψ ψ| is generated by exploiting dynamical typicality [53] from a Haar-distributed random Hilbert space vector |φ according to |ψ ∝ (1 + γA)|φ (B2) with γ = 1. For the perturbation V , we draw a matrix at random from the Gaussian unitary ensemble (GUE), so that D V (E) is a semicircle of radius 2 and consequentlŷ where J 1 (x) denotes the Bessel function of the first kind of order 1. Once again, we thus know all parameters entering the theoretical prediction (21) exactly and can compare the analytics with numerical simulations without any fit parameters. Turning to Fig. 4 for this comparison, we find that the agreement is virtually flawless, irrespective of the magnitude of either τ or δ. In the absence of any correlations (apart from those between A and ρ T necessary for nonequilibrium initial conditions), the small-τ deviations thus disappear, too.