Modification of quantum many-body relaxation by perturbations exhibiting a banded matrix structure

We investigate how the observable relaxation behavior of an isolated quantum many-body system is modified in response to weak-to-moderate perturbations within a nonperturbative typicality framework. A key role is played by the so-called perturbation profile, which characterizes the dependence of the perturbation matrix elements in the eigenbasis of the unperturbed Hamiltonian on the difference of the corresponding energy eigenvalues. In particular, a banded matrix structure is quantitatively captured by a perturbation profile which approaches zero for large energy differences. The temporal modification of the relaxation is linked to the perturbation profile via a nonlinear integral equation, which admits approximate analytical solutions for sufficiently weak and strong perturbations, and for which we work out a numerical solution scheme in the general case. As an example, we consider a spin lattice model with a pronounced banded matrix structure, and we find very good agreement of the numerics with our analytical predictions without any free fit parameter.


I. INTRODUCTION
Despite their microscopic chaoticity [1,2] the macroscopically observable behavior of isolated quantum manybody systems is often surprisingly regular. For instance, it is by now well-established that these systems generically equilibrate and usually even thermalize [2][3][4], and that the approach to equilibrium quite often follows a rather simple and direct route. Understanding how this dynamics emerges from a microscopic description, however, is still a theoretical challenge attracting considerable attention recently. A particularly interesting question in this context is how the observable relaxation behavior of a given system is modified under the influence of reasonably weak perturbations, linking, for example, analytically tractable simple systems (e.g. noninteracting, integrable) to generic ones (e.g. interacting, nonintegrable).
Characterizing the response of a given system to a perturbation is a recurrent problem in many areas of physics. Arguably the standard approach is to expand the pertinent equations of motion in terms of the perturbation strength and to solve the resulting hierarchy of simplified equations iteratively. Unfortunately, such a strategy is doomed to failure in the case of quantum systems with many degrees of freedom. Due to their extremely dense energy spectra, the concomitant small denominators of a perturbative expansion limit its applicability to extremely short time scales much smaller the observed relaxation times. While there is strong evidence that related concepts like Fermi's golden rule and linear response theory can describe many-body dynamics in certain scenarios [5][6][7], this somewhat surprisingly holds despite the many-body character and not because of it.
Here we tackle the question of how a many-body systems responds to perturbations by "nonperturbative" methods, namely a typicality approach that aims to extract and separate the macroscopically relevant perturbation characteristics from the huge number of micro-scopic degrees of freedom. Our starting point is an isolated many-body quantum system described by a timeindependent reference Hamiltonian H, and prepared in some initial state far from equilibrium. Provided that we know the observable relaxation dynamics of this unperturbed reference system, we ask how the behavior is changed when adding a weak-to-moderate perturbation λV . In other words, the system still starts from the same initial state, but now evolves in time according to the perturbed Hamiltonian One situation that could be modeled by such an approach is an unperturbed system composed of two isolated subsystems at equilibrium, which are then coupled sufficiently weakly via the perturbation V and relax to a new, joint equilibrium state. Another interesting scenario arises when the reference system H is integrable, in which case one can often calculate the unperturbed behavior analytically. In particular, integrable systems usually still equilibrate (just like the nonintegrable ones), meaning that expectation values of experimentally relevant observables approach a constant value and stay there for most of all later times. However, these integrable systems (unlike the nonintegrable ones) may not thermalize, i.e., equilibrium expectation values are not described by the pertinent thermodynamic equilibrium ensemble and call for extensions like generalized Gibbs ensembles instead [8][9][10][11]. Adding a small integrability-breaking perturbation commonly leads to "prethermalization" [4,6,[12][13][14][15], meaning that the system still follows the unperturbed (nonthermalizing) behavior for quite some time before eventually departing towards the associated thermal state. As a third example, more generally, one may think of the unperturbed system as some system for which the relaxation dynamics happens to be known, and ask for the behavior when changing some parameter of the Hamiltonian (e.g. a "quantum quench" [3,11,16]). Basing our analysis on previous results from Ref. [17], we recap those findings in Secs. II and III. More pre-cisely, we introduce the considered classes of systems and formulate the key assumptions of our theory in Sec. II, and establish the announced theoretical prediction of the many-body response in Sec. III. A crucial role is played by the resolvent (z − H λ ) −1 averaged over an ensemble of perturbed Hamiltonians H λ , whose computation we expound in Sec. IV. These results are then used in Sec. V to make the prediction from Sec. III explicit and to compare it to numerical examples for random-matrix and spin models. Finally, we summarize and discuss our results in Sec. VI.

II. SCOPE AND PREREQUISITES
Before presenting our main result, we introduce the setting and collect several key assumptions about the physical situations we aim to describe (see also Supplemental Material of Ref. [17] for further technical details).
The isolated many-body quantum system of reference is described by a time-independent Hamiltonian H = ν E ν |ν ν| and is prepared in some (pure or mixed, and generally far from equilibrium) initial state with density operator ρ(0). According to textbook quantum mechanics, the state at any later time is then given by ρ(t) = e −iHt ρ(0)e iHt ( = 1). Of primary interest to us are the time-dependent expectation values A ρ(t) := Tr[ρ(t)A] of self-adjoint operators A which model some experimentally or theoretically relevant observable, such as (sums of) local and/or few-body operators [2][3][4]. Similarly, the time-evolved state of the perturbed system with Hamiltonian H λ from (1) is given by ρ λ (t) = e −iH λ t ρ(0)e iH λ t .
Overall, the main objective of our present work is to establish quantitative predictions for the perturbed dynamics A ρ λ (t) based on the unperturbed behavior A ρ(t) and some essential characteristics of the perturbation V .
Regarding the systems under study, the following four key assumptions will be taken for granted hereafter: (i) The system should exhibit a well-defined macroscopic energy, implying that the initial state ρ(0) (and hence also ρ(t) at any later t) only significantly populates levels E ν ∈ I within a macroscopically small energy window I := [E, E + ∆]. In particular, it is assumed that the density of states (DOS) is approximately constant throughout I, D(E) ≈ ε −1 with the mean level spacing ε. At the same time, the system's many-body character entails that the window I is still microscopically large in the sense that the number of levels contained in I is exponentially large in the system's degrees of freedom [18]. We emphasize that the initial state does not necessarily define the window I. On the contrary, the window I can be a to some extent arbitrary interval with the two prerequisites that it exhibits an approximately constant DOS and contains all E ν with nonnegligible level populations ν|ρ(0)|ν .
(ii) The perturbation should be sufficiently weak so as to leave the thermodynamic properties of the system basically unchanged. Notably, phase transitions induced by the perturbation are thus ruled out. Recalling from textbook statistical mechanics that the DOS in (2) is related to the Boltzmann entropy S(E) via D(E) = e S(E)/kB S ′ (E) with Boltzmann's constant k B , this assumption particularly implies that also the DOS of the perturbed H λ remains approximately constant with mean level spacing ε [cf. assumption (i)]. Due to the generic level repulsion of interacting many-body systems [1], the spectrum of H λ is typically indeed rather stiff, meaning that the individual eigenvalues exhibit very fast fluctuations upon variation of λ, while their density only changes very slowly [19].
We remark that such negligible changes of the thermodynamic properties do not rule out interesting and nontrivial changes of the relaxation dynamics, notably if the unperturbed Hamiltonian is in some sense special (e.g., integrable, commuting with A or ρ(0), ...), see also the examples below Eq. (1) and in Sec. V.
(iii) The perturbation should be sufficiently strong so that it significantly mixes a large number of unperturbed levels. Denoting by |m λ the eigenvectors of H λ in (1), this is to say that the overlaps between the unperturbed and perturbed eigenvectors should extend across a scale Γ with Γ ≫ ε, i.e., U mν should be nonnegligible (in a coarse-grained sense, see below) as long as |E m − E ν | Γ (see also Eq. (7)). On the other hand, note that assumptions (i) and (ii) practically require Γ ≪ ∆ [20], where ∆ is the width of the energy window I from (2). The extreme density of levels of typical many-body systems [see below Eq. (2)] still leaves room for a large range of parameters λ such that ε ≪ Γ ≪ ∆. In particular, we can and will take for granted that the number of levels N v := Γ/ε that get mixed by the perturbation is still exponentially large in the system's degrees of freedom f [17], i.e., Without going into the details, we remark that perturbations which do not satisfy the requirement ε ≪ Γ turn out (as one might have expected) to actually be so weak that they do not notably modify the unperturbed relaxation on any reasonable time scale. Incidentally, the same behavior will also be correctly reproduced by our final results. In this sense, the requirement ε ≪ Γ is not really indispensable. So far, these considerations have been very general and did not exploit any more specific properties of the actual system at hand. To make any progress, it is clear that some information about the perturbation V and possibly also the observable A and initial state ρ(0) must be taken into account. The common lore of statistical physics fur-thermore suggests that, despite its microscopic complexity, the observable behavior of a many-body system can usually be described in terms of a relatively small number of macroscopic (coarse-grained) quantities, for instance some appropriately defined (local) densities. This brings us to our main assumption about the structure of admissible perturbations: (iv) On a coarse-grained level, the magnitude of the perturbation matrix elements V µν := µ|V |ν within the energy window I should only depend on the energy difference |E µ − E ν | of the coupled levels [21], i.e., where [ · · · ] loc denotes a local average over matrix elements corresponding to levels that are close to E µ and E ν in energy (see also Sec. V C for an explicit example). Put differently, the left-hand side in (5) is understood (and formally defined) analogously as when going over, e.g., from classical point particles to (local) particle densities, namely as the effective density of the perturbation's squared matrix elements (in modulus, and "local" with respect to the spectrum of H 0 ). Accordingly, σ 2 (E) in (5) is denoted as the perturbation profile, and is, by construction, a smooth [22] and slowly varying function of E (compared to the mean level spacing ε). Semiclassical arguments [23,24] as well as numerical evidence [25][26][27][28][29] suggest that a rather common feature of realistic perturbations is a so-called banded structure of the perturbation matrix V µν (see also Fig. 3 in Sec. V C below for a particular example). By definition, this means that the (coarse-grained) V µν indeed depend only on E µ − E ν and that the perturbation profile σ 2 (E) in (5) approaches zero for E → ∞. However, it should be emphasized that σ 2 (E) is also admitted to remain finite for E → ∞, i.e., the matrix V µν may but need not exhibit a banded structure [30]. Yet another common feature of many realistic perturbations is a so-called sparse matrix structure (large fraction of vanishing matrix elements V µν ), prominently arising, e.g., if the reference Hamiltonian H is noninteracting and V describes few-body interactions [24,28,31,32]. Again, our present approach is still compatible with a possibly (but not necessarily) sparse structure of V µν [the local average in (5) then must extend over many nonvanishing matrix elements].
Our next goal is to establish the key role of the perturbation profile (5) for the deviations of the perturbed expectation values A ρ λ (t) from the unperturbed A ρ(t) . The main idea is to consider not one particular V , but rather an entire ensemble of perturbations, all of which share the property (5) with the "true" perturbation of interest, but are otherwise unbiased and rather arbitrary. More precisely, apart from the trivial constraint V * µν = V νµ , we choose the matrix elements V µν to be independent random variables following a probability distribution where [ · · · ] V denotes the average over the ensemble of perturbations, and {f E (v)} E>0 is a family of probability densities on R or C with mean zero and variance σ 2 (E). Likewise, f 0 (v) is a probability density on R of vanishing mean and finite variance [22]. Note that the ensemble thus satisfies (5) in an ergodic sense, i.e., when replacing local averages [ · · · ] loc by ensemble averages [ · · · ] V .
To arrive at a prediction for the perturbed dynamics A ρ λ (t) , we first evaluate the average behavior [ A ρ λ (t) ] V over all members of the considered ensemble of perturbations. Second, we consider the deviations ξ V (t) : ] V for one particular realization from the average. It turns out [17] that the variance [ξ V (t) 2 ] V is inversely proportional to the number N v of unperturbed levels mixed by the perturbation from assumption (iii). Exploiting (4), we can therefore conclude that, for the overwhelming majority of individual perturbations in the considered ensemble, the actual behavior so that the latter in fact correctly describes the dynamics under nearly all perturbations of the ensemble for sufficiently large system sizes. Results of this kind are also commonly known as "typicality", "concentration of measure", or "ergodicity" properties [2][3][4].
Taking for granted that the perturbation profile (5) is indeed the essential quantity for deviations between the perturbed and unperturbed systems, we may expect that also the behavior of the true system of interest should follow the ensemble average. Unfortunately, it is hard to prove this for any given, concrete physical system. Nevertheless, a phenomenological justification by means of examples is possible for a variety of different models [17,33], see also Sec. V below. For the rest, we observe that the probability distribution (6) is still rather arbitrary since we only fix the first two moments of the densities f E (v). In principle and if available, additional information about the distribution of the true V µν could thus be incorporated when choosing the f E (v), but similarly as in the central limit theorem, these statistical properties turn out to be practically irrelevant, reinforcing the pivotal role of the second moment (5).
To conclude this section, we remark that the true perturbation will usually exhibit correlations (i.e., functional interdependencies) between the matrix elements V µν , which may arise, for example, due to the locality and few-body character of interactions [34,35]. Since such correlations are not accounted for in the considered perturbation ensembles, it is implicitly assumed that their effect on the dynamics is negligible. In practice, this particularly means that the reference Hamiltonian H should be sufficiently "clean" such that the individual terms constituting the perturbation V are in some sense "orthogonal" to those of H. Notably, this rules out the possibility to "reverse the roles" by defining a new reference Hamiltonian H ′ := H λ and considering a per-

