Stability of Exponentially Damped Oscillations under Perturbations of the Mori-Chain

There is an abundance of evidence that some relaxation dynamics, e.g., exponential decays, are much more common in nature than others. Recently, there have been attempts to trace this dominance back to a certain stability of the prevalent dynamics versus generic Hamiltonian perturbations. In the paper at hand, we tackle this stability issue from yet another angle, namely in the framework of the recursion method. We investigate the behavior of various relaxation dynamics with respect to alterations of the so-called Lanczos coefficients. All considered scenarios are set up in order to comply with the"universal operator growth hypothesis". Our numerical experiments suggest the existence of stability in a larger class of relaxation dynamics consisting of exponentially damped oscillations. Further, we propose a criterion to identify"pathological"perturbations that lead to uncommon dynamics.


I. INTRODUCTION
The apparent emergence of irreversibility from the underlying reversible theory of quantum mechanics is a longstanding puzzle that lacks an entirely satisfying answer to this day [1]. Over the course of the last decades, fundamental concepts like the "eigenstate thermalization hypothesis" [2][3][4] and "quantum typicality" [5][6][7] have crystallized, which give conditions under which isolated quantum systems eventually reach an equilibrium state. However, while these mechanisms ensure eventual equilibration, they make no statement in which manner the equilibrium state is actually reached, i.e., they do not narrow down the eligible routes to equilibrium. In contrast, it is evidently true that some relaxation dynamics, e.g., exponential decays, are much more common in nature than others, e.g., recurrence dynamics. There are recent attempts to explain this prevalence with the idea that some dynamics are stable versus perturbations of the Hamiltonian. These attempts include, for example, investigations based on Hamiltonian perturbations on the level of random matrices [8]. Further, there exist advances that consider an entire ensemble of permissible perturbations [9][10][11]. It is then analytically shown that weak "typical" perturbations lead to an exponential damping of the original dynamics [9]. This renders exponential decays stable since only the decay constant changes, whereas recurrence dynamics are exponentially suppressed.
With the paper at hand, we address the issue of stability of certain [classes of] relaxation dynamics in the framework of the recursion method [12,13]. Central quantities that appear within this framework are the so-called Lanczos coefficients, real numbers that characterize the complexity growth of operators over the course of time. As will be presented in Sec. II, the Lanczos coefficients can * rheveling@uos.de † jiaowang@uos.de ‡ cbartsch@uos.de § jgemmer@uos.de be interpreted as hopping amplitudes in a tight-binding model. In this manner, many physical problems, like calculating correlation functions, can practically be reduced to a one-dimensional [finite or semi-infinite] chain, which we refer to as the "Mori-chain" in the title of this paper.
In the following, we consider perturbations on the level of Lanczos coefficients and examine the resulting effect on various kinds of relaxation dynamics. The particular advantage of our approach within the recursion method framework is that all considered scenarios can be set up to directly comply with the universal operator growth hypothesis [14], which previous approaches have been lacking [8]. Said hypothesis concerns the asymptotic growth of the Lanczos coefficients and it basically states that the coefficients should eventually attain linear growth [with a logarithmic correction in one dimension]. The hypothesis is backed up by analytical as well as numerical evidence [14][15][16]. The paper at hand is structured as follows: Sec. II constitutes a preliminary section on the recursion method and the operator growth hypothesis. Afterwards, in Sec. III, the concrete strategy to study the stability of certain classes of dynamics is explained in detail. In Sec. IV we present and discuss our numerical results. We conclude in Sec. V.

II. PRELIMINARIES: RECURSION METHOD AND OPERATOR GROWTH HYPOTHESIS
In this section, we briefly recall the basics of the recursion method [12,13] as well as the universal operator growth hypothesis [14]. We consider a system described by some Hamiltonian H. An observable of interest represented by a Hermitian operator O gives rise to a corresponding autocorrelation function where O(t) = e iHt Oe −iHt is the time-dependent operator in the Heisenberg picture ( = 1). In the following, it is convenient to work directly in Liouville space, i.e., the Hilbert space of operators, and denote its elements O as states |O). Elements of the Liouville space evolve under the Liovillian L = [H, · ], i.e., |O(t)) = e iLt |O), similar to wave functions that evolve under the Hamiltonian H.
Using the Liovillian superoperator, the autocorrelation function may be written as C(t) = (O|e iLt |O). The Liouville space is equipped with an infinite-temperature , which induces a norm via ||O|| = (O|O). In the following, the central object of interest is the Liouvillian L represented in a particular basis {|O n )}, the so-called Krylov basis. The Krylov basis is routinely constructed as part of the Lanczos algorithm. In this basis, which is determined by some "seed" observable O, the representation of the Liouvillian is tridiagonal. To initialize the algorithm, we take the normalized state |O 0 ) = |O), i.e., (O|O) = 1, and set b 1 = ||LO 0 || as well as |O 1 ) = L|O 0 )/b 1 . Then we iteratively compute The tridiagonal representation of the Liouvillian in the Krylov basis {|O n )} is then given by where the Lanczos coefficients b n are real, positive numbers output by the algorithm. For our purposes, it is sufficient to assume a finite-dimensional space. Then, the algorithm halts at step n = d+1, where d is the dimension of the Liouville space, and L is a d × d-matrix. Rewriting the Heisenberg equation of motion in the Krylov basis yields where we defined ϕ n := i −n (O n |O(t)). The initial condition is given by ϕ n (0) = δ n0 and both ϕ n and b n are set to zero by convention when n is "out of bounds". The above Eq. (4) takes the form of a discrete Schrödinger equation and can be numerically solved by familiar means of, e.g., exact diagonalization or iterative schemes like Runge-Kutta and Chebyshev polynomials. As mentioned, the Lanczos coefficients b n can be interpreted as hopping amplitudes in a tight-binding model. Then, the correlation function coincides with the amplitude of the first site, i.e., C(t) = ϕ 0 (t).
For later reference, we introduce the spectral function Φ(ω) as the Fourier transform of the correlation function, i.e., It can be shown that the Lanczos coefficients b n appear in the continued fraction expansion of Φ(ω).
Consequently, there exists a (non-linear) one-to-one map between the Lanczos coefficients b n and the autocorrelation function C(t). Thus, a set of b n 's uniquely determines C(t) and vice versa. Lastly, we present the universal operator growth hypothesis as brought forth in Ref. [14]. The hypothesis concerns the asymptotic behavior of the Lanczos coefficients b n and basically states that in generic, non-integrable systems the Lanczos coefficients of local, few-body observables grow asymptotically linear, i.e., above some n the growth is given by where α > 0 and γ are real constants and o(g n ) denotes some real sequence f n with lim n→∞ |f n /g n | = 0. In the special case of a one-dimensional system, the asymptotic growth is sub-linear due to an additional logarithmic correction [14].

III. PROBING FOR STABILITY
In this section, we devise our strategy to probe the stability of classes of dynamics with respect to certain types of perturbations.

A. Dynamics and Lanczos coefficients
The preliminary section already touched on the relation between a correlation function C(t) and the corresponding Lanczos coefficients b n . As mentioned, this correspondence is one-to-one, thus, a set of b n 's uniquely determines C(t) and vice versa. However, the convergence can be quite subtle, i.e., there may be fairly similar dynamics with vastly different Lanczos coefficients. On the other side, similar Lanczos coefficients may lead to quite different dynamics.
In the numerical experiments below, we proceed as follows. First, we choose a correlation function C (t) whose stability we want to probe. We may specify C (t) as a concrete analytical function, or as a member of a class of functions, e.g., exponential decays. The exact Lanczos coefficients corresponding to C (t) are denoted by b n . These coefficients are, in general, unknown and obtaining them is quite a difficult task. Our objective is to find Lanczos coefficients b n corresponding to a correlation function C(t) that should be practically indistinguishable from C (t) for all intents and purposes [C(t) C (t)]. These "approximate" coefficients b n may be obtained in three different ways: i. The exact coefficients b n may be analytically known. In this case b n = b n as well as C(t) = C (t).
ii. The coefficients may be obtained via an educated guess. In this case, the b n may be quite different from b n while still C(t) C (t).
iii. The coefficients may be "reverse-engineered" from the correlation function C (t). For details see App. A. In the following, we will omit the technical distinction between C(t) and C (t), since we assume that Lanczos coefficients can be found for which the two correlation functions are practically indistinguishable.