III. TYPICAL PERTURBED RELAXATION
Given the prominent role of the Hamiltonian as the generator of time evolution, it will be no surprise that the transformation matrices U mν between the eigenbases of H and H λ [see Eq. (3)] are of particular importance to relate the unperturbed and perturbed dynamics. Especially relevant turns out to be the so-called overlap distribution u(E), which describes the squared magnitude of the U mν averaged over the considered ensemble of perturbations, Due to the (approximate) constancy of the level density [assumptions (i) and (ii)] and the fact that the statistics of the V µν in (6) only depend on E µ − E ν , it follows that the statistics of the U mν from (3) must be translationally invariant in energy, and hence the second moment in (7) must only depend on the energy difference E m − E ν . Referring to Ref. [17] for the details, the typicality approach outlined below Eq. (6) then eventually yields that, for the overwhelming majority of perturbations in any admissible ensemble (6), the perturbed time evolution is given by We recall that A ρ(t) is the reference dynamics observed under the unperturbed Hamiltonian H. Furthermore, the density operatorρ appearing on the righthand side of (8) is defined via its matrix elements In other words,ρ may thus be viewed as the unperturbed diagonal ensemble associated with the initial state ρ(0) which is in addition locally washed out via the functionũ(E), arising as the convolution of u(E) with itself. According to [3,[36][37][38], this operatorρ can usually be well approximated by the microcanonical ensemble ρ mc corresponding to the pertinent energy window I from (2). Finally, the so-called response profile g(t) on the right-hand side of (8) is the Fourier transform of u(E) from (7), In particular, it can be readily verified that |g(0)| 2 = 1 and |g(t)| 2 → 0 as t → ∞. According to (8), this function g(t) thus describes how the unperturbed behavior is modified to approach the perturbed equilibrium value A ρ , i.e., it encodes the system's response to the perturbation. The remaining task is to compute the function u(E) from (7). To this end, we introduce the resolvent G(z) := (z − H λ ) −1 of H λ , which encodes the overlaps on the lefthand side of (7) as |U mν [1,39]. Since the ensemble average of G(z) can be written as with the scalar function G(z) defined in a minute, we can exploit D(E m ) ≈ ε −1 [cf. assumptions (i) and (ii)] to arrive at Finally, the above introduced ensemble-averaged resolvent G(z) itself can be obtained as the solution of the following nonlinear integral equation [17,24], In summary, the strategy to obtain a prediction for the perturbed relaxation thus is to follow the sequence of equations (8)- (11) in reverse order: First, for a given perturbation profile σ 2 (E) and perturbation strength λ, we solve Eq. (11) for G(z). Second, this gives us access to the overlap distribution u(E) via Eq. (10). Third, evaluating its Fourier transform (9) we obtain the response profile g(t), which then allows us, fourth, to predict A ρ λ (t) from the unperturbed A ρ(t) according to Eq. (8). Clearly, the first step, namely to solve the nonlinear integral equation (11), is the most demanding task. This problem is at the focus of the next section.

IV. EVALUATION OF THE ENSEMBLE-AVERAGED RESOLVENT
In this section, we will discuss solutions G(z) of Eq. (11) and the resulting overlap distributions u(E) from (10). We first consider in Sec. IV A two limiting cases for which analytical approximations will be obtained. Thereafter, we elaborate on how to solve Eq. (11) in the intermediate regime numerically using pseudospectral Chebyshev expansions [40,41]. The evaluation of the predicted dynamics and its comparison with explicit examples is deferred to the ensuing Sec. V.

A. Analytically tractable special cases
According to assumption (iv) from Sec. II, the perturbation profile σ 2 (E) from (5) is a well-behaving (continuous) function, so that the quantitȳ exists [22]. Essentially,σ thus characterizes the "intrinsic strength" of the perturbations V . As explained in Sec. II, the perturbation matrix V µν in the eigenbasis of the unperturbed Hamiltonian H is often expected to exhibit a banded structure, meaning that its perturbation profile σ 2 (E) approaches zero as E → ∞. The corresponding so-called "band width" or "perturbation range" may thus be quantified by However, in full generality we will also admit cases where σ 2 (E) does not approach zero for large E. In such a case, but also when σ 2 (E) only decays very slowly with E, the band width ∆ v will be infinite. Our first approximation starts from the observation that if the perturbation is sufficiently weak [sufficiently small λ in (1)] then also the mixing of eigenvectors between the unperturbed and perturbed Hamiltonians should be weak in the sense that the concomitant eigenvector overlaps (7) are only non-negligible for small energy differences E m − E ν of the corresponding eigenvalues. In view of (10), we therefore inspect the case that the function G(z − E) in the integrand in (11) exhibits (as a function of E, and for any preset z of later relevance) a very narrow peak compared to variations of the perturbation profile σ 2 (E). Accordingly, the integral is dominated by the region around the maximum of G(z − E) at E ≈ |z|, and we can approximate σ 2 (|E|) by its central value σ 2 (|z|). Together with D(E) ≈ ε −1 [cf. assumption (i)] we thus obtain with Exploiting once again that G(z) exhibits a very narrow peak compared to the variations of σ 2 (|z|) implies with (12) that σ 2 (|z|) ≈σ 2 for all the relevant values of |z| for which G(z) significantly deviates from zero. Furthermore, focusing in view of (10) on arguments z of the form z = x − iη with x ∈ R, the quantity C(z) in (15) assumes the same constant value C(−iη) for all z. In other words, G(z) in (14) can be written as 1/(z − c) for some constant c ∈ C. Consequently, when evaluated in the principal value sense, C(z) in (15) only depends on the sign of the imaginary part of the denominator in (14), yielding C(z) = ∓iπ for sgn(Im z) = ±1 as the only consistent solution. Altogether, we thus arrive at the approximation and with (10) we conclude that u(E) approximately assumes the Breit-Wigner form Hence Γ quantifies the peak width of u(E), and likewise for G(z). Our initial assumption that G(z) is sharply peaked thus means that σ 2 (E) must exhibit only small changes upon variations of E on the order of Γ. Viewing the perturbation strength λ as variable and all other system properties as fixed, we may thus consider (16)-(18) as a weak perturbation approximation. Importantly, this approximation is expected to apply for practically any reasonable perturbation profile σ 2 (E) provided the perturbation strengths λ are sufficiently small. In many cases, one furthermore expects that the band width (13) at the same time quantifies the scale on which σ 2 (E) exhibits notable variations, yielding as the pertinent condition for the validity of the above approximations. On the other hand, in cases where the variations of σ 2 (E) remain relatively small for arbitrary E, those approximations will actually apply to arbitrary coupling strengths λ (apart from the general restrictions in Sec. II).
Our second approximation is similar in spirit but complementary to the first one. Namely, we follow the same reasoning as before with the roles of G(z) and σ 2 (E) reversed, i.e., we now consider the case that the perturbation profile σ 2 (E) is sharply peaked compared to the variations of G(z) for arguments of the form z = E − iη. In the integrand in (11), we thus approximate G(z −E) ≈ G(z) and (as before) D(E) ≈ ε −1 , leading to Solving this algebraic equation for G(z) and observing that sgn(Im G(z)) = − sgn(Im z) due to G(z) = [(z − λV ) −1 ] V [see above (10)], we obtain Substituting into (10), we are left with the semicircular distribution where Θ(x) denotes the Heaviside step function. The condition that σ 2 (E) is sharply peaked thus means that G(z = E −iη) must exhibit only small changes upon variations of E on the order of γ. Viewing the perturbation strength λ as variable and all other system properties as fixed, we may thus consider (21)-(23) as a strong perturbation approximation. More precisely, this approximation is expected to apply for practically any reasonable perturbation profile σ 2 (E) provided the perturbation strengths λ are sufficiently large, and provided that σ 2 (E) does approach zero for large E in the first place.
In particular, this is the case if the band width ∆ v from (13) is finite. Furthermore, if ∆ v at the same time quantifies the scale on which σ 2 (E) exhibits notable variations, then the pertinent condition for the validity of the above approximations assumes the form Essentially, the overlap distribution u(E) from (7) is thus predicted to approximately assume the Breit-Wigner form (18) under the weak perturbation condition (19), and the semicircular form (23) under the strong perturbation condition (24), largely independently of any further details of the perturbation profile σ 2 (E) from (5). In the intermediate regime, characterized by Γ ≃ γ, or equivalently one thus expects a smooth crossover between these limiting cases (see also Fig. 1 below), which will depend on the detailed shape of the perturbation profile σ 2 (E), and which in general will only be tractable by numerical means.

B. Numerical treatment of the general case
Our goal is to determine the overlap distribution u(E) according to (10) by numerically solving the nonlinear integral equation (11) for largely general perturbation profiles σ 2 (E). As in the previous subsection, we thus can and will focus in (11) on arguments z of the form z = x − iη with x ∈ R and η > 0 very small. As noted below (21), the relation G( ≥ 0 for η > 0 and vice versa, i.e., the sign of the imaginary part of G(z) jumps when crossing the real line. For purely real z, in turn, this implies that the solution of (11) becomes ambiguous, depending on whether one chooses to continue from the upper or lower half-plane. Bearing in mind that the latter option is appropriate in (10), we introduce the abbreviation Exploiting (as usual) that D(E) ≈ ε −1 [cf. assumption (i)], the integral equation (11) can thus be rewritten for real-valued x as (27) with the additional constraint that Our method of choice to solve Eq. (27) numerically is an expansion in terms of Chebyshev rational functions B n (x) (n = 0, 1, . . .), which are derived from the Chebyshev polynomials of the first kind T n (x) by a compactification of the real line, Here ℓ is an arbitrary, fixed parameter that sets the scale for compactification and should roughly reflect the typical scale of the function to be expanded for optimal convergence [41]. Hence we express where the real-valued functions G R (x) and G I (x) are truncated Chebyshev series, i.e., (31) The (real-valued) coefficients G R n and G I n are then to be determined such that G + (x) from (30) satisfies (27) and (28) "as well as possible." For given expansion coefficients G := (G R 0 , G I 0 , . . . , G R M , G I M ) and x, the residual [i.e., the violation of Eq. (11)] is defined as (30) and (31). We minimize |R(G, x)| 2 by means of pseudospectral methods [40,41], requiring Re R(G, x m ) = Im R(G, x m ) = 0 for a discrete set of real-valued collocation points x m (m = 0, 1, . . . , M ). A common choice for these x m is to use the roots of the (M + 1)th Chebyshev rational function B M+1 (x), so that the pseudospectral method coincides with a spectral expansion when an optimal Gaussian quadrature rule is used to calculate inner products numerically [40,41].
Altogether, forcing Re R(G, x m ) = Im R(G, x m ) = 0 results in a set of 2(M + 1) algebraic equations for the 2(M + 1) unknown expansion coefficients G R n , G I n ∈ R. This system of equations is then solved iteratively by the Newton-Raphson method using either of the limiting distributions (18) or (23) for the first initial guess, and gradually varying λ across the intermediate regime thereafter. If the initial guess is sufficiently close to the actual solution and satisfies Im G + (x) ≥ 0, this ensures that also the finally obtained approximation will fulfill the constraint (28).
In Fig. 1, we display the so-obtained numerical solutions u(E) in (10) for different perturbation profiles σ 2 (E) and various perturbations strengths λ along with the limiting Breit-Wigner functions (18) expected for small λ and the semicircular functions (23) expected for large λ. The selected perturbation profiles are a step function, an exponential function, and a double-Breit-Wigner function, All three perturbation profiles are also shown in the insets of the left panels in Fig. 1 (23), expected for weak and strong perturbations, respectively. As predicted below (25), the crossover between these two limits occurs around λ ≃ λc ≈ 0.05. Note that the vertical axes are scaled as indicated in the top-left corner of each panel.
crossover coupling strength in (25). Moreover, the order of the Chebyshev expansions is M = 80 throughout, with the parameter ℓ varying between 0.5 and 8 [roughly optimizing the global residual (32)].
For each of the three profiles (33)- (35), the predicted crossover from the Breit-Wigner to the semicircular shape of u(E) is clearly visible as λ is increased. The intermediate regime, where neither the Breit-Wigner nor the semicircular distribution offers a satisfactory approximation, appears to be somewhat smaller for the discontinuous step profile than for the smooth exponential and double-Breit-Wigner profiles. In any case, in this intermediate regime there is a (relatively mild) dependence of u(E) on the detailed shape of σ 2 (E). It therefore seems reasonable to expect that -at least in principle -it may be possible to reconstruct from a sufficiently precisely known function u(E) the underlying perturbation profile σ 2 (E).

V. EVALUATION OF THE RELAXATION DYNAMICS AND EXAMPLES
With our above obtained results for the overlap distribution u(E) at hand, we now turn to their implications for the response profile g(t), which governs the deviations of the perturbed from the unperturbed relaxation behavior according to (8). Specifically, we will first address in Sec. V A some more general issues, while in the subsequent Secs. V B and V C, we will compare our theoretical prediction (8) with two explicit examples of randommatrix and spin models, respectively.

A. Response profile
Exploiting in (9) our usual approximation D(E) ≈ ε −1 [cf. assumption (i)], the response profile g(t) can be read- ily obtained via Fourier transformation from our analytical and numerical findings for u(E) in the previous Sec. IV. For the two analytically tractable special cases from Sec. IV A, the Fourier transformation can again be performed analytically, whereas for the numerical solutions from Sec. IV B, also the Fourier transformation is only possible by numerical means.
In the limit of weak perturbations, when u(E) assumes the Breit-Wigner form (18), one readily finds along these lines that g(t) amounts to an exponential decay, where the rate Γ is the full width at half maximum of u(E) as defined in (17). Likewise, for (moderately) strong perturbations such that u(E) takes the semicircular shape (23), its Fourier transform is where J 1 (x) is the Bessel function of the first kind of order 1, and γ as specified in Eq. (21) is the radius of the semicircle.
In the intermediate regime, our findings for u(E) imply that g(t) must exhibit a crossover between these two limiting behaviors. Calculating the Fourier transforms of the numerical solutions for u(E) from Fig. 1, we obtain the solid curves shown in Fig. 2 for |g(t)| 2 , which is the actually relevant quantity in (8). This illustrates quantitatively the expected crossover from (36) to (37) with increasing λ.
The first general conclusion is that the perturbed relaxation becomes faster with increasing λ. Quite obviously, the underlying physical reason is a corresponding broadening of u(E) with increasing λ, which in turn indicates (as expected) that an increased number of unperturbed energy levels are coupled by the perturbation according to (7).
The second general conclusion is that the functions g(t) become independent of any further details of the perturbation profile σ 2 (E) for asymptotically large or small λ, while some (rather moderate) functional dependence on σ 2 (E) remains in the intermediate regime. Again, the underlying reasons are our analogous observations for the overlap distributions u(E) in the preceding section. Though the functional dependence of g(t), and thus of the perturbed relaxation in (8), is quantitatively rather weak, it still may be possible, at least in principle, to infer the (coarse-grained) perturbation profile (5) of the specific perturbation V for some given many-body system (1) from the observable temporal relaxation via (8).

B. Random matrix example
To verify that the theoretical prediction (8) indeed describes the behavior of many-body quantum systems (provided that assumptions (i) through (iv) from Sec. II hold), we finally compare it to explicit numerical examples.
The first example is a (in some sense artificial) random matrix model that satisfies the requirements from Sec. II by construction and thus serves as a testbed for the validity of the approximations employed in the derivation of Eq. (8) (see also Ref. [17]). The reference Hamiltonian has equally spaced energy levels E ν = νε with ε = 1/512. The perturbation V is a complex Hermitian random matrix distributed according to (6) with On average, the matrices V µν are thus sparse with a fraction p of nonvanishing entries following a complex normal distribution of varianceσ 2 (E) for µ < ν, and V νµ = V * µν . For simplicity, the diagonal matrix elements V νν are sampled similarly, but with a real normal distribution for the nonvanishing entries. Consequently, the perturbation profile (5) is given by Specifically, we implemented the three perturbation pro-files (33)-(35) withσ 2 = p = 0.2 and ∆ v = 1.46 (corresponding to about 750 levels). The initial state ρ(0) = |ν 0 ν 0 | is an eigenstate of the reference Hamiltonian H from the middle of the spectrum, and we observe its survival probability or fidelity [42,43], i.e., A = ρ(0). Hence A ρ(t) = 1 for all t while A ρ = A mc ≈ 0 for a sufficiently large energy window I from (2), so that the prediction (8) reduces to In other words, recording the dynamics in this setup for one particular perturbation sampled from (38), we should exactly recover the solid curves in Fig. 2. The dashed lines in the figure represent one such example dynamics for a Hilbert space of dimension 2 14 = 16 384 and an initial eigenstate |ν 0 with ν 0 = 2 13 = 8192. The main conclusion is that the simulation results indeed agree almost perfectly with the theoretically predicted solid curves throughout the entire crossover regime.

C. Spin lattice example
Finally, we test the theoretical prediction (8) in a more realistic two-dimensional spin-1 2 model. We consider a square lattice of L×L sites as sketched in Fig. 3(a), where the reference Hamiltonian couples nearest neighbors with an isotropic spin-spin interaction, Here σ i,j := (σ x i,j , σ y i,j , σ z i,j ) with σ α i,j denoting the Pauli matrices acting on site (i, j). The perturbation adds spinflip terms between next-nearest neighbors, In all of the numerics presented here, we used L = 4 and J = 1, and we focused on the sector with vanishing total magnetization in the z-direction.
To obtain the perturbation profile (5) of V , we first fix an energy window I by choosing the central 60 % of energy levels, which comprises a total of 7722 states ranging from E = −8.8 to E = 5.8, implying a mean level spacing ε = 0.0019. Next we compute the matrix elements V µν with E µ , E ν ∈ I by diagonalizing the reference Hamiltonian H. A coarse-grained view of the resulting matrix is shown in Fig. 3(b), visualizing the bandedness of the perturbation matrix. We proceed by binning the V µν according to the energy difference E µ −E ν of the associated levels and evaluate the average of |V µν | 2 within each bin. The obtained relation between the coarse-grained |V µν | 2 and E µ − E ν is displayed as a black curve in Fig. 3(c), indicating an approximately exponential dependency. The function σ 2 (E) is then determined by fitting the exponential form (34) to the empirical distribution, yielding the red line in Fig. 3(c) withσ 2 = 0.00502 and ∆ v = 7.32. This implies a value of λ c = 0.75 for the predicted location of the crossover (25) between the exponential and Bessel-type decay characteristics (36) and (37), respectively.
As a first observable, we investigate the magnetization correlation m c in the z direction between next-nearest neighbors from the center of the lattice, One could consider these two spins at (2,2) and (3,3) as the system and all other surrounding spins as a bath. In the reference Hamiltonian H, the system spins can thus only interact via the bath, whereas the perturbation V adds a direct interaction between them. For the initial state ρ(0) = |ψ ψ|, we choose those two system spins at (2,2) and (3,3) to be in the "up" state, while the bath is supposed to be at equilibrium, which we emulate by choosing a Haar-distributed random vector in the bath's subspace. However, to ensure assumption (i) of a well-defined macroscopic energy, we finally apply a Gaussian projection Π E,∆E of mean energy E = 0 and standard deviation ∆ E = 2 to the soobtained state, simulating a macroscopic measurement of the system energy that yielded E = 0 [44][45][46]. If |φ denotes a Haar-distributed random vector on the full (zero-magnetization) Hilbert space, we thus have with σ + i,j := σ x i,j + iσ y i,j and In Fig. 4(a), we compare the observed dynamics obtained by exact diagonalization (dashed lines) with our theoretical prediction (8) (solid lines) for several perturbation strengths λ. For the theoretical prediction, we use the numerical reference dynamics (i.e., the dash-dotted black curve with λ = 0) for m c ρ(t) . The function g(t) is the Fourier transform of u(E) calculated as explained in Sec. IV B from the empirically determined approximate perturbation profile σ 2 (E), i.e., the red curve in Fig. 3(c). The so-obtained response profiles g(t) are also displayed in the inset of Fig. 4(a). Since the long-time limiting values exhibit some finite-size variations, we do not use the microcanonical value m c mc = −0.0805 (within the 60 % window I) for m c ρ , but instead compute the predicted coarse-grained diagonal ensembleρ directly as detailed below Eq. (8), making use of our solution for u(E) and the known occupations ν|ρ(0)|ν of the initial state from (44). The resulting quantitative values of m c ρ for the various perturbations strengths λ are given in the figure caption. The agreement between theory and numerics is very good despite the rather small system size and several idealizations. In particular, the assumptions of a homogeneous density of states [assumption (i)], of an exponential perturbation profile [assumption (iv) and Fig. 3(c)], and of uncorrelated matrix elements V µν [see above (6)] are all violated to some extent and are thus potential origins of the visible small deviations in Fig. 4 for short times. The fluctuations for longer times, in contrast, are likely caused predominantly by finite-size effects. We emphasize that there are no free parameters in the theoretical prediction; all ingredients in (8) were extracted directly from properties of the model (41)- (42).
As a second observable, we consider the spin-flip or hopping correlation j c between the same sites (2, 2) and (3,3) from the center of the lattice in Fig. 3(a), where σ ± i,j := σ x i,j ±iσ y i,j . For the initial state, we employ a dynamical typicality setup [47,48] to prepare the system far from equilibrium, choosing where |φ is a Haar-distributed random state as before, Π is a projector onto the central 2048 states in the zeromagnetization sector [ensuring assumption (i)], and κ is a real parameter (in the examples, we use κ = 2). A similar comparison as for m c between numerical simulations and the theoretical prediction (8) is shown for the hopping correlation j c from (46) in Fig. 4(b). In particular, the functions g(t) are the same in both panels of Fig. 4. On the other hand, in this setup j c ρ is well approximated by the thermal expectation value j c mc = 0 (by symmetry), so that we used this value throughout. Altogether, this amounts again to an entirely parameterfree prediction of the perturbed dynamics, which agrees well with the actually observed behavior.

VI. CONCLUSIONS
We investigated the response of quantum many-body systems to weak-to-moderate perturbations within a nonperturbative typicality framework. In particular, we presented a method to theoretically predict time-dependent expectation values of observables for the perturbed system from the unperturbed relaxation behavior. This prediction (8) entails that the perturbed relaxation resembles the unperturbed one, but is modified by a characteristic response profile function g(t) that pushes the system towards a coarse-grained diagonal ensemble state, which can usually be identified with the pertinent thermal state. The function g(t), in turn, is essentially determined by the perturbation profile, i.e., the locally averaged squared absolute value (5) of the perturbation's matrix elements V µν in the unperturbed basis.
For asymptotically weak perturbations, the response profile g(t) describes an exponential decay, where the decay rate corresponds to the energy scale across which the perturbation mixes unperturbed eigenstates, scaling quadratically with the perturbation strength λ. Broadly speaking, this may be understood as a nonperturbative justification of Fermi's golden rule in a many-body setting.
The nonperturbative character of our method becomes manifest as the perturbation strength is increased. Our results then predict a crossover of g(t) towards the Besseltype shape (37), whose inverse relaxation time scale γ still quantifies the mixing of energy levels, but now scales linearly with λ and additionally depends on the energy range ∆ v of the perturbation.
We verified all those theoretical predictions in an explicit example of a spin system on a 4 × 4 square lattice. Using exact diagonalization to determine the perturbation profile of the applied perturbation empirically [cf. Fig. 3(c)], the function g(t) derived from it indeed describes the actually observed perturbed dynamics remarkably well as long as the key assumptions (i) through (iv) collected in Sec. II are satisfied. Notably, the theory does not involve any free parameters, i.e., all quantities were determined first-hand from the underlying spin model. Since the perturbation profile is the only variable input for the theory, this establishes that said profile encodes the dynamical response on a fundamental level.
Then again, the correspondence between the perturbation profile and the dynamical response may in principle be exploited the other way round, too. The rapidly improving experimental capabilities to observe time-dependent expectation values of mesoscopic quantum systems may thus offer a way to probe the (coarsegrained) matrix elements of applied perturbations. A similar proposal to extract matrix structures from dynamics can also be found in the recent work [49] using periodic driving and working in the regime of weak perturbations governed by the exponential law (36). Our present approach can be considered complementary in that it avoids time-dependent manipulations and extends to significantly stronger perturbations. Given the important role of matrix elements in the energy eigenbasis for the dynamics in general and for questions of equilibration and thermalization (e.g. the eigenstate thermalization hypothesis) in particular, this sets up new possibilities to explore the underlying mechanisms by means of time series analysis.