B. Design of the Lanczos coefficients
In the following, we compare dynamics of correlation functions C A (t) and C B (t) with respect to their stability under a certain class of perturbations. Corresponding sets of Lanczos coefficients b A n and b B n may be obtained by one of the three methods mentioned in the last section. We demand that these Lanczos coefficients satisfy the following [partly physically motivated] requirements sufficiently well: i . The Lanczos coefficients b X n should reproduce the respective correlation function C X (t) to an acceptable accuracy, where X = A, B (this is sort of obvious and has already been addressed above).
ii . The Lanczos coefficients b X n should comply with the universal operator growth hypothesis, i.e., they should eventually attain linear growth.
iii . The resulting correlation functions C A (t) and C B (t) should decay on more or less the same time scale. This is to ensure a fair comparison between the two, since faster dynamics are typically less affected by perturbations than slower ones. iv . The Lanczos coefficients b A n and b B n should be similar in magnitude. In particular, we demand that the quan- We show in App. B that this sum is related to the spectral variance of the Hamiltonian. Hence, this condition fosters the notion that the two correlation functions C A (t) and C B (t) originate from different observables while the underlying Hamiltonians are quite similar. In practice, a value of close to unity is desirable.

C. Design of the perturbations
The considered perturbations are designed on the level of Lanczos coefficients. The coefficients b n corresponding to some correlation function C(t) will be slightly altered according tob where the perturbed coefficients are denoted byb n . Here, λ is the perturbation strength and v n is specified below. We show in App. C that this particular form of the perturbation yields a sensible scaling of the perturbed Hamiltonian with λ.
In general, it is an intricate problem to determine how a perturbation in the form of v n corresponds to a perturbation V on the level of Hermitian matrices. This is an issue that we do not attempt to tackle in this paper.
Here, we try to make as few assumptions as possible and model the perturbation v n as random numbers [with the only restriction of a minimal correlation length, see below]. Concretely, we set where the x k , y k are real, random numbers from a Gaussian distribution with zero mean and unit variance. They are normalized as k x 2 k + y 2 k = 1. The sum in Eq. (10) is capped at a number N f . This corresponds to excluding shorter wavelength or higher frequencies, which induces a minimal correlation length in the coefficientsb n [reddish noise]. No truncation at all, i.e., N f = d, would be equivalent to adding uncorrelated random numbers [white noise] to the coefficients b n . This choice can be justified by invoking the tight-binding model interpretation in which the Lanczos coefficients represent hopping amplitudes between neighboring sites. Simply altering the hopping amplitudes in a random, uncorrelated manner leads to localization effects similar to Anderson localization [17]. Through unsystematic varying of N f , we find that the transition from unlocalized to to localized behavior is quite sharp. Since the aim is to study the stability of certain relaxation dynamics, we choose N f as large as possible, but small enough to avoid said localization effects. In practice, this amounts to N f ≈ d/3. In Sec. IV C, we ease this restriction to identify "pathological" perturbations.

D. Assessing stability
Once the perturbation has acted on the coefficients b n and yielded the perturbed coefficientsb n , we can calculate the corresponding perturbed dynamicsC(t) as laid out in Sec. II. Now, some tools are needed in order to judge how strongly the dynamics were altered by the perturbation. In particular, it is important to assess to what extent the perturbed dynamics still falls into the original class of functions to which the unperturbed dynamics belonged. For example, we may ask ifC(t) is still exponential, given that the original dynamics C(t) was exponential [the decay constants may differ]. To this end, we attempt to fit the perturbed dynamics with a function that also described the unperturbed dynamics, e.g., we may try to fitC(t) with f (t) = Ae −µt . In practice, the fit is obtained via a standard fitting routine. Before we introduce quantifiers that measure the effect of the perturbation, some remarks on the fitting ansatz are in order. Correlation functions are symmetric in time, thus, all odd moments necessarily vanish. In particular, correlation functions have zero slope at t = 0. Further, we normalize the correlation function to unity at t = 0. However, for fitting we choose functions that lack these features, i.e., the above exponential ansatz may yield an A that slightly differs from unity. Further, the slope of an exponential at t = 0 is never zero. This "negligence" is due to the fact that we are more interested in an overall description of relaxation dynamics, rather than in the details of the short-time behavior. We assess the quality of a fit f (t) by calculating "how far off" it is from the given perturbed dynamics. Concretely, we define a measure of the "error" or "deviation" by the expression Here, the respective functions are evaluated at available points in time t n = nδt, where δt is the time step used to solve the equation of motion. The upper bound N eq corresponds to a time at which the dynamics in question has seemingly equilibrated [18].
To get a feeling of which numerical value of constitutes a good or bad fit, we refer to the exemplary numerical data displayed in Figs. 3, 4, 8 and 9. We introduce a second quantifier σ that measures how strongly the unperturbed dynamics are altered due to the perturbation in the first place.
The construction of this quantity is similar to the construction of the quantity that measures the fit quality.

IV. NUMERICAL ANALYSIS
In this section, we apply the presented strategy to certain relaxation dynamics. In particular, we investigate and compare the stability of two classes of dynamics. The first class in question is the class of damped oscillations, whose damping is due to an ordinary exponential factor. This class also includes simple exponential decays, which can be viewed as damped oscillations with zero frequency. These dynamics are ubiquitous in nature. For example, slow exponential dynamics may commonly arise whenever a system interacts weakly with an environment or whenever long-wavelength Fourier components of spatial densities of conserved quantities are considered. Further, exponentially damped oscillations are routinely observed as short-wavelength components. The possible choices for the second class of dynamics are of course manifold. Since an exhaustive analysis is impossible, we decide on a "Gaussianized" version of the first class. This is supposed to mean that the class includes Gaussian decays as well as oscillations damped by Gaussian factors. The restriction of the second class of course stifles any aspiration of generality. Thus, the following results should be viewed as mere numerical experiments that corroborate the existence of stability in the first class. Concretely, in Sec. IV A, we compare exponential and Gaussian decay. In Sec. IV B, the comparison is between two damped oscillations, one with an exponential damping factor and the other with a Gaussian damping factor. Lastly, in Sec. IV C, we identify "pathological" perturbations that destroy the considered dynamics. Throughout the next sections, we set the dimension of the Liouville space to d = 10000 and the frequency cut-off to N f = 3333 ≈ d/3. We repeat the strategy N = 1000 times by drawing random numbers x k , y k and present statistics for the quantifiers and σ.
A. Exponential vs. Gaussian decay In the first round of our stability investigation, an exponential decay competes against a Gaussian decay. Following the strategy laid out in the previous sections, the first step is to find suitable Lanczos coefficients that comply with the four conditions presented in Sec. III B. We begin with the Gaussian decay. The Lanczos coefficients that exactly correspond to Gaussian decay of the form C(t) = e −t 2 /2 are analytically known [method one in Sec. III A)], i.e., b n = √ n [19]. Nevertheless, the second condition, the compliance with the operator growth hypothesis, needs to be satisfied. To this end, we choose a cutoff-point n at which the coefficients are continued in a linear fashion. This yields the following Lanczos coefficients ["g" = "Gaussian"] b g n = √ n , n ≤ n αn + γ , n > n . The parameters α and γ are determined by the requirement of a smooth transition from square-root growth to linear growth, i.e., α = 1/(2 √ n ) and γ = √ n /2. The change from square-root to linear growth for n > n does not strongly affect the "Gaussianity" of the correlation dynamics [compare Fig. 2] such that condition one remains fulfilled. Next, coefficients for an (approximate) exponential decay need to be found. As mentioned, a correlation function can never be truly exponential, since it necessarily features zero slope at t = 0. Thus, we only require the exponential decay to be present after a short Zeno-time, which is usually exceedingly short compared to the relaxation time [20,21]. We achieve an approximate exponential decay via an educated guess [method two in Sec. III A]. Slow dynamics are characterized by a relatively small first coefficient b 1 , followed by a jump to a larger remainder of coefficients b n≥2 . Thus, we make the following ansatz for the coefficients ["e" = "exponential"] b e n = a , n = 1 αn + γ , n ≥ 2 .
The parameter a is set to 1.2, a similar magnitude as the first Gaussian coefficient b g 1 . The other parameters α and γ are the same as in the Gaussian case. Thus, condition two, i.e., compliance with the universal operator growth hypothesis, is fulfilled. The Lanczos coefficients as defined in Eq. (13) and Eq. (14) are depicted in Fig. 1. Since both sets of coefficients coincide for n > n , condition four is satisfied as the quantity q(b g n , b e n ) = 0.999993 is sufficiently close to unity. The corresponding dynamics are depicted in Fig. 2. The ansatz for the coefficients b e n in Eq. (14) indeed yields a nice exponential decay. A fit of the form Ae −λt with A = 1.02 and λ = 0.24 captures the dynamics quite well. Thus, condition one is sufficiently fulfilled for the exponential decay as well. Further, we can extract from Fig. 2 that both dynamics decay on more or less the same time scale. Therefore, we also view condition three as satisfied. If any, the Gaussian dynamics is faster and thus less prone to perturbations. Now that we have checked all four conditions presented in Sec. III B, we are ready to apply the perturbation as laid out in Sec. III C. We set λ = 0.5. Three exemplary perturbed dynamics with respective fits for the exponential case are depicted in Fig. 3. It is evident that the perturbed dynamics are still nicely described by an exponential, only the decay constant changes. For the Gaussian decay, three exemplary perturbed dynamics with respective fits are depicted in Fig. 4. Two displayed perturbed Gaussian curves feature oscillations, which can not possibly be captured by a Gaussian fit ansatz. The three depicted fits are much worse than in the exponential case. Insets show corresponding Lanczos coefficients.  The impression from these six exemplary curves is confirmed in Fig. 5, which shows histograms of deviations of the fits from the respective perturbed dynamics. The symbol Ω denotes the number of values within a bin of size 5 · 10 −4 . There is a clear division between the exponential cases (blue) and the Gaussian cases (red). The deviations of the exponential decays are much smaller than those of the Gaussian decays. In particular, the mean value¯ e = 0.002 (indicated by a blue, dashed line) in the exponential case is about twenty times smaller than in the Gaussian case¯ g = 0.042. For a visualization of this comparison see Fig. 3 and Fig. 4, whose exemplary curves feature values of close to the respective mean values. We have to be aware that this apparent stability of the exponential may be due to the possibility that exponential decays are generally less affected by our constructed perturbation than the Gaussian decays. To check this possibility and to show that this is not the case, we consider the quantifier σ, which measures the difference between the perturbed and the unperturbed dynamics. The inset of Fig. 5 depicts a scatter plot of all 1000 pairs (σ i , i ). If one decay would be consistently less affected than the other, a cluster of points to the left [small values of σ] would emerge. This is evidently not the case, as the distribution along the horizontal axis is relatively similar for both decays. The existence of the diagonal edge is expected, since there can be no fit that is "farther away" from the perturbed dynamics than the unperturbed dynamics itself, which is, of course, also an eligible candidate for fitting. Thus, we conclude that exponential decays is indeed stable with respect to perturbations [as in Eq. (10)]. In contrast, Gaussian decays seem to be quite unstable. This is the first main result of the paper at hand.

B. Damped oscillations: exponential vs. Gaussian
In the second round, we investigate and compare two types of damped oscillations, one with an exponential damping factor, the other with a Gaussian damping factor. As in the last section, first, two sets of Lanczos coefficients need to be found that satisfy all four conditions imposed in Sec. III B. We again begin with the Gaussian case. The coefficients b gdo n ["gdo" = "Gaussian damped oscillation"] that correspond to an oscillation damped by a Gaussian factor a neither analytically known, nor can we make an educated guess. This only leaves the third method, in which we "reverse-engineer" the coefficients from the dynamics itself. To this end, we choose a particular correlation function C(t) = e −t 2 /8 cos(2t) and proceed as detailed in App. A. This procedure allows us to obtain about 50 Lanczos coefficients. After that, the coefficients are continued "by hand" in a manner that first, respects the pattern exhibited by the coefficients so far, and second, eventually becomes linear. For the exponentially damped oscillation we can again make an educated guess [method two] for the coefficients b edo n ["edo" = "exponentially damped oscillation"]. We set the first two coefficients b 1 , b 2 to some values (b 1 = 2.0, b 2 = 1.6), followed by a jump to larger coefficients b n≥3 , which then grow in a linear fashion. The slope is determined by and coincides with the slope in the Gaussian case. The unperturbed Lanczos coefficients and corresponding correlation functions are depicted in Fig. 6 and Fig. 7, respectively. The "reverse-engineered and continued by hand" coefficients b gdo n indeed reproduce the given correlation function C(t) = e −t 2 /8 cos(2t). Further, the guess for b edo n yields a nice exponentially damped curve, i.e., the fit ansatz with [A = 1.04, µ = 0.57, ω = 0.32, φ = 2.19] nicely captures the dynamics. Therefore, condition one is satisfied for both cases. The two sets of coefficients were designed to eventually attain linear growth, thus condition two is fulfilled. It is evident from Fig. 7 that both dynamics decay on similar time scales, rendering condition three satisfied. Lastly, the quantity q(b gdo n , b edo n ) practically equals unity. Thus, all four conditions are met.  Next, we switch on the perturbation and set λ = 0.1. Exemplary perturbed dynamics are displayed in Fig. 8 and Fig. 9, respectively. The perturbed dynamics are fitted with the 4-parametric ansatz C(t) = Ae −µt cos(ωt − φ) in the exponential and C(t) = Ae −µt 2 cos(ωt − φ) in the Gaussian case. While the perturbed coefficients in the exponential case still lead to correlation functions that are within the class of exponentially damped oscillations [with different decay constants µ and frequencies ω], the same can not be said for the Gaussian case, where the perturbed dynamics decay too slowly and can not by captured by the fit ansatz above. The impression from the six exemplary curves is confirmed in Fig. 10, which displays histograms of the quantifier . Again, the division into two distinct peaks is evident. The respective mean values gdo = 0.026 and edo = 0.006 are indicated as dashed, vertical lines and differ by a factor of about five. For a visualization of what these values mean for the fit quality, see Fig. 8 and Fig. 9, whose curves feature values of close to the respective mean values. Thus, we conclude that the exponentially damped dynamics are quite stable and, in particular, more stable than the oscillations damped by a Gaussian. The disparity between the two is not as strong as for the decays in the previous section [however, different values of λ may not be directly comparable.] The inset of Fig. 10 again shows a scatter plot of all points (σ i , i ). The blue cluster of dots is a little more concentrated to small values of σ, indicating that the exponentially damped oscillations are a little less altered by the perturbation than their Gaussian counterpart. However, the marginal distribution over σ is still less partitioned into two peaks than the one over [which is the data in the histogram itself]. Hence, the stability of exponentially damped oscillations is still due to the nature of the dynamics and not due to a smaller effect of the perturbation on the dynamics. We conclude that exponentially damped oscillations are stable with respect to perturbations [as in Eq. (10)]. In contrast, oscillations damped by a Gaussian seem to be quite unstable. This is the second main result of the paper at hand.

C. Pathological perturbations
In this section, we lift the restriction imposed on the v n earlier, which was that frequencies were cut at N f ≈ d/3. Instead, we include all frequencies in the construction of the perturbation, which amounts to setting N f = d. This case simply corresponds to adding uncorrelated random numbers to the unperturbed coefficients b n . Otherwise, the strategy is pursued in the same manner as above. Three exemplary dynamics for exponential decays are depicted in Fig. 11 [for conciseness, we refrain from showing more exemplary data for the other cases]. As is evident, the exponential decay is completely absent and replaced by quite irregular dynamics, which do not seem to reach an equilibrium [N eq is set as large as possible [18] in these cases]. These curves hint at the presence of localization effects, which prevent the "particle" [in the tight-binding picture] to leave the first site. A histogram of the Gaussian vs. exponential data can be viewed in Fig. 12. The mean values read¯ g = 0.15 in the Gaussian case and¯ e = 0.12 in the exponential case, which is some orders of magnitude larger than before when shorter wavelengths were excluded. Further, the accumulation of points along the diagonal in the inset suggests that the curves no longer resemble their original form, which is expected judging from Fig. 11. The histogram for the oscillating cases is displayed in Fig. 13. The mean deviations read¯ gdo = 0.035 and¯ edo = 0.025. These values are comparable to the Gaussian case in Fig. 10. However, since the quantifier measures more or less the deviation "per time step", the comparison of values¯ between dynamics that do and do not equilibrate can be void of meaning. In practice, judging from all figures displaying exemplary curves, a value of about one percent corresponds to fits that "look good to the eye". Based on these observations, we specify a property of perturbations that lead to non-generic or pathological dynamics. Namely, a perturbation is "untypical", if it gives rise to Lanczos coefficients that vary non-smoothly, i.e., there is no minimal correlation length within the coefficients. Again, what this means on the level of Hermitian matrices is difficult to assess and beyond the scope of this work. This is the third main result of the paper at hand.

V. CONCLUSION
In this paper, we performed numerical experiments to probe the stability of two classes of relaxation dynamics. The first class consisted of exponentially damped oscillations, which also includes exponential decays. The second class was chosen as a Gaussian counterpart to the first class, i.e., including Gaussian decays and oscillations damped by a Gaussian. The whole strategy was formulated in the framework of the recursion method, in particular, the perturbations were constructed as an alteration of the Lanczos coefficients. Unperturbed coefficients b n and perturbation v n were chosen to satisfy certain physically motivated conditions.
The first main message of the paper at hand is that the exponential class of dynamics is relatively stable under the considered perturbations. In contrast, the Gaussian counterpart is found to be quite unstable. These findings confirm and extent upon previous results based on random matrices [8], which did not comply with the operator growth hypothesis. We want to emphasize that within this work the main focus should be put on the stability of the former class of relaxation dynamics, as these are ubiquitous in nature, examples are given at the beginning of Sec. IV. The investigation of the latter class should just be taken as an exemplary comparison. In fact, the choice of a Gaussian counterpart to the first class is quite arbitrary. Any number of relaxation dynamics could have been investigated instead. The second main message of the paper at hand concerns the nature of the perturbations themselves. Not only, but also in the context of the works on "typicality of perturbations" [9][10][11], it could be interesting to find properties of perturbations that lead to non-generic, "pathological" dynamics. Here, we identified such a criterion, namely that the perturbation should yield Lanczos coefficients whose minimal correlation length is still above some threshold value. As is evident from Fig. 11, including short correlations seems to lead to unorthodox dynamics. The displayed curves do not seem to reach an equilibrium, at least on the available time scale. Rather, localization-like effects are introduced, which cause some part of the wave function to remain on the first site. On the other hand, sufficiently smooth coefficients alter the relaxation dynamics in a "controlled yet non-trivial" manner. Perturbations in various numerical investigations based in random matrices [8,9] as well as spin lattice models [22,23] seem to naturally possess the above specified property. Hence, a more systematic investigation on the existence of a minimal correlation length in the coefficients for realistic setups could be a possible prospect for future research.